含高溫度梯度及接觸熱阻非線性熱力耦合問題的譜元法1)

史 鑫 趙建寧 楊苗苗 張家雷 劉冬歡 ,2)

* (北京科技大學(xué)數(shù)理學(xué)院,磁光電復(fù)合材料與界面科學(xué)北京市重點(diǎn)實(shí)驗(yàn)室,北京 100083)

? (北京科技大學(xué)能源與環(huán)境工程學(xué)院,北京 100083)

** (中國工程物理研究院流體物理研究所,四川綿陽 621900)

引言

對于高超音速飛行器來說,由于急劇的氣動加熱在駐點(diǎn)區(qū)域會產(chǎn)生極高的溫度梯度,并由此產(chǎn)生非線性熱力耦合問題.求解此類問題的傳統(tǒng)有限元方法往往采用低階單元,因此需要在駐點(diǎn)區(qū)域采用相當(dāng)密集的網(wǎng)格才能得到可信的結(jié)果[1].采用廣義有限元方法是求解此類特殊問題的一個選擇,O’Hara等[2-4]研究了廣義有限元方法在存在局部高溫度梯度的瞬態(tài)傳熱問題上的應(yīng)用,在相對粗糙的網(wǎng)格上通過在局部區(qū)域附加特殊插值函數(shù)的方式來獲得高精度解.文獻(xiàn)[5-6]建立了局部加熱引起的高溫度梯度及高強(qiáng)應(yīng)力應(yīng)變問題的廣義有限元方法,通過局部加密網(wǎng)格和插值函數(shù)升階的方法來捕捉場變量的局部高梯度效應(yīng).Tang 等[7]提出了一種廣義有限元方法來求解熱傳遞及熱彈性問題,通過采用高階的插值函數(shù)來保證單元邊界上溫度梯度的連續(xù)性.Sanchez 和Duarte[8]建立了高階廣義有限元方法并將其用于多尺度結(jié)構(gòu)動力學(xué)和波動問題的求解.Kim 等[9]很早就指出廣義有限元方法潛在的缺點(diǎn)是如果局部網(wǎng)格加密或者插值函數(shù)升階不當(dāng),那么會帶來更大誤差.Calabrò等[10]提出了一種求解高梯度橢圓方程的機(jī)器學(xué)習(xí)策略.Peng[11]針對高梯度問題提出了一種基于最小二乘法的新型擴(kuò)展有限元方法,可以有效提高數(shù)值解的精度但是其收斂性會下降.對于此類高溫度梯度問題也可以直接采用高階p 型有限元方法[12].p 型有限元方法一般利用Lagrange多項(xiàng)式作為插值函數(shù),且插值結(jié)點(diǎn)是等距分布的,對于高階插值存在Runge現(xiàn)象.而譜方法[13]的插值函數(shù)采用插值結(jié)點(diǎn)非等距分布的正交多項(xiàng)式,可以很好地避免Runge 現(xiàn)象,在波傳播等周期性問題中得到了廣泛的應(yīng)用[14].譜方法的一個缺點(diǎn)是只能在規(guī)則區(qū)域上進(jìn)行離散,限制了使用范圍,而實(shí)際工程結(jié)構(gòu)的求解域往往是非規(guī)則的.因此,將譜方法與有限元方法相結(jié)合而成的譜有限元法,既保留了譜方法的高精度與快速收斂性又具有有限元方法可以求解任意不規(guī)則區(qū)域的優(yōu)異適應(yīng)性,受到了研究者的廣泛關(guān)注[15-16].譜有限元法的插值函數(shù)可以基于正交多項(xiàng)式,比如Legendre 多項(xiàng)式或者Chebyshev 多項(xiàng)式,也可以基于非等距結(jié)點(diǎn)的Lagrange 多項(xiàng)式,非等距的結(jié)點(diǎn)位置往往可以選取各種正交多項(xiàng)式的零點(diǎn).Komatitsch[17]分別采用四邊形和三角形譜單元分析了彈性波在二維結(jié)構(gòu)中的傳播行為.Kudela 等[18]基于Lagrange 多項(xiàng)式建立了一維桿、梁和二維平面譜單元的計算格式,并將數(shù)值結(jié)果與有限元計算結(jié)果和實(shí)驗(yàn)測試結(jié)果比較,驗(yàn)證了譜單元方法的高效性和高精度.Liu[19]基于Chebyshev 多項(xiàng)式構(gòu)造了一種Mindlin 偽譜單元,并對板結(jié)構(gòu)進(jìn)行了靜力、動力和波傳播分析.譜元法在其他工程領(lǐng)域也得到了廣泛的應(yīng)用,比如Merzari[20]采用基于Gauss-Lobatto-Legendre 結(jié)點(diǎn)集的Lagrange 多項(xiàng)式的譜元法研究了核反應(yīng)堆堆芯內(nèi)部的渦流問題.Xu 等[21]基于第一類Chebyshev多項(xiàng)式的譜元法研究了功能梯度材料結(jié)構(gòu)熱彈性性能.劉玉柱[22]利用譜元法進(jìn)行全波形反演,提高了反演算法的計算精度與效率.Zakian 和Bathe[23]采用基于Lobatto 結(jié)點(diǎn)集的Lagrange 多項(xiàng)式的譜元法研究了結(jié)構(gòu)動力學(xué)問題,結(jié)合質(zhì)量陣可對角化的特性以及時間積分的Noh-Bathe 格式,大大提高了計算精度和效率.

