-
辐射换热是热传递的三种基本方式之一，对计算空间目标外表面温度分布和红外辐射特性有着重要影响。辐射换热主要与目标表面温度分布、光谱发射/反射特性，以及几何外形有关。由于辐射换热计算量庞大，如何提高辐射换热的计算效率就显得尤为重要[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所示，蒙特卡洛法用于计算面元间的辐射传递系数可描述为如下步骤：
图 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所示，蒙特卡洛方法中，当入射光线以能量Ein与漫反射面元相交时，面元以吸收率ε吸收能量，反射能量Eref 以随机方向进行反射，然后继续进行下一轮的射线追踪，直到光束能量低于能量阈值，追踪结束。
图 2 (a) 蒙特卡洛方法中的漫反射表面反射光线表示；(b) 快速方法中的漫反射表面反射光束表示
Figure 2. (a) The representation of reflected light by diffuse surface in Monte-Carlo method; (b) Fast method
对于改进的快速方法而言，由于漫反射表面对入射能量的反射在半空间满足朗伯余弦定律，因此可以假设反射光束不是1根随机光线，而是无数根向半空间随机反射的子光线构成的光线集，光线集的总能量为Eref，分布服从漫反射分布。若该漫反射表面的辐射传递系数已知，则该面元的反射能量能够被各个面元依辐射传递系数被吸收，射线追踪结束，可以显著减少射线追踪次数，实现快速计算。
改进后的辐射传递系数的计算方法可描述如下：
1) 假设系统中有N1个漫反射面元和N2个非漫反射面元，先利用蒙特卡洛法计算出全部N1个漫反射面元对其他面元的辐射传递系数；
2) 计算非漫反射面元的辐射传递系数步骤可分为：
① 面元光线生成方式同蒙特卡洛法步骤1)；
② 进行射线追踪和能量计算：如图3所示，若光束在进行射线追踪过程中与非漫反射面元相交（如镜反射面元j、l），则按照表面反射特性进行能量吸收和光线反射计算，此时方式与蒙特卡洛一致；若射线与漫反射面元相交（如面元p、k），则用该漫反射面元对其他面元的辐射传递系数与反射能量相乘，计算其他面元吸收的反射能量，结束射线追踪过程，避免后续的射线追踪，提高计算效率；若射线射空，处理方式与蒙特卡洛法一致。
图 3 非漫反射表面的辐射传递系数计算示意图
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−ε。
图 4 (a) 三维几何模型与面元编号；(b) 随机发射光线的位置和方向
Figure 4. (a) Illustration of 3D geometry model and its surface ID; (b) Illustration of random positions and direction vectors for emitted rays
假设面元向半空间发射的光线总数为m，则每束光线携带的初始能量为：
$$ E=\frac{(1-\rho) \sigma T^{4}}{m} $$ (1) 式中：T为面元的瞬态温度；ρ为面元反射率；σ为Steven-Boltzmann常数。
如图4(b)所示，射线在三角形内的随机发射位置P可以通过公式 (2) 确定[15]：
$$ \begin{aligned} P =\left(1-\sqrt{R_{1}}\right) P_{\mathrm{a}}+\sqrt{R_{1}}\left(1-R_{2}\right) P_{\mathrm{b}} +\left(R_{2} \sqrt{R_{1}}\right) P_{\mathrm{c}} \end{aligned} $$ (2) 式中：Pa、Pb和Pc为面元内三个顶点的坐标向量，R1和R2是[0,1]内均匀分布集的随机数。四边形由两个三角形组成，并用相同的方式产生射线的随机发射位置P。
射线的发射方向F(θ, φ) 为弧度表示的指向内壁空间的天顶角和方位角[12]，满足公式 (3)：
$$ \left\{\begin{array}{l} \theta=\arcsin \sqrt{R_{\theta}} \\ \varphi=2 \pi R_{\varphi} \end{array}\right. $$ (3) 式中：Rθ、Rφ为[0,1]内均匀分布的随机数。
-
快速方法计算流程如图5所示。在射线追踪过程中，单束光线携带的初始能量E是依面元吸收率ε吸收而逐渐衰减的。当光线进行第一次追踪时，相交面元吸收能量为εE，光线反射能量为ρE。当光线继续进行第二次追踪时，相交面元吸收能量为ρεE，反射光束能量为ρ2E…第k次追踪时，相交面元吸收能量为ρk-1εE，光线反射能量为ρkE。给定能量阈值η，若反射能量低于阈值，则视为光线能量被完全吸收，射线追踪结束[16]。正六面体内腔模型无射线射空情况。
对于镜反射面元的随机发射光线，若在追踪中与镜反射面元相交，则按反射定律计算反射；若光线经过k次镜面反射后到达漫反射面元i，则该漫反射面元吸收能量为ρk-1εE，反射能量ρkE依该面元的辐射传递系数被各个面元吸收。如漫反射面元i对面元j的辐射传递系数为Dij，则此反射射线中被面元j吸收的能量为ρkEDij，射线追踪结束。若射线在此过程中仅与镜反射面元相交，当光束能量低于阈值η时，也结束追踪过程。
-
依次完成m根射线的追踪和能量计算过程，统计此过程中各个面元的吸收能量，如面元i对面元j的辐射传递系数可表示为：
$$ D_{i j}=\frac{E_{j}}{m E} $$ (4) 式中：Dij为面元i对面元j的辐射传递系数；Ej为面元j吸收的能量；m为面元i的发射光线数；E为单根光线的初始能量。
本模型中，立方体内壁面元间的辐射传递系数是一个 6×6 的矩阵。
图 5 漫反射与镜反射表面模型辐射传递系数快速计算
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结果的平均绝对误差，计算公式如下：
表 1 立方体内表面各面元之间的辐射传递系数（蒙特卡洛方法）
Table 1. Radiative heat transfer coefficient for the inner surfaces of the cube (Monte-Carlo method)
Facet 1(AA'B'B) 2(BB'C'C) 3(CC'D'D) 4(DD'A'A) 5(A'B'C'D') 6(ABCD) Integrality 1 0.05897 0.17917 0.21741 0.18131 0.18132 0.18116 0.99934 2 0.17956 0.05924 0.18051 0.21850 0.18032 0.18122 0.99934 3 0.21776 0.18120 0.04888 0.18594 0.18307 0.18249 0.99934 4 0.18201 0.21647 0.18750 0.04936 0.18190 0.18212 0.99934 5 0.17921 0.17974 0.18220 0.18333 0.07049 0.20437 0.99934 6 0.17872 0.17926 0.18274 0.18213 0.20545 0.07105 0.99934 表 2 立方体内表面各面元之间的辐射传递系数（快速方法）
Table 2. Radiative heat transfer coefficient for the inner surfaces of the cube (Fast method)
Facet 1(AA'B'B) 2(BB'C'C) 3(CC'D'D) 4(DD'A'A) 5(A'B'C'D') 6(ABCD) Integrality MAE 1 0.05937 0.18057 0.21859 0.18222 0.17886 0.17975 0.99934 1.29227‰ 2 0.18068 0.05972 0.18135 0.21884 0.17903 0.17972 0.99934 0.93054‰ 3 0.22136 0.18010 0.04669 0.18631 0.18214 0.18324 0.99985 1.48970‰ 4 0.18218 0.22024 0.18636 0.04688 0.18098 0.18323 0.99985 1.59877‰ 5 0.17933 0.18143 0.18179 0.18329 0.06566 0.20834 0.99985 1.84531‰ 6 0.18179 0.18012 0.18301 0.18145 0.20777 0.06570 0.99985 2.09184‰ $$ {MAE}=\frac{1}{6} \sum_{i=1}^{6}\left|y_{i}-x_{i}\right| $$ (5) 式中：yi、xi分别表示第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)所示。
图 6 (a) 立方体内表面反射率ρ=0.4时蒙特卡洛方法与快速方法的计算时间对比; (b) ρ=0.8时的计算时间对比
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箭头所指位置为腔体开口位置。
图 7 L型非封闭内腔模型
Figure 7. L-shape unenclosed cavity model
内表面面元间的辐射传递系数仿真参数与前述模型相同，面元初始发射光线数量m为10万，射线能量阈值η为1‰。
表3所示为面元反射率ρ=0.4、0.8时，使用快速方法计算得到的目标各面元辐射传递系数计算误差（与蒙特卡洛方法相比）和能量吸收份额，其中第1、6、7号面元为漫反射面元。
表 3 L型腔计算误差与能量完整性（快速方法）
Table 3. MAE and integrality of radiative heat transfer coefficient for L-shape cavity model (Fast method)
Facet ρ=0.4 ρ=0.8 MAE Integrality MAE Integrality 1 1.11640‰ 0.98594 0.61132‰ 0.85933 2 1.29494‰ 0.94656 2.29617‰ 0.77797 3 0.55498‰ 0.79605 1.91049‰ 0.64609 4 1.19067‰ 0.79290 1.34883‰ 0.64261 5 6.59430‰ 0.97410 4.82850‰ 0.80405 6 0.98649‰ 0.98910 0.64184‰ 0.85692 7 0.64990‰ 0.93160 0.56992‰ 0.76133 8 1.33612‰ 0.65359 1.32864‰ 0.41103 9 5.29473‰ 0.94401 5.60488‰ 0.77703 10 1.45282‰ 0.79286 1.44651‰ 0.65081 11 1.03800‰ 0.79272 1.53688‰ 0.64366 12 2.62882‰ 0.94494 3.66704‰ 0.77457 13 5.80383‰ 0.98413 4.92681‰ 0.84846 由表3可见，同一面元反射率ρ下，由于漫反射面元1、6、7号均采用蒙特卡洛方法，其MAE数值较小；镜反射面元由于调用了漫反射面元的计算结果，存在误差累积且面元辐射传递系数之和Integrality相对蒙特卡洛法较高，MAE数值相对较大，但这些误差都在千分位上，能够保证计算精度。
此外，还可以看出，由于腔体不封闭的原因，导致各个面元的出射能量经过多次反射后有部分能量逸出，致使出射能量被吸收的比例小于1；尤其是面元8正对腔体开口，逸出能量份额最大。此外，还可以看出，反射率ρ越大，被吸收的能量之和越小。可以想象，假如某条光线经过k次反射逸出腔体射空，面元的反射率越高，逸出光线携带的能量ρkE就越大，其被面元吸收的能量就越少。
图8所示为反射率ρ=0.4和0.8时，两种方法的计算时间对比。可以看到，漫反射面元1、6、7号两种方法的计算时间基本一致，各个镜反射面元的计算时间都得到了不同程度的降低，具体情况与面元全部射线的平均反射次数有关。此外，随着面元反射率ρ的增高，镜反射面元的辐射传递系数计算效率也得到进一步提升。
图 8 (a) L型腔体内表面反射率ρ=0.4时蒙特卡洛方法与快速方法的计算时间对比; (b) ρ=0.8时的计算时间对比
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是根据面元反射率ρ 得出的，之间的关系为：
$$ {n}=\left\lceil\frac{\lg \eta}{\lg \rho}\right\rceil $$ (6) 式中：$\left\lceil\;\;\right\rceil $表示向上取整；η为能量阈值，上述实验中取值为1‰。
图9中红色曲线为蒙特卡洛方法的射线平均追踪次数随反射率ρ的变化关系。从右侧红色坐标轴可见，当ρ=0.9时，反射次数可达66次之多。
图 9 射线平均追踪次数示意
Figure 9. Illustration of the average number of ray-tracing times
对于快速方法，其非漫反射表面的发射光线的追踪次数和系统中漫反射面元的占比有关：光线击中漫反射面元，追踪进程就结束。假如射线追踪次数为k，单根射线的追踪次数的概率pk可表示为：
$$p_k= \begin{cases}\left(\dfrac{N_2-1}{N-1}\right)^{k-1}\left(\dfrac{N_1}{N-1}\right), & 1 \leqslant k< n \\ \left(\dfrac{N_2-1}{N-1}\right)^{k-1}, & k=n\end{cases}$$ (7) 式中：N、N1、N2分别为系统面元总数、漫反射面元数量、非漫反射面元数量。
因此，非漫反射面元发射光束的平均追踪次数$ \bar{n} $可表示为：
$$ \bar{n}=\sum_{k=1}^{n} p_{k} \cdot k $$ (8) 当面元反射率较大时，非漫反射面元平均追踪次数近似表示为：
$$ \bar{n}_{\max } \approx \frac{N-1}{N_{1}} $$ (9) 以正方体模型为例 (N=6)，射线追踪次数随面元反射率ρ和漫反射面元数量N1的变化关系如图9中黑色曲线所示。
可知，对于非漫反射面元，射线追踪次数随面元反射率ρ缓慢增加，最大反射次数如公式 (9) 所述。而漫反射面元的射线追踪次数随面元反射率ρ急剧升高，为快速方法提升计算效率创造了条件。
表4所示为采用两种方法计算得到的非漫反射面元的射线平均追踪次数比值（快速方法/蒙特卡洛方法）随面元反射率ρ、漫反射面元占总面元的比值Cd的变化关系矩阵。如果单次射线追踪的时间为定值，则追踪次数之比即为镜反射面元的计算耗时之比。
表 4 非漫反射面元的理论追踪次数比值（快速方法/蒙特卡洛方法）
Table 4. Ratio of average-tracking-times for a single ray emitted from the non-diffuse surface calculated with Fast method compared to Monte-Carlo method
ρ Cd 1/6 2/6 3/6 4/6 5/6 0.1 0.813 0.653 0.520 0.413 0.333 0.2 0.672 0.461 0.330 0.250 0.200 0.3 0.615 0.397 0.277 0.208 0.167 0.4 0.520 0.307 0.208 0.156 0.125 0.5 0.446 0.248 0.167 0.125 0.100 0.6 0.341 0.178 0.119 0.089 0.071 0.7 0.247 0.125 0.083 0.063 0.050 0.8 0.161 0.081 0.054 0.040 0.032 0.9 0.076 0.038 0.025 0.019 0.015 从表4可知，对于非漫反射表面，快速方法相对于蒙特卡洛方法，其计算优势随着面元反射率ρ的增加而提高，并随着系统中漫反射面元占比的增加而提升。
2) 辐射传递系数计算总时间
尽管快速方法能够显著提高非漫反射面元的射线追踪效率，但是由于漫反射面元的辐射传递系数计算时间与蒙特卡洛方法持平，系统整体计算耗时提升受漫反射面元占比影响显著。假如单次射线追踪时间为t，则快速方法和蒙特卡洛计算总时间比τ如公式 (10) 所示。
$$ \begin{split} \tau & =\frac{N_1 m n t+N_2 m \bar{n} t}{N m n t} = \frac{N_1}{N}+\left(1-\frac{N_1}{N}\right) \frac{\bar{n}}{n} \\ & =C_d+\left(1-C_d\right) C_n \end{split} $$ (10) 式中：Cd为漫反射面元占模型总面元的比例因子；Cn为非漫反射面元使用两种方法的平均追踪次数之比（快速方法/蒙特卡洛方法）。
表5所示为快速方法和蒙特卡洛方法计算全部面元的辐射传递系数的耗时比值随面元反射率ρ、漫反射面元占比Cd的变化关系矩阵。从表中可以看出，随着反射率ρ的增加，快速方法的计算效率在也在提高，但计算时间比值不会小于系统中漫反射面元的占比。例如，第3列中，计算时间比值最小为0.513，趋近但大于漫反射面元的占比3/6。这是因为无论快速方法中非漫反射面元的计算时间如何降低，计算总耗时中，漫反射面元的那部分计算时间与蒙特卡洛方法相比是持平的。
表 5 全部面元辐射传递系数理论计算总时间之比（快速方法/蒙特卡洛方法）
Table 5. Ratio of consumed time for all the surfaces using Fast method compared with Monte-Carlo method
ρ Cd 1/6 2/6 3/6 4/6 5/6 0.1 0.844 0.769 0.760 0.804 0.889 0.2 0.727 0.641 0.665 0.750 0.867 0.3 0.679 0.598 0.638 0.736 0.861 0.4 0.600 0.538 0.604 0.719 0.854 0.5 0.539 0.499 0.583 0.708 0.850 0.6 0.451 0.452 0.560 0.696 0.845 0.7 0.373 0.417 0.542 0.688 0.842 0.8 0.301 0.387 0.527 0.680 0.839 0.9 0.230 0.359 0.513 0.673 0.836 同时可以看出，随着漫反射面元占比的增加，虽然从上节中分析得到镜反射面元的计算效率得到提升，但由于加入了漫反射面元的计算时间，使得整体的计算时间呈上升趋势，仅在反射率较小时（例如ρ=0.1~0.5时）有所波动。
为验证理论分析的正确性，以6面元立方体模型为例，当面元反射率ρ=0.8 时，通过改变漫反射面元的数量，分别使用两种方法计算内壁面元间的辐射传递系数，获取计算耗时比值，并与理论值进行对比，计算结果如表6所示。可见，上述理论分析结果是可信的。
表 6 快速方法与蒙特卡洛方法计算时间比值的理论值与实测值 (ρ=0.8)
Table 6. Theoretical and measured values of the computation time ratio between Fast method and Monte-Carlo method when ρ=0.8
Cd 1/6 2/6 3/6 4/6 5/6 Theoretical 0.301 0.387 0.527 0.680 0.839 Measured 0.331 0.384 0.524 0.695 0.846 -
文中针对蒙特卡洛方法计算辐射传递系数计算时间过长的问题，提出了一种计算漫反射灰体表面与非漫反射灰体表面间的辐射传递系数的快速方法，并通过案例和理论分析对计算精度和计算效率进行了验证和分析。结果表明：快速方法在保证计算精度的同时，能够有效提高非漫反射面元的辐射传递系数计算效率；面元的反射率越大、系统中漫反射面元的占比越高，非漫反射面元的计算效率优势越明显。同时，系统的辐射传递系数计算总时间随反射率的增大而降低，并受制于系统中漫反射面元的占比。该方法为估计漫反射与非漫反射表面间的辐射传递系数提供了更多的选择。
Fast calculation of radiative heat transfer coefficient between diffuse and non-diffuse surfaces
-
摘要: 针对蒙特卡洛方法计算目标表面间辐射传递系数耗时过长的问题，提出了一种计算漫反射灰体与非漫反射灰体表面间辐射传递系数的快速方法。首先，利用蒙特卡洛方法在漫反射表面生成发射光束集并进行射线追踪，计算漫反射表面对目标表面的辐射传递系数；然后，在非漫反射表面生成发射光束并进行射线追踪，如果射线在多次反射过程中击中任意漫反射表面，则利用反射光束的漫反射特性，直接调用漫反射表面的辐射传递系数对反射能量进行分配，并停止射线追踪。对正六面体模型内表面（两个漫反射面元、四个镜反射面元）的辐射传递系数进行仿真，结果表明，当表面反射率ρ=0.4时，镜反射表面的辐射传递系数计算时间仅为蒙特卡洛方法的1/3；ρ=0.8时，相应的计算时间占比低于1/10，并且保持良好的计算精度。理论推导并分析了影响该方法计算速度的模型参数，结果表明：表面反射率越高，系统中漫反射表面占比越大，快速方法中单束光线的平均追踪次数较蒙特卡洛方法越少，非漫反射表面的辐射传递系数计算优势越明显。Abstract:
Objective Radiative heat transfer is one of the three basic modes of heat transfer and has an important impact on the study of the temperature distribution and infrared radiation characteristics of the outer surface of a space target. For solving the radiation heat transfer between incomplete gray-diffuse surfaces system (both diffuse surfaces and non-diffuse surfaces in the model), there is usually a lack of corresponding analytical solutions. The Monte-Carlo method has the advantage of good calculation accuracy, but it has the disadvantage of long calculation time. In order to solve this problem, this paper proposes a representation of the ray reflection energy of the diffuse surface and a calculation method for the radiative transfer coefficient of the incomplete diffuse surfaces system based on the diffuse reflection characteristics of the diffuse surface. This reduces multiple ray tracing before the ray energy threshold is reached and improves the calculation speed. Methods A method for expressing the reflected energy of diffuse surfaces is proposed. When the Monte-Carlo method is used to conduct ray tracing, if a ray hits a diffuse surface, the reflection energy of the ray is defined as the diffuse emission beam set in the upper half space of the surface using the diffuse reflection characteristics of the diffuse surface (Fig.2(b)). In addition, modifications were made to the Monte-Carlo tracking process. First, the radiative transfer coefficient of the diffuse surface was calculated using the Monte-Carlo method. Then, when calculating the radiative transfer coefficient of diffuse surfaces, if the light emitted from the surface intersects with the diffuse surface, the radiative transfer coefficient of the other surfaces is multiplied by the reflected energy, and the reflected energy absorbed by the other surfaces is calculated to end the ray tracing process, avoiding subsequent multiple ray tracing and achieving fast calculation of diffuse surfaces to improve the computational efficiency of the entire system (Fig.5). Results and Discussions Using the cube model and assuming that the No.1 and No.2 surfaces are diffuse and the rest of the surfaces are specular (Fig.4(a)), the radiative transfer coefficient from each surfaces to the other surfaces (including itself) are calculated using the Monte-Carlo method and the fast algorithm proposed in this paper, and the results of the calculations are shown (Tab.1-2). It can be seen that both of them have the same accuracy, but as for the computation time, the new method is more efficient due to the significant reduction of the number of ray tracing times in Monte-Carlo (Fig.6). In addition, the 13-facet L-type unenclosed cavity model is used as proof to show that the method is also applicable in complex models. Finally, taking the cube model as an example, the advantages of the fast method compared with the Monte-Carlo method are analyzed from the theoretical point of view. For the emitted beam on a non-diffuse surface, the average tracing times of its rays are much smaller than that of the Monte-Carlo method, and the higher the reflectivity of the surface element, the more significant the computational advantages are. For example, if the energy threshold is 0.001, the number of diffuse reflective surfaces in the model is 2. When the surface reflectance is 0.4, the calculation time of the non-diffuse reflective surfaces in the fast calculation method is 0.307 times that of the Monte-Carlo method. When the reflectance of the surface is increased to 0.8, the calculation time of the non-diffuse reflective surface is only 0.081 of that of the Monte-Carlo method, which is more than ten times higher. Conclusions Aiming at the problem of the traditional Monte-Carlo method in calculating the radiative transfer coefficient between diffuse surfaces and diffuse surfaces for a long time, a fast calculation method is proposed. Firstly, the realization principle of the method and the difference with Monte-Carlo method are introduced, and then the cube model and L-shape unenclosed cavity model are used to compare the calculation results and calculation time of the fast method to Monte-Carlo method, which illustrates the advantages of fast method over Monte-Carlo method in terms of the calculation efficiency, and then the coefficient affecting the calculation advantages of the method are illustrated through the theoretical analysis, and finally the outlook on the next step of how to improve the calculation time of the diffuse surface elements is proposed. -
图 1 蒙特卡洛法计算辐射传递系数示意
Figure 1. Monte-Carlo method for calculating radiative heat transfer coefficient
图 2 (a) 蒙特卡洛方法中的漫反射表面反射光线表示；(b) 快速方法中的漫反射表面反射光束表示
Figure 2. (a) The representation of reflected light by diffuse surface in Monte-Carlo method; (b) Fast method
图 3 非漫反射表面的辐射传递系数计算示意图
Figure 3. Fast calculation method for calculating the radiative heat transfer coefficient of non-diffuse surface element i
图 4 (a) 三维几何模型与面元编号；(b) 随机发射光线的位置和方向
Figure 4. (a) Illustration of 3D geometry model and its surface ID; (b) Illustration of random positions and direction vectors for emitted rays
图 5 漫反射与镜反射表面模型辐射传递系数快速计算
Figure 5. Fast calculation of radiative heat transfer coefficient for model with diffuse and specular surfaces
图 6 (a) 立方体内表面反射率ρ=0.4时蒙特卡洛方法与快速方法的计算时间对比; (b) ρ=0.8时的计算时间对比
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
图 8 (a) L型腔体内表面反射率ρ=0.4时蒙特卡洛方法与快速方法的计算时间对比; (b) ρ=0.8时的计算时间对比
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 立方体内表面各面元之间的辐射传递系数（蒙特卡洛方法）
Table 1. Radiative heat transfer coefficient for the inner surfaces of the cube (Monte-Carlo method)
Facet 1(AA'B'B) 2(BB'C'C) 3(CC'D'D) 4(DD'A'A) 5(A'B'C'D') 6(ABCD) Integrality 1 0.05897 0.17917 0.21741 0.18131 0.18132 0.18116 0.99934 2 0.17956 0.05924 0.18051 0.21850 0.18032 0.18122 0.99934 3 0.21776 0.18120 0.04888 0.18594 0.18307 0.18249 0.99934 4 0.18201 0.21647 0.18750 0.04936 0.18190 0.18212 0.99934 5 0.17921 0.17974 0.18220 0.18333 0.07049 0.20437 0.99934 6 0.17872 0.17926 0.18274 0.18213 0.20545 0.07105 0.99934 
下载: 导出CSV
表 2 立方体内表面各面元之间的辐射传递系数（快速方法）
Table 2. Radiative heat transfer coefficient for the inner surfaces of the cube (Fast method)
Facet 1(AA'B'B) 2(BB'C'C) 3(CC'D'D) 4(DD'A'A) 5(A'B'C'D') 6(ABCD) Integrality MAE 1 0.05937 0.18057 0.21859 0.18222 0.17886 0.17975 0.99934 1.29227‰ 2 0.18068 0.05972 0.18135 0.21884 0.17903 0.17972 0.99934 0.93054‰ 3 0.22136 0.18010 0.04669 0.18631 0.18214 0.18324 0.99985 1.48970‰ 4 0.18218 0.22024 0.18636 0.04688 0.18098 0.18323 0.99985 1.59877‰ 5 0.17933 0.18143 0.18179 0.18329 0.06566 0.20834 0.99985 1.84531‰ 6 0.18179 0.18012 0.18301 0.18145 0.20777 0.06570 0.99985 2.09184‰
下载: 导出CSV
表 3 L型腔计算误差与能量完整性（快速方法）
Table 3. MAE and integrality of radiative heat transfer coefficient for L-shape cavity model (Fast method)
Facet ρ=0.4 ρ=0.8 MAE Integrality MAE Integrality 1 1.11640‰ 0.98594 0.61132‰ 0.85933 2 1.29494‰ 0.94656 2.29617‰ 0.77797 3 0.55498‰ 0.79605 1.91049‰ 0.64609 4 1.19067‰ 0.79290 1.34883‰ 0.64261 5 6.59430‰ 0.97410 4.82850‰ 0.80405 6 0.98649‰ 0.98910 0.64184‰ 0.85692 7 0.64990‰ 0.93160 0.56992‰ 0.76133 8 1.33612‰ 0.65359 1.32864‰ 0.41103 9 5.29473‰ 0.94401 5.60488‰ 0.77703 10 1.45282‰ 0.79286 1.44651‰ 0.65081 11 1.03800‰ 0.79272 1.53688‰ 0.64366 12 2.62882‰ 0.94494 3.66704‰ 0.77457 13 5.80383‰ 0.98413 4.92681‰ 0.84846
下载: 导出CSV
表 4 非漫反射面元的理论追踪次数比值（快速方法/蒙特卡洛方法）
Table 4. Ratio of average-tracking-times for a single ray emitted from the non-diffuse surface calculated with Fast method compared to Monte-Carlo method
ρ Cd 1/6 2/6 3/6 4/6 5/6 0.1 0.813 0.653 0.520 0.413 0.333 0.2 0.672 0.461 0.330 0.250 0.200 0.3 0.615 0.397 0.277 0.208 0.167 0.4 0.520 0.307 0.208 0.156 0.125 0.5 0.446 0.248 0.167 0.125 0.100 0.6 0.341 0.178 0.119 0.089 0.071 0.7 0.247 0.125 0.083 0.063 0.050 0.8 0.161 0.081 0.054 0.040 0.032 0.9 0.076 0.038 0.025 0.019 0.015
下载: 导出CSV
表 5 全部面元辐射传递系数理论计算总时间之比（快速方法/蒙特卡洛方法）
Table 5. Ratio of consumed time for all the surfaces using Fast method compared with Monte-Carlo method
ρ Cd 1/6 2/6 3/6 4/6 5/6 0.1 0.844 0.769 0.760 0.804 0.889 0.2 0.727 0.641 0.665 0.750 0.867 0.3 0.679 0.598 0.638 0.736 0.861 0.4 0.600 0.538 0.604 0.719 0.854 0.5 0.539 0.499 0.583 0.708 0.850 0.6 0.451 0.452 0.560 0.696 0.845 0.7 0.373 0.417 0.542 0.688 0.842 0.8 0.301 0.387 0.527 0.680 0.839 0.9 0.230 0.359 0.513 0.673 0.836
下载: 导出CSV
表 6 快速方法与蒙特卡洛方法计算时间比值的理论值与实测值 (ρ=0.8)
Table 6. Theoretical and measured values of the computation time ratio between Fast method and Monte-Carlo method when ρ=0.8
Cd 1/6 2/6 3/6 4/6 5/6 Theoretical 0.301 0.387 0.527 0.680 0.839 Measured 0.331 0.384 0.524 0.695 0.846
下载: 导出CSV
-
[1] 曹飞飞, 吉洪湖, 于明飞, 等 . 低发射率材料涂敷区域对排气系统壁温和红外特性的影响 [J]. 红外与激光工程,2020 ,49 (10 ):20190131 . Cao Feifei, Ji Honghu, Yu Mingfei, et al. Effects of low emissivity material coating site on wall temperature and infrared characteristics of exhaust system [J]. Infrared and Laser Engineering, 2020, 49(10): 20190131. (in Chinese)[2] 汪洪源, 陈赟 . 天基空间目标红外动态辐射特性建模与仿真 [J]. 红外与激光工程,2016 ,45 (5 ):0504002 . doi: 10.3788/irla201645.0504002 Wang Hongyuan, Chen Yun. Modeling and simulation of infrared dynamic characteristics of space-based space targets [J]. Infrared and Laser Engineering, 2016, 45(5): 0504002. (in Chinese) doi: 10.3788/irla201645.0504002[3] Downer J, Kabir M, Xu J. The design and development of a smart multilayer coating with variable emissivity capability for space vehicle thermal control systems[C]//22nd IEEE Intersociety Conference on Thermal and Thermomechanical Phenomena in Electronic Systems (ITherm), 2023: 1-12. [4] 刘硕, 朱希安, 厉夫兵 . 弹头包络目标辐射传递系数的蒙特卡洛计算 [J]. 战术导弹技术,2018 (3 ):44 -52 . Liu Shuo, Zhu Xi'an, Li Fubing. Monte Carlo method for calculating radiative heat transfer coefficient of warhead envelope target [J]. Tactical Missile Technology, 2018(3): 44-52. (in Chinese)[5] Gebhart B. Surface temperature calculations in radiant surroundings of arbitrary complexity-for gray, diffuse radiation [J]. International Journal of Heat and Mass Transfer, 1961, 3(4): 341-346. doi: 10.1016/0017-9310(61)90048-5 [6] 杨世铭, 陶文铨. 传热学[M]. 北京: 高等教育出版社, 2006. [7] 卞伯绘. 辐射换热的分析与计算[M]. 北京: 清华大学出版社, 1988. [8] 张伟清, 宣益民, 韩玉阁 . 单元表面间辐射传递系数的新型计算方法 [J]. 宇航学报,2005 ,26 (1 ):77 -80+85 . Zhang Weiqing, Xuan Yimin, Han Yuge. New method of radiative transfer coefficient between unit surfaces [J]. Journal of Astronautics, 2005, 26(1): 77-80+85. (in Chinese)[9] He F, Shi J, Zhou L, et al. Monte Carlo calculation of view factors between some complex surfaces: rectangular plane and parallel cylinder, rectangular plane and torus, especially cold-rolled strip and W-shaped radiant tube in continuous annealing furnace [J]. International Journal of Thermal Sciences, 2018, 134: 465-474. doi: 10.1016/j.ijthermalsci.2018.05.050 [10] Narayanaswamy A. An analytic expression for radiation view factor between two arbitrarily oriented planar polygons [J]. International Journal of Heat & Mass Transfer, 2015, 91: 841-847. [11] Li H, Wang S, Yuan D, et al. Electromagnetic-thermal simulation of cold crucible considering surface-to-surface heat radiation[C]//IEEE International Conference on Applied Superconductivity and Electromagnetic Devices (ASEMD), 2015: 368-369. [12] Howell J R, Perlmutter M. Monte Carlo solution of thermal transfer through radiant media between gray walls [J]. Journal of Heat Transfer, 1964, 86(1): 116-122. doi: 10.1115/1.3687044 [13] Knoll A, Wald I, Parker S, et al. Interactive isosurface ray tracing of large octree volumes[C]//IEEE Symposium on Interactive Ray Tracing, 2006: 115-124. [14] Gro M, Lojewski C, Bertram M, et al. Fast implicit KD-trees: accelerated isosurface ray tracing and maximum intensity projection for large scalar fields[C]//Proceedings of the Ninth IASTED International Conference on Computer Graphics and Imaging, 2007: 67-74. [15] Osada R, Funkhouser T, Chazelle B, et al. Matching 3D models with shape distributions[C]//Proceedings International Conference on Shape Modeling and Applications, 2001: 154-166. [16] 厉夫兵, 冷俊敏, 李红莲 . 诱饵球镜面反射内壁的辐射传递系数计算 [J]. 现代防御技术,2017 ,45 (3 ):193 -199+214 . doi: 10.3969/j.issn.1009-086x.2017.03.030 Li Fubing, Leng Junmin, Li Honglian. Radiative heat transfer coefficient calculation of inflatable sphere with specular inner surface [J]. Modern Defence Technology, 2017, 45(3): 193-199, 214. (in Chinese) doi: 10.3969/j.issn.1009-086x.2017.03.030 -
点击查看大图
计量
- 文章访问数: 34
- HTML全文浏览量: 5
- PDF下载量: 23
- 被引次数: 0
