基于MTLAB的3节点三角形单元求解平面弹性问题

1、基于MTLAB的3节点三角形单元求解平面弹性问题一 引言 本程序相比其他参考资料的程序，在操作方面得到大大简化，只需要一次输入相关力学参量、单元总数、节点总数、节点坐标、单元节点编号，利用此程序就能够很快得到所有单元的单元刚度矩阵，以及系统的整体刚度矩阵。为了便于后面计算单元的应力，在计算单元刚度矩阵时，MATLAB程序输出了弹性系数矩阵D和几何矩阵B。 此外，本程序在引入位移边界条件求解刚度方程时，采用了化零置1法，避免了手工化简整体刚度矩阵的麻烦。当然，这在一定程度上加重了计算机的负担，增加了计算机的计算时间，但是这个延长的时间相对于手工化简刚度矩阵所需的时间来说，是微不足道的。程序需要提

2、前输入的参数有：（1） 节点坐标node_gcoord（2） 单元节点编号element_node_number（3） 节点总数node_number（4） 单元总数element_number（5） 弹性模量E（6） 泊松比NU（7） 厚度t（8） 平面问题性质指示参数ID，其中1为平面应力问题，2为平面应变问题该程序主要包括以下部分：1. 单元刚度矩阵程序：用于计算各单元的单元刚度矩阵；2. 整体刚度矩阵程序：用于组装整体刚度矩阵；3. 位移及支反力求解程序：用于求解未知位移和未知支反力4. 单元应力计算程序：用于计算各单元的单元应力二 问题描述为了验证程序的正确性，我们选取了参考文献1中

3、的一个简单实例进行验证。问题如下：如图所示为一矩形薄平板，在右端部受集中力F=100kN的作用，材料常数为：弹性模量E=1e7Pa，泊松比NU=1/3，板的厚度t=0.1m。试按平面应力问题计算各节点位移及支座反力。三 问题求解（1） 结构的离散化与编号1） 节点描述 2) 单元描述 （2） 计算单元刚度矩阵KE(:,:,1) = 1.0e+006 * 0.2813 0 0 0.1875 -0.2813 -0.1875 0 0.0938 0.1875 0 -0.1875 -0.0938 0 0.1875 0.3750 0 -0.3750 -0.1875 0.1875 0 0 1.1250 -0

4、.1875 -1.1250 -0.2813 -0.1875 -0.3750 -0.1875 0.6563 0.3750 -0.1875 -0.0938 -0.1875 -1.1250 0.3750 1.2188KE(:,:,2) = 1.0e+006 * 0.2813 0 0 0.1875 -0.2813 -0.1875 0 0.0938 0.1875 0 -0.1875 -0.0938 0 0.1875 0.3750 0 -0.3750 -0.1875 0.1875 0 0 1.1250 -0.1875 -1.1250 -0.2813 -0.1875 -0.3750 -0.1875 0.65

5、63 0.3750 -0.1875 -0.0938 -0.1875 -1.1250 0.3750 1.2188（3） 组装整体刚度矩阵KK = 1.0e+006 * 0.6563 0.3750 -0.3750 -0.1875 -0.2813 -0.1875 0 0 0.3750 1.2188 -0.1875 -1.1250 -0.1875 -0.0938 0 0 -0.3750 -0.1875 0.6563 0 0 0.3750 -0.2813 -0.1875 -0.1875 -1.1250 0 1.2188 0.3750 0 -0.1875 -0.0938 -0.2813 -0.1875 0

6、 0.3750 0.6563 0 -0.3750 -0.1875 -0.1875 -0.0938 0.3750 0 0 1.2188 -0.1875 -1.1250 0 0 -0.2813 -0.1875 -0.3750 -0.1875 0.6563 0.3750 0 0 -0.1875 -0.0938 -0.1875 -1.1250 0.3750 1.2188（4） 求解位移及支反力delta = 0.1877 -0.8992 -0.1497 -0.8422 0 0 0 0F = 1.0e+004 * 0 -0.5000 0 -0.5000 -2.0000 -0.0702 2.0000 1.

7、0702（5） 计算单元应力stress(:,:,1) = 1.0e+005 * -0.8419 -0.2806 -1.5791stress(:,:,2) = 1.0e+004 * 8.4187 -2.8953 -4.2094四 结果分析可以看出，计算得到的单元1的应力分量为x =-84190Pa,y=-28060Pa,xy=-157910Pa,单元2的应力分量为x =84187Pa,y=-28953Pa,xy=-42094Pa.此结果与参考文献1上给出的结果完全吻合。五 程序优点 很多参考资料中的程序在求解单元刚度矩阵时，往往需要人工进行重复操作，例如一个系统具有n个单元，那么就需要人为的运

