Structure-preserving operators for thermal-nonequilibrium hydrodynamics

TL;DR

This paper introduces a structure-preserving numerical scheme for radiation hydrodynamics that maintains energy conservation and accurately captures shock discontinuities by discretizing conservation laws with algebraic transformations.

Contribution

The authors develop a novel discretization strategy that preserves the mathematical structure of conservation laws at the discrete level, ensuring physical fidelity in simulations.

Findings

Global conservation errors are within round-off level.

The scheme accurately reproduces shock jump conditions.

Flexible to various numerical methods like central differences and Runge-Kutta.

Abstract

Radiation hydrodynamics simulations based on the one-fluid two-temperature model may violate the law of energy conservation because the governing equations are expressed in a nonconservative formulation. Here, we maintain the important physical requirements by employing a strategy based on the key concept that the mathematical structures associated with the conservative and nonconservative equations are preserved, even at the discrete level. To this end, we discretize the conservation laws and transform them via exact algebraic operations. The proposed scheme maintains the global conservation errors within the round-off level. In addition, a numerical experiment concerning the shock tube problem suggests that the proposed scheme well agrees with the jump conditions at the discontinuities regulated by the Rankine-Hugoniot relationship. The generalized derivation allows us to employ arbitrary central difference, artificial dissipation, and Runge-Kutta methods.

Equations

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{u})=0,$$

$$\frac{\partial \rho \mathbf{u}}{\partial t}+\nabla \cdot (\rho \mathbf{uu}+p_{\mathrm{i}}+p_{\mathrm{e}})=0,$$

$$\frac{\partial \rho e_{\mathrm{i}}}{\partial t}+\nabla \cdot (\rho e_{\mathrm{i}}\mathbf{u})+p_{\mathrm{i}}\nabla \cdot \mathbf{u}=0,$$

$$\frac{\partial \rho e_{\mathrm{e}}}{\partial t}+\nabla \cdot (\rho e_{\mathrm{e}}\mathbf{u})+p_{\mathrm{e}}\nabla \cdot \mathbf{u}=0,$$

$$\frac{\partial \left(\rho e_{\mathrm{i}}+\rho e_{\mathrm{e}}+\frac{1}{2}\rho |\mathbf{u}{|}^2\right)}{\partial t}+\nabla \cdot \left(\rho e_{\mathrm{i}}\mathbf{u}+\rho e_{\mathrm{e}}\mathbf{u}+\frac{1}{2}\rho |\mathbf{u}{|}^2\mathbf{u}+p_{\mathrm{i}}\mathbf{u}+p_{\mathrm{e}}\mathbf{u}\right)=0.$$

$$\frac{\partial \rho }{\partial t}+\rho \nabla \cdot \mathbf{u}+\mathbf{u}\cdot \nabla \rho =0.$$

$$\frac{{\rho }_j^{n+1}-{\rho }_j^n}{\Delta t}+\frac{⟨\rho u{⟩}_{j^{+}}^n-⟨\rho u{⟩}_{j^{-}}^n}{\Delta x}=0,$$

$$\frac{{\rho }_j^{n+1}-{\rho }_j^n}{\Delta t}+{\rho }_j^n\frac{⟨u{⟩}_{j^{+}}^n-⟨u{⟩}_{j^{-}}^n}{\Delta x}+u_j^n\frac{⟨\rho {⟩}_{j^{+}}^n-⟨\rho {⟩}_{j^{-}}^n}{\Delta x}=0,$$

$$\sum _{j=1}^N\frac{{\rho }_j^{n+1}-{\rho }_j^n}{\Delta t}+\sum _{j=1}^N\frac{⟨\rho u{⟩}_{j^{+}}^n-⟨\rho u{⟩}_{j^{-}}^n}{\Delta x}=0,$$

$$\sum _{j=1}^N\frac{{\rho }_j^{n+1}-{\rho }_j^n}{\Delta t}+\sum _{j=1}^N{\rho }_j^n\frac{⟨u{⟩}_{j^{+}}^n-⟨u{⟩}_{j^{-}}^n}{\Delta x}+\sum _{j=1}^N{u}_j^n\frac{⟨\rho {⟩}_{j^{+}}^n-⟨\rho {⟩}_{j^{-}}^n}{\Delta x}=0,$$

$$\sum _{j=1}^N\frac{{\rho }_j^{n+1}-{\rho }_j^n}{\Delta t}+\frac{⟨\rho u{⟩}_{N+1/2}^n-⟨\rho u{⟩}_{1/2}^n}{\Delta x}=0.$$

$$\frac{f^{n+1}{g}^{n+1}-f^n{g}^n}{\Delta t}$$

$$=f^{n+1}\frac{g^{n+1}-g^n}{\Delta t}+\frac{f^{n+1}-f^n}{\Delta t}{g}^n,$$

$$\frac{f^{n+1}{g}^{n+1}-f^n{g}^n}{\Delta t}=\frac{f^{n+1}+f^n}{2}\frac{g^{n+1}-g^n}{\Delta t}+\frac{f^{n+1}-f^n}{\Delta t}\frac{g^{n+1}+g^n}{2}.$$

$$\frac{f^{n+1}/g^{n+1}-f^n/g^n}{\Delta t}$$

$$\frac{{\rho }_j^{n+1}-{\rho }_j^n}{\Delta t}+\frac{⟨\rho u{⟩}_{j^{+}}^n-⟨\rho u{⟩}_{j^{-}}^n}{\Delta x}=0,$$

