Ray Effects and False Scattering in Improved Discrete Ordinates Method

: The improved discrete ordinates method (IDOM) developed in our previous paper is extended to solve radiative transfer in three-dimensional radiative systems with anisotropic scattering medium. In IDOM, radiative intensities in a large number of new discrete directions are calculated by direct integration of the conventional discrete ordinates method (DOM) results, and radiative heat ﬂux is obtained by integrating radiative intensities in these new discrete directions. Ray effects and false scattering, which tend to compensate each other, are investigated together in IDOM. Results show that IDOM can mitigate both of them effectively with high computation efﬁciency. Finally, the effect of scattering phase function on radiative transfer is studied. Results of radiative heat ﬂux at boundaries containing media with different scattering phase functions are compared and analyzed. This paper indicates that the IDOM can overcome the shortages of the conventional DOM well while inheriting its advantages such as high computation efﬁciency and easy implementation.


Introduction
Efficient and precise solutions for radiative heat transfer are essential for thermal analysis in high-temperature processes, such as combustion, solar energy utilities, etc. Many numerical methods are employed to solve the radiative transfer equation (RTE). Among these, the discrete ordinates method (DOM) is widely used due to its high computation efficiency and easy implementation [1,2]. However, the accuracy of the DOM is affected by ray effects and false scattering [3][4][5], which are related to the angular and spatial discretization, respectively.
The most straightforward approaches for mitigating ray effects and false scattering are to increase discrete directions and refine spatial grids, respectively. Unfortunately, these approaches are computationally expensive due to increased computation time and memory. To use more accurate spatial discretization schemes is another common way to reduce false scattering [6][7][8]. The performance of several schemes including the upwind, central, hybrid, exponential, and high-resolution SMART schemes have been comprehensively compared by Liu et al. [6] to predict radiative intensity as well as radiative heat flux. Jessee and Fiveland [7] applied several bounded, high-resolution schemes to DOM, and the investigated schemes showed obvious improvement in accuracy over the standard step scheme. A later paper by Coelho [8] showed that bounded skew high-resolution schemes had better accuracy compared with the standard high-resolution schemes, yet with higher computational requirements. Besides, the double rays method is developed by Li et al. [9] to reduce false scattering error. In order to treat collimated irradiation, Li additionally used a hybrid spatial differencing scheme to reduce false scattering in arbitrarily specified discrete directions [10].
In terms of ray effects, in addition to increase discrete directions, to modify direction layout is another common remedy [11][12][13][14][15]. Numerical experiments have been conducted by Li et al. [11] showing that ray effects can be reduced by modifying direction layouts, and three corresponding countermeasures were proposed in their paper. Tencer [12] proposed a ray effects mitigation technique by averaging the computed results of various direction layouts, which are generated by arbitrarily rotating the reference layout. Quadrature rotation was also applied to solve time-dependent transport and results showed that ray effects were reduced significantly, including for small numbers of quadrature points [13]. Recently, goal-oriented angular adaptive algorithm was also utilized to mitigate ray effects [14,15]. However, neither increasing discrete directions nor modifying direction layout is a partially effective remedy for ray effects [16,17]. Several modifications to the discrete ordinates scheme were also proposed to eliminate ray effects. Lathrop [16] converted the discrete ordinates equation into a spherical harmonic-like equation to eliminate ray effects. A modification of the discrete ordinates named artificial scattering was proposed to mitigate ray effects [18]. Results demonstrated that artificial scattering yielded better accuracy than the traditional DOM with the same number of discrete ordinates. A modified discrete ordinates method (MDOM) has also been developed to mitigate ray effects caused by non-uniform emission of boundaries [19,20]. Coelho [21] extended the MDOM to reduce ray effects originated from non-uniform emission of the boundaries as well as the media. A similar concept was also applied to YIX method to mitigate ray effects [22]. Huang et al. [23] developed an improved discrete ordinate method (IDOM) based on the discrete ordinates scheme with infinitely small weight (DOS+ISW) concept [24]. Radiative intensities in a large number of new discrete directions are calculated by direct integration of the conventional DOM results. Several cases with isotropic scattering media were studied. Results showed that the IDOM can mitigate ray effects effectively with negligible computation time increased comparing with the conventional DOM.
In the papers mentioned above, ray effects and false scattering are resolved separately with different strategies. However, it is reported that there is an interaction between these two errors and that they tend to compensate each other [21,25,26]. That is to say, reducing one error while keeping the other unchanged may make the accuracy of final results worse. Then, in order to evaluate the performance of strategies applied to the DOM, ray effects and false scattering need to be considered at the same time. However, to the best of our knowledge, there is no previous work proposing or discussing a strategy that mitigates ray effects and false scattering concurrently.
In this paper, the IDOM is extended to solve radiative transfer in three-dimensional systems with anisotropic scattering medium. Ray effects, as well as false scattering, in the IDOM are investigated. First, the description of the IDOM is given. Then, ray effects and false scattering in the IDOM are studied by comparing the calculated radiative heat flux with the benchmark solution provided by the reverse Monte Carlo (RMC) method. Finally, the effect of scattering phase function on radiative transfer is studied by the IDOM.

