Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
Skip to main content

An Asymptotic-Preserving IMEX Method for Nonlinear Radiative Transfer Equation

  • Published:
Journal of Scientific Computing Aims and scope Submit manuscript

Abstract

We present an asymptotic preserving method for the radiative transfer equations in the framework of \(P_N\) method. An implicit and explicit numerical scheme is proposed to solve the \(P_N\) system based on the order analysis of the expansion coefficients of the specific intensity, where the order of each expansion coefficient is derived by the Chapman-Enskog method. The coefficients at higher-order are treated explicitly while those at lower-order are treated implicitly in each equation of the \(P_N\) system. Energy inequality is proved for this numerical scheme. Several numerical examples validate the efficiency of this scheme in both optically thick and thin regions.

This is a preview of subscription content, log in via an institution to check access.

Access this article

Subscribe and save

Springer+ Basic
$34.99 /Month
  • Get 10 units per month
  • Download Article/Chapter or eBook
  • 1 Unit = 1 Article or 1 Chapter
  • Cancel anytime
Subscribe now

Buy Now

Price excludes VAT (USA)
Tax calculation will be finalised during checkout.

Instant access to the full article PDF.

Fig. 1
Fig. 2
Fig. 3
Fig. 4
Fig. 5
Fig. 6
Fig. 7
Fig. 8
Fig. 9
Fig. 10
Fig. 11
Fig. 12

Similar content being viewed by others

Data Availability

Enquiries about data availability should be directed to the authors.

