羅愛玲 景運革

摘要：使用彈性力學(xué)解決扭轉(zhuǎn)問題實際上就是求解偏微分方程組，很難得到解析解，有時甚至得不到解析解。為了克服彈性力學(xué)解決扭轉(zhuǎn)問題的缺陷，采用有限元法來解決扭轉(zhuǎn)問題，首先以三節(jié)點的三角形單元劃分網(wǎng)格，并對應(yīng)力函數(shù)進行插值，構(gòu)造了可用于描述各個子域的場函數(shù)。然后利用最小余能原理推導(dǎo)出了扭轉(zhuǎn)問題的泛函，通過求解泛函極值，得到了單元剛度矩陣，最后用Matlab編寫了對應(yīng)的程序用于模擬有限元計算過程。數(shù)值算例表明，隨著網(wǎng)格的細(xì)化，數(shù)值解越來越精確，只要網(wǎng)格劃分得當(dāng)，有限元的解能夠較好地逼近解析解。

關(guān)鍵詞：扭轉(zhuǎn);應(yīng)力函數(shù);有限元法;三角形單元

中圖分類號：TH133文獻標(biāo)志碼：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 引言

桿結(jié)構(gòu)是工程問題中經(jīng)常使用的構(gòu)件。在航天領(lǐng)域，飛船太陽能帆板的可展折骨架采用鉸鏈連桿機構(gòu);在機械工程領(lǐng)域，起重機、塔吊使用梁和鋼桁架來增加其剛度和強度，汽車傳動軸、機床的光桿及測量儀器的傳動裝置都是采用桿結(jié)構(gòu)來傳遞運動和力。桿結(jié)構(gòu)已經(jīng)滲透到生產(chǎn)及生活中的各個領(lǐng)域，是一種十分重要的構(gòu)件。

桿件的受力形式主要有拉、壓、剪切、彎曲、扭轉(zhuǎn)以及彎扭組合，其中扭轉(zhuǎn)問題普遍存在于工程實際中。對于對動力要求不大的汽車多采用圓截面軸，若需要傳遞更大的扭矩則需要采用橢圓截面或矩形截面軸，軸的抗扭剛度和內(nèi)部的應(yīng)力將直接決定汽車的安全性能。井鉆的鉆桿在工作時也要承受巨大的扭矩，有時鉆井深度可達幾千米，鉆桿一旦斷裂，鉆頭很可能無法取出，將造成巨大的經(jīng)濟損失。車床上的光桿和絲桿也要承受扭轉(zhuǎn)載荷，其抗扭剛度將間接影響到車床的切削精度。很多橋梁尤其是跨度比較大的懸索橋，在風(fēng)速較大時容易因渦振而發(fā)生扭轉(zhuǎn)，倘若橋身抗扭剛度不夠，很可能造成橋身斷裂，橋毀人亡。因此，為了計算其內(nèi)部的應(yīng)力應(yīng)變、抗扭剛度、極限載荷、破壞條件和壽命，預(yù)防工程事故的發(fā)生，對扭轉(zhuǎn)桿進行力學(xué)分析是十分必要的。

近幾年，國內(nèi)外學(xué)者對扭轉(zhuǎn)問題的研究取得了巨大的進步，黃宏等[1]研究了受軸力和扭矩荷載作用的構(gòu)件并分析了構(gòu)件的受力性能。董云霞等[2]分析了箱梁的扭轉(zhuǎn)問題，對于扭轉(zhuǎn)單元，形函數(shù)采用三階多項式，根據(jù)能量法和變分原理得到了單元剛度矩陣和等效結(jié)點載荷。朱鑫等[3]對行李箱蓋鉸鏈桿的剛度進行了有限元分析，并進行剛度試驗，對行李箱蓋鉸鏈桿的結(jié)構(gòu)進行了優(yōu)化。楊浩等[4]推導(dǎo)了結(jié)構(gòu)整體剛度矩陣與單元剛度矩陣各元素之間的關(guān)系用于計算桿梁組合結(jié)構(gòu)，并建立了基于關(guān)聯(lián)表的結(jié)構(gòu)整體分析模型。許晶等[5]以Vlasov扭轉(zhuǎn)理論為基礎(chǔ)，構(gòu)建了壓扭桿的位移場函數(shù)，建立了與控制方程等效的泛函。利用求泛函極值得出了解析型壓扭桿單元列式，并推導(dǎo)了用于桿件內(nèi)力分析的單元剛度矩陣[5]。

