任意多連通截面扭轉(zhuǎn)應(yīng)力函數(shù)的溫度場比擬解法

何 琳，王家林

(重慶交通大學(xué)土木建筑學(xué)院，重慶400074)

在涉及桿件的扭轉(zhuǎn)問題時，求解扭件的應(yīng)力、應(yīng)變狀態(tài)及抗扭剛度十分重要，它是強度、剛度分析和設(shè)計的基礎(chǔ)。對于機械設(shè)備中的軸類扭件，像橫截面為正三角形、圓、橢圓等形狀簡單、規(guī)則的特殊桿件扭轉(zhuǎn)問題，可以采用圣維南逆解法或半逆解法獲得解析解［1］。關(guān)于扭轉(zhuǎn)問題多連通的變分法，林鴻蓀也早有討論。在土木結(jié)構(gòu)工程中，大跨徑橋梁主梁常采用薄壁箱形截面，因主梁截面的抗扭參數(shù)直接影響到結(jié)構(gòu)的抗風(fēng)穩(wěn)定性，需要準(zhǔn)確進行計算［2］。對于橫截面形狀復(fù)雜的桿件，為多連通域或不規(guī)則截面，截面的抗扭參數(shù)卻難以求解，從而使其成為設(shè)計中的難題。許多研究者針對具體的工程問題，采用簡化處理，給出各種近似的解決方法［3-7］。

筆者根據(jù)任意截面桿件扭轉(zhuǎn)問題與二維穩(wěn)態(tài)熱傳導(dǎo)問題在控制方程及邊界條件形式上的相似性，提出了一種利用穩(wěn)態(tài)溫度場分析計算任意多連通截面扭轉(zhuǎn)應(yīng)力函數(shù)的直接比擬有限元方法。通過算例比較，表明此法可以避免“化復(fù)為單［7］”求解應(yīng)力函數(shù)的缺陷，更便捷地解決各種截面的扭轉(zhuǎn)問題。該方法理論概念明確、實用、精度高，具有通用性，為進一步計算任意復(fù)雜截面桿件的扭轉(zhuǎn)剛度和應(yīng)力等問題提供了新的有效手段。

1 桿件扭轉(zhuǎn)問題的應(yīng)力函數(shù)理論

1.1 等直桿扭轉(zhuǎn)問題的控制微分方程

根據(jù)Prandt l應(yīng)力函數(shù)理論，柱體扭轉(zhuǎn)問題的域內(nèi)控制微分方程和邊界條件為［8］:

在域內(nèi):

沿邊界:

式中:Φ(x，y)為Prandt l應(yīng)力函數(shù)。

對于多連通截面，由于有多個邊界，在各邊界的取值可以不同，以Γ0表示截面的外邊界，以Γi(i=1，…，n)表示各孔洞對應(yīng)的內(nèi)邊界。則在各邊界Γi(i=0，1，…，n)上，有:

由于應(yīng)力函數(shù)增加或者減少一個常數(shù)，應(yīng)力分量不受影響，因此不管是單連通截面還是多連通截面，均可令外邊界Γ0上Φ=0(或其它常數(shù))。對于多連通截面，理論上也可令應(yīng)力函數(shù)在任意一個邊界上取0或其它常數(shù)。

1.2 位移單值條件

對于多連通區(qū)域截面，由于只能指定一個邊界上的應(yīng)力函數(shù)為已知常數(shù)，其余邊界上的應(yīng)力函數(shù)屬于未知常數(shù)，其相應(yīng)的式(3)不能組成完備的定解條件，還需要補充位移單值條件才能構(gòu)成完備的定解條件。

等直桿在扭轉(zhuǎn)過程中，橫截面的投影面發(fā)生剛性轉(zhuǎn)動，即截面上的位移為:

式中:α為單位長度扭轉(zhuǎn)角;w(x，y)反映了截面的翹曲情況。

對于位移w=w(x，y)，存在位移單值條件:

