基于改进蒙特卡洛法求解辐射传热问题的计算精度评价

李国军 钟佳琪 李定勇 王晓东

摘要：针对传统蒙特卡洛法计算辐射传输耗时问题，提出了一种改进蒙特卡洛方法，通过比例迭代累加法来求解反射及散射能量，从而大幅减少了计算时间.引入直接评价方法，以包含参与性介质的密闭系统（方形和圆形为例）为例，分析了网格密度、发射能束数及物性参数对改进蒙特卡洛法计算精度的影响.当光学厚度为0.005时，采用改进蒙特卡洛方法求得方腔及圆形腔的表面微元辐射通量的相对均方根误差值分别为0.0025及0.0023，而采用传统蒙特卡洛方法时对应误差分别为0.0080及0.0037.可见，相同计算条件下，改进蒙特卡洛方法对辐射换热问题具有更高的精度.进一步研究了追踪能束数对计算误差的影响规律，给出了计算误差与追踪能束数拟合关系式，为计算该类问题的能束数选取提供了支撑.

关键词：蒙特卡洛法;评价;计算精度;计算效率

中图分类号：TK124

文献标志码：A

能源动力、航空航天等领域的传热过程通常以辐射传热为主，大多数情况该类问题难于进行解析求解而多借助于数值方法进行求解，故而对该类问题的数值计算方法研究具有重要意义[1-5].蒙特卡洛方法（Monte Carlo method，MCM）可以精确地处理光谱特性、非均匀介质、各向异性散射及复杂几何形状等复杂辐射计算问题，已经成为解决辐射传热的主要数值计算方法之一[6-9].MCM是一种统计模拟方法，其数值结果的精度随着抽样能束数目的增加而提高，而计算时间随着追踪能束数的增加而增加.为兼顾计算精度及计算时间，笔者提出了一种改进MCM[10]，该方法可以在只进行一次抽样情况下完成辐射换热求解，且其计算精度及计算效率高于传统MCM.随着MCM在辐射传热数值模拟计算应用范围的不断扩展，如何定量分析和评估其计算结果的误差及精度已成为关注的焦点，建立公认的数值误差分析和精度评估方法成为MCM主要研究内容之一.Siegel等[11]首次用不确定度来估计MCM计算辐射传热问题的统计误差.Planas Almazan[12]运用MCM射线路径轨迹量化了混合网交换公式固有的统计误差.Plotnikov和Shkarupa[13]应用直接模拟MCM求解稀薄气体动力学问题的统计误差.阮立明等[14]利用辐射交换因子的守恒性和互易性的检验来评估MCM的计算误差，得到一种间接评价MCM计算精度的方法.Yarahmadi等[15]提出一种改进平均相对不确定评价方法，验证由于表面温度和发射率的不确定性而引起的净热流密度不确定度的新表达式.Wang等人[16-17]提出了一种直接定量评价MCM精确度的评价方法.

为评价改进MCM计算精度，本文拟采用改进MCM分别对漫灰表面、参与性各项同性介质的方形及圆形封闭腔体内的辐射传热问题进行研究，建立表面和体积微元辐射热通量的误差计算模型，得出其最小误差与能束数的函数关系.采用直接定量评价方法开展基于辐射热通量的改进MCM计算误差的分析和精度评价，研究网格密度及采样能束数变化对改进MCM求解辐射传热计算精度及效率的影响.

1改进蒙特卡洛法简介

应用区域法求解辐射换热时，当表面段和体积段的参数确定后，各段之间的辐射传递因子也可通过计算获得.若假定微元发射的全部能量到达其他各段能量的比例与该微元反射能量到达其他各段能量的比例相同，与发射能量份额及反射能量的份额无关，则采用MCM求解辐射传递因子时，对全部微元发射的能束只需进行一次采样追踪，以确定微元段发射能量到达其他微元段的比例.当表面微元反射时，将反射能束按发射能束处理，且到达其他段的比例在之前已经确定，无需再计算，散射情况也类似处理.该思路即为改进MCM，具体求解方法如下.

