Solving the Riemann Problem for Realistic Astrophysical Fluids

TL;DR

This paper introduces new methods for solving the Riemann problem with realistic equations of state in astrophysical fluids, enabling more accurate and efficient modeling of complex astrophysical flows.

Contribution

It develops both exact and approximate Riemann solvers for general EoS, guaranteeing solution existence and uniqueness under monotonicity, and compares two thermodynamic calculation methods.

Findings

Exact and approximate Riemann solutions match for perfect gas EoS.

Analytic EoS solver is faster than interpolation-based methods.

Numerical results show similar accuracy for both EoS solving methods.

Abstract

We present new methods to solve the Riemann problem both exactly and approximately for general equations of state (EoS) to facilitate realistic modeling and understanding of astrophysical flows. The existence and uniqueness of the new exact general EoS Riemann solution can be guaranteed if the EoS is monotone regardless of the physical validity of the EoS. We confirm that: (1) the solution of the new exact general EoS Riemann solver and the solution of the original exact Riemann solver match when calculating perfect gas Euler equations; (2) the solution of the new Harten-Lax-van Leer-Contact (HLLC) general EoS Riemann solver and the solution of the original HLLC Riemann solver match when working with perfect gas EoS; and (3) the solution of the new HLLC general EoS Riemann solver approaches the new exact solution. We solve the EoS with two methods, one is to interpolate 2D EoS tables by the bi-linear interpolation method, and the other is to analytically calculate thermodynamic variables at run-time. The interpolation method is more general as it can work with other monotone and realistic EoS while the analytic EoS solver introduced here works with a relatively idealized EoS. Numerical results confirm that the accuracy of the two EoS solvers is similar. We study the efficiency of these two methods with the HLLC general EoS Riemann solver and find that analytic EoS solver is faster in the test problems. However, we point out that a combination of the two EoS solvers may become favorable in some specific problems. Throughout this research, we assume local thermal equilibrium.

Tables

Table 1. Table 1: Test 1 is a shock tube test. We choose the ratio of specific heat γ = 1.05 𝛾 1.05 \gamma=1.05 and 1.667 1.667 1.667 while keep the mean atomic weight μ = 1.3 𝜇 1.3 \mu=1.3 . The density, velocity and temperature of the left state are listed from the second column to the fourth column while the rightest three columns list the variables of the right state.

ID	$t_{\max }$	${\rho }_L$	$u_L$	$T_L$	${\rho }_R$	$u_R$	$T_R$
	[s]	[g$\cdot$cm-3]	[km$\cdot$s-1]	[K]	[g$\cdot$cm-3]	[km$\cdot$s-1]	[K]
1	0.02	${10}^{-13}$	$0$	$3000$	${10}^{-15}$	$0$	$300$

Table 2. Table 2: Simulation run-times ( t max subscript 𝑡 max t_{\rm{max}} ) and initial left and right states of the Riemann problem tests. μ = 1.3 𝜇 1.3 \mu=1.3 .

ID	$t_{\max }$	${\rho }_L$	$v_L$	$T_L$	${\rho }_R$	$v_R$	$T_R$
	[s]	[g$\cdot$cm-3]	[km$\cdot$s-1]	[K]	[g$\cdot$cm-3]	[km$\cdot$s-1]	[K]
2	0.012	${10}^{-13}$	0	$20000$	${10}^{-15}$	0	$300$

Table 3. Table 3: Run-times ( t max subscript 𝑡 max t_{\rm max} ) and initial left and right states of the four tests with pure hydrogen gas. The tests evolve a strong-shock, asymmetric two-shock, symmetric two-rarefaction, and asymmetric two-shock-wave respectively. More detailed descriptions are in the text.

ID	$t_{\max }$	${\rho }_L$	$v_L$	$T_L$	${\rho }_R$	$v_R$	$T_R$
	[s]	[g$\cdot$cm-3]	[km$\cdot$s-1]	[K]	[g$\cdot$cm-3]	[km$\cdot$s-1]	[K]
2	0.012	${10}^{-13}$	0	$20000$	${10}^{-15}$	0	$300$
3	0.100	$2\times {10}^{-15}$	40	$800$	${10}^{-15}$	-60	$800$
4	0.008	$2\times {10}^{-11}$	-30	$15000$	$2\times {10}^{-11}$	30	$15000$
5	0.070	$1.2\times {10}^{-15}$	55	$500$	${10}^{-15}$	-65	$500$

Table 4. Table 4: CPU time of the three HLLC Riemann solvers of test 2-5 in Table 3 in second. We show the CPU time of the original HLLC Riemann solvers, the HLLC general EoS Riemann solver with interpolation method, and the HLLC general EoS Riemann solver with analytic EoS solver. The ratio of the CPU time of the two HLLC general EoS Riemann solvers to the original HLLC Riemann solver are shown below each solver.

	test2	test3	test4	test5
$t_o$ [s]	$2.250\times {10}^{-1}$	$2.546\times {10}^0$	$1.750\times {10}^{-1}$	$2.186\times {10}^0$
$t_i$ [s]	$1.729\times {10}^1$	$8.849\times {10}^1$	$2.318\times {10}^1$	$6.175\times {10}^1$
speed	76.84	34.76	132.5	28.25
$t_a$ [s]	$5.472\times {10}^0$	$5.335\times {10}^1$	$1.022\times {10}^1$	$2.846\times {10}^1$
speed	24.32	20.95	58.40	13.01

Equations

$$ϵ=\frac{p}{\rho (\gamma -1)},$$

$$a=\sqrt{\Gamma p/\rho },$$

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

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

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

$$E$$

$$=\frac{\rho u^2}{2}+e_{\mathrm{g}},$$

$$ϵ=ϵ(\rho ,p)\text{or}p=p(\rho ,ϵ),$$

$${\left(\frac{\partial p}{\partial \rho }\right)}_{ϵ}>0$$

$${\left(\frac{\partial p}{\partial ϵ}\right)}_{\rho }>0$$

$${\mathbf{V}}_t+\mathbf{A}(\mathbf{V}){\mathbf{V}}_x=0,$$

$$\mathbf{V}={\left[\frac{1}{\rho },u,p\right]}^T$$

$$\mathbf{A}(\mathbf{V})=\left[\begin{array}{ccc}u& -\frac{1}{\rho }& 0\\ 0& u& \frac{1}{\rho }\\ 0& \rho a^2& u\end{array}\right],$$

$${\mathbf{l}}^1\cdot d\mathbf{V}=\frac{dp}{\rho a}-du=0,$$

$${\mathbf{l}}^2\cdot d\mathbf{V}=\frac{dp}{(\rho a{)}^2}+d\left(\frac{1}{\rho }\right)=0,$$

$${\mathbf{l}}^3\cdot d\mathbf{V}=\frac{dp}{\rho a}+du=0.$$

$$[u]=\mp \frac{[p]}{J_{\mathrm{S}}},$$

$$\frac{[p]}{J_{\mathrm{S}}^2}=-\left[\frac{1}{\rho }\right],$$

$$[ϵ]=-\frac{p_{\mathrm{S}}+p_{\ast }}{2}\left[\frac{1}{\rho }\right],$$

$$J_{\mathrm{S}}^2[ϵ]=\frac{1}{2}[p^2].$$

$$p_{\mathrm{exact}}({\mathbf{W}}_{\mathrm{L}},{\mathbf{W}}_{\mathrm{R}},\Gamma )$$

$$p_{\min \ast }$$

$$p_{\max \ast }$$

$$J_{\mathrm{S},\min }=0,$$

$$J_{\mathrm{S},\max }=10\times {\rho }_{\bar{\mathrm{S}}}(|v_{\bar{\mathrm{S}}}|+a_{\bar{\mathrm{S}}}),$$

$$\Delta u(p_{\ast })=f_{\mathrm{L}}(p_{\ast })-f_{\mathrm{R}}(p_{\ast })$$

$$ϵ$$

$$\frac{1}{{\rho }_{\mathrm{S}\ast }}$$

$$J_{\mathrm{S}}^2[ϵ]$$

$$W\left({ϵ}_{\mathrm{S}\ast }-ϵ\right)=\frac{p_{\ast }^2-p^2}{2},$$

$$\frac{d\text{LHS}}{dW}={ϵ}_{\ast }(W)-ϵ+W{\left(\frac{\partial {ϵ}_{\mathrm{S}\ast }}{\partial {\rho }_{\mathrm{S}\ast }}\right)}_{p_{\ast }}\frac{d{\rho }_{\mathrm{S}\ast }}{dW},$$

$$\frac{d{\rho }_{\mathrm{S}\ast }}{dW}=-{\rho }_{\mathrm{S}\ast }^2\frac{p_{\ast }-p}{W^2}<0,$$

$$\displaystyle\left(\frac{\partial\epsilon}{\partial\rho}\right)_{p}=-\left(\frac{\partial\epsilon}{\partial p}\right)_{\rho}\bigg{/}\left(\frac{\partial\rho}{\partial p}\right)_{\epsilon}=\epsilon_{\rho}\equiv\epsilon_{\rho_{\rm S*}}<0.$$