The Improved Discrete Ordinates Method
In this section, a brief description of the conventional discrete ordinates method is given firstly. Then, the detailed description of the improved discrete ordinates method and the comparison of these two methods in calculation accuracy will follow.

Conventional Discrete Ordinates Method
In a three-dimensional radiative system filled with emitting, absorbing and scattering gray medium, discrete ordinate equations of the DOM with N discrete angular directions are written as where I is radiative intensity, β and S are the extinction coefficient and the source function of the medium, respectively. The source function and the boundary condition can be given as where Φ is the scattering phase function, ω and ε are the scattering albedo of the medium and the emissivity of the wall, respectively. w is the weight of a specific discrete direction. A general control volume p in three-dimensional radiative system is shown in Figure 1. It has six faces labeled as w, e, s, n, f, b. The finite difference form of Equation (1) can be obtained by integrating it over a volume element, where A x = ∆y∆z, A y = ∆x∆z, A z = ∆x∆y, are face areas in different directions, V p = ∆x∆y∆z is the element volume. The volume-averaged intensity, I p,i , is expressed as the weighted average of the surface-averaged intensities, where γ x , γ y , γ z are differential weighting factors, and 0.5 ≤ γ x , γ y , γ z ≤ 1. Substitute Equation (5) into Equation (4), and the volume-averaged intensity is expressed as where S p,i is the volume-averaged radiative source function, Solving Equations (5) and (6) over all N directions and all the elements, radiative intensities at all the volumes and surfaces in N direction is updated. The procedure is repeated until the convergence criteria are met. Then, radiative intensities in N directions at all the volumes and surfaces are obtained. Heat flux at any position in the radiative system can be calculated by integrating the known radiative intensities. For instance, heat flux in z direction at a point can be obtained by integrating radiative intensities at this point as

Improved Discrete Ordinates Method
IDOM is based on the conventional DOM. Radiative intensity results calculated in the conventional DOM are needed in a follow-up calculation of IDOM. Therefore, the first step of IDOM is utilizing the conventional DOM to obtain radiative intensities in N directions at all volumes and surfaces. Then, the second step begins by generating M new directions with finite weights w sk , k = 1,2,..., M. For these M new directions, the discrete ordinate equations can be rewritten as where s is the geometric length. If source function S k , k = 1,2, . . . , M is known, Equation (9) can be solved by direct piecewise integration. Take the k-th new direction at point B and the control volume p for example, as shown in Figure 2, the ray AB intersects the control volume p at the line segment CD, then the integral form of Equation (9) along CD can be given as where ∆s is the geometric length of the line segment CD, β p is the extinction coefficient, and S p,k is the source function of control volume p in k-th new direction. For control volume p, the source function S p,k can be solved by rewriting Equation (2) as For the starting point A of the ray AB, the boundary intensity at point A can be solved by rewriting Equation (3) as where w i is the weight of i discrete direction in the conventional DOM, I i is radiative intensity in i discrete direction and is already known in the first step calculation, i.e., the conventional DOM calculation. With all the source functions and boundary conditions expressed by Equations (11) and (12) are solved, radiative intensity in k-th new direction at point B can be calculated by integrating Equation (9) from point A to point B, where N k and N p is the number of control volumes that the ray AB and CB passes through, respectively. Similarly, radiative intensities in all M directions at point B, I B,k (k = 1, 2, . . . , M), can be obtained by repeating this procedure for each of the M discrete coordinate equations. Then, the heat flux in z direction at point B can be obtained by integrating the radiative intensities in M new directions as