1.1微元发射能束直接到达其他微元段比例

为方便计算，将表面微元按顺序依次命名为1，2，...，Ns，体积微元安排在表面微元之后为Ns+1，Ns+2，...，Ns+Ng，其中Ng和Ns分别为体积和表面微元总数，则总微元数为Ns+Ng.则热交换场中微元i发射能束直接到达微元j比例为

式中：kUi，（jk=1，2，3，...）表示第i个微元发射的能量第k次循环到达第j个微元的能束数，Ni为第i段发射的总能束数，当k=1时表示直接到达.

1.2微元发射能束到达其他微元段的总能束数

第i个微元发射的能量最终到达第j个微元的能量Ui，j分为直接到达被吸收的能量及经k次反射后到达被吸收的∑能量的累加，则

式中：Ui，j为第i微元段发射能束到达第j微元段能束总数，其中包含直接到达与反射到达情况，ε为微元段黑度.

对体积微∑元j=Ns+1，Ns+2，...，Ns+Ng有

式中：ω为散射反照率，U达第j微元段能束总数，其中包含直接到达与散射到达情况.

当介质为各项同性散射介质时，

式中：K。為含粒子介质系的衰减系数，Km为含粒子介质系的粒子的衰减系数，ω表示散射介质的散射反照率.

首先确定一个随机数R，如果R》ω，则能束被吸收，否则被散射.

确定散射方向是MCM研究含粒子系辐射传递的关键.本文仅计算各向同性散射，已知散射相函数的归一化条件.

对于各项同性散射的散射相函数表达为

对于各项同性散射，即b=0时

1.3计算终止条件

当微元段i发射的能束经过多次的反射或散射后剩余能量逐渐减少，当剩余能量与发射能量的比值满足式（10）时，式中ξ为无穷小量，则认为计算精度已经∑满足要求.

2直接定量评价法

将辐射体系划分为M个面元和N个体积微元，则表面和体积微元净辐射热通量

因计算过程中存在误差，则其真值可以描述为计算值与计算误差之和，即

式中：Fm→i为表面微元m对表面微元i的辐射交换因子;Fn→j为体积微元n对体积微元j的辐射交换因子;ε是表面微元i的黑度;κa是气体光学厚度;A、V分别为表面微元的面积和体积微元的体积;q、q分别为表面微元i、体积微元j上的净辐射热通量;q00算值;Δqai、Δqvj为计算误差.

在等温辐射平衡状态下，净辐射通量理论上为0，则微元段吸收的辐射能量与发射的辐射能量理论上绝对相等，从而得到表面微元无因次方程为

由式（11）（13）（15），可得表面微元i的无因次净辐射热流与精确值的偏差

同理可得，体积微元j的无因次净辐射热流与精确值的偏差

则表面微元和体积微元的辐射热通量相对误差为

表面和体积微元辐射通量的相对均方根误差表示为

式中：Q0ai、Q0vj分别为表面微元i、体积微元j上吸收无因次热流的计算值;ΔQai、ΔQvj分别为表面微元i、体积微元j上吸收无因次热流的计算误差.

3结果与讨论

3.1网格密度对计算误差的影响以漫灰表面、参与性各向同性介质的方形及圆

形封闭腔体内的辐射传热问题为例，如图1所示，将边长为L的方腔离散按长度方向Nx和高度方向Ny均匀划分网格;将半径为R的圆形腔体沿径向Nr及周向Nz划分为均匀网格.本节中引入光学厚度τg=κa⋅（LNx）并分析网格划分密度对计算误差影响.计算设定条件为：设方腔离散网格数为N=N，=5，10，20，40，圆形腔体离散网格数为N=4，14，28，40，N=2.微元发射能束N=10，表面发射率s=0.5，介质的散射反照率ω=0.5.为减小伪随机数的随机性对计算结果的影响，本文采用2次计算结果及其平均值进行研究.