$$\frac{(\rho u{)}_j^{n+1}-(\rho u{)}_j^n}{\Delta t}+\frac{⟨\rho u^2+p{⟩}_{j^{+}}^n-⟨\rho u^2+p{⟩}_{j^{-}}^n}{\Delta x}=0,$$

$$\frac{(\rho e+\frac{1}{2}\rho u^2{)}_j^{n+1}-(\rho e+\frac{1}{2}\rho u^2{)}_j^n}{\Delta t}+\frac{⟨\rho eu+\frac{1}{2}\rho u^3+pu{⟩}_{j^{+}}^n-⟨\rho eu+\frac{1}{2}\rho u^3+pu{⟩}_{j^{-}}^n}{\Delta x}=0,$$

$$\frac{(\rho u^2{)}_j^{n+1}-(\rho u^2{)}_j^n}{\Delta t}$$

$$=\frac{{\rho }_j^{n+1}+{\rho }_j^n}{2{\rho }_j^{n+1}{\rho }_j^n}\frac{\{(\rho u{)}^2{\}}_j^{n+1}-\{(\rho u{)}^2{\}}_j^n}{\Delta t}-\frac{\{(\rho u{)}^2{\}}_j^{n+1}+\{(\rho u{)}^2{\}}_j^n}{2{\rho }_j^{n+1}{\rho }_j^n}\frac{{\rho }_j^{n+1}-{\rho }_j^n}{\Delta t},$$

$$=\frac{{\rho }_j^{n+1}+{\rho }_j^n}{2{\rho }_j^{n+1}{\rho }_j^n}\left\{(\rho u{)}_j^{n+1}+(\rho u{)}_j^n\right\}\frac{(\rho u{)}_j^{n+1}-(\rho u{)}_j^n}{\Delta t}-\frac{\{(\rho u{)}^2{\}}_j^{n+1}+\{(\rho u{)}^2{\}}_j^n}{2{\rho }_j^{n+1}{\rho }_j^n}\frac{{\rho }_j^{n+1}-{\rho }_j^n}{\Delta t},$$

$$=\frac{{\rho }_j^{n+1}+{\rho }_j^n}{2{\rho }_j^{n+1}{\rho }_j^n}\left\{(\rho u{)}_j^{n+1}+(\rho u{)}_j^n\right\}\left(-\frac{⟨\rho u^2+p{⟩}_{j^{+}}^n-⟨\rho u^2+p{⟩}_{j^{-}}^n}{\Delta x}\right)$$

$$+\frac{\{(\rho u{)}^2{\}}_j^{n+1}+\{(\rho u{)}^2{\}}_j^n}{2{\rho }_j^{n+1}{\rho }_j^n}\frac{⟨\rho u{⟩}_{j^{+}}^n-⟨\rho u{⟩}_{j^{-}}^n}{\Delta x}.$$

$$\frac{(\rho e{)}_j^{n+1}-(\rho e{)}_j^n}{\Delta t}+\frac{⟨\rho eu{⟩}_{j^{+}}^n-⟨\rho eu{⟩}_{j^{-}}^n}{\Delta x}+\frac{⟨pu{⟩}_{j^{+}}^n-⟨pu{⟩}_{j^{-}}^n}{\Delta x}-\frac{{\rho }_j^{n+1}+{\rho }_j^n}{4{\rho }_j^{n+1}{\rho }_j^n}\{(\rho u{)}_j^{n+1}+(\rho u{)}_j^n\}\frac{⟨p{⟩}_{j^{+}}^n-⟨p{⟩}_{j^{-}}^n}{\Delta x}=$$

$$-\frac{({\rho }^2{u}^2{)}_j^{n+1}+({\rho }^2{u}^2{)}_j^n}{4{\rho }_j^{n+1}{\rho }_j^n}\frac{⟨\rho u{⟩}_{j^{+}}^n-⟨\rho u{⟩}_{j^{-}}^n}{\Delta x}+\frac{{\rho }_j^{n+1}+{\rho }_j^n}{4{\rho }_j^{n+1}{\rho }_j^n}\{(\rho u{)}_j^{n+1}+(\rho u{)}_j^n\}\frac{⟨\rho u^2{⟩}_{j^{+}}^n-⟨\rho u^2{⟩}_{j^{-}}^n}{\Delta x}$$

$$-\frac{1}{2}\frac{⟨\rho u^3{⟩}_{j^{+}}^n-⟨\rho u^3{⟩}_{j^{-}}^n}{\Delta x}.$$

$$\frac{\partial (\rho e)}{\partial t}+\frac{\partial (\rho eu)}{\partial x}+p\frac{\partial u}{\partial x}=0\mathrm{or}\frac{\partial (\rho e)}{\partial t}+\frac{\partial (\rho eu)}{\partial x}+\frac{\partial (pu)}{\partial x}-u\frac{\partial p}{\partial x}=0.$$

