二維功能梯度材料熱傳導格點型有限體積法研究

龔京風， 徐宗著， 宣領寬， 江致遠， 鄭文利

(武漢科技大學 汽車與交通工程學院， 湖北 武漢 430065)

功能梯度材料(functionally graded materials，FGM)被廣泛應用于航空航天、核反應堆、內燃機等領域。由于FGM常被應用于高溫高壓的惡劣工作環境，FGM結構內的溫度場、應力場變化劇烈，其熱力性能得到了廣泛關注。劉文光等[1]研究了熱流密度、陶瓷體積分數指數和熱環境等參數對功能梯度圓柱殼熱傳導的影響；Fu等[2]分析了在彈性基礎上，考慮非線性熱環境的多孔功能梯度材料圓柱殼的熱聲響應；許新等[3]基于Eluer-Bernoulli梁理論和單向耦合的熱傳導理論，研究了FGM微梁的熱彈性阻尼；梅靖等[4]研究了考慮厚度方向穩態溫度分布的石墨烯增強功能梯度復合材料板條熱彈性問題。

Steinberg等[5]提出為承受復雜的溫度場，要求材料參數性能在2個甚至3個方向變化，因此，研究二維FGM熱傳導問題是有必要的。Rahul等[6]在一階剪切變形理論基礎上，對二維溫度變化作用下的雙向功能梯度圓板進行了自由振動分析。Tang等[7]研究了具有橫向和縱向位移耦合的雙向功能梯度梁的非線性濕熱動力學。這些研究都考慮多方向功能梯度材料的溫度場影響，研究熱環境下的二維FGM結構性能，準確預測其溫度場是前提。

許楊健等[8]推導了復雜邊界條件下的二維FGM板的穩態溫度場解析解；蔣科[9]在離散域采用數值法，非離散域采用解析法的混合方法，研究了FGM結構的熱傳導特性。解析方法雖然能獲得精確解，但難以應用于復雜工程問題，而數值模擬方法具有更廣泛的適用性。Sladek等[10]采用無網格局部Petrov-Galerkin(meshless local Petrov-Galerkin,MLPG)方法分析了非均勻各向異性FGM三維軸對稱結構的瞬態熱傳導問題；Ching等[11]采用MLPG方法求解了功能梯度梁的瞬態熱彈性變形。MLPG方法計算量較大，處理結構復雜網格量大的熱機械問題存在一定困難。仝國軍等[12]采用有限元方法(finite element methed,FEM)研究了非均勻溫度場下變物性二維FGM板的瞬態熱傳導分布問題。針對FGM熱傳導問題，Gong等[13-14]提出非結構時域有限體積法(finite volume methed,FVM)，該方法將四邊形單元處理為雙線性單元，有效避免由于將四邊形單元作為線性單元處理的數值振蕩問題；隨后進一步發展該方法，提出交錯格點型有限體積法(cell-vetex finite volume methed,CV-FVM)，用于研究活塞、功能梯度多孔材料等熱傳導問題[15]。當前，基于FEM的數值計算方法在固體熱傳導分析領域占據主導地位。考慮到FVM在流體力學里的廣泛應用及其在固體力學中的成功應用，為了發展統一的多物理場數值求解方法，擬基于FVM進一步探究其在FGM熱傳導問題中的應用。

為了研究二維FGM材料熱傳導問題，本文進一步發展交錯CV-FVM，將材料物性作為空間坐標的函數，考慮全域的物性參數變化?；诰€性三角形單元和雙線性四邊形單元，建立二維FGM熱傳導數值求解方法，并驗證其正確性，討論其適用性；同時，在處理復雜FGM熱傳導問題，考慮大溫度梯度區域難以預知情況時，探究CV-FVM在不加密該區域網格情況下的數值穩定性。本文通過研究求解二維FGM問題的CV-FVM解法，為開發擁有自主知識產權的CAE分析軟件奠定基礎。

1 熱傳導數值方法

1.1 控制方程

在控制體M上建立導熱方程，控制體體積為V，邊界面積為S。根據能量守恒定律建立熱傳導方程：

(1)

式中：ρ、c、T分別是密度、比熱容和溫度；q為熱流密度；n為控制體邊界面的外法線矢量；Q為熱源在單位體積內產生的熱量。

由傅里葉定律可知-q=k▽T，則可將式(1)寫成：

(2)

式中k為導熱系數，考慮各項同性問題，k為標量。

1.2 數值離散

1.2.1 控制體的建立

對于二維問題，CV-FVM依次將圍繞節點的相鄰單元中心、邊中點連接，構造控制體。采用三角形單元或四邊形單元離散計算域。圖1為局部網格示意圖，空心點表示單元中心或邊中點，實心點為單元節點。

圖1 二維CV-FVM離散單元示意Fig.1 Schematic diagram of two-dimensional CV-FVM discrete unit

