Bidirectional monte carlo method for thermal radiation transfer in participating medium

Abstract: A bidirectional Monte Carlo (BDMC) method based on reversibility of bundle trajectory and reciprocity of thermal radiative energy exchange was developed to solve radiative heat transfer in absorbing and scattering medium. Two types of sampling models were introduced into the Monte Carlo (MC) simulation, namely the equivalent sampling and the weight sampling, respectively. Mathematical formula for the sampling models and the statistical calculation of sampling bundles were derived. Furthermore, the reciprocity error correlation of radiative exchange factors between the BDMC method and the traditional Monte Carlo (TMC) method were demonstrated and analyzed. Radiative heat transfer in a two-dimensional rectangular domain with absorbing and scattering media was solved by using both the BDMC method and the TMC method. Radiative exchange factors and radiative equilibrium temperature profiles predicted by the BDMC method were compared with those predicted by the TMC method. The performance parameter P, defined to evaluate the performance of MC methods, was computed and compared between the BDMC and the TMC methods. The results showed the superiority of BDMC method compared with the TMC method for radiative heat transfer, in addition, the weight sampling was proved to be more flexible than the equivalent sampling in the BDMC method.


Introduction
Radiative heat transfer in participating medium is described by radiative transfer equation (RTE), which is an integro-differential equation. Many methods have been developed to solve the RTE, the common used numerical methods such as the discrete ordinate method (DOM) and the finite volume method (FVM), rely on both spatial discretization and angular discretization, which are difficult to solve the RTE with high accuracy. Monte Carlo (MC) method is a stochastic statistical method based on the physical processes, it has been applied to solve radiative heat transfer in various participating media [1][2][3][4][5][6]. Moreover, the results predicted by MC method can often be treated as benchmark solutions due to its high accuracy of solution [7,8].
However, for the statistical nature of the MC method, the high computational cost is still a considerable disadvantage of MC simulation. Many attempts have been made in order to improve the computational efficiency of the method [9]. For example, several sampling approaches were developed to improve the speed of convergence, such as the importance sampling method [10], the rejection sampling method [11], the differential sampling method [12] and the weight-equivalent sampling method [13]. In addition, parallel computing technique was introduced into the MC simulation to improve computation efficiency evidently [1,9]. Howell [9] analyzed the advantages of various programming strategies of the MC method for radiative heat transfer in absorbing and scattering medium. All of these improvements were mainly mathematical efforts and hardware improvement of computer, which are universal to the MC simulations.
It has been noticed that the physical feature of thermal radiation transfer can also provide possibility to improve the MC simulation of radiative heat transfer in participating medium. For example, Walters [14] developed a reverse method based on the reciprocity principle for radiative heat transfer in a generalized enclosure containing an absorbing, emitting and scattering medium, the reverse method was proved to be efficient. Cherkaoui et al. [15,16] developed a net exchange Monte Carlo method based on a net-exchange formulation, provided an efficient way of systematically fulfilling the reciprocity principle, the computing time was proved much smaller than the conventional Monte Carlo approach. Lataillade [17] and his cooperators applied the net exchange Monte Carlo approach for radiative heat transfer in optically thick medium with spectral dependent radiative properties. Eymet et al. [18] extended this method to absorbing, emitting, and scattering media. Tessé et al. [19] improved the forward Monte Carlo (FMC) method based on the reciprocity principle, the method was used for radiative transfer in real gases, and proved to be a better choice for optically thick or nearly isothermal media compared with the forward Monte Carlo method. In addition, the reverse Monte Carlo (RMC) methods, based on the reversibility of radiative transfer trajectory, has been developed to solve the radiative transfer in absorbing and scattering media [20][21][22][23][24]. Kovtanyuk et al. [25] presented a recursive algorithm based on modification of Monte Carlo method, the modified method was used to solve the coupled radiative and conductive heat transfer in an absorbing and scattering medium, and was proved to be more accurate. Soucasse et al. [26] proposed a Monte Carlo formulation for radiative transfer in quasi-isothermal media which consists in directly computing the difference between the actual radiative field and the equilibrium radiative field at the minimum temperature in the medium.
In the present study, a bidirectional Monte Carlo (BDMC) method based on reversibility of bundle trajectory and reciprocity of radiative energy exchange was developed to solve thermal radiation transfer in absorbing and scattering medium. Two types of sampling models for MC simulation were presented, namely the equivalent sampling and the weight sampling. The equivalent sampling was chosen for the uniform mesh while the weight sampling was more suitable to the non-uniform mesh. The bidirectional information of tracing a sampling bundle was utilized by the BDMC method, the solution precision or efficiency can be evidently improved. Radiative heat transfer in a two-dimensional rectangular domain with absorbing and scattering media was solved by the BDMC method and the TMC method, respectively. The radiative exchange factors and the temperature profiles were investigated, in addition, the performance parameter defined by Howell [9] was also calculated to evaluate the two MC methods.

