Analysis of the multi-spectral energy bundle method: finding the number of sub-bundles for the minimum computational time

The multi-spectral energy bundle (MSB) is a method proposed to reduce the computational time of Monte Carlo simulations in media with spectrally dependent properties. In this method the typical energy bundles of the Monte Carlo are treated as multi-spectral by assuming that they are composed of a set of monochromatic sub-bundles. The number of sub-bundles in each bundle affects the computational time in such a way that there exists an optimum number of sub-bundles that leads to the minimum time for the computation. However, due to the stochastic nature of the Monte Carlo, the optimum number of sub-bundles is not easy to be exactly found. This paper proposes a methodology which allows to accurately determine the number of sub-bundles for the minimum computational time, and apply the proposed methodology to 25 cases consisting of a one-dimensional slab of gas composed of H2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document}O and non-emitting/absorbing species, where the FSK was used to consider spectral dependence of properties. The results showed that for the majority of the considered cases, 3 sub-bundles lead to either the minimum or a computational time quite close to the minimum, and always provides reduction of computational time. Thus, 3 sub-bundles is indicated for the analyzed class of problems.


3
(FSK), presented by Modest and Zhang [7,8] is a transformation of variables, which allows one to take into account at once several spectral intervals in which the absorption coefficient assumes specific same values. The method is exact for homogeneous gases, that is, it provides results as accurate as the line-by-line. However, when the state of the gas varies with location, the correlated-k assumption must to be invoked, that leads to some deviation of the FSK results with respect to the line-by-line solution. The absorption-line-blackbody distribution function (ALBDF) [9], employed in the spectral-line-based-weighted-sum-of-gray gases model (SLW) [10] is closely related to the cumulative k-distribution function of the FSK. In the case in which the number of gray gases approaches infinity, the method is called exact SLW [11] and provides solutions with line-byline accuracy for homogeneous media, likewise the FSK.
The Monte Carlo is a powerful method which can deal with some difficulties in a relatively easy manner as compared to other numerical techniques [2,3]. The method is well suited to consider irregular geometries, directional and spectral surface properties and anisotropic scattering of radiation. However, the Monte Carlo is associated with high computational cost, especially when spectral properties of molecular gases are considered in detail. The Monte Carlo can be applied to line-by-line computations [12] or with spectral models. Wang et al. [13] combined the Monte Carlo with the FSK; Maurente et al. [14,15], applied the Monte Carlo along with the ALBDF.
The multi-spectral energy bundle method (MSB) improves the computational efficiency of the Monte Carlo method applied to media with spectrally dependent properties [12,16,17]. In this method the Monte Carlo energy bundles are treated as multi-spectral and composed of a set of monochromatic sub-bundles. The reduction on computational time depends on the number of sub-bundles, and as demonstrated by Maurente and França [16] and Maurente and Freitas [17], there exists an optimum number of subbundles, which provides the minimum computational time. Thus, the efficiency of the MSB depends on this optimum number of sub-bundles. However, due to statistical oscillations inherent to the Monte Carlo, it is difficult to exactly determine the number of sub-bundles which will be associated with the minimum computational time.
In this research, a methodology which allows one to find with accuracy the optimum number of sub-bundles is proposed and applied to 25 test cases, all consisting of a gas slab confined between two cold walls. For all cases, the gas was composed of H 2 O and non-participating species. However, the distance between walls and/or H 2 O concentration was different from case to case. The employed spectral model was the FSK.
The next section explains some concepts about the MSB method and its application along with the FSK model; following, in Sect. 3, the proposed methodology to determine the optimum number of sub-bundles is presented.

The multi-spectral energy bundle formulation
The formulation of the multi-spectral energy bundle (MSB) used in this research was presented by Maurente and França [16]. The following Sect. 2.1 presents equations for emission of multi-spectral energy bundles; Sect. 2.2 treats the absorption process, while the last Sect. 2.3 shows how to apply the MSB along with the FSK spectral model.

