趙蘭磊，陳文家，何彥忠，李鷺揚＊，黃 鑫

（1.揚州大學(xué)機(jī)械工程學(xué)院，江蘇揚州225127；2.江蘇揚力集團(tuán)精密機(jī)床研究所，江蘇揚州225127）

復(fù)雜截面扭轉(zhuǎn)問題的孔洞擬填充法分析

趙蘭磊1，陳文家1，何彥忠2，李鷺揚1＊，黃 鑫1

（1.揚州大學(xué)機(jī)械工程學(xué)院，江蘇揚州225127；2.江蘇揚力集團(tuán)精密機(jī)床研究所，江蘇揚州225127）

為解決復(fù)雜截面應(yīng)力、變形問題的求解，以柱體扭轉(zhuǎn)問題的應(yīng)力函數(shù)及多連域轉(zhuǎn)化為單連域的方法為基礎(chǔ)，應(yīng)用有限元－邊界元法，提出了“孔洞擬填充”的方法，并利用有限元分析軟件ANSYS Workbench以偏心輪為例模擬求解扭轉(zhuǎn)問題.結(jié)果證明該方法可以適用于任何截面的桿件.

復(fù)雜截面；應(yīng)力函數(shù)；多連域；擬填充；有限元；扭轉(zhuǎn)

在車輛以及機(jī)械設(shè)備中普遍使用承受扭矩載荷的部件（如偏心輪，偏心齒輪）及采用非圓截面（如矩形、三角形、橢圓形等）的軸，這些承扭件的橫截面往往為復(fù)雜的多連域，對其扭轉(zhuǎn)問題通常難以求得解析解，從而使其強(qiáng)度和剛度計算成為設(shè)計中的難點.目前處理復(fù)雜截面的扭轉(zhuǎn)問題所采用的方法有變分解法、邊界元法［1－3］、虛邊界元法［4］、分離變量法［5］、復(fù)變函數(shù)法［6］、多連域轉(zhuǎn)化為單連域法［7］、有限元法［8－9］等.如何借助計算機(jī)和現(xiàn)代分析手段分析承扭桿的結(jié)構(gòu)問題，并精確求解這些扭桿件的應(yīng)力、應(yīng)變狀態(tài)及抗扭剛度，具有一定的研究價值.本文以有限元－邊界元法及多連域轉(zhuǎn)化為單連域的方法為基礎(chǔ)，基于虛單元、截面的拓?fù)鋬?yōu)化［10］、輪廓的孔洞填充算法［11］、無網(wǎng)格有限元法［8］、間斷有限元法［9］、任意復(fù)雜薄壁截面自由扭轉(zhuǎn)常數(shù)的計算［12］、多單元扭轉(zhuǎn)常數(shù)的分析［13］、任意截面形狀的鋼筋混凝土桿的扭轉(zhuǎn)剛度［14］，分析多單體梁的扭轉(zhuǎn)、變形和翹曲［15］，提出了“孔洞擬填充”的概念，應(yīng)用“孔洞擬填充”思想來處理復(fù)雜截面的扭轉(zhuǎn)問題.

圖1 任意多連域截面Fig.1 Arbitrary multipleconnected cross sections

1 復(fù)雜截面扭轉(zhuǎn)問題的理論以及實例分析

1.1 有限元－邊界元法

根據(jù)彈性力學(xué)中在柱體周圍側(cè)面（曲面）上的應(yīng)力關(guān)系及柱體扭轉(zhuǎn)問題的性質(zhì)，可以得出應(yīng)力在邊界上的關(guān)系式為τxzcos（n，x）＋τyzcos（n，y）＝0（n是側(cè)面的外法線方向），也可以寫成

由式（1）可得知應(yīng)力函數(shù)在邊界上一定是一個常數(shù).對于一般情況下的多連域截面扭轉(zhuǎn)問題，其圖形如圖1所示，則多連域截面的扭矩可以寫為

式中G為剪切彈性模量；Φ為應(yīng)力函數(shù)；θ為單位長度的扭轉(zhuǎn)角.

對式（2）進(jìn)行積分計算，并運用格林公式將面積的積分轉(zhuǎn)化為沿曲線的積分，可以簡化為

由式（1）已得知應(yīng)力函數(shù)在邊界上是一個常數(shù)，由于通常情況下取截面最外邊界Γ0（截面外邊界）上的應(yīng)力函數(shù)Φ0＝0，而其他的內(nèi)邊界Γ1，…，Γn在每一邊界上應(yīng)力函數(shù)Φ都是常數(shù)，但這個常數(shù)不一定相等，故取此常數(shù)為k（待定的常數(shù)），即Φi＝ki，其中i＝1，2，3，…，n.假設(shè)S與截面內(nèi)一個邊界Γi重合，所以式（3）可以變?yōu)?/p>