$$J_{\mathrm{S},\max }^2-J_{\mathrm{S},\min }^2<2{\delta }_{\mathrm{machine}}$$

$$p_{\max \ast }-p_{\min \ast }<2{\delta }_{\mathrm{machine}},$$

$${\lambda }_{\mathrm{S}}$$

Full text

††journal: Journal of Computational Physics

1 Introduction

The need to incorporate complex equations of state (EoS) in hydrodynamics is becoming increasingly important in computational astrophysics. The state-of-the-art for many numerically intensive research problems in such diverse sub-fields as accretion disks, binary mergers, star formation, novae, and supernovae, is approaching the state where more accurate EoS relating fluid state variables (e.g., density, pressure, and temperature) are needed. Although the assumption of a perfect gas EoS with a fixed constant ratio of specific heats is now prevalent in simulations with Riemann solvers, such constancy is usually not fully justified, especially where ionization or phase transitions occur (e.g., hydrogen ionization transition). In a perfect gas EoS, the specific internal energy

[TABLE]

where $p,\rho$, and $\gamma$ are the pressure, density and ratio of the specific heats (i.e. $\gamma=C_{p}/C_{V}$), respectively. $C_{p}$ is the specific heat capacity at constant pressure and $C_{V}$ is the specific heat capacity at constant volume. The adiabatic sound speed $a$ of a general fluid can be calculated as

[TABLE]

where $\Gamma=\left(\partial\ln p/\partial\ln\rho\right)_{s}$ and the subscript $s$ denotes that this derivative is taken at constant specific entropy. If the fluid is a perfect gas $\Gamma=\gamma$. For a realistic gas $\gamma$ and $\Gamma$ may not be the same (see Appendix A for details).

A further complication for general EoS is binding/ionization energy. For example, it takes $13.6$ eV to ionize a hydrogen atom, and the ionization happens around $T\approx 8000$ K (dependent on $\rho$), corresponding to an average thermal kinetic energy of $\approx 1.0$ eV. This makes the binding energy significantly larger than the translational energy. Similar situations occur during the disassociation of molecular hydrogen and the ionization of \chHe and \chHe^+. Neglecting these binding energies and other endothermic processes would lead to wrong gas temperature, causing significant errors in radiative transfer and chemical/nuclear reaction calculations. The hydrogen ionization instability powering outbursts observed in the accretion disks of white dwarf systems called dwarf novae [1] is one prominent example. Refs. [2, 3, 4] show that using an EoS, which accurately models the hydrogen ionization transition is essential for understanding dwarf novae.

The 1D general EoS hydrodynamic Euler equations that we are interested in solving are,

[TABLE]

where

[TABLE]

and an EoS is required to close the system of equations

[TABLE]

where $\rho,u,p,E,e_{\rm g}$ and $\epsilon$ are the density, velocity, pressure, total energy density, internal energy, and specific internal energy, respectively. Local thermal equilibrium (LTE) is a requirement for all hydrodynamic models, and we assume the LTE condition is satisfied throughout this paper. As noted before, $\Gamma=a^{2}\rho/p$ where $a$ is the adiabatic sound speed. We will use $\Gamma$ in our general EoS Riemann solvers.

Some recent work has aimed to devise efficient, high fidelity general EoS Riemann solvers [5, 6] with the EoS in assumed analytic form. Nonetheless, it is vital to test the accuracy of such general EoS Riemann solvers. The idea of general EoS Riemann solver of the Harten-Lax-van Leer-Contact (HLLC) type has been mentioned by E. F. Toro in [7] (the last remark in Section 10.4.2, Chapter 10), however, there has been a lack of systematic assessment of the accuracy and efficiency of the HLLC type general EoS Riemann solvers.

In this paper, we devise and examine a new efficient HLLC [8] general EoS Riemann solver that can handle realistic gases whose $C_{V}$, $C_{p}$ and $\Gamma$ may vary drastically when the temperature changes. The rapidly varying thermal dynamic variables may become stiff to the Euler equations. Therefore an approximate Riemann solver cannot stand alone. To justify the accuracy of the HLLC general EoS Riemann solver, we design and program a new exact general EoS Riemann solver that can work with non-convex EoS and compare the solution of the HLLC general EoS Riemann solver to the exact Riemann solver. We point out that the exact general EoS Riemann solver is the true solution with the EoS we use. Our exact general EoS Riemann solver only needs two monotone conditions,

[TABLE]

discussed in [9, 10]. An additional convex condition $(\partial^{2}p/\partial\rho^{2})_{s}>0$ was required by [11], however we find that this is not a requirement of the new exact general EoS Riemann solver discussed in this paper (see Section 2.3 for further discussion).

In Section 2, we present the exact general EoS Riemann solver. We introduce the application of the exact general EoS Riemann solver to the Godunov scheme [12] in Section 3. The novel HLLC general EoS Riemann solver will be discussed in Section 4. Section 5 will focus on the numerical results generated by the exact and the HLLC general EoS Riemann solvers. In particular, we will establish the correspondence of the exact general EoS Riemann solver to the perfect gas Riemann solver, and the HLLC general EoS Riemann solver to the original HLLC Riemann solver. We obtain the solutions from the two general EoS Riemann solvers with a Godunov scheme [12] and compare them to the exact solution. Section 5 also contains our efficiency and convergence study of the HLLC general EoS Riemann solver. We conclude in Section 6.

2 The exact general EoS Riemann solver

Much of this section is a review of the seminal work [11], however we utilize a more robust root finding algorithm to calculate the correct middle state pressure $p_{\rm*}$ in Section 2.3, and numerically integrate the ODEs describing the simple/rarefaction waves. We present the derivation here briefly for completeness. Readers are encouraged to read the original paper [11].

2.1 The quasi-linear hyperbolic system in non-conservative form

The one-dimensional Euler equations can be written in the non-conservative form,

[TABLE]

where,

[TABLE]

[TABLE]

and $a$ is the adiabatic sound speed. The (Jacobian) matrix $\mathbf{A}$ has the left and right eigenvectors ($\mathbf{l}^{1},\mathbf{r}^{1}$),($\mathbf{l}^{2},\mathbf{r}^{2}$), and ($\mathbf{l}^{3},\mathbf{r}^{3}$), which correspond to the left, middle (contact discontinuity) and right waves respectively. Figure 1 shows the solution structure of a rightward shock-tube Riemann problem. A key condition in [11] is that $p_{\rm L*}=p_{\rm R*}=p_{\rm*}$ across the contact discontinuity, as previously noted by others including [13]. In appendix B, we prove the veracity of this condition for the general EoS Euler equations.

The characteristic equations are obtained by setting $\mathbf{l}^{i}\cdot d\mathbf{V}=0$ for $i=1,2,3$:

[TABLE]

For the right rarefaction wave, the solutions of $\rho,\ u$ and $p$ are calculated by Equation 14 and 15 while for the left rarefaction wave, the solutions of $\rho,\ u$ and $p$ are calculated by Equation 15 and 16. For the left and right shock waves, the pre-shock and post-shock conditions are determined by the Rankine-Hugoniot relation:

[TABLE]

where $[q]=q_{\rm*}-q_{\rm S}$ and S=L,R. $J_{\rm S}=\mp\rho_{S}(\lambda_{S}-u_{S})$ is the mass flux across the shock wave. The shock front speed is expressed by $\lambda_{\rm S}=u_{\rm S}\mp\frac{J_{\rm S}}{\rho_{\rm S}}$. Substituting Equation 19 into Equation 18 and eliminating $\rho$, we obtain

[TABLE]

In Section 2.3.2, we will use Equation 20 to calculate $J_{\rm S}$.

2.2 Uniqueness of solution and the p-u phase plot

[9, 10] have proved that Equation 9 and 10 are sufficient but not necessary to the uniqueness of solution of Euler equation with a general EoS. We assume the gas we work with satisfy these two conditions. More in depth discussion of the EoS we use can be found in Appendix C which contains a discussion about these two conditions.

In Appendix B, Equation 87 and 89 prove that the contact discontinuity solution of the Euler equations with a general EoS should be $u_{\rm{L*}}=u_{\rm{R*}}=u_{\rm{*}}$ and $p_{\rm{L*}}=p_{\rm{R*}}=p_{\rm{*}}$. This enables us to follow [14] and think about the solution as the intersection of curves in $p\textnormal{-}u$ phase space (See Figure 2). The initial state $p_{\rm S}$ and $u_{\rm S}$ is a point on the $p\textnormal{-}u$ phase plot. The rarefied or shocked state $(p_{\rm{S*}},u_{\rm{S*}})$ can be calculated by Equation 14-19 and generate a continuous curve on the $p\textnormal{-}u$ phase space for each of the two values of $S$ (L, R). The following two propositions will be helpful to the understanding of the $p\textnormal{-}u$ phase plot and the uniqueness of a solution:

Proposition 1**.**

$dp/du<0$* for the left initial state.*

Proof: if $p>p_{\rm L}$, there is a shock wave on the left side, use Equation 17 and we can get $dp/du=-J_{\rm L}<0$. If $p<=p_{\rm L}$, there is a rarefaction wave on the left side, use Equation 16 and we can get $dp/du=-1/(\rho a)<0$.

