蒙特卡罗法求解二维各向异性导热-辐射耦合传热

Monte Carlo Method for Solving Two-Dimensional Coupled Anisotropic Conduction-Radiation Heat Transfer

摘要

为了探究二维各向异性材料中的导热-辐射耦合传热特性，基于随机游走的基本思想，在统一离散网格的基础上，开发了一种求解二维各向异性材料的导热-辐射耦合传热的计算方法，并验证了计算方法的准确性，分析了各向异性导热系数对传热特性的影响.结果表明，对于二维平板结构，当存在各向异性导热系数时，会导致平板温度场向相应的各向异性方向发生偏移；当k_{\mathrm{12}}^{*}与k_{\mathrm{22}}^{*}比值相同时，平板的各向异性基本保持不变；而当k_{\mathrm{22}}^{*}保持不变而k_{\mathrm{12}}^{*}增加时，会导致平板内温度传导的方向发生改变，导致平板内温度的均匀性降低.

Abstract

In order to study the coupled thermal conductivity-radiation heat transfer characteristics of two-dimensional anisotropic materials， a computational method was developed for solving the coupled conduction-radiation heat transfer of two-dimensional anisotropic materials based on the basic idea of random walk for a uniform discrete mesh，and the accuracy of the computational method was verified. The effects of anisotropic thermal conductivity coefficient on the heat transfer characteristics were analyzed. The results show that for a two-dimensional plate， an anisotropic thermal conductivity coefficient results in a shift of the temperature field of the plate towards the corresponding anisotropy. The anisotropy of the plate is basically unchanged when the ratio of k_{\mathrm{12}}^{*} to k_{\mathrm{22}}^{*} is the same. The increase of k_{\mathrm{12}}^{*} leads to the change of the direction of the temperature conduction in the plate when the k_{\mathrm{22}}^{*} remains constant， which makes the temperature homogeneity decrease in the plate.

Graphical abstract

关键词

各向异性 / 蒙特卡罗法 / 统一网格 / 耦合传热 / 随机游走算法

Key words

anisotropy / Monte Carlo method / uniform grid / coupled heat transfer / random walk algorithm

导热-辐射耦合传热在半透明材料、太阳能利用、航空航天、加热炉等领域具有广泛的应用前景^［1-4］.各向异性材料因其具有方向性的传热特性和良好的力学性能被广泛应用^［5］，例如多层不同材料的结构体^［6］、高温炉内的观测镜头等.因此，对各向异性材料中的导热-辐射耦合传热过程进行研究具有十分重要的意义.导热-辐射耦合传热过程主要由能量方程和辐射传递方程的耦合求解得到.目前，针对耦合传热问题的数值计算方法主要包括两种类型：一种是采用两种方法分别求解能量方程和辐射传递方程，两者结合来求解导热-辐射耦合传热问题^［7-8］；另一种则是深入探究同一计算方法求解导热-辐射耦合传热方程^［9-10］.其中采用有限体积法^［11］（FVM）求解能量方程和辐射传递方程，与建立统一格子玻尔兹曼框架^［12］（LBM）分析耦合传热问题的方式应用十分广泛.对于某些特定的工作条件，每种方法都有其自身的优点和缺点^［13］.在上述方法中，辐射和导热部分都需要进行复杂的数值计算或大量的迭代过程.实际应用涉及可变的物理性质和复杂的边界条件，导致难以获得合理的数值解.

蒙特卡罗法^［14］（MCM）作为一种随机性方法，因其求解各种复杂形状问题的高精度和灵活性以及耦合求解能力强，已成为目前应用最广泛的数值方法之一.Howell等^［15］将蒙特卡罗法应用于计算灰体辐射热传输问题.现如今通过消除表面反射采样等方式^［16-17］，蒙特卡罗法的计算效率得到了明显的提升.蒙特卡罗法因其随机性，被广泛用于模拟电磁能、电子、声子在吸收和散射介质中的传输过程^［18］，使得计算结果更加符合实际条件.Haji-Sheikh等^［19］基于蒙特卡罗法研究了一种可应用于求解热传导问题的随机游走方法.Kowsary等^［20-21］运用有限差分格式的蒙特卡罗法求解了含有各向异性材料的稳态和瞬态热传导问题，研究了各向异性热传导蒙特卡罗法的计算条件及计算效率.Naeimi等^［22］采用有限体积法对随机游走蒙特卡罗法进行了优化，并将其应用于求解复杂几何结构和不同边界条件下的传热问题.

