## 1 Introduction

The equation of radiation hydrodynamics describes the radiative transport through a fluid with coupled momentum and energy exchange Radiation-1 ; Radiation-2 ; Radiation-3 . Its applications are mainly in high-temperature hydrodynamics, including gaseous stars in astrophysics, supernova explosions, combustion phenomena, reentry vehicles, fusion physics and inertial confinement fusion. The importance of thermal radiation increases as the temperature is raised in the above problems. Such as for the moderate temperature, the role of radiation is primarily one of transporting energy by radiative process. But for the higher temperature, the energy and momentum densities of the radiation field may become comparable to or even dominates the corresponding fluid quantities.

In the case of the zero diffusion limit radiative-GKS-0 ; radiative-GKS-1 ; radiative-GKS-2 , the equation of radiation hydrodynamics can be written into a nonlinear hyperbolic system of conservation laws. But for the more complicated equilibrium-diffusion limit Radiation-Diffusion-1 , another nonlinear diffusion term for radiative heat transfer should be added. Due to the highly non-linearity of this radiation diffusion terms, it becomes more challenge to design a high-order and robust numerical method. For the scales of characteristic time between the radiation and hydrodynamics are different by several orders of magnitude, and it usually requires the radiation part to be solved implicitly to guarantee the numerical stability. There are many numerical method for the radiation hydrodynamics in equilibrium-diffusion limit. With the operator splitting method, the Godunov schemes were proposed for the hyperbolic part and an implicit scheme is proposed for the radiative heat transfer Radiation-Diffusion-1 ; Radiation-Diffusion-2 ; Radiation-Diffusion-3 . The only second-order accuracy can be achieved in space and time for both the equilibrium diffusion and streaming limit, and it is also capable of computing radiative shock solutions accurately. Furthermore, to achieve the high-order accuracy, one-dimension implicit-explicit (IMEX) Lagrangian high-order scheme was developed in Radiation-Diffusion-4 . The essentially non-oscillatory (ENO) ENO method is used for the advection and radiation diffusion term to obtain the high spatial accuracy, and the strong stability preserving method is used for high order temporal accuracy. The more work on the RHEs’ computation can be found in Radiation-Diffusion-5 ; Radiation-Diffusion-6 ; Radiation-Diffusion-7 .

In the last decades, the gas-kinetic scheme (GKS) and based on the Bhatnagar-Gross-Krook (BGK) model BGK-1 ; BGK-2 have been developed systematically for the computations from low speed flows to supersonic ones GKS-Xu1 ; GKS-Xu2 . The gas-kinetic scheme is based on an analytical integral solution of the BGK equation, and gas distribution function at a cell interface provides a multi-scale evolution process from the kinetic particle transport to the hydrodynamic wave propagation. With the two-stage fourth-order method for Lax-Wendroff type flow solvers GRP-high-1 ; GRP-high-2 , the high-order gas-kinetic schemes were developed GKS-high-1 ; GKS-high-2 . The high-order scheme not only reduces the complexity of computation, but also improves the accuracy of the numerical solution. Most importantly, the robustness is as good as the second-order shock capturing scheme. Furthermore, with the discretization of particle velocity space, a unified gas-kinetic scheme (UGKS) has been developed for the flow study in entire Knudsen number regimes from rarefied to continuum ones UGKS-Xu1 ; UGKS-Xu2 ; UGKS-Xu3 . Recently, the UGKS is extended to solve radiative transfer system with both scattering and absorption/emission effects radiative-UGKS-1 ; radiative-UGKS-2 ; radiative-UGKS-3 . The asymptotic preserving (AP) property can be accurately recovered. For the equation of radiation hydrodynamics, a multi-scale scheme is developed, in which GKS is used for the compressible inviscid flow and UGKS is used for the non-equilibrium radiative transfer radiative-UGKS-4 . Due to the possible large variation of fluid opacity in different regions, the transport of photons through the flow system is simulated by the multi-scale scheme.