Description and formulation
Radiative heat transfer in a two-dimensional (2-D) rectangular domain with absorbing, emitting and/or scattering medium was investigated. Figure 1 shows the 2-D rectangular geometry as well as the coordinate system. The four walls were assumed to be diffuse and gray. Radiative properties, such as the absorption coefficient, the scattering coefficient were assumed to be constant. The rectangular domain was divided into M × N = MN cells, any of which was depicted by V ij (medium cell ) or A ij (wall cell), wherein the subscripts I [1,M] and j [1,N]. Monte Carlo method has robust

The TMC method and radiative exchange factor
In the present study, the radiative exchange factor , ij kl RD   was introduced to decouple the solution of radiative transfer from that of the temperature profile. It was defined as a fraction of the spectral radiative energy emitted from cell V ij (or A ij ), that is absorbed by cell V kl (or A kl ) [27,28].
is tenable according to the definition of radiative exchange factors. Note that the radiative exchange factor depends only on the system geometry and the radiative properties distribution of the medium [29]. The reciprocity relation between a couple of radiative exchange factors was given by [30].
As Eq 6, after solving the radiative exchange factors, the net radiative energy exchange ij kl As Eq 7, where, is the Stefan-Boltzmann constant (= 5.67×10 -8 W/m 2 K 4 ), E bλ (T) is the spectral blackbody emissive power at temperature T.
In the MC simulation, As Eq 8, the number ij kl N  can be counted from the MC simulation, however, only an approximate value of ij kl N  can be obtained because of the pseudo randomicity [32] in the sampling process. In fact, the reciprocity correlation shown in Eq 6 cannot be strictly satisfied, and random errors always exist. In order to obtain more accurate results, one need to increase the sampling bundles ij N , which results in the increasing computing cost. In the TMC method, the propagation trajectories of ij N sampling bundles for cell ij V are tracked and counted to get the value of ij kl N  for any cell kl V . Similarly, the number kl ij N  is obtained after tracing another kl N propagation trajectories of the sampling bundles for cell kl V . Wherein, only the forward propagation information of a trajectory is used, the number of sampling bundles for cell ij V is ij N . According to the reversibility of light propagation, for a tracing trajectory from kl V to ij V , a bundle from ij V can also propagate to kl V along the reverse direction of the tracing trajectory. The information of the reverse direction can be used in the MC simulation, the bidirectional information can be obtained from the forward tracing, which results in the development of the present BDMC method.

The BDMC method and the equivalent sampling
In order to simplify the description of the BDMC method, the formula in the following sections were given for gray medium with gray walls. The basic idea of BDMC is to calculate the radiative exchange factor by the reversibility of light beams. According to the reversibility of light beams, the transmission path of the bundles received by cells can also be used as the transmission path of the bundles transmitted by cells. In this way, the effective information of sampling bundles if fully utilized without increasing the number of sampling bundles and calculation amount. Theoretically, the statistical sampling bundles in calculation can be doubled. And it can better satisfy the relationship of the reciprocal of radiative exchange factor. Figure 1 shows the schematic diagram of radiative transfer in a rectangular domain with participating medium in the BDMC method, it can be noticed that there are ij N tracing trajectories departing from cell ij V , and where, ij kl RD   is used to denote the radiative exchange factor in the BDMC method. The information of the bidirectional trajectories has been used, as is shown in Eq 11, it indicates that the effective sampling bundles are increased, while the computing cost remains unchanged. Therefore, more accurate results can be predicted from Eq 11 than that from Eq 8 under the same computing cost.
It should be pointed out that Eqs 9 and 11 are valid only if the sampling bundles for cell V ij and any other cells are equivalent. According to the reversibility of light beams, the transmission path of the bundles received by cells can also be used as the transmission path of the bundles transmitted by cells. In this way, the effective information of sampling bundles if fully utilized without increasing the number of sampling bundles and calculation amount. In other words, the forward and reverse bundles in the BDMC method should have the same contribution for radiative exchange factors calculation. To use Eq 11, the equivalent sampling was introduced, in which the number of sampling bundles for cell (i, j) was determined by where 0 N is the sampling density for a reference case defined as where, r N is the number of sample bundles of the reference cell r V , the subscript r refers to the reference cell.