Proposition 2**.**

$dp/du>0$* for the right initial state.*

Proof: similar to Proposition 1.

Figure 2 shows an example $p\textnormal{-}u$ phase plot of a rightward shock with realistic EoS and a tabulated EoS solver. The initial state is the same as test 2 in Table 2. The left and right initial states are on the p-axis in Figure 2 The intersection of the two curves is marked with a red star and it is the solution of $p_{*}$ and $u_{*}$. Let us name the curve of points that can be connected to the left state via a rarefaction wave be $f_{\rm L}(p)=u_{\rm L}$ and the curve of points that can be connected to the right state via a shock wave be $f_{\rm R}(p)=u_{\rm R}$. The numerical method in Section 2.3 is used to find the value of these two equations and is the main topic of the exact general EoS Riemann solver.

2.3 Numerical method

The more robust root finding algorithm we propose is a two iterations bisection method111Other bracketing root finding techniques (e.g. the Brent-Dekker method) could also be used. solver. In [11], Secant and Newton’s method were proposed to find the solution of $p_{\rm*}$. Secant and Newton’s method have better convergence rate but they do not guarantee a correct solution. For example, when using these two methods to solve $f(x)=0$ and $f^{\prime}(x)>0$ everywhere, if $f^{\prime\prime}(x)<0$ happens near the root, these two methods may not converge to a correct solution. In [11], a convex condition, $(\partial^{2}p/\partial\rho^{2})_{s}>0$, is required to circumvent this situation. However, this condition may not be satisfied at all time in realistic gas, specifically, during the phase transition, making the secant and Newton’s method incorrect. In Figure 16, we show the adiabatic sound speed $a=\left(\frac{\partial p}{\partial\rho}\right)_{s}$ of a realistic gas and emphasize the sound speed change during two phase transitions to exemplify the complexity of $a$ during phase transitions.

In the following of this section, we divide the procedure for finding the middle state pressure $p_{\rm*}$ in general EoS Riemann problem into five major steps.

2.3.1 Estimation of $p_{\rm min*}$ and $p_{\rm max*}$

There is no unique way of determining the bracketing pressure values $p_{\rm min*}$ and $p_{\rm max*}$; here we describe our procedure. First, we need to know the maximum and minimum $\Gamma$ of the EoS table. We employ the exact perfect gas iterative Riemann solver described in Chapter 4 of [7] to estimate the lower and upper bound of $p_{\rm*}$. For simplicity, the solution of the exact perfect gas iterative Riemann solver is represented by a function as:

[TABLE]

Where $\mathbf{W}_{\rm S}=[\rho_{\rm S},u_{\rm S},p_{\rm S}]^{T}$ is the vector of primitive (fluid) variables of the left or right state. The readers should note that Equation 21 is actually a procedure with logical control and iteration. Equation 21 gives the exact $p_{\rm*}$ of the middle states with the given left, right states and the $\Gamma$ for perfect gas. We are slightly abusing this procedure here because $\Gamma\neq\gamma$ for realistic gas. However, we only use them for an initial guess of the lower and upper bound in which $p_{\rm*}$ may reside in and it turns out to be useful. The lower and upper bound of $p_{\rm*}$ are estimated by,

[TABLE]

We multiply or divide the guessed pressure by a factor of two to ensure that the solution is indeed bracketed.

2.3.2 Estimate $J_{\rm S,min}$ and $J_{\rm S,max}$ for $p_{\rm*}>p_{\rm S}$ (shock present)

When $p_{\rm*}>p_{\rm S}$, there is a left or right shock. For each estimate of $p_{\rm*}>p_{\rm S}$ in the range $[p_{\rm min*},p_{\rm max*}]$ during the $p_{*}$ bisection iteration (see Section 2.3.3), we estimate the corresponding range of $J_{\rm S}$ ($[J_{\rm S,min}$, $J_{\rm S,max}]$) and initiate the second bisection loop. We use the chosen $p_{\rm*}$ to estimate $J_{\rm S,min}$ as

[TABLE]

since the mass flux could be close to [math] according to the definition of $J_{\rm S}$, and

[TABLE]

where $\bar{\rm S}$ is the opposite state of $\rm S$ and $a_{\rm\bar{S}}$ is the adiabatic sound speed of the $\bar{S}$ state. The choice of a large $J_{\rm S,max}$ is heuristic because we are presently interested only in the correctness of the exact general EoS Riemann solver, not its efficiency. This means that we only require that the root ($J_{\rm S}$) is bracketed. Practically, one may not be able to apply the exact general EoS Riemann solver to 2D and 3D calculations even with an optimized algorithm.

2.3.3 Using the bisection method to determine $p_{\rm*}$ and $J_{\rm S}$

In this section, we show that our bisection method ($J_{\rm S}$ iteration) guarantees to find one and only one solution of $J_{\rm S}$.

Let $u_{\rm L*}$ and $u_{\rm R*}$ be the velocity at the left and right side of the contact discontinuity calculated with any given $p_{\rm*}$. If $p_{\rm*}$ exists within the range $[p_{\rm min*},p_{\rm max*}]$, then the expression $f(p_{\rm*})=u_{\rm L*}-u_{\rm R*}$, should have different signs at $p_{\rm min*}$ and $p_{\rm max*}$, otherwise the solution has not been properly bracketed. $\Delta u$ can be defined as,

[TABLE]

where $f_{\rm L}$ corresponds to the curve of left state in the $p\textnormal{-}u$ phase space and $f_{\rm R}$ corresponds to the curve of right state in the $p\textnormal{-}u$ phase space (Figure 2). Put another way, inserting the correct $p_{\rm*}$ will yield $u_{\rm L*}=u_{\rm R*}$ and $\Delta u(p_{\rm*})=0$. We find the correct $p_{\rm*}$ in the range $[p_{\rm min*},p_{\rm max*}]$ by first picking a $p_{\rm*}$ in the range $[p_{\rm min*},p_{\rm max*}]$.

If $p_{\rm*}>p_{\rm S}$, then there is a shock wave on the $S$ side. Find the correct $J_{\rm S}$ in the range of $[J_{\rm S,min},J_{\rm S,max}]$ using the bisection method. The correct $J_{\rm S}$ should satisfy Equations 8, 18, and 20 simultaneously. For convenience, we list and rearrange the equations here again,

[TABLE]

Equation 28 relates $\rho_{\rm S*}$ to $J_{\rm S}^{2}$. Substitute Equations 27 and 28 into Equation 29, to obtain a single equation for $J_{\rm S}^{2}$. For convenience, let us define $W\equiv J_{\rm S}^{2}\geq 0$. Then

[TABLE]

where $\epsilon_{\rm S*}=\epsilon_{\rm*}(W)$ is implicitly a function of $W$ (via Equations 27-29). The right hand side of Equation 30 is a positive constant in the iteration of $J_{\rm S}$. The derivative of the left hand side is

[TABLE]

from Equation 28

[TABLE]

and from Equation 9 and 10,

[TABLE]

Thus, $d\text{LHS}/dW\geq 0$ and only equals 0 when $W=J_{\rm S}^{2}=0$, $\text{LHS}=0$ when $W=J_{\rm S}^{2}=0$, so there is exactly one solution of $W$ to Equation 30. In Section 2.3.2, we have chosen the $J_{\rm S,max}$ to be large enough to ensure the solution $J_{\rm S}\leq J_{\rm S,max}$.

To apply the bisection root finding algorithm to Equation 30, we pick a $J_{\rm S}$ in the range of $[J_{\rm S,min},J_{\rm S,max}]$. We then calculate $\rho_{\rm*}$ using Equation 28 and calculate $\epsilon_{\rm*}$ using Equation 27 for the choice of $J_{\rm S}$ and substitute these results into Equation 30. If LHS$\neq$RHS, we must iterate the bounds of $J_{\rm S}$ until the bracketing range narrows in on a solution. A solution of $J_{\rm S}$ is thought to be found if:

[TABLE]

where $J_{\rm S,max}$ and $J_{\rm S,min}$ are the updated upper and lower bound of the bracket. $\delta_{\rm{machine}}$ is the machine precision unit.

After a proper $J_{\rm S}$ has been found, $u_{\rm{S*}}$ can be calculated by Equation 17.

If $p_{\rm*}\leq p_{\rm S}$, there is a rarefaction wave on the $S$ side. We then simply solve the ODEs 14 and 15 for a right rarefaction wave and 15 and 16 for a left rarefaction wave. The definite integration starts at $p_{\rm S}$ and ends at $p_{\rm*}$. We then calculate $u_{\rm S*}$.

After obtaining $u_{\rm L*}$ and $u_{\rm R*}$ with the given $p_{\rm*}$, we compare $u_{\rm L*}-u_{\rm R*}$ to the value of 0, update $p_{\rm min*}$ or $p_{\rm max*}$ and increment $p_{\rm*}$ according to bisection method.

The iteration is stopped when $p_{\rm min*}$ and $p_{\rm max*}$ are sufficiently close. In practice, we define “sufficiently close” as,

