Adaptability analysis of radiative transport diffusion approximation in planar-graded-index media

Abstract

To study the applicability of diffusion approximation method, the influence of the optical thickness, scattering albedo, and wall emissivity of the examined medium on the accuracy of the diffusion approximation is analyzed. The accuracy is determined by comparison with results obtained by using the Monte Carlo method. By calculating the dimensionless radiation heat flux distributions at the bottom of the medium for various parameter combinations, the influence of each optical parameter is analyzed. In addition, the adaptability of the diffusion approximation method to the case of a gradient-index medium is examined. The results show that in the homogeneous-refractive-index medium and the gradient-index medium, the accuracy of diffusion approximation method is improved by an increase in the optical thickness or the scattering albedo and by a decrease in the wall emissivity. Moreover, the maximum relative errors are primarily distributed at the medium boundary and center. In the case of increased optical thickness or decreased wall emissivity, the error of the diffusion approximation method becomes stable and remains at a small value.

Introduction

Radiation transfer, as a basic transfer mode, occupies an important position in the field of industrial applications. The research is not only limited to the field of solar radiation transmission ¹ and high-temperature gas radiation ² but also the micro-scale radiative heat transfer, ³ spacecraft thermal control design, ⁴ and radiation inverse problem method to identify the temperature field of target ⁵ and so on. Obviously, the study of radiative transfer is still of great significance.

At present, the calculation of the radiative transfer equation mainly depends on different numerical methods, such as Monte Carlo method (MCM), ⁶ discrete transfer method, ⁷ discontinuous finite element method, ⁸ natural element method, ⁹ and so on. However, these methods have different applicable conditions. For instances, the MCM can deal with all kinds of complex problems, the complexity of the simulation calculation increases proportionately with the complexity of the problem and the amount of calculation increases greatly when dealing with the radiation of non-uniform medium. The discrete transfer method is considered to be very convenient for dealing with homogeneous and inhomogeneous medium radiation, but it is not easy to converge when solving the energy equation. The discontinuous finite element method has strong adaptive ability and high accuracy, but the computation is still slightly larger.

The P1 diffusion approximation (DA) method, which is based on spherical harmonics, eliminates the angle coordinates and also reduces the dimension processing, which greatly simplifies this problem. In recent years, the DA method is mainly applied to the field of medical imaging because of its fast computation. Ruan ¹⁰ used diffusion equation to solve the forward problem in the shape diffusion optical tomography technology. Furthermore, Xu et al. ¹¹ investigated total time-resolve diffusion reflectance from a semi-infinite turbid media based on the DA method. For the research on the gradient-index medium, to our knowledge, Liu and colleagues^12,13 have derived the DA equation of steady-state radiative transfer in multidimensional gradient-index medium, and the application conditions of the DA method are briefly described, so that the degree of influence of each optical parameter has not been further revealed.

In this work, the impact of τ_L, ε, and ω on the accuracy of the DA method, for media having a homogeneous refractive index and a gradient index distribution, is analyzed and compared. By calculating the distribution of the radiation heat flux at the bottom of the medium, the degree of influence of each optical parameter is determined. Note that, in order to simplify this work, the two-dimensional (2D) steady-state radiative transfer equation is considered.

Method

The radiative transfer equation for participating media can be expressed as

$$\begin{array}{c}{n}^2\frac{d}{\mathrm{ds}}\left(\frac{I(r,\mathrm{\Omega })}{n^2}\right)=n^{2}k{I}_b(r)-(k+\sigma s)I(r,\mathrm{\Omega })\hfill \\ +\frac{{\sigma }_s}{4\pi }\underset{4\pi }{\int }I(r,\mathrm{\Omega }\prime )\mathrm{\Phi }(\mathrm{\Omega }\prime ,\mathrm{\Omega })d\mathrm{\Omega }\prime \hfill \end{array}$$

(1)

where $I(r,\mathrm{\Omega })$ is the radiation intensity, n is the refractive index of the medium, s is the geometric path length, k is the absorbing coefficient, ${\sigma }_s$ is the scattering coefficient, $I_b(r)$ is the blackbody radiation intensity, and $\mathrm{\Phi }(\mathrm{\Omega }\prime ,\mathrm{\Omega })$ is the scattering phase function.