Accuracy Analysis of IDOM
Ray effects and false scattering are two shortcomings of the DOM and cause inaccuracy of radiative intensity results as well as radiative heat flux which is calculated by integrating radiative intensity as shown in Equation (8).
In order to analyze the accuracy of IDOM, the source function in Equation (11) is rewritten as the sum of emission source function S p,k,em and in-scattering source function S p,k,sca , S p,k = S p,k,em + S p,k,sca Similarly, the boundary intensity in Equation (12) is rewritten as the sum of boundary emission intensity I w,k,em and boundary reflection intensity I w,k,re , Then, radiative intensity in k-th new direction at point B can be rewritten as the sum of two terms, where I d B,k is the direct part of radiative intensity contributed by emission of the media and the boundary, I s B,k is the indirect part of radiative intensity contributed by in-scattering of the media and reflection of the boundary. Correspondingly, radiative heat flux calculated by IDOM can also be expressed as the sum of direct heat flux q d IDOM and indirect heat flux As indicated by Equations (15)- (26), the direct part of heat flux is calculated by integrating radiative intensity contributed by media and boundary emission, which can be accurately evaluated when the number of discrete directions (M) generated in the second step of IDOM is large. For the calculation of indirect heat flux in IDOM, since indirect radiative intensity is calculated by path integrating of in-scattering source function and boundary reflection intensity rather than sweeping over all the volume elements in the conventional DOM, the error is accordingly reduced. Overall, ray effects and false scattering are significantly relieved in IDOM, and this will be demonstrated by numerical experiments given in Section 3.

Model Description
The IDOM and the conventional DOM are applied to solve radiative transfer in a three-dimensional cubic radiative system as shown in Figure 2. The side length of the cube is L. Six boundaries are all black, and the enclosed medium is gray with uniform extinction coefficient β. Three different scattering phase functions as described by Equation (27) where θ s is the angle between the directions of incoming and scattered radiative ray. Different values of A 1 , equal to 0, −1, and 1, which represent isotropic, backward and forward scattering, are considered. Three cases with different radiative conditions as listed in Table 1 are examined. In order to examine ray effects caused by non-uniform boundary emission, Case 1 is a radiative system with only one boundary emission which is similar to cases used in refs. [20,22]. Case 2 with a stepwise medium emission which is similar to cases in ref. [22] is used to study ray effects caused by non-uniform medium emission. A radiative system with uniform medium emission is considered in Case 3.  For all calculations in this paper, the medium is discretized into NX × NY × NZ uniform volume elements. NX, NY and NZ are set to be the same value for the cubic systems. Six boundaries are divided in to 6 × NX × NX uniform surface elements. For the angular quadrature in the conventional DOM and the N original directions in IDOM (first step calculation), LSO angular quadrature scheme is employed. For the M new directions in IDOM (second step calculation), the spherical ring angular quadrature scheme [27] is employed and abbreviated as SR. Therefore, IDOM S8 + SR10 means that the N original discrete directions and weights are based on the LSO S8 scheme, and the M new discrete directions and weights are based on the spherical ring angular quadrature scheme with 10 levels. All calculations are performed on a personal computer with an Intel Dual-Core 3.40 GHz CPU, 8.0 GB memory, 64-bit operating system. In all calculations, a single CPU core is employed. Program is written in FORTRAN language.