Emission of multi-spectral energy bundles
In the traditional Monte Carlo, in which the bundles are monochromatic, the energy carried by a bundle emitted from a uniform volume of medium is where N is the total number of bundles emitted per unity time and Q e is the rate of emission of energy from the volume, which can be computed as where I b,η is the blackbody spectral intensity given by the Planck distribution, κ η is the spectrally dependent absorption coefficient, η is the wavenumber and is the solid angle. Thus, the fraction of energy emitted in the wavenumber η around dη is Spectrally dependent absorption coefficient of a medium constituted of 20% H 2 O and N 2 , at 1000 K and 1 atm, computed using spectroscopic data from HITRAN2012 [1] where the numerator is the amount of energy emitted in the wavenumber η, per unity time, and the denominator is the rate of emission of energy, that is, the total energy emitted per unit time.
Assuming that the energy is emitted in N energy bundles per unity time, all carrying the same amount of energy (Q e /N), the fraction P(η)dη may be thought as the probability of emission of an energy bundle with wavenumber η. Therefore, from the probability and statistical theory, P(η) is the probability density function. Consequently, the cumulative distribution function of probabilities is where the asterisk denotes a dummy variable of integration.
The cumulative distribution function given by Eq. (4) is equivalent to the probability that any given energy bundle has a wavenumber between −∞ and η. However, since a negative wavenumber has no physical meaning, Eq. (4) is equivalent to the probability that the bundle has a wavenumber between 0 and η. Thereafter, the probability that a bundle has a wavenumber between 0 and +∞ (when η → ∞) is equal to 1.
Cumulative distribution functions are used in Monte Carlo simulations of radiation heat transfer to determine quantities associated with energy bundles, as direction of emission, wavenumber and location of absorption [2,3]. In the case of the example, the wavenumber of an energy bundle can be determined by first choosing a random value between 0 and 1 to R, and then obtaining η(R) after inversion of Eq. (4).
In traditional Monte Carlo simulations, Eq. (4) is used to determine the wavenumber of the monochromatic bundles. In the case of the MSB, Eq. (4) still applies; however, it is used to compute the wavenumber of the sub-bundles.
The energy of the multi-spectral bundle will be the summation of the energy of its sub-bundles, that is: where b η,o,j is the energy of the j-th sub-bundle and N η is the total number of sub-bundles that compose the multispectral bundle. Thus, in order that the energy carried by the multi-spectral bundle be equal to that of a monochromatic bundle, the energy of each sub-bundle must be given by the following equation:

Absorption process for MSB
There are two methods of simulating the absorption process: one is assuming that a bundle is entirely absorbed in a certain point of the medium at random; the other is considering that the energy is gradually absorbed as the bundle travels along a certain path in the medium. According to [13], the second method is the most computationally efficient. For this reason, only the formulation for continuous absorption will be presented here.
In traditional Monte Carlo simulations the fraction of the energy of a monochromatic bundle that is absorbed as the bundle travels a path s in a uniform absorbing medium is given by This equation still applies in the case of the MSB method; however, it is used to compute the fraction of the energy of each sub-bundle that is absorbed as the multi-spectral bundle travels the path s. It follows that the energy absorbed in a volume of medium as a multi-spectral bundle travels a path s within this volume is where p η,j is given by Eq. (7), and the subscript j is the index of the sub-bundle in the multi-spectral bundle. The subscript o indicates that b η,o,j is evaluated before the bundle travels the path s. Consequently, the energy remaining in each sub-bundle after the bundle travels the distance s is The energy remaining in the bundle is given by the summation of the energy remaining in each sub-bundle: As the bundle enters the next volume of medium, the energy carried by each sub-bundle is updated to b η,o,j = b η,j , and the procedure described above is repeated. The energy of the multi-spectral bundle is gradually depleted along its path until either the bundle carries an insignificant amount of energy or is entirely absorbed by a surface boundary. The MSB can be applied with different spectral models. The next section presents the MSB formulation for applications along with the FSK.

The multi-spectral energy bundle with the FSK method
In the present study, only the basic concepts and equations are going to be presented for a uniform volume of medium. Modest [3] presents a comprehensive explanation about the FSK method, and details about the application of the Monte Carlo along with the FSK or ALBDF can be found in [13][14][15].
The full spectrum k-distribution is given by where δ(k − κ η ) is the Dirac-delta function and I b is the total blackbody intensity, given by The cumulative k-distribution is computed from the fullspectrum k-distribution as where the asterisk indicates dummy variable of integration.
Physically, the cumulative k-distribution, g(k), gives the fraction of the blackbody energy in the spectral intervals where the absorption coefficient, κ η , is lower than k. Thus, it is clear that as k varies from zero to the maximum value of κ η , g(k) varies from zero to unity.
In FSK Monte Carlo computations, Eq. (4) can no longer be used for randomly determining the wavenumber of the bundle and its respective absorption coefficient, since this equation is suited only for line-by-line computations. Maurente et al. [14] proposed a cumulative distribution function of probabilities which is directly related to the cumulative k-distribution (or to the ALBDF) as In order to determine the absorption coefficient of an energy bundle in FSK Monte Carlo simulations, a number is randomly chosen for R g in an interval between 0 and 1. Then, g(k) is set equal to this number and subsequently a corresponding k(g) is found for the absorption coefficient. In such a formulation, the amount of energy in the bundles varies with the absorption coefficient. However, Eqs. (1), (5) and (6), for emission of multi-spectral energy bundles, are still valid, provided that Q e be replaced with (4πkVI b ) , in Eqs. (1) and (6). Similarly, for simulating the absorption of multi-spectral energy bundles, Eqs. (7)-(10) can be applied, but using the absorption coefficient k of the FSK instead of the spectrally dependent absorption coefficient, k η . Therefore, κ η should be replaced with k in Eq. (7).
The FSK cumulative k-distribution functions, g(k), are calculated from the line-by-line absorption coefficients obtained from databases such as HITRAN and HITEMP [1,18]. Thus it is necessary to precalculate them before performing radiation heat transfer computations, otherwise the FSK would become even more computationally expensive than line-by-line calculations. Pearson et al. [19] presented useful correlations, computed from data obtained with HITEMP2010 [18], which provides the g(k) functions for performing FSK computations in media constituted by H 2 O and non-participating species. These correlations were employed in this work to obtain the necessary g(k).