值得注意的是,對于工程結(jié)構(gòu)中經(jīng)常遇到的存在高溫度梯度現(xiàn)象,此前的研究已經(jīng)證明了高階單元在捕捉局部高梯度現(xiàn)象中的優(yōu)勢[1],與此同時,固體界面的接觸熱阻也會導(dǎo)致界面上溫度的跳變,這可以看作是一種廣義的高溫度梯度.本文建立可用于求解含高溫度梯度及接觸熱阻的結(jié)構(gòu)非線性熱力耦合問題的譜元法計算格式,并通過數(shù)值算例驗(yàn)證本文方法和程序的有效性及精度,為解決此類非線性熱力耦合問題提供了一種新的高效高精度方法.

1 熱力耦合問題的譜元法格式

這里采用順序耦合的方法建立熱力耦合問題的譜元法格式,分別建立熱傳導(dǎo)問題和熱應(yīng)力問題的計算格式,通過數(shù)值迭代的方法得到熱力耦合問題的解.

1.1 熱傳導(dǎo)問題

邊界為 ?Ω的區(qū)域 Ω 上的熱傳導(dǎo)問題的控制方程為

其中,? 為梯度算子,k是對稱的熱傳導(dǎo)系數(shù)矩陣,T為溫度,q0是熱源函數(shù).問題的Dirichlet 邊界條件和Neumann 邊界條件分別為

其中,?ΩD和 ?ΩN分別為Dirichlet 邊界和Neumann邊界,n為邊界 ?Ω的單位外法線向量,q為熱流密度向量,分別為給定的溫度值和熱流密度值.

在單元域內(nèi)方程(1)的等效積分弱形式為

這里,ve為任意的權(quán)函數(shù).

進(jìn)一步考慮到

其中,熱流密度為qe=-ke?Te.將式(4)代入式(3)可得原問題的等效積分形式為

其中,ne為單元e的邊界 ?Ωe上的外法線.

方程(5)中出現(xiàn)了單元邊界 ?Ωe上的積分項(xiàng),如果該單元邊界是結(jié)構(gòu)的外邊界,那么熱流可以由給定熱流、對流換熱或輻射換熱等方式給出.而對于不存在接觸熱阻的單元內(nèi)邊界,該積分項(xiàng)會在組裝總體有限元方程組時自動消除,此時熱量平衡方程的等效積分形式為

對于存在接觸熱阻R的單元內(nèi)邊界,采用間斷有限元方法的處理思想[24],將接觸界面處的熱流密度由接觸熱阻的定義式給出

其中Te和Teb分別是界面兩側(cè)的溫度.

與常規(guī)有限元方法一樣,譜元法需要對單元內(nèi)的場變量進(jìn)行插值近似

本文采用Lagrange 型插值函數(shù),其m階一維形式為

