xs=diff(x,s); xt=diff(x,t); ys=diff(y,s); yt=diff(y,t); J=xs*yt-ys*xt; for i=1:8 Ns(i)=diff(N(i),s); Nt(i)=diff(N(i),t); end for i=1:8 g(i)=yt*Ns(i)-ys*Nt(i); h(i)=xs*Nt(i)-xt*Ns(i); end B=sym(zeros(3,16)); for i=1:8 B(1,2*i-1)=B(1,2*i-1)+g(i); B(2,2*i)=B(2,2*i)+h(i); B(3,2*i-1)=B(3,2*i-1)+h(i); B(3,2*i)=B(3,2*i)+g(i); end Bnew=simplify(B); Jnew=simplify(J); BD=transpose(Bnew)*D*Bnew/Jnew; r=int(int(BD,t,-1,1),s,-1,1); z=hou*r; k=double(z); 〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓 function K=kk(no,k,K) % no表示相应单元在整体坐标系中的编号,k为该单元的刚度矩阵 % K返回结构的整体刚度矩阵 for i=1:8 for j=1:8 K(2*no(i)-1,2*no(j)-1)=K(2*no(i)-1,2*no(j)-1)+k(2*i-1,2*j-1); K(2*no(i)-1,2*no(j))=K(2*no(i)-1,2*no(j))+k(2*i-1,2*j); K(2*no(i),2*no(j)-1)=K(2*no(i),2*no(j)-1)+k(2*i,2*j-1); K(2*no(i),2*no(j))=K(2*no(i),2*no(j))+k(2*i,2*j); end end (责任编辑:泉水) |