以節點1為例，構造的控制體為圖1中虛線圍成的區域O1-L2-O2-L3-O3-L4-O4-L1-O1，該區域即圍繞節點1的控制體M，其中，L1-O1-L2表示節點1在相鄰單元1內的控制體邊界面，節點2～7皆為節點1的相鄰節點。

采用交錯網格的思想，將待解溫度定義在網格節點，并假設溫度在控制體內均勻分布；將材料參數、體積熱源定義在單元中心，假設其在網格單元內均勻分布。由于材料物性的空間分布，在控制體內，物性參數非均勻，從而將物性參數的空間分布引入離散方程。

1.2.2 導熱方程的離散

為了方便方程的離散，將式(2)改寫為：

(3)

式中nx、ny分別表示控制體邊界面的法向矢量在x、y方向的分量。

方程(3)右端第1項借用FEM的思想，采用形函數對該項進行離散：

(4)

由于體積源項定義在單元中心并假設在單元內均勻分布，因此方程(3)右端第2項可離散為：

(5)

方程(3)左端為一階時間項，采用向后差分的方式離散：

(6)

式中：t為當前時刻;Δt為時間步長。

對于邊界上的節點，相應的導熱方程離散需要考慮恒溫邊界SD，熱流密度邊界SN和換熱邊界SR的影響：

SD:T=TB

(7)

(8)

(9)

式中：TB為邊界溫度;qB為熱流密度;hB為對流換熱系數;Tf為環境溫度。

對于邊界SD上的節點，采用式直接給定溫度。對于邊界SN和邊界SR上的節點，將式(8)和式(9)代入式(4)得：

(10)

將方程的最終離散形式整理為:

(11)

式中：ABi表示與節點相鄰的第i個邊界網格面的面積矢量；nN表示與當前節點相鄰的位于SN上的邊界網格面的個數;nR表示與當前節點相鄰的位于SR上的邊界網格面的個數。

本文只考慮各項同性且無內熱源的穩態問題，因此，式(11)可以簡化為：

(12)

1.3 程序設計

為實現上述數值算法，采用C++語言編寫程序，圖2為程序流程圖。

圖2 程序流程Fig.2 Program flow chart

為了考慮FGM熱傳導問題，在物性參數初始化時，根據其變化規律，基于單元中心坐標計算網格單元內的物性參數。為了考慮非均勻溫度邊界條件，在初始化時，根據邊界面單元中心坐標計算非均勻邊界參數。同時，計算程序不采用迭代法求解，而采用直接法求解。

將離散后的控制方程整理為式的形式：

AX=B

(13)

(14)

(15)

(16)

式中：e表示節點總數；A中的元素φlm為：

(17)

式中：下標l和m皆表示節點編號；P表示當前節點l的相鄰單元所包含節點編號中，存在相鄰節點m的單元數。采用文獻[17]中的置大數法處理恒溫邊界。

2 熱傳導方法驗證

2.1 網格適用性的驗證

計算內徑0.08 m，外徑0.1 m的圓環熱傳導問題。圓環內外兩側皆為恒溫邊界，內側溫度T1=0 ℃，外側溫度T2=1 ℃。圓環的導熱系數沿徑向呈指數函數變化：

(18)

式中：k0=17 W/(m·℃)；λ為分布參數，取λ=50 m-1。

分別采用三角形單元、四邊形單元、混合網格單元劃分計算域，網格尺寸為0.005 m。計算得到的圓環徑向溫度值(見表1)與文獻[18]解析解吻合良好。

表1 不同網格單元在徑向上的溫度值Table 1 The temperature values of different grid cells in the radial direction

圖3為基于不同的網格計算得到的圓環溫度分布云圖，其徑向方向的溫度變化皆與解析解相匹配，說明CV-FVM對于線性三角形單元、雙線性四邊形單元和二者的混合單元皆具有良好的適用性，不存在數值振蕩問題。

圖3 不同網格單元溫度分布云圖Fig.3 Cloud map of temperature distribution in different grid cells

2.2 對二維FGM熱傳導問題的適用性

計算圖4所示的FGM板。正方形板邊長a=0.01 m。上邊界B1給定對流換熱邊界條件，下邊界B2給定熱流密度邊界條件，左右邊界B3、B4給定恒溫邊界：

圖4 功能梯度平板Fig.4 Functional gradient plate

(19)

計算域采用1×10-3m的四邊形網格單元進行劃分。

假設導熱系數沿x和y方向按指數函數規律雙向變化，對流換熱系數沿x方向按指數規律分布，即：

k(x,y)=k0ecx+dy,h(x)=h0ecx

(20)

式中：k0為x=y=0 m處的導熱系數;h0為x=0 m處的換熱系數;c、d為導熱系數梯度分布參數，當c=d=0 m-1時，表示材料參數均勻分布，退化為普通材料。

考慮對流換熱邊界溫度場的解析解為[19]:

T(x,y)=Tw+

(21)

其中：

Kn=4n2π2+a2c2

Ln=1-cos(nπ)eac/2

分別計算2種條件下二維FGM板的熱傳導問題。

