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

k=k+H(j)*H(l)*B'*D*B*det(J); end end

G=zeros(8,2*n); %初始总刚变换矩阵 G(1,2*elem(i,1)-1)=1; %总刚变换矩阵 G(2,2*elem(i,1))=1; G(3,2*elem(i,2)-1)=1; G(4,2*elem(i,2))=1; G(5,2*elem(i,3)-1)=1; G(6,2*elem(i,3))=1; G(7,2*elem(i,4)-1)=1; G(8,2*elem(i,4))=1;

K=K+G'*k*G; %总体刚度矩阵合成 end KK=K;

b=size(w,1); for i=1:b

K(2*w(i)-1,2*w(i)-1)=1e20; K(2*w(i),2*w(i))=1e20;

end %置大数法 f=zeros(2*n,1); %初始载荷矩阵 f(10)=-10e3; %加载荷10kN U=K\f; %节点位移

for i=1:m %将每个单元各个节点位移集

合

u(:,i)=[U(2*elem(i,1)-1);U(2*elem(i,1));U(2*elem(i,2)-1);U(2*elem(i,2));U(2*elem(i,3)-1);U(2*elem(i,3));U(2*elem(i,4)-1);U(2*elem(i,4))]; end

for i=1:m %求单元应力 syms Ks Et x y I1 I2 a b B; e=[1,1;-1,1;-1,-1;1,-1]; for j=1:4

N(j)=0.25*(1+e(j,1)*Ks)*(1+e(j,2)*Et); end x=0;y=0; for j=1:4

x=x+N(j)*X(elem(i,j),1); y=y+N(j)*X(elem(i,j),2);