Ray Effects in IDOM
In this section, the IDOM and DOM are applied to calculate radiative heat flux distributions along the diagonal line of the top boundary for all three cases. Ray effects in both methods are investigated by considering the reverse Monte Carlo method with 2 million energy bundles as the benchmark solution. For all cases considered in this section, NX = 20 is applied to spatial discretization. Figure 3 shows radiative heat flux along the diagonal line of the top boundary in Case 1. In this case, the enclosed medium is cold and purely scattering, and the optical thickness βL is 1.0. The bottom boundary is hot with unity emissive power, while the other boundaries are all cold. As shown in the figure, the difference between heat flux calculated by DOM S8 and RMC is obvious because of ray effects. Comparing with the RMC method, there is still obvious difference even a higher order scheme (DOM S12) is applied. Results calculated by IDOM are in better agreement with the benchmark solution. Take results in Figure 3b for example, the average errors of heat flux by DOM S8 and DOM S12 are about 18%, and the average errors for IDOM S8 + SR10 (M = 530) and S8 + SR20 (M = 2058) reduce to 4.3% and 2.0%, respectively. It shows that the IDOM can reduce ray effects effectively both in isotropic and anisotropic scattering medium. Heat flux at the top boundary in Case 2 is displayed in Figure 4. In this case, the enclosed medium is absorbing and scattering with optical thickness βL equal to 2.0. The bottom half of the medium is hot with a unity emission, while the top half and all boundaries are cold. As displayed in the figure, both DOM S8 and DOM S12 suffer from obvious ray effects due to the non-uniform emission of the medium. Their calculation results are significantly larger than the benchmark solution, and the average errors are larger than 8.2%. In contrast, heat flux calculated by IDOM S8 + SR10 is in good agreement with the benchmark solution and the average error is less than 1.4%. The above results once again show that IDOM can effectively reduce the error of ray effects and significantly improve the accuracy of radiative heat flux.  Figure 5 shows heat flux along the diagonal line of the top boundary in Case 3. All the conditions are the same as those in Case 2, except that the medium has uniform emission in this case. As displayed in the figure, DOM with both S8 and S12 has the same and acceptable accuracy since ray effects is not obvious in the uniform medium. Obvious improvement of the IDOM results can still be observed and the heat flux curves by the IDOM confirm excellently with those of the RMC method. In DOM, the results keep unchanged while the angular discrete scheme changes from S8 to S12. This means the ray effects are not the cause of the error for DOM in this case, and false scattering may be the main reason and this will be investigated in the following section. Furthermore, as displayed in Figure 5d, results of both RMC and IDOM show that heat flux with forward scattering medium is smaller than that with backward scattering medium in the middle of the boundary, while is larger at the edge of the boundary. At first glance, the trend is contradictory to our initial expectation that the forward scattering medium benefits for energy transfer and should have larger radiative heat flux at the boundary. This interesting result, which is related to the effect of scattering phase function on radiative transfer, is to be studied in Section 3.4.

False Scattering in IDOM
As shown in previous section, the IDOM reduces ray effects effectively since radiative intensity in a large number of discrete directions are obtained and used to evaluate radiative heat flux. However, false scattering and ray effects usually have impacts on each other, reducing one of them may increase the other one [21,25,26]. Then, the false scattering in IDOM is to be examined in this section.
In order to decouple the errors caused by ray effects and false scattering on final results, Case 3 with uniform emission medium is considered in which ray effects is not obvious as shown in Figure 5. The enclosed medium is isotropic scattering, and the optical thickness βL is equal to 2.0. Heat flux along the diagonal line of the top boundary with different spatial discrete number NX = 10 and 40 are depicted in Figure 6. From the figure, results calculated by DOM S8 and DOM S12 agree well with each other for all the spatial discrete numbers, which shows ray effects has no influence on results accuracy in these cases. Meanwhile, as the number of NX increases, results of DOM have better accuracy. This means the false scattering is the main cause of the error, and this error is reduced by refining spatial grids.   When the optical thickness is thick, as shown in Figure 7c, although the IDOM has better accuracy than the DOM, it also suffers from false scattering with low spatial resolution (NX = 10). Refining the spatial resolution to NX = 20, the IDOM significantly reduces the false scattering error, and its calculation results are in good agreement with the benchmark solution, while the DOM still has obvious false scattering error, as shown in the Figure 7d. These results indicate that, if the radiative system studied has a thick optical thickness, it is also necessary to refine spatial grids in IDOM to reduce the false scattering error. Meanwhile, the false scattering error in IDOM is much less than that in the conventional DOM if the same spatial grids are employed.

Effect of Scattering Phase Functions on Radiative Transfer
Scattering behavior is one of the key factors affecting radiative transfer. Radiative rays change their original propagating directions by scattering in the medium, and the directional distribution of scattered rays is described by the scattering phase function. In this section, the effect of scattering phase function on radiative energy transfer is analyzed by the IDOM. Radiative heat flux results calculated by IDOM S8 + SR20 with discrete number NX = 20, which have been shown to have good accuracy in Section 3.2, are compared for media with different scattering phase functions. Figure 8 shows heat flux along the diagonal line on boundary of Case 1. As displayed in the figure, heat flux at the top boundary with forward scattering medium is larger than that with backward scattering medium, while the trend of heat flux at the bottom boundary influenced by scattering phase function is reverse. The explanation for these results is straightforward. The forward scattering medium scatters more radiative energy forward; thus, heat flux at boundaries located in the upward direction of emission plane (z = 0) has bigger values, and at boundaries located in the downward direction has smaller values. Heat flux curves along the diagonal line at boundaries in Case 2 are presented in Figure 9. At the top boundary in Figure 9a, heat flux with forward scattering medium is larger than that with backward scattering medium, and this trend is the same as that in Case 1. However, at the bottom boundary in Figure 9b, heat flux at the center area with forward scattering medium is smaller than that with backward scattering medium, while an opposite trend is observed at the edge area. This is because radiative energy received by the bottom boundary is affected by both the backward scattering part (of the energy emitted upward by the source) and the forward scattering part (of the energy emitted downward by the source). These two parts are influenced by the scattering phase functions in an opposite way, and the trend of radiative heat flux at the bottom boundary depends on the overall effects of these two parts. If we revisit the results in Figure 5d, the same trend as in Figure 9b can be observed and can be explained similarly. These results show that unlike the case with emitting boundary, the effect of scattering phase functions on radiative transfer within systems with emitting medium is complicated and needs to be considered case by case.

