A Nonlinear Moment Model for Radiative Transfer Equation

The radiative transfer equation (RTE) depicts the motion of photons and their interaction with the background medium. It has lots of applications, such as radiation astronomy [35], optical imaging [26, 40], neutron transport in reactor physics [37, 16], light transport in atmospheric radiative transfer [31] and heat transfer [27]. Due to the integro-differential form and the high-dimensionality of RTE, how to develop efficient methods to solve it numerically is an important but challenging topic. So far, the commonly used numerical methods can be categorized into two groups: the probabilistic methods, like the direct simulation Monte Carlo (DSMC) method [21, 3, 24, 14, 1], and the deterministic methods [4, 28, 39, 25, 13, 15, 36, 2, 17, 20, 19], such as the discrete ordinates method ($S_{N}$) [4, 28, 39], the moment methods [25, 13, 15, 2, 20, 19] and etc.

The discrete ordinates method ($S_{N}$) is one of the most popular numerical methods to simulate the RTE, which solve the RTE along with a discrete set of angular directions from a given quadrature set. However, the $S_{N}$ model assumes that the particles can only move along the directions in the quadrature set, thus once the coordinate system is rotated, the results of the $S_{N}$ model can be different. The lack of rotational invariance results in numerical artifacts, known as ray effects [28].

In order to reduce the complexity of the RTE, the moment method focuses on the evolution of a finite number of moments of the specific intensity, which avoids the high-dimensionality of directly solving the RTE. Since the governing equation of a lower order moment commonly contains higher order moments, the moment system is often not automatically closed. Hence one has to take a moment closure to close the moment system. A practical method for the moment closure is to construct an ansatz to approximate the specific intensity. The pioneer works in moment method include the spherical harmonics method ($P_{N}$) [37] and the maximum entropy method ($M_{N}$) [29, 15, 36]. The $P_{N}$ model constructs the ansatz using spherical harmonic polynomials. It can be regarded as a polynomial expansion of the specific intensity around the equilibrium, which is a constant function. One of the flaws is that the resulting system may lead to nonphysical oscillations, or even worse, negative particle concentration [6, 7, 34]. The $M_{N}$ model constructs the ansatz using the principle of maximum entropy, as the maximum entropy closure for Boltzmann equation [29, 15]. Unfortunately, no explicit expression of the moment closure for $M_{N}$ model can be given when the order $N\geq 2$. To implement the model numerically, one has to solve an ill-conditioned optimization problem to obtain an approximate moment closure. This almost prohibits the application of the $M_{N}$ model.

Recently, a nonlinear moment model (called the ${M\!P}_{\!N}$ model) was proposed in [20] for the RTE in slab geometry. This model takes the ansatz of the $M_{1}$ model (the first order $M_{N}$ model) as the weight function, then constructs the ansatz by expanding the specific intensity around the weight function in terms of orthogonal polynomials in the velocity variables. Numerical examples in [20] demonstrated a quite promising performance as an improved approximation of the intensity in comparison of the $P_{N}$ model. The ${M\!P}_{\!N}$ model was further improved in [19] by a globally hyperbolic regularization following the framework developed in [8, 9, 10, 18]. We note that the regularization in [19] is a subtle modification of the work in [10] instead of a direct application. Otherwise, the resulting system may change the $M_{1}$ model, which leads to a wrong higher-order Eddington approximation. Eventually, the ${H\!M\!P}_{\!N}$ model was proposed in [19] with not only global hyperbolicity, but also a physical higher-order Eddington approximation.

