高可樂, 王克用

(上海工程技術(shù)大學(xué) 機(jī)械與汽車工程學(xué)院, 上海 201620)

近年來,Trefftz有限元法因兼有傳統(tǒng)有限法和邊界元法的諸多優(yōu)點而備受關(guān)注,其最早可追溯到1926年Trefftz對Laplace方程的求解.1977年Jirousek等[1]在對薄板體彎曲問題的研究中正式提出雜交Trefftz有限元法(HT-FEM).與傳統(tǒng)有限元法不同,雜交Trefftz有限元法采用雙位移場插值模式,即單元域內(nèi)場和輔助網(wǎng)線場.單元域內(nèi)場采用控制方程的齊次解構(gòu)建形函數(shù),使得非常稀疏網(wǎng)格和相對少量自由度情況下仍能獲得高精度解.此外,僅含邊界積分的有限元格式顯著提高了Trefftz單元的抗畸變能力.目前,Trefftz有限元法已成功應(yīng)用于位勢問題[2]、平面彈性問題[3]、接觸問題[4-5]、軸對稱問題[6-8]和復(fù)合材料問題[9-12].

汽輪機(jī)、火箭、航空發(fā)動機(jī)等裝備中的機(jī)械結(jié)構(gòu)常由于工作溫度變化引起熱荷載和機(jī)械荷載雙重作用,對其進(jìn)行熱彈性分析十分重要.一般熱彈性問題中溫度場和變形場是耦合的,從數(shù)學(xué)角度講熱彈性問題是非齊次方程求解問題.借助Trefftz有限元法求解時,作為體力項的溫度荷載會添加域積分到單元剛度方程中,使得雜交基本解有限元法僅含邊界積分的優(yōu)勢消失.繼區(qū)域離散法、雙重互易法之后,特解法作為消除域積分的有效方法之一,由Qin等[13]引入到雜交Trefftz有限元法中.此后,Wang等[14]應(yīng)用徑向基函數(shù)獲得特解研究正交各向異性位勢問題;Wang等[15]用特解法求解極小曲面問題;劉博等[16]應(yīng)用特解法分析軸對稱Poisson方程問題的求解;Zhou等[6-8]利用基本解構(gòu)造單元域內(nèi)插值函數(shù),提出雜交基本解有限元法(HFS-FEM),并分析軸對稱位勢問題以及軸對稱彈性力學(xué)問題.本文基于雜交基本解有限元格式分析熱彈性問題,通過應(yīng)用特解法以消除其中的域積分.

1 控制方程及特解理論

1.1 控制方程和邊界條件

線彈性問題的控制方程可寫成

(λ+μ)uj,ji+μui,jj+bi=0

(1)

式中:u為位移;bi為非齊次項;λ和μ分別為拉梅常數(shù)和剪切彈性模量.逗號表示求導(dǎo),i,j=x,y,z為笛卡爾空間直角坐標(biāo)系的3個維度.考慮邊界條件為

(2)

(3)

式中:t為表面力;Γ=Γu∪Γt,Γu∩Γt=?;Γ為求解域Ω的邊界;Γu和Γt分別為Dirichlet邊界條件和Neumann邊界條件,變量上方橫線表示給定值.

在邊界條件(2)和(3)下,式(1)的解可表達(dá)為齊次解和特解的疊加形式.其中,齊次解uh滿足

(4)

相應(yīng)的邊界條件修改為

(5)

(6)

特解不受邊界條件限制,只需滿足關(guān)系式

(7)

1.2 離心載荷下的特解

為分析軸對稱問題,首先在笛卡爾直角坐標(biāo)系下給出離心力分量并求出相應(yīng)位移特解,然后轉(zhuǎn)化為圓柱坐標(biāo)系下的特解形式.

假設(shè)軸對稱物體繞z軸以角速度ω旋轉(zhuǎn),且物體密度為ρ,則體力b分量可寫成

