有限元法在扭轉桿計算中的應用

羅愛玲 景運革

摘要：使用彈性力學解決扭轉問題實際上就是求解偏微分方程組，很難得到解析解，有時甚至得不到解析解。為了克服彈性力學解決扭轉問題的缺陷，采用有限元法來解決扭轉問題，首先以三節點的三角形單元劃分網格，并對應力函數進行插值，構造了可用于描述各個子域的場函數。然后利用最小余能原理推導出了扭轉問題的泛函，通過求解泛函極值，得到了單元剛度矩陣，最后用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-），男，博士，教授，研究領域為粗糙集理論與粒計算。 （編輯：刁少華）