另外，有限元法是一種重要的數(shù)值分析方法，在各種領(lǐng)域都表現(xiàn)出了獨特的優(yōu)勢，已被廣大學(xué)者所普遍接受。楊權(quán)等[6]采用工程測試的方法對傾動機構(gòu)進行傾動力矩的測量，將扭力桿裝置的作用力作為邊界條件，通過有限元分析的方法計算了扭力桿裝置的應(yīng)力。在 Bin 等[7]的研究中，三維有限元法被用以分析發(fā)動機連桿的各項參數(shù)。Zhang等[8]采用有限元法分析了電機的輸入和輸出性能。為了驗證該方法的有效性，將其應(yīng)用于棒形超聲電機的優(yōu)化設(shè)計中，設(shè)計結(jié)果表明優(yōu)化后與電機性能相關(guān)的設(shè)計指標(biāo)有明顯改善。Orwoll 等[9]用基于臨床 QCT影像的有限元模型計算髖關(guān)節(jié)強度，并研究了該方法預(yù)測老年人髖部骨折的效果，證明了 QCT掃描的有限元法生物力學(xué)分析對預(yù)測男性髖部骨折是有效的。石曉燕等[10]依據(jù)上海某深基坑工程，基于桿系有限元法，分別討論了每道支撐位置、擋土結(jié)構(gòu)厚度及其入土深度對板式擋土結(jié)構(gòu)內(nèi)力、位移的影響。劉健[11]對位置有限元的基本理論與方法進行系統(tǒng)地研究，以工程中常用的梁單元為研究對象，構(gòu)造了新型的三結(jié)點等截面梁單元和變截面梁單元，用于結(jié)構(gòu)的靜態(tài)及動態(tài)響應(yīng)分析。綜合國內(nèi)外學(xué)者對扭轉(zhuǎn)問題的研究來看，抗扭剛度和剪應(yīng)力是扭轉(zhuǎn)桿的兩個重要參數(shù)，如何便捷高效的求解是個值得深入研究的問題。因此本文采用有限元法對柱體扭轉(zhuǎn)問題進行分析。首先，以三節(jié)點的三角形單元對截面進行網(wǎng)格劃分，之后對應(yīng)力函數(shù)進行插值，構(gòu)造出各個子域上的應(yīng)力場函數(shù)。利用最小余能原理推導(dǎo)出與扭轉(zhuǎn)問題控制方程等價的泛函，通過對泛函求極值，得到單元剛度矩陣和單元等效結(jié)點載荷，將單元剛度矩陣按位置組裝程總體剛度矩陣，單元等效結(jié)點載荷組裝成總體節(jié)點載荷向量，最后用Matlab編寫對應(yīng)的程序用于模擬有限元計算過程。

1 扭轉(zhuǎn)問題基本方程

扭轉(zhuǎn)問題本是一種三維問題，但本文只考慮自由扭

轉(zhuǎn)，即兩端無軸向載荷。此時截面上的應(yīng)力、應(yīng)變以及位移分布情況完全相同，只需分析任意一個截面即可，這樣，就把三維問題簡化為二維問題。由材料力學(xué)假設(shè)可知，截面上除了切應(yīng)力其他量均為0，因此需要以應(yīng)力表示其他物理量。

1.1 應(yīng)力表示的基本方程

對于一等截面直桿，兩端作用大小相等方向相反的扭矩，兩端無約束，屬于自由扭轉(zhuǎn)。由彈性力學(xué)可知，應(yīng)力的物理方程可描述為：

式中：τxz和τzy是截面上的切應(yīng)力;σx 、σy 、σz 和τxy分別為不同方向的應(yīng)力分量，對于非圓截面的扭轉(zhuǎn)問題，則不同方向的應(yīng)力分量值為0，σx =σy =σz =0，τxy =0。

根據(jù)Θ=σx +σy +σz ，扭轉(zhuǎn)相容方程為：

由于材料在變形過程中是連續(xù)的，不應(yīng)出現(xiàn)“重疊”或“撕裂”，應(yīng)力解法還必須滿足用應(yīng)力表示的變形協(xié)調(diào)方程：

式中： l 、m 分別為側(cè)面單位法向量的方向余弦。