综上所述，蒙特卡罗法多被用来单独求解辐射传递问题或导热问题，很少有研究将其应用于导热-辐射耦合传热研究.而关于各向异性材料的研究，大多只针对各向异性热传导或者各向异性散射介质的热辐射传输这两方面进行深入探究，对于各向异性材料耦合传热的研究较少.因此，本文针对二维各向异性材料，基于蒙特卡罗法的思想，运用热平衡法建立统一的随机游走导热-辐射耦合传热计算模型，深入研究了二维平板内各向异性对耦合传热的影响.

1 数学模型

1.1 随机游走原理

应用蒙特卡罗法求解耦合传热问题，需要构造相应的随机游走概率模型.对于稳态传热过程，二维导热-辐射耦合传热能量方程为

k_{\mathrm{11}}\frac{{\partial }^{\mathrm{2}}T}{\partial x^{\mathrm{2}}}+k_{\mathrm{22}}\frac{{\partial }^{\mathrm{2}}T}{\partial y^{\mathrm{2}}}+k_{\mathrm{12}}\frac{{\partial }^{\mathrm{2}}T}{\partial x\partial y}+k_{\mathrm{21}}\frac{{\partial }^{\mathrm{2}}T}{\partial x\partial y}+Q=\mathrm{0}.(1)

式中：T为温度，K；k_{\mathrm{11}}，k_{\mathrm{12}}，k_{\mathrm{21}}，k_{\mathrm{22}}为各个方向的导热系数，W/（m·K）；Q为单位体积辐射热源，W/m³，表示为辐射热通量密度的散度.

Q=-\mathrm{d}\mathrm{i}\mathrm{v}{\vec{\mathit{q}}}_{\mathrm{r}}.(2)

式中，{\vec{\mathit{q}}}_{\mathrm{r}}为辐射热通量密度，W/m².

对于粒子的随机游走过程，考虑变换域中的二维几何模型，质点在二维坐标系的离散节点间根据随机移动规则进行游走，通过不断碰撞产生不规则的迁移运动，最终到达边界并获取边界信息.因此，本文采用热平衡法对式（1）中各项进行离散处理.如图1所示，根据Fourier定律可知，W侧节点通过界面传导到节点M（i，j）的热流量{\Psi }_{\mathrm{W}}为

{\Psi }_{\mathrm{W}}=k_{\mathrm{11},\mathrm{W}}∆y\frac{T_{i-\mathrm{1},j}-T_{i,j}}{∆x}.(3)

式中，k_{\mathrm{11},\mathrm{W}}为W侧导热系数，W/（m·K）.参照式（3）可以得出如图1所示的E，S，N方向热流量方程.

对于式（1）中的混合偏导数，正向节点的离散格式为

{\Psi }_{\mathrm{N}\mathrm{E}}=k_{\mathrm{12},\mathrm{N}\mathrm{E}}\left(∆x\frac{T_{i+\mathrm{1},j+\mathrm{1}}-T_{i+\mathrm{1},j}}{∆y}-∆y\frac{T_{i,j+\mathrm{1}}-T_{i,j}}{∆x}\right).(4)

由于本文主要研究各向异性对导热-辐射耦合传热的影响，因此考虑对耦合方程进行适当简化，即假设各向异性导热系数在各自方向上不随温度变化，得到微元节点M （i， j）的能量守恒方程为

{\Psi }_{\mathrm{W}}+{\Psi }_{\mathrm{E}}+{\Psi }_{\mathrm{N}}+{\Psi }_{\mathrm{S}}+{\Psi }_{\mathrm{N}\mathrm{E}}+{\Psi }_{\mathrm{q}}=\mathrm{0}.(5)

