Fast calculation of radiative heat transfer coefficient between diffuse and non-diffuse surfaces

辐射换热是热传递的三种基本方式之一，对计算空间目标外表面温度分布和红外辐射特性有着重要影响。辐射换热主要与目标表面温度分布、光谱发射/反射特性，以及几何外形有关。由于辐射换热计算量庞大，如何提高辐射换热的计算效率就显得尤为重要^[1-4]。

辐射传递系数是从一个表面发出的能量中，最终被另一个表面吸收的能量所占的份额^[5]。已知表面间的辐射传递系数就可以快速计算出物体间的辐射换热量。黑体和灰体是研究辐射换热的理想化物体^[6]。对于黑体表面，辐射传递系数即辐射角系数^[7]（角系数是从一个表面发出的能量中，直接到达另一个表面的能量所占的份额）。

对于灰体表面，根据光线反射形式的不同，可分为漫反射、镜面反射、以及不同类型组合而成的复杂反射。对于漫反射灰体，已知其辐射角系数就可以快速计算出其辐射传递系数，张伟清通过构建漫反射表面的辐射传递系数与角系数间之间的关系实现立方体内表面辐射传递系数的快速求解^[8]。由于自然界中的大部分物体，近似意义上都可以等效为漫反射灰体，因此通常采用计算其表面间角系数来求解其辐射换热。角系数的计算通常可以采用蒙特卡洛法、单位球法、半立方体法等方法。例如Fei He利用蒙特卡洛法计算冷扎带与W形辐射管之间的辐射角系数^[9]；Arvind 利用单位球法计算出任意定向多边形的辐射角系数的解析表达式^[10]；Li Hailin 采用直接面积积分法和半立方法分别计算冷坩埚悬浮炉面元间的热辐射^[11]。

对于系统中同时存在漫反射和非漫反射表面组成的混合表面模型，无法通过迭代角系数来计算辐射传递系数，而几乎只能使用蒙特卡洛法 (Monte-Carlo method, MCM) 计算辐射传递系数^[12]。蒙特卡洛法是辐射传递系数计算的“金标准”，该方法需要目标表面生成大量随机光线并进行多次射线追踪，计算结果准确性高，时间消耗严重。为了提高蒙特卡洛射线追踪的计算效率，许多基于GPU的技术被开发出来，如八叉树 (Octree)、KD树^[13-14]等。然而，对于大量面元，计算量仍然相当大，计算效率仍需进一步提高。文中针对蒙特卡洛方法计算辐射传递系数耗时过长的问题，提出了一种计算漫反射与非漫反射表面间辐射传递系数的快速方法。

蒙特卡洛法作为一种随机模拟方法，本质在于通过大量重复实验，利用统计学获得计算结果。如图1所示，蒙特卡洛法用于计算面元间的辐射传递系数可描述为如下步骤：

Figure 1.  Monte-Carlo method for calculating radiative heat transfer coefficient

1）对于任意发射面元i，在面元上随机位置、以随机方向生成多根射线，其中每根射线携带相同的初始能量；

2）进行射线追踪。根据发射能量光束的角度和位置计算相交面元和其交点，若射线与非漫反射表面相交（如镜反射面元j、l），根据镜面面元吸收率计算射线被吸收能量，并根据镜反射定律在该点求解射线反射方向；若与漫反射表面相交（如面元p、k），根据表面吸收率计算面元吸收的射线能量，并在交点位置随机生成反射方向；如果射线与除自身之外的任意面元都不相交，则光线射空，该射线携带的能量被舍弃，不计入任何面元吸收的能量之内。

3）重复过程2），直到射线能量低于追踪阈值或射空，结束射线追踪过程；

4）依次完成面元i上所有射线的追踪，统计此过程中每个面元吸收的能量占面元i发射光束总能量的比例，即为面元i对其他面元的辐射传递系数；

5）依次完成所有面元的辐射传递系数计算。