Encouraged by the elegant mathematical structure and the promising numerical performance of the ${H\!M\!P}_{\!N}$ model for RTE in slab geometry, we in this paper try to extend the method to derive the ${H\!M\!P}_{\!N}$ model for 3D problems. The steps of the extension are clear while there are still numerous difficulties. Fortunately, the 3D $M_{1}$ model is explicit, which allow us to construct the ansatz to approximate the specific intensity using again the weighted polynomials with the weight function is the ansatz of the $M_{1}$ model. To construct the function space of the weighted polynomials, we need to give the orthogonal polynomial basis with respect to the weight function. For the slab geometry [20], this can be implemented by a simple Gram-Schmidt orthogonalization. For the 3D case, we have to use quasi-orthogonal polynomials rather than orthogonal polynomials. Otherwise, it can be extremely involving to accomplish the calculation, which makes further analysis to the resulted model prohibited. We propose a procedure to make a quasi Gram-Schmidt orthogonalization to obtain the quasi-orthogonal polynomials in explicit expressions. This provides us a 3D ${M\!P}_{\!N}$ model in explicit formation, which can be mathematically analyzed. To achieve global hyperbolicity, we still adopt the method in [19] to regularize the 3D ${M\!P}_{\!N}$ model. Quite smoothly a globally hyperbolic 3D ${H\!M\!P}_{\!N}$ model is eventually attained with a list of fabulous mathematical natures inherited from its 1D counterpart for the slab geometry. The resulted model is rotational invariant, with wave speeds not greater than that of light, spectral approximation accuracy, and correct higher-order Eddington approximation. We carry out preliminary numerical simulating using an abruptly splitting scheme to validate the new 3D ${H\!M\!P}_{\!N}$ model. Some numerical examples on typical problems are presented with satisfactory performance.

The rest of this paper is arranged as follows. In Section 2, we briefly introduce the moment methods for RTE, and review how the ${M\!P}_{\!N}$ and the ${H\!M\!P}_{\!N}$ model in slab geometry were derived in [20, 19]. In Section 3, we derive the 3D ${M\!P}_{\!N}$ model and prove that the model is rotational invariant. The hyperbolic regularization is applied to give the ${H\!M\!P}_{\!N}$ model in Section 4. The model is analyzed in detail therein. In Section 5, we introduce the numerical scheme to carry out numerical simulations and present some numerical examples. The paper is then ended with a short conclusion remarks.

To model radiative transfer, the governing equation is a time-dependent equation of the specific intensity $I$ as

where $c$ is the speed of light, and the specific intensity $I=I(t,\bm{x};\bm{\Omega},\nu)$ depends on time $t\in\mathbb{R}^{+}$, the spatial coordinate of the photon $\bm{x}\in\mathbb{R}^{3}$, the velocity direction $\bm{\Omega}\in\mathbb{S}^{2}$ and the frequency $\nu\in\mathbb{R}^{+}$. In this paper, our study omits the independent variable $\nu$ that $I$ is a function of $t$, $\bm{x}$ and $\bm{\Omega}$ only. The right hand side ${\mathcal{S}}(I)$ denotes the actions by the background medium on the photons. A form of ${\mathcal{S}}(I)$ adopted commonly was given in [5, 33] as

where $a$ is the radiation constant, and $s=s(t,\bm{x})$ is an isotropic external source of radiation. The scattering coefficient $\sigma_{s}$, the absorption coefficient $\sigma_{a}$, and the material temperature $T(t,\bm{x})$ depend on time $t$ and the spatial position $\bm{x}$. The total opacity coefficient is $\sigma_{t}=\sigma_{a}+\sigma_{s}$.

In case that the problems slab geometry and spherical symmetric geometry are considered, the 3D RTE (2.1) can be simplified to 1D problem. Precisely, in the slab geometry, the specific intensity depends only upon the single spatial coordinate $z$ and the single angular coordinate $\arccos\mu$, the angle between $\bm{\Omega}$ and the $z$-axis. Then the specific intensity becomes $I=I(z,\mu)$, and (2.1) is simplified as

The spherical geometry with perfect symmetry is a slightly more complicated case. The RTE, where the specific intensity depends upon only on the distance from the origin $r=\|\bm{x}\|$, and the angular variable $\arccos\mu$, which is the angle between $\bm{\Omega}$ and $\bm{x}$. In this case, $I=I(r,\mu)$, and the 3D RTE (2.1) is simplified as

In [2, 20, 19, 30], for slab geometry and spherical symmetric geometry, some nonlinear moment models had been derived with global hyperbolicity and promising performance in handling problems with fair extreme specific intensity functions. The major aim of this paper is to develop models for 3D problems with similar techniques.