In this paper, a high-order gas-kinetic scheme is proposed for the equation of radiation hydrodynamic in the equilibrium-diffusion limit. Based on the zeroth-order Chapman-Enskog expansion, the hydrodynamic part of radiation hydrodynamic equation can be obtained from the modified BGK equation with modified equilibrium state. The radiation diffusion term is considered as source term. Since the time scales of radiation diffusion and hydrodynamic part are different and it will make the time step of an explicit scheme very small. Thus, an implicit-explicit (IMEX) scheme is developed for solving the radiation hydrodynamic equation, in which the fluid advection term is treated explicitly and the radiation diffusion is treated implicitly. For the hydrodynamic part, the gas-kinetic solver with the modified equilibrium state is used to solve the compressible flow equations. The nonlinear Newton-GMRES method GMRES-1 ; GMRES-2 is used to deal with the radiation diffusion. To achieve the temporal accuracy, the two-stage third-order temporal discretization is developed, which is an extension of two-stage method was developed for Lax-Wendroff type flow solvers. To achieve the spatial accuracy, the classical weighted essentially non-oscillatory (WENO) WENO-JS ; WENO-Z method reconstruction is used. With the two-stage temporal discretization and WENO reconstruction, a reliable framework was provided for equation of radiation hydrodynamics. Various numerical experiments are carried out to validate the performance of current scheme.

This paper is organized as follows. In Section 2, the equation of radiation hydrodynamics and corresponding BGK model are introduced. The high-order gas-kinetic for RHE is presented in Section 3. Numerical examples are included in Section 4 and the last section is the conclusion.

## 2 Equation of radiation hydrodynamics and BGK model

### 2.1 Equation of radiation hydrodynamics

The equation of radiation hydrodynamics (RHE) describes the motion of flows under a radiation field. It consists of the Euler equations coupling with the radiation momentum and energy sources and the radiation-transport equation, and can be given as follows

where , , and are the density, velocity, total energy and pressure of the matter, respectively. is the radiation momentum source, is the radiation energy source, is the speed of light, is radiation intensity, is the photon direction of flight and is the angle and frequency dependent radiation source representing the radiation-matter interaction. The radiation-transport equation is essentially the conservation of the photon number, which reveals the relationship between the photon free transport and radiation-matter interaction, i.e. photon emission, photon absorption and photon scatter. The source terms and can be written as the zeroth and first order frequency-integrated angular moments of , respectively,

where is radiation energy density, is radiation flux and is radiation-pressure. The equilibrium-diffusion approximation imposes four basic assumptions to simplify RHE, i.e. the photon mean-free-path is small compared to the size of the absorption-dominated system, the matter-radiation system is in thermal equilibrium, the radiation flux is diffusive, and the radiation pressure is isotropic. It is also assumed that the radiative temperature and the fluid temperature are equal and the gas is radiatively opaque so that the equilibrium diffusion will be dealt with. With the assumptions above, the radiation energy density, radiation flux and radiation-pressure can be simplified as

where is temperature, is diffusion constant and is radiation constant representing the ratio of the radiation energy to the material energy. The time-derivative of radiation flux in the total-momentum can be dropped in accordance with the diffusion approximation. Therefore, the equation of radiation hydrodynamics in the equilibrium-diffusion limit can be written as

(1) |

where the total energy and the total pressure are given as

The polytropic ideal gas is considered, and the equation of state is given by

where is the specific internal energy, is the heat capacity at constant volume and is the specific heat ratio.

### 2.2 BGK model

In this paper, a high-order gas-kinetic scheme will be presented for two-dimensional flows. To recover the macroscopic equation of radiation hydrodynamics Eq.(1), the modified BGK equation BGK-1 ; BGK-2 can be written as

(2) |

where is the gas distribution function, is the equilibrium distribution, is the particle velocity and is the collision time. For the equation of radiation hydrodynamics, a modified equilibrium state function is introduced radiative-GKS-2

where for two dimensional system, is the macroscopic velocity, the internal variables are defined as and

and the internal degrees of freedom of

and satisfyThe parameters can be given by

The collision term also satisfies the compatibility condition

