Solution of Linear Radiative Transfer Equation in Hollow Sphere by Diamond Difference Discrete Ordinates and Decomposition Methods

Schramm M.; Ladeia C. A.; Fernandes J. C. L.;

DOI 10.5433/1679-0375.2024.v45.51961

Citation Semin., Ciênc. Exatas Tecnol. 2024, v. 45: e51961

Received: 29 November 2024 Received in revised for: 22 December 2024 Accepted: 26 December 2024 Available online: 28 December 2024

Abstract:

In this article, we present a methodology to solve radiative transfer problems in spherical geometry without other forms of heat exchange. We use a decomposition method based on the Adomian formulations, together with a diamond difference scheme and a trapezoidal rule to approximate the integral part of the solution. The algorithm is simple, highly reproducible, and can be easily adapted to further problems or geometries. Also, we demonstrate its consistency and showed that using an analytical solution with a trapezoidal rule improves the order of convergence compared to using the finite difference method. These considerations are necessary for future applications in more complex cases. The numerical results are compared with some classical and recent cases in the literature, along with a simplified version of a complete (fully coupled with heat exchange) case.

Keywords: radiativetransfer, spherical geometry, decomposition method, diamond difference, discrete ordinates

Introduction

Extensive research on the radiative transfer equation in spherical geometry has been developed over the last decades (Howell et al., 2016; Stamnes et al., 2017). The present discussion considers the radiative transfer problem in hollow spheres. Solutions found in the literature are typically determined by numerical methods; see, for instance, Abulwafa, 1993; Ladeia et al., 2020; Sghaier, 2013.

There are few recent papers about solving the radiative transfer equation. Xu et al., 2023 used an approximate technique for remote sensing and Li et al., 2020 used a moment-based method to solve the non-linear transport equation, but there is little effort in showing numerical analysis like consistency and convergence.

In this work, we will present a hybrid semi-analytical methodology with a focus on formalism, thus standing out from other works in the literature that do not present such necessary formalism for iterative, recursive, and discretization methods. This methodology is characterized by the combination of the methods of discrete ordinates and diamond difference in the treatment of the angular variable and a decomposition method based on Adomian, 1988 to solve the resulting system of ordinary differential equations.

Many researchers have been using the decomposition technique for many fields of science and technology, for example, Allahviranloo, 2005; Haq et al., 2018; Wazwaza & El-Sayed, 2001. For instance, in transport theory, Ladeia et al., 2019 used the modified decomposition method for nonlinear radiative conductive transfer problem in solid sphere and Segatto et al., 2017 solved the multi-group neutron transport equation in slab geometry using a combination of the Adomian decomposition method and the discrete ordinates technique.

In the present discussion, we demonstrate the consistency of the solution obtained from the recursive scheme. Also, we show the order of convergence of part of our methodology and compare with the finite difference method. Furthermore, we report cases with numerical solutions that are compared with data in the literature (Abulwafa, 1993; Ladeia et al., 2020). Last, we make concluding remarks about the achievements in developing this research together with future perspectives.

Materials and methods

We begin considering the radiative transfer equation in spherical geometry for hollow spheres (Ozisik, 1973), $$\begin{gathered} \label{1} \frac\mu{r^2} \frac \partial{\partial r} \left[ \vphantom{Q_Q^Q} r^2 I(r,\mu) \right] + \frac1r \frac \partial {\partial \mu } \left[ \vphantom{Q_Q^Q} \left(1-\mu^2 \right)I(r,\mu) \right] + I(r,\mu) \\ = \left( \vphantom{Q_Q^Q} 1- \omega\left(r\right)\right) I_b\left(T\right) + \frac {\omega\left(r \right)} 2 \int_{-1}^1 p\left(\mu,\mu^\prime \right) I\left(r , \mu^\prime\right) \text{d}\mu^\prime\, , \end{gathered}$$ where $R_1 \leq r \leq R_2$ is the optical space variable and $-1 \leq \mu\leq 1$ is the direction cosine variable. $R_1$ and $R_2$ are the inner and outer spherical surfaces radii, respectively. Further, $I(r,\mu)$ is the radiation intensity, $I_b\left(T\right)$ is the black body radiation for temperature $T$, $\omega$ is the single scattering albedo and $p\left(\mu,\mu^\prime \right)$ is the phase function. According to Chandrasekhar, 1950, $p\left(\mu,\mu^\prime \right)$ may be expressed in terms of Legendre

polynomials, $$\begin{aligned} p\left(\mu,\mu^\prime \right) = \sum _{\ell=0}^{L} \beta_\ell \mathcal P_\ell \left(\mu\right) \mathcal P_\ell \left(\mu^\prime\right) \, , \end{aligned}$$ where $\beta_\ell$ are the expansion coefficients of the Legendre polynomials and $\mathcal P_\ell$ is the $\ell$-th Legendre polynomial. Here, $L$ refers to the degree of anisotropy. The boundary conditions of equation (1) are $$I(R_1, \mu) =\epsilon_1 I_{b1}(T) - 2\rho_1 \int_{-1}^0 I\left(r,\mu^\prime\right) \mu^\prime \,\text{d} \mu^\prime,$$ for $0<\mu\leq 1$ and $$I(R_2, \mu) = \epsilon_2 I_{b2}(T) + 2\rho_2 \int_0^1 I\left(r,\mu^\prime\right) \mu^\prime \,\text{d} \mu^\prime,$$ for $-1\leq\mu<0$, where $\epsilon_1$ and $\epsilon_2$ are emissivities of the inner and outer surfaces, respectively. In the same way, $\rho_1$ and $\rho_2$ are the diffusive reflectivities for the inner and outer surfaces, respectively. $I_{b1}(T)$ and $I_{b2}(T)$ are the black body radiations for inner and outer surfaces in temperature $T$, respectively. In this paper, we consider standard units ($I$, $I_b$, $I_{b1}$ and $I_{b2}$ are in in W $^{-2}$ sr$^{-1}$, $p$ is in sr$^{-1}$, and the other parameters are adimensional) and, for comparison purposes, we show only numerical values of the parameters.