At first, we define the moments of the 3D specific intensity. Let $\alpha\in\mathbb{N}^{3}$ be a 3D multi-index, i.e. $\alpha=(\alpha_{1},\alpha_{2},\alpha_{3})^{T}$, $\alpha_{1},\alpha_{2},\alpha_{3}\in\mathbb{N}$. We define a function of $t$, and $\bm{x}$, denoted by ${\langle I\rangle}_{\alpha}(t,\bm{x})$, as

where $\bm{\Omega}^{\alpha}=\Omega_{1}^{\alpha_{1}}\Omega_{2}^{\alpha_{2}}\Omega_{3}^{\alpha_{3}}$. We call that ${\langle I\rangle}_{\alpha}$ is the $\alpha$-th moment of the specific intensity $I$.

Notice that $\bm{\Omega}\in\mathbb{S}^{2}$ implies $\|\bm{\Omega}\|=1$, thus one has

where $e_{d}$ represents the multi-index whose $d$-th index is 1, and the else two indexes are 0. Therefore, we only need to consider these moments ${\langle I\rangle}_{\alpha}$, $\alpha\in{\mathcal{I}}$, where ${\mathcal{I}}$ is a set of multi-indexes, defined by

The order of the multi-index $\alpha$ is defined as $|\alpha|=\sum_{d=1}^{3}\alpha_{d}$, and we denote that ${\mathcal{I}}_{N}\triangleq\left\{\alpha:\alpha\in{\mathcal{I}},|\alpha|\leq N\right\}$. It is clear that once $\left\{{\langle I\rangle}_{\alpha}:\alpha\in{\mathcal{I}}_{N}\right\}$ is determined, one can obtain $\left\{{\langle I\rangle}_{\alpha}:\alpha\in\mathbb{N}^{3},|\alpha|\leq N\right\}$. This allows us to discuss ${\langle I\rangle}_{\alpha}$ for $\alpha\in{\mathcal{I}}$ only to derive reduced models.

Multiplying (2.1) by $\bm{\Omega}^{\alpha}$, and taking the integration with respect to $\bm{\Omega}$ over $\mathbb{S}^{2}$, one can have

In order to derive a moment model for (2.1), we first truncate the system by discarding all the governing equations of high order moments ${\langle I\rangle}_{\alpha}$, where $|\alpha|>N$, for a given integer $N\in\mathbb{N}$. The truncated moment system is

However, the governing equations of ${\langle I\rangle}_{\alpha}$, $|\alpha|=N$, involve three $N+1$ order moments ${\langle I\rangle}_{\alpha+e_{d}}$ for $d=1,2,3$, thus the truncated system (2.8) is not closed. Therefore, we need to determine all these moments, $\left\{{\langle I\rangle}_{\alpha},\alpha\in{\mathcal{I}}_{N+1}\right\}$ to make this truncated system (2.8) closed.

We divide the moments in $\left\{{\langle I\rangle}_{\alpha},\alpha\in{\mathcal{I}}_{N+1}\right\}$ into two parts: lower-order moments (also referred as known moments later on), and higher-order moments (also referred as unknown moments later on).

The aim of the so-called moment closure to this system is to approximate the unknown moments as functions of known moments, saying to give a formulation as

To achieve this goal, a practical approach is to construct an ansatz for the specific intensity. Precisely, let $E_{\alpha}$, $\alpha\in{\mathcal{I}}_{N}$, be the known moments for a certain unknown specific intensity $I$. Then one may propose an expression $\hat{I}(\bm{\Omega};E_{\alpha},\alpha\in{\mathcal{I}}_{N})$, called an ansatz to approximate $I$, such that

Often we require that $\hat{I}$ is uniquely determined by the consistency relations (2.11). With $\hat{I}$ given, the higher-order moments of $I$ are then approximated by the higher-order moments of $\hat{I}$, i.e.,

Therefore, one may take the closed moment system

as the reduced model to approximate the original RTE, where $E_{\alpha+e_{d}}$ are functions of $E_{\alpha}$, $\alpha\in{\mathcal{I}}_{N}$, defined in (2.10) and (2.12).