$$\frac{(\rho e_{\mathrm{s}}{)}_j^{n+1}-(\rho e_{\mathrm{s}}{)}_j^n}{\Delta t}+\frac{⟨\rho e_{\mathrm{s}}u{⟩}_{j^{+}}^n-⟨\rho e_{\mathrm{s}}u{⟩}_{j^{-}}^n}{\Delta x}+\frac{⟨p_{\mathrm{s}}u{⟩}_{j^{+}}^n-⟨p_{\mathrm{s}}u{⟩}_{j^{-}}^n}{\Delta x}-\frac{{\rho }_j^{n+1}+{\rho }_j^n}{4{\rho }_j^{n+1}{\rho }_j^n}\{(\rho u{)}_j^{n+1}+(\rho u{)}_j^n\}\frac{⟨p_{\mathrm{s}}{⟩}_{j^{+}}^n-⟨p_{\mathrm{s}}{⟩}_{j^{-}}^n}{\Delta x}$$

$$=-\frac{({\rho }^2{u}^2{)}_j^{n+1}+({\rho }^2{u}^2{)}_j^n}{8{\rho }_j^{n+1}{\rho }_j^n}\frac{⟨\rho u{⟩}_{j^{+}}^n-⟨\rho u{⟩}_{j^{-}}^n}{\Delta x}+\frac{{\rho }_j^{n+1}+{\rho }_j^n}{8{\rho }_j^{n+1}{\rho }_j^n}\{(\rho u{)}_j^{n+1}+(\rho u{)}_j^n\}\frac{⟨\rho u^2{⟩}_{j^{+}}^n-⟨\rho u^2{⟩}_{j^{-}}^n}{\Delta x}$$

$$-\frac{1}{4}\frac{⟨\rho u^3{⟩}_{j^{+}}^n-⟨\rho u^3{⟩}_{j^{-}}^n}{\Delta x}.$$

$$\frac{{\rho }_j^{n+1}-{\rho }_j^n}{\Delta t}+\frac{⟨\rho u{⟩}_{j^{+}}^n-⟨\rho u{⟩}_{j^{-}}^n}{\Delta x}=0,$$

$$\frac{(\rho u{)}_j^{n+1}-(\rho u{)}_j^n}{\Delta t}+\frac{⟨\rho u^2+p{⟩}_{j^{+}}^n-⟨\rho u^2+p{⟩}_{j^{-}}^n}{\Delta x}=\frac{⟨A{⟩}_{j^{+}}^n-⟨A{⟩}_{j^{-}}^n}{\Delta x},$$

$$\frac{(\rho e+\frac{1}{2}\rho u^2{)}_j^{n+1}-(\rho e+\frac{1}{2}\rho u^2{)}_j^n}{\Delta t}+\frac{⟨\rho eu+\frac{1}{2}\rho u^3+pu{⟩}_{j^{+}}^n-⟨\rho eu+\frac{1}{2}\rho u^3+pu{⟩}_{j^{-}}^n}{\Delta x}=\frac{⟨B{⟩}_{j^{+}}^n-⟨B{⟩}_{j^{-}}^n}{\Delta x},$$

$$\frac{(\rho e_{\mathrm{s}}{)}_j^{n+1}-(\rho e_{\mathrm{s}}{)}_j^n}{\Delta t}+\frac{⟨\rho e_{\mathrm{s}}u{⟩}_{j^{+}}^n-⟨\rho e_{\mathrm{s}}u{⟩}_{j^{-}}^n}{\Delta x}+\frac{⟨p_{\mathrm{s}}u{⟩}_{j^{+}}^n-⟨p_{\mathrm{s}}u{⟩}_{j^{-}}^n}{\Delta x}-\frac{{\rho }_j^{n+1}+{\rho }_j^n}{4{\rho }_j^{n+1}{\rho }_j^n}\{(\rho u{)}_j^{n+1}+(\rho u{)}_j^n\}\frac{⟨p_{\mathrm{s}}{⟩}_{j^{+}}^n-⟨p_{\mathrm{s}}{⟩}_{j^{-}}^n}{\Delta x}$$

$$=-\frac{{\rho }_j^{n+1}+{\rho }_j^n}{8{\rho }_j^{n+1}{\rho }_j^n}\{(\rho u{)}_j^{n+1}+(\rho u{)}_j^n\}\frac{⟨A{⟩}_{j^{+}}^n-⟨A{⟩}_{j^{-}}^n}{\Delta x}+\frac{⟨B_{\mathrm{s}}{⟩}_{j^{+}}^n-⟨B_{\mathrm{s}}{⟩}_{j^{-}}^n}{\Delta x}-\frac{({\rho }^2{u}^2{)}_j^{n+1}+({\rho }^2{u}^2{)}_j^n}{8{\rho }_j^{n+1}{\rho }_j^n}\frac{⟨\rho u{⟩}_{j^{+}}^n-⟨\rho u{⟩}_{j^{-}}^n}{\Delta x}$$

$$+\frac{{\rho }_j^{n+1}+{\rho }_j^n}{8{\rho }_j^{n+1}{\rho }_j^n}\{(\rho u{)}_j^{n+1}+(\rho u{)}_j^n\}\frac{⟨\rho u^2{⟩}_{j^{+}}^n-⟨\rho u^2{⟩}_{j^{-}}^n}{\Delta x}-\frac{1}{4}\frac{⟨\rho u^3{⟩}_{j^{+}}^n-⟨\rho u^3{⟩}_{j^{-}}^n}{\Delta x}.$$

$$\frac{\partial \rho }{\partial t}+\frac{\partial (\rho u)}{\partial x}=0,$$

$$\frac{\partial (\rho u)}{\partial t}+\frac{\partial (\rho u^2+p_{\mathrm{i}}+p_{\mathrm{e}})}{\partial x}=0,$$