式中，{\Psi }_{\mathrm{q}}为节点处的辐射热流量.

将式（3）~式（5）代入式（1），可得

T_{i,j}=P_{\mathrm{W}}{T}_{i-\mathrm{1},j}+P_{\mathrm{S}}{T}_{i,j-\mathrm{1}}+P_{\mathrm{E}}{T}_{i+\mathrm{1},j}+P_{\mathrm{N}}{T}_{i,j+\mathrm{1}}+P_{\mathrm{N}\mathrm{E}}{T}_{i+\mathrm{1},j+\mathrm{1}}+Q_{\mathrm{d}\mathrm{s}}.(6)

式中，节点离散方程的温度系数及离散源项如表1所示，其中P为温度系数总函数，表示为

P=\frac{\mathrm{2}{k}_{\mathrm{11}}}{∆x^{\mathrm{2}}}+\frac{\mathrm{2}{k}_{\mathrm{22}}}{∆y^{\mathrm{2}}}-\frac{k_{\mathrm{12}}+k_{\mathrm{21}}}{∆x∆y}.(7)

由表1可知，温度系数的总和为1，因此满足作为随机移动概率值的条件.值得注意的是，表 1需要保证P值非负，同时考虑结构网格的传热过程，可以得到参数设定需满足的条件为

(8)

其中，表1中的温度系数即为随机游走系数，表示随机游走质点向不同方向移动的概率.对于各向异性热传导过程中质点的随机移动，需要生成（0，1）区间上的均匀随机数R，并根据如下概率区间进行移动：

(9)

随机游走是通过一组随机概率系数P_m 来决定质点的随机移动方向.如图2所示，质点每行走一步都要从一个节点根据该概率系数移动到下一个节点，且质点的每一步移动都是相互独立的.

本文考虑适用范围广泛的Dirichlet边界条件，即当质点移动到边界时随机游走终止，质点被边界吸收.节点M（x，y）的蒙特卡罗温度估计即为质点根据概率系数随机移动后获取到所有边界信息的算术平均值.当质点的数量m足够多时，得到的数值解基本接近真实值.因此，节点M（x，y）的温度为

T_{x,y}={\sum }_m{P}_m{T}_m+Q_{\mathrm{d}\mathrm{s}}.(10)

1.2 蒙特卡罗法求解辐射源项

对于本文所研究的二维平板内的导热-辐射耦合传热，可直接利用导热过程中所划分的质点随机移动网格.假设平板内部存在N个体元节点和4个面元节点，并且在应用蒙特卡罗法求解的过程中需要引入辐射传递因子^［23］K（i，j）的概念：微元体i辐射的能量经过系统内其他微元一次或多次的折射、反射后到达微元体j的份额.

若微元为体积微元\mathrm{d}V：

K^{*}\left(i,j\right)=\mathrm{4}{n}^{\mathrm{2}}{\kappa }_{\mathrm{a}}\mathrm{d}{V}_{i}K(i,j).(11)

若微元为面积微元\mathrm{d}A：

K^{*}\left(i,j\right)=\epsilon \mathrm{d}{A}_{i}K(i,j).(12)

式中：K^*（i，j）为辐射全交换面积；n为介质折射率；{\kappa }_{\mathrm{a}}为吸收系数，m^-1；ε为表面发射率.

能量束的发射主要分为两种情况：体积微元和面积微元.在二维矩形平板内，体积微元的光束发射方向为

{\theta }_{\mathrm{v}}=\mathrm{2}\mathrm{\pi }{R}_{\theta }.(13)

式中：θ_v为体积微元射线的发射方向；R_θ 为（0，1）区间内均匀生成的发射随机数.

值得注意的是，虽然上述提到表面为4个面元节点发射，但考虑表面随机微元发射，且当发射能量束数目足够多时，可近似表示为整体壁面发射.表面随机发射位置为

\left.\begin{array}{c}{L}_{\mathrm{R}x}=L_x{R}_{\mathrm{A}},\\ L_{\mathrm{R}y}=L_y{R}_{\mathrm{A}}.\end{array}\right\}(14)

