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

龔京風(fēng)， 徐宗著， 宣領(lǐng)寬， 江致遠(yuǎn)， 鄭文利

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

功能梯度材料(functionally graded materials，F(xiàn)GM)被廣泛應(yīng)用于航空航天、核反應(yīng)堆、內(nèi)燃機(jī)等領(lǐng)域。由于FGM常被應(yīng)用于高溫高壓的惡劣工作環(huán)境，F(xiàn)GM結(jié)構(gòu)內(nèi)的溫度場、應(yīng)力場變化劇烈，其熱力性能得到了廣泛關(guān)注。劉文光等[1]研究了熱流密度、陶瓷體積分?jǐn)?shù)指數(shù)和熱環(huán)境等參數(shù)對功能梯度圓柱殼熱傳導(dǎo)的影響；Fu等[2]分析了在彈性基礎(chǔ)上，考慮非線性熱環(huán)境的多孔功能梯度材料圓柱殼的熱聲響應(yīng)；許新等[3]基于Eluer-Bernoulli梁理論和單向耦合的熱傳導(dǎo)理論，研究了FGM微梁的熱彈性阻尼；梅靖等[4]研究了考慮厚度方向穩(wěn)態(tài)溫度分布的石墨烯增強(qiáng)功能梯度復(fù)合材料板條熱彈性問題。

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

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

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

1 熱傳導(dǎo)數(shù)值方法

1.1 控制方程

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

(1)

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

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

(2)

式中k為導(dǎo)熱系數(shù)，考慮各項同性問題，k為標(biāo)量。

1.2 數(shù)值離散

1.2.1 控制體的建立

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

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

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

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

1.2.2 導(dǎo)熱方程的離散

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

(3)

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

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

(4)

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

(5)

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

(6)

式中：t為當(dāng)前時刻;Δt為時間步長。

對于邊界上的節(jié)點，相應(yīng)的導(dǎo)熱方程離散需要考慮恒溫邊界SD，熱流密度邊界SN和換熱邊界SR的影響：

SD:T=TB

(7)

(8)

(9)

式中：TB為邊界溫度;qB為熱流密度;hB為對流換熱系數(shù);Tf為環(huán)境溫度。

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

(10)

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

(11)

式中：ABi表示與節(jié)點相鄰的第i個邊界網(wǎng)格面的面積矢量；nN表示與當(dāng)前節(jié)點相鄰的位于SN上的邊界網(wǎng)格面的個數(shù);nR表示與當(dāng)前節(jié)點相鄰的位于SR上的邊界網(wǎng)格面的個數(shù)。

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

(12)

1.3 程序設(shè)計

為實現(xiàn)上述數(shù)值算法，采用C++語言編寫程序，圖2為程序流程圖。

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

為了考慮FGM熱傳導(dǎo)問題，在物性參數(shù)初始化時，根據(jù)其變化規(guī)律，基于單元中心坐標(biāo)計算網(wǎng)格單元內(nèi)的物性參數(shù)。為了考慮非均勻溫度邊界條件，在初始化時，根據(jù)邊界面單元中心坐標(biāo)計算非均勻邊界參數(shù)。同時，計算程序不采用迭代法求解，而采用直接法求解。

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

AX=B

(13)

(14)

(15)

(16)

式中：e表示節(jié)點總數(shù)；A中的元素φlm為：

(17)

式中：下標(biāo)l和m皆表示節(jié)點編號；P表示當(dāng)前節(jié)點l的相鄰單元所包含節(jié)點編號中，存在相鄰節(jié)點m的單元數(shù)。采用文獻(xiàn)[17]中的置大數(shù)法處理恒溫邊界。

2 熱傳導(dǎo)方法驗證

2.1 網(wǎng)格適用性的驗證

計算內(nèi)徑0.08 m，外徑0.1 m的圓環(huán)熱傳導(dǎo)問題。圓環(huán)內(nèi)外兩側(cè)皆為恒溫邊界，內(nèi)側(cè)溫度T1=0 ℃，外側(cè)溫度T2=1 ℃。圓環(huán)的導(dǎo)熱系數(shù)沿徑向呈指數(shù)函數(shù)變化：

(18)

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

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

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

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

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

2.2 對二維FGM熱傳導(dǎo)問題的適用性

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

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

(19)

計算域采用1×10-3m的四邊形網(wǎng)格單元進(jìn)行劃分。

假設(shè)導(dǎo)熱系數(shù)沿x和y方向按指數(shù)函數(shù)規(guī)律雙向變化，對流換熱系數(shù)沿x方向按指數(shù)規(guī)律分布，即：

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

(20)

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

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

T(x,y)=Tw+

(21)

其中：

Kn=4n2π2+a2c2

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

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

算例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數(shù)值解及解析解對比，如圖5所示。本文計算得到的溫度曲線與FEM解及解析解吻合良好。

圖5 第1種情況結(jié)果對比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計算得到的特征點溫度與文獻(xiàn)中的解析解對比并計算誤差，見表2和表3。可以看出，相同網(wǎng)格尺寸下，本文CV-FVM計算得到的結(jié)果與FEM解的求解精度相當(dāng)。同時，不同非均勻邊界條件下，基于CV-FVM計算得到的二維FGM板的熱傳導(dǎo)結(jié)果都與解析解吻合良好，驗證了本文CV-FVM對不同邊界條件的適用性。

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

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

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

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

2.3 CV-FVM的優(yōu)勢

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

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

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

3 結(jié)論

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

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

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