式(5)對任意環(huán)線均成立，Γ+表示沿環(huán)線的逆時針方向。

將式(5)應(yīng)用于各邊界，得到多連通截面上Prandtl應(yīng)力函數(shù)的補充定解條件:

對于設(shè)置了應(yīng)力函數(shù)為已知常數(shù)的邊界，盡管式(6)仍然成立，但是由于位移單值條件自動得到滿足，不必將對于該邊界的式(6)作為定解條件。

綜上所述，對于任意連通截面，Prandt l應(yīng)力函數(shù)的完備定解條件如下。

1)在截面內(nèi)，有:

2)在各邊界 Γi(i=0，1，…，n)上Φ = Φi，但是只能在某一個Γk上設(shè)置應(yīng)力函數(shù)的值:

3)對于Γk以外的邊界Γi(i≠k)，有:

對單連通問題，由于只有一個邊界，在外邊界上取Φ=0，式(8)、式(9)自動得到滿足。

2 穩(wěn)態(tài)熱傳導(dǎo)理論

均質(zhì)材料的二維穩(wěn)態(tài)熱傳導(dǎo)問題所滿足的控制方程［9］為:

式中:T為域內(nèi)的溫度函數(shù);ρ為材料密度;Q為物體內(nèi)部的熱源密度;k為材料的熱傳導(dǎo)系數(shù)。

導(dǎo)熱偏微分方程的定解條件有以下3類。

第1類邊界條件是物體邊界上的溫度函數(shù)已知:

式中:Γ為物體邊界，Γ方向為逆時針方向;ˉT是已知溫度值。

第2類邊界條件是物體邊界上的熱流密度q已知:

第3類邊界條件是已知與物體接觸的流體介質(zhì)溫度Tf和換熱系數(shù)α:

3 兩類問題的比較

比較扭轉(zhuǎn)問題控制方程(7)和二維穩(wěn)態(tài)熱傳導(dǎo)問題的控制方程(10)不難發(fā)現(xiàn):式(10)當(dāng)取k=1、ρ=1和Q=2時與式(7)的形式完全相同。扭轉(zhuǎn)問題邊界條件式(8)類似于熱傳導(dǎo)第1類邊界條件式(11)，扭轉(zhuǎn)問題邊界條件式(9)類似熱傳導(dǎo)第2類邊界條件式(12)。

于是，基于兩種問題在控制方程和定解條件方面的相近性，可采用比擬概念，使用通用有限元軟件的穩(wěn)態(tài)熱傳導(dǎo)模塊計算對應(yīng)的溫度場來獲得相應(yīng)扭轉(zhuǎn)問題的應(yīng)力函數(shù)。在有限元軟件ABAQUS中，具體的比擬方法為:

1)取二維穩(wěn)態(tài)熱傳導(dǎo)問題的導(dǎo)熱系數(shù)k=1、材料密度ρ=1、物體內(nèi)部的熱源密度Q=2，使二維穩(wěn)態(tài)熱傳導(dǎo)問題的控制方程與扭轉(zhuǎn)問題完全一致，此時扭轉(zhuǎn)應(yīng)力函數(shù)Φ對應(yīng)于溫度函數(shù)Τ。

2)對于截面的外邊界，設(shè)置其溫度為0。

3)對于各內(nèi)邊界，任取邊界環(huán)線上一點為主控節(jié)點，通過綁定方式將環(huán)線上其它節(jié)點設(shè)置為主控節(jié)點的從節(jié)點，使得整條環(huán)線具有相同的溫度，從而滿足式(3)的條件。

4)多連通截面的位移單值邊條件式(9)對應(yīng)于熱傳導(dǎo)問題的第2類邊界條件式(13)，從溫度場角度，即是令各內(nèi)邊界上所有流入的熱流值等于對應(yīng)邊界所圍成的面積的兩倍，在ABAQUS中的實現(xiàn)方法為:對于各內(nèi)邊界環(huán)線，任取環(huán)線上一點(一般可取主控節(jié)點)設(shè)置集中熱流2Ai。

