导热-辐射耦合传热的多尺度分析和数值模型

2021-10-20童自翔李明佳李冬

航空学报 2021年9期

童自翔，李明佳，李冬

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所示，其温度场计算结果能够反映细观结构内的波动现象，因此所建立的多尺度模型能够被应用于更为复杂的复合材料结构内部传热过程的预测。