罗爱玲 景运革

摘要：使用弹性力学解决扭转问题实际上就是求解偏微分方程组，很难得到解析解，有时甚至得不到解析解。为了克服弹性力学解决扭转问题的缺陷，采用有限元法来解决扭转问题，首先以三节点的三角形单元划分网格，并对应力函数进行插值，构造了可用于描述各个子域的场函数。然后利用最小余能原理推导出了扭转问题的泛函，通过求解泛函极值，得到了单元刚度矩阵，最后用Matlab编写了对应的程序用于模拟有限元计算过程。数值算例表明，随着网格的细化，数值解越来越精确，只要网格划分得当，有限元的解能够较好地逼近解析解。

关键词：扭转;应力函数;有限元法;三角形单元

中图分类号：TH133文献标志码：A文章编号：1009-9492（2021）11-0125-04

A Finite Element Method for Calculation of Torsional Rods

Luo Ailing ，Jing Yunge

（Mechanical and Electrical Engineering Department， Yuncheng College， Yuncheng， Shanxi 044000， China）

Abstract： In fact， when the solution of torsional problems was got by elastic mechanics， partial differential equations are needed to solve， which is very difficult to get the analytical solutions. To overcome these shortcomings， the finite element method （FEM） was used to solve the problem. Triangular element with three nodes was used to divide the mesh， and stress function was interpolated to construct the field function which can be used to describe each subdomain. Based on the principle of minimum residual energy， the functional of torsion problem was derived， then the element stiffness was obtained by solving the functional extremum. Finally， corresponding program was written to simulate the calculation process. Examples show that the numerical solutions become more and more accurate with the refinement of the meshes. The solution acquired by FEM can approach the analytical solution well as long as the meshes were divided properly.

Key words： torsion; prandtl function; FEM; triangular element

0 引言

杆结构是工程问题中经常使用的构件。在航天领域，飞船太阳能帆板的可展折骨架采用铰链连杆机构;在机械工程领域，起重机、塔吊使用梁和钢桁架来增加其刚度和强度，汽车传动轴、机床的光杆及测量仪器的传动装置都是采用杆结构来传递运动和力。杆结构已经渗透到生产及生活中的各个领域，是一种十分重要的构件。

杆件的受力形式主要有拉、壓、剪切、弯曲、扭转以及弯扭组合，其中扭转问题普遍存在于工程实际中。对于对动力要求不大的汽车多采用圆截面轴，若需要传递更大的扭矩则需要采用椭圆截面或矩形截面轴，轴的抗扭刚度和内部的应力将直接决定汽车的安全性能。井钻的钻杆在工作时也要承受巨大的扭矩，有时钻井深度可达几千米，钻杆一旦断裂，钻头很可能无法取出，将造成巨大的经济损失。车床上的光杆和丝杆也要承受扭转载荷，其抗扭刚度将间接影响到车床的切削精度。很多桥梁尤其是跨度比较大的悬索桥，在风速较大时容易因涡振而发生扭转，倘若桥身抗扭刚度不够，很可能造成桥身断裂，桥毁人亡。因此，为了计算其内部的应力应变、抗扭刚度、极限载荷、破坏条件和寿命，预防工程事故的发生，对扭转杆进行力学分析是十分必要的。

近几年，国内外学者对扭转问题的研究取得了巨大的进步，黄宏等[1]研究了受轴力和扭矩荷载作用的构件并分析了构件的受力性能。董云霞等[2]分析了箱梁的扭转问题，对于扭转单元，形函数采用三阶多项式，根据能量法和变分原理得到了单元刚度矩阵和等效结点载荷。朱鑫等[3]对行李箱盖铰链杆的刚度进行了有限元分析，并进行刚度试验，对行李箱盖铰链杆的结构进行了优化。杨浩等[4]推导了结构整体刚度矩阵与单元刚度矩阵各元素之间的关系用于计算杆梁组合结构，并建立了基于关联表的结构整体分析模型。许晶等[5]以Vlasov扭转理论为基础，构建了压扭杆的位移场函数，建立了与控制方程等效的泛函。利用求泛函极值得出了解析型压扭杆单元列式，并推导了用于杆件内力分析的单元刚度矩阵[5]。

