Solution of Pure Scattering Radiation Transport Equation (RTE) Using Finite Difference Method (FDM)

Abstract

Radiative transfer is the physical phenomenon of energy transfer in the form of electromagnetic radiation. The propagation of radiation through a medium is affected by absorption, emission, and scattering. Radiative Transfer Equation (RTE) have been applied in a many subjects including optics, astrophysics, atmospheric science, remote sensing, etc. Analytic solutions for RTE exist for simple cases, but, for more realistic media with complex multiple scattering effects, numerical methods are required. In the RTE, six different independent variables define the radiance at any spatial and temporal point. By making appropriate assumptions about the behavior of photons in a scattering medium, the number of independent variables can be reduced. These assumptions lead to the diffusion theory (or diffusion equation) for photon transport. In this work, the diffusive form of RTE is discretized, using a Forward-Time Central-Space (FTCS) Finite Difference Method (FDM). The results reveal the radiance penetration according to Beer-Lambert law.

1 Introduction

Radiative transfer is the physical phenomenon of energy transfer in the form of electromagnetic radiation. Radiative transfer has applications in a wide variety of subjects including optics, astrophysics, oceanography, atmospheric science, remote sensing, infra-red imaging, etc. [3, 11, 16].

The propagation of radiation through a medium is affected by absorption, emission, and scattering processes [3]. Absorption is the process by which the energy of the photons is transferred to particles such as electrons present in its transmission medium. The transfer of momentum raises the localized kinetic energy of the electrons and hence the temperature. A common daily life example of absorption is that on a sunny day we can feel the warmth of the sun, even if the surrounding temperature is quite low.

Emission is the process through which radiation transfers in the form of waves/particles. Emissions may originate from a point source, such as a bulb filament or a spark, or from a surface such as an ionized tube, neon bulb, etc.

Scattering is the deviation of the radiation waves/particles from their original path. Scattering occurs as a result of particle-particle collisions. It can be categorized as specular, such as in the case of mirror or diffused reflections [2].

Radiative transfer can be expressed mathematically in the form of radiation transport equation [4]. The equation poses a challenge for those wishing to obtain a definite solution considering its partial differential nature. Therefore, various studies have employed different methodologies. For example, some researchers have used the Monte-Carlo method to solve the radiation transport models [6]. Similarly, other numerical techniques such as the discrete-ordinate method have also been employed to find the solution [8].

This paper focuses on solving the pure scattering radiation transport equation using the finite difference method. This methodology has previously been used to solve the heat equation and simulate an infra-red signature [9, 10, 14].

2 Radiation Transport Equation

The equation of radiative transfer [13, 17] is given mathematically as shown in Eq. (1),

(1)

where $I_{v}$ the is spectral radiance of electromagnetic waves, $c$ is the speed of light, $\hat{\varOmega }$ is the vectorial position of a solid angle (polar angle $\theta$ and azimuthal angle $\varphi$), $\varOmega$ is a solid angle, $k_{v,s}$ is the scattering opacity of the medium, $k_{v,a}$ is the absorption opacity of the medium, $j_{v}$ is the emission coefficient of the medium and $t$ is the time variable.

The energy contents of electromagnetic waves can be calculated as shown in Eq. (2),

$$\begin{aligned} dE_v = I_v\left( \mathbf {r},\hat{\mathbf {n}},t \right) \cos (\theta )\ dv\ dA\ d\varOmega \ dt \end{aligned}$$

(2)

where $E_{v}$ is the radiation energy, $\theta$ is an angle that the unit vector $\hat{\mathbf {n}}$ makes with the normal of elemental area $dA$, positioned at $\mathbf {r}$, and $v$ is the frequency. This is illustrated in Fig. 1.

Fig. 1.

Radiance energy originating from differential area $dA$ at vectorial position $\mathbf {r}$ within a solid angle $d\varOmega$.

Similarly, radiance intensity and energy flux can be written as shown in Eqs. (3) and (4):

$$\begin{aligned} \phi _{v}\left( \mathbf {r},t \right) = \int _{4\pi }^{\ }{I_{v}\left( \mathbf {r},\hat{\mathbf {n}},t \right) }\ d\varOmega \end{aligned}$$

(3)