4 計算實例

下面給出單連通、復(fù)連通和多箱薄壁桿件3個實例的扭轉(zhuǎn)分析，通過數(shù)值解與理論解的比較，充分驗證了本方法的有效性與實用性。計算數(shù)據(jù)還驗證了此法避免了文獻［8］的缺陷，即在輸入過大的熱傳導(dǎo)系數(shù)(如1010以上)時會帶來計算上的不收斂現(xiàn)象，而在輸入的熱傳導(dǎo)系數(shù)較小時(如105)，填充空洞域后，會帶來填充域內(nèi)溫度的不等值現(xiàn)象。

4.1 圓環(huán)截面直桿扭轉(zhuǎn)

扭桿截面如圖1(a)，外徑R1=2 m，內(nèi)徑R2=1 m。有限元計算網(wǎng)格如圖1(b)。應(yīng)力函數(shù)結(jié)果如表1。內(nèi)外邊界處，理論值和模擬計算數(shù)據(jù)完全一致，從A—B間，最大誤差為0.000 526%。本實例結(jié)果表明，通過直接比擬熱場分析來求解扭轉(zhuǎn)問題的應(yīng)力函數(shù)精確度非常高，且無需考慮填充材料熱傳導(dǎo)系數(shù)取值的問題。

圖1 圓環(huán)形扭桿截面、沿A—B的應(yīng)力函數(shù)Fig.1 Cross-section，finite element mesh and stress function along path A—B of annular torsion bar

表1 圓環(huán)形扭桿截面A—B應(yīng)力函數(shù)對比Table 1 Comparison of stress function along path A—B of annular torsion bar

4.2 矩形閉口薄壁桿件

如圖2，矩形閉口薄壁桿截面，截面外邊界長0.1 m，寬0.05 m，壁厚 0.005 m。采用直接比擬法求出的應(yīng)力函數(shù)在截面內(nèi)邊界上完全一致，應(yīng)力函數(shù)結(jié)果為1.570 5E-4，在相同的網(wǎng)格劃分情況下，采用化復(fù)為單比擬解法，當(dāng)填充材料的熱傳導(dǎo)系數(shù)取1.0E+10時，在截面的內(nèi)邊界上，各結(jié)點的應(yīng)力函數(shù)結(jié)果為1.574 8E-4。兩種比擬方法沿A—B的應(yīng)力函數(shù)值(如圖3)，除內(nèi)邊界稍差2.0E-8外，越往外邊界計算結(jié)果越趨一致，最大誤差值(在內(nèi)邊界處)為0.001 27%。

圖2 矩形閉口薄壁桿截面Fig.2 Cross-section of closed thin-walled rectangular bar

圖3 沿A—B的應(yīng)力函數(shù)(×10－6)Fig.3 Stress function along path A—B(×10－6)

對圖2的薄壁桿截面，根據(jù)薄壁柱體扭轉(zhuǎn)的理論［8］可得到內(nèi)壁應(yīng)力函數(shù)值為 1.526 79E-4，與直接比擬法得到的結(jié)果1.570 5E-4誤差2.78%，誤差的原因在于薄壁柱體扭轉(zhuǎn)中假設(shè)了薄壁沿厚度方向各處切應(yīng)力相等，即應(yīng)力函數(shù)法向梯度為常數(shù)。另一方面，誤差較小表明此截面用薄壁柱體扭轉(zhuǎn)理論處理具有可靠的精度。

4.3 多箱薄壁截面扭轉(zhuǎn)

如圖4，多箱矩形閉口薄壁桿截面，箱室具有不同的壁厚，其中 AB厚8 mm，CD厚12 mm，EF厚16 mm，GH及外壁厚10 mm。分別采用直接比擬法、化復(fù)為單的比擬解法(填充材料的熱傳導(dǎo)系數(shù)取1.0E+10)進行計算，表2為在截面3個箱室內(nèi)邊界上的應(yīng)力函數(shù)結(jié)果。通過計算表明，兩種計算方法的結(jié)果非常接近，應(yīng)力函數(shù)值最大誤差0.052 7%。同樣地，文中方法不存在化復(fù)為單法存在熱傳導(dǎo)系數(shù)取值上的計算缺陷。