Error analysis and weight sampling
The bidirectional counting of sampling bundles employed in the BDMC method is equivalent to double the sampling number if the random error does not exist. In fact, any of the Monte Carlo methods always accompanied by random error. For the BDMC method, that leads to 2 ij ij N N   and the energy of a bundle is not strictly equal to other cells. For any cell, it can be easily demonstrated by Eq 11, and the radiative exchange factors calculated by the BDMC method satisfied the conservation relation strictly. But the reciprocity relation between a couple of radiative exchange factors cannot be satisfied strictly. According to Eq 11, the following equation can be obtained , , According to Eq 8, for gray medium, (15) replace N ij and N kl according to Eq 12, in addition, considering Eqs 4 and 5, then Eq 15 can be transformed into , , where, Considering the conservation relation therefore, , , Because of the randomicity of the reciprocity error ij kl   , the value of The reciprocity error for a couple of radiative exchange factors predicted by the BDMC method is always smaller than those predicted by the TMC method.
For some cases, the equivalent sampling may be inconvenient or inaccurate. For example, if the differences of apparent radiation characteristics and/or cell volume are very large among different cells, the equivalent sampling will result in a very small sampling number for a cell, and may be a very large sampling number for another. The former is not expected for the statistical analysis, while the latter increase computational cost.
The weight sampling was introduced into the BDMC method to avoid the above shortages of the equivalent sampling. For weight sampling, the number of sampling bundles for a cell was determined only based on the statistical request, but the weight of a sampling bundle was taken into account in the final statistical calculation. If the number of sampling bundles for cell V ij was N ij , then, its bundle weight W ij was defined as where N r is the number of sample bundles for reference cell V r . Then, the radiative exchange factor Eq 11 should be substituted by Eq 24 if the weight sampling is employed.

Numerical analysis on the BDMC method
Radiative transfer in a rectangular domain with absorbing, emitting and isotropic scattering gray medium was solved separately by the BDMC method and the TMC method. The radiative exchange factors, the radiative equilibrium temperature profiles, and the performance parameter defined by Farmer and Howell [1,9], predicted by the BDMC method were compared with those predicted by the TMC method. In addition, the performance of the equivalent sampling and the weight sampling in the MC simulation were also examined. The radiative properties of the medium considered were supposed to be uniform and depicted by constant extinction coefficient  e , scattering albedo  , and refractive index n = 1.0. The length and width of the rectangular domain were  Table 1, where, x c = y c were the width and length of the largest volume cell in the non-uniform mesh system, and x c = y c , the parameter K = 1.04, M = N = 21, the mesh sizes in the non-uniform mesh system were expressed as The sampling bundles needed in the MC simulation should be estimated as the results were calculated based on statistical computation. Figures 2 and 3 show the radiative equilibrium temperature profiles at dimensionless locations of Y = 0.3 and Y = 0.5 predicted by the TMC method with different sampling bundles of N b = 10 4 , 10 5 , and 10 6 , respectively. Figure 2 shows the results for absorbing medium, while the results shown in Figure 3 were for absorbing and isotropic scattering medium with scattering albedo of = 0.5. For both the two cases, the wall temperatures were imposed separately as T WE = T WW = 1000 K and T WN = T WS = 1500 K, the emissivity of the walls were assumed to be constant and equal to 0.5, the optical thickness = k e L = 15. For both the absorbing and absorbing-scattering medium, the predicted temperature profiles were found to be close to each other for different sampling bundles of N b = 10 4 , 10 5 , and 10 6 , moreover, the temperature difference for N b = 10 5 and N b = 10 6 were less than 0.5%, therefore, the results were considered to reach convergence for sampling bundles N b = 10 5 in the present study, and the sampling bundles of N b = 10 5 was also chosen as the reference sampling bundles N r in the BDMC method.

