導熱-輻射耦合傳熱的多尺度分析和數值模型

童自翔，李明佳，李冬

1. 西安交通大學 人居環境與建筑工程學院，西安 710049

2. 西安交通大學 能源與動力工程學院 熱流科學與工程教育部重點實驗室，西安 710049

復合材料及多孔材料在高溫條件下的傳熱特性預測在實際工程中有重要應用價值。例如，航空航天熱防護結構的精細化設計，需要準確預測纖維增強復合材料、氣凝膠復合隔熱材料等熱防護材料的高溫傳熱特性[1-3]。新型高溫太陽能吸熱器的性能預測，需要掌握陶瓷泡沫多孔材料在高熱流密度條件下的導熱、對流和輻射耦合傳熱規律[4]。

國家數值風洞工程項目針對中國CFD軟件系統自主化的瓶頸，在轉捩與湍流模型、多相多介質計算模型、多物理場耦合計算模型、高階精度計算格式等多個方面獲得了學術進展[5]。針對復合熱防護材料高溫傳熱特性的高效預測模型，作為數值風洞中壁面熱效應預測的子模型之一，也是國家數值風洞工程的重要組成部分。

復合材料高溫傳熱特性的研究存在兩點困難。首先，由于材料使用于高溫和高熱流密度的條件，熱輻射成為材料的主要傳熱途徑之一。導熱-輻射耦合的傳熱過程涉及到導熱方程和輻射傳輸方程，其分析相對于純導熱過程更為復雜。其次，上述復合材料雖然在細觀尺度具有周期性，但其單個周期內材料的幾何結構復雜，組分物性變化劇烈。變化的物性會造成材料局部溫度的非均勻分布，從而影響材料實際性能。例如，在復合材料的氧化和燒蝕過程中，材料局部的非均勻溫度和物性會影響局部燒蝕速率，產生非均勻的燒蝕形貌[6]。在高溫熱透波材料中，溫度的非均勻可能引起局部材料介電常數的變化，從而影響透波性能。因此，建立能夠反映材料細觀特征的多尺度傳熱模型，對指導復合材料的實際應用具有重要意義。

部分復合材料傳熱特性的研究集中于材料等效導熱系數的獲取。通過重構能夠反映材料細觀特征的表征單元，采用理論分析和數值模擬等方法，建立材料等效導熱系數與材料組分體積分數、幾何特征和物性等參數之間的關系[2-3]。此類模型通常針對特定的材料結構，難以具有通用性，且無法反映表征單元內傳熱過程與宏觀傳熱過程之間的關聯。

基于多尺度漸進展開的均勻化方法，通過將溫度場分解為宏觀平均溫度場與細觀的溫度漲落，利用分析的方法建立細觀與宏觀尺度傳熱過程的關聯，并獲得等效傳熱參數[7-8]。該方法具有嚴格的數學基礎，且適用于不同的材料結構，已被廣泛應用于周期結構復合材料的傳熱和傳質問題的研究[9]。

目前傳熱過程的均勻化模型研究主要集中于導熱問題，部分研究考慮了導熱-輻射耦合傳熱問題。例如，Liu和Zhang[10]采用均勻化方法預測多孔材料導熱過程的等效導熱系數，并通過單胞內的輻射計算，補充等效輻射傳熱系數。Allaire和El Ganaoui[11]針對多孔材料的導熱-輻射耦合傳熱問題開展了多尺度漸進分析研究，其控制方程依然為導熱方程，而在材料內部孔洞壁面采用熱輻射邊界條件。Yang等[12-14]針對類似問題開展了進一步的研究工作，改進了輻射邊界條件中人為假設的小尺度參數，并將原模型拓展到二階、非穩態及三空間尺度的模擬。Haymes和Gal[15]利用此類模型研究了多孔墻體材料的傳熱問題。該模型被進一步推廣到包含導熱、對流和輻射的更復雜的傳熱過程[16-17]。在這些模型中，熱輻射只作為模型內部邊界條件作用于傳熱過程，并未考慮孔洞內部的傳熱及輻射傳輸過程。另一方面，Cao等[18-19]研究了帶有正比于溫度四次方的輻射體積熱源項的導熱方程，建立了多尺度漸進分析模型，應用于導熱-輻射耦合傳熱問題，此類模型中同樣未考慮輻射傳輸過程。