Many existing models can be regarded as consequences using this moment closure approach. For example, the $P_{N}$ model [25], the $M_{N}$ model [29, 15], the positive $P_{N}$ model [22], the $B_{2}$ model [2], and the ${M\!P}_{\!N}$ model [20] are in this fold. The ${M\!P}_{\!N}$ model we proposed in [20], and then improved in [19] as the ${H\!M\!P}_{\!N}$ model, is limited for problems in slab geometry, where it exhibits satisfactory numerical performance for some standard benchmarks. To extend the method therein to 3D RTE, below we first briefly review the methods in [20, 19] to derive models in slab geometry to clarify our idea.

The ${M\!P}_{\!N}$ model derived in [20] for RTE in slab geometry is based on a method to combine the $P_{N}$ model and the $M_{N}$ model, which was implemented by expanding the specific intensity around the ansatz of the $M_{1}$ model in terms of orthogonal polynomials. The ansatz of the $M_{1}$ model in slab geometry is

where $\varepsilon$ and $c_{0}$ are determined by the $0$-th moment $E_{0}$ and $1$-th moment $E_{1}$, formulated as

Then the specific intensity is approximated by a weighted polynomial, with the weight function

The function space of the weighted polynomials is

Then the ansatz of the ${M\!P}_{\!N}$ model is written as

where $\Phi^{[c_{0}]}_{i}(\mu)=\phi^{[c_{0}]}_{i}(\mu)\omega^{[c_{0}]}(\mu)$, $i=0,1,\cdots,N$, are the basis functions, $\phi^{[c_{0}]}_{i}(\mu)$ are orthogonal polynomials with respect to the weight function, and $f_{i}$ are the expansion coefficients.

The orthogonal polynomials $\phi^{[c_{0}]}_{k}(\mu)$ can be calculated by a simple Gram-Schmidt orthogonalization, formulated as

where ${\mathcal{K}}_{j,k}=\int_{-1}^{1}\mu^{j}\phi^{[c_{0}]}_{k}(\mu)\omega^{[c_{0}]}(\mu)\,\mathrm{d}\mu$, calculated by

Furthermore, by

one can determine the coefficients $f_{k}$, and the ansatz $\hat{I}(\mu;E_{0},E_{1},\cdots E_{N})$. Finally, we have the moment closure as

Meanwhile, in the viewpoint of orthogonal projection, if we define the orthogonal projection to function space $\mathbb{H}^{[c_{0}]}_{N}$,

then the ${M\!P}_{\!N}$ moment system can be written as

The weight function (2.16) permits the ${M\!P}_{\!N}$ model to approximate a strongly anisotropic distribution with very high accuracy. In [19], the ${M\!P}_{\!N}$ model was further improved by a hyperbolic regularization, which provides global hyperbolicity. To achieve the required hyperbolicity, the model reduction framework in [10, 18] suggests adding one more projection between the operators $\mu\cdot$ and $\dfrac{\partial{\cdot}}{\partial{z}}$ in (2.24) to regularize the ${M\!P}_{\!N}$ model to be globally hyperbolic, and the resulting model is

An interesting point observed in [19] is that this regularization changes the ${M\!P}_{\!N}$ model when $N=1$. It is definitely inappropriate that the $M_{1}$ model is changed by the regularization. In order to fix this defect, another weight function is introduced as

and then a new function space is defined as

In this subspace, one can define the orthogonal polynomials $\tilde{\phi}^{[c_{0}]}_{k}(\mu)$ and the basis function $\tilde{\Phi}^{[c_{0}]}_{k}(\mu)$, $0\leq k\leq N$. Furthermore, a new projection $\tilde{\mathcal{P}}_{N}$ is defined as

The new hyperbolic regularization in [19] adds one more projection $\tilde{\mathcal{P}}_{N}$ to give the ${H\!M\!P}_{\!N}$ model formulated as

It was revealed that the ${H\!M\!P}_{\!N}$ model enjoys some desired properties, such as

In this section, we adopt the strategy introduced in Section 2 to derive a ${M\!P}_{\!N}$ type model for the 3D RTE at first.