$$\frac{\partial (\rho e_{\mathrm{i}}+\rho e_{\mathrm{e}}+\frac{1}{2}\rho u^2)}{\partial t}+\frac{\partial (\rho e_{\mathrm{i}}u+\rho e_{\mathrm{e}}u+\frac{1}{2}\rho u^3+p_{\mathrm{i}}u+p_{\mathrm{e}}u)}{\partial x}=0.$$

$$\frac{1}{4}⟨f{⟩}_{j-1/2}+⟨f{⟩}_{j+1/2}+\frac{1}{4}⟨f{⟩}_{j+3/2}=\frac{3}{2}(f_j+f_{j+1}),$$

$$⟨f_1+f_2⟩=⟨f_1⟩+⟨f_2⟩,$$

Full text

††journal: Journal of Computational Physics

1 Introduction

Radiation hydrodynamics (RHD)[1] is one of the major techniques employed in laser-plasma simulations. In RHD, a neutral charge is assumed for a fluid composed of ions and electrons. This assumption allows large-scale plasma simulations to be performed, as the grid interval is not limited by the Debye length, unlike the conventional particle-in-cell (PIC) method. This type of simulation has been employed in laser-plasma simulations that model, for example, the implosion dynamics of inertial confinement fusion (ICF)[2].

One of the most simplified RHD simulations combines radiative transfer and one-temperature hydrodynamics, which are employed by the FastRad3D code[3]. However, thermal nonequilibrium typically exists between the ions and electrons of a laser plasma. The incident laser is absorbed by the inverse bremsstrahlung process, and the energy is deposited on the electrons. Thus, the one-fluid two-temperature (1F2T) model has been employed in order to investigate plasma hydrodynamics more accurately [4, 5, 6]. The governing equations of the 1F2T model excluding viscous and heat conduction effects are expressed in the form

[TABLE]

where $\rho$ is the mass density, $\mathbf{u}=^{\mathrm{T}}[u,v,w]$ is the flow velocity, $e$ is the specific internal energy, and $p$ is the pressure. The subscripts “$\mathrm{i}$” and “$\mathrm{e}$” denote the ions and electrons, respectively. Note that nonhydrodynamic energy transports such as radiative transfer are neglected. These nonconservative equations have been used in many RHD simulations [7, 8].

Three main methods of solving the 1F2T model exist. The first (hereafter referred to as “nonconservative method”) discretizes the governing equations directly with the Lagrangian or arbitrary Lagrangian–Eulerian (ALE) grids[9]. This method has been adopted for use in the RHD codes known as HYDRA[10], DRACO[11, 12], and PINOCO[13], which are specialized for ICF implosions. However, because the nonconservative equations are discretized directly using the artificial dissipation, the law of energy conservation is not satisfied, especially near the discontinuities such as a shock wave. For example, it has been reported that Helios-CR[14], a one-dimensional RHD code, has an error of several percent with regard to the law of energy conservation in the typical cases. The second discretization technique is based on approximate Riemann solvers and methods of this kind are employed by CRASH[15], FLASH[16], RAGE[17], and RAICHO[18], which are primarily implemented in astrophysical simulations. Note that some of these codes incorporate adaptive mesh refinement (AMR) for the performance of high-resolution simulations, and cross-code comparisons are also conducted[19]. These codes are designed to solve the Euler equation and the hydrodynamic heating is artificially divided between the ions and electrons in proportion to the pressure ratio. Therefore, the results are not guaranteed to converge to the exact solutions, because the governing equations are not discretized directly. However, the energy conservation is satisfied even in the discrete form.

The most recently developed method (hereafter, “conservative method”) is implemented in the ASTER[20] code developed at the Laboratory of Laser Energetics, University of Rochester. This code solves the conservation laws of mass, momentum, and total energy, along with the energy equation given in Eq. (4). Thus, the conservation law of total energy is satisfied automatically and the numerical solutions converge to the exact solutions, provided the solutions do not possess any discontinuities. However, relatively high-intensity lasers ($\sim$10${}^{15}\mathrm{\ W/cm^{2}}$) are usually employed in typical laser-plasma investigations[21, 22, 23, 24], which yield strong shock waves. The internal energies of the ions and electrons may not be accurately determined near these shock waves. This problem has a critical impact on RHD simulations, because the emission and absorption coefficients of radiative transfer are strongly dependent on the electron internal energy.

In this paper, we propose a structure-preserving scheme that can overcome the above problems. The proposed scheme should achieve the following requirements: (i) solve the partially nonconservative equations (1)–(4), (ii) maintain the error of energy conservation at the round-off level, and (iii) get a good agreement with exact solutions including discontinuities. In Section 2, we derive the conservative scheme through a discussion using the discrete mathematics. The mathematical keys for the exact derivation are the product and quotient rules in the discrete calculus; therefore, these tools are also introduced. The numerical methods employed in this study are described in Section 3. Furthermore, in Section 4, the results of a numerical experiment on an extended shock-tube problem are discussed, and the superiority of the proposed scheme over the existing nonconservative and conservative methods is demonstrated. The accuracy of the proposed scheme is verified in Section 5 via a linear advection problem concerning entropy waves. Section 6 presents the conclusion of this article.

2 Mathematical background of the proposed approach

2.1 Causes of energy conservation law violation

Equations (1)–(4) can be mathematically transformed into the law of energy conservation by using the product and quotient rules,

[TABLE]