綜上可見，現有研究中對于綜合考慮導熱方程和輻射傳輸方程的多尺度模型仍鮮有報道。因此，本文將針對導熱-輻射耦合傳熱問題開展多尺度漸進分析研究，基于耦合的導熱方程和輻射傳輸方程，建立多尺度傳熱特性計算模型，進一步拓展傳熱問題多尺度方法的理論模型和應用范圍，為復合材料高溫傳熱特性的高效準確預測提供模型方法。

1 導熱-輻射耦合傳熱模型

導熱-輻射耦合傳熱模型的控制方程包括導熱方程和輻射傳輸方程，其形式為[20]

(1)

(2)

式中：Λ為能夠反映材料各向異性傳熱的導熱系數矩陣；T(x)為坐標x處的溫度；qr為輻射熱流密度；I(x,Ω)為坐標x處沿方向Ω的輻射強度；n為沿方向Ω的單位向量；β、α和σ分別為介質的衰減、吸收和散射系數，并有關系β=α+σ；σB為Stefan-Boltzmann常數；Φ(Ω,Ω′)為介質散射相函數，表示由方向Ω′散射到方向Ω上的輻射增強部分。式(1)為包含輻射源項的導熱方程，式(2) 為輻射傳輸方程。熱流密度散度與輻射強度之間的關系為

(3)

對于模型的邊界條件，本文假設導熱方程為溫度邊界條件，而輻射傳輸方程為發射率為1的黑體邊界，即在邊界處有：

(4)

對于具有周期性復雜結構的多孔材料或復合材料，由于細觀尺度組分物性的變化，可假設與組分物性相關的Λ、β、α、σ和Ω等參數是y=x/ε的函數，并具有周期性，其中ε是與材料表征單元尺度相對應的一個小周期參數。

2 導熱-輻射耦合傳熱模型的均勻化方法

針對導熱方程和輻射傳輸方程，本節將采用多尺度均勻化方法建立其宏觀尺度等效模型。通過引入坐標y=x/ε，可將溫度場和輻射強度場漸近展開為ε的多項式形式：

T(x)=T0(x,y)+εT1(x,y)+ε2T2(x,y)+…

(5)

I(x,Ω)=I0(x,y,Ω)+εI1(x,y,Ω)+…

(6)

相應的梯度算子可以表示為

(7)

將式(5)～式(7)代入式(1)～式(3)，并整理ε各階次數所對應的項，可得到ε各階的方程。

首先，在ε-2階只有對應的導熱方程：

(8)

從式(8)可得T0僅為宏觀坐標x的函數：

T0=T0(x)

(9)

其次，在ε-1階存在相應的導熱和輻射方程：

(10)

(11)

根據式(11)可知I0也僅為x和Ω的函數：

I0=I0(x,Ω)

(12)

利用式(9)，可將式(10)簡化為

(13)

根據式(13)的形式，可假設T1為

(14)

式中：N(y)為表征單元內定義的周期向量函數。將式(14)代入式(13)，即可得到N的控制方程：

(15)

最后，ε0階所對應的導熱和輻射方程為

(16)

(17)

將式(16)和式(17)在表征單元內積分，并利用T1、T2和I1的周期性，最終得到宏觀尺度的均勻化導熱-輻射耦合傳熱方程：

(18)

(19)

(20)

式中：Y為表征單元；|Y|為表征單元的體積。等效吸收、衰減系數及散射系數和相函數的乘積均為其在表征單元內的體積平均值，即

(21)

(22)

(23)

3 導熱-輻射耦合多尺度模型計算方法

3.1 多尺度模型計算流程

第2節所建立的導熱-輻射耦合多尺度模型的計算流程如圖1所示，其具體描述如下。

圖1 多尺度模型計算流程Fig.1 Calculation procedure of multiscale model

1) 針對周期性復合材料，建立宏觀尺度和細觀表征單元尺度模型并劃分網格。

2) 在細觀表征單元尺度，假設函數N(y)在表征單元內的平均值為零并具有周期性，基于式(15) 求解N(y)的各分量。

3) 利用所求得的N(y)，通過式(20)計算宏觀等效導熱系數，并利用式(21)～式(23)計算等效輻射衰減和吸收系數及散射系數與相函數之積。