傳統(tǒng)的一維Lagrange 插值函數(shù)的插值結(jié)點(diǎn)在其插值區(qū)間上是等距均布的,而譜單元插值函數(shù)的結(jié)點(diǎn)是非等距分布的,這些非等距結(jié)點(diǎn)一般為Lobatto結(jié)點(diǎn)集或第二類Chebyshev 結(jié)點(diǎn)集[14].Lobatto 結(jié)點(diǎn)集即為Lobatto 多項(xiàng)式的零點(diǎn),Lobatto 多項(xiàng)式是對Legendre 多項(xiàng)式進(jìn)行求導(dǎo)得到的.第二類Chebyshev結(jié)點(diǎn)集即第二類Chebyshev 多項(xiàng)式的零點(diǎn),其m+1個結(jié)點(diǎn)的坐標(biāo)為:ξi=cos[(i-1)π/m] .基于不同結(jié)點(diǎn)集的前五階Lagrange 插值函數(shù)圖像如圖1 所示.值得注意的是,也可以直接采用各種正交多項(xiàng)式比如Legendre 多項(xiàng)式或Chebyshev 多項(xiàng)式作為插值函數(shù).

圖1 基于不同結(jié)點(diǎn)集的Lagrange 插值函數(shù)Fig.1 Lagrange interpolation function with different interpolation points

二維Lagrange 插值函數(shù)可以利用ξ 和η 方向上的一維插值函數(shù) φ(ξ)與 φ(η)的張量積得到,例如二階二維Lagrange 插值函數(shù)可寫為

其中,域內(nèi)熱傳導(dǎo)系數(shù)矩陣為

由單元內(nèi)邊界上接觸熱阻引起的與當(dāng)前單元e及其鄰居單元eb相關(guān)的系數(shù)矩陣分別為

由單元域內(nèi)熱源項(xiàng)引起的等效熱載荷為

由給定熱流qn的邊界條件引起的等效載荷列向量為

對單元有限元方程(11)進(jìn)行組裝,即可得到結(jié)構(gòu)溫度場分析的譜元法方程組

引入強(qiáng)制溫度邊界條件即可對方程組(16)進(jìn)行求解,經(jīng)過后處理即可得到結(jié)構(gòu)的溫度場和熱流密度場.

1.2 熱應(yīng)力問題

結(jié)構(gòu)熱應(yīng)力問題的平衡方程可寫成

問題的Dirichlet 邊界條件和Neumann 邊界條件分別為

這里,σ 為應(yīng)力向量,f為體力向量,u為位移向量,p為邊界力向量.

考慮熱變形的本構(gòu)方程為

其中,De是彈性矩陣,T0是計算熱變形的參考溫度,εe為應(yīng)變向量,αe是熱膨脹系數(shù)向量.

與推導(dǎo)溫度場控制方程的等效弱形式類似,力平衡方程的等效積分弱形式可寫為

這里,we為單元位移場分析的權(quán)函數(shù).

與溫度場分析時在單元接觸界面需要引入表征非完全熱接觸的接觸熱阻相對應(yīng),應(yīng)力場分析時需要在接觸界面上引入位移的不可穿透條件.這里采用共節(jié)點(diǎn)的處理方法.

在典型單元內(nèi)位移場和應(yīng)變場可分別表示成

位移場插值函數(shù)的推導(dǎo)與溫度場類似,此處不再贅述.與溫度場不同的是,對于平面二維問題,單元節(jié)點(diǎn)位移有兩個基本未知量,因此位移場的插值函數(shù)為

幾何矩陣為

取權(quán)函數(shù)we=,并將幾何方程和本構(gòu)方程代入平衡方程,則位移場分析的單元有限元方程為

其中,域內(nèi)剛度矩陣為

由體力引起的載荷列向量為

由溫度場引起的載荷列向量為

由力邊界條件引起的載荷列向量

對單元有限元方程(24)進(jìn)行組裝,即可得到位移場分析的譜元法方程組

引入強(qiáng)制位移邊界條件即可對方程組(29)進(jìn)行求解,經(jīng)過后處理即可得到結(jié)構(gòu)的位移場和應(yīng)力場.

1.3 譜元法的數(shù)值實(shí)現(xiàn)