式中：L_x，L_y 为矩形腔的边长，m；L_Rx，L_Ry 为随机确定的发射位置；R_A为区间（0，1）内均匀生成的发射随机数.

根据光束在传播过程中的连续衰减规律，对其进行随机取样以确定随机发射距离r.光束的传播距离为

r=-\frac{\mathrm{1}}{{\kappa }_{\mathrm{a}}}\mathrm{l}\mathrm{n}(\mathrm{1}-R_{\mathrm{r}}).(15)

式中：r为能量束发射距离；R_r为区间（0，1）内均匀生成的发射随机数.

最终，根据随机发射方向和距离确定是否到达对应的节点微元.当能量束到达节点微元时，若考虑散射介质，根据吸收随机化的原理，需要进行随机取样来判断能量束是否被吸收.如果随机数R_ij ≥ε，则其被吸收；否则被散射.当能量束发射后并未到达节点而到达表面边界时，需要根据表面发射率判断能量束是被吸收还是被反射.类似地，如果随机数R_ij ≤ε，则其被吸收；否则被反射.

针对本文所考虑的二维物理模型，利用辐射全交换面积可得到辐射热源项：

Q_{i,j}=\sigma {\sum }_{k=\mathrm{0}}^{N_X+\mathrm{1}}{\sum }_{m=\mathrm{0}}^{N_Y+\mathrm{1}}{K}^{*}{\left(k,m,i,j\right)}_M\times \left(A_{M,T_{k,m}}{T}_{k,m}^{\mathrm{4}}-A_{M,T_{i,j}}{T}_{i,j}^{\mathrm{4}}\right)\left(\mathrm{1}\right.\le i\le N_X,\mathrm{1}\le j\le N_Y).(16)

式中：Q_{i， j} 为辐射热源，W；K^*（k，m，i，j） _M 是其他控制体（k，m）发射，控制体M（i，j）吸收的辐射传递因子；σ是斯蒂芬-玻尔兹曼常数，W/（m²·K⁴）；A_{M,T_{k,m}}是节点（k，m）的温度T_k，m 对控制体M的交换面积，m².

综上所述，对于各向异性导热-辐射耦合传热计算，首先质点在所划分的温度节点上进行移动，移动过程中会生成一个随机数，而后根据移动概率系数和随机数确定移动方向.每一个离散的质点在移动过程中都会携带一份离散的辐射热源，在平板内进行随机移动，直至到达边界.图3为耦合质点随机移动流程图.

2 结果与讨论

2.1 计算模型验证

由于目前很少有文献研究各向异性介质的导热-辐射耦合传热特性，因此本文针对各向同性介质的导热-辐射耦合传热进行了研究，将计算结果与离散传递法^［7］（DTM）以及格子玻尔兹曼法^［12］（LBM）的结果进行对比，以验证计算模型的正确性.

文献中所用物理模型为二维矩形，其中L_x =L_y =1 m，矩形各个壁面均为黑体壁面，内部介质为吸收、发射的各向同性介质，即k_{\mathrm{12}}=k_{\mathrm{21}}，k_{\mathrm{11}}=k_{\mathrm{22}}.平板N，W，E面壁面温度设定为T₀，平板S面壁面温度设定为T_R=2T₀.考虑折射率n=1，衰减系数β=1 m^-1，壁面发射率ε=1.0.表2展示了不同导热-辐射参数N_{\mathrm{c}\mathrm{r}}=k\beta /(\mathrm{4}\sigma T_{\mathrm{R}}^{\mathrm{3}})下，中轴线特定位置（x/L_x =0.5，y/L_y ）的稳态无量纲温度分布T/T_R.图4给出了导热-辐射参数N_cr=0.01，0.1，1.0时本文方法与DTM以及LBM在x/L_x =0.5位置沿y方向无量纲温度分布的比较.

由表2和图4可以看出，本文所提出方法的计算结果与两种文献方法得到的结果均吻合较好：与DTM结果相比的平均相对偏差为0.79%，最大相对偏差为3.41%；与LBM结果相比的平均相对偏差为0.18%，最大相对偏差为2.00%.这说明本文方法对于求解二维导热-辐射耦合传热问题具有很高的准确性.