[TABLE]

When condition 35 is satisfied, $\left|u_{\rm{R*}}-u_{\rm{L*}}\right|$ may not be less than $2\delta_{\rm{machine}}$. Different problems may also give different gaps. This is because $p_{*}$ and $u_{*}$ are in general converging at different rate, it is impractical to find a constant $\delta$ and force their gaps to become smaller than $\delta$ at the same iteration. In addition, the nature of bisection method does not allow further numerical calculation if condition 35 is satisfied therefore we would just stick to this convergence criterion. Calculation shows that $\left|u_{\rm{R*}}-u_{\rm{L*}}\right|<10^{-5}$ cm$\cdot$s*-1* for all of our tests and the gap decreases as $\left|p_{\rm max*}-p_{\rm min*}\right|$ decreases.

2.3.4 Calculation of $\lambda_{\rm S}$ and $u_{\rm*}$

To calculate the three wave speeds discussed in Section 2.1 for a shock wave, we need to utilize the $J_{\rm S}$ computed in the previous section:

[TABLE]

If there is a rarefaction wave, then

[TABLE]

where $a_{\rm S}$ is the adiabatic sound speed of the S state. Here $u_{\rm*}$ can be calculated by solving ODEs 14 and 15 for the right rarefaction wave and ODEs 15 and 16 for the left rarefaction wave, where the bounds of the integration are $p_{\rm S}$ and $p_{\rm*}$. In any given Riemann problem, the value $u_{\rm*}$ calculated by the left and right waves must be equal (as proved in Appendix B), otherwise the previous steps of our solution procedure have not been followed properly (e.g. finding $p_{\rm*}$ or $J_{\rm S}$).

2.3.5 Sample results

The solution ($\mathbf{W}(x,T)$) of a Riemann problem could be one of the three types: (i) two shock waves; (ii) two rarefaction waves: (iii) and one shock wave with one rarefaction wave. As before the left wave speed is $\lambda_{\rm L}$ and the right wave speed is $\lambda_{\rm R}$. The ratio of coordinate $x$ to given time $t$ is related to the different waves and state vectors $\mathbf{W}(x,T)=\mathbf{W}$ as follows:

For the left shock wave,

If $x/t<\lambda_{\rm L}$, then $\mathbf{W}=\mathbf{W}_{\rm L}$

• 2.

If $\lambda_{\rm L}\leq x/t<u_{\rm*}$, then $\mathbf{W}=\mathbf{W}_{\rm L*}$

For the left rarefaction wave,

If $x/t<\lambda_{\rm L}$, then $\mathbf{W}=\mathbf{W}_{\rm L}$

• 2.

If $\lambda_{\rm L}\leq x/t<u_{\rm*}-a_{\rm*}$, integrate the ODEs 15 and 16 until $x/t=u-a$ to get the correct $\mathbf{W}$ in the fan region.

• 3.

If $u_{\rm*}-a_{\rm*}\leq x/t<u_{\rm*}$, then $\mathbf{W}=\mathbf{W}_{\rm L*}$

For the right shock wave,

If $\lambda_{\rm R}<x/t$, then $\mathbf{W}=\mathbf{W}_{\rm R}$

• 2.

If $u_{\rm*}<x/t\leq\lambda_{\rm R}$, then $\mathbf{W}=\mathbf{W}_{\rm R*}$

For the right rarefaction wave,

If $\lambda_{\rm R}<x/t$, then $\mathbf{W}=\mathbf{W}_{\rm R}$

• 2.

If $u_{\rm*}+a_{\rm*}\leq x/t<\lambda_{\rm R}$, integrate 14 and 15 until $x/t=u+a$ to get the correct $\mathbf{W}$ in the fan region.

• 3.

If $u_{\rm*}\leq x/t<u_{\rm*}+a_{\rm*}$, then $\mathbf{W}=\mathbf{W}_{\rm R*}$

We present a flow chart of the exact general EoS Riemann solver in Figure 3 to make it easier to follow.

2.4 Tabulated EoS solver for the exact general EoS Riemann solver

In the exact general EoS Riemann solver, we constantly need to solve

[TABLE]

in Equation 14-19. Therefore, we prepare tables with of $\log a(\log\rho,\log p)$ and $\log\epsilon(\log\rho,\log p)$ and use bi-linear interpolation method to interpolate the tables to get the value of $a$ and $\epsilon$. Bi-linear interpolation method is easy to implement and also preserve the monotonicity of the table. The tables we use are rectangular and the domain of the tables are $\log\rho\in[-15,-10]$ g$\cdot$cm*-3and $\log p\in[-5.38,2.62]$ erg$\cdot$cm-3*. The resolution of the domain is $501\times 701$ including the endpoints.

We can calculate $\epsilon,\ T$ and $\mu$ by interpolating $\log\epsilon(\log\rho,\log p)$, $\log T(\log\rho,\log p)$, and $\mu(\log\rho,\log p)$ tables with bi-linear interpolation method. The number density $n_{i}$ of all species can be calculated after getting the value of $\mu$.

3 Applying the exact general EoS Riemann solver to the Godunov scheme

In this section, we discuss how to apply the exact general EoS Riemann solver to the Godunov scheme to first order. For an introduction to the Godunov scheme, we refer readers to standard computational fluid mechanics books, for example [15, 7]. The quantities that are needed in the Godunov scheme are the flux and the fastest wave speed at the interface. In Section 2, we have solved the exact general EoS Riemann problem. We are able to calculate the $\mathbf{W}$ at any point and get the nonlinear wave speeds. We can easily compare the wave speeds and pick the fastest one. To calculate the flux, we should set $dx/dt=x/t=0$ in Section 2.3.5 and calculate the corresponding $\mathbf{W}$ and the flux through the flux function

[TABLE]

4 The Harten-Lax-van Leer-Contact (HLLC) general EoS Riemann solver

4.1 A brief review of the original HLLC Riemann solver

We assume that readers are familiar with the original HLLC Riemann solver [8] and only briefly list the main steps here. The original HLLC Riemann solver differs from the HLL [16] Riemann solver by dividing the middle region into two states separated by the middle wave (i.e. the contact discontinuity). Figure 4 shows a solution structure of the original HLLC Riemann solver, where $\lambda_{\rm L}$ and $\lambda_{\rm R}$ are the speeds of the left and the right waves, and $u_{\rm*}$ is the middle wave speed. The main difference between the exact and HLLC Riemann Solvers is that HLLC assumes that the left and right waves are always shocks. The wave speeds $\lambda_{\rm L}$ and $\lambda_{\rm R}$ need to be calculated before calculating the middle wave speed:

[TABLE]

where

[TABLE]

Equation 42-44 is referred as the “hybrid estimate” in [8]. We briefly describe the estimation of $p_{*}$ in Section 4.2.1.

Additional constraints are placed by the the Rankine-Hugoniot conditions (Equations 17-19) which can be rewritten as follows (Section 10.4.1, Chapter 10 of [7]):

[TABLE]

where $\mathbf{U}_{\rm S}=[\rho_{\rm S},\rho_{\rm S}u_{\rm S},E_{\rm S}]^{T}$ are the conservative quantities, and $\lambda_{\rm L},\lambda_{\rm R}$ and $u_{\rm*}$ are three waves speeds as before. The original HLLC Riemann solver also reintroduces the contact discontinuity

[TABLE]

to HLL type Riemann solvers. With the addition of these constraints, the middle wave speed can be calculated

[TABLE]

enabling the computation of the two middle states,

[TABLE]

4.2 Algorithms of the HLLC general EoS Riemann solver

Ideally, the HLLC general EoS Riemann solver should reduce to the original HLLC Riemann solver when supplied with EoS tables of perfect gas. Consequently, we attempted to minimize the changes when expanding the HLLC method to general EoS. For completeness, we also outline the whole procedure here.

4.2.1 Estimate $p_{*}$

Unlike the perfect gas, realistic gas has variable $\Gamma$. When the $\Gamma_{\rm L}$ of the left and $\Gamma_{\rm R}$ of the right state are very different, it may be natural to use one of the $\Gamma$. We have experimented with $\Gamma=\min{(\Gamma_{\rm L},\Gamma_{\rm R})}$, $\Gamma=\max{(\Gamma_{\rm L},\Gamma_{\rm R})}$, $\Gamma=(\Gamma_{\rm L}+\Gamma_{\rm R})/2$, the upwind or downwind $\Gamma$, and both $\Gamma_{\rm L}$ and $\Gamma_{\rm R}$. We find the difference in them is not big in terms of the solution of $\rho,\ u,\ p,\ T,\ e_{\rm g}$, etc. The reason for this could be that the solution of $\Gamma$ is already different from the exact solution (see Figure 8-11) because $\Gamma$ is very sensitive to $\rho$ and $T$. If $\Gamma$ is not accurate, any manipulation of $\Gamma$ would not be effective. In this paper, we use both $\Gamma$ to calculate $p_{*}$. The algorithm that we use to calculate $p_{*}$ is a small variation of the adaptive non-iterative Riemann solver in Section 9.5.2, Chapter 9 of [7]. It consists of three sub-solvers: the primitive variable, two-rarefaction, and two-shock Riemann solvers. We substitute all the sound speed with the adiabatic sound speed calculated from the EoS solvers, which means we are using both $\Gamma_{\rm L}$ and $\Gamma_{\rm R}$. In the two-rarefaction Riemann solver, one $\Gamma$ is still needed, in this situation, we use $\Gamma=\min(\Gamma_{\rm L},\Gamma_{\rm R})$.

