平面四节点等参单元matlab实现(6)

end

J1=jacobian([x;y],[Ks;Et]); J=J1';

for j=1:4

I1=diff(N(j),Ks); I2=diff(N(j),Et); C=(J^-1)*[I1;I2]; a=C(1); b=C(2);

B(1,2*j-1)=a; B(1,2*j)=0; B(2,2*j-1)=0; B(2,2*j)=b; B(3,2*j-1)=b;

B(3,2*j)=a; %以上同前面部分为得到B阵 end

D=(E/(1-NU*NU))*[1,NU,0;NU,1,0;0,0,(1-NU)/2]; w=D*B*u(:,i);

w1=subs(w,{Ks,Et},{1,1}); %求单元上各节点的应力 sigma1(:,i)=double(w1); w2=subs(w,{Ks,Et},{-1,1}); sigma2(:,i)=double(w2); w3=subs(w,{Ks,Et},{-1,-1}); sigma3(:,i)=double(w3); w4=subs(w,{Ks,Et},{1,-1}); sigma4(:,i)=double(w4); end

c=[24,29,47,58,78,79,137,149,160,166,186]'; %如截图选取圆半径方向的节点号 d=[166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184

,185]';%圆周方向选择的节点号 e=size(c,1); f=size(d,1);

sigmar=zeros(3,e); %r相同角度不同节点应力矩阵 sigmat=zeros(3,f); %角度不同r不同节点应力矩阵 msigmar=zeros(1,e); %半径方向节点处的mises应力 msigmat=zeros(1,f); %圆周方向节点处的mises应力

for i=1:e %绕节点平均 g=0;

for j=1:m %判断节点在单元的那个位置并加上相应的应

力值