2.2 各向异性热传导对耦合传热特性的影响

考虑相同的物理模型，忽略其他外在因素对传热过程的影响，即二维平板边界长度L_x =L_y =1 m；W，N，S，E面均为黑体壁面；设定S面温度为1 000 K，W，N，E面温度均为500 K.定义无量纲温度为T/T_R，本文以T_R=1 000 K为参考温度.设定导热-辐射参数N_cr=0.1，其他参数设置与2.1节保持一致.

材料的各向异性主要体现在各个方向的导热系数不同，并且在相互垂直的方向上体现得更为明显.本文以k_{\mathrm{11}}为参考导热系数，即导热-辐射参数N_{\mathrm{c}\mathrm{r}}=k_{\mathrm{11}}\beta /(\mathrm{4}\sigma T_{\mathrm{R}}^{\mathrm{3}})，定义k_{\mathrm{12}}^{*}=k_{\mathrm{12}}/k_{\mathrm{11}}为相对导热系数.根据互易性关系可知k_{\mathrm{12}}=k_{\mathrm{21}}，因此本文主要针对导热系数k_{\mathrm{12}}^{*}和k_{\mathrm{22}}^{*}，研究各向异性对整个传热过程的影响.

首先，本文研究导热系数的变化对垂直方向温度梯度的影响，考虑在k_{\mathrm{12}}^{*}=\mathrm{0}时，k_{\mathrm{22}}^{*}=0.5，0.75，1.0和1.5对平板传热过程的温度变化情况.计算所得在x/L_x =0.5位置沿y方向的无量纲温度分布曲线如图5所示.由图可知，随着k_{\mathrm{22}}^{*}的增加，温度曲线越来越平缓.当k_{\mathrm{12}}^{*}保持不变时，改变k_{\mathrm{22}}^{*}会影响温度在N，S方向的传导能力，导热系数的增加使得两个方向的传导能力随之增强.当时，辐射在导热-辐射耦合传热中的比重增加，温度的非线性在不断增强.当k_{\mathrm{22}}^{*}=\mathrm{0.5}时，平板高温区的下降趋势更加明显，并且平板的低温区更加靠近冷边界.当k_{\mathrm{22}}^{*}>\mathrm{0.75}时，随着k_{\mathrm{22}}^{*}的增加，曲线的下降趋势基本保持一致，说明热传导过程在导热-辐射耦合传热中的比重增大.

各向异性导热系数的变化会导致热流传递方向发生改变，影响整个传热过程的热传导能力，使得导热与辐射在传热过程中所占比重发生变化，进而对导热-辐射耦合传热产生影响.根据式（6）~式（9）可知，k_{\mathrm{12}}^{*}和k_{\mathrm{22}}^{*}的取值会影响质点在平板内随机移动的概率.因此，在考虑式（8）中各向异性导热系数设定的基础上，主要研究两种类型的情况：一是当k_{\mathrm{12}}^{*}与k_{\mathrm{22}}^{*}比值相同时，各向异性对平板温度场的影响；二是当k_{\mathrm{22}}^{*}保持恒定时，k_{\mathrm{12}}^{*}的变化对平板温度场的影响.综上所述，本文考虑了如表3所示的4种条件进行数值计算.图6为这4种条件下平板的无量纲温度分布图.

对于图6a和图6c，当k_{\mathrm{12}}^{*}与k_{\mathrm{22}}^{*}比值相同时，随着导热系数整体的增加，平板内温度梯度变化增大，但两种工况下热流传递方向趋势基本相同.图7为平板y/L_y =0.5位置沿x方向变化的无量纲温度曲线.从图7a中可以更加直观地看出，两条曲线温度差异最大出现在x/L_x =0.8的位置，最大相对偏差可达7.01%.并且工况3的温度曲线整体都要高于工况1的温度曲线.这说明当各向异性导热参数保持等量变化时，导热系数的增加会影响整体传热能力，而对传导方向的变化影响较小.同时辐射传热在较高温度区间的影响较为明显，从而导致温差随着温度的升高而加大.

