A Positivity-Preserving Finite Volume Scheme for Nonequilibrium Radiation Diffusion Equations on Distorted Meshes

Abstract

In this paper, we propose a new positivity-preserving finite volume scheme with fixed stencils for the nonequilibrium radiation diffusion equations on distorted meshes. This scheme is used to simulate the equations on meshes with both the cell-centered and cell-vertex unknowns. The cell-centered unknowns are the primary unknowns, and the element vertex unknowns are taken as the auxiliary unknowns, which can be calculated by interpolation algorithm. With the nonlinear two-point flux approximation, the interpolation algorithm is not required to be positivity-preserving. Besides, the scheme has a fixed stencil and is locally conservative. The Anderson acceleration is used for the Picard method to solve the nonlinear systems efficiently. Several numerical results are also given to illustrate the efficiency and strong positivity-preserving quality of the scheme.

1. Introduction

The subject capacity of entropy is huge. In thermodynamics, entropy can refer to a physical quantity that can be measured by the change of heat. In many computational problems, the nonequilibrium radiation diffusion equations play a very important role. These real physical applications include the inertial confinement fusion, Z-pinch experiments, and astrophysical problems. When the thermodynamic equilibrium between the radiation field and the material is not reached, a set of coupled radiation diffusion and material conduction equations are needed to simulate the transfer and exchange of energy. With the multiple materials, strongly nonlinear and tightly coupled with the problems, the accurate numerical method is essential to simulate the diffusion processes of these applications.

In recent years, various numerical algorithms [1,2,3] were proposed to obtain reliable numerical solutions of nonequilibrium radiation diffusion equations. The Jacobian-free Newton–Krylov method is given in [4,5] to solve the equations. In [6], the preconditioned Jacobian-free Newton–Krylov methods are considered, and it investigates in [7] a minor improvement to the operator-splitting preconditioner. The time integration methods are presented in [8,9,10], and their efficiency and accuracy are further considered in [11]. The second-order time discretization method for the coupled multidimensional flux-limited nonequilibrium radiation diffusion and material conduction equations is studied in [12]. Moreover, two different time-step control methods were given in [13]. Nevertheless, these methods cannot be used to solve the equations on distorted meshes. The radiation diffusion problems are often combined with other physical processes. For example, in Lagrangian radiation hydrodynamics, the solution is based on hydrodynamic mesh, which is often distorted by fluid motion. Therefore, the constructed numerical method should be able to simulate the radiation diffusion problem on severely distorted meshes.

Positivity-preserving is one of the key requirements for constructing the discrete scheme of nonequilibrium radiation diffusion equations. Against the background of heat transfer, if the discrete scheme does not satisfy this property, it may violate the entropy constraint of the second law of thermodynamics, and it will affect the numerical accuracy of the scheme and bring spurious oscillations. In [14,15], the positivity-preserving finite volume schemes for nonequilibrium radiation diffusion problems are constructed, and the positivity-preserving property of the schemes is also derived. With the applied interpolation method, these methods are only applicable to the problems with geometric restrictions. Moreover, the interpolation algorithms are not positivity-preserving. In [16], two finite volume element methods are presented. It was proved that one of them is monotonic, and another is positivity-preserving under some postprocessing techniques. However, these methods only can be used for the meshes with geometric restrictions. In [17], the authors propose a positivity-preserving finite point method for the nonequilibrium radiation diffusion equations. However, this method is not conservative, and it can only be used for the uniform mesh. In [18], a unified gas-kinetic scheme (UGKS) for a coupled system of radiative transport and material heat conduction with different diffusive limits was constructed, which also only considers uniform mesh.

In this paper, we propose a new positivity-preserving finite volume scheme for the nonequilibrium radiation diffusion equations based on [19,20] This numerical scheme is fixed stencil, local conservative, and positivity-preserving. The numerical scheme has the characteristics of fixed stencil, local conservative, positivity-preserving and so on. The cell-vertexes are used to define auxiliary unknowns. Therefore, vertex interpolation algorithm will be used to obtain the value of the auxiliary unknown. With the distorted meshes, the existing vertex positivity-preserving interpolation algorithms usually have significant accuracy loss. However, this scheme dose not need the interpolation method to be a positivity-preserving one. Besides, the decomposition of normal vectors on unstructured meshes is highly efficient. In practical applications, when the classical Picard iteration method solves the final nonlinear algebraic system, its convergence rate may be very slow. To speed up the convergence of nonlinear iterations, the Anderson acceleration method is used here.

The outline of this paper is as follows. In Section 2, some notations of mesh and the nonequilibrium radiation diffusion equations are presented. In Section 3, we present the construction of the new positivity-preserving finite volume scheme. Some theoretical analysis of discrete schemes is in Section 4. Then, in Section 5 we present some numerical experiments to illustrate the features of the scheme. We end our presentation in Section 6 with some conclusions.

2. Model and Notations

We consider a system of multimaterial nonequilibrium radiation diffusion coupled with material energy balance equation, which is defined in the domain $\mathsf{\Omega }\in {\mathrm{R}}^2$ with the reflection boundary conditions. The equations are written as:

$$\begin{array}{c}\hfill \frac{\partial E}{\partial t}-\nabla ·(D\nabla E)={\sigma }_a(T^4-E),\end{array}$$

(1)

$$\begin{array}{c}\hfill \frac{\partial T}{\partial t}-\nabla ·(\mathsf{\Lambda }\nabla T)={\sigma }_a(E-T^4),\end{array}$$

(2)

where E is defined as the radiation energy density; D is the radiation diffusion coefficient; T is defined as the material temperature; $\mathsf{\Lambda }$ is the corresponding material conduction coefficient; and t is represented as time.

The energy exchange is controlled by the photon absorption cross-section:

$${\sigma }_a\left(T\right)=\frac{Z^3}{T^3},$$

(3)

where Z represents the atomic mass number, and its value is related to the material.

We first define the radiation diffusion coefficient without flux limiter as:

$$D=\frac{1}{3{\sigma }_a}.$$

(4)

However, in regions with strong radiation energy gradients, the model may be unphysical, with the propagation velocity of a radiation wave front in a vacuum faster than the speed of light. To avoid this unphysical phenomenon, we add a limiting term to the diffusion coefficient and adopt such a flux-limited diffusion coefficient:

$$D=\frac{1}{3{\sigma }_a+\frac{|\nabla E|}{E}}.$$

(5)

The material conduction coefficient $\mathsf{\Lambda }$ is taken as the following form:

$$\mathsf{\Lambda }=c_0{T}^{5/2},$$

(6)

where the constant $c_0$ is chosen to be $0.01$.

Throughout this article, we employ the following notations and define the discretization of a finite volume scheme on $\mathsf{\Omega }$ as $\mathcal{D}=(\mathcal{M},\mathcal{E},\mathcal{O},\mathcal{P})$, where