Methodology to determine the optimum number of sub-bundles
Monte Carlo results are usually presented within a confidence interval, which decreases as the number of bundles is increased. As described by Maurente and França [16] and Maurente and Freitas [17], when multi-spectral bundles are used, the necessary number of bundles to find the results within a certain prescribed level of uncertainty tends to decrease as the number of sub-bundles is increased. On the other hand, the computational time per bundle increases linearly with the number of sub-bundles. where obtained by Maurente and França [16], for line-byline MSB computations. Observing Fig. 3, one can notice that the optimum number of sub-bundles is in between 3 and 6, but cannot be exactly identified. This is due to statistical oscillations inherent to Monte Carlo results. In order to identify the optimum number of sub-bundles the statistical oscillations must be reduced. However, following the procedures adopted by Maurente and França [16] and Maurente and Freitas [17], it would require extremely high computational time.
In this paper a methodology to accurately determine the optimum number of sub-bundles is proposed. It was developed based on one-dimensional problems consisting of gas slabs confined between two parallel cold black walls. The step-by-step procedure is described as follows: Step 1: Initially, an index i = 1 and a number of N p = 500 energy bundles are assigned.
Step 2: The number of bundles emitted in iteration i (starting with i = 1) is computed as N i = N p i.
Step 3: For each iteration, i, the heat transfer problem is solved M = 20 times, generating 20 samples for the wall heat flux, each different one from the other due to statistical oscillations.
Step 4: Using the following statistical relations [20], the wall heat flux is presented within a confidence interval. First the mean heat flux to the wall is computed as where q w,m is the sample for the wall heat flux, and the subscript m is the index of the sample. Following, the standard deviation is computed: and the confidence interval can be obtained as where z is a function which depends on the desired confidence. This function can be obtained from tables [20]. In the computations accomplished in this work the considered confidence is 99%. Finally, the wall heat flux can be presented within a confidence interval as Step 5: Steps 3 and 4 are repeated N s = 100 times, each resulting in a confidence interval somewhat different one from the other due to statistical oscillations. Then, these 100 values for the confidence interval are used to compute an average confidence interval as The N p = 500, M = 20, and N s = 100 were selected by numerical experimentation. If the averaged confidence interval falls below the prescribed convergence criterion (for example, a confidence interval of 1%), the computation stops, and the associated number of sub-bundles and computational time are saved. On the other hand, if the confidence interval is higher than the prescribed value, the simulation index is updated to i = i + 1, and the procedure starts again from step 2.

Results
The methodology described in Sect. 3 was initially applied to the following case.

Case 1
This case considers a one-dimensional slab of gas 1 m thick constituted of 20% of H 2 O at 1600 K, which is confined between two cold parallel black walls. Figure 4 shows the number of bundles necessary to obtain the wall heat flux within a confidence interval of 1.5% with 99% of confidence. As one can see, the statistical oscillations are considerably smaller using the proposed methodology as compared with that appearing in Fig. 2. Similar reductions on statistical oscillations were obtained also for the computational time per bundle, which increases linearly with the number of sub-bundles, as shown in Fig. 5.
Such reduction on the statistical oscillations allows one to exactly determine the optimum number of sub-bundles, which was 3 for this case, as can be seen in Fig. 6. However, quite close computational time savings would be obtained using 4 sub-bundles. Thus, either 3 or 4 sub-bundles would provide the effective best result for the reduction on the computational time.
In addition to the confidence interval of 1.5%, other confidence intervals were considered, namely, 0.8, 1 and 1.2%. It was not verified any significant influence of the confidence interval on the optimum number of sub-bundles, which was always either 3 or 4, both providing quite close reduction on the computational time.
Based on the results one can conclude that 3 sub-bundles minimize the computational time for the case in which the MSB method is applied with the FSK to simulate the heat transfer in a slab of medium 1 m thick composed of 20% H 2 O and non-participating species at 1600 K. It would be useful to find a number of sub-bundles that reduces the computational time to values close to the minimum in general, which could be applied to different cases. Aiming at finding such a number of sub-bundles, simulations were performed in 23 cases, in addition to case 1.