由以上原理可知，蒙特卡洛法计算面元间辐射传递系数时，需要对大量射线进行多次射线追踪，如果系统面元反射率较高，计算时间与反射率呈指数倍数增长，计算效率较低。

需要指出的是，文中，漫反射表面和镜反射表面指的是反射方向的不同，在光谱上仍被视为灰体表面，即表面光谱发射率ε_λ在全谱段内为常数。

根据漫反射表面随机反射的原理，提出一种入射光束在漫反射表面反射能量的表示方法。如图2所示，蒙特卡洛方法中，当入射光线以能量E_in与漫反射面元相交时，面元以吸收率ε吸收能量，反射能量E_ref 以随机方向进行反射，然后继续进行下一轮的射线追踪，直到光束能量低于能量阈值，追踪结束。

Figure 2.  (a) The representation of reflected light by diffuse surface in Monte-Carlo method; (b) Fast method

对于改进的快速方法而言，由于漫反射表面对入射能量的反射在半空间满足朗伯余弦定律，因此可以假设反射光束不是1根随机光线，而是无数根向半空间随机反射的子光线构成的光线集，光线集的总能量为E_ref，分布服从漫反射分布。若该漫反射表面的辐射传递系数已知，则该面元的反射能量能够被各个面元依辐射传递系数被吸收，射线追踪结束，可以显著减少射线追踪次数，实现快速计算。

改进后的辐射传递系数的计算方法可描述如下：

1) 假设系统中有N₁个漫反射面元和N₂个非漫反射面元，先利用蒙特卡洛法计算出全部N₁个漫反射面元对其他面元的辐射传递系数；

2) 计算非漫反射面元的辐射传递系数步骤可分为：

① 面元光线生成方式同蒙特卡洛法步骤1)；

② 进行射线追踪和能量计算：如图3所示，若光束在进行射线追踪过程中与非漫反射面元相交（如镜反射面元j、l），则按照表面反射特性进行能量吸收和光线反射计算，此时方式与蒙特卡洛一致；若射线与漫反射面元相交（如面元p、k），则用该漫反射面元对其他面元的辐射传递系数与反射能量相乘，计算其他面元吸收的反射能量，结束射线追踪过程，避免后续的射线追踪，提高计算效率；若射线射空，处理方式与蒙特卡洛法一致。

Figure 3.  Fast calculation method for calculating the radiative heat transfer coefficient of non-diffuse surface element i

③同蒙特卡洛法步骤4）；

④同蒙特卡洛法步骤5）。

建立6面元立方体空腔模型，计算立方体内壁面元间的辐射传递系数，对比辐射传递系数计算精度和计算时间，说明快速方法相对于蒙特卡洛方法的优势。为简化计算并不失一般性，假设模型中只存在漫反射和镜反射两种面元类型。

如图4(a)所示建立参考坐标系XYZ，立方体模型及面元编号设置如下，1号面元AA'B'B，2号面元BB'C'C，3号面元CC'D'D，4号面元DD'A'A，5号面元为顶面A'B'C'D'，6号面元为底面ABCD。设定面元吸收率为ε，反射率ρ=1−ε。

Figure 4.  (a) Illustration of 3D geometry model and its surface ID; (b) Illustration of random positions and direction vectors for emitted rays

假设面元向半空间发射的光线总数为m，则每束光线携带的初始能量为：

式中：T为面元的瞬态温度；ρ为面元反射率；σ为Steven-Boltzmann常数。

如图4(b)所示，射线在三角形内的随机发射位置P可以通过公式 (2) 确定^[15]：

式中：P_a、P_b和P_c为面元内三个顶点的坐标向量，R₁和R₂是[0,1]内均匀分布集的随机数。四边形由两个三角形组成，并用相同的方式产生射线的随机发射位置P。

射线的发射方向F(θ, φ) 为弧度表示的指向内壁空间的天顶角和方位角^[12]，满足公式 (3)：

式中：R_θ、R_φ为[0,1]内均匀分布的随机数。