以上系數(shù)矩陣的數(shù)值計算一般采用Gauss 積分法[25].對于采用正交多項(xiàng)式作為插值函數(shù)的譜元法,其系數(shù)矩陣計算時也可以利用插值函數(shù)的正交性而采用顯式解析積分方案,有時能得到對角化的質(zhì)量陣,對于瞬態(tài)問題來說,這是一個非常有用的特性[26].對于不同的正交多項(xiàng)式插值函數(shù)可以選擇不同的積分方法,如果選擇Lobatto 多項(xiàng)式作為插值函數(shù),Gauss-Lobatto-Legendre 積分方法計算效率最高,而如果選擇Chebyshev 多項(xiàng)式作為插值函數(shù),Gauss-Legendre 積分方法計算效率最高.

本文對于存在應(yīng)力相關(guān)接觸熱阻的熱力耦合問題采用順序耦合的方法進(jìn)行迭代求解,即在初始假設(shè)的接觸熱阻條件下進(jìn)行溫度場分析,考慮到材料熱導(dǎo)率的溫度相關(guān)性,在溫度場分析時也需要采用直接迭代法進(jìn)行迭代求解,得到結(jié)構(gòu)溫度場之后進(jìn)行結(jié)構(gòu)位移場應(yīng)力場分析,得到結(jié)構(gòu)應(yīng)力場之后對接觸熱阻和接觸狀態(tài)進(jìn)行更新,然后再次進(jìn)行溫度場分析,待溫度場和位移場同時滿足收斂條件時迭代終止.整個的計算流程如下:

(1)輸入幾何、材料、載荷、邊界條件及控制參數(shù);

(2)初始化界面接觸熱阻R;

(3)基于設(shè)定的插值函數(shù)階次,利用式(9)的張量積形式確定溫度場插值函數(shù);

(4)基于式(12)～式(15)得到單元有限元方程的系數(shù)矩陣,根據(jù)常規(guī)有限元方法組裝總體有限元方程的方法,得到總體譜元法方程式(16);

(5)求解式(16)得到結(jié)構(gòu)溫度場,考慮到材料屬性的溫度相關(guān)性,反復(fù)迭代步驟4～5 得到收斂的溫度場;

(7)基于結(jié)構(gòu)應(yīng)力場更新界面接觸熱阻,如果接觸熱阻達(dá)到收斂準(zhǔn)則,則進(jìn)入下一步,否則回到步驟4,進(jìn)行新一輪迭代計算;

(8)輸出結(jié)果,計算結(jié)束.

2 數(shù)值算例

本節(jié)給出三個數(shù)值算例來驗(yàn)證本文所建立的求解含高溫度梯度及接觸熱阻的非線性熱力耦合問題的譜元法的可行性及算法精度.

2.1 算例1

圖2 一維熱傳遞問題譜元法數(shù)值解與精確解的對比Fig.2 Comparisons of numerical results with spectral element method and exact solution for one-dimensional heat transfer problem

從圖2 可以看出,譜單元可以用較少的單元數(shù)得到與精確解一致的結(jié)果,且能充分描述域內(nèi)溫度的劇烈變化.與此同時,和相同總自由度數(shù)目的經(jīng)典線性單元相比,在捕捉溫度峰值等細(xì)節(jié)上譜單元的精度更高.圖3 進(jìn)一步給出了不同階次譜單元結(jié)果的收斂速度,這里采用的溫度場2-范數(shù)誤差定義為

其中,TSEM為譜元法計算得到的節(jié)點(diǎn)溫度列向量,Texact為相同維度的節(jié)點(diǎn)溫度列向量精確解.

圖3 結(jié)果表明對于一維熱傳導(dǎo)問題,基于譜單元方法結(jié)果的收斂速度和計算精度都要高于經(jīng)典線性單元,在總求解自由度數(shù)比較小的情況下,插值階次越高則精度越好.而當(dāng)總自由度數(shù)比較大時,溫度場的2-范數(shù)誤差并不總是最小的,這是因?yàn)樽杂啥容^高時基于各種插值階次的計算精度都比較高,此時數(shù)值誤差占主導(dǎo),而插值函數(shù)和插值方法的誤差處于次要地位.

圖3 基于不同插值階次譜單元的計算結(jié)果的收斂速度Fig.3 Convergency rate of numerical results with spectral element of different interpolation orders

2.2 算例2

