張 凱, 王克用, 齊東平

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

0 引 言

雜交Trefftz有限元法(HT-FEM)最早由Jirousek和Leon[1]在對(duì)薄板體彎曲問(wèn)題的研究中提出,該方法融合了傳統(tǒng)有限元法(FEM)和邊界元法(BEM)的諸多優(yōu)點(diǎn),并且摒棄了它們的一些缺點(diǎn)[2-4]．對(duì)于某些物理問(wèn)題,雜交Trefftz有限元法很難得到相應(yīng)的完備解,且需要選擇合適的Trefftz項(xiàng)數(shù)來(lái)獲得預(yù)期結(jié)果．為克服這一缺點(diǎn),近年來(lái),基于基本解的雜交有限元法(HFS-FEM)得到了廣泛的關(guān)注[4-6]．與雜交Trefftz有限元法相似,雜交基本解有限元法也假設(shè)兩個(gè)獨(dú)立的溫度場(chǎng):單元域內(nèi)溫度場(chǎng)和單元邊界網(wǎng)線場(chǎng)．不同之處在于,后者利用控制方程基本解的線性組合來(lái)近似單元域內(nèi)場(chǎng),而前者采用完備解．引入修正變分泛函將上述兩個(gè)場(chǎng)關(guān)聯(lián)起來(lái),導(dǎo)出僅含邊界積分的有限元列式,降低了求解維度,從而減少了計(jì)算量,稀疏網(wǎng)格下也能獲得理想的精度[6-11]．顯然,該方法幾乎繼承了雜交Trefftz有限元法的所有優(yōu)點(diǎn),且規(guī)避了難以構(gòu)造Trefftz函數(shù)的問(wèn)題．到目前為止,雜交基本解有限元法已成功應(yīng)用于熱傳導(dǎo)問(wèn)題[4]、平面彈性問(wèn)題[12]、熱彈性問(wèn)題[13]、軸對(duì)稱問(wèn)題[7-8,14]和裂紋問(wèn)題[15]等．

由于基本解的奇異性[16],適當(dāng)數(shù)量的源點(diǎn)需布置在單元域外,以避免基本解涉及奇異積分問(wèn)題．根據(jù)以往的研究,源點(diǎn)的布局主要有兩種方式:一種是假設(shè)與單元邊界形狀相似的虛擬邊界(偏置邊界),另一種是在單元域外假設(shè)一個(gè)圓形虛擬邊界[17]．虛擬邊界偏移量和源點(diǎn)數(shù)目的選擇直接影響計(jì)算結(jié)果的準(zhǔn)確性和穩(wěn)定性,需要根據(jù)問(wèn)題具體分析,目前尚無(wú)準(zhǔn)確的數(shù)值供選擇．上述雜交基本解有限元法均采用第一種虛擬邊界求解問(wèn)題,而對(duì)其他源點(diǎn)布局未作闡述．在求解扭轉(zhuǎn)彈性問(wèn)題[18]中,通過(guò)對(duì)比上述兩種邊界方式下的計(jì)算結(jié)果,認(rèn)為當(dāng)源點(diǎn)布局方式與單元域形狀相似時(shí),可得到良好的結(jié)果．

在問(wèn)題求解中,一種源點(diǎn)的方案可能不適用于所有的問(wèn)題,當(dāng)其失效時(shí),可以考慮其他源點(diǎn)布局方式．為了探究不同源點(diǎn)對(duì)雜交基本解有限元法的影響,本文結(jié)合上述兩種虛擬邊界以及雙重虛擬邊界[19]分析了熱傳導(dǎo)問(wèn)題．通過(guò)兩個(gè)算例驗(yàn)證了不同源點(diǎn)布局方式的可行性,探討了源點(diǎn)數(shù)目和偏移量對(duì)計(jì)算精度的影響,以及虛擬邊界參數(shù)對(duì)單元?jiǎng)偠染仃嚄l件數(shù)的影響,同時(shí)也分析了不同源點(diǎn)布局下的收斂性．此外,還說(shuō)明了當(dāng)網(wǎng)格畸變程度較大時(shí),哪種源點(diǎn)布局更適合求解．