References

  1. Brunner, T. A.: Forms of approximate radiation transport. Sandia report, (2002)

  2. Crestetto, A., Crouseilles, N., Dimarco, G., Lemou, M.: Asymptotically complexity diminishing schemes (ACDS) for kinetic equations in the diffusive scaling. J. Comput. Phys. 394, 243–262 (2019)

    Article  MathSciNet  MATH  Google Scholar 

  3. Densmore, J.: Asymptotic analysis of the spatial discretization of radiation absorption and re-emission in Implicit Monte Carlo. J. Comput. Phys. 230(4), 1116–1133 (2011)

    Article  MathSciNet  MATH  Google Scholar 

  4. Densmore, J., Park, H., Wollaber, A., Rauenzahn, R., Knoll, D.: Monte Carlo simulation methods in moment-based scale-bridging algorithms for thermal radiative-transfer problems. J. Comput. Phys. 284, 40–58 (2015)

    Article  MathSciNet  MATH  Google Scholar 

  5. Di, Y., Fan, Y., Kou, Z., Li, R., Wang, Y.: Filtered hyperbolic moment method for the Vlasov equation. J. Sci. Comput. 79(2), 969–991 (2019)

    Article  MathSciNet  MATH  Google Scholar 

  6. Discacciati, N., Hesthaven, J., Deep, R.: Controlling oscillations in high-order Discontinuous Galerkin schemes using artificial viscosity tuned by neural networks. J. Comput. Phys. 409, 109304 (2020)

    Article  MathSciNet  MATH  Google Scholar 

  7. Fleck, J., Cummings, J.: An implicit Monte Carlo scheme for calculating time and frequency dependent nonlinear radiation transport. J. Comput. Phys. 8(3), 313–342 (1971)

    Article  MathSciNet  MATH  Google Scholar 

  8. Gentile, N.: Implicit Monte Carlo diffusion-an acceleration method for Monte Carlo time-dependent radiative transfer simulations. J. Comput. Phys. 172(2), 543–571 (2001)

    Article  MATH  Google Scholar 

  9. Hammer, H., Park, H., Chacón, L.: A multi-dimensional, moment-accelerated deterministic particle method for time-dependent, multi-frequency thermal radiative transfer problems. J. Comput. Phys. 386, 653–674 (2019)

    Article  MathSciNet  MATH  Google Scholar 

  10. Hou, T., Li, R.: Computing nearly singular solutions using pseudo-spectral methods. J. Comput. Phys. 226(1), 379–397 (2007)

    Article  MathSciNet  MATH  Google Scholar 

  11. Jang, J., Li, F., Qiu, J., Xiong, T.: Analysis of asymptotic preserving dg-imex schemes for linear kinetic transport equations in a diffusive scaling. SIAM J. Numer. Anal. 52(4), 2048–2072 (2014)

    Article  MathSciNet  MATH  Google Scholar 

  12. Jang, J., Li, F., Qiu, J., Xiong, T.: High order asymptotic preserving DG-IMEX schemes for discrete-velocity kinetic equations in a diffusive scaling. J. Comput. Phys. 281, 199–224 (2015)

    Article  MathSciNet  MATH  Google Scholar 

  13. Jin, S., Levermore, C.: The discrete-ordinate method in diffusive regimes. Transp. Theory Stat. Phys 20(1–2), 413–439 (1991)

    Article  MathSciNet  MATH  Google Scholar 

  14. Jin, S., Levermore, C.: Fully discrete numerical transfer in diffusive regimes. Transp. Theory Stat. Phys 22(6), 739–791 (1993)

    Article  MathSciNet  MATH  Google Scholar 

  15. Jin, S., Pareschi, L., Toscani, G.: Uniformly accurate diffusive relaxation schemes for multiscale transport equations. SIAM J. Numer. Anal. 38(3), 913–936 (2000)

    Article  MathSciNet  MATH  Google Scholar 

  16. Kershaw, D.: Flux limiting nature’s own way. Technical Report UCRL-78378, Lawrence Livermore National Laboratory, Livermore, CA, (1976)

  17. Klar, A.: An asymptotic-induced scheme for nonstationary transport equations in the diffusive limit. SIAM J. Numer. Anal 35(6), 1073–1094 (1998)

    Article  MathSciNet  MATH  Google Scholar 

  18. Koch, R., Krebs, W., Wittig, S., Viskanta, R.: Discrete ordinates quadrature schemes for multidimensional radiative transfer. J. Quant. Spectrosc. Ra. 53(4), 353–372 (1995)

    Article  Google Scholar 

  19. Los Alamos National Laboratory. An implicit Monte Carlo code for thermal radiative transfer: Capabilities, development, and usag. LA-14195-MS, 2000

  20. Laboure, V., McClarren, R., Hauck, C.: Implicit filtered \(P_N\) for high-energy density thermal radiation transport using discontinuous galerkin finite elements. J. Comput. Phys. 321, 624–643 (2016)

    Article  MathSciNet  MATH  Google Scholar 

  21. Laiu, M., Frank, M., Hauck, C.: A positive asymptotic-preserving scheme for linear kinetic transport equations. SIAM J. Sci. Comput. 41, A1500–A1526 (2019)

    Article  MathSciNet  MATH  Google Scholar 

  22. Larsen, A., Morel, J.: Asymptotic solutions of numerical transport problems in optically thick, diffusive regimes. J. Comput. Phys. 69(2), 283–324 (1987)

    Article  MathSciNet  MATH  Google Scholar 

  23. Larsen, A., Morel, J.: Asymptotic solutions of numerical transport problems in optically thick, diffusive regimes ii. J. Comput. Phys. 83(1), 212–236 (1989)

    Article  MathSciNet  MATH  Google Scholar 

  24. Larsen, E., Kumar, A., Morel, J.: Properties of the implicitly time-differenced equations of thermal radiation transport. J. Comput. Phys. 238, 82–96 (2013)

    Article  MathSciNet  MATH  Google Scholar 

  25. Lathrop, K., Garlson, B.: Discrete ordinates angular quadrature of the neutron transport equation. Los Alamos Scientific Laboratory, (1965)

  26. Lemou, E., Mieussens, L.: A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit. SIAM J. Sci. Comput. 31, 334–368 (2010)

    Article  MathSciNet  MATH  Google Scholar 

  27. Lewis, E., Miller, W.: Computational Methods in Neutron Transport. United States, (1993)

  28. Li, W., Liu, C., Zhu, Y., Zhang, J., Xu, K.: Unified gas-kinetic wave-particle methods iii: Multiscale photon transport. J. Comput. Phys. 408, 109280 (2020)

    Article  MathSciNet  MATH  Google Scholar 

  29. Maginot, P., Ragusa, J., Morel, J.: High-order solution methods for grey discrete ordinates thermal radiative transfer. J. Comput. Phys. 327, 719–746 (2016)

    Article  MathSciNet  MATH  Google Scholar 

  30. Mathews, K.: On the propagation of rays in discrete ordinates. Nucl. Sci. Eng. 132, 155–180 (1999)

    Article  Google Scholar 

  31. McClarren, R., Evans, T., Lowrie, R., Densmore, J.: Semi-implicit time integration for PN thermal radiative transfer. J. Comput. Phys. 227(16), 7561–7586 (2008)

    Article  MathSciNet  MATH  Google Scholar 

  32. McClarren, R., Hauck, C.: Robust and accurate filtered spherical harmonics expansions for radiative transfer. J. Comput. Phys. 229(16), 5597–5614 (2010)

    Article  MathSciNet  MATH  Google Scholar 

  33. McClarren, R., Hauck, C.: Simulating radiative transfer with filtered spherical harmonics. Phys. Lett. A 374(22), 2290–2296 (2010)

    Article  MATH  Google Scholar 

  34. McClarren, R., Holloway, J., Brunner, T.: On solutions to the \(P_n\) equations for thermal radiative transfer. J. Comput. Phys. 227(5), 2864–2885 (2008)

    Article  MathSciNet  MATH  Google Scholar 

  35. Mieussens, L.: On the asymptotic preserving property of the unified gas kinetic scheme for the diffusion limit of linear kinetic models. J. Comput. Phys. 253, 138–156 (2013)

    Article  MathSciNet  MATH  Google Scholar 

  36. Morel, J., Wareing, T., Lowrie, R., Parsons, D.: Analysis of ray-effect mitigation techniques. Nucl. Sci. Eng. 144, 1–22 (2003)

    Article  Google Scholar 

  37. Olson, G.: Second-order time evolution of \(P_N\) equations for radiation transport. J. Comput. Phys. 228(8), 3072–3083 (2009)

    Article  MathSciNet  MATH  Google Scholar 

  38. Park, H., Knoll, D., Rauenzahn, R., Wollaber, A., Densmore, J.: A consistent, moment-based, multiscale solution approach for thermal radiative transfer problems. Transp. Theory Stat. Phys. 41(3–4), 284–303 (2012)

    Article  MATH  Google Scholar 

  39. Peng, Z., Cheng, Y., Qiu, J., Li, F.: Stability-enhanced AP IMEX-LDG schemes for linear kinetic transport equations under a diffusive scaling. J. Comput. Phys. 415, 109485 (2020)

    Article  MathSciNet  MATH  Google Scholar 

  40. Seibold, B., Frank, M.: Starmap code. website. http://www.math.temple.edu/~seibold/research/starmap

  41. Seibold, B., Frank, M.: Starmap-a second order staggered grid method for spherical harmonics moment equations of radiative transfer. ACM T. Math. Software (TOMS) 41(1), 4 (2014)

    MathSciNet  MATH  Google Scholar 

  42. Shi, Y., Song, P., Sun, W.: An asymptotic preserving unified gas kinetic particle method for radiative transfer equations. J. Comput. Phys. 420, 109687 (2020)

    Article  MathSciNet  MATH  Google Scholar 

  43. Sun, W., Jiang, S., Xu, K.: An asymptotic preserving unified gas kinetic scheme for gray radiative transfer equations. J. Comput. Phys. 285(15), 265–279 (2015)

    Article  MathSciNet  MATH  Google Scholar 

  44. Sun, W., Jiang, S., Xu, K.: An asymptotic preserving implicit unified gas kinetic scheme for frequency-dependent radiative transfer equations. Int. J. Numer. Anal. Mod. 15(1–2), 134–153 (2018)

    MathSciNet  MATH  Google Scholar 

  45. Tang, M., Wang, L., Zhang, X.: Accurate front capturing asymptotic preserving scheme for nonlinear gray radiative transfer equation. SIAM J. Sci. Comput. 43(3), B759–B783 (2021)

    Article  MathSciNet  MATH  Google Scholar 

  46. Warsa, J., Wareing, T., Morel, J.: Krylov iterative methods and the degraded effectiveness of diffusion synthetic acceleration for multidimensional \(S_N\) calculations in problems with material discontinuities. Nucl. Sci. Eng. 147, 218–248 (2004)

    Article  Google Scholar 

  47. Xiong, T., Jang, J., Li, F., Qiu, J.: High order asymptotic preserving nodal discontinuous Galerkin IMEX schemes for the BGK equation. J. Comput. Phys. 284, 70–94 (2015)

    Article  MathSciNet  MATH  Google Scholar 

  48. Xiong, T., Sun, W., Shi, Y., Song, P.: High order asymptotic preserving discontinuous Galerkin methods for gray radiative transfer equations. arXiv:2011.14090, (2020)

  49. Yan, B., Jin, S.: A successive penalty-based asymptotic-preserving scheme for kinetic equations. SIAM J. Sci. Comput. 35(1), A150–A172 (2013)

    Article  MathSciNet  MATH  Google Scholar 

