b1=yj-ym; b2=ym-yi; b3=yi-yj; c1=xm-xj; c2=xi-xm; c3=xj-xi;
%计算各单元的几何矩阵B
B(:,:,loopi)=1/(2*A(loopi))*[b1 0 b2 0 b3 0;0 c1 0 c2 0 c3;c1 b1 c2 b2 c3 b3]; if ID==1 E=E; NU=NU; elseif ID==2 E=E/(1-NU^2); NU=NU/(1-NU); end
%形成平面问题的弹性系数矩阵D
D=E/(1-NU^2)*[1 NU 0;NU 1 0;0 0 (1-NU)/2]; %形成单元的刚度矩阵
KE(:,:,loopi)=B(:,:,loopi)'*D*B(:,:,loopi)*A(loopi)*t; end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3.形成整体刚度矩阵
%%%%%%%%%%%%%%%%%%%triangle2d3nodeassembly%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function KK=triangle2d3nodeassembly(KE,node_number,element_number,element_node_number) %--------------------------------------------------
%KK=triangle2d3nodeassembly(KE,node_number,element_number,element_node_number) %该函数用于组装整体刚度矩阵 %KK返回系统的整体刚度矩阵 %KE单元刚度矩阵 %node_numbre节点总数 %element_number单元总数
%element_node_number单元的节点编号
%-------------------------------------------------- KK=zeros(2*node_number,2*node_number); for loopi=1:element_number i=element_node_number(loopi,1); j=element_node_number(loopi,2); m=element_node_number(loopi,3); dof(1)=2*i-1; dof(2)=2*i; dof(3)=2*j-1; dof(4)=2*j; dof(5)=2*m-1; dof(6)=2*m;
第6页 共9页
k=KE(:,:,loopi); for n1=1:6 for n2=1:6
KK(dof(n1),dof(n2))=KK(dof(n1),dof(n2))+k(n1,n2); end end end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4.求解未知位移和未知支反力
%%%%%%%%%%%%%%%%%%%%%%triangle2d3noded_f%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [delta,F]=triangle2d3noded_f(KK,deltas,FS,node_number) %----------------------------------------------------------------
%该函数用于引入位移边界条件化简整体刚度矩阵和节点载荷列阵，并且求出未知位移和支反力 Tlta返回求解出的节点位移列阵 %F返回求解出的节点力列阵 %KK原整体刚度矩阵
Tltas节点位移列阵，里面含有未知节点位移 %FS节点力列阵，里面含有未知节点力 %node_number节点总数
%---------------------------------------------------------------- KK1=KK;
for loopi=1:2*node_number if (deltas(loopi,1)==0) KK1(:,loopi)=0; KK1(loopi,:)=0; KK1(loopi,loopi)=1; FS(loopi,1)=0; end end KKS=KK1;
delta=double(KKS\\FS); F=KK*delta;
%%%%%%%%%%%%%%%%%%%%%%%%%end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5.计算单元应变和单元应力
%%%%%%%%%%%%%%%%%%%%%%%%%triangle2d3nodestrainstress%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function
[strain,stress]=triangle2d3nodestrainstress(D,B,delta,element_node_number,element_number) %------------------------------------------------------------
%[strain,stress]=triangle2d3nodestrainstress(D,B,delta,element_node_number,element_number) %该函数用于计算各单元的应变和应力 %strain返回各单元的应变 %stress返回各单元的应力 %D弹性系数矩阵 %B单元的几何矩阵
第7页 共9页
Tlta单元的位移列阵 %element_number单元总数
%------------------------------------------------------------ for loopi=1:element_number
dof1=element_node_number(loopi,1); dof2=element_node_number(loopi,2); dof3=element_node_number(loopi,3);
delta1=[delta(2*dof1-1,1);delta(2*dof1,1);delta(2*dof2-1,1);delta(2*dof2,1);delta(2*dof3-1,1);delta(2*dof3,1)];
strain(:,:,loopi)=B(:,:,loopi)*delta1; stress(:,:,loopi)=D*strain(:,:,loopi); end
%%%%%%%%%%%%%%%%%%%%%%%%%%end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
第8页 共9页
参考文献
[1]王勖成 编著.有限单元法.北京：清华大学出版社，2003
[2]徐斌，高跃飞，余龙 编著.MATLAB有限元结构动力学分析与工程应用.北京：清华大学出版社，2009 [3]江见鲸，何放龙，何益斌，陆新征 编著.有限元法及其应用.北京：机械工业出版社，2006 [4]曾攀 编著.有限元基础教程.北京：高等教育出版社，2009
[5]李人宪 编著.有限元法基础（第2版）.北京：国防工业出版社，2004
第9页 共9页