1.2 應(yīng)力函數(shù)表示的基本方程

假設(shè)應(yīng)力函數(shù)為φ0時，單位長度扭轉(zhuǎn)角α=1，則實際應(yīng)力函數(shù)描述為：

式中：函數(shù)φ（x，y）為普朗特應(yīng)力函數(shù)。

2 扭轉(zhuǎn)問題的有限元分析

本節(jié)給出有限元法的相關(guān)概念，并給出有限元法的流程。

2.1 形函數(shù)

有限元分析最重要的一部就是離散，選用何種單元進行離散影響著最終求解精度。常見的平面單元有三節(jié)點三角形單元、六節(jié)點三角形單元、四節(jié)點矩形單元、八節(jié)點矩形單元以及多邊形單元等。就單個單元來說，三節(jié)點三角形單元逼近效果較差，但其對邊界形狀的適應(yīng)性較好，能夠減少離散誤差，而且隨著網(wǎng)格的細(xì)化，計算誤差將越來小，只要設(shè)置合適的網(wǎng)格密度，即可得到比較精確的結(jié)果。因此，本文采用三節(jié)點的三角形單元對求解域進行離散。三角形單元如圖1所示。

i、j、k 分別為三角形單元的3個頂點，依次按逆時針順序排列， （xi ，yi）、 （xj ，yj）、 （xm ，ym）分別為3個節(jié)點的坐標(biāo)，φi、φj 、φm 分別為3個節(jié)點已知的應(yīng)力函數(shù)值，假設(shè)應(yīng)力函數(shù)的形式如下：

φ=a1+a2 x +a3y （5）

2.2 結(jié)構(gòu)離散化

針對具體截面形狀，需要根據(jù)其應(yīng)力分布情況合理布置單元的數(shù)量和疏密程度，從幾何上來說就是用一張連接各節(jié)點的網(wǎng)格代替截面。單元和單元之間靠節(jié)點連接，因此各個單元的物理量也是通過節(jié)點來傳遞，這就需要將節(jié)點號、節(jié)點坐標(biāo)和單元號聯(lián)系起來，形成節(jié)點索引矩陣。單元數(shù)量不多或自然離散時可以手動劃分網(wǎng)格，然后對單元和節(jié)點進行編號并計算出各個節(jié)點的坐標(biāo)。但當(dāng)邊界較復(fù)雜，需要較多的單元來逼近幾何形狀時，使用程序來劃分網(wǎng)格較為妥當(dāng)，這樣可以很方便地對網(wǎng)格進行細(xì)化處理，而且能保證單元信息的輸入準(zhǔn)確無誤。

2.3 單元分析

得到單元信息后，可得扭轉(zhuǎn)問題的控制方程和邊界條件描述如下：

其等價的泛函描述為：

2.4 單元分析

采用文獻[1]中直接剛度法進行單元組裝，組裝過程如下：首先建立一個行數(shù)和列數(shù)均為總節(jié)點數(shù)的方陣，其行號和列號均就是總體剛度矩陣中的節(jié)點號，各單元剛度矩陣中的元素按其對應(yīng)的節(jié)點號在總體剛度矩陣中對號入座，同位置元素累加。同理，單元等效載荷中的元素按其對應(yīng)的節(jié)點號在總體載荷向量中對號入座，同位置元素累加。

式中：K 為總體剛度矩陣;f為總體載荷陣列。

2.5 施加邊界條件并求解

根據(jù)已經(jīng)得到了所有節(jié)點的平衡方程，但還不能聯(lián)立解出節(jié)點應(yīng)力函數(shù)，因為此時的總體剛度矩陣為奇異矩陣，方程沒有唯一解，只有引入約束后才能解出唯一的節(jié)點應(yīng)力函數(shù)。本文采用直接縮減法進行邊界條件的引入?？傮w平衡方程劃分為：

2.6 有限元法的流程圖

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

3 扭轉(zhuǎn)問題仿真實驗分析

為了驗證有限元法計算扭轉(zhuǎn)問題的有效性，在計算機上對扭轉(zhuǎn)實際問題進行了模擬仿真實驗，實驗所用的平臺是 Windows 7操作系統(tǒng)，并利用Matlab 7.0軟件編寫程序代碼對具體算例進行計算，并將計算結(jié)果與解析解作對比分析。實驗所用的實例如例1所示。