考慮正方形域上的二維傳熱問題,正方形邊長為6 m,其材料熱導(dǎo)率為1 W/mK,四條邊上保持恒定溫度為0 K.其溫度場的精確解為T(x,y)=,相應(yīng)地,內(nèi)生熱率可由溫度場的精確解求得.圖4 給出了基于Lobatto 結(jié)點(diǎn)集的五階譜單元得到的結(jié)構(gòu)溫度場與精確解的對比.圖5 進(jìn)一步給出了沿著x=1 和y=0 兩條路徑上溫度分布的數(shù)值解與精確解的對比.

圖4 二維熱傳遞問題基于譜元法的溫度場與精確解的對比Fig.4 Comparisons of numerical results with spectral element method and exact solution of the temperature field for two-dimensional heat transfer problem

從圖4 和圖5 可以看出,無論是定性比較還是定量比較,基于譜元法得到的數(shù)值解都和精確解吻合得很好.圖6 進(jìn)一步給出了采用等距結(jié)點(diǎn)的經(jīng)典高階有限元與采用Lobatto 和Chebyshev 非等距結(jié)點(diǎn)集的譜單元法得到的數(shù)值結(jié)果的收斂速度對比,從中可以明顯看出,基于三種不同結(jié)點(diǎn)集的計算結(jié)果的收斂速度相差不大,但是在相同自由度的條件下,基于非等距結(jié)點(diǎn)集譜元法的精度更高,而且基于Lobatto 結(jié)點(diǎn)集的結(jié)果精度比基于Chebyshev 結(jié)點(diǎn)集的精度還要更好一些.

圖5 沿x=1 和y=0 兩條路徑溫度分布的數(shù)值解與精確解對比Fig.5 Comparisons of numerical results with spectral element method and exact solution for temperature distribution along x=1 and y=0

圖6 基于不同結(jié)點(diǎn)集譜單元法的計算結(jié)果的收斂速度Fig.6 Convergency rate of numerical results with spectral element with different interpolation points

2.3 算例3

兩塊長度均為1 m 且寬度均為0.2 m 的矩形板沿長度方向?qū)釉谝黄?組合結(jié)構(gòu)的左端給定溫度273 K,右端施加熱流1000 W/m2.左右兩端均固定.材料為45#鋼,其彈性模量、泊松比、熱膨脹系數(shù)和熱導(dǎo)率均隨溫度變化[27].兩塊板的接觸界面存在應(yīng)力相關(guān)的接觸熱阻R=A0(σ/σ0)-2.5,其中A0=0.01 m2K/W,σ0=70 MPa,σ 為界面法向接觸擠壓應(yīng)力,單位為MPa.

基于8 個矩形和四邊形譜單元(有限元網(wǎng)格參見圖7)計算得到的沿結(jié)構(gòu)長度方向y=0 上的溫度分布和x方向的位移分布與參考解的對比如圖8 所示.這里參考解基于ANSYS 軟件,采用了經(jīng)過收斂性測試的100 000 個四節(jié)點(diǎn)平面應(yīng)力單元.全場的溫度分布和沿x方向位移、應(yīng)力云圖分別如圖9～圖11 所示.

圖7 矩形和四邊形譜單元網(wǎng)格劃分Fig.7 Spectral element with rectangular and quadrilateral mesh

圖8 基于譜單元的計算結(jié)果與參考解的對比Fig.8 Comparisons of numerical solutions of spectral element with reference solution

圖9 基于三階單元的結(jié)構(gòu)溫度場Fig.9 Numerical solutions of temperature field with third order element

圖10 基于三階單元的x 方向位移場Fig.10 Numerical solutions of displacement field along x-coordinate with third order element

圖11 基于三階譜單元的x 方向應(yīng)力場Fig.11 Numerical solutions of stress field along x-coordinate with third order element

可以看出,無論是溫度場還是水平方向的位移場和應(yīng)力場,本文基于譜元法得到的計算結(jié)果都與參考解吻合得很好,由于兩板接觸界面上存在應(yīng)力相關(guān)的接觸熱阻,因此這是一個非線性問題,需要進(jìn)行迭代求解,由于接觸熱阻的存在,使得界面左右兩側(cè)出現(xiàn)了明顯的溫度跳變,如圖8(a)所示.與此同時,無論是矩形譜單元還是任意四邊形譜單元,計算結(jié)果都與參考解吻合得非常好,這表明本文方法在求解非線性熱力耦合問題時具有良好的精度,相比傳統(tǒng)譜方法只能在矩形單元上求解,譜元法很好地融合了譜方法和有限元方法的優(yōu)點(diǎn),可以有效地在不規(guī)則網(wǎng)格上求解,因此譜單元可以很好地應(yīng)用在幾何形狀復(fù)雜的實(shí)際工程問題中.