由于伪随机数的随机性，不同的离散网格密度的计算结果随机误差均值表现出不同的噪声.以图2（a）（b）（c）（d）和图2（e）（f）（g）（h）所示，离散网格密度较小则误差波动明显，离散网格密度越大，结果越稳定，并且离散网格密度的增加使改进MCM计算结果的随机性和任意性最小化到更小的程度.从图2可以看出，方形网格数取Nx=Ny=40，圆形网格数取Nz=40，Nr=2即可满足精度条件.

3.2辐射传热量的精确值和计算值

以漫灰表面、参与性各向同性介质的方形及圆形封闭腔体内的辐射传热问题为例，将方形离散成均匀表面微元和体积微元，其中Nx=Ny=40.将圆形腔体沿径向及周向分别划分为均匀网格，即Nz=40，Nr=2.网格划分如图1所示.

当微元发射能束数N=105，表面发射率ε方=0.8，ε圆=0.5，介质的散射反照率ω=0.5，气体光学厚度分别为τg=0.0005，τg=0.5，τg=50时，表面微元无因次热通量计算结果分别如图3所示.

由图3可知，光学厚度不同時，分别采用改进MCM及MCM计算方腔及圆形腔内表面微元无因次热通量的计算值与真实值误差较小.当τg=0.005时，采用改进MCM求解得到的方腔及圆形腔的相对均方根误差Ea分别为0.0025及0.0023（MCM求解得到的方腔及圆形腔的相对均方根误差Ea分别为0.0080及0.0037）;当τg=0.5时，其改进MCM求解的Ea分别为0.0029及0.0028（MCM求解的Ea分别为0.0141及0.0116）;当τg=50时，其改进MCM求解的Ea分别为0.0032及0.0031（MCM求解的Ea分别为0.0270及0.0238）.由此可以看出，改进MCM对方腔及圆形腔具有相同适用性，在相同条件下其求解精度高于MCM.

图4给出了气体光学厚度分别为τg=50，τg=0.5，τg=0.0005时方腔及圆形腔体积微元辐射热通量计算值与真实值关系.由图可知，当τg=0.0005时，用改进MCM求解得到的方腔及圆形腔的相对均方根误差Ev分别为0.00138及0.0268（MCM求解的Ev分别为0.0403及0.0978）;当τg=0.5时，方腔及圆形腔的相对均方根误差Ev分别为0.0142及0.0512（MCM求解的Ev分别为0.0403及0.1207）;当τg=50时，方腔及圆形腔的相对均方根误差Ev分别为0.4997及0.1462（MCM求解的Ev分别为1.0018及0.4782）.可见，当其他条件不变时，随着光学厚度增加，MCM及改进MCM计算结果误差有增大趋势.

3.3微元发射能束数与计算误差关系

设方形及圆形的表面发射率ε=0.5，介质的散射反照率ω=0.5，当微元发射能束数N分别为N=1000，3000，5000，10000，100000，500000时，求解表面和体积微元无因次热通量最小误差值.

表1给出了Ea，min、Ev，min与微元发射能束数之间的关系.可以看出，随着微元发射能束数增加，最小误差逐渐减小.对圆形封闭腔，无因次热通量最小误差值Ea，min与微元发射能束数N的关系可拟合为如式（22）所示.

由图5可知，计算精度随着微元能束数的增加而提高，但是当达到一定的微元能束数时，能束数继续增加对计算精度提升逐渐不明显.通常对于燃烧模拟，辐射传热计算中可接受的Ea，min误差水平要求约为0.01，所以对于用改进MCM求解辐射传热的方形和圆形算例分别需要微元能束数N=5000和N=10000即为满足精度要求.通过式（23）可以获得不同最小误差要求情况下需要的微元能束数目.

4结论

本文阐述了基于改进MCM进行求解辐射换热问题的思路，并引入直接定量评价方法对不同光学厚度情况下的圆形腔及方腔内参与性介质的辐射热通量计算精度进行了评价.研究结果如下：

1）改进MCM对方腔及圆形腔内辐射问题具有适用性，在相同计算条件下，其直接评价结果明显高于采用改进MCM进行计算的精度.如对于光学厚度为0.005时，对方腔及圆形腔，采用改进MCM的Ea值分别0.0025及0.0023而采用MCM时Ea分别为0.0080及0.0037.