快速方法计算流程如图5所示。在射线追踪过程中，单束光线携带的初始能量E是依面元吸收率ε吸收而逐渐衰减的。当光线进行第一次追踪时，相交面元吸收能量为εE，光线反射能量为ρE。当光线继续进行第二次追踪时，相交面元吸收能量为ρεE，反射光束能量为ρ²E…第k次追踪时，相交面元吸收能量为ρ^k-1εE，光线反射能量为ρ^kE。给定能量阈值η，若反射能量低于阈值，则视为光线能量被完全吸收，射线追踪结束^[16]。正六面体内腔模型无射线射空情况。

对于镜反射面元的随机发射光线，若在追踪中与镜反射面元相交，则按反射定律计算反射；若光线经过k次镜面反射后到达漫反射面元i，则该漫反射面元吸收能量为ρ^k-1εE，反射能量ρ^kE依该面元的辐射传递系数被各个面元吸收。如漫反射面元i对面元j的辐射传递系数为D_ij，则此反射射线中被面元j吸收的能量为ρ^kED_ij，射线追踪结束。若射线在此过程中仅与镜反射面元相交，当光束能量低于阈值η时，也结束追踪过程。

依次完成m根射线的追踪和能量计算过程，统计此过程中各个面元的吸收能量，如面元i对面元j的辐射传递系数可表示为：

式中：D_ij为面元i对面元j的辐射传递系数；E_j为面元j吸收的能量；m为面元i的发射光线数；E为单根光线的初始能量。

本模型中，立方体内壁面元间的辐射传递系数是一个 6×6 的矩阵。

Figure 5.  Fast calculation of radiative heat transfer coefficient for model with diffuse and specular surfaces

假设模型4(a) 中1号(AA'B'B)、2号面元(BB'C'C)为漫反射表面，其余四个面元为镜反射表面，每个面元随机发射光线数目m为10万，面元反射率ρ为0.4，射线追踪的截止能量阈值η为1‰。

表1所示为蒙特卡洛方法计算面元间的辐射传递系数的计算结果，表格中每一行为该面元对其他面元（包括其自身）的辐射传递系数；列“Integrality”意为完整性，即该面元对其他面元（含自身）的辐射传递系数的总和。表2为快速方法的计算结果，其中，漫反射面元1、2号的辐射传递系数采用蒙特卡洛方法计算，其余面元调用了漫反射面元的计算结果；列“MAE(Mean Absolute Error)”是表2结果相较于表1结果的平均绝对误差，计算公式如下：

式中：y_i、x_i分别表示第i个面元的蒙特卡洛计算结果和快速计算结果。

由于模型的辐射传递系数难以获得解析解，因此使用蒙特卡洛的计算结果作为参考值，来衡量快速方法的计算精度。

对比表1和表2，由于1、2号面元为漫反射面元，均采用蒙特卡洛方法计算，计算结果没有很大区别，MAE为1‰左右。表2中3~6号面元的辐射传递系数计算结果与表1几乎一致，误差主要集中在千分位上，MAE在该例中为2‰左右，误差较大的项主要集中在镜反射面元对自身的辐射传递系数上。其原因主要如下：一方面，快速方法在计算过程中需要调用漫反射面元的辐射传递系数以加速计算。由于蒙特卡洛方法计算得到的辐射传递系数不是“真值”，存在一定的误差，在调用过程中存在着误差累积，使快速方法的误差略大于蒙特卡洛方法；另一方面，由于面元对自身的辐射传递系数较小，使得面元辐射出去的光线中，“击中自身”的数量偏少，在光线总数量不充分大的情况下，容易产生较大的偏差。

对能量完整性进行分析，可知快速方法的镜反射面元的完整性更高。蒙特卡洛射线追踪中，能量阈值取值不等于0，射线能量低于阈值后即停止追踪，因此完整性只能达到0.99934。而使用快速方法计算镜反射面元的辐射传递系数时，若射线击中漫反射面元，则调用漫反射面元的辐射传递系数，剩余能量被完全定义为该面元上的漫发射光束集，因此对剩余能量的处理更加完全，能量完整性更高，达到0.99985，提高了能量的利用率。