Examination of radiative exchange factors
In this section, the radiative exchange factors predicted by the BDMC method were compared with those predicted by the TMC method, the numerical simulations were implemented using both uniform and non-uniform grids, meanwhile, the equivalent sampling and the weight sampling were also employed separately. The predicted radiative exchange factors always satisfy the conservation relation strictly, therefore, only the reciprocity for radiative exchange factors was examined in the present study. The relative error of the reciprocity for a couple of radiative exchange factors was defined as 2 100% ij kl kl ij where, ij kl RD  refers not only to the modified form of the radiative exchange factors given by Eq 8 for the TMC method, but also to those given by Eq 11 (BDMC method, equivalent sampling) and by Eq 24 (BDMC method, weight sampling). Table 2 shows the modified radiative exchange factors predicted by the BDMC method and the TMC method employing the equivalent sampling and the uniform grids as well as the corresponding relative error D r defined by Eq 25. The emissivity of the walls were taken as = 0.5, the optical thicknesses of the rectangle medium along the length and width direction were 1.0   , the reference sampling bundles was taken as 5 10 r N  . It indicates that the reciprocity error for a couple of radiative exchange factors predicted by the BDMC method was at least an order of magnitude less than those predicted by the TMC method with the same sampling bundles. The maximum relative error for all the elements considered did not exceed 0.58% in the BDMC method, while it reached 30.1% in the TMC method.   Table 3. Comparison of radiative exchange factors predicted by the BDMC method with the equivalent sampling and the weight sampling using the non-uniform grids.  (11,11) (1,11) 6.730  10 -4 5.000  10 -4 29.5 9.413  10 -4 9.414  10 -4 0.01 (11,11) Table 3 shows the radiative exchange factors predicted by the BDMC method using the equivalent sampling and the weight sampling in a non-uniform mesh system. The computing parameters were 0.5   , , and 5 10 r N  . The reciprocity error for a couple of radiative exchange factors predicted by the BDMC method with the equivalent sampling was unacceptable for the non-uniform mesh, it indicates that the equivalent sampling was not recommended to solve radiative transfer in the BDMC method with non-uniform grids. If there is obvious difference of heat capacity between cells, the equivalent sampling will lead to big difference of sampling number for the same cells. The calculation accuracy will come under influence if the equivalent sampling is applied. However, the BDMC method with the weight sampling can predict quite satisfactory results for all the cells considered. Table 4 shows the radiative exchange factors predicted by the BDMC method and the TMC method employing the weight sampling and the non-uniform grids. It can be seen that the reciprocity error in the TMC method was larger than those in the BDMC method for all the cells considered. The BDMC method showed greatly superior to the TMC method when the weight sampling and the non-uniform grids was used. Table 4. Comparison of radiative exchange factors by the BDMC method and the TMC method using the weight sampling and the non-uniform grids.

Examination of radiative equilibrium temperature profiles
The advantage of the BDMC method was further verified by comparing the predicted temperature profiles. For gray medium, the radiative equilibrium temperature ij T of cell (i,j) satisfies the energy equation Utilizing the predictions of radiative exchange factors, the radiative equilibrium temperature profile can be solved iteratively from Eq 26 or Eq 27. In the iteration, the iteration is stopped by setting the residuals. The number of iterations is not certain, which is related to the residuals.   Figure 4 shows the radiative equilibrium temperature profiles for an absorbing, emitting, and non-scattering medium predicted by the TMC method and the BDMC method employing the equivalent sampling and the uniform grids. The wall temperatures were imposed separately as T WE = T WW =1000 K and T WN = T WS =1500 K, the emissivity of the walls were assumed to be constant and equal to 0.5, the optical thickness = k e L x = k e L y = 15. First, the BDMC predictions were very close to the TMC predictions, it indicates that the development of the present BDMC for predicting temperature profiles is correct, in addition, the BDMC predictions were shown smoother than those predicted by the TMC method, the BDMC method converged faster than the TMC method. Similar superior of the BDMC method to the TMC method can be seen from Figure 5, where the weight sampling and the non-uniform grids were employed and the medium was absorbing and isotropic scattering with scattering albedo of = 0.5. Figure 6 shows the comparison of temperature profiles predicted by the BDMC method with non-uniform grids employing different sampling models of the equivalent sampling and the weight sampling. The temperature profiles predicted by using the weight sampling was smooth, while the predictions using equivalent sampling model contains some small rectangles, this may due to the fact that the number of the sampling bundles for different cells employing equivalent sampling were significantly different, which lead to large random error and affected the temperature profiles. Figure 7 shows the radiative equilibrium temperature profiles predicted by the TMC method and the BDMC method with a relative small number of weight sampling bundles using the non-uniform grids. The temperature profiles predicted by the TMC method became unacceptable as the reference sampling bundles decreased to N r = 10 3 , while the BDMC method can still predict better results, the BDMC method converged faster than the TMC method. Therefore, the BDMC method would be a better choice if the sampling bundles for radiation computation is limited or efficient computation is required.