Physically, the law of energy conservation should be maintained by solving Eqs. (1)–(4). However, when these equations are solved in discrete form, the energy conservation is often violated. Sometimes, this error can be fatal to the simulations. The conservation law is violated because the energy equations (3) and (4) are expressed in nonconservative formulation. To explain why the conservation laws are violated in the discrete form, we convert Eq. (1) into nonconservative form using the product rule, as an example. Thus, Eq. (1) becomes

[TABLE]

The forward-time central-space method is then applied to Eqs. (1) and (6) to obtain the discretized equations in the finite difference method (FDM). Thus,

[TABLE]

where “$\langle\rangle$” denotes an arbitrary interpolation operator at the half points ($j^{+}=j+1/2,j^{-}=j-1/2$), and $n$ and $j$ are the indices of time and space, respectively. Note that the spatial derivatives in the $y$ and $z$ directions are omitted for simplicity. Of course, Eqs. (7) and (8) correspond to Eqs. (1) and (6), respectively. We take the summation of Eqs. (7) and (8) over the computational domain in order to discuss the global conservation in the discrete system; thus, we obtain

[TABLE]

where $N$ is the number of grids. The advection terms of Eq. (9) are cancelled, yielding

[TABLE]

Therefore, the total mass of the discrete system does not change if the periodic [$\langle\rho u\rangle_{1/2}=\langle\rho u\rangle_{N+1/2}$] or Neumann [$\langle\rho u\rangle_{1/2}=\langle\rho u\rangle_{N+1/2}=0$] boundaries are applied. In contrast, Eq. (10) violates the conservation laws, as the advection terms are not generally cancelled.

The above example is the most simple explanation of why the conservation laws are violated in the discrete form, and of why the mathematical structure of the product rule should be preserved, even at the discrete level. Below, we discuss our idea of how to construct a structure-preserving scheme.

2.2 Product and quotient rules in discrete form

The nonconservative energy equation is obtained from the conservation laws of mass, momentum, and energy. Although the derivation is performed using the product and quotient rules, these formulae may violate the global conservation in the discrete form. Therefore, we introduce the discrete product and quotient rules.

First, we introduce the product rule in discrete form:

[TABLE]

where $f$ and $g$ are arbitrary functions depending on $t$. Algebraic operations performed on the numerator, which add and subtract $f^{n+1}g^{n}$, correspond to the same approach as that used to prove the original formula in the differential form. The other forms of the product rule are obtained in a similar manner, where

[TABLE]

The quotient rule in discrete form is also obtained via the same strategy, yielding

[TABLE]

The quotient rule of Eq. (14) is mathematically equivalent to the product rule of Eq. (13); however, we must expand the kinetic energy term as a function of the mass and momentum. The quotient rule is useful when a conservative variable appears in the denominator of a function.

2.3 Derivation of structure-preserving scheme

In this subsection, we derive a scheme which exactly satisfies the law of energy conservation with solving the internal energy of ions and electrons. It is the structure-preserving scheme that the mathematical structure associating the conservative formulation with the nonconservative one is strictly maintained. The Euler equation is discretized using the FDM approach, such that

[TABLE]

where $p=p(\rho,e)$ is the total pressure. As shown in Section 2.1, the conservation laws of mass, momentum, and energy are strictly satisfied in the system. Furthermore, although we employ the Euler explicit method here, the scheme can be extended to the Runge–Kutta (RK) methods. Proof of this is given in A. The multidimensional description is given in B.

The structure-preserving scheme is obtained by expanding Eq. (17). The time derivative of the kinetic energy is expanded in the form

[TABLE]

Here, the quotient and product rules of Eqs. (14) and (13), respectively, are used in the derivation. Some readers may feel that this expansion is unnecessary, as both $\rho^{n+1}_{j}$ and $u^{n+1}_{j}$ have already been obtained in Eqs. (15) and (16). However, this operation is necessary as it clarifies the contributions of the ion and electron pressures on the right hand side (RHS) of Eq. (18). This clarification is required in order to prove that the proposed scheme qualifies as the discretized equations of Eqs. (3) and (4). In order to obtain the nonconservative formulation, Eq. (18) is substituted into Eq. (17), yielding

[TABLE]

This equation provides important information regarding the energy conservation in the 1F2T model. The left hand side (LHS) of Eq. (19) is a discretized formulation of the nonconservative energy equation

[TABLE]

The rest of the terms on the RHS are error terms that are cancelled out when $\Delta t\to 0$ and $\Delta x\to 0$. The information on this error is lost in Eq. (20); therefore, it is difficult to reconstruct the error terms using intuitive discretizations. Furthermore, note that Eq. (20) is the sum of Eqs. (3) and (4), which correspond to the first law of thermodynamics for ions and electrons. This relationship must be satisfied even at the discrete level; therefore, Eq. (19) provides a constraint condition that regulates the law of energy conservation in the discrete 1F2T model.

Finally, the structure-preserving scheme for the 1F2T model can be written in the form

[TABLE]

The terms incorporating the internal energy ($e$) or pressure ($p$) in Eq. (19) can be easily separated because of the physical requirements. However, the other terms have no explicit restrictions for partition. Equation (21) distributes the error terms equally among the ions and electrons; hence, the symmetry of Eqs. (3) and (4) is retained. The Rankine–Hugoniot relationship and the law of equipartition can only be reproduced using this approach. This will be mentioned in the later verification. Note that the extension to the multitemperature model is straightforward, because of the law of equipartition.