图6(a)所示为面元的反射率ρ=0.4（发射率ε=0.6）时，两种方法获取辐射传递系数的计算耗时。保持其他参数不变，将面元的反射率提高至ρ=0.8（ε=0.2），两种方法的计算时间对比如图6(b)所示。

Figure 6.  (a) Comparison of computational time between Monte-Carlo method and Fast method when the inner surface reflectivity of the cube is ρ=0.4; (b) ρ=0.8

对比两图可知，当面元反射率ρ=0.8时，由于反射率的提高，射线的平均反射次数显著提高，蒙特卡洛方法的计算时间也随之增加。但对于快速方法，当计算镜反射面元的随机发射光线时，只要射线击中漫反射面元，便可以结束追踪，计算耗时受反射率影响不显著。在图6(b)示例中，镜反射面元的辐射传递系数计算耗时不足蒙特卡洛方法的1/10。面元反射率越高，该方法的计算优势越明显。

图7所示为一个相对复杂的单端开口L型腔模型。模型共13个面元，面元编号和法矢量方向在图中已标出。其中，第1、6、7号面元（已置灰）为漫反射表面，其余面元为镜反射表面，Y轴正方向Open Cavity箭头所指位置为腔体开口位置。

Figure 7.  L-shape unenclosed cavity model

内表面面元间的辐射传递系数仿真参数与前述模型相同，面元初始发射光线数量m为10万，射线能量阈值η为1‰。

表3所示为面元反射率ρ=0.4、0.8时，使用快速方法计算得到的目标各面元辐射传递系数计算误差（与蒙特卡洛方法相比）和能量吸收份额，其中第1、6、7号面元为漫反射面元。

由表3可见，同一面元反射率ρ下，由于漫反射面元1、6、7号均采用蒙特卡洛方法，其MAE数值较小；镜反射面元由于调用了漫反射面元的计算结果，存在误差累积且面元辐射传递系数之和Integrality相对蒙特卡洛法较高，MAE数值相对较大，但这些误差都在千分位上，能够保证计算精度。

此外，还可以看出，由于腔体不封闭的原因，导致各个面元的出射能量经过多次反射后有部分能量逸出，致使出射能量被吸收的比例小于1；尤其是面元8正对腔体开口，逸出能量份额最大。此外，还可以看出，反射率ρ越大，被吸收的能量之和越小。可以想象，假如某条光线经过k次反射逸出腔体射空，面元的反射率越高，逸出光线携带的能量ρ^kE就越大，其被面元吸收的能量就越少。

图8所示为反射率ρ=0.4和0.8时，两种方法的计算时间对比。可以看到，漫反射面元1、6、7号两种方法的计算时间基本一致，各个镜反射面元的计算时间都得到了不同程度的降低，具体情况与面元全部射线的平均反射次数有关。此外，随着面元反射率ρ的增高，镜反射面元的辐射传递系数计算效率也得到进一步提升。

Figure 8.  (a) Comparison of computational time between Monte-Carlo method and Fast method when the inner surface reflectivity of the L-shape cavity is ρ=0.4; (b) ρ=0.8

下面从理论角度对快速方法的计算效率进行分析和说明。为简化分析过程，此处假定目标为封闭模型，并且光线击中各个面元的概率相同。

1) 射线平均追踪次数

对于蒙特卡洛法而言，单根射线的追踪次数n是根据面元反射率ρ 得出的，之间的关系为：

式中：$\left\lceil\;\;\right\rceil $表示向上取整；η为能量阈值，上述实验中取值为1‰。

图9中红色曲线为蒙特卡洛方法的射线平均追踪次数随反射率ρ的变化关系。从右侧红色坐标轴可见，当ρ=0.9时，反射次数可达66次之多。