按照應(yīng)力函數(shù)求解多連域截面桿扭轉(zhuǎn)問題的邊界元法，可得邊值問題的關(guān)系式

因為多連域截面內(nèi)的應(yīng)力分量為τxy，τyz，因此對應(yīng)的普朗特應(yīng)力函數(shù)為

由式（9）可以得到某一點處的應(yīng)力τ為

由此可見，用邊界元法求解多連域截面桿扭轉(zhuǎn)問題，須將多連域截面的邊值問題轉(zhuǎn)化為在相同區(qū)域內(nèi)同一離散格式下進(jìn)行求解.首先解出Φ和θ，然后代入式（8），再與式（9）、（10）聯(lián)合求解，即可確定多連域截面桿在扭矩作用下的應(yīng)力問題.

依據(jù)上述方法，以偏心輪為例，借助于有限元－邊界元法對其進(jìn)行分析計算.已知：R1，R2，R3，R4；截面承受的扭矩為M；各點對應(yīng)的位置坐標(biāo)見圖2.根據(jù)圖2及其相關(guān)參數(shù)確定各坐標(biāo)點之間的關(guān)系式和矩陣關(guān)系式，然后確定孔洞的面積，利用公式Ai＝Ωjdxdy對孔洞部分進(jìn)行面積的求解，將所求面積及公式Φ＝代入式（7），可以算出待定常數(shù)ki，再回代到公式Φ＝，即可求得應(yīng)力函數(shù)Φ的值；將所得Φ值和由式（2）所得單位長度的扭轉(zhuǎn)角θ代入式（9），（10），可求得偏心輪某點的應(yīng)力值.

圖2 偏心輪Fig.2 The eccentric wheel

1.2 “孔洞擬填充”的思想

隨著等值面數(shù)的增加以及區(qū)域Ω（截面域）形狀的多變，結(jié)果的求解也將變得復(fù)雜起來，原因是由于多連域及其內(nèi)部的邊界條件存在，使得結(jié)果的求解處理難以進(jìn)行.究竟采用什么方法才可以有效避免多連域和邊界條件的影響？為此，考慮到有限單元法.因為其可以把計算域劃分為有限個互不重疊的單元，在每個單元內(nèi)選擇一些合適的節(jié)點作為求解函數(shù)的插值點，從而把復(fù)雜問題簡單化；但是將復(fù)雜截面劃分成有限個小單元后，由于不同的小單元具有的材料性能可以不相同，那么能否將復(fù)雜截面中的孔洞進(jìn)行填充，然后轉(zhuǎn)變?yōu)閱芜B域問題求解？據(jù)此出現(xiàn)了“孔洞擬填充”的思想，其實質(zhì)性的原理是將截面內(nèi)邊界Γi所圍成的孔洞部分視作由密度和剪切彈性模量極小的虛擬材料填充而成，然后對此部分虛擬域進(jìn)行單元劃分，再根據(jù)應(yīng)力函數(shù)的連續(xù)性和位移單值條件，將多連域轉(zhuǎn)化為單連域，即將運算轉(zhuǎn)變?yōu)樵谕膺吔绂?所圍成的整個二維面域Ω0上進(jìn)行求解計算.

根據(jù)應(yīng)力函數(shù)理論，可知多連域截面扭轉(zhuǎn)問題的控制方程（歐拉方程）和邊界條件的一般形式為

其中由于是假設(shè)擬填充孔洞，而不同單元具有的材料特性可能不同，故將剪切彈性模量G看作是x，y的函數(shù).當(dāng)采用有限單元法時，將G在每個小單元內(nèi)視為常數(shù)，孔洞區(qū)域的G可取有材料區(qū)域的10－10倍［7］159.若能求解出應(yīng)力函數(shù)Φ，則再根據(jù)式（10）即可得出某一點的應(yīng)力值τ0.為了便于Φ的計算，利用變分原理可以進(jìn)行邊值問題求解，故采用Galerkin法的等效積分從泛函方面進(jìn)行分析求解，由此得到的多連域截面扭轉(zhuǎn)問題的泛函表達(dá)式為

取式（13）的極值條件，由駐值條件δΠ（Φ）＝0可得到關(guān)系式

式中Ω0為截面外邊界Γ0所包圍的區(qū)域.由于有限元法實際上是Ritz法的推廣，所以根據(jù)Ritz法，由公式δΠ（Φ）＝0可得到的近似解為

由式（15）可得單元矩陣關(guān)系式

式中n表示每個小單元節(jié)點個數(shù)；Φe表示單元的應(yīng)力函數(shù)且滿足強(qiáng)制性邊界條件；Νi為試探函數(shù)（形函數(shù)）.