• $\mathcal{M}=\left\{K\right\}$ denotes a family of partitions of the domain $\mathsf{\Omega }$ into nonoverlapping mesh cells, and $\overline{\mathsf{\Omega }}={\cup }_{K\in \mathcal{M}}\overline{K}$. For $K\in \mathcal{M}$, $\partial K$, $h_K$ and $|K|$ denote the cell boundary, diameter (the maximum distance between any two points in K) and measure, respectively. Besides, $h={max}_{K\in \mathcal{M}}{h}_K$ is the mesh sizes;
• $\mathcal{E}=\left\{\sigma \right\}$ is a finite family of disjoint edges in $\overline{\mathsf{\Omega }}$ such that for $\sigma \in \mathcal{E}$, $\sigma$ is a line segment whose measure is defined as $|\sigma |$. Let ${\mathcal{E}}^{int}=\mathcal{E}\cap \mathsf{\Omega }$ and ${\mathcal{E}}^{ext}=\mathcal{E}\cap \partial \mathsf{\Omega }$. For $K\in \mathcal{M}$, there exists a subset ${\mathcal{E}}_K$ of $\mathcal{E}$ such that $\partial K={\cup }_{\sigma \in {\mathcal{E}}_K}\overline{\sigma }$. ${\mathit{n}}_{K,\sigma }$ denotes the unit vector normal to $\sigma$ outward to K;
• $\mathcal{O}=\{{\mathit{x}}_K,K\in \mathcal{M}\}$ is a set of points defined as cell centers, where ${\mathit{x}}_K\in K$;
• $\mathcal{P}={\cup }_{K\in \mathcal{M}}{\mathcal{P}}_K$ is also a set of point, where ${\mathcal{P}}_K=\{{\mathit{x}}_{K,1},{\mathit{x}}_{K,2},\cdots ,{\mathit{x}}_{K,n_K}\}$ represents the set of vertices of cell K, where ${\mathit{x}}_{K,i}(1\le i\le n_K)$ are oriented in a counter clockwise direction. $n_K$ is the number of vertex for cell K.

Let this problem time discrete with the uniform time step $\Delta t$. $E_K$ and $T_K$ are defined as the approximate solutions of E and T at the cell center ${\mathit{x}}_K$ as the primary variables, respectively. The approximation of solutions E and T at the vertex point ${\mathit{x}}_{K,i}$ are defined as the auxiliary variables and denoted as $E_{K,i}$ and $T_{K,i}$.

According to the idea of finite volume framework, the above Equations (1) and (2) are integrated into each control volume K, and we can then obtain:

$$\begin{array}{c}\hfill {\int }_K\frac{\partial E}{\partial t}\mathrm{d}\mathit{x}+\sum _{\sigma \in {\mathcal{E}}_K}{\mathcal{F}}_{E,K,\sigma }={\int }_K{\sigma }_a(T^4-E)\mathrm{d}\mathit{x},\\ \hfill {\int }_K\frac{\partial T}{\partial t}\mathrm{d}\mathit{x}+\sum _{\sigma \in {\mathcal{E}}_K}{\mathcal{F}}_{T,K,\sigma }={\int }_K{\sigma }_a(E-T^4)\mathrm{d}\mathit{x},\end{array}$$

where ${\mathcal{F}}_{E,K,\sigma }$ and ${\mathcal{F}}_{T,K,\sigma }$ are the fluxes $-{\int }_{\sigma }(D\nabla E)·{\mathit{n}}_{K,\sigma }\mathrm{d}s$ and $-{\int }_{\sigma }(\mathsf{\Lambda }\nabla T)·{\mathit{n}}_{K,\sigma }\mathrm{d}s$, respectively. $F_{E,K,\sigma }$ (resp.$F_{T,K,\sigma }$) is the one-sided flux that approximates ${\mathcal{F}}_{E,K,\sigma }$ (resp.${\mathcal{F}}_{T,K,\sigma }$) using only the information of cell K.

3. Methods

In this section, we will make use of four steps to construct a new positivity-preserving finite volume scheme.

3.1. Construction of One-Sided Flux

In this article, we consider that the cells $K\in \mathcal{M}$ are star-shaped. The geometric center ${\mathit{x}}_K$ is the average of the cell vertices and taken as the cell center. Under the linear preserving approach [20,21], we can use the cell center ${\mathit{x}}_K$ and the interpolation points ${\mathcal{P}}_K$ to obtain an approximation of the one-sided flux. The algebraic form of the flux is as follows:

$${\mathbf{F}}_{E,K}={\mathbb{A}}_{E,K}\left(E_K{\mathbf{I}}_K-{\mathbf{E}}_K\right),{\mathbf{F}}_{T,K}={\mathbb{A}}_{T,K}\left(T_K{\mathbf{I}}_K-{\mathbf{T}}_K\right),$$

(7)

where ${\mathbf{F}}_{E,K}=\left(F_{E,K,\sigma },\sigma \in {\mathcal{E}}_K\right)$ and ${\mathbf{E}}_K=\left(E_{K,i},{\mathit{x}}_{K,i}\in {\mathcal{P}}_K\right)$ have $n_K$ components, ${\mathbb{A}}_{E,K}$ is an $n_K\times n_K$ cell matrix, and ${\mathbf{I}}_K$ is a vector whose components are all 1, and its size is $n_K$. The ${\mathbf{F}}_{T,K}$, ${\mathbf{T}}_K$ and ${\mathbb{A}}_{T,K}$ are defined similarly. Now the construction procedure of the cell matrix ${\mathbb{A}}_{E,K}$ will be presented. The cell matrix ${\mathbb{A}}_{T,K}$ can be obtained by the same procedure.

For any $\sigma \in {\mathcal{E}}_K$, there are two endpoints ${\mathit{x}}_{K,i},{\mathit{x}}_{K,i+1}$ of $\sigma$ as in Figure 1. The following decomposition can be obtained:

$$|\sigma |{\mathit{n}}_{K,\sigma }={\alpha }_{K,\sigma ,i}\left({\mathit{x}}_{K,i}-{\mathit{x}}_K\right)+{\alpha }_{K,\sigma ,i+1}\left({\mathit{x}}_{K,i+1}-{\mathit{x}}_K\right),$$

(8)

where

$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill {\alpha }_{K,\sigma ,i}& =\hfill & \hfill \frac{|\sigma |{\mathit{n}}_{K,\sigma }^{\prime }\mathcal{R}({\mathit{x}}_{K,i+1}-{\mathit{x}}_K)}{{({\mathit{x}}_{K,i}-{\mathit{x}}_K)}^{\prime }\mathcal{R}({\mathit{x}}_{K,i+1}-{\mathit{x}}_K)},\\ \hfill {\alpha }_{K,\sigma ,i+1}& =\hfill & \hfill \frac{|\sigma |{\mathit{n}}_{K,\sigma }^{\prime }\mathcal{R}({\mathit{x}}_{K,i}-{\mathit{x}}_K)}{{({\mathit{x}}_{K,i+1}-{\mathit{x}}_K)}^{\prime }\mathcal{R}({\mathit{x}}_{K,i}-{\mathit{x}}_K)}.\end{array}\end{array}$$

(9)

$\mathcal{R}$ denotes an operator that rotates a vector clockwise to its normal direction, ${\mathit{n}}_{K,\sigma }^{\prime }$ denotes the transposition of vector ${\mathit{n}}_{K,\sigma }$. Based on the decomposition (8), the linearity-preserving one-sided flux of radiation energy is given as follows:

$$F_{E,K,\sigma }=D_{K,\sigma }[{\alpha }_{K,\sigma ,i}(E_K-E_{K,i})+{\alpha }_{K,\sigma ,i+1}(E_K-E_{K,i+1})],\sigma \in {\mathcal{E}}_K.$$

(10)

The linearity-preserving one-sided flux of the material temperature T is defined similarly:

$$F_{T,K,\sigma }={\mathsf{\Lambda }}_{K,\sigma }[{\alpha }_{K,\sigma ,i}(T_K-T_{K,i})+{\alpha }_{K,\sigma ,i+1}(T_K-T_{K,i+1})],\sigma \in {\mathcal{E}}_K.$$

(11)

The following conditions need to be satisfied to make the discrete scheme to remain the positivity-preserving property:

$${\alpha }_{K,\sigma ,i}+{\alpha }_{K,\sigma ,i+1}>0,\sigma \in {\mathcal{E}}^{int},$$

(12)

Theorem 1.

If K is a star-shaped polygon with respect to its cell-centered ${\mathit{x}}_K$. With the coefficients of one-sided flux $F_{E,K,\sigma }$, $F_{T,K,\sigma }$ given by (10) and (11), the inequality (12) will be satisfied.

Proof. 

The star-shaped cell K in the mesh:

$$\begin{array}{c}\hfill {({\mathit{x}}_{K,i}-{\mathit{x}}_K)}^{\prime }\mathcal{R}({\mathit{x}}_{K,i+1}-{\mathit{x}}_K)=-{({\mathit{x}}_{K,i+1}-{\mathit{x}}_K)}^{\prime }\mathcal{R}({\mathit{x}}_{K,i}-{\mathit{x}}_K)=|\sigma |d_{K,\sigma },\end{array}$$

can be obtained by considering its geometric relationship.

Then, we deduce from (9) that:

$${\alpha }_{K,\sigma ,i}+{\alpha }_{K,\sigma ,i+1}=\frac{|\sigma |{\mathit{n}}_{K,\sigma }^{\prime }\mathcal{R}({\mathit{x}}_{K,i+1}-{\mathit{x}}_{K,i})}{{({\mathit{x}}_{K,i}-{\mathit{x}}_K)}^{\prime }\mathcal{R}({\mathit{x}}_{K,i+1}-{\mathit{x}}_K)}=\frac{|\sigma |}{d_{K,\sigma }},$$

(13)

where $\mathcal{R}({\mathit{x}}_{K,i+1}-{\mathit{x}}_{K,i})=|\sigma |n_{K,\sigma }$. Thus, the proof is completed.    □

Remark 1.

We directly define $D_{K,\sigma }={\left(3\frac{z_K^3}{T_{K,\sigma }^3}\right)}^{-1}$ with the $T_{K,\sigma }=\frac{d_{K,\sigma }}{d_{K,\sigma }+d_{L,\sigma }}{T}_K+\frac{d_{L,\sigma }}{d_{K,\sigma }+d_{L,\sigma }}{T}_L$, to represent the radiation diffusion coefficient without flux limiter.

The material conduction coefficient is taken as ${\mathsf{\Lambda }}_{K,\sigma }=c_0{\left(T_{K,\sigma }\right)}^{5/2}$. For the radiation diffusion coefficient with flux limiter, we take:

$$\begin{array}{c}\hfill D_{K,\sigma }=\frac{1}{3\frac{z_K^3}{T_{K,\sigma }^3}+\frac{|\nabla E_{K,\sigma }|}{E_{K,\sigma }}},\end{array}$$

where $\nabla E_{K,\sigma }$ is an approximation of $\nabla E$ on triangle $▵{\mathit{x}}_K{\mathit{x}}_{K,i}{\mathit{x}}_{K,i+1}$ and $E_{K,\sigma }=\frac{d_{K,\sigma }}{d_{K,\sigma }+d_{L,\sigma }}{E}_K+\frac{d_{L,\sigma }}{d_{K,\sigma }+d_{L,\sigma }}{E}_L$. ${\mathsf{\Omega }}_{K,\sigma }$ denotes the triangle $▵{\mathit{x}}_K{\mathit{x}}_{K,i}{\mathit{x}}_{K,i+1}$ (see Figure 1).

Then, by the Green formulas:

$$\begin{array}{cc}\hfill {\int }_{{\mathsf{\Omega }}_{K,\sigma }}\nabla Ed\mathit{x}& ={\int }_{\partial {\mathsf{\Omega }}_{K,\sigma }}E\mathit{n}ds\hfill \\ \hfill & \approx \frac{1}{2}(E_{K,i}+E_{K,i+1})|\sigma |{\mathit{n}}_{K,\sigma }+\frac{1}{2}(E_K+E_{K,i+1})|{\mathit{x}}_{K,i+1}{\mathit{x}}_K|{\mathit{n}}_{{\mathit{x}}_{K,i+1}{\mathit{x}}_K}\hfill \\ \hfill & +\frac{1}{2}(E_K+E_{K,i})|{\mathit{x}}_K{\mathit{x}}_{K,i}|{\mathit{n}}_{{\mathit{x}}_K{\mathit{x}}_{K,i}}\hfill \\ \hfill & =\frac{1}{2}(E_K-E_{K,i+1})|{\mathit{x}}_K{\mathit{x}}_{K,i}|{\mathit{n}}_{{\mathit{x}}_K{\mathit{x}}_{K,i}}+\frac{1}{2}(E_K-E_{K,i})|{\mathit{x}}_{K,i+1}{\mathit{x}}_K|{\mathit{n}}_{{\mathit{x}}_{K,i+1}{\mathit{x}}_K}.\hfill \end{array}$$

Here, ${\mathit{n}}_{{\mathit{x}}_K{\mathit{x}}_{K,i}}$ (resp.,${\mathit{n}}_{{\mathit{x}}_{K,i+1}{\mathit{x}}_K}$) denotes the unit outward normal to ${\mathit{x}}_K{\mathit{x}}_{K,i}$ (resp., ${\mathit{x}}_{K,i+1}{\mathit{x}}_K$). It holds that:

$$\begin{array}{c}\hfill \frac{1}{2}|{\mathit{x}}_K{\mathit{x}}_{K,i}||{\mathit{x}}_{K,i+1}{\mathit{x}}_K|sin{\theta }_K\nabla E_{K,\sigma }\\ \hfill =\frac{1}{2}|{\mathit{x}}_K{\mathit{x}}_{K,i}||{\mathit{x}}_{K,i+1}{\mathit{x}}_K|\left(\frac{E_K-E_{K,i+1}}{|{\mathit{x}}_K{\mathit{x}}_{K,i+1}|}{\mathit{n}}_{{\mathit{x}}_K{\mathit{x}}_{K,i}}+\frac{E_K-E_{K,i}}{|{\mathit{x}}_K{\mathit{x}}_{K,i}|}{\mathit{n}}_{{\mathit{x}}_{K,i+1}{\mathit{x}}_K}\right).\end{array}$$

Then, we have:

$$\begin{array}{c}\hfill \nabla E_{K,\sigma }=\frac{1}{sin{\theta }_K}\left(\frac{E_K-E_{K,i+1}}{|{\mathit{x}}_K{\mathit{x}}_{K,i+1}|}{\mathit{n}}_{{\mathit{x}}_K{\mathit{x}}_{K,i}}+\frac{E_K-E_{K,i}}{|{\mathit{x}}_K{\mathit{x}}_{K,i}|}{\mathit{n}}_{{\mathit{x}}_{K,i+1}{\mathit{x}}_K}\right).\end{array}$$

3.2. A Unique Definition of the Facet Flux

For the interior edge $\{\sigma \in {\mathcal{E}}_K\cap {\mathcal{E}}_L\}$, the $F_{E,K,\sigma }$ and $F_{T,K,\sigma }$ are defined as (10) and (11). With the continuity of normal flux component, a linear combination of the one-sided fluxes is defined as follows:

$${\tilde{F}}_{\vartheta ,K,\sigma }={\mu }_{\vartheta ,K,\sigma }{F}_{\vartheta ,K,\sigma }-{\mu }_{\vartheta ,L,\sigma }{F}_{\vartheta ,L,\sigma },{\tilde{F}}_{\vartheta ,L,\sigma }={\mu }_{\vartheta ,L,\sigma }{F}_{\vartheta ,L,\sigma }-{\mu }_{\vartheta ,K,\sigma }{F}_{\vartheta ,K,\sigma },$$

