平面单元应力求解

E=210000000;
v=0.28;
t=0.01;
%%%%%%各单元节点坐标向量,各行分别表示x1 y1 x2 y2 x3 y3
GDYJDZB=[0,1,0,2/3,1/3,2/3;
0,1,1/3,2/3,1/3,1;
1/3,1,1/3,2/3,2/3,2/3;
1/3,1,2/3,2/3,2/3,1;
2/3,1,2/3,2/3,1,2/3;
2/3,1,1,2/3,1,1;
0,2/3,0,1/3,1/3,1/3;0,2/3,1/3,1/3,1/3,2/3;
1/3,2/3,1/3,1/3,2/3,1/3;1/3,2/3,2/3,1/3,2/3,2/3;
2/3,2/3,2/3,1/3,1,1/3;2/3,2/3,1,1/3,1,2/3;
0,1/3,0,0,1/3,0;0,1/3,1/3,0,1/3,1/3;
1/3,1/3,1/3,0,2/3,0;1/3,1/3,2/3,0,2/3,1/3;
2/3,1/3,2/3,0,1,0;2/3,1/3,1,0,1,1/3];
%%%%%%位移编号,各行分别为u1 v1 u2 v2 u3 v3
WYBH=[0,0,0,0,7,8;0,0,7,8,1,2;
1,2,7,8,9,10;1,2,9,10,3,4;
3,4,9,10,11,12;3,4,11,12,5,6;
0,0,0,0,13,14;0,0,13,14,7,8;
7,8,13,14,15,16;7,8,15,16,9,10;
9,10,15,16,17,18;9,10,17,18,11,12;
0,0,0,0,19,20;0,0,19,20,13,14;
13,14,19,20,21,22;13,14,21,22,15,16;
15,16,21,22,23,24;15,16,23,24,17,18];
%%%%%%求各单元单元刚度矩阵与集成整体刚度矩阵
K=zeros(24);
for i=1:1:18
AA=[1,GDYJDZB(i,1),GDYJDZB(i,2);1,GDYJDZB(i,3),GDYJDZB(i,4);1,GDYJDZB(i,5),GDYJDZB(i,6)]
A=0.5*det(AA);
abc=det(AA)*inv(AA);%求矩阵AA各元素的代数余子式,得到ai bi ci
b=abc(:,2);c=abc(:,3)
B=1/det(AA)*[b(1),0,b(2),0,b(3),0;0,c(1),0,c(2),0,c(3);c(1),b(1),c(2),b(2),c(3),b(3)];
D=E/(1-v^2)*[1,v,0;v,1,0;0,0,(1-v)/2];
k=B'*D*B*t*A;%单元刚度矩阵
P=find(WYBH(i,:));
CC=WYBH(i,P);
K(CC,CC)=K(CC,CC) k(P,P);%整体刚度矩阵
end
WLXL=[0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]';%外力向量,规定:向右、向上为正方向
WY=K\WLXL; %求解位移,或者 WY=inv(K)*WLXL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DYJDL=[];DYYL=[];DYYB=[];
for i=1:1:18
U=[0,0,0,0,0,0]';
AA=[1,GDYJDZB(i,1),GDYJDZB(i,2);1,GDYJDZB(i,1),GDYJDZB(i,4);1,GDYJDZB(i,5),GDYJDZB(i,6)];
A=0.5*det(AA);
abc=det(AA)*inv(AA);%求矩阵AA各元素的代数余子式,得到ai bi ci
b=abc(:,2);c=abc(:,3);
B=1/det(AA)*[b(1),0,b(2),0,b(3),0;0,c(1),0,c(2),0,c(3);c(1),b(1),c(2),b(2),c(3),b(3)];
D=E/(1-v^2)*[1,v,0;v,1,0;0,0,(1-v)/2];
k=B'*D*B*t*A;%单元刚度矩阵
P=find(WYBH(i,:));
CC=WYBH(i,P);
U(P)=WY(CC);
yb=B*U;
DYYB=[DYYB yb];
yl=D*B*U;
DYYL=[DYYL yl];
f=k*U;
DYJDL=[DYJDL f];%各节点的力应平衡
end
%%%%%%%%%%%%%%节点处应力修正(取各单元均值)
jdljdy=[1,2,0,0,0,0;2,3,4,0,0,0;4,5,6,0,0,0;4,0,0,0,0,0;
1,7,8,0,0,0;1,2,3,8,9,10;3,4,5,10,11,12;5,6,12,0,0,0;
7,13,14,0,0,0;7,8,9,14,15,16;9,10,11,16,17,18;11,12,18,0,0,0;
13,0,0,0,0,0;13,14,15,0,0,0;15,16,17,0,0,0;17,18,0,0,0,0];
JDYL=[];
for i=1:1:16
P=find(jdljdy(i,:));
CC=jdljdy(i,P);
x=numel(P);
jdyl=sum(DYYL(:,CC),2)/x;
JDYL=[JDYL jdyl]
end