數(shù)值計算時基于臺式電腦Intel (R)Core (TM)i7-10750H CPU@2.60 GHz,由于對非線性問題需要迭代求解,因此計算時間較長.從表1 可以看出,無論是溫度場還是位移場的相對誤差,同樣的總自由度(104)時高階譜單元結(jié)果比經(jīng)典線性單元精度更高用時更少,甚至高階譜單元用更少的總自由度(104和170)得到了比線性單元采用更多自由度(210)更高的計算精度.與此同時,自由度總數(shù)為170 的高階譜單元方法計算用時25.144 s,而自由度總數(shù)為104 的經(jīng)典線性單元則用時61.557 s.這充分表明,本文提出的高階譜單元方法,不僅利用較少的總自由度數(shù)即可得到高精度的計算結(jié)果,而且計算時間比耗費(fèi)更多自由度的經(jīng)典線性單元還要少,也即高階譜單元兼有高精度和高效率的特點(diǎn),在存在高溫度梯度及接觸熱阻的非線性熱力耦合問題上具有廣闊的應(yīng)用前景.

表1 不同網(wǎng)格及插值階次譜元法的計算精度及計算時間Table 1 Computational accuracy and CPU time of spectral element method with different mesh grid and interpolation orders

進(jìn)一步為了展示本文譜方法在三維熱力耦合問題上的應(yīng)用,保持本算例結(jié)構(gòu)的長度和寬度不變,將厚度方向上的尺寸取為和結(jié)構(gòu)寬度一樣,材料屬性、邊界條件、載荷情況和界面參數(shù)等也保持不變.基于本文譜方法得到的計算結(jié)果與商用軟件ANSYS 結(jié)果的對比如圖12～圖15 所示.

圖12 結(jié)構(gòu)溫度場Fig.12 Temperature field

圖13 結(jié)構(gòu)x 方向位移場Fig.13 Displacement field along x-coordinate

圖14 結(jié)構(gòu)y 方向位移場Fig.14 Displacement field along y-coordinate

圖15 結(jié)構(gòu)x 方向應(yīng)力場Fig.15 Stress field along x-coordinate

可以看出,針對三維熱力耦合問題本文譜方法得到的結(jié)果和商用軟件得到的結(jié)果吻合得很好,從而驗(yàn)證了本文譜方法在三維問題上的可行性和精度.

3 結(jié)論

基于Lobatto 及Chebyshev 非等距插值結(jié)點(diǎn)建立了Lagrange 型譜單元插值函數(shù),并將其應(yīng)用于存在高溫度梯度及接觸熱阻的非線性熱力耦合問題的求解,建立了相關(guān)譜元法的計算格式及計算流程,通過數(shù)值算例驗(yàn)證了本文建立的譜元法的計算精度及計算效率.數(shù)值結(jié)果表明:

(1)基于非等距結(jié)點(diǎn)的譜元法可以用較少的自由度捕捉到結(jié)構(gòu)域內(nèi)溫度劇烈變化的趨勢,計算精度和效率都比經(jīng)典線性有限元法高;

(2)任意四邊形單元的譜元法適用于幾何形狀復(fù)雜的實(shí)際工程問題,克服了傳統(tǒng)譜方法只能使用矩形單元的弱點(diǎn),適應(yīng)性更強(qiáng);

(3)對于含高溫度梯度及接觸熱阻的非線性熱力耦合問題,譜元法可以高效高精度地給出結(jié)構(gòu)溫度場及熱變形熱應(yīng)力,具有廣闊的工程應(yīng)用前景.

值得注意的是,在本文研究中,雖然對場變量的插值函數(shù)采用了高階格式,但是對單元幾何場,采用了場變量插值函數(shù)一致的插值函數(shù),存在一定的幾何離散誤差,下一步對幾何場可以嘗試采用非均勻有理NURBS 張量積為基函數(shù)進(jìn)行近似[28],有望進(jìn)一步將譜元法的優(yōu)勢發(fā)展到極致.