基于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） 结构的离散化与编号
第1页 共9页
1） 节点描述 2) 单元描述
节点编号节点坐标x节点坐标y121节点i220单元编号30112400
23（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.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.2188
KE(:,:,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.6563 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 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
第2页 共9页
节点j节点m3421 （4） 求解位移及支反力 delta =
0.1877 -0.8992 -0.1497 -0.8422 0 0 0 0 F =
1.0e+004 *
0 -0.5000 0 -0.5000 -2.0000 -0.0702 2.0000 1.0702
（5） 计算单元应力 stress(:,:,1) =
1.0e+005 *
-0.8419 -0.2806 -1.5791
stress(:,:,2) =
1.0e+004 *
8.4187 -2.8953 -4.2094
第3页 共9页
四． 结果分析
可以看出，计算得到的单元1的应力分量为σx =-84190Pa,σy=-28060Pa,
τxy=-157910Pa,单元2的应力分量为σx =84187Pa,σy=-28953Pa,τxy=-42094Pa.此结果与参考文献[1]上给出的结果完全吻合。
五． 程序优点
很多参考资料中的程序在求解单元刚度矩阵时，往往需要人工进行重复操作，例如一个系统具有n个单元，那么就需要人为的运行n次程序，并且每次都需要手动更改所需参数来求解单元刚度矩阵。在组装整体刚度矩阵时，也需要人为的运行n次程序，并且每次都需要手动更改所需参数来组装整体刚度矩阵。这样一来，使得操作极为不便，浪费大量时间。
本程序相比其他参考资料的程序，在操作方面得到大大简化，只需要一次输入相关力学参量、单元总数、节点总数、节点坐标、单元节点编号，利用此程序就能够很快得到所有单元的单元刚度矩阵，以及系统的整体刚度矩阵。为了便于后面计算单元的应力，在计算单元刚度矩阵时，MATLAB程序输出了弹性系数矩阵D和几何矩阵B。
此外，本程序在引入位移边界条件求解刚度方程时，采用了化零置1法，避免了手工化简整体刚度矩阵的麻烦。当然，这在一定程度上加重了计算机的负担，增加了计算机的计算时间，但是这个延长的时间相对于手工化简刚度矩阵所需的时间来说，是微不足道的。
第4页 共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_gcoord返回单元的节点坐标 %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,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为平面应变问题 %------------------------------------------------------------------ for loopi=1:element_number xi=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]); %计算各三角形单元的面积
第5页 共9页
