軸對稱熱彈性問題雜交基本解Trefftz有限元分析

高可樂, 王克用

(上海工程技術大學 機械與汽車工程學院, 上海 201620)

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

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

1 控制方程及特解理論

1.1 控制方程和邊界條件

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

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

(1)

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

(2)

(3)

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

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

(4)

相應的邊界條件修改為

(5)

(6)

特解不受邊界條件限制,只需滿足關系式

(7)

1.2 離心載荷下的特解

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

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

(8)

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

(9)

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

(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ù)滿足關系式

(14)

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

(15)

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

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

2.1 雙位移場插值模式

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

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

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

(16)

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

(17)

LDLTNe=0

(18)

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

(19)

(20)

其中

(21)

根據(jù)彈性力學位移應變關系及本構關系,相應應力寫成微分算子形式為

(22)

(23)

其中

(24)

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

(25)

2.2 修正變分泛函

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

(26)

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

(27)

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

(28)

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

(29)

其中

(30)

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

(31)

2.3 剛體運動項恢復

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

3 數(shù)值算例

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

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

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

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

3.2 溫度荷載下厚壁圓筒

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

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

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

(32)

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

(33)

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

(34)

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

(35)

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

表1 徑向應力Table 1 Radial stress Pa

表2 周向應力Table 2 Circumferential stress Pa

4 結 語

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