(8)

文獻(xiàn)[14]分析了控制方程(1)的特解形式,本文給出笛卡爾空間直角坐標(biāo)系下一個特解為

(9)

將式(9)轉(zhuǎn)化為圓柱坐標(biāo)系下形式,則位移分量特解可寫成

(10)

1.3 溫度荷載下的特解

單獨溫度荷載下,控制方程(1)中非齊次項可寫成

bi=βT,i

(11)

將式(11)代入式(1)中,可得到位移分量表示的熱彈性平衡微分方程為

(λ+μ)uj,ji+μui,jj=-βT,i

(12)

其中

β=α(3λ+2μ)

(13)

式中:α為熱膨脹系數(shù);T為溫度.

式(12)的特解可通過熱彈性位移勢函數(shù)Φ(x,y,z)確定.該函數(shù)滿足關(guān)系式

(14)

一旦找到符合條件的熱彈性位移勢,則位移特解可由式(14)計算得到.為獲得上述熱彈性位移勢,將式(14)代入式(12),化簡后得到熱彈性位移勢與溫度場的關(guān)系式為

(15)

式(15)為標(biāo)準(zhǔn)Poisson方程形式.若溫度場已知,根據(jù)式(14)和(15)即可求出位移特解,進(jìn)而通過位移應(yīng)變關(guān)系和本構(gòu)關(guān)系求出應(yīng)力特解.需要注意的是,若溫度場位移勢為簡單函數(shù)形式,熱彈性位移勢可通過式(15)直接獲得解析式.若溫度場解析解比較復(fù)雜甚至不存在,則需采用近似函數(shù)(如徑向基函數(shù)或多項式函數(shù))求解近似熱彈性位移勢[16].

2 雜交基本解Trefftz有限元格式

2.1 雙位移場插值模式

雜交基本解有限元法采用雙位移場插值模式,如圖1所示.單元域內(nèi)場插值函數(shù)通過精確滿足控制方程的截斷完備解或基本解來構(gòu)造,而輔助網(wǎng)線場用以保證單元間的連續(xù)性,其插值函數(shù)采用常規(guī)方法建立.

圖1 軸對稱體與8節(jié)點環(huán)狀單元Fig.1 Axisymmetric body and 8-node annular element

由于控制方程是非齊次的,單元域內(nèi)場需將特解考慮進(jìn)去,可寫成

(16)

式中:下標(biāo)e代表單元內(nèi)變量.輔助網(wǎng)線場不受非齊次項影響,其特解部分包含在單元節(jié)點自由度列陣de中,可寫成

(17)

LDLTNe=0

(18)

式中:L為微分算子矩陣;D為彈性矩陣,其形式為

(19)

(20)

其中

(21)

根據(jù)彈性力學(xué)位移應(yīng)變關(guān)系及本構(gòu)關(guān)系,相應(yīng)應(yīng)力寫成微分算子形式為

(22)

(23)

其中

(24)

輔助網(wǎng)線場可通過自然坐標(biāo)系ξ∈[-1,1]來構(gòu)造.對圖1所示的8節(jié)點四邊形單元,每條邊布置3個節(jié)點,可構(gòu)造二次網(wǎng)線函數(shù)為

(25)

2.2 修正變分泛函

單元域內(nèi)場和輔助網(wǎng)線場之間的聯(lián)系是通過修正變分泛函實現(xiàn)的.控制方程為齊次時,通過高斯散度定理可去除變分泛函中的域積分.然而,控制方程為非齊次時,變分泛函中的域積分無法直接去除.為解決這個問題,本文只構(gòu)造原問題齊次情形的變分泛函,暫時排除邊界條件中因非齊次項誘發(fā)的特解部分.整個求解域?qū)?yīng)的變分泛函可寫成所有單元泛函的疊加形式為

(26)

考慮式(18),對式(26)第一項應(yīng)用高斯散度定理可得

(27)

