最近写有有限元作业，用matlab实现了整个的过程，顺便熟悉一下有限元的整体解题思路，感觉这确实是一种很好方法。下面对书上的一道题做仿真吧，最后结果如下。
题目：如图所示薄板

忽略自重，在力F的作用下，求各节点位移和各单元应力，已知

假定该问题为平面应力问题。

 这是一个正方形的薄板，上面两个顶点处被吊起来，在中心处受到一个向下的力。现分成四个单元（下面的两个顶点编号为1,2，上面个两个为3,4，中间的5）

执行程序得到如下结果

1节点的位移:（u，v） = （-1.42857e-006, -7.14286e-006）
2节点的位移:（u，v） = （1.42857e-006, -7.14286e-006）
3节点的位移:（u，v） = （0, 0）
4节点的位移:（u，v） = （0, 0）
5节点的位移:（u，v） = （0, -8.57143e-006）
第1单元:
应力为（750000，3.75e+006，-2.25e+006）
应变为（3.57143e-006，1.78571e-005，-2.14286e-005）
第2单元:
应力为（1.5e+006，-1.5e+006，6.98492e-010）
应变为（7.14286e-006，-7.14286e-006，1.01644e-020）
第3单元:
应力为（750000，3.75e+006，2.25e+006）
应变为（3.57143e-006，1.78571e-005，2.14286e-005）
第4单元:
应力为（0，9e+006，0）
应变为（0，4.28571e-005，0）

还可以实现各顶点移动的动画图，本来想做个gif，但发现getframe截不了屏，郁闷！！！只能在这贴几张图了

这个是不是就是ANSYS的有限元分析啊！

代码部分如下，显示动画的代码麻烦点，就不贴了，从代码也可以看到，好的数据结构的重要性：

main.m文件：

clc
clear all
close all

% 数据
choise = 1;
switch choise
    case 1
        % 作业第3题
        verts = [0 0;
            2 0;
            2 2;
            0 2;
            1 1;]*0.2;
        eleindex = [1 5 4;
            2 5 1;
            3 5 2;
            4 5 3];
        vertdisind = [1 1 1 1 0 0 0 0 0 1]';
        R = [0 0 0 0 0 0 0 0 0 -30000]';
    case 2
        % 书上例题
        verts = [0 2;
            0 1;
            1 1;
            0 0;
            1 0;
            2 0];
        eleindex = [5 2 4;
            6 3 5;
            2 5 3;
            3 1 2];
end
figure
patch('Faces', eleindex, 'Vertices', verts, 'FaceColor', 'r'); % 绘制数据
axis equal
axis off

% 构造FEM
hFEM.vertcount = size(verts, 1);
hFEM.verts = verts;
hFEM.elecount = size(eleindex, 1);
hFEM.eleindex = eleindex;
hFEM.mu = 0;
hFEM.E = 210e9;
hFEM.t = 0.01;
hFEM = FEM(hFEM);

% 求解有限元方程得各节点位移
hFEM = FEMSolve(hFEM, vertdisind, R);
for k = 1:hFEM.vertcount
    fprintf('%d节点的位移:（u，v） = （%g, %g）\n', k, hFEM.vertdis([2*k-1, 2*k]));
end

% 求各单元应力、应变
for k = 1:hFEM.elecount
    ind(1:2:5) = 2*hFEM.eleindex(k, :)-1;
    ind(2:2:6) = 2*hFEM.eleindex(k, :);
    stress = hFEM.eletruct(k).S * hFEM.vertdis(ind);
    strain = hFEM.eletruct(k).B * hFEM.vertdis(ind);
   
    fprintf('第%d单元:\n应力为（%g，%g，%g）\n', k, stress);
    fprintf('应变为（%g，%g，%g）\n', strain);
end

FEM.m文件：

% 计算FEM分析的各种矩阵
function hFEM = FEM(hFEM)
% 处理每个单元
for i = 1:hFEM.elecount
    hFEM.eletruct(i) = FEMele(hFEM, i);
end

% 计算合成刚度矩阵
hFEM.K = FEMCalcCombiningK(hFEM);
end

function hFEMele = FEMele(hFEM, ind)
hFEMele.ele = hFEM.verts(hFEM.eleindex(ind, :), :);
hFEMele.B = FEMCalcB(hFEMele);
hFEMele.D = FEMCalcD(hFEM);
hFEMele.S = FEMCalcS(hFEMele);
hFEMele.K = FEMCalcK(hFEM, hFEMele);
end

% 计算几何矩阵
function B = FEMCalcB(hFEMele)
ele = hFEMele.ele;
A(:, 2:3) = ele;
A(:, 1) = 1;
invA = inv(A);

B(1, 1:2:5) = invA(2, :);
B(2, 2:2:6) = invA(3, :);
B(3, 1:2:5) = invA(3, :);
B(3, 2:2:6) = invA(2, :);
end

% 计算应力矩阵
function D = FEMCalcD(hFEM)
mu = hFEM.mu;
E = hFEM.E;

D = [1 mu 0;
    mu 1 0;
    0 0 (1-mu)/2]*E/(1-mu^2);
end

% 计算应力变换矩阵
function S = FEMCalcS(hFEMele)
S = hFEMele.D*hFEMele.B;
end

% 计算单元刚度矩阵
function K = FEMCalcK(hFEM, hFEMele)
triarea = CalcTriangleArea(hFEMele.ele);
K = hFEM.t*triarea*hFEMele.B'*hFEMele.D*hFEMele.B;
end

% 计算三角形面积
function area = CalcTriangleArea(ele)
A(:, 2:3) = ele;
A(:, 1) = 1;
area = 0.5*det(A);
end

% 计算合成矩阵
function K = FEMCalcCombiningK(hFEM)
vertcount = hFEM.vertcount;
elecount = hFEM.elecount;
eleindex = hFEM.eleindex;

eleK = zeros(6, 6, elecount);
for k = 1:elecount
    eleK(:, :, k) = hFEM.eletruct(k).K;
end

K = zeros(vertcount*2);
for i = 1:vertcount
    for j = 1:vertcount
        for k = 1:elecount
            indi = find(eleindex(k, :)==i);
            indj = find(eleindex(k, :)==j);
            if ~isempty(indi) && ~isempty(indj)
                K(2*i-1:2*i, 2*j-1:2*j) = K(2*i-1:2*i, 2*j-1:2*j) + eleK(2*indi-1:2*indi, 2*indj-1:2*indj, k);
            end
        end
    end
end
end

FEMSolve.m文件：

% 求解有限元方程
function hFEM = FEMSolve(hFEM, vertdisind, R)
K = hFEM.K;
hFEM.vertdis = SolveFEM(vertdisind, R, K);
for i = 1:hFEM.elecount
    hFEM.eletruct(i).vertdis(1:2:5) = hFEM.vertdis(2*hFEM.eleindex(i, :) - 1); % u
    hFEM.eletruct(i).vertdis(2:2:6) = hFEM.vertdis(2*hFEM.eleindex(i, :)); % v
end
end

function vertdis = SolveFEM(vertdisind, R, K)
vertdis = vertdisind;

ind = find(vertdisind ~= 0);
Krule = K(ind, ind);
Rrule = R(ind);
vertdis(ind) = Krule\Rrule;
end