Performance parameter analysis
Farmer and Howell introduced the performance parameter P to evaluate various MC methods or strategies [1,9], the performance parameter was defined as where, t is the CPU time spent by the concerned MC simulation and 2  is the variance of the results. A good method or strategy for the MC simulation tends to minimize the performance parameter P. The variance 2  is given by It has been discussed that the acceptable numerical results can be obtained if the sampling bundles N b in the MC simulation reached 10 5 . In order to compute the performance parameters for the TMC method and the BDMC method, the "exact solution" of the radiative exchange factors should be firstly obtained. The approximate "exact solution" can be obtained by increasing the sampling bundles of the MC simulation due to the high precision of the MC method. To do so, the reference sampling bundles were taken as N r = 10 6 and N r = 10 8 , respectively, the predicted results for the two cases were found to be nearly unchanged, and the maximum relative error for any of the radiative exchange factors correspond to different cells was less than 0.1% as the reference sampling bundles N r increased from 10 6 to 10 8 . Therefore, the results employing N r = 10 8 were treated as the "exact solution" in the present study.
For each exercise, only the initial random number and sampling order were changed, therefore, the CPU time spent by each exercise should be nearly unchanged. The total exercise times was set as H =10. The performance parameter of the BDMC method and the TMC method were written as P BDM and P TM , respectively, while the performance increment P  was defined as  Table 5 reports the performance parameter for the TMC method and the BDMC method using the equivalent sampling and the uniform grids. It indicates that the performance parameter for the TMC method was always larger than that for the BDMC method, in addition, the performance increment decrease with the increasing sampling bundles. Table 6 shows the performance parameter for the TMC method and the BDMC method employing the weight sampling and the non-uniform grids. Similar conclusions can be drawn as those employing the equivalent sampling and the uniform grids. Table 7 shows the performance parameter for the BDMC method with equivalent sampling and the weight sampling using the non-uniform grids. The performance parameter for BDMC with the equivalent sampling was found to be larger than that with the weight sampling, this indicate that the weight sampling was more suitable for non-uniform grids in the BDMC simulations. Table 5. Performance comparison of the TMC method and the BDMC method using the equivalent sampling and the uniform grids.   Table 7. Performance comparison of the equivalent sampling and the weight sampling in the BDMC method using the non-uniform grids.

Conclusions
The BDMC method along with the equivalent sampling and the weight sampling were developed to solve radiative transfer in absorbing and scattering medium. The BDMC method approximately doubles the information from the bundle tracing to possess a high efficiency by making the best use of the reversibility of the bundle trajectory. The formula for the BDMC method and the corresponding error analysis were derived and presented. Radiative heat transfer in a two-dimensional rectangular domain of absorbing and/or scattering medium were solved by the TMC method and the BDMC method, respectively. The reciprocity of the radiative exchange factors, the radiative equilibrium temperature profiles, and the performance parameter predicted by the two MC methods were examined and compared. The results showed that (1) the BDMC method can greatly improve the reciprocity satisfaction of radiative exchange factors, which is helpful for temperature profile solution, (2) the BDMC method was more accurate or efficient than the TMC method, (3) in the BDMC simulations, the weight sampling was found to be more flexible than the equivalent sampling.