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

专栏导读

作者简介：工学博士，高级工程师，专注于工业软件算法研究
本文已收录于专栏：《有限元编程从入门到精通》本专栏旨在提供 1.以案例的形式讲解各类有限元问题的程序实现，并提供所有案例完整源码；2.单元类型包含：杆单元，梁单元，平面三角形单元，薄板单元，厚板单元，壳单元，四/六面体实体单元，金字塔单元等；3.物理场问题涉及：力学、传热学、电磁学及多物理场耦合等问题的稳态（静力学）和瞬态（动力学）求解。专栏旨在帮助有志于有限元工业软件开发的小伙伴，快速上手有限元编程，在案例中成长，摆脱按部就班填鸭式教学。
【所有专栏文章均提供视频教程和源码】《有限元编程从入门到精通30讲》视频教程与源码获取地址： https://www.bilibili.com/video/BV1kP4y1d7Zo
欢迎订阅专栏，订阅用户可私聊进入有限元编程交流群（知识交流、问题解答），并获赠丰厚的有限元相关学习资料（教材、源码、视频课）
专栏订阅地址：有限元编程从入门到精通_suoge223的博客-CSDN博客

文章目录

专栏导读

文章目录

一、写在文前

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

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

2、铁木辛柯梁 Timoshenko Beam

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

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

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

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

四、我的Matlab有限元编程精品课

一、写在文前

梁在工程中应用广泛，是重要的结构构件。从几何上看，梁是任意截面形状的承受横向力的杆状结构，与杆的区别仅在于二者承受的载荷不同。在梁结构中，不同的梁固接在一起，既能传递力，又能传递力矩。本文针对二维和三维空间梁结构的matlab有限元编程进行讲解，涉及的梁单元类型有欧拉梁单元和铁木辛柯梁单元。重点讲解二者的基本力学假定、适用范围、对应的三大类方程的建立、有限元离散方程的建立（包括形函数、刚度矩阵推导等）以及通过Matlab编程的实现上述两类梁单元静力分析求解和模态分析求解，获得梁结构的位移、剪力、弯矩图，以及模态各阶频率和振型，并探讨了铁木辛柯梁单元剪切自锁问题，并且对比了欧拉梁、铁木辛柯梁（完全积分）、铁木辛柯梁（减缩积分）的计算结果。

【本案例源码链接】：Matlab梁单元有限元编程 | 铁木辛柯梁 | 欧拉梁 | Matlab源码 | 理论文本 |深梁 |短梁

【本案例视频教程】：(试看)铁木辛柯梁的matlab有限元编程/欧拉-伯努利梁/刚度矩阵/剪切自锁/减缩积分/高斯积分/等参单元/剪切变形/弯曲/仿真/结构力学/数值计算/微分方程_哔哩哔哩_bilibili

【《Matlab有限元编程从入门到精通》合集课】：Matlab有限元编程从入门到精通/后处理/刚度矩阵/等参单元/高斯积分/Hammer积分/非线性/静力/动力/传热/三维/二维/梁板壳/实体/平面_哔哩哔哩_bilibili

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

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

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

2、铁木辛柯梁 Timoshenko Beam

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

最终的运动学方程

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

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

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

接下来基于欧拉-伯努利梁理论推导平面纯弯梁的有限元方程，该理论适用于细长梁。在局部坐标系下，平面纯弯梁的每个节点有两个自由度（DOF），分别是方向的挠度和平面内的转角。因此，每个纯弯梁单元有4个自由度。当然也可以将轴向自由度考虑进去，直梁单元拉压与弯曲的变形和应力状态相互之间不耦合，因此可以直接叠加轴向刚度，相当于轴力杆单元和梁单元的简单叠加，就是将二者的刚度矩阵的叠加。本文只讲解4自由度的梁单元有限元方程的推导。（1）形函数考虑如图1所示的梁单元，单元长度2a，节点编号为 1和2。x方向沿梁的轴线方向，局部坐标系的坐标原点在梁单元的中点。为了推导梁单元的有限元方程，首先需要推导梁单元的形函数。由于梁单元有4个自由度，所以需要有4个形函数。在推导形函数时，通常会使用一种称为自然坐标系( natural coordinate 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程序验证对比了欧拉梁单元、铁木辛柯梁单元（完全积分）、铁木辛柯梁单元（减缩积分）三种单元对不同跨高比的梁结构的计算结果。本案例已收录至《Matlab有限元编程从入门到精通20讲》强烈推荐学习者订阅，本课程已上架仿真秀和b站课堂。。

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

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

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

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

四、我的Matlab有限元编程精品课

以上就是笔者关注Matlab梁单元有限元编程铁木辛柯梁与欧拉梁对比分析，该内容已经收录在我的原创视频课程里面《Matlab有限元编程从入门到精通20讲》强烈推荐学习者订阅。本课程已上架仿真秀和b站课堂。

本课程为matlab有限元编程专题课，并配套的全部案例源码和讲义PPT/理论文本，旨在以案例的形式讲解各类有限元问题程序实现及算法原理，并提供完整Matlab源码供大家练习，案例源码均包含前后处理模块和求解器模块；单元类型包含：杆单元，梁单元，平面三角形单元，板壳单元，四/六面体实体单元等；物理场问题涉及：静力学、动力学、传热学、材料非线性、几何非线性、接触非线性等求解，内容丰富，不断更新完善。。
此外，笔者为所有订阅用户提供知识圈答疑服务和VIP用户交流群。并附赠课程相关资料等（仿真秀/b站平台支持自行开具电子发票）。
大家学习Matlab进行有限元编程过程中，强烈推荐大家学习以下资料（在文章末尾点赞或再看，截图发到仿真秀公众号，回复 SimPC ，我会24小时发给你以下资料）。