对于图6b和图6d，两种工况下的温度分布基本相似，但由于各向异性增强，导热方向发生了变化.从图7b可以看出，当k_{\mathrm{12}}^{*}=\mathrm{0.25}时，整体无量纲温度曲线比k_{\mathrm{12}}^{*}=\mathrm{0.5}0时的曲线更加平缓对称.这主要是因为当k_{\mathrm{22}}^{*}保持不变时，随着k_{\mathrm{12}}^{*}的增加，平板内节点向其他节点的温度传导份额发生了不同的变化，导致平板内温度均匀性降低.同时，导热系数的增加，对应平板传导热量的能力随之增强.从图中可以看出，两条曲线在x/L_x =0.4处出现交点，并且k_{\mathrm{12}}^{*}=\mathrm{0.5}0对应的无量纲温度曲线在最高点时（x/L_x =0.8）略高于k_{\mathrm{12}}^{*}=\mathrm{0.25}的情况，最大相对偏差可达1.65%.

3 结 论

1） 通过将采用DTM与LBM求解导热-辐射耦合传热得到的计算结果与本文方法进行对比分析发现：相较于DTM，最大相对偏差为3.41%，平均相对偏差为0.79%；相较于LBM，最大相对偏差为2.00%，平均相对偏差为0.18%.结果表明：本文方法对于求解导热-辐射耦合传热具有较高的精度.

2） 二维平板传热过程中，k_{\mathrm{12}}^{*}的变化会导致热量传导能力向非垂直方向发生偏移，k_{\mathrm{22}}^{*}的改变会影响平板在N、S方向上的热传导能力，并且过小的k_{\mathrm{22}}^{*}会导致辐射在耦合传热过程中所占比重不断提升.因此，在选择各向异性材料时，根据所指定的温度条件合理地选择材料，可以优化导热和辐射的耦合传热过程.

3） 当k_{\mathrm{12}}^{*}与k_{\mathrm{22}}^{*}比值相同时，平板的各向异性基本保持不变.随着导热系数的增加，热传导能力会沿相同的方向增强.辐射的参与使得温差在较高温度区间达到最大，最大相对偏差可达7.01%.因此，尽管不同材料的各向异性可能存在相似性，但导热系数的大小对耦合传热过程具有直接影响.过小的导热系数会导致导热和辐射的比重发生变化.

4） 当k_{\mathrm{22}}^{*}保持不变时，k_{\mathrm{12}}^{*}的增加会导致平板内温度梯度发生变化，从而导致平板内温度的均匀性降低，但同时也会提高平板的热传导能力.在高温环境下，由于导热系数的改变，可能会导致温度出现一定的偏差.在本文所设定的计算条件下，最大相对偏差为1.65%.因此，在选择具有方向性传热特性的材料时，需要考虑各向异性对热传导能力的影响.

参考文献

[1] Ge W A， Zhao C Y， Wang B X. Thermal radiation and conduction in functionally graded thermal barrier coatings. Part I： experimental study on radiative properties［J］. International Journal of Heat and Mass Transfer， 2019， 134： 101-113.

[2] Ge W A， Zhao C Y， Wang B X. Thermal radiation and conduction in functionally graded thermal barrier coatings. Part Ⅱ： experimental thermal conductivities and heat transfer modeling［J］. International Journal of Heat and Mass Transfer， 2019， 134： 166-174.

[3] Wellele O， Orlande H R B， Ruperti J N， et al. Coupled conduction-radiation in semi-transparent materials at high temperatures［J］. Journal of Physics and Chemistry of Solids， 2006， 67 （9/10）：2230-2240.

[4] Yi Z， Su Z G， Li G J， et al. Development of a double model slab tracking control system for the continuous reheating furnace［J］. International Journal of Heat and Mass Transfer， 2017， 113： 861-874.

[5] Liu D L， Chen H L. Numerical verification of a nonlocal discrete model for anisotropic heat conduction problems［J］. International Journal of Thermal Sciences， 2023， 191： 108360.