例1：假設(shè)有一截面為等邊三角形的桿，邊長為 a =2 m ，截面端部受扭矩 M =6×103 N ? m 的作用，剪切彈性模量為 G =80 GPa，求截面上的切應(yīng)力和抗扭剛度。

根據(jù)相關(guān)文獻可知截面的應(yīng)力函數(shù)解析解計算公式為：

通過有限元法對上述實例進行模擬實驗，實驗過程如下。

（1） 網(wǎng)格劃分如圖3所示。

（2） 得到的節(jié)點應(yīng)力函數(shù)為離散點，對 x 軸上對應(yīng)點的應(yīng)力函數(shù)值進行插值，并與解析法進行對比，結(jié)果如圖4所示。從圖4可知有限元解和解析法的精度基本一致，說明有限元解是有效的。

（3） 由式（4） 可知，對應(yīng)力函數(shù)求偏導(dǎo)即可得到切應(yīng)力，對應(yīng)力函數(shù)求導(dǎo)，得到 x 軸上各點的切應(yīng)力τyz，結(jié)果如圖5所示。從圖5可知抗扭剛度有限元解 D =8.322 N ? m-1，解析解 D =8.327 N ? m-1。

4 結(jié)束語

本文對有限元方法計算扭轉(zhuǎn)桿截面上的切應(yīng)力以及桿件的抗扭剛度的優(yōu)點進行總結(jié)，并將對扭轉(zhuǎn)桿復(fù)雜計算進行深入研究。

（1） 應(yīng)力函數(shù)的使用很大程度上減小了求解偏微分方程組的難度，尤其是在劃分網(wǎng)格時減少了總自由度的數(shù)量，使得總體剛度矩陣的維數(shù)減小了一半，提高了計算效率。

（2） 相比于翹曲函數(shù)法，應(yīng)力函數(shù)法的邊界條件更為簡單，有利于約束條件的引入。

（3） 對于抗扭剛度的求解，直接令單位扭轉(zhuǎn)角為1，求出應(yīng)力函數(shù)后直接對其進行積分，此方法比翹曲函數(shù)法的積分公式簡單一些，結(jié)果也非常精確。

（4） 本文提出的有限元法僅僅是對扭轉(zhuǎn)桿中一個簡單的截面進行分析，但是在實際應(yīng)用中截面往往是復(fù)雜的，因此，如何提高算法的計算效率是將來繼續(xù)深入研究的問題。

參考文獻：

[1]黃宏， 郭曉宇， 陳夢成.圓中空夾層鋼管混凝土壓扭構(gòu)件有限元分析[J].建筑結(jié)構(gòu)學(xué)報， 2013， 34（1）：50-56.

[2]董云霞.箱梁橋扭轉(zhuǎn)與畸變應(yīng)力分析及程序設(shè)計[D].湘潭：湖南科技大學(xué)，2011.

[3]朱鑫，寒冬桂， 劉芳，等.乘用車行李箱蓋鉸鏈桿的剛度有限元分析與結(jié)構(gòu)優(yōu)化[J].機械工程，2019，57（2）：21-25.

[4]楊浩， 羅帥，邢國然，等.桿梁組合結(jié)構(gòu)的有限元分析[J].工程力學(xué)，2019，36（1）：154-157.

[5]許晶，夏文忠，王宏志，等.考慮大位移影響的解析型壓扭桿單元[J].工程力學(xué)，2019，36（4）：44-51.

[6]楊權(quán)，李小標(biāo)，何康.基于工程測試的托圈扭力桿裝置有限元應(yīng)力分析[J].平頂山學(xué)院學(xué)報，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]石曉燕，賀永勝，王啟睿，等.基于干系有限元法板式擋土結(jié)構(gòu)內(nèi)力位移的影響分析[J].水利與建筑工程學(xué)報，2019，17（6）：188-193.

[11]劉健.彈性梁幾何非線性問題的位置有限元法及其應(yīng)用研究[D].濟南：山東大學(xué)，2012.

作者簡介：

羅愛玲（1970-），女，山西永濟人，大學(xué)本科，實驗員，研究領(lǐng)域為智能計算。

景運革（1970-），男，博士，教授，研究領(lǐng)域為粗糙集理論與粒計算。 （編輯：刁少華）