matlab 有限元分析

clc;
clc;
t=1;%单元厚度
E=2e11;%弹性模量
u=0.3;%泊松比
F1=50;F2=100;F3=25;F4=50;%外力
D=(E/(1-u^2))*[1 u 0;u 1 0;0 0 (1-u)/2];
%各结点坐标
no1=[0 -0.5];no2=[0.5 -0.5];no3=[1 -0.5];no4=[1.5 -0.5];no5=[2 -0.5];
no6=[2 0];no7=[2 0.5];no8=[1.5 0.5];no9=[1 0.5];no10=[0.5 0.5];
no11=[0 0.5];no12=[0 0];no13=[0.5 0];no14=[1 0];no15=[1.5 0];
%各个单元对应的结点坐标，顺序为 i,j,m
ele1=[1 2 12];ele2=[12 2 13];ele3=[12 13 11];ele4=[11 13 10];
ele5=[2 3 13];ele6=[13 3 14];ele7=[13 14 10];ele8=[10 14 9];
ele9=[3 4 14];ele10=[14 4 15];ele11=[14 15 9];ele12=[9 15 8];
ele13=[4 5 15];ele14=[15 5 6];ele15=[15 6 8];ele16=[8 6 7];
%各个单元对应的结点 i,j,m
eleno1=[no1;no2;no12];eleno2=[no12;no2;no13];eleno3=[no12;no13;no11];
eleno4=[no11;no13;no10];
eleno5=[no2;no3;no13];eleno6=[no13;no3;no14];eleno7=[no13;no14;no10];
eleno8=[no10;no14;no9];
eleno9=[no3;no4;no14];eleno10=[no14;no4;no15];eleno11=[no14;no15;no9];
eleno12=[no9;no15;no8];
eleno13=[no4;no5;no15];eleno14=[no15;no5;no6];eleno15=[no15 ;no6 ;no8];
eleno16=[no8;no6;no7];
%各个单元劲度矩阵
for ele =1:16%单元从 1 到 16
for i =1:3
eval(['x' num2str(i) '=eleno' num2str(ele) '(i,1)']);
%i,j,m 坐标顺序用 1,2,3 代替
eval(['y' num2str(i) '=eleno' num2str(ele) '(i,2)']);

end
Aa=[1,x1,y1;1,x2,y2;1,x3,y3];
A=det(Aa)/2;%各个单元的面积 A
b1=y2-y3;b2=y3-y1;b3=y1-y2;
c1=x3-x2;c2=x1-x3;c3=x2-x1;
B1=[b1,0;0,c1;c1,b1];
B2=[b2,0;0,c2;c2,b2];
B3=[b3,0;0,c3;c3,b3];
B=[B1,B2,B3]/(A*2);
k=B'DB.*t.*A;
eval(['k' num2str(ele) '=k']);%单元 1 到单元 16 的劲度矩阵
eval(['BB' num2str(ele) '=B']);%单元 1 到单元 16 的 B 矩阵
k=mat2cell(k,[2,2,2],[2,2,2]);
for j=1:3
for l=1:3
eval(['k' num2str(ele) num2str(j) num2str(l) '=k{j,l}']);
%第一个数表示单元第二和第三个数表示位置
end
end
end
K=zeros(30,30);
K=mat2cell(K,[2,2,2,2,2,2,2,2,2,2,2,2,2,2,2],[2,2,2,2,2,2,2,2,2,2,2,2,2,2,2]);
for i=1:15
for j=1:15
for e=1:16

eval(['elem=ele' num2str(e)]);
for m=1:3
if elem(m)==i
for n=1:3
if elem(n)==j
eval(['kk=k' num2str(e) num2str(m) num2str(n)]);
K{i,j}=K{i,j}+kk;
end
end
end
end
end
end
end
%结点 x,y 方向的受力
F=[F1;0;
F2;0;
F2;0;
F2;0;
0;F