Download references

Acknowledgements

We thank Prof. Ruo Li from PKU, Prof. Zhenning Cai from NUS, Prof. Tao Xiong from XMU, Prof. Kailiang Wu from SUSTech, Dr. Zhichao Peng from MSU and Prof. Jiequan Li, Prof. Wenjun Sun, Dr. Yi Shi from IAPCM for their valuable suggestions.

Funding

Weiming Li is partially supported by the Science Challenge Project (No. TZ2016002) and the National Natural Science Foundation of China (12001051). Peng Song is partially supported by the Science Challenge Project (No. TZ2016002), the CAEP foundation (No.CX20200026). The work of Yanli Wang is partially supported by Science Challenge Project (No. TZ2016002) and the National Natural Science Foundation of China (Grant No. 12171026, U1930402 and 12031013).

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Weiming Li.

Ethics declarations

Conflict of interest

The authors have not disclosed any competing interest.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Appendix

Appendix

1.1 The Gray Approximation of the Radiative Transfer Equations for 1D Angle Problem and Related \(P_N\) Equations

The time-dependent gray approximation of the radiative transfer equations [31] in a one-dimensional planar geometry medium has the form as

$$\begin{aligned}&\frac{\epsilon ^2}{c} \dfrac{\partial {I}}{\partial {t}} +\epsilon \mu \dfrac{\partial {I}}{\partial {x}} = \sigma \left( \frac{1}{2}a c T^4 - I\right) , \qquad x \in [0, L], \end{aligned}$$
(A.1a)
$$\begin{aligned}&\epsilon ^2 C_{v} \dfrac{\partial {T}}{\partial {t}} = \sigma \left( \int _{-1}^1 I \,\mathrm {d}\mu - ac T^4\right) , \end{aligned}$$
(A.1b)

where \(I(t, x, \mu )\) is the specific intensity of radiation, \(\mu = \cos \theta \in [-1, 1]\) is the internal coordinate associated with the angle \(\theta \in [0, \pi ]\). T(t, x) is the material temperature, and \(\sigma \) is the absorption opacity. Moreover, the one-dimensional form of (2.10) is reduced into

$$\begin{aligned} \frac{\epsilon ^2}{c} \dfrac{\partial {I}}{\partial {t}} +\epsilon \mu \dfrac{\partial {I}}{\partial {x}} = \sigma \left( \frac{1}{2} \int I \,\mathrm {d}\mu - I\right) , \qquad x \in [0, L]. \end{aligned}$$
(A.2)

For (A.1), the basis function for the \(P_N\) method is the Legendre polynomials. The moments are defined as

$$\begin{aligned} I_l = \int _{-1}^1 P_l(\mu ) I(t, x, \mu ) \,\mathrm {d}\mu , \qquad l = 0, \cdots M, \end{aligned}$$
(A.3)

where \(P_l\) is the Legendre polynomial. Then, the \(P_N\) equations for (A.1) are

$$\begin{aligned} \begin{aligned}&\frac{\epsilon ^2}{c}\dfrac{\partial {\bar{\mathbf{I}}}}{\partial {t}} + \epsilon \mathbf{B}^{\mathrm{low }} \dfrac{\partial {\bar{\mathbf{I}}}}{\partial {x}} + \epsilon {\mathbf{B}}^{\mathrm{up}} \dfrac{\partial {{\bar{\mathbf{I}}}}}{\partial {x}} = -\sigma {\bar{\mathbf{I}}} + \sigma a c T^4 e_1, \\&\epsilon ^2 C_{v} \dfrac{\partial {T}}{\partial {t}} = \sigma \left( I_0 - ac T^4\right) , \end{aligned} \end{aligned}$$
(A.4)

where \(\bar{\mathbf{I}}= (I_0, I_1, \cdots , I_M)\) and \(e_1 = (1, 0, \cdots , 0)^T\). \(\mathbf{B}^{\mathrm{low}}\) and \(\mathbf{B}^{\mathrm{up}}\) are triangular matrix with the non-zero entries as

$$\begin{aligned} \mathbf{B}^{\mathrm{low}}(i+1, i) = \frac{i}{2 i+1}, \qquad \mathbf{B}^{\mathrm{up}}(i, i+1) = \frac{i}{2(i-1)+1}, \qquad i = 1, \cdots M. \end{aligned}$$
(A.5)

