基于MATLAB的有限元结构分析

浮 涛

（郑州市交通规划勘察设计研究院，河南 郑州 450000）

0 引言

有限元法是当今工程分析中应用最为广泛的数值分析方法。在以往的结构有限元分析编程实践中，通常采用BASIC、FORTRAN、C、C++语言等进行编程，相比较而言，MATLAB 被视为有效快捷的工程分析工具，其优点是可将编程计算、数据分析及数据图形化成功地集成在工作环境中，被广泛应用于多领域。使用者可编译M 文件或者运用软件自身提供的强大工具箱来进行动态仿真，可求解出一系列工程问题，且求出的数值结果具有可视化、图形功能完善等特点。

在求解普通结构的工程分析中，使用较多的二维有限元是平面梁单元、平面桁架单元和平面刚架单元，三种单元都具有两种坐标系，分别为总体坐标系和局部坐标系。

平面梁单元与平面刚架单元的参数有E、I、A、L（E代表弹性模量，I代表截面惯性矩，A代表截面面积，L代表单元长度），每个平面梁单元与平面刚架单元均有2个节点，并且相对于总体坐标系的X轴正向逆时针的倾斜角为θ。其中平面梁单元有4 个自由度（每个节点有2 个自由度即1 个位移自由度和1 个转角自由度）；平面刚架单元有6个自由度（每个节点有3个自由度即2个位移自由度和一个转角自由度）[1-3]。

平面桁架单元的参数有E、A、L（E代表弹性模量，A代表截面面积，L代表单元长度），每个平面桁架单元有2个节点，并且相对于总体坐标系的X轴正向逆时针的倾斜角为θ。平面桁架单元有4 个自由度（每个节点有2个自由度即2个位移自由度）。

上述平面单元均约定位移向上为正，转角逆时针方向为正。

本文采用MATLAB 进行有限元分析的方法，并对某平面桁架结构、平面刚架结构进行实例分析研究。

1 MATLAB有限元分析

1.1 MATLAB有限元分析的基本步骤

结构的离散化，写出单元刚度矩阵，集成整体刚度矩阵，引入边界条件，求解方程组，结果处理。

1.2 编写M函数文件

1.2.1 平面桁架单元的M函数文件

（1）PTEStiffness（E,A,L,θ）—计算平面桁架单元（E,A,L,θ）的单刚矩阵。

（2）PTAssemble（K,k,i,j）—将连接节点i、j的平面桁架单元的单刚矩阵k集成到整刚矩阵K。

（3）PTEForce（E,A,L,θ,u）—计算单元节点力矢量。（4）PTEStress（E,L,θ,u）—计算单元应力。

1.2.2 平面刚架单元的M函数文件

（5）PFEStiffness（E,A,I,L,θ）—计算平面刚架单元(E,A,I,L,theta)的单刚矩阵。

（6）PFAssemble（K,k,i,j）—将连接节点i、j的平面刚架单元的单刚矩阵k集成到整刚矩阵K。

（7）PFEForces（E,A,I,L,θ,u）—计算平面刚架单元的节点力矢量。

（8）PFEADiagram（f,L）—绘制节点力矢量为f和长度为L的单元轴力图。

（9）PFESDiagram（f,L）—绘制节点力矢量为f和长度为L的单元剪力图。

（10）PFEMDiagram（f,L）—绘制节点力矢量为f和长度为L的单元弯矩图。

在MATLAB软件集成环境中，若进行平面梁单元结构分析，使用（1）～（4）这4个M函数文件，根据实际结构给出函数所需要的参数直接调用即可；若进行平面桁架单元结构分析，使用（5）～（10）这6个M函数文件，根据实际结构给出函数所需要的参数直接调用即可[4]。

MATLAB 软件的使用方法，本文不再详细叙述，用MATLAB 软件写出矩阵方程和在MATLAB 集成环境中编写并使用M函数文件请参阅文献[5]。

2 实例分析

2.1 平面桁架结构

