三维铁木辛柯梁单元Matlab有限元编程案例分析
导读：上一篇《Matlab梁单元有限元编程：铁木辛柯梁VS欧拉梁讲解》针对二维和三维空间梁结构的matlab有限元编程进行讲解，涉及的梁单元类型有欧拉梁单元和铁木辛柯梁单元。今天笔者继续讲解三维铁木辛柯梁单元Matlab有限元编程技术，本文案例通过matlab实现了三维铁木辛柯梁单元和弹性支撑单元的有限元编程。
一、引言
梁单元为三维铁木辛柯梁单元，采用 Timoshenko 梁理论，通过截面剪切系数考虑梁截面的剪切变形影响。三维铁木辛柯梁单元是由3 个节点构成的直梁单元，其示意图见图 1-1。其中，o-xyz 为梁单元局部坐标系;节点i、j为单元物理节点，用于确定单元边界，为必选节点:节点k为单元梁截面主轴方向节点，用于确定梁截面主轴Z 的指向，为可选节点，当不设置 3 号节点时，梁的截面主轴将按照梁截面绕梁轴线的转角进行确认。三维铁木辛柯梁单元每个节点包含 6 个位移自由度，分别为沿单元局部坐标x、y、z轴的平动自由度 u、v、w，以及绕单元局部坐标
轴的转动自由度.
图1 铁木辛科梁单元自由度示意图
程序实现的案例为列车轮轴的静力分析，将车轴简化如图2所示。将轮对主轴简化为圆截面的空间三维Timoshenko梁单元如图3所示，将车轮简化为轴承支撑单元。车轴为阶梯轴，按照截面大小不同，将车轴划分为几个不同截面面积的梁单元组合，另外为了得到轮轴过盈配合连接处精确的位移，在车轮车轴连接处建立多个节点，连接处节点单元均使用轴承支撑单元，即弹性支撑单元。因此本代码还涉及弹性支撑单元的有限元编程。
图2
图3轮对Timoshenko梁模型
二、三维铁木辛柯梁几何方程与物理方程
三维Timoshenko 梁理论用于描述当梁的横截面尺寸与其长度相比不小时短梁的行为，具体包括轴向拉压、横向剪切、扭转和弯曲的三维 Timoshenko 梁单元的理论公式。
2.1 几何方程
如图1所示的一个长度为l的两节点梁单元，每个节点六个自由度，在局部坐标系下，梁的轴向与x轴重合。单元横截面上任意点的位移向量d=[U V W]可以表示为
（1）
式中u为梁轴线处沿轴向位移，v,w为轴线处横向两个方向的位移，分别由弯曲和截切变形两部分组成
（2）
横向两个方向的位移v,w与弯曲和剪切变形的关系为
（3）
和
时xy和xz平面的剪切应变，
和
为两个方向上的弯曲变形引起的转角。因此铁木辛柯梁的应变与位移的关系为
（4）
请注意，扭转过程中的轴向翘曲位移被忽略。
2.2 物理方程
根据弹性力学胡克定律，以及铁木辛柯梁理论对于梁横向剪应变的假定，应力应变关系如下
（5）
三、三维铁木辛柯梁形函数
基于最小位能原理的考虑剪切影响的梁单元和不考虑剪切影响的梁单元相同，仍以挠度w和转角为结点参数，但在泛函中引人了剪切应变能的影响。在引人剪切影响的方法上有以下两种方案: 1.在经典梁单元基础上引入剪切变形影响；2.挠度和截面转动各自独立插值。后者会引起剪切自锁现象，具体可参考徐斌的《MATLAB有限元结构动力学分析与工程应用》中关于二维铁木辛柯梁的理论，这样构造的单元是C0型单元，因为挠度和转动采用同阶插值表达式，所以剪切应变中dw/dx和两项是不同阶的。本理论模型采用方法1，在经典的欧拉梁理论的而基础上引入剪切变形影响，因此梁单元仍属于C1型单元。因为并非独立插值，而是由导出。
三维铁木辛柯梁轴向和扭转变形不与其他变形耦合，因此他们的形函数，采用一维线性拉个朗日插值，具体如下
（6）
其中
是无量纲的轴向参数坐标。Xy平面上弯曲变形的形函数通过下述公式推导
（7）
沿单元的剪切应变假定为一常数
。因此弯曲变形引起的转角为
（8）
同样弯矩与曲率的关系为
（9）
由剪切应变引起的剪切力为
（10）
其中k为剪切修正系数，考虑了剪切应力沿着梁截面A分布的不均匀性；E是弹性模量；G是剪切模量，
是沿z轴的转动惯量。根据平衡方程，弯矩与剪切力存在下述关系
（11）
联立公式9和公式10可以得到剪切应变的表达式
（12）
同时形函数还应满足下列边界条件
（13）
写成矩阵形式为
（14）
上式可简写为
d=Aa （15）
通过上式可求得a，并将
和ξ=x/l带入公式7，并经过简化可得
（16）
同理
可以写成
（17）
用同样的方法可以推导处对于xz平面的弯曲变形的形函数
（18）
从上式18可以看出，
（19）
为铁木辛科梁单元任意点的位移（轴线处）；
为铁木辛柯梁单元节点i和j处的位移。
四、单元刚度矩阵
铁木辛柯梁单元的总势能表达式可写为
（20）
将上一节位移的离散表达式带入上述势能方程中，利用最小势能原理可以得到单元局部坐标系下的刚度矩阵的表达式
（21）
和
均布荷载的等效节点力表达式为
（23）
五、局部坐标与整体坐标的转换
上述公式均在局部坐标系下进行推导，在进行整体结构分析时，需要将上述局部坐标系下的矩阵转换到全局坐标系。首先对于全局坐标系下的位移矩阵和局部坐标系下的位移矩阵存在下述转换关系。
全局坐标系下的位移矩阵为
局部坐标系下的位移矩阵为
转换矩阵为
（24）
因此全局坐标系和局部坐标系下的位移、力、刚度矩阵的转换关系为
（25）
接下来讨论如何计算T矩阵，首先确定全局坐标系下的梁轴向向量的分量I,给定梁两个节点的全局坐标(x1,y1,z1,x2,y2,z2)，可以确定T矩阵的第一行元素为
（26）
，由图1所示的节点k来定义。，因此xz平面由向量i和向量k’确定，则xz的法向量
,所以T矩阵的第二行元素为
（27）
图1坐标系中xy平面的法向量
。所以T矩阵第三行元素为
（28）
六. 轴承弹性支座的有限元求解原理
在许多工程应用中，梁都是由弹性构件进行支承的，比如本案例中的车轴通常由滚珠、滚柱或轴颈轴承进行支承。单排滚珠轴承可以看做是：在每个轴承处都有一个节点，并将轴承刚度加到单元刚度矩阵中对应垂直自由度的对角位置上。
滚珠轴承简化为弹性支承梁
对于弹簧支撑的边界条件，可将之看作是混合单元问题，对于本求解模型，看作是梁单元和弹簧单元组合而成的模型。模型本身的边界条件为x节点处受到弹簧的支撑作用和外载荷作用，转换为混合单元后，模型就变为弹簧与地面连接处的自由度为0。
整体刚度矩阵组装分为两步，第一步是组装梁单元，遍历梁单元数量，对单元刚度矩阵按照局部编码和全局编码的原则，组装到整体刚度矩阵中；第二步是组装弹簧单元，遍历弹簧单元的数量，对单元刚度矩阵按照局部编码和全局编码的原则，组装到整体刚度矩阵中，最终的整体刚度矩阵就是模型的整刚。
对于弹簧刚度，这里做一个补充。一维情况下的弹簧刚度矩阵为：
二维情况下，弹簧刚度分为x向刚度和y向刚度，弹簧刚度矩阵为：
对于本次案例，三维铁木辛柯梁，只考虑平动自由度的弹簧刚度，其刚度矩阵为
其中
七、我的Matlab有限元编程精品课
以上就是笔者关注三维铁木辛柯梁单元Matlab有限元编程分析，该内容已经收录在我的原创视频课程里面《Matlab有限元编程从入门到精通30讲》强烈推荐学习者订阅（持续加餐内容）。
可开电子发票，赠送答疑专栏
提供vip群交流，课程可反复回看
识别下方二维码，立即试看