To obtain a solution, we use the discrete ordinates diamond difference technique in the angular variable $I_m(r)= I(r, \mu_m)$ (Chandrasekhar, 1950). This technique is based on the angular variable discretization in an enumerable set of angles or equivalently their direction cosines, and it is often called $\mathrm{S}_M$ for $M$ angles. Here, the discretized direction cosines are the discrete ordinates $\mu_m$ for $m=1,2,\dots,M$, and $w_m$ are the weights in a quadrature rule for integrals over $\left[-1,1\right]$. Using the diamond difference algorithm (Lewis & Miller, 1984), the discretized derivative term with term with respect to $\mu$ writes as $$\left[ \frac\partial{\partial \mu} \left[ \left(1-\mu^2\right) I\left(r,\mu\right) \right] \right]_{\mu=\mu_m} = \frac{ \alpha_{m+1/2}I_{m+1/2}-\alpha_{m-1/2}I_{m-1/2} }{w_m} \,,$$ where $\alpha_{m\pm1/2}$ are the angular differencing coefficients, obtained by the recursion formulas $$\begin{aligned} \alpha_{1/2} & = 0 \,, \label{5} \\ \alpha_{m+1/2} & = \alpha_{m-1/2} - 2 \mu_mw_m \,, \label{6} \\ \label{7} I_{m+1/2} &= 2I_m - I_{m-1/2}\,, \end{aligned}$$ for $m=1,2,\dots,M$. This generates a recursive set of differential equations, and to start the recursion, $I_{1/2}$ is required. Here, we define $I_{1/2}\left(r\right) = I\left(r,-1\right)$. Expanding both derivatives of equation (1) and evaluating the whole equation at $\mu=-1$ we get $$- \frac {\text{d}}{\text{d} r} I_{1/2} + I_{1/2} = \left( \vphantom{Q_Q^Q} 1-\omega\left(r\right)\right) I_b\left(T\right) + \frac{\omega\left(r\right)}{2} \sum_{m^\prime=1}^M w_{m^\prime} p\left(-1,\mu_{m^\prime}\right) I_{m^\prime} \,, \label{8}$$ whose solution starts the recursion in this diamond difference scheme. Then, taking equations (1), (3) and (4), expanding derivatives with respect to $r$, using the diamond difference approximation in the derivative terms with respect to $\mu$, evaluating the equations in $\mu=\mu_m$ and using the recursion formula in equation (7) we get

$$\begin{gathered} \label{9} \mu_m \frac{\text{d}}{\text{d}r} I_m + \frac{2\mu_m}{r} I_m + \frac{2\alpha_{m+1/2}}{r w_m} I_m - \frac{\alpha_{m+1/2} + \alpha_{m-1/2}}{r w_m} I_{m-1/2} + I_m \\ = \left(1 - \omega(r)\right) I_b(T) + \frac{\omega(r)}{2} \sum_{m^\prime=1}^M w_{m^\prime} p\left(\mu_m, \mu_{m^\prime}\right) I_{m^\prime} \end{gathered}$$ for $m=1,2,\dots,M$, $$\begin{aligned} I_m\left(R_1\right) = \epsilon_1 I_{b1}(T) - 2 \rho_1 \sum_{m^\prime=1}^{M/2} w_{m^\prime} \mu_{m^\prime} I_{m^\prime}\left(R_1\right) \end{aligned}$$ for $m=M/2+1,M/2+2,\dots,M$, and $$\begin{aligned} I_m\left(R_2\right) = \epsilon_2 I_{b2}(T) + 2 \rho_2 \sum_{m^\prime=M/2+1}^{M} w_{m^\prime} \mu_{m^\prime} I_{m^\prime}\left(R_2\right) \end{aligned}$$ for $m=1,2,\dots,M/2$. Also, only even values are valid choices for $M$, so there is no integer $m$ such that $\mu_m=0$. Furthermore, without loss of generality, we are are choosing $\mu_m$ and $w_m$ to be the abscissas and weights of a $M$-th order Gauss-Legendre quadrature rule.

To solve this set of equations from equation (9) to equation (8) we use the decomposition method (Adomian, 1988). To this end, first we write (9) as $$\begin{gathered} \label{12} \frac{\text{d}}{ \text{d} r} I_m + \left( \frac2r + \frac{2\alpha_{m+1/2}} {r \mu_m w_m} + \frac{1}{\mu_m} \right) I_m = \frac{\left( \vphantom{Q_Q^Q} 1-\omega\left(r\right)\right) I_b\left(T\right)}{\mu_m} \\ + \frac{\omega\left(r\right)}{2\mu_m} \sum_{m^\prime=1}^M w_{m^\prime} p\left(\mu_m,\mu_{m^\prime}\right) I_{m^\prime} + \frac{\alpha_{m+1/2}+\alpha_{m-1/2}} {r \mu_m w_m} I_{m-1/2} \,. \end{gathered}$$ Writing equation (12) in implicit form we get $$\begin{aligned} & \frac {\text{d}}{\text{d} r}I_m + \left(\frac2r + \frac{2\alpha_{m+1/2}} {r \mu_m w_m} + \frac{1}{\mu_m} \right) I_m = \Psi_m\left(r\right)\,, \end{aligned}$$ where $$\Psi_m\left(r\right) = \frac{\left( \vphantom{Q_Q^Q} 1-\omega\left(r\right)\right) I_b\left(T\right)}{\mu_m} + \frac{\omega\left(r\right)}{2\mu_m} \sum_{m^\prime=1}^M w_{m^\prime} p\left(\mu_m,\mu_{m^\prime}\right) I_{m^\prime} + \frac{\alpha_{m+1/2}+\alpha_{m-1/2}} {r \mu_m w_m} I_{m-1/2}\,,$$ whose solution is $$\begin{gathered} \label{13} I_m\left(r\right) = \exp\left( - \int_{c}^r \left(\frac2\tau + \frac{2\alpha_{m+1/2}} {\tau \mu_m w_m} + \frac{1}{\mu_m} \right) \text{d} \tau \right) \\ \times \left[ \vphantom{\int_{c}^r} I_m\left(c\right) + \int_{c}^r \exp\left( \int_{c}^\tau \left(\frac2\eta + \frac{2\alpha_{m+1/2}} {\eta \mu_m w_m} + \frac{1}{\mu_m} \right) \text{d} \eta \right) \Psi_m\left(\tau\right) \text{d} \tau \right] \end{gathered}$$ where $c$ is an arbitrary point in $\left[R_1,R_2\right]$.