1 問(wèn) 題 描 述

1.1 控制方程與邊界條件

對(duì)于穩(wěn)態(tài)熱傳導(dǎo)問(wèn)題,一般二維區(qū)域Ω的Laplace控制方程為

(1)

考慮Dirichlet和Neumann邊界條件:

(2)

(3)

式中,u和q分別代表未知溫度和熱流,字母上方的橫線代表已知邊界值,k為導(dǎo)熱系數(shù),Γ=Γu∪Γq為求解域Ω圍成的整個(gè)邊界,nx1和nx2分別代表邊界Γ上任意點(diǎn)外法線方向余弦,x=[x1,x2]T．

1.2 基本解及源點(diǎn)布局

對(duì)于均質(zhì)各向同性材料,Laplace方程(1)的基本解[4,13]為

(4)

類型1 單元域外假想一個(gè)與單元邊界形狀相似的虛擬邊界,源點(diǎn)ysj(j=1,2,…,ns)位置可由基本解法確定:

ysj=xb+λ(xb+xc),

(5)

式中,λ為無(wú)量綱參數(shù),決定著源點(diǎn)到單元邊界的距離,xb為單元邊界上的場(chǎng)點(diǎn),xc為單元形心．圖1(a)為一個(gè)典型的單元源點(diǎn)分布．

類型2 單元域外假想一個(gè)半徑為R的圓形虛擬邊界,其形式如圖1(b)所示,源點(diǎn)和場(chǎng)點(diǎn)滿足下列方程式:

(6a)

(6b)

類型3 雙層虛擬邊界,與類型2類似,在單元域外假想兩個(gè)半徑分別為R1和R2的圓,其形式如圖1(c)所示,源點(diǎn)分別布置在兩個(gè)圓上:

圖1 源點(diǎn)配置方式Fig. 1 Configurations of source points

(7a)

(7b)

式中,θ為單元邊界節(jié)點(diǎn)的角度．為了便于區(qū)分,假設(shè)R1

為了更好地對(duì)源點(diǎn)布局進(jìn)行比較分析,構(gòu)造了以無(wú)量綱參數(shù)λ為變量的關(guān)系式:

(8)

式中,dmax為單元邊上點(diǎn)到單元形心的最大距離,d為雙層虛擬邊界之間的距離,即d=R2-R1=αR2,若無(wú)特殊說(shuō)明,α=0.1．

2 雜交有限元列式

2.1 假設(shè)溫度場(chǎng)

類似于Trefftz有限元法,雜交基本解有限元法將所考慮的問(wèn)題域劃分為一系列小的單元,采用兩個(gè)獨(dú)立的插值模式(單元域內(nèi)溫度場(chǎng)和輔助網(wǎng)線溫度場(chǎng))．圖2為一個(gè)典型的四節(jié)點(diǎn)單元．單元域內(nèi)溫度場(chǎng)由基本解的線性組合構(gòu)成,保證了單元域內(nèi)變量場(chǎng)的計(jì)算精度; 單元邊界定義的輔助網(wǎng)線溫度場(chǎng)確保了單元間的連續(xù)性,一般采用常規(guī)有限元法的形函數(shù)[7]．

1) 單元域內(nèi)溫度場(chǎng)

(9)

式中,ue為單元域內(nèi)的溫度,Ne為基本解Ne(x,yj)=u*(x,y)的線性組合構(gòu)成的行向量,ns為單元域外源點(diǎn)個(gè)數(shù),ce為待定參數(shù)構(gòu)成的列向量,Ωe為單元邊界Γe圍成的單元區(qū)域．

利用式(9),單元邊界Γe上的外法向熱流qe為

(10)

式中

(11)

2) 輔助網(wǎng)線溫度場(chǎng)

(12)