由公式δΠ（Φ）＝0也可以得出關(guān)系式Π＝（Π／Φ1）δΦ1＋（Π／Φ2）δΦ2＋…＋（Π／Φnl）.δΦnl＝0，其中Φ1，Φ2，…，Φn相互獨立，由此可以確定每個小單元的應(yīng)力函數(shù)Φi；所以，對每個小單元應(yīng)用變分原理，利用式（15）和（16）以及Π／Φ＝0可以得到KΦ＝R，其中K表示整個截面的剛度矩陣，它的每一個單元的剛度矩陣都反映了一定的剛度特性，K＝∑eKe，Ke＝，其中Ge表示單元剪切彈性模量；R表示整個截面應(yīng)力函數(shù)的列陣，R＝∑eRe，Re＝2θ∫Ω0NTΦ，其中Φ表示整個截面的節(jié)點應(yīng)力函數(shù)的列陣.此時應(yīng)力函數(shù)Φ可求出，將其代入式（10）即可求得某點的應(yīng)力值.

用孔洞擬填充方法求解，相對于邊界元法，求解偏心輪上某點的應(yīng)力函數(shù)比較簡單，其求解目的是將偏心輪轉(zhuǎn)化為無孔洞的承扭桿件的截面.已知偏心輪的剪切彈性模量為G，將其上的孔洞用剪切彈性模量10－10G的材料進(jìn)行填充，求解過程為：根據(jù)上述分析已知每個小單元的的應(yīng)力函數(shù)Φi為常數(shù)，由于試探函數(shù)Νi事先給定，所以利用KΦ＝R，式（11）以及變分原理式可以求出應(yīng)力函數(shù)Φ，將所得Φ代入式（12），可得到單位長度扭轉(zhuǎn)角θ，然后將求解得到的Φ和θ代入式（9），（10），即可求解出偏心輪某點的應(yīng)力值.

2 基于ANSYS Workbench的實例分析

采用Solidworks建立偏心輪三維模型，圖3為一般的求解模型圖，圖4為虛擬填充模型圖，并將其分別保存為Parasolid格式.將偏心輪的模型圖導(dǎo)入到有限元分析軟件ANSYS Workbench中，進(jìn)行加載、約束、定義分析類型、分析選項、載荷數(shù)據(jù)和載荷選項等步驟，然后開始有限元求解，其中對圖3采用有限元法求解，圖4采用孔洞擬填充法求解.在進(jìn)行分析求解時，偏心輪的材料選擇碳素鋼，2種情況下偏心輪所具有的基本參數(shù)（密度ρ、彈性橫量E、泊松比μ，屈服應(yīng)力σs、破壞應(yīng)力σu）如表1所示.

表1 偏心輪基本參數(shù)Tab.1 Basic parameters of the eccentric with hole and dummy filling

對圖3、圖4兩種情況下的偏心輪進(jìn)行網(wǎng)格劃分，其圖形如圖5～6所示.

圖3 偏心輪的三維模型圖Fig.3 The 3D model of the eccentric wheel

圖4 孔洞擬填充偏心輪的三維模型圖Fig.4 The 3Dmodel of the eccentric wheel with hole and dummy filling

圖5 偏心輪網(wǎng)格Fig.5 Eccentric wheel grid

圖6 孔洞擬填充偏心輪網(wǎng)格Fig.6 The grid of the eccentric wheel with hole and dummy filling

對以上偏心輪的模型加載求解，經(jīng)求解處理后偏心輪在2種情況下所得到的應(yīng)力、位移以及應(yīng)變的分析結(jié)果如圖7～10所示.

圖7 應(yīng)力云圖Fig.7 Stress nephogram

圖8 孔洞擬填充偏心輪的應(yīng)力云圖Fig.8 The stress nephogram of the eccentric with hole and dummy filling

圖9 位移云圖Fig.9 Displacement nephogram

圖10 孔洞擬填充偏心輪的位移云圖Fig.10 The displacement nephogram of the eccentric wheel with hole and dummy filling

對偏心輪采用2種方法計算后所得數(shù)據(jù)見表2，求解后偏心輪的變形曲線和應(yīng)力曲線如圖11～12所示.

通過以上計算對所得云圖以及表2的數(shù)據(jù)進(jìn)行對比分析可知：偏心輪和孔洞擬填充偏心輪在2種工況下對應(yīng)的應(yīng)力曲線圖和變形曲線圖形狀基本相同，最大應(yīng)力值和最大變形值出現(xiàn)的位置相同，并且最大應(yīng)力值都小于許用應(yīng)力值142MPa，從而滿足偏心輪本身所具有的強(qiáng)度要求.由此說明偏心輪在2種工況下的計算結(jié)果具有相似性.利用“孔洞擬填充”思想可以有效處理復(fù)雜截面扭轉(zhuǎn)問題的應(yīng)力值，求解其變形量，從而充分驗證了本文方法的準(zhǔn)確性.