To compute the last integral of equation (13) over $\tau$ a recursive scheme a quadrature rule is necessary. The chosen rule applies over some key values of $I_m\left(r\right)$ at the abscissas $r_\imath$, that are the assigned values in a running code. We conveniently chose the trapezoidal rule, as it does not require any kind of interpolation for the evaluation of the integral part. Thus, we discretize the interval $\left[R_1,R_2\right]$ in $N$ intervals $$\begin{aligned} \Delta r & = \frac{R_2- R_1}N \,, \label{14} \\ r_\imath & = R_1 + \imath \Delta r \,, \label{15} \end{aligned}$$ for $\imath=0,1,2,\dots,N$. Note that $r_0=R_1$ and $r_N=R_2$, and this discretization is made only for evaluation of the integral.

Depending on the sign of $\mu_m$, the $I_m^\imath$ are sequentially updated for crescent $\imath$ ($\mu>0$) or decrescent $\imath$ ($\mu<0$), but always for crescent $m$, according to the diamond difference scheme. Considering this, the $r$ discretization and the trapezoidal rule for the last integral over $\tau$, we write equations (10), (11) and (13) in the final form for the diamond difference scheme ($I_m^\imath= I_m\left(r_\imath\right)$ and $\Psi_m^\imath= \Psi_m\left(r_\imath\right)$), $$\begin{aligned} I_m^N &= \epsilon_2 I_{b2}(T) + 2\rho_2 \sum_{m^\prime=M/2+1}^M w_{m^\prime} \mu_{m^\prime} I_{m^\prime}^N \\ \label{17} I_m^\imath & = B_m^\imath \left( \vphantom\sum I_m^{\imath+1} - \Psi_m^{\imath}\frac{\Delta r}2 \right) - \Psi_m^{\imath+1} \frac{\Delta r}2 \end{aligned}$$ where $$\begin{aligned} B_m^\imath = \exp\left( \int_{r_\imath}^{r_{\imath+1}} \left(\frac2r + \frac{2\alpha_{m+1/2}} {r \mu_m w_m} + \frac{1}{\mu_m} \right) \text{d} r \right) \end{aligned}$$ for $m=1,2,\dots,M/2$ ($\mu<0$) and $\imath=N-1,N-2,\dots,0$ (decreasing $\imath$); and $$\begin{aligned} I_m^0 &= \epsilon_1 I_{b1} (T) - 2\rho_1 \sum_{m^\prime=1}^{M/2} w_{m^\prime} \mu_{m^\prime} I_{m^\prime}^0 \\ \label{19} I_m^\imath &= B_m^\imath \left( \vphantom\sum I_m^{\imath} + \Psi_m^{\imath}\frac{\Delta r}2 \right) + \Psi_m^{\imath-1} \frac{\Delta r}2 \end{aligned}$$ where $$\begin{aligned} B_m^\imath = \exp\left( - \int_{r_{\imath-1}}^{r_\imath} \left(\frac2r + \frac{2\alpha_{m+1/2}} {r \mu_m w_m} + \frac{1}{\mu_m} \right) \text{d} r \right) \end{aligned}$$ for $m=M/2+1,M/2+2,\dots,M$ ($\mu>0$) and $\imath=1,2,\dots,N$ (increasing $\imath$). Both in equations (17) and (19) we compute $\Psi_m^\imath$ as $$\Psi_m^\imath = \frac{\left( \vphantom{Q_Q^Q} 1-\omega\left(r_\imath\right)\right) I_b\left(T\right)}{\mu_m} + \frac{\omega\left(r_\imath\right)}{2\mu_m} \sum_{m^\prime=1}^M w_{m^\prime} p\left(\mu_m,\mu_{m^\prime}\right) I_{m^\prime}^\imath + \frac{\alpha_{m+1/2}+\alpha_{m-1/2}} {r_\imath \mu_m w_m} I_{m-1/2}^\imath$$ for $\imath=0,1,\dots,N$ and $m=1,2,\dots,M$.

To compute the last term of equation (21) we use equation (7), as it is valid for all $r$ (hence, for all $r_\imath$), $$\begin{aligned} I_{m+1/2}^\imath = 2I_m^\imath - I_{m-1/2}^\imath \,. \end{aligned}$$

Using the same procedures as in equation (17), the discretized equations of $I_{1/2}^\imath$ (equations (4) and (8) with $\mu=-1$) are $$\begin{aligned} I_{1/2}^N &= \epsilon_2 I_{b2}(T) + 2\rho_2 \sum_{m^\prime=M/2+1}^M w_{m^\prime} \mu_{m^\prime} I_{m^\prime}^N \,, \\[0.15cm] I_{1/2}^\imath &= B_{1/2}^\imath \left( \vphantom\sum I_{1/2}^{\imath+1} - \Psi_{1/2}^{\imath}\frac{\Delta r}2 \right) - \Psi_{1/2}^{\imath+1} \frac{\Delta r}2 \,, \end{aligned}$$

where $$\begin{aligned} B_{1/2}^\imath &= \exp\left( - \int_{r_\imath}^{r_{\imath+1}} \text{d} r \right) \,, \\[0.15cm] \label{23} \Psi_{1/2}^\imath &= - \left( \vphantom{Q_Q^Q} 1-\omega\left(r_\imath\right)\right) I_b\left(T\right) - \frac{\omega\left(r_\imath\right)}2\sum_{m^\prime=1}^M w_{m^\prime} p\left(-1,\mu_{m^\prime}\right) I_{m^\prime}^\imath \end{aligned}$$ for $\imath=N-1,N-2,\dots,0$ (decreasing $\imath$, as $\mu=-1<0$).

Although it is possible to solve this linear algebraic system using an iterative method, we chose to use a decomposition method based on Adomian, 1988. It presents a standard algorithm for power-like non-linearities, which is usual in radiative transfer modeling.