Based on P1 DA method, I can be rewritten in terms of the incident radiation G, and radiation heat flux q can be written as

$$I(r,\mathrm{\Omega })=\frac{1}{4\pi }[G(r)+3q(r)·\hat{s}]$$

(2)

Combined with formulas (1) and (2), the steady-state radiative transfer equation can be obtained

$$\begin{array}{c}-\nabla ·[\frac{1}{\beta [1-\frac{A_1\omega }{3}]}\nabla G(r)]+\nabla \left[\frac{1}{\beta [1-\frac{A_1\omega }{3}]}\frac{\nabla n^2}{n^2}\right]\hfill \\ +3\beta (1-\omega )[G(r)-4\pi n^2{I}_b(r)]=0\hfill \end{array}$$

(3)

The boundary condition are expressed as

$$G=\frac{2-\epsilon }{\epsilon }2\mathrm{q}·n_w\pi n^2{I}_{\mathrm{bw}}$$

(4)

where $n_w$ is the outer normal vector of the wall.

Using the conventional Galerkin finite element method, finally we can obtain

$$\begin{array}{c}\sum _i\{\underset{V}{\int }\nabla {\mathrm{\Gamma }}_i(r)·\nabla {\mathrm{\Gamma }}_j(r)\mathrm{dV}+\underset{V}{\int }(\nabla \ln n^2(r)·\nabla {\mathrm{\Gamma }}_i(r)){\mathrm{\Gamma }}_j(r)\mathrm{dV}\hfill \\ +\underset{V}{\int }(k[\beta (3-A_1\omega )]+{\nabla }^2\ln n^2(r)){\mathrm{\Gamma }}_i(r){\mathrm{\Gamma }}_j(r)\mathrm{dV}\hfill \\ -\underset{\partial V}{\int }(\frac{\nabla n^2}{n^2}·n_w+\frac{3}{2}\frac{{\epsilon }_w}{2-{\epsilon }_w}\beta \prime ){\mathrm{\Gamma }}_i(r){\mathrm{\Gamma }}_j(r)\mathrm{dA}\}{G}_i\hfill \\ =\underset{V}{\int }4\pi [\beta (3-A_1\omega )]k{n}^2{I}_b(r){\mathrm{\Gamma }}_j(r)\mathrm{dV}\hfill \\ -\underset{\partial V}{\int }\frac{3}{2}\frac{{\epsilon }_w}{2-{\epsilon }_w}\beta \prime 4\pi n^2{I}_{\mathrm{bw}}{\mathrm{\Gamma }}_j(r)\mathrm{dA}\hfill \end{array}$$

(5)

where ${\mathrm{\Gamma }}_i$ is the shape function, V and A represent the computing domains and the area, respectively, and $\beta$ is the extinction coefficient.

The discrete formula (5) can be written in a matrix form

$$K_{\mathrm{ji}}{G}_i=H_j$$

(6)

Each matrix element is defined as

$$\begin{array}{l}{K}_{ji}={\displaystyle \sum _i\{{\displaystyle \underset{V}{\int }\nabla {\mathrm{\Gamma }}_i(r)}\cdot \nabla {\mathrm{\Gamma }}_j(r)dV+{\displaystyle \underset{V}{\int }(\nabla \ln n^2(r)}\cdot \nabla {\mathrm{\Gamma }}_i(r)){\mathrm{\Gamma }}_j(r)dV}\\ +{\displaystyle \underset{V}{\int }(k[\beta (3-A_1\omega )]+{\nabla }^2\ln n^2(r)){\mathrm{\Gamma }}_i(r)}{\mathrm{\Gamma }}_j(r)dV\\ -{\displaystyle \underset{\partial V}{\int }\left(\frac{\nabla n^2}{n^2}\cdot n_w+\frac{3}{2}\frac{{\epsilon }_w}{2-{\epsilon }_w}\beta \prime \right){\mathrm{\Gamma }}_i(r)}{\mathrm{\Gamma }}_j(r)dA\}\end{array}$$

(7)