In the above process of estimation, the only possible values are $\Gamma_{\rm L}$ and $\Gamma_{\rm R}$. Therefore, in the case of a perfect gas, the estimated $\Gamma$ will always be the constant value (i.e., $\Gamma=\Gamma_{\rm L}=\Gamma_{\rm R}=\gamma$). This guarantees that the HLLC general EoS Riemann solver can be reduced to the original HLLC Riemann solver when fed with an EoS table with constant $\Gamma$. We will demonstrate this with two numerical examples in Section 5.2.

We choose this adaptive non-iterative Riemann solver because it is simple and fast. We would admit this choice is quite heuristic at this stage in the context of general EoS. Other choices of Riemann solver are also available [17]. We find that using the Roe average Riemann solver [17] gives similar results to the one introduced here, but it is slower in our implementation as we need to solve the EoS to get the sound speed of the middle state.

4.2.2 Calculate $\lambda_{\rm S}$, $u_{\rm*}$, $\mathbf{U}_{\rm S*}$ and the flux

Calculate $p_{\rm*}$ with the procedure described in Section 4.2.1.

• 2.

Calculate $\lambda_{\rm L}$ or $\lambda_{\rm R}$ from Equation 42 and 43. Use the left and right adiabatic sound speeds that are interpolated from the EoS tables or calculated by the analytic EoS solver.

• 3.

Calculate $u_{\rm*}$ from Equation 50 with the given $\lambda_{\rm L}$ and $\lambda_{\rm R}$.

• 4.

The $U_{\rm S*}$ can be calculated by Equation 51.

• 5.

Calculate the intermediate flux from Equation 45 or 47 with all other calculated quantities.

Figure 5 shows the main steps to get the flux at the interface of the HLLC general EoS Riemann solver. After getting the flux, it is straightforward to apply it to any conservative schemes (for example, Godunov scheme) and calculate the conserved quantities. In designing the algorithm of the HLLC general EoS Riemann solver, we intentional keep the change to a low level to achieve a faster speed. The changes in our procedure compared to the original HLLC Riemann solver include an estimation the $\Gamma$ before the calculations and interpolation of the adiabatic sound speed. These relatively small changes result in numerically accurate solutions as we show if the following section.

4.3 EoS solvers for the HLLC general EoS Riemann solver

We have implemented two types of EoS solvers for the HLLC type of general EoS Riemann solver. The tabulated EoS solver uses the interpolation method described in Section 2.4. In addition, we need to be able to calculate

[TABLE]

Equation 52 is needed in calculating the left and right adiabatic sound speed. Equation 53 will be used to convert the conservative quantities to the primitive quantities. We use bi-linear interpolation method to interpolate $\Gamma(\log\rho,\log p)$ and $\log p(\log\rho,\log\epsilon)$. The domain of the table $(\log\rho,\log\epsilon)$ is $\log\rho\in[-15,-10]$ g$\cdot$cm*-3and $\log\epsilon\in[9.79,13.33]$ erg$\cdot$g-1*with a resolution of $501\times 701$ including endpoints.

The analytic EoS solver is discuss in detail in Appendix C.

5 Numerical results

5.1 The correspondence between the exact general EoS Riemann solver and the exact perfect gas Riemann solver

The perfect gas Riemann solver should be a special case of the general EoS Riemann solver. When using the EoS tables of a perfect gas, the general EoS Riemann solver should give the same solution as the original exact ideal gas Riemann solver. To test this aspect, we create EoS tables of two kinds of perfect gas differ in $\gamma=\Gamma=1.05,1.667$, but the same mean molecular weight $\mu=1.3$. The initial left and right states are listed in Table 1. The left and right states are separated by origin. We solve the EoS by interpolation method described in Section 2.4. The results of both exact Riemann solvers are shown in Figure 6. Red circles and blue dots indicate the solution from the exact general EoS Riemann solver and the solution from the exact perfect gas Riemann solver, respectively. The solutions from the two exact Riemann solvers are very close, confirming that the exact general EoS Riemann solver can be reduced to an exact perfect gas Riemann solver when a perfect gas EoS is used. Figure 6 also shows that low $\gamma$ perfect gas has much higher compressibility compared to the high $\gamma$ perfect gas.

5.2 The correspondence between the HLLC general EoS Riemann solver and the original HLLC Riemann solver

For the same reason, as stated in Section 5.1. The HLLC general EoS Riemann solver should reduce to the original HLLC Riemann solver when the EoS tables of a perfect gas are used. The EoS is solved with the interpolation method described in Section 4.3. We carry out two tests that have the solution structure of a shock and rarefaction waves with $\gamma=1.05$ and $1.667$. The left and right states are listed in Table 2. The coordinate that separates the initial left and right states is $x=0.5$ km, and the simulation domain is $[0,1]$ km. We use the first order Godunov scheme [12] and let $\text{CFL}=0.7$. There are 100 cells in our simulation and Figure 7. From Figure 7, we can tell that the solutions from the two HLLC Riemann solvers also produce very similar results, confirming that the HLLC general EoS Riemann solver reduces to the original HLLC Riemann solver for a perfect gas.

5.3 Comparison of the exact solution to solutions of the HLLC general EoS Riemann solver with interpolation EoS solver and with analytic EoS solver

In this section, we use pure hydrogen gas with two binding energies (\chH2 disassociation and \chH ionization) as an example of realistic gas. Therefore, four species are considered (\chH2, \chH, \chH^+, and \che^-, see Appendix C for more information). We compare three kinds of numerical solvers (with different EoS schemes) on four tests using our pure hydrogen gas EoS: the exact solution calculated by the algorithm in Section 2, the HLLC general EoS Riemann solver with an interpolated EoS (discussed in Section 4.3), and the HLLC general EoS Riemann solver with an analytic EoS solver. The analytic EoS solver solves the EoS by calculating thermodynamic variables at run-time and is explicitly described in Appendix C. We show that the solution of the HLLC general EoS Riemann solver with the two EoS solvers approaches the solution of the exact solution as the $L^{1}$ norm of their difference decreases monotonically with the increasing resolution. Explicitly, we calculate,

[TABLE]

where $f_{N}$ is the solution of any specific field, $N$ is the number of points. $f_{\rm exact}$ is the solution of exact general EoS Riemann solver of any specific field, and $L=1$ km is the length of the simulation domain. In this research, we choose the number of cells $N=100,200,400,800$ (we shall use the term “resolution” interchangeably) and the density field to calculate the corresponding $L^{1}$ norm. We designed four tests to demonstrate that the HLLC general EoS Riemann solver can handle the 1D hydrodynamic problem with multiple phase transitions well, which are listed in Table 3. The coordinate that separates the initial left and right states for all these four tests is $x=0.5$ km and the simulation domain is $[0,1]$ km. We use the first order Godunov scheme and let $\text{CFL}=0.7$.

Test 2 is a strong shock tube test (same as the initial condition in Table 2). Its left state consists of fully ionized hydrogen, and its right state consists of purely molecular hydrogen. Test 3 is an asymmetric two-shock-waves test. Its left and right gases are all in a purely molecular state, but the shocked gas would be (partially) ionized. Test 4 is a symmetric two-rarefaction-wave test. Its left and right states are all fully ionized. The rarefied gas in the middle region will recombine and form atomic hydrogen. Test 5 is also an asymmetric two-shock-wave test. Its difference from test 3 is that the shocked gas will be more ionized.

When using Riemann solvers with Godunov scheme, we can get the solution of $\rho$, $v$, and $\epsilon$ after each step. The tabulated EoS solver calculates $p,\ T,\ \Gamma$, and $\mu$ with the bi-linear interpolation method discussed in Section 4.3. With $\mu$ and $\rho$, $n_{\ch{H2}},\ n_{\ch{H},\ n_{\ch{H+}}}$, and $n_{\ch{e-}}$ can be calculated analytically. We also calculate

[TABLE]

as the concentration percentage of each hydrogen species. The analytic EoS solver first calculates $T$ by solving the inverse problem of $T=T(\rho,\epsilon)$, then it calculate $p,\ \Gamma$ and all species number density $n_{i}$ by solving the forward problem of $z=z(\rho,T)$ where $z$ should be replaced by the interested variable. $\mu$ is calculated by $\mu=\rho/(m_{\mu}\sum_{i}n_{i})$ where $m_{\mu}$ is the atomic mass unit.

The solutions produced by the three aforementioned Riemann solvers with $N=100$ for these four tests are shown in Figure 8, 9, 10, and 11. The arrangement of the panels in all four figures are the same. The solutions of the tabulated EoS solver and the analytic EoS solver have little difference.

