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

龔京風(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)的溫度場(chǎng)、應(yīng)力場(chǎng)變化劇烈，其熱力性能得到了廣泛關(guān)注。劉文光等[1]研究了熱流密度、陶瓷體積分?jǐn)?shù)指數(shù)和熱環(huán)境等參數(shù)對(duì)功能梯度圓柱殼熱傳導(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ù)雜的溫度場(chǎng)，要求材料參數(shù)性能在2個(gè)甚至3個(gè)方向變化，因此，研究二維FGM熱傳導(dǎo)問題是有必要的。Rahul等[6]在一階剪切變形理論基礎(chǔ)上，對(duì)二維溫度變化作用下的雙向功能梯度圓板進(jìn)行了自由振動(dòng)分析。Tang等[7]研究了具有橫向和縱向位移耦合的雙向功能梯度梁的非線性濕熱動(dòng)力學(xué)。這些研究都考慮多方向功能梯度材料的溫度場(chǎng)影響，研究熱環(huán)境下的二維FGM結(jié)構(gòu)性能，準(zhǔn)確預(yù)測(cè)其溫度場(chǎng)是前提。

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

為了研究二維FGM材料熱傳導(dǎo)問題，本文進(jìn)一步發(fā)展交錯(cuò)CV-FVM，將材料物性作為空間坐標(biāo)的函數(shù)，考慮全域的物性參數(shù)變化?；诰€性三角形單元和雙線性四邊形單元，建立二維FGM熱傳導(dǎo)數(shù)值求解方法，并驗(yàn)證其正確性，討論其適用性；同時(shí)，在處理復(fù)雜FGM熱傳導(dǎo)問題，考慮大溫度梯度區(qū)域難以預(yù)知情況時(shí)，探究CV-FVM在不加密該區(qū)域網(wǎng)格情況下的數(shù)值穩(wěn)定性。本文通過研究求解二維FGM問題的CV-FVM解法，為開發(fā)擁有自主知識(shí)產(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ù)，考慮各項(xiàng)同性問題，k為標(biāo)量。

1.2 數(shù)值離散

1.2.1 控制體的建立

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

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

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

采用交錯(cuò)網(wǎng)格的思想，將待解溫度定義在網(wǎng)格節(jié)點(diǎn)，并假設(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項(xiàng)借用FEM的思想，采用形函數(shù)對(duì)該項(xiàng)進(jìn)行離散：

(4)

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

(5)

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

(6)

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

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

SD:T=TB

(7)

(8)

(9)

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

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

(10)

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

(11)

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

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

(12)

1.3 程序設(shè)計(jì)

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

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

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

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

AX=B

(13)

(14)

(15)

(16)

式中：e表示節(jié)點(diǎn)總數(shù)；A中的元素φl(shuí)m為：

(17)

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

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

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

計(jì)算內(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)格單元?jiǎng)澐钟?jì)算域，網(wǎng)格尺寸為0.005 m。計(jì)算得到的圓環(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)格計(jì)算得到的圓環(huán)溫度分布云圖，其徑向方向的溫度變化皆與解析解相匹配，說明CV-FVM對(duì)于線性三角形單元、雙線性四邊形單元和二者的混合單元皆具有良好的適用性，不存在數(shù)值振蕩問題。

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

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

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

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

(19)

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

假設(shè)導(dǎo)熱系數(shù)沿x和y方向按指數(shù)函數(shù)規(guī)律雙向變化，對(duì)流換熱系數(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í)，表示材料參數(shù)均勻分布，退化為普通材料。

考慮對(duì)流換熱邊界溫度場(chǎng)的解析解為[19]:

T(x,y)=Tw+

(21)

其中：

Kn=4n2π2+a2c2

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

分別計(jì)算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。將計(jì)算得到的x=1×10-3m和x=5×10-3m處的溫度與FEM數(shù)值解及解析解對(duì)比，如圖5所示。本文計(jì)算得到的溫度曲線與FEM解及解析解吻合良好。

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

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

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

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

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

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

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

圖6 c=d=0時(shí)等溫線Fig.6 Isotherm diagram when c=d=0

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

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

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