We did not consider non-linearities of any kind in this paper. However we are preparing to tackle these non-linear problems in the future, hence the choice of a decomposition method.

This method consists in expanding the unknowns in infinite series, $$\begin{gathered} \label{24} I_m^\imath = \sum_{\jmath=0}^\infty \left[ \vphantom\sum \mathcal{I}_m^\imath \right]_\jmath \,, \\ \label{25} \Psi_m^\imath = \sum_{\jmath=0}^\infty \left[ \vphantom\sum \Psi_m^\imath \right]_\jmath \,, \end{gathered}$$ for $m=1/2,1,2,3,\dots,M$ and $\imath = 0,1,\dots,N$ and making a recursive set of equations where the heterogeneities are considered only when evaluating $\left[ \vphantom\sum \mathcal{I}_m^\imath \right]_0$. For computational purposes, we truncate these series when $\jmath$ reaches some $J$, when a stop criterion is satisfied.

Substituting equations (24) and (25) in equations (16) to (23) and organizing the terms in a recursive set of equations, we have $$\begin{gathered} \label{26} \left[ \vphantom\sum \Psi_{1/2}^\imath \right]_0 = - \left( \vphantom{Q_Q^Q} 1-\omega\left(r_\imath\right)\right) I_b\left(T\right)\,, \\ \label{27} \left[ \vphantom\sum \Psi_{1/2}^\imath \right]_\jmath = - \frac{\omega\left(r_\imath\right)}2 \sum_{m^\prime=1}^M w_{m^\prime} p\left(-1,\mu_{m^\prime}\right) \left[ \vphantom\sum \mathcal{I}_{m^\prime}^\imath \right]_{\jmath-1} \,, \\ \label{28} \left[ \vphantom\sum \Psi_m^\imath \right]_0 = \frac{\left( \vphantom{Q_Q^Q} 1-\omega\left(r_\imath\right)\right) I_b\left(T\right)}{\mu_m} + \frac{\omega\left(r_\imath\right)}{2\mu_m} \sum_{m^\prime=1}^{m-1} w_{m^\prime} p\left(\mu_m,\mu_{m^\prime}\right) \left[ \vphantom\sum \mathcal{I}_{m^\prime}^\imath \right]_0 + \frac{\alpha_{m+1/2}+\alpha_{m-1/2}} {r_\imath \mu_m w_m} \left[ \vphantom\sum \mathcal{I}_{m-1/2}^\imath \right]_0 \,, \end{gathered}$$ $$\begin{gathered} \label{29} \left[ \vphantom\sum \Psi_m^\imath \right]_\jmath = \frac{\omega\left(r_\imath\right)}{2\mu_m} \left( \sum_{m^\prime=1}^{m-1} w_{m^\prime} p\left(\mu_m,\mu_{m^\prime}\right) \left[ \vphantom\sum \mathcal{I}_{m^\prime}^\imath \right]_{\jmath} \right. \\ \left. + \sum_{m^\prime=m}^{M} w_{m^\prime} p\left(\mu_m,\mu_{m^\prime}\right) \left[ \vphantom\sum \mathcal{I}_{m^\prime}^\imath \right]_{\jmath-1} \right) + \frac{\alpha_{m+1/2}+\alpha_{m-1/2}} {r_\imath \mu_m w_m} \left[ \vphantom\sum \mathcal{I}_{m-1/2}^\imath \right]_\jmath \,, \end{gathered}$$ $$\begin{gathered} \label{30} \left[ \vphantom\sum \mathcal{I}_m^{N} \right]_0 = \epsilon_2 I_{b2}(T) \,,\\ \label{31} \left[ \vphantom\sum \mathcal{I}_m^{N} \right]_\jmath = 2 \rho_2 \sum_{m^\prime=M/2+1}^M w_{m^\prime} \mu_{m^\prime} \left[ \vphantom\sum \mathcal{I}_{m^\prime}^{N} \right]_{\jmath-1} \,, \\ \label{32} \left[ \vphantom\sum \mathcal{I}_m^0 \right]_0 = \epsilon_1 I_{b1}(T) - 2 \rho_1 \sum_{m^\prime=1}^{M/2} w_{m^\prime} \mu_{m^\prime} \left[ \vphantom\sum \mathcal{I}_{m^\prime}^0 \right]_0 \,,\\ \label{33} \left[ \vphantom\sum \mathcal{I}_m^0 \right]_\jmath = - 2 \rho_1 \sum_{m^\prime=1}^{M/2} w_{m^\prime} \mu_{m^\prime} \left[ \vphantom\sum \mathcal{I}_{m^\prime}^0 \right]_\jmath \,, \end{gathered}$$ $$\begin{gathered} \label{34} \left[ \vphantom\sum \mathcal{I}_m^\imath \right]_\jmath = B_m^\imath \left( \left[ \vphantom\sum \mathcal{I}_m^{\imath+1} \right]_\jmath - \left[ \vphantom\sum \Psi_m^{\imath} \right]_\jmath \frac{\Delta r}2 \right) - \left[ \vphantom\sum \Psi_m^{\imath+1} \right]_\jmath \frac{\Delta r}2 \,, \\ \label{35} \left[ \vphantom\sum \mathcal{I}_m^\imath \right]_\jmath = B_m^\imath \left( \left[ \vphantom\sum \mathcal{I}_m^{\imath-1} \right]_\jmath + \left[ \vphantom\sum \Psi_m^{\imath} \right]_\jmath \frac{\Delta r}2 \right) + \left[ \vphantom\sum \Psi_m^{\imath-1} \right]_\jmath \frac{\Delta r}2 \,, \\ \label{36} \left[ \vphantom\sum \mathcal{I}_{m+1/2}^\imath \right]_\jmath = 2 \left[ \vphantom\sum \mathcal{I}_{m}^\imath \right]_\jmath - \left[ \vphantom\sum \mathcal{I}_{m-1/2}^\imath \right]_\jmath \,. \end{gathered}$$ Equations (26) to (36) are valid for different ranges of $\imath$, $m$ and $\jmath$. Those ranges can be checked in Table 1.

The steps of this recursive solver are listed below.