2.4 Shock capturing method

When solving the hydrodynamic field with the discontinuities, the shock capturing method is required to obtain the entropy solutions. Here, we derive the structure-preserving scheme including artificial dissipation terms. The discretized Euler equation Eqs. (15)–(17) is modified as

[TABLE]

where $A$ and $B$ are the artificial dissipations required to capture the discontinuities. The nonconservative equations about the internal energy of ions and electrons are derived by the same way in Sec. 2.3:

[TABLE]

Note that the first term on the RHS about $A$, which is associated with the viscosity, is equally separated between ions and electrons since the viscosity transfers the momentum whose contributions of ions and electrons are inseparable.

2.5 Summary of the proposed approach

The summary of the proposed approach toward construction of the structure-preserving scheme is as follows:

Discretize the conservation laws of mass, momentum, and total energy Eqs. (26)–(28) so that the conservation laws are automatically satisfied even in the discrete level. 2. 2.

Expand the time derivative of kinetic energy in Eq. (28) using the discrete product and quotient rules to obtain the nonconservative energy equation about $\rho(e_{\mathrm{i}}+e_{\mathrm{e}})$ with error terms. 3. 3.

Separate the contributions of ions and electrons in the discretized nonconservative energy equation. The terms with the pressure $p_{\mathrm{s}}$ and the specific internal energy $e_{\mathrm{s}}$ are separated by the physically accurate way. The rest terms only including the density $\rho$ and the velocity $u$ are mathematically equally shared by assuming the law of equipartition. 4. 4.

Solve the discretized equations about the density $\rho$, the momentum $\rho u$, the internal energy of ions $\rho e_{\mathrm{i}}$, and that of electrons $\rho e_{\mathrm{e}}$.

[TABLE]

The above methodology is also available if the artificial dissipation terms are included for the purpose of the shock capturing method.

3 Numerical implementations

In this study, the fourth-order Padé-type interpolation[25] is utilized, with

[TABLE]

where $f$ is the flux of an arbitrary conservative equation. This interpolation operator has the linearity as follows:

[TABLE]

where $\alpha$ is an arbitrary constant. Time integration is performed using the third-order total variation diminishing (TVD) RK method[26].

Artificial dissipations are modeled using the bulk viscosity $\beta$ and the thermal conductivity $\kappa_{\mathrm{s}}$ for ions and electrons, where

[TABLE]

These transport coefficients are modeled in a similar manner to the localized artificial diffusivity (LAD) scheme[27, 28]

[TABLE]

where $H$ is the Heaviside function, $a$ is the sound speed, and $r=4$, according to the typical LAD usage. $C_{\beta}$ and $C_{\kappa}$ are the nondimensional parameters of the LAD scheme and fixed to $5$ and $2$, respectively. High-order derivatives obtained from these compact schemes have noisy profiles; thus, the obtained $\beta$ and $\kappa_{\mathrm{s}}$ should be smeared using the appropriate truncated Gaussian blur[29]. The first and fourth derivatives are derived using the fourth-order compact schemes:

[TABLE]

Note that tridiagonal matrices appearing in the compact schemes are solved using the Thomas algorithm.

Hydrodynamic simulations with compact differences are often coupled with low-pass filters[30, 31], owing to the stabilization of numerical dispersion especially at high wavenumbers. Low-pass filters are usually applied to the conservative variables. The 1F2T model only has three conservation laws of mass, momentum, and energy although the number of governing equations Eqs. (1)–(4) is four. Therefore, the filtering scheme should be applied to the 1F2T model carefully. We derive a constraint condition that should be preserved through the filtering operations:

[TABLE]

where the superscript “$\star$” denotes the filtered quantities and “$\lceil\ \rfloor$” is the filtering operator. The momentum and total energy composed of the filtered primitive variables are

[TABLE]

where we assume that the filtering operator also has the linearity. Therefore, the conservation laws are maintained when the filtering operators are used like Eqs. (39)–(42). In this study, an eighth-order compact filter[32] is used in the final step of the RK integration to stabilize the numerical dispersion, such that

[TABLE]

where $b_{0}=-\frac{35}{256}$, $b_{1}=\frac{21}{256}$, $b_{2}=-\frac{7}{256}$, and $b_{3}=\frac{1}{256}$. The filtering parameter $\alpha_{\mathrm{f}}$ is set to 0.495.

4 Verification via shock tube problem

The shock tube problem is employed to verify the effects of the conservative and nonconservative schemes, as this problem includes discontinuous solutions, i.e., the shock wave and contact discontinuity[33]. Here, we extend the shock tube problem to the 1F2T model, so that the proposed scheme can be verified using the exact solution.

We assume that a diaphragm separating the high- and low-pressure sections of a tube bursts at $t=0$ and $x=0.5$. The computational domain is $0\leq x\leq 1$ and the number of grids is 201. The initial conditions are almost identical to the original Sod’s problem: $\rho_{\mathrm{L}}=1$, $\rho_{\mathrm{R}}=0.125$, $u_{\mathrm{L}}=u_{\mathrm{R}}=0$, $e_{\mathrm{i,L}}=1.5$, and $e_{\mathrm{i,R}}=e_{\mathrm{e,L}}=e_{\mathrm{e,L}}=1$, where the subscripts “$\mathrm{L}$” and “$\mathrm{R}$” denote the areas on the left and right sides of the diaphragm, respectively. The quantity $e_{\mathrm{i}}+e_{\mathrm{e}}$ is identical to the original Sod’s problem. Therefore, these conditions simply add the temperature nonequilibrium between the ions and electrons to the original problem. Note that $e$ is associated with $p$ by the thermally and calorically ideal equations of state (EOS), with the ratio of specific heat $\gamma=1.4$. The boundaries satisfy the Neumann condition such that the summations of the mass and energy are fixed to the initial values.

