Matlab梁单元有限元编程：铁木辛柯梁VS欧拉梁讲解

7月前浏览4244

本文摘要（由ai生成）：

这篇文档主要讲解了二维和三维空间梁结构的 Matlab 有限元编程，包括欧拉梁单元和铁木辛柯梁单元。作者首先介绍了两种梁单元的基本理论，然后详细推导了它们的有限元离散方程，包括形函数、应变矩阵和刚度矩阵。最后，作者通过自编有限元 Matlab 程序，对比分析了不同单元对不同跨高比梁结构的计算结果，并探讨了铁木辛柯梁单元的剪切自锁问题。

导读：2月6日，仿真秀2023年度用户好评的“受用户好评的仿真好内容”名单公示。由仿真秀专栏作者SimPC创作的《Matlab有限元编程从入门到精通20讲：快速获得各典型有限元案例的Matlab代码》，荣获仿真秀2023年度“受用户好评的仿真好内容”，可喜可贺！诚邀用户关注和分享。

以下是正文

一、写在文前

梁在工程中应用广泛，是重要的结构构件。从几何上看，梁是任意截面形状的承受横向力的杆状结构，与杆的区别仅在于二者承受的载荷不同。在梁结构中，不同的梁固接在一起，既能传递力，又能传递力矩。

本文针对二维和三维空间梁结构的matlab有限元编程进行讲解，涉及的梁单元类型有欧拉梁单元和铁木辛柯梁单元。重点讲解二者的基本力学假定、适用范围、对应的三大类方程的建立、有限元离散方程的建立（包括形函数、刚度矩阵推导等）以及通过Matlab编程的实现上述两类梁单元静力分析求解和模态分析求解，获得梁结构的位移、剪力、弯矩图，以及模态各阶频率和振型，并探讨了铁木辛柯梁单元剪切自锁问题，并且对比了欧拉梁、铁木辛柯梁（完全积分）、铁木辛柯梁（减缩积分）的计算结果。

二、欧拉梁&铁木辛柯梁基本理论

1、欧拉-伯努利梁 Euler-Bernoulli Beam

欧拉-伯努利梁适用的前提条件是发生小变形、线弹性范围内、材料各向同性、等截面，其变形特征在于只有弯曲形变、横截面没有产生切应变；梁受力发生变形时，横截面依然为一个平面，且始终垂直于中性轴。受力方向垂直于中性轴。仅一个独立的变量v，即，垂直方向的位移。由于它忽略了切应变的效果，计算出的梁的变形量低于现实中梁的变形量；因此适用于细长梁。其平衡方程、几何方程、物理方程为：

2、铁木辛柯梁 Timoshenko Beam

铁木辛柯梁适用的前提条件是发生小变形、线弹性、各向同性的材料、等截面，其变形特征在于梁产生弯曲变形和梁的横截面产生切应变；梁受力发生变形时，横截面依然为一个平面，但不再垂直于中性轴。存在两个变量独立的变量，v垂直方向的位移和 θ横截面旋转角。由于它考虑了切应变的效果，计算出的梁的变形量，接近于现实梁的变形量；因此适用于短梁。其平衡方程、几何方程、物理方程为（公式参考https://zhuanlan.zhihu.com/p/134198227）：

最终的运动学方程

值得一提的是，两种梁模型也分别对应了2维的Kirchhoff板和Mindlin板模型，类似的思路可以直接平移到板理论去分析2D问题。

三、欧拉梁&铁木辛柯有限元离散

1、欧拉-伯努利梁有限元离散

接下来基于欧拉-伯努利梁理论推导平面纯弯梁的有限元方程，该理论适用于细长梁。在局部坐标系下，平面纯弯梁的每个节点有两个自由度（DOF），分别是方向的挠度和平面内的转角。因此，每个纯弯梁单元有4个自由度。当然也可以将轴向自由度考虑进去，直梁单元拉压与弯曲的变形和应力状态相互之间不耦合，因此可以直接叠加轴向刚度，相当于轴力杆单元和梁单元的简单叠加，就是将二者的刚度矩阵的叠加。本文只讲解4自由度的梁单元有限元方程的推导。

（1）形函数

考虑如图1所示的梁单元，单元长度2a，节点编号为 1和2。x方向沿梁的轴线方向，局部坐标系的坐标原点在梁单元的中点。为了推导梁单元的有限元方程，首先需要推导梁单元的形函数。由于梁单元有4个自由度，所以需要有4个形函数。在推导形函数时，通常会使用一种称为自然坐标系( naturalcoordinate system )的无量纲局部坐标系，以便将坐标系的取值范围变换到[0, 1]或者[-1, 1]之间。本例中自然坐标系的取值范围为[-1, 1]。

图1 局部坐标系下的梁单元

自然坐标系与局部坐标系的关系为

为了推导自然坐标系下的四个形函数，我们可以将挠度写成

的三次多项式

写作矩阵的形式

在小变形情况下，转角

可由挠度微分得到

引入边界条件：

可得

求解

之后可以得到

其中，N是形函数矩阵

（2）应变矩阵和刚度矩阵

欧拉-伯努利梁理论中，梁的横截面始终垂直于中性轴。将2.1.1节的用形函数表示的位移表达式带入到1.1节的几何方程中，得

得到应变矩阵后，可以求出单元刚度矩阵

扩展到三维空间梁单元，自由度按以下顺序排列：

首结点 1轴平动、2轴平动、3轴平动、绕1轴旋转、绕2轴旋转、绕3轴旋转；末结点 1轴平动、2轴平动、3轴平动、绕1轴旋转、绕2轴旋转、绕3轴旋转，共12个自由度。因此对应的刚度矩阵为

2、铁木辛柯梁有限元离散

参考《Matlab有限元结构动力学分析与工程应用-徐斌》梁单元的介绍，Timoshenko 梁单元中，横向位移v和转角是独立的，可各自独立插值，有

式中，N1和N2为线性形函数

根据1.1节的几何方程和物理方程，将上述位移表达式带入之后，可得

因此梁的刚度矩阵为

3、铁木辛柯梁的剪切锁死现象

上述的 Timoshenko梁刚度矩阵用于扁梁（高宽比小）时会产生剪切锁死现象。原因在于v和采用了同阶插值造成截切刚度过大。为了避免这种现象，在建立单元刚度矩阵的时候可采用减缩积分的方法建立。笔者通过自编有限元matlab程序验证对比了欧拉梁单元、铁木辛柯梁单元（完全积分）、铁木辛柯梁单元（减缩积分）三种单元对不同跨高比的梁结构的计算结果。

（1）L=5000,H=700,B=300的梁的计算结果

（2）L=5000,H=2700,B=300的梁的计算结果

（3）L=5000,H=7,B=300的梁的计算结果

可见，对于（1）所示以弯曲变形为主的梁结构，三种单元均能给出一个较为准确的结果；对于（2）所示的深梁，剪切变形不可忽略，导致欧拉梁与铁木辛柯梁的计算结果相差较大，而且两种铁木辛柯梁均给出较精确的结果；对于（3）所示的扁梁，以受弯变形为主，因此欧拉梁和减缩积分的铁木辛柯梁均能给出合理结果，但是由于完全积分的铁木辛柯梁会发生严重的剪切自锁现象，因此与欧拉梁和减缩积分的铁木辛柯梁的结果相差较大。