如图1 的桁架，已知各杆E及A均相同，且E=210GPa，A=1E-2m2，L1=400mm，L2=300mm，L3=500mm，L4=400mm，求解各杆件内力。

图1 桁架结构图

根据M函数文件PTEStiffness（）求出单元的刚度矩阵k1、k2、k3、k4，该结构有四个节点，所以整体刚度矩阵为8×8。由于该结构有四个桁架单元，故需要四次调用函数文件PTAssemble（），求得整体刚度矩阵。

被约束的自由度有：U1x=U1y=U2y=U4x=U4y=0。则活动自由度编号为3，5，6。活动自由度对应的节点载荷F3=20000N，F5=0，F6=-25000N，采用高斯消去法进行求解，最后求得节点2的水平位移3.81e-3（mm）,节点3的水平位移7.94e-4（mm），节点3的竖向位移3.125e-3（mm）。建立结构的节点位移矢量U,然后计算结构的节点力矢量F。由此可求得节点1的水平支反力-15.833kN（方向向左），垂直支反力3.125kN（方向向上），节点2的水平支反力20kN（方向向右），垂直支反力21.875kN（方向向上），节点4的水平支反力4.167kN（方向向左），垂直支反力0，满足力的平衡条件。

建 立 单 元 位 移 矢 量u1、u2、u3和u4，然 后 调 用PTEForces（）求解各杆轴力，调用PTEStress（）求解各杆应力，见表1所示。

表1 各单元轴力及应力[4]

2.2 平面刚架结构

如图2的刚架，已知各杆E、I及A均相同，且A=1E-2m2，I=3E-5m4，L1=6m，L2=8m，L3=6m，E=210GPa，求节点2和3的位移和转角，并采用MATLAB绘出单元2的剪力图和弯矩图。

图2 刚架结构图

根据M函数文件PFEStiffness（）求出单元的刚度矩阵k1、k2、k3，该结构有四个节点，所以整体刚度矩阵为12×12。由于该结构有三个刚架单元，故需要三次调用函数文件PFAssemble（），求得整体刚度矩阵。

被约束的自由度有：U1x=U1y=R1z=0；U4x=U4y=R4z=0。则活动自由度编号为4，5，6，7，8，9。活动自由度对应的节点载荷F4=-30000N，F5=0，F6=0，F7=0，F8=0，F9=-20000kN·m，采用高斯消去法进行求解，最后求得节点2 的水平位移-0.0611m（向左），竖向位移0，转角位移0.0078rad（逆时针）；节点3 的水平位移-0.0610m（向左），竖向位移0，转角位移0.0043rad（逆时针）。

建立结构的节点位移矢量U,然后计算结构的节点力矢量F。由此可求得节点1 的水平支反力13.187kN（方向向右），垂直支反力7.158kN（方向向上），弯矩为47.753kN·m（顺时针）节点4 的水平支反力16.813kN（方向向右），垂直支反力7.158kN（方向向下），弯矩为54.983kN·m（顺时针），满足力的平衡条件。

表2 建立单元位移矢量u1、u2和u3，然后调用PFE Forces（）函数求出单元矢量f1、f2和f3。

表2 单元轴力、剪力、弯矩值[4]（单位：kN）

最后调用PFEADiagram（）、PFESDiagram（）以及PFEMDiagram（）可分别绘出单元的轴力图、剪力图和弯矩图。单元2 的轴力图、剪力图和弯矩图见图3、图4、图5。

图3 单元2轴力图（单位：N）

图4 单元2剪力图（单位：N）

图5 单元2弯矩图（单位：N）

3 结束语

本文以平面桁架结构和平面刚架结构为例，利用MATLAB 进行有限元结构分析，经实践证明，该方法简便、可行、实用。本文虽仅从平面的角度为例验证其可行性，但其原理和方法可推广到其他单元，如板壳、实体单元的刚度矩阵，甚至三维空间和非线性刚度矩阵之中。同时，该方法可为相关工程人员继续深入研究提供一定参考依据。