Here we show the exact solution of the 1F2T shock tube problem. The governing equations can be described by the following quasi-linear form:

[TABLE]

This is a hyperbolic system which has the eigenvectors $\mathbf{k}$ and the corresponding eigenvalues $\lambda$:

[TABLE]

where $a=\sqrt{\gamma(p_{\mathrm{i}}+p_{\mathrm{e}})/\rho}$ is the sound speed. The first and fourth eigenmodes denote the pressure waves while the second and third ones are the entropy waves for ions and electrons, respectively. According to the eigenstructure, the ratio of internal energy $e_{\mathrm{s}}/(e_{\mathrm{i}}+e_{\mathrm{e}})$ only changes at the contact surface. The exact solutions for each internal energies are easily derived from the original Sod’s solution with this fact.

Figure 1 shows a time history of the errors with respect to the global conservation of energy. As mentioned in Section 1, the nonconservative method obviously violates the law of energy conservation, because this method discretizes the nonconservative equations. On the other hand, the conservative and proposed methods keep the errors within the machine zero level of the double-precision floating-point numbers ($\sim$2$~{}\times 10^{-16}$). This is because the law of energy conservation is discretized directly using these techniques. The behaviors of these errors differ from each other; however, this discrepancy has no physical meaning, as these are the round-off rather than truncation errors. Therefore, the conservative and proposed schemes are proven to be conservative. Moreover, the round-off errors of these schemes do not linearly accumulate in the numerical experiments. This is primarily because we employ the finite-volume-method (FVM) like interpolation and filtering schemes[25, 32]. Note that this is the key to maintain the global conservation errors at the round-off level.

Figures 2 and 3 are spatial profiles of $e_{\mathrm{i}}$ and $e_{\mathrm{e}}$, respectively. The proposed scheme well agrees with the exact solutions including the discontinuities. However, profiles obtained by the existing nonconservative and conservative methods clearly deviate from the exact solutions between the shock wave and contact surface. The exact solutions are associated with the Rankine–Hugoniot relationship; thus, they are influenced by the local principles such as the conservation laws. Therefore, Figures 2 and 3 indicate that the nonconservative and conservative schemes violate some of these principles. For the nonconservative scheme, $e_{\mathrm{i}}$ and $e_{\mathrm{e}}$ remain identical between the discontinuities, while $e_{\mathrm{i}}+e_{\mathrm{e}}$ is lower than the exact solution. In other words, the nonconservative scheme maintains the law of equipartition but violates the law of energy conservation at the discrete level. In contrast, the conservative scheme reproduces the spatial profile of $e_{\mathrm{i}}+e_{\mathrm{e}}$ but cannot maintain the equilibrium of the ions and electrons between the discontinuities. This finding suggests that the conservative scheme maintains the law of energy conservation, but the law of equipartition is violated near the discontinuities. Note that the law of equipartition is strongly related to the symmetry of Eqs. (3) and (4). Thus, this symmetry should be preserved at the discrete level. Both the nonconservative scheme, which directly discretizes the nonconservative equations, and the proposed scheme determine $e_{\mathrm{i}}$ and $e_{\mathrm{e}}$ with the symmetric formulation. In contrast, although the energy equation of Eq. (4) is discretized directly by the conservative scheme, $e_{\mathrm{i}}$ is indirectly obtained from the other discretized equations; this explains why the existing conservative scheme violates the law of equipartition. The derivation of the proposed scheme is unique; the mathematical structures of the governing equations must be preserved in order to maintain the important physical principles.

5 Accuracy verification

The spatial accuracy is assessed to check whether the formal accuracy is reproduced by the proposed scheme. The 1F2T hydrodynamic equations are decomposed into an eigenstructure comprised of two pressure waves and the entropy waves of the ions and electrons. We introduce this test problem to examine the spatial accuracy of the 1F2T model:

[TABLE]

where $x\in[0,1]$. The initial profiles are depicted in Fig. 4. Again, the EOS is thermally and calorically ideal. These conditions eliminate the pressure waves which correspond to the nonlinear field. Therefore, this is a linearized problem involving ion and electron entropy waves, whose solutions are easily obtained via the hyperbolic solvers. A periodic condition is applied to the boundaries, so that the initial conditions are recovered at $t=1$. The time interval $\Delta t=10^{-5}$ is sufficiently small to make the temporal error negligible.

Figure 5 is a log-log graph to illustrate the accuracy verification. The results show that $\rho$, $e_{\mathrm{i}}$, and $e_{\mathrm{e}}$ have the fourth-order formal accuracy of the compact scheme (Eq. (29)). Note that $e_{\mathrm{i}}$ and $e_{\mathrm{e}}$ are not discretized separately; the discretizations are only applied to the conservation laws of mass, momentum, and energy. On the other hand, the spatial accuracy of $u$ cannot be observed via this numerical experiment, because $u$ and the static pressure $p_{\mathrm{i}}+p_{\mathrm{e}}$ are initially constant in the entire region. Equation (16) is simplified by assuming $u^{n}_{j}=U=\mathrm{const.}$ and $p^{n}_{j}=\mathrm{const.}$, such that