Computation Cost Comparison
For conventional DOM, its computation time mainly depends on the number of control volumes, discrete directions, and the number of iterations. Take the cubic radiative system in Figure 2 for example, assuming that the number of control volumes is N s 3 , the number of discrete directions is N and the number of iterations is N i , the computation time of the conventional DOM is roughly proportional to N s 3 × N 2 × N i . For IDOM, in addition to the computation time spent in the first step, which is the same as that of the conventional time, extra computation time to make the integral calculation of Equation (13) in the second step is needed. On average, the number of integrated control volumes in each direction is around N s . Assuming the number of new discrete directions is M and heat flux at N q points is to be calculated, the additional computation time of IDOM is roughly proportional to N s × M × N × N q . Since N s × M × N × N q is much smaller than N s 3 × N 2 × N i in most cases, the additional computation time of IDOM is almost negligible compared with the computation time of the conventional DOM.
The above analysis can also be proved by calculation times of different methods listed in Table 2. The computation time is about specific cases considered in Section 3.2. Take Case 1 with isotropic scattering medium as listed in the first line of the table for instance, calculation time consumed by DOM increases from 1.22 s to 5.88 s with angular discretization number increases from 80 (S8) to 168 (S12), while the average errors of heat flux by DOM S8 and DOM S12 are about 16% as shown in Figure 3a. For IDOM S8 + SR10, its calculation time is only about 0.08 s more than that of DOM S8, while the average error of heat flux reduces to about 3.8%. This error in IDOM can be further reduced by using more new directions (such as IDOM S8 + SR20) with negligible increased computation time. There are similar conclusions in the other cases. That is to say, IDOM can provide results with much higher accuracy while keeping computation time almost the same as the conventional DOM.
Furthermore, since the second step in IDOM is merely an integral calculation, the memory usage for new variables in IDOM is limited and negligible compared with that of the conventional DOM.

Conclusions
The IDOM is extended to calculate radiative heat flux in three-dimensional radiative systems with anisotropic scattering medium. Radiative heat flux calculated in IDOM is divided into the direct part contributed by emission of the media and the boundaries and the indirect part contributed by in-scattering of the media and reflection of the boundaries. Theoretical analysis indicates that IDOM has higher accuracy than the conventional DOM.
Ray effects, together with false scattering in IDOM, were investigated in this paper by considering the reverse Monte Carlo method as the benchmark solution. In the considered radiative system with only one boundary emission, which suffers obvious ray effects, the average error of heat flux was about 16%, and the calculation time was about 1.22 s for conventional DOM. For IDOM, the average error of heat flux was less than 4% and the calculation time was only 0.08 s more than that of the conventional DOM. More numerical experiments show that the IDOM could reduce ray effects effectively both in isotropic and anisotropic scattering medium with negligible increased computation time compared with the conventional DOM. Meanwhile, in a radiative system with uniform medium emission of which the accuracy is not influenced by ray effects, the calculation error of conventional DOM was larger than 3% while that of IDOM is less than 1%, when the same spatial resolution was applied in both methods. This numerical experiment proves that the false scattering error in IDOM was much smaller than that in the conventional DOM. In a word, the IDOM can overcome the shortages of the conventional DOM well, while inheriting its advantages such as high computation efficiency and easy implementation.
The IDOM was also applied to investigate the effect of scattering phase functions on radiative transfer. Results show that the trend of radiative heat flux with different scattering phase functions is simple for cases of emitting boundaries, while it is complicated for cases of emitting media.