(14)

where ${\mu }_{\vartheta ,K,\sigma }$, ${\mu }_{\vartheta ,L,\sigma }$ are the parameters

$${\mu }_{\vartheta ,K,\sigma }+{\mu }_{\vartheta ,L,\sigma }=1,\vartheta =E,T.$$

(15)

Obviously, we have:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\tilde{F}}_{\vartheta ,K,\sigma }+{\tilde{F}}_{\vartheta ,L,\sigma }=0,\vartheta =E,T,\sigma \in {\mathcal{E}}^{int}.\end{array}\end{array}$$

(16)

By taking (10) into the first equation of (14), we will obtain:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\tilde{F}}_{E,K,\sigma }=\left({\mu }_{E,K,\sigma }{D}_{K,\sigma }\sum _{i=1}^{n_K}{\alpha }_{K,\sigma ,i}\right)E_K-\left({\mu }_{E,L,\sigma }{D}_{L,\sigma }\sum _{i=1}^{n_L}{\alpha }_{L,\sigma ,i}\right)E_L\\ \hfill +{\mu }_{E,L,\sigma }{a}_{E,L,\sigma }-{\mu }_{E,K,\sigma }{a}_{E,K,\sigma },\end{array}\end{array}$$

(17)

where:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill a_{E,K,\sigma }=D_{K,\sigma }\sum _{i=1}^{n_K}{\alpha }_{K,\sigma ,i}{E}_{K,i},a_{E,L,\sigma }=D_{L,\sigma }\sum _{i=1}^{n_L}{\alpha }_{L,\sigma ,i}{E}_{L,i}.\end{array}\end{array}$$

Based on (15) and (17), we define:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\mu }_{E,K,\sigma }=\frac{|a_{E,L,\sigma }|+ϵ}{|a_{E,K,\sigma }|+|a_{E,L,\sigma }|+2ϵ},\end{array}\end{array}$$

(18)

and ${\mu }_{E,L,\sigma }=1-{\mu }_{E,K,\sigma }$.

The set

$$\begin{array}{c}\hfill B_{E,\sigma }=\left\{\begin{array}{c}{\mu }_{E,L,\sigma }{a}_{E,L,\sigma }-{\mu }_{E,K,\sigma }{a}_{E,K,\sigma },\sigma \in {\mathcal{E}}_K\cap {\mathcal{E}}_L,\hfill \\ -a_{E,K,\sigma },\sigma \in {\mathcal{E}}_K\cap {\mathcal{E}}^{ext},\hfill \end{array}\right.\end{array}$$

(19)

and

$$\begin{array}{c}\hfill \begin{array}{c}\hfill B_{E,\sigma }^{+}=\frac{|B_{E,\sigma }|+B_{E,\sigma }}{2},B_{E,\sigma }^{-}=\frac{|B_{E,\sigma }|-B_{E,\sigma }}{2}.\end{array}\end{array}$$

Then, (17) can be rewritten as:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\tilde{F}}_{E,K,\sigma }=A_{E,K,\sigma }{E}_K-A_{E,L,\sigma }{E}_L+B_{E,\sigma }^{ϵ},\end{array}\end{array}$$

(20)

where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill A_{E,K,\sigma }& ={\mu }_{E,K,\sigma }{D}_{K,\sigma }\sum _{i=1}^{n_K}{\alpha }_{K,\sigma ,i}+\frac{B_{E,\sigma }^{+}}{E_K+ϵ},\hfill \\ \hfill A_{E,L,\sigma }& ={\mu }_{E,L,\sigma }{D}_{L,\sigma }\sum _{i=1}^{n_L}{\alpha }_{L,\sigma ,i}+\frac{B_{E,\sigma }^{-}}{E_L+ϵ},\hfill \end{array}\end{array}$$

(21)

and

$$\begin{array}{c}\hfill \begin{array}{c}\hfill B_{E,\sigma }^{ϵ}=\frac{B_{E,\sigma }^{+}ϵ}{E_K+ϵ}-\frac{B_{E,\sigma }^{-}ϵ}{E_L+ϵ}.\end{array}\end{array}$$

Here, $ϵ$ denotes a positive number with the machine precision.

Finally, ignoring $B_{E,\sigma }^{ϵ}$ in (20) and by (14), a new definition of the unique edge flux is obtained:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\tilde{F}}_{E,K,\sigma }& =A_{E,K,\sigma }{E}_K-A_{E,L,\sigma }{E}_L,\hfill \\ \hfill {\tilde{F}}_{E,L,\sigma }& =A_{E,L,\sigma }{E}_L-A_{E,K,\sigma }{E}_K.\hfill \end{array}\end{array}$$

(22)

For the face $\sigma \in {\mathcal{E}}_K\cap {\mathsf{\Gamma }}_D$, we can simply define:

$${\tilde{F}}_{E,K,\sigma }=F_{E,K,\sigma }=\left(\sum _{i=1}^{n_K}{\alpha }_{K,\sigma ,i}\right)E_K+B_{E,\sigma }^{+}-B_{E,\sigma }^{-},$$

where $B_{E,\sigma }^{+}$ can be treated in a way similar to the interior face, while $B_{E,\sigma }^{-}$ moves to the right of the final finite volume equations.

For the face $\sigma \in {\mathcal{E}}_K\cap {\mathsf{\Gamma }}_N$, we can obtain the corresponding ${\tilde{F}}_{E,K,\sigma }$ by integrating the Neumann boundary on $\sigma$.

The one-sided flux of the material temperature T can be defined similarly.

3.3. Interpolation of the Auxiliary Variables

To make the finite volume scheme cell-centered, we use interpolation procedure to eliminate the intermediate vertex variables in the flux expression. We express the intermediate variables as linear combinations of the primary variables. As shown in Figure 2, we might as well assume that the interior cell vertex ${\mathit{x}}_{\upsilon }$ is surrounded by cells $K_i$$(1\le i\le \gamma )$. Let ${\overline{\mathit{x}}}_{{\sigma }_i}$ is the midpoint of ${\sigma }_i$, and ${\mathit{x}}_{K_i}$ denotes the cell center of $K_i$. The interior vertex ${\mathit{x}}_{\upsilon }$ is the intersection of edges ${\sigma }_i$ and ${\sigma }_{i+1}$. Considering the geometry of the mesh, we now that $K_{\gamma +1}=K_1$, ${\sigma }_{\gamma }={\sigma }_0$, etc. In this section, we let $E_{\upsilon }$ represent the unknowns defined at the vertex of ${\mathit{x}}_{\upsilon }$. $D_{K_i}$ and $E_{K_i}$ are the values of D and E on the primary unknown at ${\mathit{x}}_{K_i}$, respectively. $S_{i,j}$ represents the area of triangle $▵{\mathit{x}}_{K_i}{\mathit{x}}_{\upsilon }{\overline{\mathit{x}}}_{{\sigma }_{i+j-1}}(j=1,2)$.

For $i=1,2,\cdots ,\gamma$ and $j=1,2$, we define:

$$\begin{array}{c}\hfill {\xi }_{i,j}=\frac{D_{K_i}}{2S_{i,j}}{|{\overline{\mathit{x}}}_{{\sigma }_{i+j-1}}-{\mathit{x}}_{\upsilon }|}^2,\\ \hfill {\overline{\xi }}_{i,j}=\frac{D_{K_i}}{2S_{i,j}}{({\overline{\mathit{x}}}_{{\sigma }_{i+j-1}}-{\mathit{x}}_{\upsilon })}^{\prime }({\mathit{x}}_{\upsilon }-{\mathit{x}}_{K_i}),\end{array}$$