4) 求解宏觀尺度導熱-輻射耦合傳熱方程式(18) 和式(19)，獲得宏觀尺度T0及I0。

5) 最終，一階精度的溫度場可以通過以下關系式計算得到：

(24)

該多尺度模型通過表征單元尺度的計算獲得宏觀等效導熱系數以及輻射傳輸方程中的等效吸收、散射和衰減系數，從而計算宏觀平均的溫度場和輻射強度場。同時，表征單元尺度的計算結果能夠為平均溫度場提供細觀尺度的修正量T1，從而進一步提升模型的計算精度。在第4節中將通過具體的案例展示多尺度模型對計算精度的提升效果。

在第2節的模型推導中引入了小周期參數ε來區分宏觀和細觀尺度。需要強調的是，當宏觀和細觀網格確定后，ε的取值并不影響計算結果。例如，若取一新的ε1，其值比ε縮小C倍，即ε1=ε/C。則有y1=Cy，其中y為參數ε所對應的細觀尺度坐標，而y1為參數ε1所對應的細觀尺度坐標。從表征單元內周期函數N的計算公式式(15) 可以看到，相應條件下計算出的N1=CN，從而推導出ε1N1=εN。因此，根據式(24)所重構出的溫度場不受到ε的具體取值的影響。實際計算中，可將ε取為細觀網格尺寸和宏觀網格尺寸的比值。

3.2 多尺度模型數值方法

從2.1節所示多尺度模型計算流程可知，模型需要求解式(15)和式(18)所對應的帶有源項的擴散方程，以及式(19)所對應的輻射傳輸方程。本文將采用有限容積方法(FVM)與格子Boltzmann方法(LBM)相結合的數值方法，利用FVM求解導熱方程，并利用LBM求解輻射傳輸方程[21]。

本文采用正方形結構化網格開展計算。在FVM方法中，通過將控制方程在控制體內積分，并將界面處通量表示為相鄰控制體節點的差分形式，獲得離散代數方程，通過三對角矩陣算法等迭代求解。

對輻射傳輸方程式(19)，考慮其非穩態形式：

(25)

式中：c為輻射傳播速度。假設輻射強度被離散為對應于離散速度ei的輻射強度Ii，則可將方程沿輻射強度傳播方向在δt時間內積分，得到[21-22]

(26)

其中：ci=|ei|為離散速度ei的大小；Φi,l為離散形式的散射相函數；wl為對應的權系數。式(26)與LBM的演化方程具有一致性，其中離散輻射強度Ii類比于LBM中的速度分布函數，而ei類比于LBM中的離散速度[23]。取eiδt與劃分的網格相一致，則離散輻射強度Ii將在δt時間內從當前網格節點遷移到相鄰的網格節點。因此，可通過LBM方法的碰撞和遷移步驟對式(26)進行迭代求解。

離散速度ei的選取影響著輻射方程求解的準確性。研究表明，采用30離散速度方向的模型模擬結果與離散坐標方法的模擬結果有較好的一致性[20]。其中，e1～6的方向為(1, 0, 0)及其反方向和置換，e7～30的方向為(2, 1, 1)及其反方向和置換，其在二維平面上的速度投影如圖2所示。由于LBM需要保證離散方向上的輻射強度在1個時間步長內遷移到相鄰網格節點，因此不同ei的速度大小會有不同，即輻射強度在不同方向上的傳播速度有差異。由于輻射傳播速度為光速，且本文主要考慮穩態問題，因此該速度差異對計算結果沒有影響[24]。

圖2 二維平面離散速度投影示意圖Fig.2 Sketch of projections of discrete velocities on 2D plane

最后，需要針對以上30離散速度方向確定其相應的權系數wi。基于離散速度的對稱性，采用以下關系式確定wi的數值：

(27)

式中：μi為ei沿x軸(坐標系的3個坐標軸分別為x、y、z)的方向余弦。求解得到w1～6=0.140 068 0π，w7～30=0.131 649 7π。

4 模型應用及驗證

4.1 問題描述