圖2 典型四節(jié)點(diǎn)單元和單元邊形函數(shù)Fig. 2 A typical 4-node element and the element side shape functions

(13)

2.2 修正變分泛函

整個(gè)求解域Ω的雜交泛函Πm為所有單元泛函Πme的疊加,即Πm=∑eΠme．為將上述兩個(gè)溫度場(chǎng)關(guān)聯(lián)在一起,單元泛函可表示為

(14)

(15)

將式(9)、(10)和(12)代入式(15),有

(16)

式中

(17)

為保證單元間連續(xù)性,未知向量ce由單元節(jié)點(diǎn)自由度向量de表示,對(duì)泛函(16)應(yīng)用駐值原理分別消去ce和de,即

(18)

(19)

2.3 恢復(fù)剛體運(yùn)動(dòng)項(xiàng)

ue=c0+Nece,

(20)

(21)

進(jìn)一步得知

(22)

式中,ne為單元的節(jié)點(diǎn)個(gè)數(shù),一旦通過(guò)整體剛度矩陣求得節(jié)點(diǎn)的溫度,向量ce即可通過(guò)式(18)確定,然后c0可由式(22)求得,最后可通過(guò)式(20)求出單元域內(nèi)任一點(diǎn)的溫度ue．

3 數(shù) 值 算 例

為驗(yàn)證不同源點(diǎn)布局方式下本文方法的可行性,考慮了兩個(gè)數(shù)值算例:正方形板和偏心環(huán)空的熱傳導(dǎo)．若無(wú)特別說(shuō)明,導(dǎo)熱系數(shù)取為k=1 W/(m·K)．為了從量化角度理解計(jì)算結(jié)果的準(zhǔn)確性,引入了任意變量f的平均相對(duì)誤差(ARE):

(23)

式中,Nt為測(cè)試點(diǎn)的數(shù)目,(fHFS-FEM)i和(fref)i分別為雜交基本解有限元法和ABAQUS在點(diǎn)i的數(shù)值解．

3.1 正方形板

第一個(gè)算例考慮的是正方形板,其幾何尺寸和邊界條件如圖3所示．板的上下表面均施加Dirichlet邊界條件,右表面施加Neumann邊界條件．為研究該方法的收斂性,將正方形板分別離散為2×2,4×4,6×6和8×8網(wǎng)格．

圖4 不同源點(diǎn)數(shù)目的四節(jié)點(diǎn)單元Fig. 4 A 4-node element with different numbers of source points

圖5 無(wú)量綱參數(shù)λ對(duì)計(jì)算精度的影響 圖6 無(wú)量綱參數(shù)λ對(duì)矩陣He條件數(shù)的影響 Fig. 5 Effects of dimensionless parameter λ on Fig. 6 Effects of dimensionless parameter λ on the computation accuracy the condition number of matrix He

為研究本文方法對(duì)網(wǎng)格畸變的不敏感性,圖7給出了5種4×4網(wǎng)格畸變方案,由無(wú)量綱參數(shù)ψ=e/l決定網(wǎng)格的畸變程度．圖8給出了溫度的平均相對(duì)誤差隨著畸變程度的變化情況,盡管網(wǎng)格畸變程度較大,該方法仍能得到良好的計(jì)算結(jié)果,這主要得益于修正變分泛函僅涉及邊界積分．從圖9可以看出,采用不同源點(diǎn)布局方式的計(jì)算結(jié)果收斂的趨勢(shì)是相同的,而采用雙重虛擬邊界收斂的速度更快一些．表1和表2分別給出了若干點(diǎn)的溫度和沿著x2方向的熱流分量,盡管個(gè)別點(diǎn)熱流誤差較大,但仍在工程實(shí)際容許范圍內(nèi),并且隨著網(wǎng)格的加密,精度會(huì)逐漸改善．

圖7 網(wǎng)格畸變方案Fig. 7 Mesh distortion schemes