and

$$\begin{array}{c}\hfill {\eta }_{i,j}=\frac{D_{K_i}}{2S_{▵{\mathit{x}}_{\upsilon }{\overline{\mathit{x}}}_{{\sigma }_i}{\overline{\mathit{x}}}_{{\sigma }_{i+1}}}}{({\overline{\mathit{x}}}_{{\sigma }_{i+1}}-{\overline{\mathit{x}}}_{{\sigma }_i})}^{\prime }({\overline{\mathit{x}}}_{{\sigma }_{i+j-1}}-{\mathit{x}}_{\upsilon }).\end{array}$$

Then, based on the interpolation algorithm LPEW2, we obtain the following interpolation formula:

$$\begin{array}{c}\hfill E_{\upsilon }=\sum _{i=1}^{\gamma }{\omega }_i{E}_{K_i},\end{array}$$

where

$$\begin{array}{c}\hfill {\omega }_i=\frac{{\overline{\omega }}_i}{{\sum }_{k=1}^{\gamma }{\overline{\omega }}_k}\mathrm{and}{\overline{\omega }}_k=\frac{{\eta }_{k-1,1}-{\eta }_{k,2}}{{\overline{\xi }}_{k-1,2}+{\overline{\xi }}_{k,1}}{\xi }_{k,1}+\frac{{\eta }_{k,1}-{\eta }_{k+1,2}}{{\overline{\xi }}_{k,2}+{\overline{\xi }}_{k+1,1}}{\xi }_{k,2}.\end{array}$$

The auxiliary variables of unknown T in flux expressions also can be computed as the above formulations.

Remark 2.

With flux discrete (10) and (11), the auxiliary variables on boundary may be used. The vertex values on the boundary can be obtained using the interpolation algorithm in [20]. For the flux-limited diffusion coefficient $D_{K_i}$, we have the discrete $\nabla E$ on cell $K_i$:

$$\begin{array}{cc}\hfill \nabla E|K_i\approx \frac{1}{|K_i|}{\int }_{K_i}\nabla Ed\mathit{x}& =\frac{1}{|K_i|}\sum _{\sigma \in {\mathcal{E}}_{K_i}}{\int }_{\sigma }E{\mathbf{n}}_{K_i,\sigma }ds\hfill \\ \hfill & \approx \frac{1}{2|K_i|}\sum _{\sigma \in {\mathcal{E}}_{K_i}}|\sigma |(E_{K_i,j}+E_{K_i,j+1}){\mathbf{n}}_{K_i,\sigma }.\hfill \end{array}$$

3.4. The Finite Volume Scheme

To obtain the complete discretization scheme of the Equations (1) and (2), we use backward Euler to discretize time. Based on the definitions of ${\tilde{F}}_{E,K,\sigma }$ and ${\tilde{F}}_{T,K,\sigma }$, we formulate the general cell-centered positivity-preserving finite volume scheme as follows:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \frac{E^{n+1}-E^n}{\Delta t}m\left(K\right)+\sum _{\sigma \in {\epsilon }_K}{\tilde{F}}_{E,K,\sigma }^{n+1}={\sigma }_{a,K}^{n+1}({\left(T_K^{n+1}\right)}^4-E_K^{n+1})m\left(K\right),\\ \hfill \frac{T^{n+1}-T^n}{\Delta t}m\left(K\right)+\sum _{\sigma \in {\epsilon }_K}{\tilde{F}}_{T,K,\sigma }^{n+1}={\sigma }_{a,K}^{n+1}(E_K^{n+1}-{\left(T_K^{n+1}\right)}^4)m\left(K\right),\end{array}\end{array}$$

(23)

where ${\sigma }_{a,K}^{n+1}=\frac{z_K^3}{{\left(T_K^{n+1,k}\right)}^3}$, $m\left(K\right)$ is the measure of cell K, and ${\left(T_K^{n+1}\right)}^4\approx {\left(T_K^{n+1,k}\right)}^3{T}_K^{n+1,k+1}$.

The nonlinear scheme (23) can be expressed in the following matrix form:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \mathbb{M}\left({\mathbf{U}}^{n+1}\right){\mathbf{U}}^{n+1}=\mathbf{F}\left({\mathbf{U}}^n\right),\end{array}\end{array}$$

(24)

where ${\mathbf{U}}^{n+1}$ denotes the vector unknowns of cell-centroid, the right-hand side vector $\mathbf{F}\left({\mathbf{U}}^n\right)$ is obtained by the known value, and $\mathbb{M}\left(\mathbf{U}\right)$ is the coefficient matrix.

This nonlinear algebraic system can be solved by the Picard iterative method (Algorithm 1).

Algorithm 1 Picard method.

1: Choose a small positive value ${\epsilon }_{non}$ and take ${\mathbf{U}}^{n+1,0}={\mathbf{U}}^n\ge 0$. 2: for$k=1,2,...$do 3:     Solve the linear system $$\begin{array}{c}\hfill \begin{array}{c}\hfill \mathbb{M}\left({\mathbf{U}}^{n+1,k}\right){\mathbf{U}}^{n+1,k+1}=\mathbf{F}\left({\mathbf{U}}^n\right)\end{array}\end{array}$$ (25) to obtain ${\mathbf{U}}^{n+1,k+1}$. 4:     Use $$\begin{array}{c}\hfill \begin{array}{c}\hfill \parallel \mathbb{M}\left({\mathbf{U}}^{n+1,k+1}\right){\mathbf{U}}^{n+1,k+1}-\mathbf{F}\left({\mathbf{U}}^{n+1,k+1}\right)\parallel \le {\epsilon }_{non}\parallel \mathbb{M}\left({\mathbf{U}}^n\right){\mathbf{U}}^n-\mathbf{F}\left({\mathbf{U}}^n\right)\parallel \end{array}\end{array}$$ (26) as the convergence criterion. 5: end for

To speed up the convergence of nonlinear iterations, we employ the Anderson acceleration for the Picard iteration, which is formalized as the following Algorithm 2. The parameter m is a fixed positive integer, which affects the efficiency of the iteration. The solution of constrained minimization problem (27) is shifted to the solution of a saddle-point problem [22].

Algorithm 2 Anderson acceleration of Picard method.

1: Give a small positive value ${\epsilon }_{non}$, and define ${\mathbf{U}}^{n+1,0}={\mathbf{U}}^n\ge 0$. 2: Apply two Picard iterations (25) and set ${\tilde{\mathbf{U}}}^{n+1,k}={\mathbf{U}}^{n+1,k}$, $\delta {\mathbf{U}}^{n+1,k}={\tilde{\mathbf{U}}}^{n+1,k}-{\mathbf{U}}^{n+1,k-1}$, $k=1,2$. 3: for$k=2,...$do 4:     Set $m_k=\min (m,k)$. 5:     Solve the minimization problem $$\begin{array}{c}\hfill \begin{array}{c}\hfill \min \parallel \sum _{i=1}^{m_k}{\alpha }_i\delta {\mathbf{U}}^{n+1,k-m_k+i}\parallel \mathrm{with}\sum _{i=1}^{m_k}{\alpha }_i=1.\end{array}\end{array}$$ (27) 6:     Set $${\mathbf{U}}^{n+1,k+1}=\sum _{i=1}^{m_k}{\alpha }_i{\tilde{\mathbf{U}}}^{n+1,k-m_k+i}.$$ 7:     Let ${\mathbf{U}}^{n+1,k+1}:={\mathbf{U}}^{n+1,k+1}-\min \{{\mathbf{U}}^{n+1,k+1},0\}\mathbf{e}$, where $\mathbf{e}$ is a vector with all entries equal to 1. 8:     Set ${\mathbf{U}}^{n+1}={\mathbf{U}}^{n+1,k+1}$, while the criterion (26) is satisfied. 9:     Solve the linear system $$\mathbb{M}\left({\mathbf{U}}^{n+1,k+1}\right){\tilde{\mathbf{U}}}^{n+1,k+1}=\mathbf{F}\left({\mathbf{U}}^n\right)$$ and set $\delta {\mathbf{U}}^{n+1,k+1}={\tilde{\mathbf{U}}}^{n+1,k+1}-{\mathbf{U}}^{n+1,k+1}$. 10: end for

