胡 寬

(同濟大學土木工程學院，上海 200082)

0 引言

扭轉(zhuǎn)是結(jié)構(gòu)構(gòu)件一種常見的受力狀態(tài)。在針對桿件的扭轉(zhuǎn)分析中，抗扭慣性矩是基本且重要的一個力學特性，對于簡單的截面形式(如矩形，圓形，圓環(huán))以及鋼結(jié)構(gòu)等薄壁截面，均有較為成熟的計算方法。但是對于非規(guī)則的多連通或形狀復(fù)雜的截面，難以得到抗扭慣性矩的解析解[1]。Prandtl曾提出了薄膜比擬的實驗方法，并設(shè)計了實驗設(shè)備，但是操作復(fù)雜、測量精度不高。有學者曾提出有限差分法(FDM)，該計算方法簡單而有效，在單連通截面的分析中應(yīng)用廣泛，但不適合復(fù)聯(lián)通截面[2]。另有研究者針對具體的工程問題，進行簡化處理，給出近似解，但不具有通用性[3-7]。在當前計算機有限元分析與智能化設(shè)計的潮流下，有必要探索出一種適用性強、可程序化的方法。本文從彈性力學理論出發(fā)，基于薄膜理論的思想、有限元方法，提出一種適用于任意復(fù)聯(lián)通截面的高精度計算方法。選取三個算例，證明該方法的通用性與準確性。

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

1.1 應(yīng)力函數(shù)方程

由應(yīng)力函數(shù)理論可得[8]，對于圖1所示的任意復(fù)連通截面形狀的桿件在扭轉(zhuǎn)作用下的應(yīng)力函數(shù)方程與邊界條件應(yīng)滿足：

(1)

其中，φ為應(yīng)力函數(shù);G為材料剪切模量;θ為單位長度扭轉(zhuǎn)角;Γi為截面內(nèi)外邊界;Ci為未知常數(shù)。應(yīng)力函數(shù)與桿件承受扭矩M應(yīng)滿足:

(2)

其中，Ω為截面域(除去孔洞);Ωi為截面的孔洞。

1.2 薄膜比擬方法

假設(shè)有一片張緊的均勻無重力薄膜，受到垂直薄膜平面的均勻壓力?？傻肹9]，薄膜的垂度函數(shù)z和薄膜變形后與邊界平面包圍的體積V關(guān)系滿足:

(3)

(4)

其中,q為薄膜受到的均勻壓力;FT為薄膜內(nèi)張力;Γ為薄膜邊界;Ω為薄膜平面域。

若令q/FT=2Gθ，可注意到，薄膜垂度函數(shù)與扭轉(zhuǎn)作用下的桿件截面應(yīng)力函數(shù)滿足同樣的微分方程與邊界條件，因此必然具有相同且唯一的解[9]，即φ≡z。同理，扭矩M必然等于薄膜與邊界平面包圍體積的兩倍，即M=2V。在物理含義上，應(yīng)力函數(shù)的求解問題可以轉(zhuǎn)化為薄膜垂度函數(shù)的計算問題。引入有限單元法的分析思想，借助網(wǎng)格劃分算法，即可以得到應(yīng)力函數(shù)的數(shù)值解，進而計算抗扭慣性矩。

2 扭轉(zhuǎn)問題的有限元理論

有限單元法分析的一個重要步驟是結(jié)構(gòu)的離散化。本文選取Delaunay三角化方法，該方法可以適應(yīng)各種邊界、多連通域等復(fù)雜情況，將任意截面形狀離散化為三角形網(wǎng)格。

2.1 基于三角形薄膜單元的有限元列式

根據(jù)式(1)，應(yīng)用伽遼金法和變分基本原理，可得:

(5)

令轉(zhuǎn)角θ=-1，假設(shè)單元形函數(shù)矩陣N，應(yīng)變矩陣B，單元垂度自由度Ze，則在單元內(nèi)部Ωe有:

(6)

針對三角形單元，采用一次位移插值形函數(shù)，由上式和單元平衡方程KeZe=Pe，推導(dǎo)單元剛度矩陣和荷載列向量為:

(7)

Pe=2GA[1/3 1/3 1/3]T

(8)

其中，A為三角形面積。剛度矩陣中元素為：