表2 偏心輪的最大應(yīng)力值和變形量Tab.2 Maximum stress value and deformation of the eccentric wheel

［1］ 歐貴寶.多連域截面桿扭轉(zhuǎn)問題的邊界元法［J］.哈爾濱船舶工程學(xué)院學(xué)報，1988，9（4）：469－475.

［2］ SAPOUNTZAKIS E J，TSIPIRAS V J.Composite bars of arbitrary cross section in nonlinear elastic nonuniform torsion by BEM［J］.J Eng Mech，2009，135（12）：1354－1367.

［3］ SAPOUNTZAKIS E J，DIKAROS I C.Non－linear flexural－torsional dynamic analysis of beams of arbitrary cross section by BEM［J］.Int J Nonlinear Mech，2011，46（5）：782－794.

［4］ 楊冬升，凌靜.虛邊界元多連通區(qū)域邊界離散及邊界條件分析［J］.佳木斯大學(xué)學(xué)報：自然科學(xué)版，2011，29（2）：161－163.

［5］ 程耀芳，賴遠(yuǎn)明，王云峰.任意截面形狀直桿扭轉(zhuǎn)問題的解析解［J］.蘭州鐵道學(xué)院學(xué)報，1998，17（2）：13－16.

［6］ 尹昌言.含多個縱向圓孔的圓軸的扭轉(zhuǎn)［J］.東北重型機(jī)械學(xué)院學(xué)報，1989，13（3）：10－18.

［7］ 李永奇，張衛(wèi)紅.模擬求解任意截面桿扭轉(zhuǎn)的化復(fù)為單有限元方法［J］.蘭州理工大學(xué)學(xué)報，2006，32（2）：158－161.

［8］ IDELSOHN S R，ONATE E，CALVO N，et al.The meshless finite element method［J］.Int J Numer Methods Eng，2003，58（6）：893－912.

［9］ 康瀾，張其林，王忠全，等.任意復(fù)雜薄壁截面自由扭轉(zhuǎn)常數(shù)的數(shù)值計算方法［J］.中南大學(xué)學(xué)報：自然科學(xué)版，2011，42（5）：1437－1441.

［10］ KIM Y Y，KIM T S.Topology optimization of beam cross sections［J］.Int J Solids Struct，2000，37（3）：477－493.

［11］ 張德才，周春光，周強(qiáng)，等.基于輪廓的孔洞填充算法［J］.吉林大學(xué)學(xué)報：理學(xué)版，2011，49（1）：82－86.

［12］ HANSBO P，LARSON M G.Discontinuous Galerkin methods for incompressible and nearly incompressible elasticity by Nitsche’s method［J］.Comput Methods Appl Mech Eng，2002，191（17／18）：1895－1908.［13］ LUBARDA V A.On the torsion constant of multicell profiles and its maximization with respect to spar position［J］.Thin Wall Struct，2009，47（6／7）：798－806.

［14］ KIM J H，KIM Y Y.Thin－walled multicell beam analysis for coupled torsion，distortion，and warping deformations［J］.J Appl Mech，2001，68（2）：260－269.

［15］ LI Zhao－xia，KO J M，NI Y Q.Torsional rigidity of reinforced concrete bars with arbitrary sectional shape［J］.Finite Elem Anal Des，2000，35（4）：349－361.

圖11 偏心輪的應(yīng)力曲線Fig.11 Eccentric wheel stress curve

圖12 偏心輪的變形曲線Fig.12 Eccentric wheel deformation curve

Abstract：By the finite element／boundary element method，on the base of torsion problems of the cylinder stress function and converting multiple－connected regions into simple－connected region method，this paper proposes a“hole and dummy filling”method to solve the problem of stress and strain of complicated section.Using eccentric wheel as an example，the problem of torsion is analyzed by ANSYS Workbench finite element software.The study proves that the method is applicable to bars with arbitrary cross－section.

Keywords：complicated section；stress function；multiple connected regions；dummy filling；finite element；torsion

（責(zé)任編輯 賈慧鳴）

Analysis of the complicated section torsional problem with“hole and dummy filling”method

ZHAO Lan－lei1，CHEN Wen－jia1，HE Yan－zhong2，LI Lu－yang1＊，HUANG Xin1
（1.Sch of Mech Engin，Yangzhou Univ，Yangzhou 225127，China；2.Precis Mach Tool Res Inst of Jiangsu Yangli Group，Yangzhou 225127，China）

TH 123.4

A

1007－824X（2012）02－0047－05

2012－01－12

揚州市－揚州大學(xué)科技合作基金（YZ2010154）

＊聯(lián)系人，E－mail：meliluyang＠126.com