將式(27)代入式(26),原泛函可簡化為僅含邊界積分的形式為

(28)

將式(16)、式(17)、式(22)和式(23)代入式(28),可得

(29)

其中

(30)

對式(29)應(yīng)用兩次駐值原理,可得待定系數(shù)列陣和單元剛度方程為

(31)

2.3 剛體運(yùn)動項恢復(fù)

將上述單元剛度方程組裝成總體剛度方程,采用乘大數(shù)法引入位移邊界條件即可求得所有節(jié)點的齊次解,進(jìn)一步與獲得的特解疊加便可求得節(jié)點處的位移全解.需要注意的是,單元域內(nèi)任一點處的位移還需恢復(fù)剛體位移項,具體原因和恢復(fù)方法參見文獻(xiàn)[7-9],本文因篇幅所限不再復(fù)述.

3 數(shù)值算例

3.1 離心載荷作用下帶孔圓球

為分析離心載荷下物體內(nèi)部位移分布,給出帶圓柱孔的球體結(jié)構(gòu),如圖2(a)所示.由于對稱性,只分析1/4結(jié)構(gòu),圓柱孔半徑r0=1,圓球半徑為r0+1.本例分析中,彈性模量E=1,泊松比v=0.3,邊界條件及網(wǎng)格剖分如圖2(b)所示.采用雜交基本解有限元法求解時,需要注意單元域外源點位置.由于這些源點位置與網(wǎng)格單元大小有關(guān),為避免奇異性,子午面內(nèi)的源點不得配置在z軸左側(cè).圖3給出雜交基本解有限元和ABAQUS位移計算結(jié)果,二者吻合良好.

圖2 帶孔圓球及其網(wǎng)格剖分Fig.2 Perforated sphere and mesh configuration

圖3 總位移云圖Fig.3 Total displacement nephogram

3.2 溫度荷載下厚壁圓筒

本節(jié)對厚壁圓筒內(nèi)熱應(yīng)力的分布情況進(jìn)行分析.厚壁圓筒尺寸如圖4所示.厚壁圓筒內(nèi)徑a=1,外徑b=a+1,高h(yuǎn)=0.5,內(nèi)外表面為自由邊界條件,上下表面徑向固定.彈性模量和泊松比與3.1節(jié)中算例相同.

對于熱傳導(dǎo)問題,若厚壁圓筒內(nèi)表面溫度變化為t1=Ta,t2=0,外表面溫度變化為零,厚壁圓筒的溫度場解析解可表示為

圖4 厚壁圓筒及其網(wǎng)格劃分Fig.4 Thick-walled cylinder and mesh configuration

(32)

將式(32)代入式(15),并考慮圓柱坐標(biāo)系,可得方程為

(33)

對式(33)兩次積分,可得熱彈性位移勢為

(34)

求得熱彈性位移勢后,位移特解可通過式(14)獲得,應(yīng)力特解通過位移應(yīng)變關(guān)系和本構(gòu)關(guān)系推得.為方便對比,給出此溫度場下的應(yīng)力解析解為

(35)

解析解和雜交基本解有限元解見表1和表2.兩者對比,取相同精度時結(jié)果十分吻合.

表1 徑向應(yīng)力Table 1 Radial stress Pa

表2 周向應(yīng)力Table 2 Circumferential stress Pa

4 結(jié) 語

針對軸對稱熱彈性問題,引入熱彈性位移勢得到一組位移和應(yīng)力特解.若溫度場解析解比較復(fù)雜甚至不存在,可根據(jù)熱彈性位移勢導(dǎo)出Poisson方程進(jìn)而求得近似特解.將特解與雜交基本解有限元列式結(jié)合獲得齊次解,利用線性疊加原理即可求得全解.本文方法求解過程理論清晰,編程簡易,由特解法去除因溫度荷載引起的域積分后,雜交基本解有限元法的固有優(yōu)點得以保持.