為驗證本文所建立的導熱-輻射耦合多尺度模型的有效性，本節將利用該模型模擬如圖3所示的周期性復合結構內的導熱-輻射耦合傳熱問題，其計算區域由基體材料及填充其間的10×10陣列組成。本文主要考慮2種結構，一種是由SiO2氣凝膠基體材料及填充其間的顆粒狀TiO2遮光劑組成，其計算區域大小為0.8 mm×0.8 mm，TiO2顆粒直徑為20 μm，因此其表征單元為80 μm×80 μm的包含單個TiO2顆粒的正方形結構，TiO2顆粒的體積分數為4.9%。另一種為8 mm×8 mm的SiC陶瓷泡沫材料，其基體為SiC陶瓷骨架，內部為邊長0.7 mm 的正方形空氣孔隙，其表征單元為0.8 mm×0.8 mm的包含單個孔隙的正方形結構，材料孔隙率為77%。假設SiC、空氣、SiO2氣凝膠和TiO2顆粒的導熱系數、輻射吸收和散射系數為700 K條件下的常物性，其具體數值如表1所示。其中假設SiC、空氣和TiO2均無散射效應，同時忽略空氣對輻射的吸收效應。另外，假設各組分均為各向同性介質，即令散射相函數Ф≡1。計算中不考慮不同材料界面處的折射和散射效應，而在不同材料所占據的計算單元內，直接采用相應的物性參數。對導熱問題，其界面處的等效導熱系數為相鄰網格導熱系數的調和平均值。計算區域的上下邊界為周期邊界條件，左右邊界為定溫邊界條件。對SiO2氣凝膠復合材料，左右邊界分別為800 K和600 K，對SiC陶瓷泡沫，左右邊界分別為1 100 K 和300 K。

表1 SiC、空氣、SiO2氣凝膠及TiO2的物性參數Table 1 Physical properties of SiC, air, SiO2 aerogel and TiO2

圖3 計算區域示意圖Fig.3 Sketch of computational domain

模擬中，首先將整體計算區域劃分為800×800的計算網格，基于不同區域的材料物性，直接求解導熱-輻射耦合傳熱方程式(1)和式(2)，作為數值模擬的參照結果。在多尺度模擬中，將整體計算區域劃分為100×100的計算網格，采用平均物性參數計算T0和I0，將表征單元劃分為80×80的網格開展表征單元問題的求解。計算結果收斂的判據為2次相鄰迭代中的溫度場及輻射強度場的相對誤差之和小于10-8。

4.2 計算結果對比

首先針對SiO2氣凝膠復合材料開展計算。其表征單元內N(y)的計算結果如圖4所示。由于材料表征單元結構具有各向同性特征，其Nx和Ny分量具有對稱性。基于N(y)和式(20)，可計算出材料的宏觀等效導熱系數為0.029 1 W/(m·K)。同時，式(21)～式(23)計算得出的材料等效吸收及散射系數分別為3.71 cm-1及0.236 cm-1。

圖4 SiO2氣凝膠復合材料的N(y)計算結果Fig.4 Calculation results of N(y) for SiO2 aerogel composites

基于以上等效參數，利用宏觀平均導熱-輻射耦合傳熱方程，計算得出材料宏觀等效溫度場，并結合式(24)重構材料多尺度溫度場。多尺度計算得到的溫度場與完整計算得到的溫度場如圖5所示。

圖5 SiO2氣凝膠復合材料的溫度場計算結果Fig.5 Calculation results of temperature fields for SiO2 aerogel composites

同理，采用多尺度方法計算SiC陶瓷泡沫材料的導熱-輻射耦合傳熱問題。其表征單元內N(y)的計算結果如圖6所示，對應的等效導熱系數為14.24 W/(m·K)，等效輻射吸收系數為84.0 cm-1。同時，采用多尺度模型求解得到的溫度場與完整求解的溫度場對比結果如圖7所示。

圖6 SiC陶瓷泡沫的N(y)計算結果Fig.6 Calculation results of N(y) for SiC ceramic foam

從圖5和圖7可以看到，基于本文所建立的多尺度方法計算得出的溫度場與完整求解得出的溫度場具有較好的一致性。在圖8中進一步給出了沿y=L/2的溫度分布情況，其中L為計算區域邊長。可見多尺度模擬結果與完整計算結果符合良好，且能夠體現出溫度場在不同介質中的非均勻變化。