圖4 三箱薄壁截面示意Fig.4 Cross-section of closed thin-walled rectangular bar with three boxes

表2 內(nèi)壁應(yīng)力函數(shù)值對比Table 2 Results comparison of stress function at inwall

5 結(jié)論

筆者通過比較二維穩(wěn)態(tài)熱傳導(dǎo)問題和扭轉(zhuǎn)問題控制方程及邊界條件的相似性，提出了多連通任意截面直桿扭轉(zhuǎn)問題的應(yīng)力函數(shù)直接求解法，并利用有限元分析軟件ABAQUS的穩(wěn)態(tài)溫度場分析功能，模擬求解了幾種截面的扭轉(zhuǎn)問題。計算結(jié)果表明，該方法簡便實用，能精確地計算任意復(fù)雜截面桿件的扭轉(zhuǎn)應(yīng)力函數(shù)，不受截面形狀的限制，可以避免化復(fù)為單求解扭轉(zhuǎn)問題存在的材料特性取值缺陷，為進一步準(zhǔn)確計算扭轉(zhuǎn)問題的剛度和應(yīng)力等問題提供了新的有效手段。

［1］程昌鈞，王穎堅.彈性力學(xué)［M］.北京:高等教育出版社，1999:235-282.

［2］Cheng S H，LAu D T，Cheung M S.Comparison of numerical techniques for 3D flutter analysis of cable-stayed bridges［J］.Computers＆ Structures，2003，81(32):2811-2822.

［3］李金伯，何西泠，余國城.復(fù)雜形體軸的扭轉(zhuǎn)剛度［J］.工程機械，1998，29(4):11-13.Li Jinbai，He xileng，Yu guocheng.Torsional rigidity of shafts of complex profile［J］.Construction Machinery and Equipment，1998，29(4):11-13.

［4］程耀芳，賴遠明，王云峰.任意截面形狀直桿扭轉(zhuǎn)問題的解析解［J］.蘭州鐵道學(xué)院學(xué)報，1998，17(2):13-16.Cheng Yaofang，Lai Yuanming，Wang Yunfeng.Analytical solution of torsion problems of arbitrary shape cross section bars［J］.Journal of Lanzhou Railway University，1998，17(2):13-16.

［5］劉悅藏，范慕輝.凸多邊形截面桿扭轉(zhuǎn)問題的數(shù)值解法［J］.河北工業(yè)大學(xué)學(xué)報，1999，28(6):73-75.Liu Yuecang，F(xiàn)an Muhui.The numerical method for solving the torsional shafts with convex polygon cross sections［J］.Journal of Hebei University of Technology，1999，28(6):73-75.

［6］李瑞選.多孔直桿扭轉(zhuǎn)問題的邊界元法［J］.華東理工大學(xué)學(xué)報，1995，21(5):640-647.Li Ruixia.Boundary element methods for elastic torsion of multihole bars［J］.Journal of East China University of Science and Technology，1995，21(5):640-647.

［7］李永奇，張衛(wèi)紅.模擬求解任意截面桿扭轉(zhuǎn)的化復(fù)為單有限元方法［J］.蘭州理工大學(xué)學(xué)報，2006，32(2):158-161.Li Yongqi，Zhang Weihong.A special simulation solution for torsional bars with arbitrary cross-section using finite element method［J］.Journal of Lanzhou University of Technology，2006，32(2):158-161.

［8］錢偉長，葉開沅.彈性力學(xué)［M］.北京:科學(xué)出版社，1956:143-152.

［9］孔祥謙.有限單元法在傳熱學(xué)中的應(yīng)用［M］.北京:科學(xué)出版社，1986:1-7.