1. Input problem and numerical data: $\omega\left(r\right)$, $I_b\left(r\right)$, $I_{b1}\left(r\right)$, $I_{b2}\left(r\right)$, $\beta_\ell$ for $\ell=0,1,\dots,L$, $R_1$, $R_2$, $\epsilon_1$, $\epsilon_2$, $\rho_1$, $\rho_2$, $M$ and $N$.

2. Pre-processing:

   1. Compute $p\left(\mu,\mu'\right)$ as in equation (2); $\Delta r$ as in equation (14) and then $r_\imath$ for $i=0,1,\dots,N$ as in equation (15); $\mu_m$ and $w_m$ for $m=1,2,\dots,M$ using the $M$-th order Gauss-Legendre quadrature.

   2. Compute $\alpha_{m\pm1/2}$ as in equations (5) and (6) for $m=1,2,\dots,M$.

   3. (Recommended, but unnecessary) Compute and have variables assigned to the coefficients $\omega\left(r_\imath\right) w_m p\left(-1,\mu_{m'}\right)/2$, $\omega\left(r_\imath\right) w_m p\left(\mu_m,\mu_{m'}\right)/2\mu_m$, $\left(\alpha_{m+1/2}+\alpha_{m-1/2}\right)/r_\imath \mu_m w_m$, $\mu_m w_m$ and $\Delta r/2$ for $\imath=0,1,\dots,N$ and $m,m'=1,2,\dots,M$.

   4. Compute $B_m^\imath$ for $\imath=0,1,\dots,N$ and $m=1/2,1,2,3,\dots,M$ as in equations (18), (20) and (22). Note that the $B_m^\imath$ are computed analitically.

3. First terms. For $j=0$,

   1. Compute $\left[\Psi_{1/2}^\imath\right]_0$ as in equation (26) for $\imath=0,1,\dots,N$.

   2. Compute $\left[\mathcal I_{1/2}^N\right]_0$ as in equation (30) for $m=1/2$.

   3. Sequentially, for $\imath=N-1,N-2,\dots,0$, compute $\left[\mathcal I_{1/2}^\imath\right]_0$ as in equation (34) for $m=1/2$.

   4. Sequentially, for $m=1,2,\dots,M/2$,

      1. Compute $\left[\Psi_{m\vphantom/}^\imath\right]_0$ as in equation (28) for $\imath=0,1,\dots,N$.

      2. Compute $\left[\mathcal I_{m\vphantom/}^N\right]_0$ as in equation (30).

      3. Sequentially, for $\imath=N-1,N-2,\dots,0$, compute $\left[\mathcal I_{m\vphantom/}^\imath\right]_0$ as in equation (34).

      4. Compute $\left[\mathcal I_{m+1/2}^\imath\right]_0$ as in equation (36) for $\imath=0,1,\dots,N$.

   5. Sequentially, for $m=M/2+1,M/2+2,\dots,M$,

      1. Compute $\left[\Psi_{m\vphantom/}^\imath\right]_0$ as in equation (28) for $\imath=0,1,\dots,N$.

      2. Compute $\left[\mathcal I_{m\vphantom/}^0\right]_0$ as in equation (32).

      3. Sequentially, for $\imath=1,2,\dots,N$, compute $\left[\mathcal I_{m\vphantom/}^\imath\right]_0$ as in equation (35).

      4. Compute $\left[\mathcal I_{m+1/2}^\imath\right]_0$ as in equation (36) for $\imath=0,1,\dots,N$.

4. Other recursions. For $\jmath=1,2,\dots,J$,

   1. Compute $\left[\Psi_{1/2}^\imath\right]_\jmath$ as in equation (27) for $\imath=0,1,\dots,N$.

   2. Compute $\left[\mathcal I_{1/2}^N\right]_\jmath$ as in equation (31) for $m=1/2$.

   3. Sequentially, for $\imath=N-1,N-2,\dots,0$, compute $\left[\mathcal I_{1/2}^\imath\right]_\jmath$ as in equation (34) for $m=1/2$.

   4. Sequentially, for $m=1,2,\dots,M/2$,

      1. Compute $\left[\Psi_{m\vphantom/}^\imath\right]_\jmath$ as in equation (29) for $\imath=0,1,\dots,N$.

      2. Compute $\left[\mathcal I_{m\vphantom/}^N\right]_\jmath$ as in equation (31).

      3. Sequentially, for $\imath=N-1,N-2,\dots,0$, compute $\left[\mathcal I_{m\vphantom/}^\imath\right]_\jmath$ as in equation (34).

      4. Compute $\left[\mathcal I_{m+1/2}^\imath\right]_\jmath$ as in equation (36) for $\imath=0,1,\dots,N$.

   5. Sequentially, for $m=M/2+1,M/2+2,\dots,M$,

      1. Compute $\left[\Psi_{m\vphantom/}^\imath\right]_\jmath$ as in equation (29) for $\imath=0,1,\dots,N$.

      2. Compute $\left[\mathcal I_{m\vphantom/}^0\right]_\jmath$ as in equation (33).

      3. Sequentially, for $\imath=1,2,\dots,N$, compute $\left[\mathcal I_{m\vphantom/}^\imath\right]_0$ as in equation (35).

      4. Compute $\left[\mathcal I_{m+1/2}^\imath\right]_0$ as in equation (36) for $\imath=0,1,\dots,N$.

5. Compute an approximation of $I_{m}^\imath$ as in equation (24), adding up to $\jmath = J$ for $m=1,2,\dots,M$ and $\imath=0,1,\dots,N$.

It is relatively simple to show the consistency of the decomposition by making the residual term go to zero as the number of terms in the sum of equation (24) increases. Upon substituting the decomposition in equation (24) truncated in the $\left(J+1\right)$-th term ($\jmath=J$) in equations (9), (10) and (11) evaluated at $\left(r,\mu\right)=\left(r_\imath,\mu_m\right)$ for all $\imath=0,1,\dots,N$ and $m=1,2,\dots,M$ we have an approximation, as the computation of each term in the series follows equations (26) to (36) make many of them cancel each other out, except for a remaining quantity we are calling the residual term. The residual terms of every node $\left(r_\imath,\mu_m\right)$, denoted as $\left[ \varepsilon_m^\imath \right]_J$, are computed as $$\begin{gathered} \label{37} \left[ \varepsilon_m^{N} \right]_J = - 2 \rho_2 \sum_{m^\prime=M/2+1}^M w_{m^\prime} \mu_{m^\prime} \left[ \mathcal{I}_{m^\prime}^{N} \right]_J \,, \\ \label{38} % i=0,1,\dots,N-1; m=1,2,\dots,M/2 \left[ \varepsilon_m^\imath \right]_J = - \frac{\Delta r}{2} \left( B_m^\imath \sum_{m^\prime=m}^{M} w_{m^\prime} p\left(\mu_m,\mu_{m^\prime}\right) \left[ \mathcal{I}_{m^\prime}^{\imath+1} \right]_J + \sum_{m^\prime=m}^{M} w_{m^\prime} p\left(\mu_m,\mu_{m^\prime}\right) \left[ \mathcal{I}_{m^\prime}^\imath \right]_J \right) \,, \\ \label{39} \left[ \varepsilon_m^0 \right]_J = 0 \,, \\ \label{40} \left[ \varepsilon_m^\imath \right]_J = \frac{\Delta r}{2} \left( B_m^\imath \sum_{m^\prime=m}^{M} w_{m^\prime} p\left(\mu_m,\mu_{m^\prime}\right) \left[ \mathcal{I}_{m^\prime}^{\imath-1} \right]_J + \sum_{m^\prime=m}^{M} w_{m^\prime} p\left(\mu_m,\mu_{m^\prime}\right) \left[ \mathcal{I}_{m^\prime}^\imath \right]_J \right) \,. \end{gathered}$$ Here, equation (37) is valid for $m=1,2,\dots,M/2$, equation (38) is valid for $i=0,1,\dots,N-1$ and $m=1,2,\dots,M/2$, equation (39) is valid for $m=M/2+1,M/2+2,\dots,M$ and equation (40) is valid for $i=1,2,\dots,N$ and $m=M/2+1,M/2+2,\dots,M$. This set of equations can be written in matrix form, $$\varepsilon_J = C \mathbf{\mathcal I}_J \,,$$ where $\varepsilon_J$ is the vector of $\left[ \varepsilon_m^\imath \right]_J$, $\mathbf{\mathcal I}_J$ is the vector of $\left[ \mathcal I_m^\imath \right]_J$ and $C$ is the associate matrix. Taking the maximum norm of equation (41) we get $\left\| \varepsilon_J \right\|_\infty = \left\| C \mathbf{\mathcal I}_J \right\|_\infty \leq \left\| C\right\|_\infty \left\|\mathbf{\mathcal I}_J\right\|_\infty \,$.

As $C$ does not vary with $J$, one may infer that $\left\| \varepsilon_J \right\|_\infty$ is majored by a constant scale of $\left\|\mathbf{\mathcal I}_J\right\|_\infty$. In other words, if $\left\|\mathbf{\mathcal I}_J\right\|_\infty\to0$, then $\left\| \varepsilon_J \right\|_\infty \to0$, hence the recursive system is consistent, given that the sum in equation (24) is convergent, which is out of the scope of this paper.

At this point we discuss with a little more detail the order of convergence of the proposed methodology. We focus on the use of the analytical solution with numerical integration instead of a standard finite difference method and we show a gain in the order of convergence, at least locally. To avoid extensive and repetitive mathematical formulation, let $a_m\left(r\right) = \dfrac2r + \dfrac{2\alpha_{m+1/2}}{r\mu_mw_m} + \dfrac1{\mu_m} \,$,

so equation (13) is written in a simpler version, $$\frac{\text{d}I_m}{\text{d}r} + a_m\left(r\right)I_m = \Psi_m\left(r\right)\,.$$

Its analytical solution is $$I_m\left(r\right) = \exp\left(-\int_{c}^r a_m\left(\tau\right)\text{d}\tau\right) \left[ I_m\left(c\right) + \int_{c}^r \exp\left(\int_{c}^\eta a_m\left(\tau\right)\text{d}\tau\right) \Psi_m\left(\eta\right)\text{d}\eta \right] \,.$$

Using the trapezoidal rule with the order of convergence term and setting $r=r_\imath$ and $c=r_{\imath+1}$ and $m\leq M/2$ yields

$$\begin{gathered} I_m\left(r_\imath\right) = \exp\left(\int_{r_\imath}^{r_{\imath+1}} a_m\left(\tau\right)\text{d}\tau\right) \left[ \vphantom{\int_{r_0}^r} I_m\left(r_0\right) \right. \\ \left. - \frac{\Delta r}2 \left( \exp\left(-\int_{r_\imath}^{r_{\imath+1}} a_m\left(\tau\right)\text{d}\tau\right) \Psi_m\left(r\right) + \Psi_m\left(r_0\right) \right) + O\left(\Delta r^3\right) \right] \,, \end{gathered}$$ so the order of convergence term in computing $I_m\left(r_\imath\right)$ would be $\exp\left( \int_{r_\imath}^{r_{\imath+1}} a_m\left(\tau\right)\text{d}\tau \right) O\left(\Delta r^3\right)\,$.

It can be shown that the integral is always negative, so the exponential is always less than one, but as we make $\Delta r\to0$, we conclude that the order of convergence of the present methodology gets closer to $O\left(\Delta r^3\right)$. Applying a finite difference scheme in equation (42) and carrying out the order of convergence term $O\left(\Delta r\right)$ will result in $$I_m\left(r_\imath\right) = \frac{I_m\left(r_{\imath+1}\right)-\Psi_m\left(r_\imath\right)\Delta r} {1-a_m\left(r_\imath\right)\Delta r} + \frac{O\left(\Delta r^2\right)} {1-a_m\left(r_\imath\right)\Delta r} \,.$$ Making $\Delta r\to0$, we conclude that the order of convergence of this finite difference scheme gets closer to $O\left(\Delta r^2\right)$, hence the presented methodology shows an improvement in this aspect. We have shown this approach for $m\leq M/2$, however the same conclusion is obtained for $m>M/2$.

Results and discussion

The presented methodology was implemented in a Python script, ran an a domestic computer with some literature inputs, separated into two cases. For comparison purposes, we define the backward and forward radiation fluxes as $$\begin{aligned} q^-\left(r_i\right) &= \int_{-1}^0 I\left(r,\mu\right) \mu \,\text{d} \mu = \sum_{m^\prime=1}^{M/2} w_{m^\prime} \mu_{m^\prime} I_{m^\prime}^i \,, \\ q^+\left(r_i\right) &= \int_0^1 I\left(r,\mu\right) \mu \,\text{d} \mu = \sum_{m^\prime=M/2+1}^M w_{m^\prime} \mu_{m^\prime} I_{m^\prime}^i \,, \end{aligned}$$ where the sums are the integral approximations using the $M$-th order quadrature rule. In both cases, the coefficients $\beta_{\ell}$ in equation (2) are evaluated according to Table 2 (Ladeia et al., 2020; Abulwafa, 1993). Besides, if $\beta_1$ is positive or negative, we say we have a forward or backward scattering phase function, respectively. Also, if $\beta_\ell=0$ for $\ell\geq1$ we say it is an isotropic scattering phase function (Petty, 2006). Nonetheless, the authors chose, for simplicity, $N=1280$ and $M=24$ for all cases, and the stop criterion is that when the maximum value of $\left|\left(\vphantom\sum \mathcal I_m^\imath \right)_\jmath \right|$ for all $\imath=0,1,\dots,N$ and $m=1,2,\dots,M$ is less or equal than $10^{-6}$, the recursion stops.

In the following, we will demonstrate that our methodology can be applied to several sets of parameters from the literature, considering inhomogeneous source terms and anistotropic scattering in Case 1 and to sets of parameters that are commonly found in radiation transfer problems with non-linear coupling with diffusion problems in Case 2.

Case 1

In the Case 1, the results were computed and compared with Abulwafa, 1993; Ladeia et al., 2020, where considers several subcases about inhomogeneous hollow spheres, with numerical values $R_1=1$ and $R_2=2$ for the inner and outer radii. Also, the fixed parameters are $\epsilon_1 = \epsilon_2 = 0.75$, $\rho_1 = \rho_2 = 0.25$, $I_{b1}\left(T\right) = 0$, $I_{b2}\left(T\right)=4/3$.

The subcases differ in the combinations of $\left(1-\omega\left(r\right)\right) I_b\left(T\right)$, $\omega\left(r\right)$ and $p\left(\mu,\mu^\prime\right)$ given by equation (2), where the values of $\beta_\ell$ are given by Table 2 and the formulas for the different $\omega\left(r\right)$ are those in Table 3, for a total of 99 subcases.

In the description of these subcases, we use two auxiliary constants $F$ and $H$, computed as

$$\begin{aligned} F & = \frac{\left(R_2\right)^4-\left(R_1\right)^4} {\left(R_2\right)^3-\left(R_1\right)^3} \,. \end{aligned}$$

and $$\begin{aligned} \label{43} H & = \frac{\left(R_2\right)^5-\left(R_1\right)^5} {\left(R_2\right)^3-\left(R_1\right)^3} \,. \end{aligned}$$

The values of $q^+\left(R_2\right)$ and $q^-\left(R_1\right)$ for each subcase were compared with the data from the references Abulwafa, 1993 and Ladeia et al., 2020, denoted as REF1 and REF2, respectively, as well as with the results of the current methodology, referred to as PM. The results are presented in Tables 4, 5 and 6 and listed the largest values for each combination of $(1-\omega(r))I_b(T)$ and scattering, as described next.

0.15cm

0.15cm

0.15cm

To simplify the comparison of results in Case 1, we also computed the relative distances (RD) for the results in Tables 4, 5 and 6 using the equation (44)

$$\begin{aligned} \label{44} & \text{RD} = \left|\dfrac{\text{value}~{\text{reference}} - \text{value}~{\text{PM}}}{\text{value}~{\text{reference}}}\right| \,, \end{aligned}$$

and listed the largest values for each combination of $\left(1-\omega\left(r\right)\right)I_b\left(T\right)$ and scattering type in Table 7 for both references and for which $\omega\left(r\right)$ and partial flux it happened. We displayed the values of relative distances in percentage, rounded up to the third decimal place.

In other words, we show in Table 7 the largest values of the relative distance we may find compared to each reference column in Tables 4, 5 and 6. For example, the first line in Table 7 indicates that the relative distance from the presented methodology to the reference Abulwafa, 1993 for $\left(1-\omega\left(r\right)\right)I_b\left(T\right)=0.0$ and forward scattering is no larger than \(0.975\\), that we find for the outgoing flux in $R_1$ for $\omega\left(r\right)=\omega_6\left(r\right)$, from Table 3. We omitted the other relative distances we computed due to the lack of space.

We showed in Table 7 that the results show fairly good agreement to the references. The maximum relative distance among all subcases is about 4.2%, and the vast majority is below 1.0%. 13 out of 18 subcases present the maximum relative distance when $\omega\left(r\right) = \omega_7\left(r\right)$ for the outgoing flux at $R_2$, thus establishing a mode. Even this mode had all relative distances below 1%.

0.25cm

$^*$$H$ from equation (43).

Case 2

In the Case 2, we consider three subcases of homogeneous hollow spheres (constant values for $\omega$).

All subcases consider $R_1=1$, $R_2=2$, $\epsilon_1=1$, $\rho_1=0$, $p\left(\mu,\mu^\prime\right)=1$, $I_b\left(T\right)=0$, $I_{b1}\left(T\right)=1$ and $I_{b2}\left(T\right)=1$. The remaining subcases are presented in Table 8.

0.15cm

These are usual sets of parameters in radiative transfer problems coupled with heat conduction. We present the outgoing fluxes as results for these subcases in Table 9.

0.15cm

This Case 2 has only 3 subcases with isotropic scattering, no external source term and constant values for $\omega$, thus being far simpler than Case 1, however it represents many usual parameters for radiation transfer cases coupled with heat conduction, where the black body radiations are modeled as proportional to the temperature to the fourth power.

As the focus of this paper is to present a consistent methodology for linear cases, we adapted those problems removing the nonlinear part, so there are no references to compare the results with. In this context the present approach is a first step in this direction.

We were unable to compare our results with another methodology due to the lack of publications with these set of parameters.

Case 1 is about testing the current methodology to some cases with anysotropy and heterogeneity and compare the results with the literature. Case 2 is about showing the preparation to the application of the decomposition method to a nonlinear case in a future work. The results demonstrate that the proposed methodology is robust and consistently applicable across all cases anlyzed in this study.

Conclusions

The radiative transfer problems in spherical geometry are usually complex and difficult to solve in general cases. In this work we solved the linear problem for a hollow sphere, including an anisotropic case using an algorithm that consists of a combination of methods in the literature. The mentioned algorithm combines a decomposition method, where the terms are computed using a diamond difference scheme, hence avoiding dealing with complex and costly algebraic operations. This combination produced a robust and efficient algorithm to solve the problems here specified.

The cases we chose to present applications of the methodology involve some usual parameters in transport problems, like anisotropy, dependence of $r$ in parameters and diffuse-reflective boundary conditions. It is noteworthy that despite the phase function as in equation (2) showing that the anisotropy may be modeled as a sum of orthogonal polynomials, this is not required by the presented methodology. In fact, our methodology does not imply any restrictions to the construction of the phase function, single scattering albedo or the black body radiation term, as shown in the previous sections.

The presented methodology has only one local approximation in the integral from the trapezoidal rule. We showed that our methodology has a superior order of convergence compared to classical finite difference method.

Ladeia et al., 2020 used a finite volumes method, which has a local order of convergence of $O\left(\Delta r^2\right)$ like the finite difference method, and Abulwafa, 1993 used a Galerkin method, whose order of convergence we were unable to determine due to the lack of dedicated information in the paper. So the presented methodology converges locally faster than the methodology in Ladeia et al., 2020.

This was possible because there were no approximations in the computations of $B_m^i$ and $\Psi_m^i$. It is also possible to increase the order of convergence to $O\left(\Delta r^4\right)$ or higher by using higher order quadrature rules like Simpson or other Newton-Cotes formulas. However the schemes would require more algebraic work.

The order of convergence for the recursive system, the discretization of the angular variable, and the whole composite trapezoidal rule (over $R_1\leq r \leq R_2$) are basically the same as in the references.

We also presented studies on the method’s consistency and showed that upon the construction of the recursive system, the method is automatically consistent, at least for the linear case. The convergence and stability analysis are left for future work, in which we also aim to improve the decomposition to consider insertion of non-linearity temperature coupling with conduction effects.

Author contributions

M. Schramm, participated in: formal analysis, methodology, editing and writing - original draft. C. A. Ladeia participated in: methodology, supervision, writing-revision and editing. J. C. L. Fernandes participated in: supervision and writing-revision.

Conflicts of interest

The authors declare no conflict of interest.

References

Abulwafa, E. M. (1993). Radiative-transfer in a linearly-anisotropic spherical medium. 165--175. https://doi.org/10.1016/0022-4073(93)90057-O

Adomian, G. (1988). A review of the decomposition method in applied mathematics. 501--544. https://doi.org/10.1016/0022-247X(88)90170-9

Allahviranloo, T. (2005). The Adomian decomposition method for fuzzy system of linear equations. 553--563. https://doi.org/10.1016/j.amc.2004.02.020

Chandrasekhar, S. (1950). Radiative Transfer. Oxford University Press.

Haq, F., Shah, K., Rahman, G. & Shahzad, M. (2018). Numerical solution of fractional order smoking model via Laplace Adomian decomposition method. 1061--1069. https://doi.org/10.1016/j.aej.2017.02.015

Howell, J. R., Menguc, M. P. & Siegel, R. (2016). Thermal Radiation Heat Transfer. CRC Press.

Ladeia, C. A., Bodmann, B. E. J. & Vilhena, M. T. (2019). On the integro-differential radiative conductive transfer equation: A modified decomposition method. Birkhauser. 197--210. https://doi.org/10.1007/978-3-030-16077-7_16

Ladeia, C. A., Schramm, M. & Fernandes, J. C. L. (2020). A simple numerical scheme to linear radiative transfer in hollow and solid spheres. 21--30. https://doi.org/10.5433/1679-0375.2020v41n1p21

Lewis, E. E. & Miller, W. F., Jr. (1984). Computational Methods of Neutron Transport. John Wiley \& Sons.

Li, R., Li, W. & Zheng, L. (2020). A nonlinear three-moment model for radiative transfer in spherical symmetry. 285--299. https://doi.org/10.1016/j.matcom.2019.11.004

Modest, M. F. (2013). Radiative Heat Transfer. Academic Press.

Ozisik, M. (1973). Radiative Transfer and Interaction with Conductions and Convection. John Wiley \& Sons Inc..

Petty, G. W. (2006). A First Course in Atmospheric Radiation. Sundog Publishing.

Segatto, C. F., Tomaschewski, F. K., Barros, R. C. & Vilhena, M. T. (2017). On the solution of the SN multigroup kinetics equations in multilayer slabs. 229--236. https://doi.org/10.1016/j.anucene.2017.02.016

Sghaier, T. (2013). Study of radiation in spherical media using moment method. 25--32. https://doi.org/10.11648/j.ajpa.20130101.15

Stamnes, K., Thomas, G. E. & Stamnes, J. J. (2017). Radiative Transfer in the Atmosphere and Ocean. Cambridge University Press.

Tsai, J. R., Osizik, M. N. & Santarelli, F. (1989). Radiation in spherical symmetry with anisotropic scattering and variable properties. 187--199. https://doi.org/10.1016/0022-4073(89)90082-4

Wazwaza, A. M. & El-Sayed, S. M. (2001). A new modification of the Adomian decomposition method for linear and nonlinear operators. 393--405. https://doi.org/10.1016/S0096-3003(00)00060-6

Xu, F., He, X., Jin, X., Cai, W., Bai, Y., Wang, D., Gong, F. & Zhu, Q. (2023). Spherical vector radiative transfer model for satellite ocean color remote sensing. 11192--11212. https://doi.org/10.1364/OE.483221