In Figure 8, the shocked gas is compressed, and the solution structure resembles a low $\gamma$ perfect gas. Since the thickness of the compressed gas is small, at $N=100$, there are only two cells inside the compressed region. The solution of the HLLC solvers gives much “diffused” density profiles. Nonetheless, the position of the shock is well captured, and the solution of $u,\ p,\ T$, and $e_{\rm g}$ are similar to the exact solution. $\Gamma$ is smooth but different from the exact solution. We think the $\Gamma(\log\rho,\log p)$ table is well resolved as $d\log p=0.011$ which correspond to $2.56\%$ change in pressure. This point is also pronounced in the analytic EoS solver as the convergence criterion requires that the change of $T$ should be less than $0.001\%$, or $dT<1$ K at $T=10000$ K. Therefore, the EoS solvers are less likely to be the cause of this defect. There could be some fundamental limitation in the current HLLC algorithm or improvement of the HLLC Riemann solver may be made in the future to improve the accuracy of the solution. $\chi_{\ch{H2}},\ \chi_{\ch{H}}$, and $\chi_{\ch{H+}}$ are related to the calculation of $\mu$ in the tabulated EoS solver or $T$ in the analytic EoS solver. Both approximate solutions roughly caption the position of the ionization front. The sensitive relation of number density to $\mu$ or $T$ and the diffusive nature of the Godunov scheme lead to the difference between the approximate solutions and the exact solution.

In Figure 9, the position of the middle wave in the approximate solution of $\rho$ and $T$ does not match the exact solution very well. There is a little bump in the solution of $e_{\rm g}$ while $p$ does not have. This situation cannot happen in perfect gas since $e_{\rm g}$ and $p$ has a linear relation in perfect gas. The bump in $e_{\rm g}$ could be related to some defect in the HLLC solver. $\Gamma$ does not give accurate solution and the reason could be the same as test 2. The solutions of $\chi_{\ch{H2}},\ \chi_{\ch{H}}$, and $\chi_{\ch{H+}}$ are more satisfactory.

What we see in Figure 10 represents the defect of Godunov scheme. The “heating” at the center of the symmetric rarefaction waves and the “dip” in density are typically found in Godunov scheme (see e.g. [18] for the magneto-hydrodynamic case and Section 6.4.2, Chapter 6 of [7]).

In Figure 11, $\rho,\ u,\ p,\ T$, and $e_{\rm g}$ are all very satisfactory. Differences can be found in the solutions of $\chi_{\ch{H}}$ between the exact and approximate solutions. This is because the gas is on the edge of becoming fully ionized therefore $\chi_{\ch{H}}$ is very sensitive to $\rho$ and $T$.

An interesting thing to note is that both test 3 and test 5 (Figures 9 and 11 respectively) consist of two strong shocks. The “accuracy” of approximate Riemann solvers seem very different. This implies that the HLLC general EoS Riemann solvers discussed in this paper may not be very accurate when dealing with partially ionized or disassociated gas, or gas that is right in phase transitions. One reason may be that we are not finding an optimal $p_{*}$ with the adaptive non-iterative Riemann solver. However, using the exact general EoS Riemann solver to calculate $p_{*}$ would be too expensive. Another possibility is that the change of $\epsilon$ during transition is too stiff. Alternatively, this error could be cause by the fact that the discontinuities spend most of their time in between cell interfaces. When this occurs the cells containing the discontinuity will contain both neutral and ionized gas and be forced to assume some intermediate phase introducing errors. Nonetheless, the HLLC general EoS Riemann solvers proposed in this paper are still available options in general EoS hydrodynamics, and these test are designed to probe the limits of our code.

We also study the efficiency of the HLLC general EoS Riemann solver with the analytic EoS solver and tabulated EoS solver. We compare them to the original HLLC Riemann solver at the same resolution of $N=800$. We calculate the solution of the nine quantities shown in Figure 8-11 with the general EoS Riemann solvers and only calculate the solution of $\rho,\ u,\ p,\ T$, and $e_{\rm g}$ with the original HLLC Riemann solver. The CPU time for the four tests in Table 3 are listed in Table 4.

The HLLC general EoS Riemann solver with analytic EoS solver and tabulated EoS solver are about 13-58 times slower and 35-132 times slower than the original HLLC Riemann solver, respectively. Since test 2 has one rarefaction wave and test 4 has two rarefaction waves, Table 4 tells us that the HLLC general EoS Riemann solver is less efficient in calculating rarefaction waves. We further conclude some advantages and disadvantages of each general EoS Riemann solvers below:

If the span of the EoS domain is big in $(\log\rho,\log p)$, the tabulated EoS Riemann solver becomes less accurate or slower. This is because the table needs to be enlarged to maintain the same resolution (or accuracy) and can slower down the interpolation process significantly which is mainly due to the increase of the memory accessing time. 2. 2.

If the gas consists of \chH2, \chH, \chH+, \chHe, \chHe+, \chHe^2+, and \che-, the analytic EoS solver needs to solve two equilibrium system and each one of them may contain a quartic equation. This can slow down the analytic EoS solver significantly. 3. 3.

If the EoS is hard to be approximated by the analytic EoS solver, for example, the rotational and vibrational degrees of freedom in \chH2, and the problem is sensitive to these two degrees of freedom, the tabulated EoS solver becomes only option. 4. 4.

When using the tabulated EoS solver, one has to estimated the range of $\log\rho$ and $\log p$ and choose the suitable EoS tables before solving the problem. An over-estimation could reduce the efficiency and an under-estimation may lead to unexpected outcome. The analytic EoS solver does not have this limitation.

With the above conclusions, we propose that a combination of the analytic EoS and tabulated EoS solvers may be more favorable. We have not implemented other high order interpolation methods in the tabulated EoS solver, but we think it very much worth investigating as long as the interpolation method preserves the monotonicity of the tables.

At the end of this section, we plot the normalized value $\log_{2}\{\delta_{N}(\rho)/\delta_{100}(\rho)\}$ of each test for both HLLC general EoS Riemann solver with tabulated and analytic EoS solvers in Figure 12. The black line represents the linear convergence rate, which is the highest convergence rate one can achieve because the exact solutions have discontinuities (except for the two rarefaction waves). Both solvers show a clear trend of decreasing $\log_{2}\{\delta_{N}(\rho)/\delta_{100}(\rho)\}$ with increasing resolution $N$ but their convergence rate is less than the linear convergence rate. The worst convergence rate happens in test 3 that has two strong shocks with partially ionized middle states.

5.4 A 2D simulation

Figure 13 shows a shock test in 2D. We use the HLLC general EoS Riemann solver with analytic EoS solver. The simulation domain is $[-1,1]\times[-1,1]$ km2 with a fixed resolution of $200\times 200$. We place a circle with $\rho=10^{-13}$ g$\cdot$cm*-3and $T=20000$ K at the origin and its radius is $0.2$ km at $t=0$ s. The ambient has an initial condition of $\rho=10^{-15}$ g$\cdot$cm-3*and $T=300$ K. In this 2D simulation, CFL$=0.7$ and we confirm that $(v_{x}/dx+v_{y}/dy)\Delta t<1$ at all time in the domain. $v_{x}$ and $v_{y}$ in the x and y direction. $dx=dy=0.01$ km is the grid size and $\Delta t$ is the time step. The 2D integration method used in this simulation is symmetric Strang splitting [19].

Each row in Figure 13 shows $\rho,\ \chi_{\ch{H2}},\ \chi_{\ch{H},}$ and, $\chi_{\ch{H+}}$ at different time. $\chi=0$ is represented by the white area as we divide the EoS into six regions and set the number density of some species to be [math] in some region. This simulation shows an outward shock and a inward rarefaction waves. The center of the “fireball” cools down and $\ch{H}$ appears. The shocked gas is highly compressed due to the geometric effect and the higher compressibility of realistic gas.

6 Conclusions and discussions

Following the same principle of the original HLLC Riemann solver, we devise an HLLC general EoS Riemann solver that can solve realistic gas dynamics with drastically varying $C_{V}$ and $C_{P}$. We intentionally keep the modification of the original HLLC Riemann solver to a low extent to achieve a fast speed.

To examine the Riemann problem for a realistic fluid and determine the accuracy of the HLLC general EoS Riemann solver, we implemented an exact general EoS Riemann solver based on [11]. We improved the robustness of the root-finding algorithm. The existence and uniqueness of the solution are guaranteed if the EoS tables and the interpolated ones satisfy the two monotone condition in Equation 9-10. The new exact general EoS Riemann solver can be used even when a convex condition $\left(\frac{\partial^{2}p}{\partial\rho^{2}}\right)_{s}>0$ is not satisfied. We applied the exact general EoS Riemann solver to the Godunov scheme in Section 3. To examine the correctness of our new algorithm, we compared our exact general EoS Riemann solver to the exact ideal gas Riemann solver for shock tube tests of perfect gases with $\gamma=\Gamma=1.05$ and $\gamma=\Gamma=1.667$. The solutions are comparable, confirming that our exact general EoS Riemann reduces to the original exact ideal gas Riemann solver when using the corresponding perfect gas EoS. We also compared our HLLC general EoS Riemann solver to the original HLLC Riemann solver by running tests with perfect gas EoS, resulting in comparable solutions. Therefore, our HLLC general EoS Riemann solver reduces to the original HLLC Riemann solver in the case of a perfect gas.