$$\begin{array}{c}\multicolumn{1}{c}{H_j=\underset{V}{\int }4\pi [\beta (3-A_1\omega )]k{n}^2{I}_b(r){\mathrm{\Gamma }}_j(r)\mathrm{dV}}\\ -\underset{\partial V}{\int }\frac{3}{2}\frac{{\epsilon }_w}{2-{\epsilon }_w}\beta \prime 4\pi n^2{I}_{\mathrm{bw}}{\mathrm{\Gamma }}_j(r)\mathrm{dA}\hfill \end{array}$$

(8)

Results and discussion

To determine the performance and accuracy of the presented DA method, we examine radiative heat transfer problems in a rectangle with opaque and diffuse surfaces. The boundaries are cold (0°K), and the medium temperature is maintained at 1000°K. For the following numerical study, the medium is uniformly divided into 200 × 200 units. Depending on the τ_L, ε, ω, and n of the material, different results can be obtained.

Several sample problems are tested to explore the accuracy and computation efficiency of P1 DA for different τ_L, ε, and ω. In subsection “Square enclosure filled with homogeneous-refractive-index medium,” three examples regarding homogeneous-refractive-index media are presented. In subsection “Square enclosure filled with graded-index medium,” examples concerning gradient-index medium are considered.

Square enclosure filled with homogeneous-refractive-index medium

The DA algorithm is tested through comparison with the MCM. The examined medium has a fixed $n(x)=5$ . For the case of $\epsilon =1.0$ and $\omega =0.5$ , the dimensionless heat flux distributions at ${\tau }_L=1.0,5.0,\mathrm{and}10.0$ are shown in Figure 1(a)–(c), respectively. It is apparent that the DA results near the boundaries are larger than those obtained by using the MCM when τ_L is small and that the relative error is approximately 12%. Furthermore, the deviation near the center is small, with a maximum relative error of only 3.4%. With increasedτ_L $({\tau }_L=5.0)$ , the error is primarily distributed near the center of the medium, and the maximum relative error is approximately 6%. Finally, for ${\tau }_L=10.0$ , the maximum relative error is less than 5.5%. Overall, the DA method produces large errors in optically thin media, and these errors decrease as τ_L increases. Afterτ_L reaches a certain value, the relative error is no longer significantly reduced with further increases to the thickness.

Figure 1.

Dimensionless radiation heat flux distributions along the bottom wall for $\epsilon =1.0$ , $\omega =0.5$ and (a)–(c) ${\tau }_L=1.0,5.0,\mathrm{and}10.0$ , respectively.

For $\epsilon =0.5$ , ${\tau }_L=10.0$ , and different ω, the resultant dimensionless heat flux distribution are presented in Figure 2. This figure demonstrates that the maximum relative error appears at the wall in this case. Under the condition of $\omega =0$ , the maximum relative error of the DA method is 6.3%, which is higher than that obtained for $\omega =0.5$ . The reason is that high scattering albedo makes wall radiation intensity with angle distribution more uniform.

Figure 2.

Dimensionless radiation heat flux distributions along the bottom wall for $\epsilon =0.5$ , ${\tau }_L=10.0$ , and (a) $\omega =0.5$ and (b) ω = 0.

In Figure 3, the distributions of the dimensionless heat flux as obtained using both the DA method and the MCM are plotted for the case of ${\tau }_L=10.0$ , $\epsilon =0.1$ , and $\omega =0.5$ . It can be observed from Figures 1(c), 2(a), and 3 that, for the same τ_L and ω and with $\epsilon =1.0,0.5,\mathrm{and}0.1$ , respectively, the maximum relative errors are 5.5%, 3%, and 0.53%, respectively. Therefore, the DA method exhibits higher accuracy for small ε. This is because the lower the emissivity, the lower the absorption rate, and the more uniform the radiation intensity distribution.

Figure 3.

Dimensionless radiation heat flux distribution on homogeneous-refractive-index medium bottom for ${\tau }_L=10.0$ , $\epsilon =0.1$ , and $\omega =0.5$ .

Square enclosure filled with graded-index medium

Next, we consider a special layered graded-index medium, where the refractive index is modeled by

