结构有限元与matlab,有限元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