Cases 2-24
For all 23 cases the geometry is one-dimensional. It consists of a gas slab composed of H 2 O and non-participating species confined between two cold black walls. However, temperature, chemical species mole fraction and distance between the walls vary from case to case, ranging respectively from 800 to 2400 K, 5 to 100%, and 0.019 to 41.15 m. The prescribed confidence interval was 1.5% with 99% of confidence.
Similar results presented in Figs. 4, 5 and 6 were obtained for each of the 23 cases. These results are summarized in Table 1. The table also contains information about case 1.  Observing the results presented in Table 1, one can verify that the MSB has lead to reduction on the computational time in all 24 cases. The maximum reduction was about 51%, and only for 3 cases it was less than 20%. The optimum number of sub-bundles varied from 2 to 8, depending on the case.
It is useful to know in advance as to how many subbundles to use for obtaining the maximum reduction on the computational time. However, accordingly Table 1 such a number vary with mole fraction, temperature and geometric dimensions. Nevertheless, a number of sub-bundles which provides significant computational time reduction in general can be sought. Table 2 can be used to assist on determining a number of sub-bundles for general use for analysed class of problem. This table shows the computational time, averaged over all 24 cases in function of the number of sub-bundles. As one can observe, some values in the last column of Table 2 present negative sign. This indicates computational time reduction. As opposite, when the variation is positive the computational time increases. The maximum reduction in average was about 28.5%, which was obtained using either 3 or 4 sub-bundles. However, the number of 3 sub-bundles was the only one which provided reduction on the computational time for all cases, without exception. Therefore 3 sub-bundles is the number suggested for use in MSB/FSK simulations of radiative heat transfer in media composed of H 2 O and non-participating species, within the  limits of the considered test cases, which include temperatures ranging from 800 to 2400 K, and H 2 O mole fraction ranging from 5 to 100%. It is important to observe that the result of 3 sub-bundles obtained with the presented analysis is restricted to the considered problems. Nevertheless the presented procedure could be employed to sought a number of sub-bundles well suited for applications in other problems involving different chemical species, as CO 2 and CO, and/or mixtures.
Homogeneous gases were considered to analyse the effect of temperature, chemical species concentration and geometric dimensions on the optimum number of sub-bundles. In real radiation heat transfer problems homogeneous media are unlikely to occur. The suggested number of 3 sub-bundles is now employed to compute the radiative heat transfer in a nonhomogeneous gas.

Case 25
This case considers a slab of a nonhomogeneous gas, 2 m thick, whose temperature and H 2 O mole fraction vary parabolically respectively according to Eqs. (20) and (21), where temperature is given in Kelvin. Thus temperature varies from 1400 to 2400 K, while H 2 O mole fraction varies from 5 to 25%. The gas is once again confined between two black cold walls. The confidence interval considered for the wall heat flux is 1.5% with 99% of confidence.
Results were obtained using 3 sub-bundles and traditional monochromatic bundles for comparison. Figure 7 presents the divergence of the radiative heat flux, the results differer only due to statistical oscillations. Using the suggested number of 3 sub-bundles, the simulation time was 26.2% lower than that for the traditional Monte Carlo.

Conclusions
A methodology was proposed to accurately determine the optimum number of sub-bundles in MSB Monte Carlo computations. These methodology was applied to a series of test cases which consisted of slabs of homogeneous gases having as participating species H 2 O. The gas temperature and H 2 O mole fraction varied from case to case.
The analysis of the results demonstrated that the optimum number of sub-bundles varied from case to case. Nevertheless a number of sub-bundles which provides significant reduction on the computational time in average could be found. Such a number is 3 sub-bundles and have application to FSK/MSB Monte Carlo computations in gases composed of H 2 O and non-participating species, within the limits of the considered test cases. As an example, the suggested number of 3 sub-bundles was employed to a case of a nonhomogenous slab of gas and provided a reduction on the computational time of 26.2%.