1.2 Inflow Boundary Condition for PN Equations

We implement the inflow boundary condition for the \(P_N\) equations by specifying the values of coefficients of the \(P_N\) system in ghost cells. Choosing the left boundary as an example, the incoming specific intensity incident on the boundary interface is

$$\begin{aligned} I(\mu ) = I^b(\mu ), \quad \text {for}~\mu > 0. \end{aligned}$$
(A.6)

For the \(P_N\) method, the numerical boundary can be rewritten as

$$\begin{aligned} I^{\mathrm{ghost}}(\mu ) = \left\{ \begin{array}{cc} I^{b}(\mu ), &{}\mu > 0, \\ I^{i}(\mu ), &{} \mu < 0, \end{array} \right. \end{aligned}$$
(A.7)

where \(I^{i}(\mu )\) is the specific intensity at the left boundary of the area. Then, the expansion coefficient at the ghost cell is

$$\begin{aligned} I^{\mathrm{ghost}}_{l} = \int _{-1}^1 I^{\mathrm{ghost}}(\mu ) P_l(\mu )\,\mathrm {d}\mu . \end{aligned}$$
(A.8)

The implementation of the inflow boundary condition in 2D is similar in spirit to that of 1D. Supposing \({\varvec{n}}\) is the outward normal of the boundary interface, the incident specific intensity on the boundary is

$$\begin{aligned} I({\varvec{\Omega }}) = I^b({\varvec{\Omega }}), \quad \text {for}~{\varvec{\Omega }}\cdot {\varvec{n}}< 0. \end{aligned}$$
(A.9)

For the \(P_N\) method, the numerical boundary can be rewritten as

$$\begin{aligned} I^{\mathrm{ghost}}({\varvec{\Omega }}) = \left\{ \begin{array}{cc} I^{b}({\varvec{\Omega }}), &{} {\varvec{\Omega }}\cdot {\varvec{n}}< 0, \\ I^{i}({\varvec{\Omega }}), &{} {\varvec{\Omega }}\cdot {\varvec{n}}> 0, \end{array} \right. \end{aligned}$$
(A.10)

where \(I^{i}({\varvec{\Omega }})\) is the specific intensity on the interior side of the boundary interface. Thus, the expansion coefficient at the ghost cell is

$$\begin{aligned} \begin{aligned} I^{m, \mathrm ghost}_{l}&= \int _{{\mathbb {S}}^2} I^{\mathrm{ghost}}({\varvec{\Omega }}) \overline{Y^m_l}({\varvec{\Omega }}) \,\mathrm {d}{\varvec{\Omega }}\\&= \int _{{\varvec{\Omega }}\cdot {\varvec{n}}< 0} I^{b}({\varvec{\Omega }}) \overline{Y^m_l}({\varvec{\Omega }}) \,\mathrm {d}{\varvec{\Omega }}+ \int _{{\varvec{\Omega }}\cdot {\varvec{n}}> 0} I^{i}({\varvec{\Omega }}) \overline{Y^m_l}({\varvec{\Omega }}) \,\mathrm {d}{\varvec{\Omega }}\\&= \int _{{\varvec{\Omega }}\cdot {\varvec{n}}< 0} I^{b}({\varvec{\Omega }}) \overline{Y^m_l}({\varvec{\Omega }}) \,\mathrm {d}{\varvec{\Omega }}+ \sum \limits _{j=0}^M\sum \limits _{k=-j}^j I^{k, i}_j \int _{{\varvec{\Omega }}\cdot {\varvec{n}}> 0} Y^k_j({\varvec{\Omega }}) \overline{Y^m_l}({\varvec{\Omega }}) \,\mathrm {d}{\varvec{\Omega }}. \end{aligned} \end{aligned}$$
(A.11)

The integration \(\int _{{\varvec{\Omega }}\cdot {\varvec{n}}> 0} Y^k_j({\varvec{\Omega }}) \overline{Y^m_l}({\varvec{\Omega }}) \,\mathrm {d}{\varvec{\Omega }}\) does not depend on the specific numerical solutions, and is pre-computed.

Fig. 13
figure 13

The stability region for the numerical scheme (3.18) of \(P_1\) system under the condition (A.12). The x-axis is \(\beta _2 = \log _{10}(\epsilon / \Delta x)\), and the y-axis is the CFL number C. The blue region is the area where the numerical scheme is stable and the yellow region is the area where the numerical scheme is unstable

1.3 Proof of Proposition 2

In this section, the proof of the Proposition 2 is proposed here.

Proof

Following the method in [39], we will begin the Fourier analysis of (3.18) for the \(P_1\) system of the linear equation system (A.2). The result can be extended to the generalized \(P_N\) system naturally. We will first discuss two special cases where \(\xi = 0, \pi \). Therein, \(\mathbf{C}\) is reduced into a real diagonal matrix with the maximum eigenvalues equaling 1. According to the principle, the numerical scheme is stable.

Then, we study the general case by considering two scenarios according to the time step length:

  1. 1.

    \(\epsilon < \Delta x\), where

    $$\begin{aligned} \Delta t = C \Delta x^2. \end{aligned}$$
    (A.12)

    Substituting the time step length (A.12) into (3.36), we can find that \(\lambda _i, i=1,2\) are functions of C, \(\frac{\epsilon }{\Delta x}\), \(\alpha \) and \(\xi \). Introducing two variables as \(\beta _1 = \log _{10}(C)\) and \(\beta _2 = \log _{10} (\epsilon / \Delta x)\), the stability regions are plotted in Fig. 13 with fixed \(\alpha \). Here the discrete wave number \(\xi \) is uniformly taken from \([0, 2\pi ]\) with 200 samples. In this case, due to the definition of C which is the CFL number and \(\epsilon < \Delta x\), the range for \(\beta _1\) and \(\beta _2\) is changed into

    $$\begin{aligned} \beta _1< 0, \qquad \beta _2 < 0. \end{aligned}$$
    (A.13)

    Moreover, it is natural to demand that \(\epsilon< \Delta x < 0.4\). Thus, \(\alpha \) is taken uniformly from \([0, \exp (-1/0.16)]\) with 100 samples and six cases are shown in Fig. 13 due to their similar behavior. From Fig. 13, we can find that when \(\alpha = 0\), the numerical scheme is always stable. However, with the increase of \(\alpha \), the stability region is becoming smaller, especially when the CFL number C is large and the radio \(\epsilon / \Delta x\) is small. We find that when \(\alpha = \exp (-1/0.16)\), the numerical scheme is stable when \(\log _{10}(\epsilon / \Delta x) > -2.5\). Noting that when \(\alpha = \exp (-1/0.16)\), which means \(\epsilon = 0.4\), \(\log _{10}(\epsilon / \Delta x)\) is always larger than \(-2.5\) for \(\Delta x < 1\). This indicates that in the simulation of benchmark problems, the stability condition is always satisfied.

  2. 2.

    \(\epsilon > \Delta x\), where

    $$\begin{aligned} \Delta t = C \epsilon \Delta x. \end{aligned}$$
    (A.14)

    Substituting the time step length (A.14) into (3.36), we can easily find that \(\lambda _i, i = 1, 2\) are also the function of C, \(\frac{\epsilon }{\Delta x}\), \(\alpha \) and \(\xi \). Introducing the same two variables \(\beta _{i}, i = 1,2\), we plot the stability regions in Fig. 14 with fixed \(\alpha \). Here the discrete wave number \(\xi \) is uniformly taken from \([0, 2\pi ]\) with 200 samples. Since

    $$\begin{aligned} \alpha = \exp \left( -\frac{1}{\epsilon ^2}\right) , \end{aligned}$$
    (A.15)

    and assuming \(\epsilon < 1\) in the numerical test, \(\alpha \) is taken uniformly from [0, 0.5] with 100 samples. As their behavior is similar, the six cases \(\alpha = 0, 0.05, 0.1, 0. 2, 0.3\) and 0.5 are plotted here to illustrate the result. Moreover, due to the definition of C which is the CFL number, and the condition that \(\epsilon > \Delta x\), it holds that

    $$\begin{aligned} \beta _1 < 0, \qquad \beta _2 > 0. \end{aligned}$$
    (A.16)

    From Fig. 14, we can find that the numerical scheme is stable under the time step length (A.14). In the numerical tests, the upper bound of \(\epsilon \) is \(\epsilon = 10^5 \Delta x\), which is large enough for the computational parameter.

\(\square \)

Fig. 14
figure 14

The stability region for the numerical scheme (3.18) of \(P_1\) system under the condition (A.14). The x-axis is \(\beta _2 = \log _{10}(\epsilon / \Delta x)\), and the y-axis is the CFL number C. The blue region is the area where the numerical scheme is stable

1.4 Proof of Theorem 2

In this section, the proof of Theorem 2 is proposed here.

Proof

We will take \(M=2\) as an example, and it could be extended to the general case naturally. Moreover, without loss of generality, we set \(a = c = C_v = \sigma = 1\) in the proof. When \(M=2\), (4.10) is reduced into

$$\begin{aligned}&\epsilon ^2 \frac{I_{0,j}^{n+1} - I_{0,j}^n}{\Delta t} + \epsilon \frac{I_{1, j+1}^{n} - I_{1, j-1}^{n}}{2 \Delta x} - \frac{\alpha \epsilon }{2} \frac{I_{0, j+1}^n - 2 I_{0,j}^n + I_{0, j-1}^n}{\Delta x} = \left( (T_j^4)^{n+1} - I_{0,j}^{n+1}\right) , \end{aligned}$$
(A.17a)
$$\begin{aligned}&\epsilon ^2 \frac{I_{1,j}^{n+1} - I_{1,j}^n}{\Delta t} + \frac{\epsilon }{3} \frac{I_{0, j+1}^{n+1} - I_{0, j-1}^{n+1}}{2 \Delta x} +\frac{2\epsilon }{3} \frac{I_{2, j+1}^{n} - I_{2, j-1}^{n}}{2 \Delta x} - \frac{\alpha \epsilon }{2} \frac{I_{1, j+1}^n - 2 I_{1,j}^n + I_{1, j-1}^n}{\Delta x} \nonumber \\&\quad = - I_{1,j}^{n+1}, \end{aligned}$$
(A.17b)
$$\begin{aligned}&\epsilon ^2 \frac{I_{2,j}^{n+1} - I_{2,j}^n}{\Delta t} + \frac{2\epsilon }{5} \frac{I_{1, j+1}^{n+1} - I_{1, j-1}^{n+1}}{2 \Delta x} - \frac{\alpha \epsilon }{2} \frac{I_{2, j+1}^n - 2 I_{2,j}^n + I_{2, j-1}^n}{\Delta x} = - I_{2,j}^{n+1}, \end{aligned}$$
(A.17c)
$$\begin{aligned}&\epsilon ^2 \frac{T_j^{n+1} - T_j^n}{\Delta t} +\epsilon \frac{I_{0,j}^{n+1} - I_{0,j}^n}{\Delta t} + \epsilon \frac{I_{1, j+1}^{n} - I_{1, j-1}^{n}}{2 \Delta x} - \frac{\alpha \epsilon }{2} \frac{I_{0, j+1}^n - 2 I_{0,j}^n + I_{0, j-1}^n}{\Delta x} = 0. \end{aligned}$$
(A.17d)

For (A.17a), multiplying it by \(I_{0,j}^{n+1}\), we can get that

$$\begin{aligned}&\epsilon ^2 I_{0,j}^{n+1} \frac{I_{0,j}^{n+1} - I_{0,j}^n}{\Delta t} + \epsilon I_{0,j}^{n+1} \frac{I_{1, j+1}^{n} - I_{1, j-1}^{n}}{2 \Delta x} - \frac{\alpha \epsilon }{2} I_{0,j}^{n+1} \frac{I_{0, j+1}^n - 2 I_{0,j}^n + I_{0, j-1}^n}{\Delta x} \nonumber \\&\quad = I_{0,j}^{n+1} \left( (T_j^4)^{n+1} - I_{0,j}^{n+1}\right) . \end{aligned}$$
(A.18)

For (A.17b) and (A.17c), shifting it backward one time step and multiplying \(3I_{1,j}^{n}\) and \(5I_{2, j}^{n}\) respectively, we can derive that

$$\begin{aligned} \begin{aligned}&3\epsilon ^2 I_{1,j}^n \frac{I_{1,j}^{n} - I_{1,j}^{n-1}}{\Delta t} + \epsilon I_{1,j}^n \frac{I_{0, j+1}^{n} - I_{0, j-1}^{n}}{2 \Delta x} + \epsilon I_{1,j}^n \frac{I_{2, j+1}^{n-1} - I_{2, j-1}^{n-1}}{ \Delta x}\\&\quad - \frac{3\alpha \epsilon }{2} I_{1,j}^n\frac{I_{1, j+1}^{n-1} - 2 I_{1,j}^{n-1} + I_{1, j-1}^{n-1}}{\Delta x} = - 3(I_{1,j}^{n})^2, \\&5 \epsilon ^2 I_{2,j}^n\frac{I_{2,j}^{n} - I_{2,j}^{n-1}}{\Delta t} + \epsilon I_{2,j}^n \frac{I_{1, j+1}^{n} - I_{1, j-1}^{n}}{ \Delta x} \\&\quad -\frac{5 \alpha \epsilon }{2} I_{2,j}^n \frac{I_{2, j+1}^{n-1} - 2 I_{2,j}^{n-1} + I_{2, j-1}^{n-1}}{\Delta x} = - 5 (I_{2,j}^{n})^2. \end{aligned} \end{aligned}$$
(A.19)

Summing (A.18) and (A.19) over j, then it holds that

$$\begin{aligned}&\frac{\epsilon ^2}{2\Delta t} \sum _j\Big [ (I_{0,j}^{n+1})^2 - (I_{0,j}^{n})^2 + 3\left( (I_{1,j}^{n})^2 - (I_{1,j}^{n-1})^2 \right) + 5 \left( (I_{2,j}^{n})^2 - (I_{2,j}^{n-1})^2\right) \Big ]\nonumber \\&\quad + A_0 = A_1 + A_2 + A_3, \end{aligned}$$
(A.20)

where

$$\begin{aligned} A_0&= \frac{\epsilon ^2}{2\Delta t} \sum _j\Big [ (I_{0,j}^{n+1}- I_{0,j}^{n})^2 + 3(I_{1,j}^{n} - I_{1,j}^{n-1})^2 + 5(I_{2,j}^{n} - I_{2,j}^{n-1})^2 \Big ], \end{aligned}$$
(A.21a)
$$\begin{aligned} A_1&= -\frac{\epsilon }{2\Delta x}\sum _j \Big [ I_{0,j}^{n+1} (I_{1, j+1}^{n} - I_{1, j-1}^{n}) + I_{1,j}^n (I_{0, j+1}^{n} - I_{0, j-1}^{n}) + 2I_{1,j}^n (I_{2, j+1}^{n-1} - I_{2, j-1}^{n-1}) \nonumber \\&\quad + 2I_{2,j}^n (I_{1, j+1}^{n} - I_{1, j-1}^{n})\Big ], \end{aligned}$$
(A.21b)
$$\begin{aligned} A_2&= \frac{\alpha \epsilon }{2\Delta x}\sum _j \Big [ I_{0,j}^{n+1} (I_{0, j+1}^{n} - 2I_{0, j}^{n} +I_{0, j-1}^{n}) + 3I_{1,j}^n (I_{1, j+1}^{n-1} - 2 I_{1,j}^{n-1} + I_{1, j-1}^{n-1})\end{aligned}$$
(A.21c)
$$\begin{aligned}&\quad + 5I_{2,j}^n (I_{2, j+1}^{n-1} - 2I_{2, j}^{n-1} + I_{2, j-1}^{n-1})\Big ],\nonumber \\ A_3&= -\sum _j \Big [-I_{0,j}^{n+1}(T_{j}^4)^{n+1} + (I_{0,j}^{n+1})^2 + 3(I_{1, j}^n)^2 + 5(I_{2,j}^n)^2\Big ]. \end{aligned}$$
(A.21d)

Then we will begin from the approximation of \(A_1\) and \(A_2\). With some arrangement and the periodic boundary condition, \(A_1\) is changed into

$$\begin{aligned} \begin{aligned} A_1 =&-\frac{\epsilon }{2\Delta x}\sum _j \Big [ (I_{0,j}^{n+1} - I_{0,j}^{n}) (I_{1, j+1}^{n} - I_{1, j-1}^{n}) + 2(I_{1,j+1}^n - I_{1,j-1}^{n}) (I_{2, j}^{n} - I_{2, j}^{n-1}) \Big ]\\ \le&\frac{\epsilon }{2\Delta x}\sum _j \Big [ \frac{1}{2}\beta _1^2 (I_{0,j}^{n+1} - I_{0,j}^{n})^2 + \frac{1}{2 \beta _1^2} (I_{1, j+1}^{n} - I_{1, j-1}^{n})^2 + \beta _2^2 (I_{2,j}^n - I_{2,j}^{n-1})^2 \\&+ \frac{1}{\beta _2^2} (I_{1, j+1}^{n} - I_{1, j-1}^{n})^2 \Big ] \\ \le&\frac{\epsilon }{2\Delta x}\sum _j \Big [ \frac{1}{2}\beta _1^2 (I_{0,j}^{n+1} - I_{0,j}^{n})^2 + \beta _2^2 (I_{2,j}^n - I_{2,j}^{n-1})^2 + \left( \frac{2}{\beta _1^2} + \frac{4}{\beta _2^2}\right) (I_{1, j}^{n})^2 \Big ]. \end{aligned} \end{aligned}$$
(A.22)

With the estimation that

$$\begin{aligned} \begin{aligned}&\sum _j I_{0,j}^{n+1}(I_{0,j+1}^n - 2 I_{0,j}^n + I_{0,j-1}^n)\\&\qquad = \sum _j \Big [(I_{0,j}^{n+1} - I_{0,j}^n)(I_{0,j+1}^n - 2 I_{0,j}^n + I_{0,j-1}^n) + I_{0,j}^n(I_{0,j+1}^n - 2 I_{0,j}^n + I_{0,j-1}^n) \Big ] \\&\qquad \le \sum _j\Big [\frac{1}{2}\beta _3^2 (I_{0,j}^{n+1} - I_{0,j}^n)^2 + \frac{2}{\beta _3^2}(I_{0,j}^n - I_{0,j-1}^n)^2\Big ] - \sum _j (I_{0,j}^n - I_{0,j-1}^n)^2 \\&\qquad = \sum _j\Big [\frac{1}{2}\beta _3^2 (I_{0,j}^{n+1} - I_{0,j}^n)^2 + (\frac{2}{\beta _3^2} - 1)\left( I_{0,j}^n - I_{0,j-1}^n\right) ^2\Big ]. \end{aligned} \end{aligned}$$
(A.23)

Similarly, it also holds that

$$\begin{aligned} \begin{aligned}&\sum _j I_{1,j}^{n}(I_{1,j+1}^{n-1} - 2 I_{1,j}^{n-1} + I_{1,j-1}^{n-1}) = \sum _j I_{1,j}^{n-1}(I_{1,j+1}^{n} - 2 I_{1,j}^{n} + I_{1,j-1}^{n}) \\&\le \sum _j\Big [\frac{1}{2}\beta _4^2 (I_{1,j}^{n} - I_{1,j}^{n-1})^2 + 4\left( \frac{2}{\beta _4^2} - 1\right) (I_{1,j}^n)^2\Big ], \\&\sum _j I_{2,j}^{n}(I_{2,j+1}^{n-1} - 2 I_{2,j}^{n-1} + I_{2,j-1}^{n-1}) = \sum _j I_{2,j}^{n-1}(I_{2,j+1}^{n} - 2 I_{2,j}^{n} + I_{2,j-1}^{n}) \\&\le \sum _j\Big [\frac{1}{2}\beta _5^2 (I_{2,j}^{n} - I_{2,j}^{n-1})^2 + \left( \frac{2}{\beta _5^2}-1\right) (I_{2,j}^n - I_{2,j-1}^n)^2\Big ]. \end{aligned} \end{aligned}$$
(A.24)

Let

$$\begin{aligned} \beta _1^2 = \frac{2 \epsilon \Delta x}{\Delta t}- 2 \alpha , \qquad \beta _2^2 = \frac{5\epsilon \Delta x}{\Delta t}- 5 \alpha , \qquad \beta _3^2 = 2, \qquad \beta _4^2 = \frac{2 \Delta x \epsilon }{\alpha \Delta t}, \qquad \beta _5^2 = 2,\nonumber \\ \end{aligned}$$
(A.25)

then together with (A.21), (A.22), (A.23), (A.24) and (A.25), (A.20) is reduced into

$$\begin{aligned} \begin{aligned} \frac{\epsilon ^2}{2} \sum _j\Big [ (I_{0,j}^{n+1})^2 - (I_{0,j}^{n})^2 + 3\left( (I_{1,j}^{n})^2 - (I_{1,j}^{n-1})^2 \right) + 5 \left( (I_{2,j}^{n})^2 - (I_{2,j}^{n-1})^2\right) \Big ]\\ \le \beta _6\sum _j (I_{1,j}^{n})^2 + \Delta t \sum _j \Big [I_{0,j}^{n+1}(T_{j}^4)^{n+1} - (I_{0,j}^{n+1})^2 - 5(I_{2,j}^{n})^2\Big ], \end{aligned} \end{aligned}$$
(A.26)

with

$$\begin{aligned} \beta _6 = \frac{9}{10}\dfrac{\epsilon \Delta t}{\Delta x} \left( \frac{1}{\frac{\epsilon \Delta x}{\Delta t} - \alpha }\right) + \dfrac{6 \alpha ^2 (\Delta t)^2}{(\Delta x)^2}-\dfrac{6\alpha \epsilon \Delta t}{\Delta x} - 3 \Delta t. \end{aligned}$$
(A.27)

If it holds for \(\beta _6\) that

$$\begin{aligned} \beta _6 \le 0, \end{aligned}$$
(A.28)

with the time step length (3.32), then we can derive the stability result (4.16). Precisely, with (A.17a) and (A.17d), we can derive that

$$\begin{aligned} \epsilon ^2 \frac{T_j^{n+1} - T_j^n}{\Delta t} = -(T_j^4)^{n+1} + I_{0,j}^{n+1}. \end{aligned}$$
(A.29)

Multiplying (A.29) with \((T_{j}^4)^{n+1}\) and summing over j, it holds with (A.26)

$$\begin{aligned} \begin{aligned}&\sum _j \left[ \frac{\epsilon ^2}{2\Delta t}\Big ( (I_{0,j}^{n+1})^2 - (I_{0,j}^{n})^2 + 3\left[ (I_{1,j}^{n})^2 - (I_{1,j}^{n-1})^2 \right] + 5 \left[ (I_{2,j}^{n})^2 - (I_{2,j}^{n-1})^2\right] \Big ) \right. \\&\qquad \left. +\frac{\epsilon ^2}{5 \Delta t} \Big [(T_j^5)^{n+1} - (T_j^5)^{n}\Big ] \right] \le - \sum _j \left[ I_{0,j}^{n+1} - (T_{j}^4)^{n+1}\right] ^2 \le 0. \end{aligned} \end{aligned}$$
(A.30)

We derive the energy stability (4.16). The only point left is to prove (A.28), which we will be done in two cases:

  1. 1.

    \(\epsilon > \Delta x\), in which case,

    $$\begin{aligned} \Delta t= C \epsilon \Delta x. \end{aligned}$$
    (A.31)

    Substituting (A.31) into (A.27), we can deduce that

    $$\begin{aligned} \beta _6 = \frac{9C^2\epsilon ^2}{10}\left( \frac{1}{1 - \alpha C }\right) + 6\alpha \epsilon ^2C(\alpha C - 1) - 3C \Delta x \epsilon . \end{aligned}$$
    (A.32)

    Thus if

    $$\begin{aligned} 0< C < \min \left( \frac{\epsilon }{\alpha \Delta x},\frac{10 \Delta x}{3 \epsilon + 10 \Delta x \alpha } \right) , \end{aligned}$$
    (A.33)

    it holds that \(\beta _6 \le 0\).

    For the coefficients \(\beta _i^2, i = 1,\cdots 5\), it requires that \(\beta _i^2 > 0\). Thus, from (A.25), it demands that

    $$\begin{aligned} \frac{\epsilon \Delta x}{\Delta t} - \alpha > 0. \end{aligned}$$
    (A.34)

    Substituting (A.31) into (A.34), we can obtain that

    $$\begin{aligned} C < \frac{1}{\alpha }. \end{aligned}$$
    (A.35)

    Thus, the constrain on C is changed into

    $$\begin{aligned} 0< C < \min \left( \frac{1}{\alpha },\frac{10 \Delta x}{3 \epsilon + 10 \Delta x \alpha } \right) . \end{aligned}$$
    (A.36)
  2. 2.

    \(\epsilon < \Delta x\), in which case

    $$\begin{aligned} \Delta t= C \Delta x^2. \end{aligned}$$
    (A.37)

    Substituting (A.37) into (A.27), we can deduce that

    $$\begin{aligned} \beta _6 = \frac{9}{10}C^2\epsilon \Delta x^2\left( \frac{1}{\epsilon - \alpha C \Delta x}\right) + 6\alpha \Delta x C(\alpha C \Delta x - \epsilon ) - 3C \Delta x^2. \end{aligned}$$
    (A.38)

    Thus if

    $$\begin{aligned} 0< C < \min \left( \frac{\epsilon }{\alpha \Delta x},\frac{10 \epsilon }{3 \epsilon + 10 \Delta x \alpha } \right) , \end{aligned}$$
    (A.39)

    it holds that \(\beta _6 \le 0\). Similarly, we can verify that the constrain \(\beta _i^2 > 0, i = 1, \cdots 5\) will not affect the condition (A.39), then the proof is finished.

    For \(\epsilon < \Delta x\), it is always true that \(\alpha = \exp (-1/\epsilon ^2)\) is quite small, and (A.39) could be reduced into

    $$\begin{aligned} 0< C < \frac{10}{3}. \end{aligned}$$
    (A.40)

\(\square \)

1.5 Analysis of the Higher-Order Scheme

From the test of the AP property for the numerical scheme, we found that even for the IMEX3 scheme with WENO reconstruction, the convergence order is only two. Analysis of the numerical scheme shows that when solving \(T^{n+1}\), the fourth-order polynomial equation of \(T^{n+1}\) is solved, where \((T^{n+1})^{4}\) is approximated as

$$\begin{aligned} (T^{4})_{i} \approx (T_{i})^{4} \end{aligned}$$
(A.41)

instead of

$$\begin{aligned} (T)^{4}_{i} \approx \frac{\int _{x_{i-\frac{1}{2}}}^{x_{i + \frac{1}{2}}}T^{4}dx}{\Delta x}, \end{aligned}$$
(A.42)

where \(T_{i}\) is the cell average of cell i. Noting that

$$\begin{aligned} \frac{\int _{ x_{i - \frac{1}{2} }}^{ x_{i + \frac{1}{2} }} T(x) dx }{\Delta x} = T(x_{i}) + \frac{1}{24}(T(\xi _{i}))^{''}\Delta x^2,\quad \xi \in [x_{i - \frac{1}{2}}, x_{i + \frac{1}{2}}], \end{aligned}$$
(A.43)

and

$$\begin{aligned} \frac{ \int _{x_{i - \frac{1}{2}}}^{x_{i + \frac{1}{2}}} T^{4}(x) dx }{\Delta x} = T^{4}(x_{i}) + \frac{1}{24}(T^{4}(\eta _{i}))^{''}\Delta x^2,\eta _{i}\in [x_{i - \frac{1}{2}}, x_{i + \frac{1}{2}}], \end{aligned}$$
(A.44)

thus, it holds

$$\begin{aligned} \frac{\int _{x_{i - \frac{1}{2}}}^{x_{i + \frac{1}{2}}} T^{4} dx}{\Delta x} - (T_{i})^{4} = {\mathcal {O}}(\Delta x^{2}). \end{aligned}$$
(A.45)

Therefore, the convergence order of the whole numerical scheme is at most two.

Rights and permissions

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Fu, J., Li, W., Song, P. et al. An Asymptotic-Preserving IMEX Method for Nonlinear Radiative Transfer Equation. J Sci Comput 92, 27 (2022). https://doi.org/10.1007/s10915-022-01870-3

Download citation

  • Received:

  • Revised:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1007/s10915-022-01870-3

Keywords