2）研究了微元发射能束数与计算误差关系，给出了能束数与Ea关系拟合式，为采用改进MCM进行辐射换热计算提供追踪能束数选取的依据.

参考文献

[1]李国军，陈海耿，伊智.积分降重法求解辐射直接交换面积[J].东北大学学报（自然科学版），2010，31（1）：80-83.

[2] MODEST M F. Radiative heat transfer[M]. New York：Academic Press，1993：639-666.

[3] RUAN L M，TAN H P，YAN Y Y. A Monte Carlo（mc）method applied to the medium with nongray absorbing-emitting- anisotropic scattering particles and gray approximation[J].Nu⁃merical Heat Transfer，Part A：Applications，2002，42（3）： 253-268.

[4] TAN H P，SHUAI Y，XIA X L. Reliability of stray light calcula⁃tion code by the Monte Carlo method[J]. Optical Engineering， 2005，44（2）：1-2.

[5] KOVTANYUK A E，BOTKIN N D，HOFFMANN K H.Numerical simulations of a coupled radiative-conductive heat transfer model using a modified Monte Carlo method[J].International Journal of Heat and Mass Transfer，2012，55（4）：649-654.

[6] TSENG H Y，HONG Y B，WU C Y，et al.Monte Carlo simulation of radiative transfer in a medium with varying refractive index specified at discrete points[J].Applied Mathematical Modelling， 2 0 1 7 ，4 8 ：8 7 0 - 8 8 4 .

[7] WANG C H，ZHANG Y，YI H L，et al.Transient radiative trans⁃fer in two-dimensional graded index medium by Monte Carlo method combined with the time shift and superposition principle [J]. Numerical Heat Transfer，Part A：Applications，2016，69 （ 6 ）：5 7 4 - 5 8 8 .

[8]龔光彩，刘佳.空气载能辐射空调混合通风协同运行研究[J].湖南大学学报（自然科学版），2019，46（5）：148-156.

[9] KHANNA S，YADAV U. Compution of total exchange ares using Monte Carlo method[J]. Journal of Engineering Thermophysics， 2013，8（3）：1-7.

[10] LI G J，ZHONG J Q，WANG X D. An improved Monte Carlo method for radiative heat transfer in participating media[J].Jour⁃nal of Quantitative Spectroscopy and Radiative Transfer，2020， 2 5 1 ：1 0 7 0 8 1 .

[11] SIEGEL R，HOWELL J R. Thermal radiation heat transfer[M]. NewYork：Taylor&Francis，2002：23-566.

[12] PLANAS ALMAZAN P. Accuracy of Monte Carlo ray-tracing ther⁃mal radiation calculations：A practical discussion[J]. Noordwijk： European Space Agency，1997：579-591.

[13] PLOTNIKOV M Y，SHKARUPA E V.Estimation of the statisti⁃cal error of the direct simulation Monte Carlo method[J].Compu⁃tational Mathematics and Mathematical Physics，2010，50（2）： 335-344.

[14]阮立明，谭建宇，董士奎，等.辐射传递蒙特卡洛法精度分析及数值试验[J].工程热物理学报，2003，24（5）：813-816.

[15] YARAHMADI M，MAHAN J R，PRIESTLEY J. Estimation and use of the radiation distribution factor median for predicting uncer⁃tainty in the Monte Carlo ray-trace method[J]. Journal of Heat Transfer，2019，141（6）：1- 7.

[16] WANG D D，ZHOU H C.Quantitative evaluation of the computa⁃tional accuracy for the Monte Carlo calculation of radiative heat transfer[J].Journal of Quantitative Spectroscopy and Radiative T r a n s f e r ，2 0 1 9 ，2 2 6 ：1 0 0 - 1 1 4 .

[17] WANG D D，CHENG Q，ZHUO X S，et al. Effects of surface emissivity and medium scattering albedo on the computational ac⁃curacy of radiative heat transfer by MCM[J].Journal of Quantita⁃tive Spectroscopy and Radiative Transfer，2020，240：106712.