另外，有限元法是一种重要的数值分析方法，在各种领域都表现出了独特的优势，已被广大学者所普遍接受。杨权等[6]采用工程测试的方法对倾动机构进行倾动力矩的测量，将扭力杆装置的作用力作为边界条件，通过有限元分析的方法计算了扭力杆装置的应力。在 Bin 等[7]的研究中，三维有限元法被用以分析发动机连杆的各项参数。Zhang等[8]采用有限元法分析了电机的输入和输出性能。为了验证该方法的有效性，将其应用于棒形超声电机的优化设计中，设计结果表明优化后与电机性能相关的设计指标有明显改善。Orwoll 等[9]用基于临床 QCT影像的有限元模型计算髋关节强度，并研究了该方法预测老年人髋部骨折的效果，证明了 QCT扫描的有限元法生物力学分析对预测男性髋部骨折是有效的。石晓燕等[10]依据上海某深基坑工程，基于杆系有限元法，分别讨论了每道支撑位置、挡土结构厚度及其入土深度对板式挡土结构内力、位移的影响。刘健[11]对位置有限元的基本理论与方法进行系统地研究，以工程中常用的梁单元为研究对象，构造了新型的三结点等截面梁单元和变截面梁单元，用于结构的静态及动态响应分析。综合国内外学者对扭转问题的研究来看，抗扭刚度和剪应力是扭转杆的两个重要参数，如何便捷高效的求解是个值得深入研究的问题。因此本文采用有限元法对柱体扭转问题进行分析。首先，以三节点的三角形单元对截面进行网格划分，之后对应力函数进行插值，构造出各个子域上的应力场函数。利用最小余能原理推导出与扭转问题控制方程等价的泛函，通过对泛函求极值，得到单元刚度矩阵和单元等效结点载荷，将单元刚度矩阵按位置组装程总体刚度矩阵，单元等效结点载荷组装成总体节点载荷向量，最后用Matlab编写对应的程序用于模拟有限元计算过程。

1 扭转问题基本方程

扭转问题本是一种三维问题，但本文只考虑自由扭

转，即两端无轴向载荷。此时截面上的应力、应变以及位移分布情况完全相同，只需分析任意一个截面即可，这样，就把三维问题简化为二维问题。由材料力学假设可知，截面上除了切应力其他量均为0，因此需要以应力表示其他物理量。

1.1 应力表示的基本方程

对于一等截面直杆，两端作用大小相等方向相反的扭矩，两端无约束，属于自由扭转。由弹性力学可知，应力的物理方程可描述为：

式中：τxz和τzy是截面上的切应力;σx 、σy 、σz 和τxy分别为不同方向的应力分量，对于非圆截面的扭转问题，则不同方向的应力分量值为0，σx =σy =σz =0，τxy =0。

根据Θ=σx +σy +σz ，扭转相容方程为：

由于材料在变形过程中是连续的，不应出现“重叠”或“撕裂”，应力解法还必须满足用应力表示的变形协调方程：

式中： l 、m 分别为侧面单位法向量的方向余弦。

1.2 应力函数表示的基本方程

假设应力函数为φ0时，单位长度扭转角α=1，则实际应力函数描述为：

式中：函数φ（x，y）为普朗特应力函数。

2 扭转问题的有限元分析

本节给出有限元法的相关概念，并给出有限元法的流程。

2.1 形函数

有限元分析最重要的一部就是离散，选用何种单元进行离散影响着最终求解精度。常见的平面单元有三节点三角形单元、六节点三角形单元、四节点矩形单元、八节点矩形单元以及多边形单元等。就单个单元来说，三节点三角形单元逼近效果较差，但其对边界形狀的适应性较好，能够减少离散误差，而且随着网格的细化，计算误差将越来小，只要设置合适的网格密度，即可得到比较精确的结果。因此，本文采用三节点的三角形单元对求解域进行离散。三角形单元如图1所示。

i、j、k 分别为三角形单元的3个顶点，依次按逆时针顺序排列， （xi ，yi）、 （xj ，yj）、 （xm ，ym）分别为3个节点的坐标，φi、φj 、φm 分别为3个节点已知的应力函数值，假设应力函数的形式如下：

φ=a1+a2 x +a3y （5）

2.2 结构离散化

针对具体截面形状，需要根据其应力分布情况合理布置单元的数量和疏密程度，从几何上来说就是用一张连接各节点的网格代替截面。单元和单元之间靠节点连接，因此各个单元的物理量也是通过节点来传递，这就需要将节点号、节点坐标和单元号联系起来，形成节点索引矩阵。单元数量不多或自然离散时可以手动划分网格，然后对单元和节点进行编号并计算出各个节点的坐标。但当边界较复杂，需要较多的单元来逼近几何形状时，使用程序来划分网格较为妥当，这样可以很方便地对网格进行细化处理，而且能保证单元信息的输入准确无误。

2.3 单元分析

得到单元信息后，可得扭转问题的控制方程和边界条件描述如下：

其等价的泛函描述为：

2.4 单元分析

采用文献[1]中直接刚度法进行单元组装，组装过程如下：首先建立一个行数和列数均为总节点数的方阵，其行号和列号均就是总体刚度矩阵中的节点号，各单元刚度矩阵中的元素按其对应的节点号在总体刚度矩阵中对号入座，同位置元素累加。同理，单元等效载荷中的元素按其对应的节点号在总体载荷向量中对号入座，同位置元素累加。

式中：K 为总体刚度矩阵;f为总体载荷阵列。