kij=bibj+cicj
bi=xj-xm
ci=yj-ym

(9)

其中，xi，yi為三角形單元的節(jié)點坐標;i，j，m均為下標輪換，如i→j，j→m，m→i。

2.2 邊界條件的處理

對于單聯(lián)通截面，可以令所有邊界上的節(jié)點位移為零來求解方程。本論文采用對總體剛度矩陣對角元素改一法處理強制邊界條件[10]。當截面存在內(nèi)環(huán)時，應(yīng)該給內(nèi)環(huán)上所有節(jié)點增加位移耦合條件，并將環(huán)內(nèi)以及環(huán)上全部荷載移加到環(huán)上某一節(jié)點，然后令其他環(huán)山節(jié)點位移等于該節(jié)點。由上式可知，環(huán)內(nèi)荷載總和等于2GAi，其中，Ai為內(nèi)環(huán)面積，可用任意多邊形的面積計算公式獲得[11]。

施加節(jié)點位移耦合條件時，應(yīng)將從節(jié)點在總體剛度矩陣中對應(yīng)的行列元素附加到主節(jié)點對應(yīng)行列上，然后給從節(jié)點添加強制邊界條件。求解結(jié)束后，再將主節(jié)點位移賦值給從節(jié)點。

3 計算實例

基于C#編程語言實現(xiàn)該計算方法，并針對三個算例進行測試。

對于規(guī)則截面(算例1)使用理論公式進行正確性驗證；對于非規(guī)則工程用截面(算例2,算例3)，參考文獻[2]的思路，使用ANSYS有限元分析軟件建立實體懸臂模型進行驗證。選取實體懸臂梁長度L1=50 m，使用Solid95單元，懸臂端連接一個梁單元，長度為L2=1 m，使用Beam188單元。實體單元與梁單元交界處設(shè)置剛域(Regid Region)。剛域僅約束X，Y方向(截面平面內(nèi))的平動自由度，以避免箱梁的翹曲效應(yīng)的影響。在梁單元端部施加扭矩T=1 000 kN·m，通過計算實體梁在單位長度內(nèi)的扭轉(zhuǎn)角，推算抗扭慣性矩。

驗證結(jié)果表明，該計算方法對于任意復(fù)聯(lián)通截面的扭轉(zhuǎn)慣性矩結(jié)果是正確的，且具有相當高的計算精度。算例驗證如下。

3.1 圓環(huán)截面

計算外徑為D=10，壁厚為t的圓環(huán)的抗扭慣性矩。其中t分別取0.2,0.5,1三個值。對截面使用Delaunay三角網(wǎng)格化的結(jié)果見圖2，計算結(jié)果與理論值對比見表1。

表1 計算結(jié)果與誤差

壁厚t節(jié)點數(shù)單元數(shù)抗扭慣性矩計算值理論值誤差/%0.2256256147.8147.9-0.040.5128128336.9337.6-0.201.096128576.3579.6-0.58

3.2 單箱室箱型截面

選取截面為箱梁，幾何尺寸見圖3，三角網(wǎng)格化結(jié)果見圖4。

在ANSYS中建立實體梁模型見圖5，該方法計算結(jié)果與ANSYS計算結(jié)果對比見表2。

表2 計算結(jié)果與誤差

節(jié)點數(shù)單元數(shù)抗扭慣性矩/cm4計算值A(chǔ)NSYS誤差/%7321 0533.329E+93.461E+9-3.82

3.3 多箱室箱型截面

選取截面為多箱室箱梁，幾何尺寸見圖6，該方法計算結(jié)果與ANSYS計算結(jié)果對比見表3。

表3 計算結(jié)果與誤差

節(jié)點數(shù)單元數(shù)抗扭慣性矩/cm4計算值A(chǔ)NSYS誤差/%8591 2031.136E+91.160E+9-2.04

4 結(jié)語

基于彈性力學應(yīng)力函數(shù)方程和薄膜比擬方法，可以將任意截面抗扭慣性矩計算問題轉(zhuǎn)化為一維有限元分析問題。

經(jīng)過驗證，該計算方法適用于復(fù)聯(lián)通截面，具有較高的計算精度，可以滿足工程需求。通過在C#編程語言上的實現(xiàn)，說明該方法適合程序化應(yīng)用。