$$\begin{aligned} \psi _{v}\left( \mathbf {r},t \right) = \int _{4\pi }^{\ }{\hat{\mathbf {n}}\mathbf {\ }I_{v}\left( \mathbf {r},\hat{\mathbf {n}},t \right) }\ d\varOmega \end{aligned}$$

(4)

where $\phi _{v}$ and $\psi _{v}$ are radiance intensity and energy flux respectively. Each of these terms has units of W/m².

Analytical solutions to the radiative transfer equation (RTE) exist for simple cases, but, for more realistic mediums and complex multiple scattering effects, computer simulations are required.

3 Methodology

There are six independent variables defining the radiance at any spatial and temporal point in the radiation transport equation. These variables are the $x$, $y$, $z$ coordinates from a reference presented in the form of vector $\mathbf {r}$, two dimensional vectorial position of a solid angle $\hat{\varOmega }$ (polar angle $\theta$ and azimuthal angle $\varphi$) and time variable $t$. By making appropriate assumptions about the behavior of the photons in the scattering medium, the number of independent variables can be reduced. These assumptions lead to the diffusion theory (and diffusion equation) for photon transport. Two assumptions, which permit the application of diffusion theory are:

• 1. Relative to scattering events, there are very few absorption events. Likewise, after numerous scattering events, few absorption events will occur, and the radiance will become nearly isotropic. This assumption is sometimes called directional broadening.

• 2. In a primarily scattering medium, the time for substantial current density change is much longer than the time to traverse one transport mean free path. Thus, over one transport mean free path, the fractional change in current density is much less than unity. This property is sometimes called temporal broadening.

It should be noted that both of these assumptions require a high-albedo (predominantly scattering) medium. The diffusion approximation is limited to systems where reduced scattering coefficients are much larger than their absorption coefficients and have a minimum layer thickness of the order of a few transport mean free paths.

From the diffusion approximation, we can write Eq. (5),

$$\begin{aligned} I_{v}\left( \mathbf {r},\hat{\mathbf {n}},t\right) = \frac{1}{4\pi }\phi _{v}\left( \mathbf {r},t \right) + \ \frac{3}{4\pi }\ \psi _{v}\left( \mathbf {r},t\right) .\hat{\mathbf {n}}\ \end{aligned}$$

(5)

by substituting Eq. (5) in Eq. (1), we get Eq. (6),

$$\begin{aligned} \frac{\partial \phi _{v}\left( \mathbf {r},t \right) }{c\ \partial t}\ + k_{v,a}\phi _{v}\left( \mathbf {r},t \right) \ - \nabla .\psi _{v}\left( \mathbf {r},t \right) = j_{v}\left( \mathbf {r},t \right) \ \end{aligned}$$

(6)

by applying Fick’s Law [7], we get Eq. (7),

$$\begin{aligned} \psi _{v}\left( \mathbf {r},t \right) = - \frac{\nabla \phi _{v}\left( \mathbf {r},t \right) }{3\left( \left( 1 - g \right) k_{v,s} + k_{v,a} \right) } = - D\nabla \phi _{v}\left( \mathbf {r},t \right) \ \end{aligned}$$

(7)

where g is the anisotropy of the medium and D is the diffusion coefficient.

By substituting Eq. (7) in Eq. (6), we get Eq. (8),

$$\begin{aligned} \frac{\partial \phi _{v}\left( \mathbf {r},t \right) }{c\ \partial t} + k_{v,a}\phi _{v}\left( \mathbf {r},t \right) \ - D\nabla ^{2}\phi _{v}\left( \mathbf {r},t \right) = j_{v}\left( \mathbf {r},t \right) \ \end{aligned}$$

(8)

Assuming hypothetically that there are zero absorption and zero emission (pure scattering medium), we can simplify Eq. (8) to Eq. (9),

$$\begin{aligned} \frac{\partial \phi _{v}\left( \mathbf {r},t \right) }{c\ \partial t} = D\nabla ^{2}\phi _{v}\left( \mathbf {r},t \right) = D\left( \frac{\partial ^{2}\phi _{v}\left( \mathbf {r},t \right) }{\partial x^{2}} + \frac{\partial ^{2}\phi _{v}\left( \mathbf {r},t \right) }{\partial y^{2}} + \frac{\partial ^{2}\phi _{v}\left( \mathbf {r},t \right) }{\partial z^{2}} \right) \ \end{aligned}$$