In Section 5.3, we compare the solution of the HLLC general EoS Riemann solver with two different EoS solvers to the exact solution calculated by the algorithm in Section 2. One EoS solver is called the tabulated EoS solver; it calculates the thermodynamic variables by interpolating 2D tables with a bi-linear interpolation method (Section 4.3). Another EoS solver is called the analytic EoS solver; it calculates the thermodynamic variables by solving EoS at run-time (see Appendix C for details).

We find reasonably consistent results between the approximate and exact Riemann solvers for the shock-tube and two-shock tests. The HLLC solver can resolve the contact discontinuity, which is different from the HLL solver [20]. The location and the strength of the shock can be well captured. Although discrepancy is found in the two-rarefaction-wave test, we think defect is rooted in the Godunov scheme. By comparing two tests (test 3 and 5) that have two strong shock waves with different ionization level, we infer that the HLLC general EoS Riemann solver still has difficulty in resolving shocked gases that is in the middle of a phase transition. The efficiency of the new HLLC general EoS Riemann solver with both EoS solvers are compared to the original HLLC Riemann solver. We find that the analytic EoS solver is roughly 20 times slower than the original HLLC solver, while the tabulated EoS solver is roughly 60 times slower than the original HLLC solver. However, we point out that these multiples are not constant depending on the problem being solved in Section 5.3. In short, tabulated EoS solver becomes slower in problems with a large range of $\log\rho$ or $\log p$, and analytic EoS solver becomes slower in problems with many species in equilibrium. A combination of the two EoS solvers may become favorable in problems that have a large range of $\log\rho$, $\log p$, and many species present at the same time. In multi-physics simulations that include self-gravity (Poisson equation), complicated radiation transfer [21, 22], etc., hydro-step only takes a fraction of the total computational time. If the temperature is critical in those simulations, the HLLC solvers proposed in this research are still useful even though they are much more expensive than the perfect gas HLLC Riemann solver. Convergence study has been performed by comparing the solutions of our HLLC general EoS Riemann solver to the exact solution. We use the $L^{1}$ norm of density to demonstrate that the HLLC Riemann solver approaches the exact solution with increasing resolution. We found that the convergence rates are below the linear convergence rate, with the worst one happens in the case of two shocks with partially ionized gas.

Acknowledgments

We thank James Stone, Kengo Tomida, Yan-Fei Jiang, Jonathan Carroll-Nellenback, Natasha Ivanova, and Xuening Bai for their useful discussions and insight generated from their work. We also thank the anonymous referee who’s feedback lead to significant improvements of this paper. ZC is grateful to the Horton Fellowship at the University of Rochester and the CITA National Postdoctoral Fellowship for financial support. MC gratefully acknowledges support from the Institute for Advanced Study, NSF via grant AST-1515763, and NASA via grant 14-ATP14-0059. AF appreciates financial support from the US Department of Energy through grant DE-SC0001063 ,the NSF through grant AST-1515648, and the Space Telescope Science Institute through grant HST-AR-12832.01-A.

Appendix A Basic thermodynamics: differentiating $\gamma$ and $\Gamma$

We use $\gamma$ to represent the ratio of specific heats, i.e. $\gamma=C_{p}/C_{V}$ where $C_{p}$ and $C_{V}$ are the constant pressure heat capacity and the constant volume heat capacity, respectively. We use $\Gamma$ to represent the adiabatic index which is defined as

[TABLE]

and the adiabatic sound speed

[TABLE]

where the subscript $s$ means constant specific entropy.

In this section, we will show that $\gamma=\Gamma$ for perfect gas but not necessarily for realistic gas.

First, we shall derive $\gamma$. A general thermal dynamic relation that consider multiple species is:

[TABLE]

where $i$ stands for different species. U,\ H,\ s,\ \mu_{i},\and $N_{i}$ ($\mu_{i}$ is not to be confused with the mean atomic weight $\mu$) are the internal energy, enthalpy, entropy, chemical potential, and particle number, respectively. By the definition of $C_{V}$ and $C_{p}$,

[TABLE]

therefore,

[TABLE]

If we further consider ideal gas with fixed species number, i.e. $dN_{i}=0$, then $pV=Nk_{b}T$ where $N$ is the total particle number. The key step in deriving $\Gamma$ is to write $dT$ in terms of $dp$ and $d\rho$ in adiabatic process,

[TABLE]

where we have used $V=1/\rho$. Substitute Equation 66 into Equation 62 and rearrange the equation,

[TABLE]

therefore,

[TABLE]

Compare Equation 64 to Equation 68, we can conclude that, for constant species number ideal gas (and perfect gas),

[TABLE]

therefore $\gamma=\Gamma$, which is a coincidence. For realistic gas, $dN_{i}\neq 0$ thus even with ideal gas assumption,

[TABLE]

becomes nontrivial and it is even more difficult to express $dN_{i}$ in terms of $dp$ and $d\rho$ since $dN_{i}$ depends on $(\rho,T)$ in general and is solved with Saha equation. The strategy we used to derive $\Gamma$ fails. In fact, $\gamma=\Gamma$ may not hold anywhere in $(\rho,T)$ for realistic gas. We will show the calculation of $\Gamma$ for a simplified pure hydrogen gas in Appendix C.

Appendix B Quasilinear Approximation

The computation of the characteristics (i.e. speeds in the Riemann fan) is typically done with the quasilinear approximation.

B.1 Conservative Form

First we start with the Euler equations in differential, conservative form:

[TABLE]

These can be expressed as

[TABLE]

where $\mathbf{U}$, $\mathbf{F}$ are the conservative variables and their associated fluxes respectively:

[TABLE]

[TABLE]

where $E$ is the total energy density, $H$ is the total specific enthalpy

[TABLE]

and $U_{i}$ are the components of $\mathbf{U}$. We can then linearize equation 72 into

[TABLE]

where $\mathbf{A}$ is the Jacobian of $\mathbf{F}$

[TABLE]

where

[TABLE]

and

[TABLE]

Side note: Our present analysis should extend to the case of a more general EoS of the form $p(\rho,\epsilon,\bm{\psi})$ as long as $p_{\rho}>0$, $p_{\epsilon}>0$, and Eqn. 79 holds, implying

[TABLE]

For example $\bm{\psi}$ could be some collection of passively advected scalars.

We can write $\mathbf{A}$ in a more manageable form by utilizing the relation below

[TABLE]

to remove all instances of $p_{\rho}$. Where $a$ is the adiabatic sound speed.

[TABLE]

The eigenvalues along with the associated left and right eigenvectors of $\mathbf{A}$ are,

[TABLE]

Surprisingly, the only difference in the right eigenvectors $\left(\mathbf{r^{(i)}}\right)$ compared to the case of a perfect gas EoS (i.e. $\gamma$-law) is $\mathbf{r}_{3}^{(2)}$, in which $H-a^{2}\rho/p_{e}$ simplifies to $u^{2}/2$.

B.2 Requirements on the EoS to keep the Euler equations hyperbolic

It is crucial to make sure that the linearized Euler equations are hyperbolic with the EoS since Riemann solvers are based on Riemann invariants and Riemann invariants exist only when the PDEs are hyperbolic. In section B.1, we have listed the eigenvalues of the Jacobian of the quasilinear system. The quasilinear system is hyperbolic if and only if all the eigenvalues are real and the Jacobian is diagonalizable. It is apparent that when $a^{2}>0$, the three eigenvalues are distinct and the quasilinear system is hyperbolic.

In equation 82, $\epsilon_{p}>0$ because the specific internal energy is usually positively related to the pressure when the density is held constant. Therefore, the requirements of the EoS to keep the Euler equations hyperbolic now reduce to

[TABLE]

In this research, the EoS that we use has $\epsilon_{\rho}<0$, thus the linearized Euler equations are hyperbolic.

B.3 Jump condition for the contact discontinuity

For quasilinear hyperbolic systems, the Riemann invariants admitted by 77 follow the following ODEs [23, 7]

[TABLE]

The middle characteristic ($k=2$) is the contact discontinuity, and is the only characteristic for which we do not need to worry about non-linear shocks. From $i=1$ and $j=2$ we get

[TABLE]

i.e. constant velocity across the contact discontinuity. Using the above in combination with $i=1$ and $j=3$ we retrieve

[TABLE]

By recalling Eqn. 82 we can further simplify this to

[TABLE]

i.e. pressure is constant across the contact discontinuity.

Appendix C Idealized Hydrogen EoS

