有限元法在扭轉(zhuǎn)桿計(jì)算中的應(yīng)用

羅愛玲 景運(yùn)革

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

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

中圖分類號(hào)：TH133文獻(xiàn)標(biāo)志碼：A文章編號(hào)：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)域，飛船太陽(yáng)能帆板的可展折骨架采用鉸鏈連桿機(jī)構(gòu);在機(jī)械工程領(lǐng)域，起重機(jī)、塔吊使用梁和鋼桁架來(lái)增加其剛度和強(qiáng)度，汽車傳動(dòng)軸、機(jī)床的光桿及測(cè)量?jī)x器的傳動(dòng)裝置都是采用桿結(jié)構(gòu)來(lái)傳遞運(yùn)動(dòng)和力。桿結(jié)構(gòu)已經(jīng)滲透到生產(chǎn)及生活中的各個(gè)領(lǐng)域，是一種十分重要的構(gòu)件。

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

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

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

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

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

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

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

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

式中：τxz和τzy是截面上的切應(yīng)力;σx 、σy 、σz 和τxy分別為不同方向的應(yīng)力分量，對(duì)于非圓截面的扭轉(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時(shí)，單位長(zhǎng)度扭轉(zhuǎn)角α=1，則實(shí)際應(yīng)力函數(shù)描述為：

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

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

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

2.1 形函數(shù)

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

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

φ=a1+a2 x +a3y （5）

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

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

2.3 單元分析

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

其等價(jià)的泛函描述為：

2.4 單元分析

采用文獻(xiàn)[1]中直接剛度法進(jìn)行單元組裝，組裝過程如下：首先建立一個(gè)行數(shù)和列數(shù)均為總節(jié)點(diǎn)數(shù)的方陣，其行號(hào)和列號(hào)均就是總體剛度矩陣中的節(jié)點(diǎn)號(hào)，各單元?jiǎng)偠染仃囍械脑匕雌鋵?duì)應(yīng)的節(jié)點(diǎn)號(hào)在總體剛度矩陣中對(duì)號(hào)入座，同位置元素累加。同理，單元等效載荷中的元素按其對(duì)應(yīng)的節(jié)點(diǎn)號(hào)在總體載荷向量中對(duì)號(hào)入座，同位置元素累加。

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

2.5 施加邊界條件并求解

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

2.6 有限元法的流程圖

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

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

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

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

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

通過有限元法對(duì)上述實(shí)例進(jìn)行模擬實(shí)驗(yàn)，實(shí)驗(yàn)過程如下。

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

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

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

4 結(jié)束語(yǔ)

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

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

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

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

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

參考文獻(xiàn)：

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

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

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

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

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

[6]楊權(quán)，李小標(biāo)，何康.基于工程測(cè)試的托圈扭力桿裝置有限元應(yīng)力分析[J].平頂山學(xué)院學(xué)報(bào)，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é)報(bào)，2019，17（6）：188-193.

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

作者簡(jiǎn)介：

羅愛玲（1970-），女，山西永濟(jì)人，大學(xué)本科，實(shí)驗(yàn)員，研究領(lǐng)域?yàn)橹悄苡?jì)算。

景運(yùn)革（1970-），男，博士，教授，研究領(lǐng)域?yàn)榇植诩碚撆c粒計(jì)算。 （編輯：刁少華）