4. Theoretical Results

In this section, we consider the theoretical results of the finite volume scheme. The positivity-preserving property of the scheme is first proved. Then, the compatibility of the discrete solution and discrete flux of the scheme is given. Finally, we give the existence of solution for nonlinear discrete system.

4.1. The Positivity-Preserving

From the construction procedure of the finite volume scheme, we have the following theorem:

Theorem 2.

Assume that the initial solution vector ${\mathbf{U}}^0\ge 0$ and the linear algebraic equations formed by Picard iteration are solved exactly. Then, we have ${\mathbf{U}}^{n+1,k+1}\ge 0$ for $n\ge 0$ and $k\ge 0$.

Proof. 

The nonlinear discrete system (24) obtained by the positivity-preserving finite volume scheme, in which the matrix $\mathbb{M}\left({\mathbf{U}}^{n+1,k}\right)$ has the following properties:

• All diagonal elements of matrix $\mathbb{M}\left({\mathbf{U}}^{n+1,k}\right)$ are greater than zero.
• All off-diagonal elements of matrix $\mathbb{M}\left({\mathbf{U}}^{n+1,k}\right)$ are not less than zero.
• The sum of each column of matrix $\mathbb{M}\left({\mathbf{U}}^{n+1,k}\right)$ is greater than zero.

Obviously, the matrix ${\mathbb{M}}^{\prime }\left({\mathbf{U}}^{n+1,k}\right)$ is an M-matrix. Therefore, it is further obtained that the matrix ${\mathbb{M}}^{-1}\left({\mathbf{U}}^{n+1,k}\right)$ is nonsingular. With the nonnegative vector $\mathbf{F}\left({\mathbf{U}}^n\right)$, one nonnegative solution ${\mathbf{U}}^{n+1,k+1}$ will be obtained. □

4.2. The Compatibility

Assumption 1.

Assume that any mesh cell K is star-shaped; that is, the ray starting from the center of cell K intersects with only one edge of the mesh. Moreover, we assume that:

$$E\left(x\right),T\left(x\right)\in C^2\left(\overline{K}\right),D,\mathsf{\Lambda }\in C^1\left(\overline{K}\right),{\sigma }_a\in C\left(\overline{K}\right)$$

Obviously, some concave meshes are still star-shaped.

Theorem 3.

Suppose that an interior vertex ${\mathit{x}}_{\upsilon }$ is surrounded by cell $K_i$ centered at ${\mathit{x}}_{K_i}$, and ${\omega }_{E,i},{\omega }_{T,i}$ are the weight coefficients of the corresponding cell center ${\mathit{x}}_{K_i}$, γ is the number of adjacent cells of vertex ${\mathit{x}}_{\upsilon }$. Under the Assumption 1, there is a constant C independent of h, so that

$$|E_{\upsilon }-\sum _{i=1}^{\gamma }{\omega }_{E,i}{E}_{K_i}|\le C{h}^2$$

$$|T_{\upsilon }-\sum _{i=1}^{\gamma }{\omega }_{T,i}{T}_{K_i}|\le C{h}^2$$

Proof. 

According to the article [20], and the derivation process of vertex weight coefficient is the linearity-preserving, it is clear that the first formula in the lemma holds, and the same holds for the second formula. □

Theorem 4.

Suppose that $E\left(x\right),T\left(x\right)\in C^2\left(\overline{\mathsf{\Omega }}\right)$, under the Assumption 1, there is a constant C independent of h, such that:

$$\frac{1}{\sigma }|{\tilde{F}}_{E,K,\sigma }-{\mathcal{F}}_{E,K,\sigma }|\le Ch$$

$$\frac{1}{\sigma }|{\tilde{F}}_{T,K,\sigma }-{\mathcal{F}}_{T,K,\sigma }|\le Ch$$

Proof. 

(10) shows that:

$${\mathcal{F}}_{E,K,\sigma }=D_{K,\sigma }[{\alpha }_{K,\sigma ,i}(E_K-E_{K,i})+{\alpha }_{K,\sigma ,i+1}(E_K-E_{K,i+1})]+O\left(h^2\right)=F_{E,K,\sigma }+O\left(h^2\right).$$

From the Theorem 3 we known that:

$$\frac{1}{\sigma }|{\mathcal{F}}_{E,K,\sigma }-F_{E,K,\sigma }|\le Ch$$

Similarly:

$$\frac{1}{\sigma }|{\mathcal{F}}_{E,L,\sigma }-F_{E,L,\sigma }|\le Ch$$

with

$${\tilde{F}}_{E,K,\sigma }=-{\tilde{F}}_{E,L,\sigma }$$

$${\tilde{F}}_{E,K,\sigma }={\mu }_{E,K,\sigma }{F}_{E,K,\sigma }-{\mu }_{E,L,\sigma }{F}_{E,L,\sigma }$$

$0<{\mu }_{E,K,\sigma },{\mu }_{E,L,\sigma }\le 0,{\mu }_{E,K,\sigma }+{\mu }_{E,L,\sigma }=1$, and $|\sigma |=O\left(h\right)$ then

$$\frac{1}{\sigma }|{\tilde{F}}_{E,K,\sigma }-{\mathcal{F}}_{E,K,\sigma }|\le Ch$$

Similarly, the second formula holds. □

4.3. The Existence of Solution

In this subsection, we only consider the existence of discrete solution for our scheme.

Theorem 5.

Given ${\phi }_{\alpha }\left(\mathit{x}\right)\ge 0$ and $g_{\alpha }\left(\mathit{x}\right)\ge 0,\alpha \in \{e,i,r\}$, there exists a solution ${\mathbf{U}}^{n+1}\ge 0$ for the system (24).

Proof. 

The system (24) can be written as ${\mathbf{U}}^{n+1}=M\left({\mathbf{U}}^{n+1}\right)$ and:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill M\left({\mathbf{U}}^{n+1}\right)={(\mathbf{D}\left({\mathbf{U}}^{n+1}\right)+\Delta t\mathbf{A}\left({\mathbf{U}}^{n+1}\right))}^{-1}\mathbf{F}\left({\mathbf{U}}^n\right).\end{array}\end{array}$$

To prove the existence of solution for system (24), we need to prove that the map M has a fixed point ${\mathbf{U}}^{n+1}\ge 0$. The ${(\mathbf{D}\left({\mathbf{U}}^{n+1}\right)+\Delta t\mathbf{A}\left({\mathbf{U}}^{n+1}\right))}^{\prime }$ is an M-matrix. Combined with $\mathbf{F}\left({\mathbf{U}}^n\right)\ge 0$, we have:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \forall {\mathbf{U}}^{n+1}\in {\mathrm{R}}^{2\mathbb{N}},{\mathbf{U}}^{n+1}\ge 0,M\left({\mathbf{U}}^{n+1}\right)\ge 0.\end{array}\end{array}$$