圖7 SiC陶瓷泡沫的溫度場計算結果Fig.7 Calculation results of temperature fields for SiC ceramic foam

圖8 y=L/2處溫度場計算結果Fig.8 Calculation results of temperature fields on y=L/2

令完整求解出的溫度場為T(i,j)，多尺度方法計算出的溫度場為Tm(i,j)，定義多尺度方法計算的溫度場誤差為

(28)

2個問題的多尺度模型計算結果誤差如表2所示，其中列舉了宏觀平均溫度場T0的誤差及由式(24)重構的多尺度溫度場的誤差。可以看到采用多尺度模型能夠獲得具有較高精度的溫度場，且通過引入高階溫度場T1，能夠有效提高模型精度。以SiO2氣凝膠復合材料的多尺度溫度場計算為例，在圖9中展現了多尺度溫度場與完整計算結果相比的絕對誤差和相對誤差分布。可以看到，溫度場局部的絕對誤差小于1 K，相對誤差小于0.13%。在計算區域的中部誤差相對較大，主要原因是等效物性的偏差導致的平均溫度分布差異。同時，在添加物位置也存在一定的誤差波動。

圖9 SiO2氣凝膠復合材料多尺度溫度場的誤差分布Fig.9 Error distribution for multiscale temperature field of SiO2 aerogel composites

表2 多尺度模型計算結果相對誤差

在表3中對比了多尺度模擬與完整模擬的計算時間。可以看到，多尺度模擬的計算時間減少了近4個數量級，能夠有效提升復合結構傳熱特性的預測效率。本文算例中假設物性參數為常數，因此表征單元內的模擬只需要進行一次，獲得等效物性參數。針對實際中隨溫度變化的物性參數，可在一定溫度范圍內分別計算表征單元模型，獲得宏觀等效參數隨溫度的變化規律，用于宏觀模擬。由于表征單元模擬的計算時間較小，多尺度模型依然能夠有效提高數值模型計算效率。

表3 多尺度模型計算時間對比

4.3 三維平紋編織復合材料傳熱模擬

作為模型在更復雜結構中的應用舉例，采用所建立的多尺度模型模擬了三維平紋編織復合材料內的傳熱過程。整體計算區域和表征單元結構如圖10所示，利用多尺度模型計算得到的溫度分布如圖11所示，其溫度場計算結果能夠反映細觀結構內的波動現象，因此所建立的多尺度模型能夠被應用于更為復雜的復合材料結構內部傳熱過程的預測。

圖10 平紋編織復合材料示意圖Fig.10 Sketch of plane woven composite

圖11 平紋編織復合材料溫度場Fig.11 Temperature field of plane woven composite

5 結 論

本文針對復合材料高溫條件下的導熱-輻射耦合傳熱問題，基于導熱方程及輻射傳輸方程，利用多尺度均勻化方法，建立了導熱-輻射耦合傳熱問題的多尺度計算模型，主要研究成果為

1) 通過多尺度分析，建立了材料周期性表征單元內的計算模型，并給出了宏觀等效導熱系數與表征單元計算模型之間的關聯。同時，分析表明宏觀等效輻射吸收、衰減、散射系數和散射相函數等可通過表征單元內的體積平均獲得。

2) 利用均勻化分析方法，結合FVM與LBM數值方法，建立溫度場多尺度計算模型。采用SiO2氣凝膠復合材料和SiC陶瓷泡沫等材料的導熱-輻射耦合傳熱過程作為算例，驗證了多尺度模型的有效性。在本文所采用的數值方法、收斂判據及常物性參數條件下，多尺度模型相對于材料整體建模，能夠顯著減少計算時間。

同時，本文所建立的多尺度模型還存在著以下局限：

1) 該模型推導的前提是復合材料具有周期性結構，實際中復合材料的微觀結構存在著一定隨機性，因此需要選取在平均意義上具有代表性的表征單元結構。

2) 模型中宏觀尺度溫度場T0滿足邊界條件，但修正項T1會使得邊界上的溫度偏離設定的溫度邊界條件，在未來研究中需要進一步考慮邊界層的修正。

3) 當添加物的特征尺寸與輻射的波長相當時，需要采用更為精細的輻射模型。

未來可通過引入溫度場的二階修正和輻射強度場的一階修正，進一步提升模型的預測精度。