Figure 9.  Illustration of the average number of ray-tracing times

对于快速方法，其非漫反射表面的发射光线的追踪次数和系统中漫反射面元的占比有关：光线击中漫反射面元，追踪进程就结束。假如射线追踪次数为k，单根射线的追踪次数的概率p_k可表示为：

式中：N、N₁、N₂分别为系统面元总数、漫反射面元数量、非漫反射面元数量。

因此，非漫反射面元发射光束的平均追踪次数$ \bar{n} $可表示为：

当面元反射率较大时，非漫反射面元平均追踪次数近似表示为：

以正方体模型为例 (N=6)，射线追踪次数随面元反射率ρ和漫反射面元数量N₁的变化关系如图9中黑色曲线所示。

可知，对于非漫反射面元，射线追踪次数随面元反射率ρ缓慢增加，最大反射次数如公式 (9) 所述。而漫反射面元的射线追踪次数随面元反射率ρ急剧升高，为快速方法提升计算效率创造了条件。

表4所示为采用两种方法计算得到的非漫反射面元的射线平均追踪次数比值（快速方法/蒙特卡洛方法）随面元反射率ρ、漫反射面元占总面元的比值C_d的变化关系矩阵。如果单次射线追踪的时间为定值，则追踪次数之比即为镜反射面元的计算耗时之比。

从表4可知，对于非漫反射表面，快速方法相对于蒙特卡洛方法，其计算优势随着面元反射率ρ的增加而提高，并随着系统中漫反射面元占比的增加而提升。

2) 辐射传递系数计算总时间

尽管快速方法能够显著提高非漫反射面元的射线追踪效率，但是由于漫反射面元的辐射传递系数计算时间与蒙特卡洛方法持平，系统整体计算耗时提升受漫反射面元占比影响显著。假如单次射线追踪时间为t，则快速方法和蒙特卡洛计算总时间比τ如公式 (10) 所示。

式中：C_d为漫反射面元占模型总面元的比例因子；C_n为非漫反射面元使用两种方法的平均追踪次数之比（快速方法/蒙特卡洛方法）。

表5所示为快速方法和蒙特卡洛方法计算全部面元的辐射传递系数的耗时比值随面元反射率ρ、漫反射面元占比C_d的变化关系矩阵。从表中可以看出，随着反射率ρ的增加，快速方法的计算效率在也在提高，但计算时间比值不会小于系统中漫反射面元的占比。例如，第3列中，计算时间比值最小为0.513，趋近但大于漫反射面元的占比3/6。这是因为无论快速方法中非漫反射面元的计算时间如何降低，计算总耗时中，漫反射面元的那部分计算时间与蒙特卡洛方法相比是持平的。

同时可以看出，随着漫反射面元占比的增加，虽然从上节中分析得到镜反射面元的计算效率得到提升，但由于加入了漫反射面元的计算时间，使得整体的计算时间呈上升趋势，仅在反射率较小时（例如ρ=0.1~0.5时）有所波动。

为验证理论分析的正确性，以6面元立方体模型为例，当面元反射率ρ=0.8 时，通过改变漫反射面元的数量，分别使用两种方法计算内壁面元间的辐射传递系数，获取计算耗时比值，并与理论值进行对比，计算结果如表6所示。可见，上述理论分析结果是可信的。

文中针对蒙特卡洛方法计算辐射传递系数计算时间过长的问题，提出了一种计算漫反射灰体表面与非漫反射灰体表面间的辐射传递系数的快速方法，并通过案例和理论分析对计算精度和计算效率进行了验证和分析。结果表明：快速方法在保证计算精度的同时，能够有效提高非漫反射面元的辐射传递系数计算效率；面元的反射率越大、系统中漫反射面元的占比越高，非漫反射面元的计算效率优势越明显。同时，系统的辐射传递系数计算总时间随反射率的增大而降低，并受制于系统中漫反射面元的占比。该方法为估计漫反射与非漫反射表面间的辐射传递系数提供了更多的选择。