(9)

where $x,\ y,\ z$ are the space dimensions. We can reduce the dimensions to a hypothetical two-dimensional space. This will further simplify the equation as shown in Eq. (10),

$$\begin{aligned} \frac{\partial \phi _{v}\left( \mathbf {r},t \right) }{c\ \partial t} = D\left( \frac{\partial ^{2}\phi _{v}\left( \mathbf {r},t \right) }{\partial x^{2}} + \frac{\partial ^{2}\phi _{v}\left( \mathbf {r},t \right) }{\partial y^{2}} \right) \ \end{aligned}$$

(10)

In order to do so, we have to discretize the equation. In this work, we will discretize the equation using a Forward-Time Central-Space (FTCS) Finite Difference Method (FDM) [1, 10, 12, 14]. This results in Eq. (11),

(11)

where subscripts i, j are integers representing the computational points in the two-dimensional space domain and the superscript represents the transient state.

A finite difference method (FDM) is a numerical method for solving differential equations such as that in Eq. (10). This method approximates the differentials with differences by discretizing the dependent variable (radiance intensity) in the independent variable domains (space and time). Each discretized value of the dependent variable is referred to as a nodal value.

Equation (11) is solved in a two-dimensional spatial domain, as shown in Fig. 2. The two-dimensional space is discretized in equally spaced quadrilaterals. Each quadrilateral is referenced in two-dimensional space via indices $i$ and $j$. Indices $i$ and $j$ refer to positions on the horizontal and vertical axes, respectively. The $n$ and $m$ refer to the maximum value of the $i$ and $j$ indices, respectively.

Fig. 2.

Radiance intensities in the two-dimensional discretized domain. Indices $i$ and $j$ refer to the nodal position of radiance intensity.

For the stability and accuracy of the FDM, it is vital to choose the correct time step value. In this work, the Courant-Friedrichs-Lewy (CFL) condition [5, 12] is used to decide the time step size. The CFL condition is given in Eq. (12),

$$\begin{aligned} 2D c \varDelta t\ \le min\left( \left( \varDelta x \right) ^{2},\left( \varDelta y \right) ^{2} \right) \ \end{aligned}$$

(12)

where $D$ is the diffusion coefficient (m²/s), $c$ is the speed of light (m/s), $\varDelta t$ is the time step size (s), and $\varDelta x$ and $\varDelta y$ are the differences in the spatial positions of the nodes (m).

In addition, initial and boundary conditions are required. In this case, boundary conditions were specified such that they resembles to an infinite space. The source radiance was introduced on a small part of the boundary. The flow chart of the method of solution is given in Fig. 3, while the values of the constants are presented in Table 1.

Fig. 3.

Flow chart of the method of solution.

Table 1. Values of constants

4 Results and Discussion

Ten cases are set with inlet radiance intensities of 10 W/m², 50 W/m², 100 W/m², 500 W/m², 1000 W/m², 5000 W/m², 10000 W/m², 50000 W/m², 100000 W/m², and 500000 W/m². The main reason for this diverse range is to understand the response of the medium and to identify its limiting behavior. It is vital to highlight here that the medium’s spectral diffusivity plays a vital role. Given study is limited to only highly diffusive mediums. Figure 4 shows the radiance intensity penetration with different inlet radiance intensities.

Fig. 4.

Radiance intensity penetration

Fig. 5.

Contour plots with various radiance intensities.

Contour plots with various inlet radiance intensities are shown in Fig. 5. Lighten regions show the radiance intensity of 100 W/m² and above in a two-dimensional space. As expected, the results clearly demonstrate that the radiance penetration has increased with the increased value of the inlet source radiance. This behavior is also in accordance with the Beer-Lambert Law [15].

5 Conclusion

The Radiation Transport Equation (RTE) can be solved using a Forward-Time Central-Space (FTCS) Finite Difference Method (FDM) in special cases such as a highly diffusive mediums (as discussed above). The results show the radiance intensity space, which reflects on the radiance penetration.

6 Future Work

This work can be further extended by varying the diffusion coefficient and observing its impact on the radiance penetration. It will help to find the validity limits of Radiation Transport Equation (RTE) solution using Forward-Time Central-Space (FTCS) Finite Difference Method (FDM).