$$n(x)=5^{*}{(1-{0.95}^{2*}{\mathrm{x}}^2)}^{0.5}$$

(9)

Thus, this medium has a variable refractive index n(x). We set ${\tau }_L=1.0$ , $\epsilon =1.0$ , and $\omega =0.5$ , and the resultant dimensionless heat flux distributions are shown in Figure 4. Hence, it is apparent that the dimensionless heat flux density is asymmetric and the peak moves to the left in the gradient-index medium, and the uniform radiation intensity distribution cause the error at the left wall to increase, while that at the right wall decrease.

Figure 4.

Dimensionless radiation heat flux distribution on graded-index medium bottom for ${\tau }_L=1.0$ , $\epsilon =1.0$ , and $\omega =0.5$ .

Figure 5(a) and (b) shows the dimensionless heat flux distributions for ${\tau }_L=5.0\mathrm{and}10.0$ , respectively. The other conditions are as above. Comparison of these results with Figure 4, indicates that, for ${\tau }_L=1.0$ , the DA method is inaccurate at the left wall, and the maximum relative error is 23%. However, with increased τ_L, the left-wall error decreases gradually. The maximum relative error of ${\tau }_L=5.0$ is approximately 5%. However, with further increases in τ_L, for example, in the case of ${\tau }_L=10.0$ , the maximum relative error remains close to 5%. Thus, a threshold optical thickness τ_L0 exists, in the case of τ _L <τ_L0, the DA method exhibits low accuracy; otherwise, for higher accuracy is achieved.

Figure 5.

Dimensionless radiation heat flux distributions along the bottom wall for $\epsilon =1.0$ , $\omega =0.5$ , and (a) ${\tau }_L=5.0$ and (b) ${\tau }_L=10.0$ .

In Figure 6, the dimensionless heat flux distributions are plotted for ${\tau }_L=10.0$ , $\epsilon =0.5$ , and different ω values of (a) 0.5 and (b) 0. It can observed that an increased value of ω can improve the DA method accuracy. That is, the maximum relative error for $\omega =0.5$ is 2.79%, whereas that for $\omega =0$ is 3.61%. Because ${\tau }_L=10.0$ in this case, the DA method has higher accuracy, and the relative error is not large. Thus, for an optically thick media, increasing ω can improve the accuracy of DA method; however, the effect of the increased ω is not significant.

Figure 6.

Dimensionless radiation heat flux distributions along the bottom wall for ${\tau }_L=10.0$ , $\epsilon =0.5$ and (a) $\omega =0.5$ and (b) $\omega =0$ .

Figure 7 shows the distributions of dimensionless heat flux under the conditions of ${\tau }_L=10.0$ , $\epsilon =0.1$ , and $\omega =0.5$ . Comparison of this figure with Figures 5(b) and 6(a) demonstrate that, for the sameτ_L and ω, the maximum relative errors of DA method are approximately 5%, 2.8%, and 2.8% for $\epsilon =1.0,0.5,\mathrm{and}0.1$ , respectively. Although the relative errors obtained from Figures 6(a) and 7 are essentially identical, some error is apparent near the medium center in Figure 6(a). In contrast, the DA method has a significant error in the small range near the left wall surface of the medium only in Figure (7), this is in good agreement with the benchmark solution in other ranges.

Figure 7.

Dimensionless radiation heat flux distributions on graded-index medium bottom for ${\tau }_L=10.0$ , $\epsilon =0.1$ , and $\omega =0.5$ .

Conclusion

The DA method, which is a rapid algorithm, was extended to solve radiative transfer problems in a 2D rectangle containing a participating medium. Then, the applicability and precision of this method for a variety of optical parameters were examined. Under appropriate boundary conditions, the DA method was found to be consistent with the MCM, and the errors were primarily distributed at the boundaries and the medium center. Comparison of the influence of the optical parameters on the DA method in case of a homogeneous-medium-refractive-index and a graded-index medium indicated the accuracy of the DA. That is, the accuracy of the DA method increase with increased optical thickness or scattering albedo and with decreased wall emissivity. When the optical thickness increases or the wall emissivity decreases to a certain value, the error of the DA method becomes stable and keeps a small value.