圖8 畸變程度ψ對(duì)計(jì)算精度影響 圖9 網(wǎng)格密度對(duì)計(jì)算精度的影響 Fig. 8 Effects of distortion parameter ψ on Fig. 9 Effects of the mesh density on the computation accuracy the computation accuracy

表1 選擇點(diǎn)處溫度u的計(jì)算結(jié)果

表2 選擇點(diǎn)處熱流分量qx2的計(jì)算結(jié)果

3.2 偏心環(huán)空的熱傳導(dǎo)

為驗(yàn)證本文方法求解曲邊問(wèn)題的有效性,第二個(gè)算例考慮偏心環(huán)空的熱傳導(dǎo)．該算例的尺寸參數(shù)和邊界條件如圖10所示,內(nèi)表面施加Dirichlet邊界條件,外表面施加Neumann邊界條件．在計(jì)算過(guò)程中,整個(gè)求解域離散為78個(gè)4節(jié)點(diǎn)四邊形單元．為便于比較,將采用300個(gè)單元的ABAQUS數(shù)值解作為參考解．

圖10 偏心環(huán)空與有限元網(wǎng)格Fig. 10 An eccentric annulus and the finite element meshes

圖11給出了不同源點(diǎn)布局下雜交基本解有限元法采用78個(gè)單元的溫度等值線圖,可以看出與ABAQUS數(shù)值結(jié)果吻合較好．為了使數(shù)據(jù)可視化,圖12和13給出了不同源點(diǎn)布局下偏心環(huán)空內(nèi)表面和外表面單元形心的周向溫度及其熱流分量qx1和qx2的計(jì)算結(jié)果,與ABAQUS數(shù)值解對(duì)比發(fā)現(xiàn),兩者結(jié)果基本吻合,進(jìn)一步說(shuō)明了本文方法在稀疏網(wǎng)格下仍能得到滿意的結(jié)果．為說(shuō)明該方法的效率,表3列出了CPU計(jì)算所需的時(shí)間,所有的計(jì)算都是在戴爾靈越15-1558(處理器類型:Intel Core i7-5500U)上進(jìn)行的．與ABAQUS相比,HFS-FEM計(jì)算所需的時(shí)間更少,這表明該方法具有更高的效率．

圖11 偏心環(huán)空溫度u的等值線圖 圖12 單元形心處溫度u沿周向的變化Fig. 11 Contour plots of temperature u in Fig. 12 Variations of temperature u at element centroids the eccentric annulus along the circumferential direction

表3 CPU時(shí)間的對(duì)比

圖13 單元形心處熱流分量qx1和qx2沿周向的變化Fig. 13 Variations of heat flux components qx1 and qx2 at element centroids along the circumferential direction

4 結(jié) 論

對(duì)于熱傳導(dǎo)問(wèn)題,本文基于雜交基本解有限元法,分別采用3種源點(diǎn)布局方式(與單元形狀相似的虛擬邊界、圓形虛擬邊界和雙重虛擬邊界)進(jìn)行求解,所得結(jié)論如下:

1) 相比于單元偏置的虛擬邊界和圓形邊界,隨著網(wǎng)格的加密,采用雙重虛擬邊界的收斂性更好．

2) 對(duì)于網(wǎng)格畸變程度較大的算例,圓形虛擬邊界更適合用來(lái)求解．

3) 隨著源點(diǎn)與邊界距離的增加,8源點(diǎn)雙重虛擬邊界下,靠近單元邊界趨于穩(wěn)定更快,并且其矩陣條件數(shù)的上下振蕩次數(shù)較?。畬?duì)于另外兩種方法,盡管矩陣條件數(shù)較小,但其數(shù)值振蕩次數(shù)較大．

本文通過(guò)兩個(gè)數(shù)值算例,與傳統(tǒng)有限元法計(jì)算結(jié)果對(duì)比．驗(yàn)證了3種源點(diǎn)布局方式的可行性,在稀疏網(wǎng)格下,雜交基本解有限元法仍能獲得較高的精度．