8、行n次程序，并且每次都需要手动更改所需参数来求解单元刚度矩阵。在组装整体刚度矩阵时，也需要人为的运行n次程序，并且每次都需要手动更改所需参数来组装整体刚度矩阵。这样一来，使得操作极为不便，浪费大量时间。本程序相比其他参考资料的程序，在操作方面得到大大简化，只需要一次输入相关力学参量、单元总数、节点总数、节点坐标、单元节点编号，利用此程序就能够很快得到所有单元的单元刚度矩阵，以及系统的整体刚度矩阵。为了便于后面计算单元的应力，在计算单元刚度矩阵时，MATLAB程序输出了弹性系数矩阵D和几何矩阵B。 此外，本程序在引入位移边界条件求解刚度方程时，采用了化零置1法，避免了手工化简整体刚度矩阵的麻烦。

9、当然，这在一定程度上加重了计算机的负担，增加了计算机的计算时间，但是这个延长的时间相对于手工化简刚度矩阵所需的时间来说，是微不足道的。附录1. 计算单元节点坐标%triangle2d3nodeeng%function element_node_gcoord=triangle2d3nodeeng(node_gcoord,element_node_number,element_number)%-%element_node_gcoord=triangle2d3nodeeng(node_gcoord,element_node_number)%该函数用于生成单元的节点坐标 %element_node_g

10、coord返回单元的节点坐标%node_gcoord节点坐标%element_node_number单元节点编号%-for loopi=1:element_number element_node_gcoord(loopi,1)=node_gcoord(element_node_number(loopi,1),1); element_node_gcoord(loopi,2)=node_gcoord(element_node_number(loopi,1),2); element_node_gcoord(loopi,3)=node_gcoord(element_node_number(loopi,

11、2),1); element_node_gcoord(loopi,4)=node_gcoord(element_node_number(loopi,2),2); element_node_gcoord(loopi,5)=node_gcoord(element_node_number(loopi,3),1); element_node_gcoord(loopi,6)=node_gcoord(element_node_number(loopi,3),2); end%end%2. 计算单元刚度矩阵%triangle2d3nodestiffness%function KE,B,D,A=triangle

12、2d3nodestiffness(E,NU,t,element_number,element_node_gcoord,ID)%-%KE,B,D=triangle2d3nodestiffness(E,NU,t,element_number,element_node_gcoord,ID)%该函数用于计算单元的刚度矩阵%KE返回单元刚度矩阵%B返回各单元的几何矩阵%D返回平面弹性系数矩阵%输入弹性模量E，泊松比NU，厚度t%输入单元节点坐标文件element_node_gcooord%输入单元总数element_number%输入平面问题性质指示参数ID，其中1为平面应力问题，2为平面应变问题%-f

13、or loopi=1:element_numberxi=element_node_gcoord(loopi,1);yi=element_node_gcoord(loopi,2);xj=element_node_gcoord(loopi,3);yj=element_node_gcoord(loopi,4);xm=element_node_gcoord(loopi,5);ym=element_node_gcoord(loopi,6);A(loopi)=1/2*det(1 xi yi;1 xj yj;1 xm ym); %计算各三角形单元的面积b1=yj-ym;b2=ym-yi;b3=yi-yj;c

14、1=xm-xj;c2=xi-xm;c3=xj-xi;%计算各单元的几何矩阵BB(:,:,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-NU2); NU=NU/(1-NU);end%形成平面问题的弹性系数矩阵D D=E/(1-NU2)*1 NU 0;NU 1 0;0 0 (1-NU)/2;%形成单元的刚度矩阵KE(:,:,loopi)=B(:,:,loopi)'*D*B(:,:,loopi)*A(loopi)*t;end%

15、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单元

16、的节点编号%-KK=zeros(2*node_number,2*node_number);for loopi=1:element_numberi=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;k=KE(:,:,loopi);for n1=1:6 for n2=1:6 KK(dof(n1),dof(n2)=KK(dof(

17、n1),dof(n2)+k(n1,n2); endendend%end%4. 求解未知位移和未知支反力%triangle2d3noded_f%function delta,F=triangle2d3noded_f(KK,deltas,FS,node_number)%-%该函数用于引入位移边界条件化简整体刚度矩阵和节点载荷列阵，并且求出未知位移和支反力%delta返回求解出的节点位移列阵%F返回求解出的节点力列阵%KK原整体刚度矩阵%deltas节点位移列阵，里面含有未知节点位移%FS节点力列阵，里面含有未知节点力%node_number节点总数%-KK1=KK;for loopi=1:2*no

18、de_number if (deltas(loopi,1)=0) KK1(:,loopi)=0; KK1(loopi,:)=0; KK1(loopi,loopi)=1; FS(loopi,1)=0; endend KKS=KK1;delta=double(KKSFS);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)%该函数用于计算各