From the discretization of flux (22), we know that most column sum of matrix $\mathbf{A}\left({\mathbf{U}}^{n+1}\right)$ is zero. Then, multiplying system (24) with a constant vector $(1,\dots ,1)$, we have:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \sum _{\alpha =e,i,r}\sum _{K\in \mathcal{M}}{\left[\mathbf{D}\left({\mathbf{U}}^{n+1}\right)M\left({\mathbf{U}}^{n+1}\right)\right]}_{\alpha ,K}\le \sum _{\alpha =e,i,r}\sum _{K\in \mathcal{M}}{\left[\mathbf{F}\left({\mathbf{U}}^n\right)\right]}_{\alpha ,K}\le C_1.\end{array}\end{array}$$

Further, we can obtain:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \underset{1\le k\le 2\mathbb{N}}{min}\left\{{\mathbf{D}}_{k,k}\left({\mathbf{U}}^{n+1}\right)\right\}\sum _{\alpha =e,i,r}\sum _{K\in \mathcal{M}}{\left[M\left({\mathbf{U}}^{n+1}\right)\right]}_{\alpha ,K}\le \sum _{\alpha =e,i,r}\sum _{K\in \mathcal{M}}{\left[\mathbf{D}\left({\mathbf{U}}^{n+1}\right)M\left({\mathbf{U}}^{n+1}\right)\right]}_{\alpha ,K}\le C_1.\end{array}\end{array}$$

Thus, we have:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \sum _{\alpha =e,i,r}\sum _{K\in \mathcal{M}}{\left[M\left({\mathbf{U}}^{n+1}\right)\right]}_{\alpha ,K}\le C_0.\end{array}\end{array}$$

Combining with $M\left({\mathbf{U}}^{n+1}\right)\ge 0$, we obtain:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\left[M\left({\mathbf{U}}^{n+1}\right)\right]}_{\alpha ,K}\le C_0,\forall K\in \mathcal{M},\alpha \in \{e,i,r\}.\end{array}\end{array}$$

Define a convex compact subset $\mathbf{C}\in {\mathrm{R}}^{2\mathbb{N}}$ as follows:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \mathbf{C}=\{{\mathbf{U}}^n\in {\mathrm{R}}^{2\mathbb{N}}|0\le T_{\alpha ,K}^n\le C_0,\forall K\in \mathcal{M},\alpha =e,i,r\}\end{array}\end{array}$$

For $\mathbf{w}\in \mathbf{C}$, we know $M\left(\mathbf{w}\right)\in \mathbf{C}$; that is, M maps $\mathbf{C}$ into itself. For each $\mathbf{w}\in \mathbf{C}$, there is only one solution for $M\left(\mathbf{w}\right)=\mathbf{F}\left({\mathbf{U}}^n\right)$. Hence, the map M is well-defined and $\mathbb{M}\left(\mathbf{w}\right)$ is an invertible matrix, and ${\mathbf{U}}^n=M\left(\mathbf{w}\right)={\left(\mathbb{M}\left(\mathbf{w}\right)\right)}^{-1}\mathbf{F}\left({\mathbf{U}}^n\right)\in \mathbf{C}$. As ${\mathbf{w}}_k\to {\mathbf{w}}_0(k\to \infty )$ in ${\mathrm{R}}^{2\mathbb{N}}$, we know that $\mathbb{M}\left({\mathbf{w}}_k\right)\to \mathbb{M}\left({\mathbf{w}}_0\right)$. Since the inverse operator is continuous, we have ${\left(\mathbb{M}\left({\mathbf{w}}_k\right)\right)}^{-1}\to {\left(\mathbb{M}\left({\mathbf{w}}_0\right)\right)}^{-1}$. It follows that $M\left({\mathbf{w}}_k\right)\to M\left({\mathbf{w}}_0\right)$ in ${\mathrm{R}}^{2\mathbb{N}}$, so the map M is continuous.

By the Brouwer’s fixed point theorem [15,21], there exists a fixed point ${\mathbf{U}}^{n+1}\in \mathbf{C}$ such that ${\mathbf{U}}^{n+1}=M\left({\mathbf{U}}^{n+1}\right)$. Thus, ${\mathbf{U}}^n$ satisfies:

$$\begin{array}{c}\hfill \mathbb{M}\left({\mathbf{U}}^{n+1}\right){\mathbf{U}}^{n+1}=\mathbf{F}\left({\mathbf{U}}^n\right).\end{array}$$

□

5. Numerical Results

In this section, we present some numerical examples to investigate the accuracy and efficiency of the new positive-preserving finite volume scheme on distorted meshes.

The GMRES linear solver is used to solve the linear systems with stopping tolerance $1.0\times {10}^{-3}$. The stopping tolerance of nonlinear iteration is taken to be ${\epsilon }_{non}=1.0\times {10}^{-8}$. The discrete solution errors and flux errors are investigated in the discrete $L_2$ norms, which are defined by the following expressions:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & {\mathrm{E}}_u={\left(\sum _{K\in \mathcal{M}}|K|{(E\left({\mathit{x}}_K\right)-E_K)}^2\right)}^{\frac{1}{2}},\hfill \\ \hfill & {\mathrm{E}}_q={\left(\sum _{\sigma \in \mathcal{E}}{S}_{\sigma }((F_{\sigma }-F_{\sigma }^{ext})/{|\sigma |)}^2/\sum _{\sigma \in \mathcal{E}}{S}_{\sigma }\right)}^{\frac{1}{2}},\hfill \end{array}\end{array}$$

where $S_{\sigma }$ is the measure associated with $\sigma$, $F_{\sigma }$ is the numerical flux, and the analytical flux $F_{\sigma }^{ext}$ is evaluated by the midpoint rule. Besides, we define the $L_2$-norm of solution to compare the numerical results on different meshes. Let:

$$\begin{array}{c}\hfill L_2(E+T)=\sqrt{\sum _{K\in \mathcal{M}}{(E+T)}^{2}m\left(K\right)}.\end{array}$$

The energy conservation error is a criterion to judge the precision of the scheme. Define the total energy:

$$\begin{array}{c}\hfill E_{total}=\sum _{K\in \mathcal{M}}(E+T)m\left(K\right).\end{array}$$

The energy conservation error is $E_{total}^{error}=|E_{total}^N-E_{total}^0|$, where $E_{total}^0$ is the total energy at initial time, and $E_{total}^N$ is the total energy at final time. Besides, the following notations are used for the numerical tests

-

itn: average number of linear iterations;

-

nitn: average number of nonlinear Picard iterations;

-

$E_{min}$: minimal value of the numerical solution E;

-

$T_{min}$: minimal value of the numerical solution T.

5.1. Accuracy Test

In this test, we will consider the parabolic equation on the unit square $\mathsf{\Omega }={[0,1]}^2$ as:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \frac{\partial E}{\partial t}-\nabla ·(D\nabla E)=f\hfill \end{array}\end{array}$$

with the full Dirichlet boundary condition. The diffusion tensor D is taken as:

$$\begin{array}{c}\hfill D=\left\{\begin{array}{c}\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),x\le 0.5,\hfill \\ \left(\begin{array}{cc}10& 3\\ 3& 1\end{array}\right),x>0.5.\hfill \end{array}\right.\end{array}$$

The exact solution is chosen to be $E=t{E}_1(x,y)$:

$$\begin{array}{c}\hfill E_1(x,y)=\left\{\begin{array}{c}1-2y^2+4xy+6x+2y,x\le 0.5,\hfill \\ -2y^2+1.6xy-0.6x+3.2y+4.3,x>0.5.\hfill \end{array}\right.\end{array}$$

This problem is run up to time $t=0.1$ with the time step $\Delta t=0.2h^2$. This problem is solved on three type meshes as in Figure 3, such that the nodes of unstructured mesh located on line $x=0.5$ are distorted only in the y-direction. In Figure 4, the solution and flux errors are graphically depicted as log–log plots of the discrete $L_2$ norm errors versus the mesh size h. The results demonstrate that the scheme has second (respectively first)-order convergence rate with the solution (resp. flux).

5.2. Results without Flux Limiter

Mousseau and Knoll present an interesting two-dimensional multimaterial test problem that has a central blast wave moving out around two square obstacles. The Equations (1) and (2) are solved on the $64\times 64$ mesh grid as Figure 3. The background material uses $Z=1$, while the obstacles use $Z=10$. These obstacles are located at:

$$\begin{array}{c}\hfill x_1=0.3125\pm 0.125,y_1=0.6875\pm 0.125,\\ \hfill x_2=0.6875\pm 0.125,y_2=0.3125\pm 0.125.\end{array}$$

All four walls are insulated with respect to radiation diffusion and material conduction:

$$\frac{\partial E}{\partial x}{|}_{x=0}=\frac{\partial E}{\partial x}{|}_{x=1}=\frac{\partial E}{\partial x}{{|}_{y=0}=\frac{\partial E}{\partial x}|}_{y=1}=0,$$

and

$$\frac{\partial T}{\partial x}{|}_{x=0}=\frac{\partial T}{\partial x}{|}_{x=1}=\frac{\partial T}{\partial x}{{|}_{y=0}=\frac{\partial T}{\partial x}|}_{y=1}=0.$$

The initial energy distribution has a Gaussian peak near the origin:

$$E(x,y)=0.001+100exp\left[-100(x^2+y^2)\right].$$

The initial material temperature is taken to be $T(x,y)=E{(x,y)}^{1/4}$. The initial radiation spreads out and flows around the obstacles. This problem is simulated with the final time $t=1.5$ and time step $\Delta t=5.0\times {10}^{-4}$. The contour of radiation energy density and material temperature at time $t=1.5$ are presented in Figure 5 and Figure 6, respectively. We can see that the contours on random mesh and triangular mesh accord with that on uniform mesh. In

Table 1. Numerical results for problem without flux limiter.

Mesh	${\mathit{E}}_{\mathit{min}}$	${\mathit{T}}_{\mathit{min}}$	${\mathit{L}}_2$ norm	nitn	itn	${\mathit{E}}_{\mathit{total}}^{\mathit{error}}$
Mesh1	1.000 $\times {10}^{-3}$	0.1778	1.204	7.121	19.424	4.571 $\times {10}^{-3}$
Mesh2	9.999 $\times {10}^{-4}$	0.1778	1.207	7.391	19.567	5.435 $\times {10}^{-12}$
Mesh3	9.594 $\times {10}^{-4}$	0.1767	1.211	7.684	20.376	1.742 $\times {10}^{-12}$

, the minimal value and $L_2$-norm of numerical solution on different mesh are similar. Moreover, the minimal values of the numerical solution E and T are positive. The energy conservation error on different mesh is almost machine precision. The average number of nonlinear iterations in one time step and the average number of linear iterations per nonlinear iteration are also shown in Table 1. These numerical results indicate that our scheme is positivity-preserving, conservative, and robust in solving this problem.

The results of Anderson acceleration for the Picard iteration method are also presented in

Table 2. Numerical results of total number of nonlinear iterations.

Mesh	$\mathit{m}=1$	$\mathit{m}=2$	$\mathit{m}=3$	$\mathit{m}=5$	$\mathit{m}=7$	$\mathit{m}=9$
Mesh1	21,363	19,506	16,872	15,783	15,471	15,471
Mesh2	22,173	20,229	17,175	15,759	15,753	15,753
Mesh3	23,052	20,937	17,298	16,241	16,239	16,239

. The total number of nonlinear iterations is reduced significantly by the Anderson acceleration. The reasonable choice for this problem is $m=5$, since the total number of nonlinear iterative times is not observably decreasing for $m=7$ and $m=9$. From these numerical results, we know that the positivity-preserving scheme is robust and accurate in simulating this problem on distorted meshes.

5.3. Results with Flux Limiter

In this numerical example, we consider simulating the Equations (1) and (2) with flux limiter on the $64\times 64$ mesh as Figure 3. The initial and boundary conditions are given as the above numerical example. This problem is run up to $t=3$, and the time step is taken as $\Delta t=5.0$ $\times {10}^{-4}$.

The contours of radiation energy density and material temperature are shown in Figure 7 and Figure 8, respectively. The contours on unstructured meshes also accord with that on uniform mesh. Moreover, no spurious oscillation can be seen on the distorted meshes. In

Table 3. Numerical results for problem with flux limiter.

Mesh	${\mathit{E}}_{\mathit{min}}$	${\mathit{T}}_{\mathit{min}}$	${\mathit{L}}_2$ norm	nitn	itn	${\mathit{E}}_{\mathit{total}}^{\mathit{error}}$
Mesh1	1.000 $\times {10}^{-3}$	0.1778	1.169	7.427	18.671	1.065 $\times {10}^{-12}$
Mesh2	9.911 $\times {10}^{-4}$	0.1779	1.170	7.549	18.285	6.337 $\times {10}^{-12}$
Mesh3	9.947 $\times {10}^{-4}$	0.1772	1.168	8.370	21.054	7.436 $\times {10}^{-12}$

, some numerical results of the radiation energy density and material temperature are presented. The energy conservation errors of this problem are almost machine-precise. The $L_2$-norm and minimal value of the numerical solution on distorted meshes closely approximate that on uniform mesh. The average number of nonlinear iterations and linear iterations illustrate the positivity-preserving scheme is efficient. These numerical results of the scheme are very close to the numerical solution presented in [14,15,16].

The performance of the Anderson acceleration for the Picard iteration method is also considered. In

Table 4. Numerical results of total number of nonlinear iterations.

Mesh	$\mathit{m}=1$	$\mathit{m}=2$	$\mathit{m}=3$	$\mathit{m}=5$	$\mathit{m}=7$	$\mathit{m}=9$
Mesh1	44,562	42,876	37,170	33,390	33,372	33,368
Mesh2	45,294	43,932	38,652	36,006	35,151	35,144
Mesh3	50,220	47,484	42,222	39,361	37,281	37,275

, the total number of nonlinear iterations with the different values of m are presented, which shows the effectiveness of Anderson acceleration method. These results illustrate that $m=7$ is a better choice for this test.

6. Discussion

We presented a positivity-preserving finite volume scheme for nonequilibrium radiation diffusion systems. The advantages of this numerical scheme are the decomposition of the normal vector and the interpolation algorithm. Moreover, it has a fixed stencil and second-order accuracy on the distorted meshes. By applying the Anderson acceleration, the number of nonlinear iterations can clearly be reduced. Besides, numerical results illustrate that this positivity-preserving scheme is an accurate method in solving the nonequilibrium radiation diffusion equations. In the future, we hope to apply the positivity-preserving finite volume scheme to more practical problems, such as three temperature radiation diffusion problem, practical multimaterial physical problems, etc.