where and .

According to the Chapman-Enskog expansion, the Euler and Navier-Stokes equations can be derived form the BGK equation GKS-Xu1 ; GKS-Xu2 . Similarly, the macroscopic equations Eq.(1) can be also derived from Eq.(2). Taking zeroth-order Chapman-Enskog expansion, i.e. with , and taking moments of BGK equation Eq.(2), we have

where

and

According to the definition of and the relation of , the hyperbolic part of macroscopic equations Eq.(1

) can be recovered and the vector form is used in the following sections

where is the conservative variable, and are the fluxes in and directions and is the source term for radiation diffusion.

## 3 High-order gas-kinetic scheme

### 3.1 Temporal discretization

Recently, based on the time-dependent flux function of the generalized Riemann problem solver (GRP) GRP-high-1 ; GRP-high-2 and gas-kinetic scheme (GKS) GKS-high-1 ; GKS-high-2 , a two-stage fourth-order time-accurate discretization was developed for Lax-Wendroff type flow solvers, particularly applied for the hyperbolic conservation laws. Considering the following time-dependent equation with the initial condition

where is an operator for spatial derivative of flux. Introducing an intermediate state at , the two-stage temporal discretization can be written as

(3) |

It can be proved that for the hyperbolic equations the two-stage time stepping method Eq.(3) provides a fourth-order time accurate solution for at . Based on the high-order spatial reconstruction WENO-JS ; WENO-Z , successes have also been achieved for the construction of high-order gas-kinetic scheme for Euler and Navier-Stokes equations GKS-high-1 ; GKS-high-2 . The two-stage method provides a reliable framework to develop high-order scheme with the implementation of second-order flux function. Most importantly, due to the use of both flux function and its temporal derivative, this scheme is robust and works perfectly from the subsonic to the hypersonic flows.

In this paper, the high-order gas-kinetic scheme will be developed for the radiation hydrodynamics equation Eq.(1) as well. For simplicity, the two-dimensional uniform mesh is used. Taking moments of the BGK equation Eq.(2) and integrating with respect to space for the cell , the semi-discrete finite volume scheme can be written as

(4) |

where is the cell averaged conservative variables over the cell . The operator for hydrodynamic part is given by

(5) |

and the operator is given by

where and are the cell size, and are the time dependent numerical fluxes at cell interfaces in and directions and is the source for radiation diffusion. Without considering the source term, the explicit scheme is used for the Euler and Navier-Stokes equations GKS-high-1 . However, the time scales of radiation diffusion and fluid advection are different. It will make the time step of an explicit scheme very small, and the explicit scheme can be only used efficiently for the mildly non-relativistic regime for RHE. To improve the efficiency, the RHE should be discretized in an implicit-explicit (IMEX) procedure, i.e. the fluid advection term is treated explicitly and the radiation component is treated implicitly. Similar with Eq.(3), introducing an intermediate state at , the two-stage temporal discretization for Eq.(4) is given as follows

(6) |

Compared with the original two-stage method Eq.(3), the trapezoid integration is used for the source terms.

It can be proved that Eq.(6) provides a third-order accurate approximation for the system with sources Eq.(4) at . Integrating Eq.(4) on the time interval , we have

To prove this proposition above, the following Taylor expansion need to be satisfied

(7) | ||||

According to Eq.(4) and Cauchy-Kovalevskaya method, the temporal derivatives can be given by

and

For the operator , we have the following expansion up to the corresponding order

(8) |

Substituting Eq.(8) into Eq.(6), we have

(9) |

and

(10) |

For the operator , we have the following expansion up to the corresponding order as well

(11) |

Substituting Eq.(3.1) and Eq.(11) into Eq.(3.1), it is easy to verify

Therefore, the two-stage method Eq.(6) provides a third-order temporal discretization for radiation hydrodynamic equations.

### 3.2 Discretization for hydrodynamic part

In the following subsections, the implementation of the hydrodynamic part and radiative part will be given for Eq.(6). For the hydrodynamic part, the numerical flux in -direction can be given by Gaussian quadrature

