任意多連通截面扭轉應力函數的溫度場比擬解法

何 琳，王家林

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

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

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

1 桿件扭轉問題的應力函數理論

1.1 等直桿扭轉問題的控制微分方程

根據Prandt l應力函數理論，柱體扭轉問題的域內控制微分方程和邊界條件為［8］:

在域內:

沿邊界:

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

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

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

1.2 位移單值條件

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

等直桿在扭轉過程中，橫截面的投影面發生剛性轉動，即截面上的位移為:

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

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

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

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

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

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

1)在截面內，有:

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

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

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

2 穩態熱傳導理論

均質材料的二維穩態熱傳導問題所滿足的控制方程［9］為:

式中:T為域內的溫度函數;ρ為材料密度;Q為物體內部的熱源密度;k為材料的熱傳導系數。

導熱偏微分方程的定解條件有以下3類。

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

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

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

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

3 兩類問題的比較

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

于是，基于兩種問題在控制方程和定解條件方面的相近性，可采用比擬概念，使用通用有限元軟件的穩態熱傳導模塊計算對應的溫度場來獲得相應扭轉問題的應力函數。在有限元軟件ABAQUS中，具體的比擬方法為:

1)取二維穩態熱傳導問題的導熱系數k=1、材料密度ρ=1、物體內部的熱源密度Q=2，使二維穩態熱傳導問題的控制方程與扭轉問題完全一致，此時扭轉應力函數Φ對應于溫度函數Τ。

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

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

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

4 計算實例

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

4.1 圓環截面直桿扭轉

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

圖1 圓環形扭桿截面、沿A—B的應力函數Fig.1 Cross-section，finite element mesh and stress function along path A—B of annular torsion bar

表1 圓環形扭桿截面A—B應力函數對比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。采用直接比擬法求出的應力函數在截面內邊界上完全一致，應力函數結果為1.570 5E-4，在相同的網格劃分情況下，采用化復為單比擬解法，當填充材料的熱傳導系數取1.0E+10時，在截面的內邊界上，各結點的應力函數結果為1.574 8E-4。兩種比擬方法沿A—B的應力函數值(如圖3)，除內邊界稍差2.0E-8外，越往外邊界計算結果越趨一致，最大誤差值(在內邊界處)為0.001 27%。

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

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

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

4.3 多箱薄壁截面扭轉

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

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

表2 內壁應力函數值對比Table 2 Results comparison of stress function at inwall

5 結論

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

［1］程昌鈞，王穎堅.彈性力學［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］李金伯，何西泠，余國城.復雜形體軸的扭轉剛度［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］程耀芳，賴遠明，王云峰.任意截面形狀直桿扭轉問題的解析解［J］.蘭州鐵道學院學報，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］劉悅藏，范慕輝.凸多邊形截面桿扭轉問題的數值解法［J］.河北工業大學學報，1999，28(6):73-75.Liu Yuecang，Fan 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］李瑞選.多孔直桿扭轉問題的邊界元法［J］.華東理工大學學報，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］李永奇，張衛紅.模擬求解任意截面桿扭轉的化復為單有限元方法［J］.蘭州理工大學學報，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］錢偉長，葉開沅.彈性力學［M］.北京:科學出版社，1956:143-152.

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