2.5 施加边界条件并求解

根据已经得到了所有节点的平衡方程，但还不能联立解出节点应力函数，因为此时的总体刚度矩阵为奇异矩阵，方程没有唯一解，只有引入约束后才能解出唯一的节点应力函数。本文采用直接缩减法进行边界条件的引入。总体平衡方程划分为：

2.6 有限元法的流程图

有限元法的流程如图2所示。

3 扭转问题仿真实验分析

为了验证有限元法计算扭转问题的有效性，在计算机上对扭转实际问题进行了模拟仿真实验，实验所用的平台是 Windows 7操作系统，并利用Matlab 7.0软件编写程序代码对具体算例进行计算，并将计算结果与解析解作对比分析。实验所用的实例如例1所示。

例1：假设有一截面为等边三角形的杆，边长为 a =2 m ，截面端部受扭矩 M =6×103 N ? m 的作用，剪切弹性模量为 G =80 GPa，求截面上的切应力和抗扭刚度。

根据相关文献可知截面的应力函数解析解计算公式为：

通过有限元法对上述实例进行模拟实验，实验过程如下。

（1） 网格划分如图3所示。

（2） 得到的节点应力函数为离散点，对 x 轴上对应点的应力函数值进行插值，并与解析法进行对比，结果如图4所示。从图4可知有限元解和解析法的精度基本一致，说明有限元解是有效的。

（3） 由式（4） 可知，对应力函数求偏导即可得到切应力，对应力函数求导，得到 x 轴上各点的切应力τyz，结果如图5所示。从图5可知抗扭刚度有限元解 D =8.322 N ? m-1，解析解 D =8.327 N ? m-1。

4 结束语

本文对有限元方法计算扭转杆截面上的切应力以及杆件的抗扭刚度的优点进行总结，并将对扭转杆复杂计算进行深入研究。

（1） 应力函数的使用很大程度上减小了求解偏微分方程组的难度，尤其是在划分网格时减少了总自由度的数量，使得总体刚度矩阵的维数减小了一半，提高了计算效率。

（2） 相比于翘曲函数法，应力函数法的边界条件更为简单，有利于约束条件的引入。

（3） 对于抗扭刚度的求解，直接令单位扭转角为1，求出应力函数后直接对其进行积分，此方法比翘曲函数法的积分公式简单一些，结果也非常精确。

（4） 本文提出的有限元法仅仅是对扭转杆中一个简单的截面进行分析，但是在实际应用中截面往往是复杂的，因此，如何提高算法的计算效率是将来继续深入研究的问题。

参考文献：

[1]黄宏， 郭晓宇， 陈梦成.圆中空夹层钢管混凝土压扭构件有限元分析[J].建筑结构学报， 2013， 34（1）：50-56.

[2]董云霞.箱梁桥扭转与畸变应力分析及程序设计[D].湘潭：湖南科技大学，2011.

[3]朱鑫，寒冬桂， 刘芳，等.乘用车行李箱盖铰链杆的刚度有限元分析與结构优化[J].机械工程，2019，57（2）：21-25.

[4]杨浩， 罗帅，邢国然，等.杆梁组合结构的有限元分析[J].工程力学，2019，36（1）：154-157.

[5]许晶，夏文忠，王宏志，等.考虑大位移影响的解析型压扭杆单元[J].工程力学，2019，36（4）：44-51.

[6]杨权，李小标，何康.基于工程测试的托圈扭力杆装置有限元应力分析[J].平顶山学院学报，2019，33（5）：40-52.

[7] Bin Z ，Lixia J ， Yongqi L . Finite Element Analysis and Structur- al Improvement of Diesel Engine Connecting Rod[C]//2010 Sec- ond International Conference on Computer Modeling and Simula- tion. IEEE， 2010.

[8] Zhang J T ， Zhu H ， Zhou S Q ， et al. Optimal design of a rod shape ultrasonic motor using sequential quadratic programming and finite element method[J]. Finite Elements in Analysis and Design， 2012， 59（2）：11-17.

[9] Orwoll E S ， Marshall L M ， Nielson C M ， et al. Finite Element Analysis of the Proximal Femur and Hip Fracture Risk in Older Men[J]. Journal of bone and mineral research： the official journal of the American Society for Bone and Mineral Research， 2009， 24（3）：475-483.

[10]石晓燕，贺永胜，王启睿，等.基于干系有限元法板式挡土结构内力位移的影响分析[J].水利与建筑工程学报，2019，17（6）：188-193.

[11]刘健.弹性梁几何非线性问题的位置有限元法及其应用研究[D].济南：山东大学，2012.

作者简介：

罗爱玲（1970-），女，山西永济人，大学本科，实验员，研究领域为智能计算。

景运革（1970-），男，博士，教授，研究领域为粗糙集理论与粒计算。 （编辑：刁少华）