有限元matlab仿真
最近写有有限元作业，用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