(12) |

where is the Gaussian quadrature point and are quadrature weights. To construct the gas distribution function at the cell interface, the integral solution of BGK equation Eq.(2) is used

(13) |

where is the initial gas distribution function, is the corresponding equilibrium state, and and are the trajectory of particles. Similar with the gas-kinetic scheme for Euler and Navier-Stokes equations, the second-order gas-kinetic solver GKS-Xu2 can be written as follows

(14) |

The coefficients in Eq.(3.2) can be determined by the reconstructed derivatives and compatibility condition

and

where the moments of the equilibrium are defined by

Compared with the Euler and Navier-Stokes equations, the procedure for RHE is a little more complicated. As an example, the spatial derivative of can be written as

(15) |

where

Taking moments of Eq.(15), we have

To obtain the connections , the equation can be given as a linear system

(16) |

where

The detailed formulation of can be found in Appendix. The system Eq.(16) seems to be under determined, but we can also find the unique solution.

According to the chain rule, the derivative of

can be also written asIt’s easily to verify that

where and are the functions without and , respectively. Comparing the coefficients of and with Eq.(15), we have

According to the definition of and , we have

where

With the solution for and , the unique solution of the linear system is given by

where

Similarly, the coefficients for temporal derivative can be also determined. Thus, the gas distribution function Eq.(3.2) and the numerical fluxes Eq.(12) can be fully constructed.

In order to utilize the two-stage temporal discretization, the temporal derivatives of the flux function need to be determined. The flux in the time interval is expanded as the following linear form

(17) |

where the coefficients and can be fully determined by solving the linear system

Similarly, the coefficients and corresponding to the flux in -direction can be constructed as well. According to Eq.(5), and its temporal derivative at can be given by

With the procedure at the intermediate state, can be constructed as well.

### 3.3 Discretization for radiative part

According to the definition of source terms, we only need to solve energy equations at two-stages for Eq.(6), which are nonlinear systems with respect to the temperature . At each stage, the energy equations can be simplified as the nonlinear system

(18) |

where is a nonlinear function from to and is the number of cells. To obtain an approximate solution of the nonlinear system, the above system is rewritten as

where is the Jacobian of . In order to avoid calculating the Jabobian directly, can be approximated by the -derivative as follows

(19) |

where is a small scalar. In fact, it is not trivial to choose , which has great influence on the stability of algorithm and improper values will cause blowing up. According to GMRES-2 , a special is given as follows

where is the machine epsilon, , , is the typical size of and the typical size of real number equals to its order of magnitude plus one. With the initial guess for Eq.(18) and for Eq.(19), the initial residual is

and its generated Krylov subspace is

According to the nonlinear Newton-GMRES method GMRES-2 , an approximate solution can be given. Besides, in order to improve computation efficiency, a restart algorithm is applied in orthogonalization process GMRES-1 . In actual calculation, the restart times can be no more than 5 times if the restart step and convergence condition are set appropriately.

## 4 Numerical tests

In this section, the one and two dimensional tests cases are provided to the correctness and robustness of current scheme. In the computation, the collision time takes

where and denote the pressure on the left and right sides of the cell interface. In smooth flow regions, the collision term takes and the gas distribution function Eq.(3.2) will reduce to

To achieve the spatial accuracy, the classical multidimensional fifth-order WENO reconstruction is adopted and more details can be found in GKS-high-1 ; WENO-Z . In order to eliminate the spurious oscillation and improve the stability, the WENO reconstruction is performed for the characteristic variables, and the detailed analysis for characteristics can be found in Radiation-Diffusion-4 . Without special statement, the gas with is used in the following numerical examples.

### 4.1 Accuracy test

The advection of density perturbation problem for Euler equations is extended to test the order of accuracy for RHE. In this test case, the initial condition is set as follows

where the computational domain is , the periodic boundary condition is applied and .

As reference, the case with and is tested. With these parameters, the RHE degenerates into the Euler equation and the exact solution can be given by

Comments

There are no comments yet.