[6] Norouzi M， Rahmani H. An exact analysis for transient anisotropic heat conduction in truncated composite conical shells［J］. Applied Thermal Engineering， 2017， 124：422-431.

[7] Mishra S C， Talukdar P， Trimis D， et al. Computational efficiency improvements of the radiative transfer problems with or without conduction： a comparison of the collapsed dimension method and the discrete transfer method［J］. International Journal of Heat and Mass Transfer， 2003， 46（16）： 3083-3095.

[8] Mishra S C， Roy H K. Solving transient conduction and radiation heat transfer problems using the lattice Boltzmann method and the finite volume method［J］. Journal of Computational Physics， 2007， 223（1）：89-107.

[9] Ma J， Sun Y S， Li B W， et al. Spectral collocation method for radiative-conductive porous fin with temperature dependent properties［J］. Energy Conversion and Management， 2016， 111： 279-288.

[10] Ma J， Sun Y S， Li B W. Spectral collocation method for transient thermal analysis of coupled conductive， convective and radiative heat transfer in the moving plate with temperature dependent properties and heat generation［J］. International Journal of Heat and Mass Transfer， 2017， 114： 469-482.

[11] Sun Y J， Zhang X B， Howell J R. Non-gray combined conduction and radiation heat transfer by using FVM and SLW［J］. Journal of Quantitative Spectroscopy and Radiative Transfer， 2017， 197： 51-59.

[12] Wei Y J， Liu X C， Zhu K Y， et al. A unified lattice Boltzmann framework for combined radiation-conduction heat transfer［J］. International Journal of Heat and Mass Transfer， 2023， 200： 123513.

[13] Sun Y S， Li X Y， Zhao J Z， et al. Investigation of transient coupled conduction and radiation heat transfer in the linearly anisotropic scattering cylindrical medium by spectral collocation method［J］. International Journal of Thermal Sciences， 2022， 172： 107308.

[14] Modest M F. Radiative heat transfer ［M］. Cambridge： Academic Press， 2003：644-679.

[15] 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.

[16] Li G J， Zhong J Q， Wang X D. An improved Monte Carlo method for radiative heat transfer in participating media［J］. Journal of Quantitative Spectroscopy and Radiative Transfer， 2020， 251： 107081.

[17] Li D Y， Li G J， Hong D L， et al. Improved Monte Carlo method for radiative heat transfer in semitransparent media with BRDF surface［J］. International Journal of Thermal Sciences， 2023， 187： 108152.

[18] Wong B T， Pinar-Mengüç M. A unified Monte Carlo treatment of the transport of electromagnetic energy， electrons， and phonons in absorbing and scattering media［J］. Journal of Quantitative Spectroscopy and Radiative Transfer， 2010， 111（3）： 399-419.

[19] Haji-Sheikh A， Sparrow E M. The solution of heat conduction problems by probability methods［J］. Journal of Heat Transfer， 1967， 89（2）： 121-130.

[20] Kowsary F， Arabi M.Monte Carlo solution of anisotropic heat conduction［J］. International Communications in Heat and Mass Transfer， 1999， 26（8）：1163-1173.

[21] Kowsary F， Irano S. Monte Carlo solution of transient heat conduction in anisotropic media［J］. Journal of Thermophysics and Heat Transfer， 2006， 20（2）： 342-346.

[22] Naeimi H， Kowsary F. Finite volume Monte Carlo （FVMC） method for the analysis of conduction heat transfer［J］. Journal of the Brazilian Society of Mechanical Sciences and Engineering， 2019， 41（6）： 260.

[23] 谈和平， 夏新林， 刘林华， 等. 红外辐射特性与传输的数值计算： 计算热辐射学［M］. 哈尔滨： 哈尔滨工业大学出版社， 2006： 157-173.

[24] Tan He-ping， Xia Xin-lin， Liu Lin-hua， et al. Numerical calculation of infrared radiation characteristics and transmission： computational thermal radiation［M］. Harbin： Harbin Institute of Technology Press， 2006： 157-173.

基金资助

辽宁省自然科学基金资助项目(2023-MSBA-059)