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

第1页共9页基于MTLAB的3节点三角形单元求解平面弹性问题一．引言本程序相比其他参考资料的程序，在操作方面得到大大简化，只需要一次输入相关力学参量、单元总数、节点总数、节点坐标、单元节点编号，利用此程序就能够很快得到所有单元的单元刚度矩阵，以及系统的整体刚度矩阵。为了便于后面计算单元的应力，在计算单元刚度矩阵时，MATLAB程序输出了弹性系数矩阵D和几何矩阵B。此外，本程序在引入位移边界条件求解刚度方程时，采用了化零置1法，避免了手工化简整体刚度矩阵的麻烦。当然，这在一定程度上加重了计算机的负担，增加了计算机的计算时间，但是这个延长的时间相对于手工化简刚度矩阵所需的时间来说，是微不足道的。程序需要提前输入的参数有：（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]中的一个简单实例进行验证。问题如下：如图所示为一矩形薄平板，在右端部受集中力F=100kN的作用，材料常数为：弹性模量E=1e7Pa，泊松比NU=1/3，板的厚度t=0.1m。试按平面应力问题计算各节点位移及支座反力。三．问题求解（1）结构的离散化与编号第2页共9页1）节点描述2)单元描述节点编号节点坐标x节点坐标y121220301400单元编号节点i节点j节点m12342321（2）计算单元刚度矩阵KE(:,:,1)=1.0e+006*0.2813000.1875-0.2813-0.187500.09380.18750-0.1875-0.093800.18750.37500-0.3750-0.18750.1875001.1250-0.1875-1.1250-0.2813-0.1875-0.3750-0.18750.65630.3750-0.1875-0.0938-0.1875-1.12500.37501.2188KE(:,:,2)=1.0e+006*0.2813000.1875-0.2813-0.187500.09380.18750-0.1875-0.093800.18750.37500-0.3750-0.18750.1875001.1250-0.1875-1.1250-0.2813-0.1875-0.3750-0.18750.65630.3750-0.1875-0.0938-0.1875-1.12500.37501.2188（3）组装整体刚度矩阵KK=1.0e+006*0.65630.3750-0.3750-0.1875-0.2813-0.1875000.37501.2188-0.1875-1.1250-0.1875-0.093800-0.3750-0.18750.6563000.3750-0.2813-0.1875-0.1875-1.125001.21880.37500-0.1875-0.0938-0.2813-0.187500.37500.65630-0.3750-0.1875-0.1875-0.09380.3750001.2188-0.1875-1.125000-0.2813-0.1875-0.3750-0.18750.65630.375000-0.1875-0.0938-0.1875-1.12500.37501.2188第3页共9页（4）求解位移及支反力delta=0.1877-0.8992-0.1497-0.84220000F=1.0e+004*0-0.50000-0.5000-2.0000-0.07022.00001.0702（5）计算单元应力stress(:,:,1)=1.0e+005*-0.8419-0.2806-1.5791stress(:,:,2)=1.0e+004*8.4187-2.8953-4.2094第4页共9页四．结果分析可以看出，计算得到的单元1的应力分量为σx=-84190Pa,σy=-28060Pa,τxy=-157910Pa,单元2的应力分量为σx=84187Pa,σy=-28953Pa,τxy=-42094Pa.此结果与参考文献[1]上给出的结果完全吻合。五．程序优点很多参考资料中的程序在求解单元刚度矩阵时，往往需要人工进行重复操作，例如一个系统具有n个单元，那么就需要人为的运行n次程序，并且每次都需要手动更改所需参数来求解单元刚度矩阵。在组装整体刚度矩阵时，也需要人为的运行n次程序，并且每次都需要手动更改所需参数来组装整体刚度矩阵。这样一来，使得操作极为不便，浪费大量时间。本程序相比其他参考资料的程序，在操作方面得到大大简化，只需要一次输入相关力学参量、单元总数、节点总数、节点坐标、单元节点编号，利用此程序就能够很快得到所有单元的单元刚度矩阵，以及系统的整体刚度矩阵。为了便于后面计算单元的应力，在计算单元刚度矩阵时，MATLAB程序输出了弹性系数矩阵D和几何矩阵B。此外，本程序在引入位移边界条件求解刚度方程时，采用了化零置1法，避免了手工化简整体刚度矩阵的麻烦。当然，这在一定程度上加重了计算机的负担，增加了计算机的计算时间，但是这个延长的时间相对于手工化简刚度矩阵所需的时间来说，是微不足道的。第5页共9页附录1.计算单元节点坐标%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%triangle2d3nodeeng%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%functionelement_node_gcoord=triangle2d3nodeeng(node_gcoord,element_node_number,element_number)%--------------------------------------------------------------------------%element_node_gcoord=triangle2d3nodeeng(node_gcoord,element_node_number)%该函数用于生成单元的节点坐标%element_node_gcoord返回单元的节点坐标%node_gcoord节点坐标%element_node_number单元节点编号%--------------------------------------------------------------------------forloopi=1:element_numberelement_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,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]=triangle2d3nodestiffness(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为平面应变问题%------------------------------------------------------------------forloopi=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([1xiyi;1xjyj;1xmym]);%计算各三角形单元的面积第6页共9页b1=yj-ym;b2=ym-yi;b3=yi-yj;c1=xm-xj;c2=xi-xm;c3=xj-xi;%计算各单元的几何矩阵BB(:,:,loopi)=1/(2*A(loopi))*[b10b20b30;0c10c20c3;c1b1c2b2c3b3];ifID==1E=E;NU=NU;elseifID==2E=E/(1-NU^2);NU=NU/(1-NU);end%形成平面问题的弹性系数矩阵DD=E/(1-NU^2)*[1NU0;NU10;00(1-NU)/2];%形成单元的刚度矩阵KE(:,:,loopi)=B(:,:,loopi)'*D*B(:,:,loopi)*A(loopi)*t;end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%3.形成整体刚度矩阵%%%%%%%%%%%%%%%%%%%triangle2d3nodeassembly%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%functionKK=triangle2d3nodeassembly(KE,node_number,element_number,element_node_number)%--------------------------------------------------%KK=triangle2d3nodeassembly(KE,node_number,element_number,element_node_number)%该函数用于组装整体刚度矩阵%KK返回系统的整体刚度矩阵%KE单元刚度矩阵%nod