[TABLE]

Here, the relationship $\langle\rho u^{2}\rangle=U\langle\rho u\rangle$ is valid at the discrete level, because of the linearity of the interpolation operators. Hence, $u^{n+1}_{j}=u^{n}_{j}$ is obtained by substituting Eq. (15) into Eq. (55). This is why the $u$ errors are independent of the grid interval, so that the error-norm always maintains the round-off level.

Note that the discretizations are not based on the FVM but, rather, on the FDM, although FVM-like schemes of interpolation and filtering are utilized in this investigation.

6 Conclusions

In this article, a structure-preserving scheme for the 1F2T hydrodynamic equations is proposed, with the aim of improving the reliability of compressible RHD simulations. The proposed scheme exactly satisfies the conservation laws and thus the Rankine–Hugoniot relationship.

The key of constructing a physically accurate scheme that satisfies the important physical principles is to maintain the mathematical structure of the governing equations, even in the discrete form. Specifically, the product rule and the symmetry of the energy equations of each species must be maintained. Therefore, the proposed approach does not discretize the energy equations of the ions and electrons directly but, rather, discretizes the energy conservation law using the FDM approach. The ion and electron energy equations in the discrete form are derived using the product rule; this is the same strategy as that used to derive the nonconservative equations in differential form. This derivation yields error terms in the energy equations, which should be separated equally in accordance with the law of equipartition.

Verification via the shock tube problem demonstrates that the proposed scheme maintains the global conservation error to within the round-off level and well agrees with the Rankine–Hugoniot relationship of the 1F2T model. In other words, the proposed scheme strictly preserves the conservation laws of mass, momentum, and energy, and the law of equipartition in the discrete form. Moreover, accuracy verification based on the linear advection of entropy waves reveals that the proposed scheme yields the formal accuracy.

Although the proposed scheme possesses the prefer features explained above, some issues remain toward practical RHD simulations, such as those considering the ICF implosion. For example, the scheme must be generalized to curvilinear coordinates for spherical[34]/cylindrical[35] implosions, non-ideal EOS, and magnetohydrodynamics (MHD) for magnetized fast ignition [36, 37]. This problem may be solved by our strategy using a previous work about MHD scheme [38]. Furthermore, our approach may be used in hypersonic hydrodynamics of re-entries. It is modeled by multitemperature hydrodynamics regarding the translational, rotational, vibrational and excitation modes [39, 40].

Acknowledgments

This work was supported by a Grant-in-Aid from the Japan Society for the Promotion of Science (JSPS) Fellows, No. 15J02622. T.S. wishes to thank Dr. Atsushi Sunahara (Purdue University) for valuable discussions on the physical background of the 1F2T model.

Appendix A Usage of Runge–Kutta method in proposed scheme

High-order time integration is important for the performance of high-resolution simulations. Here, the implementations of the RK methods in the proposed scheme are presented.

A.1 First-order RK method

The first-order RK method is identical to the Euler explicit method. The terms in Eq. (21) can be classified into the following three components:

[TABLE]

where $ADV$, $DFS$, and $ERR$ represent the advection, artificial dissipation, and error components, respectively.

A.2 Second-order RK method

Here, the proposed scheme is extended using the second-order RK method. The conservation laws of mass, momentum, and energy are discretized as follows:

[TABLE]

where “$*$” denotes the internal time-step of the second-order RK method. Equations (60)–(62) and (63)–(65) correspond to the primary and secondary RK steps, respectively. Obviously, the following energy equations are obtained from Eq. (60)–(62):

[TABLE]

In addition, the spatial-difference terms of Eq. (63)–(65) can be interpreted as arithmetic averages of “$n$” and “$*$” steps. Thus, the energy equations in the second step are expressed as

[TABLE]

A.3 Third-order TVD RK method

The implementation of the third-order TVD RK method is derived using a similar method to the second-order case. The first step is expressed as

[TABLE]

where “$\dagger$” indicates the primary internal time-step of the third-order TVD RK. The second step is given as

[TABLE]

where “$\dagger\dagger$” indicates the secondary internal time-step of the third-order TVD RK. The final step is obtained by weighting these internal steps, such that

[TABLE]

Appendix B Multidimensional description

Multidimensional scheme is required for the numerical simulations of ICF implosions because hydrodynamic instabilities such as Rayleigh–Taylor instability [41, 42] are one of the fundamental physics to determine the fusion gain. Here, we introduce a multidimensional description of the proposed approach. Before the derivation of multidimensional scheme, we define two vector differential operators in the discrete form:

[TABLE]

where $\phi$ is an arbitrary scalar function, $\mathbf{\Psi}=^{\mathrm{T}}[\psi_{x},\psi_{y},\psi_{z}]$ is an arbitrary vector function, and $(i,j,k)$ are the spatial indices over $(x,y,z)$ directions, respectively. The three-dimensional Euler equation is discretized by using these operators as follows:

[TABLE]

Note that the key of our approach in Sec. 2.3 is the expansion of the time derivative. Hence, the three-dimensional scheme is also obtained by the same way. The discretized equations corresponding to Eqs. (3) and (4) are

[TABLE]

Extension to the shock capturing scheme is self-evident.