To test our Riemann solvers, we consider a simple, but useful, EoS which contains a disassociation transition and an ionization transition. Our simple EoS describes a plasma with only four species: molecular hydrogen (\chH2), neutral hydrogen (\chH) protons/ionized hydrogen (\chH+), and electrons (\che-). Note that we do not consider the vibrational and rotational energy in \chH2 [24, 25] here because: (1) the partition function of \chH2 will be very complicated if we take these two energies into consideration and will kill a fast algorithm we will introduce hereafter, (2) these two kinds of energies are not significant compared to the binding energies, (3) the EoS is still monotone thus all the conclusions regarding the general EoS Riemann solvers hold, and (4) the ratio of orthohydrogen to parahydrogen may still be an assumption. However, we admit that it is desirable to include the vibrational and rotational energies, as this is necessary to recover the correct variation of $\Gamma$ below the molecular dissociation temperature. This added complexity is computationally prohibitive, requiring the interpolation of EoS. The EoS study here is crude in the sense that we only require it to manifest the monotonicity. We also want to present a fast algorithm in the case if we only consider the disassociation and the ionization binding energies. A more thorough study of the error of ignoring the vibrational and the rotational energies should be done in the future, but due to the page limit, we cannot present it here. We shall assume local thermal equilibrium (LTE), ideal gas law ($p=nkT$) and consider only two reactions in equilibrium

[TABLE]

To derive this EoS we need to know the relevant partition functions. $Z_{i}$ is the partition function per volume for species $i$ (henceforth partition function will be used to mean partition function per volume). The partition function can be broken down into parts

[TABLE]

These parts are: internal bound states $Z_{i}^{\rm bound}$, nuclear spin $Z_{i}^{\rm nuc}$, electron spin and angular momentum $Z_{i}^{\rm elec}$, translation $Z_{i}^{\rm tr}$, and excitation $Z_{i}^{\rm exct}$. To further simplify the problem we shall assume that

[TABLE]

While this is not technically true for \chH, it is a relatively small effect on the EoS, but the lack of rotational and vibrational states for $\ch{H}_{2}$ effects $\Gamma$. However, this EoS is meant as a simple proof of concept and not for high precision applications.

We now list the partition functions for all the species. The partition functions not explicitly described are unity. The translational partition function has the same form for all species

[TABLE]

where $m_{i}$ is the mass of the $i$th species and $h$ is the Planck constant.

The remaining non-trivial partition functions are

[TABLE]

where $\phi_{\rm dis}=7.17\times 10^{-12}$erg [26] and $\phi_{\rm ion}=2.18\times 10^{-11}$erg while all the unspecified partition functions are unity.

To compute the EoS we must first solve number conservation, Saha equations corresponding to Equations 90 and 91, and charge neutrality to get the number density of all species:

[TABLE]

respectively, where $n_{i}$ is the number density of the $i$th species. Substitute Equation 104, 105, and 106 into 103, we can obtain a quartic equation of $n_{\ch{H+}}$,

[TABLE]

Once $n_{\ch{H+}}$ is determined, the number density of other species can be obtained by solving Equation 104, 105, and 106. Fourth degree polynomials are the highest degree where an explicit solution is guarantied. However, solving Equation 107 in computer with finite accuracy will not always give satisfactory result when $q_{\rm dis}$ and $q_{\rm ion}$ are too large or too small. Besides, solving a quartic equation is much more expensive than solving a quadratic equation. To improve the computational efficiency, we need to minimize the solution of the quartic equation. Therefore, we divide the EoS in density and temperature space $(\rho,T)$ into six regions with the following division criterion: if the concentration of \chH (Equation 55-57) in any species is less than $10^{-6}$, we assume such species does not present. The division result is shown in Figure 14.

In the meantime, we get four curves in $(\rho,T)$ that delineate these six regions. They are:

$T_{1}(\rho)$ separates the \chH2 only region and the multi-species region.

• 2.

$T_{2}(\rho)$ separates the \chH2 present and non-present region.

• 3.

$T_{3}(\rho)$ separates the \chH+ present and non-present region.

• 4.

$T_{4}(\rho)$ separates the \chH+ and \che- only region and the multi-species region.

Note that $T_{2}(\rho)<T_{3}(\rho)$ in region III but $T_{2}(\rho)>T_{3}(\rho)$ in region IV. These four $(\rho,T)$ curves will be used in calculating the temperature with $(\rho,p)$ and $(\rho,\epsilon)$.

The blue color indicates the regions with single species. When $(\rho,T)$ falls into these regions, we can easily calculate the number density of that species.

If $(\rho,T)$ belongs to II, we assume $n_{\ch{H+}}=n_{\ch{e-}}=0$. We solve Equation 103 and 104. Substitute Equation 104 into 103, we can get a quadratic equation of $n_{\ch{H}}$,

[TABLE]

and the physical solution of Equation 108 is,

[TABLE]

If $(\rho,T)$ belongs to IV, we assume all four species present and we have to solve Equation 107 and then Equation 104, 105, and 106.

If $(\rho,T)$ belongs to V, we assume $n_{\ch{H2}}=0$. We solve Equation 103, 105, and 106. Substitute Equation 105 and 106 into 103, we can get a quadratic equation of $n_{\ch{H+}}$,

[TABLE]

and the physical solution of Equation 110 is,

[TABLE]

After solving $n_{i}$, the pressure and specific internal energy can be calculated by,

[TABLE]

[TABLE]

where

[TABLE]

are the internal energy per particle of each species.

We also need to calculate the adiabatic sound speed, which is needed in the HLLC general EoS Riemann solver. If the specific entropy $s$ is differentiable, the adiabatic sound speed $a$ can be calculated by,

[TABLE]

where

[TABLE]

Although there are many derivatives in the expression of $S_{T},S_{\rho},p_{T}$, and $p_{\rho}$, the derivatives of the partition functions are simply,

[TABLE]

This would not be the case if one consider the vibrational and rotational degree of freedom of \chH2 as in [24, 25].

The only derivatives that we need to be calculated are $\left(\frac{\partial\ln n_{i}}{\partial\ln T}\right)_{\rho}$ and $\left(\frac{\partial\ln n_{i}}{\partial\ln\rho}\right)_{T}$. When all species present, i.e., $(\rho,T)\in$IV, let us re-write Equation 103-105 in log form while using 106,

[TABLE]

Let the right-hand side of Equation 132 be $\mathbf{L}$. Then,

[TABLE]

and,

[TABLE]

Then,

[TABLE]

and

[TABLE]

When $(\rho,T)\in$I

[TABLE]

and all the other derivatives are set to be [math]. When $(\rho,T)\in$III,

[TABLE]

and all the other derivatives are set to be [math]. When $(\rho,T)\in$V,

[TABLE]

and all the other derivatives are set to be [math]. When $(\rho,T)\in$II or V, the derivatives of non-presenting species are set to be [math] while the derivatives of existing species can be calculated by solving a rank 2 linear system following the same step as four species.

It is also beneficial and necessary to be able to calculate temperature for a given $(\rho,p)$ or $(\rho,\epsilon)$. We design a fast algorithm to find such temperature with the EoS that satisfy:

[TABLE]

In fact, almost all ideal gas EoS satisfy these monotone conditions. We now prove that our four monotone conditions are sufficient to the Equation 9 and 10.

[TABLE]

using Equation 147, 148, and

[TABLE]

Note that when the inequality 149 becomes an equality, Equation 152 become undefined. This is because $\epsilon(\rho,T)$ is reduced to $\epsilon(T)$. For simplicity, we define this situation to be [math]. The EoS of perfect gas is a such example. On the other hand,

[TABLE]

because of Equation 148 and 150.

The steps for solving $T=T(\rho,p)$ are,

1. Step 1

With the given $\rho$, find the four temperature $T_{1}(\rho)$, $T_{2}(\rho)$, $T_{3}(\rho)$, and $T_{4}(\rho)$ by linear interpolation. Calculate the corresponding $p_{i}$ by solving the species number and ideal gas law $p_{i}=n(T_{i})kT_{i}$. 2. Step 2

If $p\leq p_{1}$, the gas is in purely molecular state and

[TABLE]

If $p_{2}\leq p\leq p_{3}$, the gas is in fully atomic state and

[TABLE]

If $p\geq p_{4}$, the gas is in fully ionized state and

[TABLE] 3. Step 3

Otherwise, find the lower ($p_{\rm l}$) and upper ($p_{\rm u}$) limit of the pressure of the region that $(\rho,p)$ belongs to. $p_{\rm l}<p<p_{\rm u}$, use bisection method to determine $T\in[T_{l},T_{u}]$ and dynamically update the lower and upper bound ($T_{\rm l}$ and $T_{\rm u}$). In this research, we define the solution as $T=(T_{\rm l}+T_{\rm u})/2$ when $2(T_{\rm u}-T_{\rm l})/(T_{\rm u}+T_{\rm l})<10^{-5}$.

The steps for solving $T=T(\rho,\epsilon)$ follow the same logic as solving $T=T(\rho,p)$.

We present a set of example EoS tables in Figure 15, 16, 17, 18, and 19. We show $a,\ \Gamma,\ \mu,\ \epsilon$, and $p$ in temperature range $T\in[10^{1},10^{6}]$ with $d\log T=0.005$ and $\rho\in[10^{-20},10^{0}]$ with $d\log\rho=0.02$, where $T$ and $\rho$ are given in units of K and g$\cdot$cm*-3*, respectively.