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

* (北京科技大學數理學院,磁光電復合材料與界面科學北京市重點實驗室,北京 100083)

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

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

引言

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

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

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

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

1.1 熱傳導問題

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

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

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

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

這里,ve為任意的權函數.

進一步考慮到

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

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

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

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

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

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

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

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

圖1 基于不同結點集的Lagrange 插值函數Fig.1 Lagrange interpolation function with different interpolation points

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

其中,域內熱傳導系數矩陣為

由單元內邊界上接觸熱阻引起的與當前單元e及其鄰居單元eb相關的系數矩陣分別為

由單元域內熱源項引起的等效熱載荷為

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

對單元有限元方程(11)進行組裝,即可得到結構溫度場分析的譜元法方程組

引入強制溫度邊界條件即可對方程組(16)進行求解,經過后處理即可得到結構的溫度場和熱流密度場.

1.2 熱應力問題

結構熱應力問題的平衡方程可寫成

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

這里,σ 為應力向量,f為體力向量,u為位移向量,p為邊界力向量.

考慮熱變形的本構方程為

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

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

這里,we為單元位移場分析的權函數.

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

在典型單元內位移場和應變場可分別表示成

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

幾何矩陣為

取權函數we=,并將幾何方程和本構方程代入平衡方程,則位移場分析的單元有限元方程為

其中,域內剛度矩陣為

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

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

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

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

引入強制位移邊界條件即可對方程組(29)進行求解,經過后處理即可得到結構的位移場和應力場.

1.3 譜元法的數值實現

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

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

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

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

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

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

(5)求解式(16)得到結構溫度場,考慮到材料屬性的溫度相關性,反復迭代步驟4～5 得到收斂的溫度場;

(7)基于結構應力場更新界面接觸熱阻,如果接觸熱阻達到收斂準則,則進入下一步,否則回到步驟4,進行新一輪迭代計算;

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

2 數值算例

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

2.1 算例1

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

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

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

圖3 結果表明對于一維熱傳導問題,基于譜單元方法結果的收斂速度和計算精度都要高于經典線性單元,在總求解自由度數比較小的情況下,插值階次越高則精度越好.而當總自由度數比較大時,溫度場的2-范數誤差并不總是最小的,這是因為自由度較高時基于各種插值階次的計算精度都比較高,此時數值誤差占主導,而插值函數和插值方法的誤差處于次要地位.

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

2.2 算例2

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

圖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 可以看出,無論是定性比較還是定量比較,基于譜元法得到的數值解都和精確解吻合得很好.圖6 進一步給出了采用等距結點的經典高階有限元與采用Lobatto 和Chebyshev 非等距結點集的譜單元法得到的數值結果的收斂速度對比,從中可以明顯看出,基于三種不同結點集的計算結果的收斂速度相差不大,但是在相同自由度的條件下,基于非等距結點集譜元法的精度更高,而且基于Lobatto 結點集的結果精度比基于Chebyshev 結點集的精度還要更好一些.

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

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

2.3 算例3

兩塊長度均為1 m 且寬度均為0.2 m 的矩形板沿長度方向對接在一起,組合結構的左端給定溫度273 K,右端施加熱流1000 W/m2.左右兩端均固定.材料為45#鋼,其彈性模量、泊松比、熱膨脹系數和熱導率均隨溫度變化[27].兩塊板的接觸界面存在應力相關的接觸熱阻R=A0(σ/σ0)-2.5,其中A0=0.01 m2K/W,σ0=70 MPa,σ 為界面法向接觸擠壓應力,單位為MPa.

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

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

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

圖9 基于三階單元的結構溫度場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 方向應力場Fig.11 Numerical solutions of stress field along x-coordinate with third order element

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

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

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

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

圖12 結構溫度場Fig.12 Temperature field

圖13 結構x 方向位移場Fig.13 Displacement field along x-coordinate

圖14 結構y 方向位移場Fig.14 Displacement field along y-coordinate

圖15 結構x 方向應力場Fig.15 Stress field along x-coordinate

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

3 結論

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

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

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

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

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