算例1下邊界B2絕熱q0=0 W/m2，Tf=100 ℃，Tw=0 ℃，k0=1×103W/(m·℃)，h0=1×106W/(m2·℃)，c=d=200 m-1。將計算得到的x=1×10-3m和x=5×10-3m處的溫度與FEM數值解及解析解對比，如圖5所示。本文計算得到的溫度曲線與FEM解及解析解吻合良好。

圖5 第1種情況結果對比Fig.5 Comparison chart of results in the first case

算例2下邊界B2給定熱流密度q0=1×108W/m2，Tf=800 ℃，Tw=300 ℃，k0=1×103W/(m·℃)，h0=1×106W/(m2·℃)，c=d=100 m-1。分別將基于CV-FVM、FEM計算得到的特征點溫度與文獻中的解析解對比并計算誤差，見表2和表3?？梢钥闯?，相同網格尺寸下，本文CV-FVM計算得到的結果與FEM解的求解精度相當。同時，不同非均勻邊界條件下，基于CV-FVM計算得到的二維FGM板的熱傳導結果都與解析解吻合良好，驗證了本文CV-FVM對不同邊界條件的適用性。

表2 第2種情況CV-FVM結果對比Table 2 Comparison of CV-FVM results in the second case

表3 第2種情況FEM結果對比Table 3 Comparison of FEM results in the second case

為驗證CV-FVM與FEM的計算效率，基于Intel Core i7-10750H CPU，內存為16 GB的計算機，將算例2的網格尺寸加密至2×10-5m，分別采用CV-FVM與FEM再次計算，求解耗時及最大內存消耗見表4。可以看出目前在求解時間上CV-FVM比FEM長，但在內存消耗上要比FEM少，這是由于當前版本的軟件還未在求解速度上進行優化，在后續研究中會逐步改進。

表4 CV-FVM與FEM計算效率對比Table 4 Comparison of computational efficiency between CV-FVM and FEM

2.3 CV-FVM的優勢

基于2.2節的算例1，分別討論c=d=0，-200 m-1時溫度的分布。圖6為c=d=0 m-1時的等溫度線圖。此時導熱系數等于k0，材料為均質材料；換熱系數等于h0，邊界條件均勻。因此，計算得到的溫度關于x=5×10-3m對稱分布。圖7為c=d=-200 m-1時的二維FGM板溫度分布，由于導熱系數沿x，y方向變化，且邊界條件沿x方向變化，等溫線不再對稱分布，上邊界溫度梯度較大?？紤]到板的兩側為T=0 ℃ 的恒溫邊界，上邊界角點處存在劇烈的溫度變化，考察此處數值計算結果，可以體現數值算法求解大梯度問題的性能。

圖6 c=d=0時等溫線Fig.6 Isotherm diagram when c=d=0

圖7 c=d=-200 m-1時等溫線Fig.7 Isotherm diagram when c=d=-200 m-1

取y=1×10-2m處的溫度分布結果做進一步分析。分別對比網格尺寸為0.25×10-3、1×10-3、2×10-3m時的FEM結果和本文CV-FVM結果，并根據公式算出y=1×10-2m處的溫度解析解。可以看到，隨著網格尺寸的增加，FEM結果在角點附近出現了不合理的數值波動現象，而本文CV-FVM始終能給出合理的結果。圖7中這種等溫線匯集在一個點上的現象在物理上不合理，數學上這種位置稱為奇點，這是由角點處存在2個溫度而造成的[20]。在處理這種情況時，由于CV-FVM中將邊界變量存放至邊界網格面中心，使用這些邊界條件變量時只需按照邊界序號依次循環取用，若該處存在2個邊界條件，則由循環順序決定采用哪個邊界。受網格精度影響，當網格尺寸較大，此處的溫度階躍變化表現為一個線性變化(見圖8)。相比于FEM，CV-FVM處理此類問題時，計算所得的結果整體可信度更高，數值更加穩定。在處理復雜熱傳導問題時，難以提前預知大溫度梯度存在的區域，采用本文CV-FVM，即使不加密該區域的網格，也能得到合理的計算結果。

圖8 y=1×10-2 m處的溫度曲線Fig.8 Temperature curve when y=1×10-2 m

3 結論

1) 本文發展的用于二維FGM熱傳導問題的CV-FVM求解方法，能夠處理線性三角形單元、雙線性四邊形單元及混合網格，適用于不同的邊界條件類型。

2) 由于采用交錯網格技術考慮物性參數的空間變化，沒有基于任何分布假設，該方法適用于具有任意材料物性分布的二維FGM熱傳導問題。

3) 在處理邊界角點存在劇烈溫差情況時，FEM需要通過加密網格來避免數值振蕩，而本文發展的CV-FVM即使在粗糙網格情況下，也能給出合理的數值結果。因此，本文發展的CV-FVM更適用于惡劣環境下、物性參數復雜變化的FGM熱傳導問題，尤其對于無法提前預知大溫度梯度區域的問題，本文方法具有明顯優勢。