Physical-constraints-preserving Lagrangian finite volume schemes for one- and two-dimensional special relativistic hydrodynamics

TL;DR

This paper develops and analyzes physical-constraints-preserving Lagrangian finite volume schemes for 1D and 2D special relativistic hydrodynamics, ensuring positivity of density and pressure, and bounded velocities, with high accuracy and effectiveness demonstrated through numerical experiments.

Contribution

It introduces PCP Lagrangian schemes with proven positivity and boundedness properties for special relativistic hydrodynamics, including high-order accuracy techniques.

Findings

Proved PCP property for first-order schemes with HLLC solver.

Developed high-order schemes using SSP time discretization and WENO reconstruction.

Numerical experiments confirm accuracy and robustness in complex RHD scenarios.

Abstract

This paper studies the physical-constraints-preserving (PCP) Lagrangian finite volume schemes for one- and two-dimensional special relativistic hydrodynamic (RHD) equations. First, the PCP property (i.e. preserving the positivity of the rest-mass density and the pressure and the bound of the velocity) is proved for the first-order accurate Lagrangian scheme with the HLLC Riemann solver and forward Euler time discretization. The key is that the intermediate states in the HLLC Riemann solver are shown to be admissible or PCP when the HLLC wave speeds are estimated suitably. Then, the higher-order accurate schemes are proposed by using the high-order accurate strong stability preserving (SSP) time discretizations and the scaling PCP limiter as well as the WENO reconstruction. Finally, several one- and two-dimensional numerical experiments are conducted to demonstrate the accuracy and the effectiveness of the PCP Lagrangian schemes in solving the special RHD problems involving strong discontinuities, or large Lorentz factor, or low rest-mass density or low pressure etc.

Tables

Table 1. Table 4.1: Example 4.1 : Errors and orders of convergence t = 0.02 𝑡 0.02 t=0.02 obtained by using the first-order scheme.

$N$	${\epsilon }_1$	Order	${\epsilon }_2$	Order	${\epsilon }_{\infty }$	Order
20	1.687e-02	-	1.729e-02	-	3.350e-02	-
40	8.356e-03	1.01	8.523e-03	1.02	1.698e-02	0.98
80	4.260e-03	0.97	4.288e-03	0.99	8.418e-03	1.01
160	2.171e-03	0.97	2.169e-03	0.98	4.272e-03	0.98
320	1.098e-03	0.98	1.092e-03	0.99	2.139e-03	1.00

Table 2. Table 4.2: Same as Table 4.1 , except for the third-order scheme.

$N$	${\epsilon }_1$	Order	${\epsilon }_2$	Order	${\epsilon }_{\infty }$	Order	${\Theta }_N$
20	9.067e-02	-	1.038e-01	-	2.317e-01	-	25.0%
40	1.785e-02	2.34	2.289e-02	2.18	6.573e-02	1.82	20.8%
80	3.373e-03	2.40	4.356e-03	2.39	1.368e-02	2.26	12.8%
160	3.776e-04	3.16	4.932e-04	3.14	1.695e-03	3.01	6.49%
320	3.306e-05	3.51	4.151e-05	3.57	1.715e-04	3.31	3.12%

Table 3. Table 4.3: Example 4.3 : Errors and orders of convergence at t = 1 𝑡 1 t=1 .

$N$	${\epsilon }_1$	Order	${\epsilon }_2$	Order	${\epsilon }_{\infty }$	Order
20	1.287e-01	-	1.848e-01	-	9.417e-01	-
40	6.317e-02	1.03	8.755e-02	1.08	4.139e-01	1.19
80	3.148e-02	1.00	4.282e-02	1.03	2.056e-01	1.01
160	1.587e-02	0.99	2.139e-02	1.00	1.009e-01	1.03
320	7.946e-03	1.00	1.065e-02	1.01	4.918e-02	1.04

Table 4. Table 4.4: Same as Table 4.3 except for the second-order scheme.

$N$	${\epsilon }_1$	Order	${\epsilon }_2$	Order	${\epsilon }_{\infty }$	Order	${\Theta }_N$
20	8.131e-02	-	1.303e-01	-	6.264e-01	-	1.52%
40	2.199e-02	1.89	3.769e-02	1.79	2.008e-01	1.64	0.38%
80	5.458e-03	2.01	1.008e-02	1.90	7.842e-02	1.36	0.19%
160	1.277e-03	2.10	2.393e-03	2.07	2.037e-02	1.94	0.082%
320	3.000e-04	2.09	5.559e-04	2.11	4.390e-03	2.21	0.015%

Equations

$$\frac{\partial\mbox{\boldmath\small$U$}}{\partial t}+\sum\limits_{\ell=1}^{d}\frac{\partial\mbox{\boldmath\small$F$}_{\ell}(\mbox{\boldmath\small$U$})}{\partial x_{\ell}}=0,$$

$$\mbox{\boldmath\small$U$}=(D,\mbox{\boldmath\small$m$},E)^{T},\quad\mbox{\boldmath\small$F$}_{\ell}=\left(Du_{\ell},\mbox{\boldmath\small$m$}u_{\ell}+p\mbox{\boldmath\small$e$}_{\ell},(E+p)u_{\ell}\right)^{T},$$

$$\frac{d}{dt}\int_{\varOmega(t)}\mbox{\boldmath\small$U$}dV+\int_{\partial\varOmega(t)}\mbox{\boldmath\small$\mathcal{F}$}_{\mbox{\boldmath\small$n$}}(\mbox{\boldmath\small$U$})ds=\mbox{\boldmath\small$0$},~{}\mbox{\boldmath\small$\mathcal{F}$}_{\mbox{\boldmath\small$n$}}(\mbox{\boldmath\small$U$})=\left(0,p\mbox{\boldmath\small$n$}^{T},p\mbox{\boldmath\small$u$}\cdot\mbox{\boldmath\small$n$}^{T}\right)^{T},$$

$$D=\rho\gamma,\quad\mbox{\boldmath\small$m$}=Dh\gamma\mbox{\boldmath\small$u$},\quad E=Dh\gamma-p.$$

$$p=(\Gamma -1)\rho e,$$

$$c_s=\sqrt{\Gamma p/(\rho h)},$$

$$c_s^2=\frac{\Gamma p}{\rho h}=\frac{\Gamma p}{\rho +\frac{p}{\Gamma -1}+p}=\frac{(\Gamma -1)\Gamma p}{(\Gamma -1)\rho +\Gamma p}<\Gamma -1<1.$$

$$E+p=D\gamma +\frac{\Gamma }{\Gamma -1}p{\gamma }^2,$$

$$\rho=\frac{D}{\gamma},\ h=1+\frac{\Gamma p}{(\Gamma-1)\rho},\ \mbox{\boldmath\small$u$}=\frac{\mbox{\boldmath\small$m$}}{Dh}.$$

$$x_{i+\frac{1}{2}}^{n+1}=x_{i+\frac{1}{2}}^n+\Delta t^n{u}_{i+\frac{1}{2}}^n,i=1,\cdots ,N.$$

$$\overline{\mbox{\boldmath\small$U$}}_{i}^{n+1}\Delta x_{i}^{n+1}=\overline{\mbox{\boldmath\small$U$}}_{i}^{n}\Delta x_{i}^{n}-\Delta t^{n}\left(\widehat{\mbox{\boldmath\small$\mathcal{F}$}}(\mbox{\boldmath\small$U$}_{i+\frac{1}{2}}^{-},\mbox{\boldmath\small$U$}_{i+\frac{1}{2}}^{+})-\widehat{\mbox{\boldmath\small$\mathcal{F}$}}(\mbox{\boldmath\small$U$}_{i-\frac{1}{2}}^{-},\mbox{\boldmath\small$U$}_{i-\frac{1}{2}}^{+})\right),$$

$$\mbox{\boldmath\small$U$}_{t}+\widetilde{\mbox{\boldmath\small$F$}}(\mbox{\boldmath\small$U$})_{\xi}=0,\ \ \mbox{\boldmath\small$U$}(\xi,0)=\left\{\begin{array}[]{ll}\mbox{\boldmath\small$U$}^{-},&\xi<0,\\ \mbox{\boldmath\small$U$}^{+},&\xi>0,\end{array}\right.$$

$$\mbox{\boldmath\small$U$}=(D,m,E)^{T},\quad\widetilde{\mbox{\boldmath\small$F$}}(\mbox{\boldmath\small$U$})=\left(D(u-s^{\ast}),m(u-s^{\ast})+p,E(u-s^{\ast})+pu\right)^{T}.$$

$$s^{-}=\min\big{(}s_{\min}(\mbox{\boldmath\small$U$}^{-}),s_{\min}(\mbox{\boldmath\small$U$}^{+})\big{)},~{}~{}s^{+}=\max\big{(}s_{\max}(\mbox{\boldmath\small$U$}^{-}),s_{\max}(\mbox{\boldmath\small$U$}^{+})\big{)},$$

$$s_{\min}(\mbox{\boldmath\small$U$})=\frac{u-c_{s}}{1-c_{s}u},~{}~{}~{}s_{\max}(\mbox{\boldmath\small$U$})=\frac{u+c_{s}}{1+c_{s}u}.$$

$$\widetilde{\mbox{\boldmath\small$F$}}^{\ast,-}=\widetilde{\mbox{\boldmath\small$F$}}(\mbox{\boldmath\small$U$}^{-})-(s^{-}-s^{\ast})\mbox{\boldmath\small$U$}^{-}-\frac{1}{\Delta t^{n}}\int_{\Delta t^{n}(s^{-}-s^{\ast})}^{0}\mbox{\boldmath\small$U$}(\xi,t_{n+1})d\xi.$$

$$\widetilde{\mbox{\boldmath\small$F$}}^{\ast,+}=\widetilde{\mbox{\boldmath\small$F$}}(\mbox{\boldmath\small$U$}^{+})-(s^{+}-s^{\ast})\mbox{\boldmath\small$U$}^{+}+\frac{1}{\Delta t^{n}}\int_{0}^{\Delta t^{n}(s^{+}-s^{\ast})}\mbox{\boldmath\small$U$}(\xi,t_{n+1})d\xi.$$

$$\mbox{\boldmath\small$U$}^{\ast,\pm}=\frac{1}{\Delta t^{n}(s^{\pm}-s^{\ast})}\int_{0}^{\Delta t^{n}(s^{\pm}-s^{\ast})}\mbox{\boldmath\small$U$}(\xi,t_{n+1})d\xi,$$

$$\widetilde{\mbox{\boldmath\small$F$}}^{\ast,\pm}=\widetilde{\mbox{\boldmath\small$F$}}(\mbox{\boldmath\small$U$}^{\pm})-(s^{\pm}-s^{\ast})(\mbox{\boldmath\small$U$}^{\pm}-\mbox{\boldmath\small$U$}^{\ast,\pm}),$$

$$D^{\ast ,\pm }(s^{\pm }-u^{\ast ,\pm })=D^{\pm }(s^{\pm }-u^{\pm }),$$

$$m^{\ast ,\pm }(s^{\pm }-u^{\ast ,\pm })=m^{\pm }(s^{\pm }-u^{\pm })+p^{\ast ,\pm }-p^{\pm },$$

$$E^{\ast ,\pm }(s^{\pm }-u^{\ast ,\pm })=E^{\pm }(s^{\pm }-u^{\pm })+p^{\ast ,\pm }{u}^{\ast ,\pm }-p^{\pm }{u}^{\pm }.$$

$$(A^{\pm }+s^{\pm }{p}^{\ast ,\pm })u^{\ast ,\pm }=B^{\pm }+p^{\ast ,\pm },$$

$$A^{\pm }=s^{\pm }{E}^{\pm }-m^{\pm },B^{\pm }=m^{\pm }(s^{\pm }-u^{\pm })-p^{\pm }.$$

$$p^{\ast ,-}=p^{\ast ,+}=:p^{\ast },u^{\ast ,-}=u^{\ast ,+}=:u^{\ast }(=s^{\ast }).$$

$$(1-s^{\pm }{s}^{\ast })p^{\ast }=s^{\ast }{A}^{\pm }-B^{\pm }.$$

$${\mathcal{C}}_0+{\mathcal{C}}_1{s}^{\ast }+{\mathcal{C}}_2(s^{\ast }{)}^2=0,$$

$${\mathcal{C}}_0=B^{+}-B^{-},{\mathcal{C}}_1=A^{-}+s^{+}{B}^{-}-A^{+}-s^{-}{B}^{+},{\mathcal{C}}_2=s^{-}{A}^{+}-s^{+}{A}^{-}.$$

$$s^{\ast }=u^{\ast }=\frac{-{\mathcal{C}}_1-\sqrt{{\mathcal{C}}_1^2-4{\mathcal{C}}_0{\mathcal{C}}_2}}{2{\mathcal{C}}_2},$$

$$p^{\ast }=\frac{s^{\ast }{A}^{-}-B^{-}}{1-s^{-}{s}^{\ast }}=\frac{s^{\ast }{A}^{+}-B^{+}}{1-s^{+}{s}^{\ast }}.$$

$$\mbox{\boldmath\small$U$}_{i+\frac{1}{2}}^{-}=\overline{\mbox{\boldmath\small$U$}}_{i}^{n},\quad\ \mbox{\boldmath\small$U$}_{i+\frac{1}{2}}^{+}=\overline{\mbox{\boldmath\small$U$}}_{i+1}^{n}.$$

Full text

Physical-constraints-preserving Lagrangian finite volume schemes for one- and two-dimensional special relativistic hydrodynamics

Dan Ling1, Junming Duan111HEDPS, CAPT & LMAM, School of Mathematical Sciences, Peking University, Beijing 100871, China. E-mail: dan$\_$ling@math.pku.edu.cn; duanjm@pku.edu.cn, Huazhong Tang*1,*222School of Mathematics and Computational Science, Xiangtan University, Hunan Province, Xiangtan 411105, P.R. China. E-mail: hztang@math.pku.edu.cn

Abstract

This paper studies the physical-constraints-preserving (PCP) Lagrangian finite volume schemes for one- and two-dimensional special relativistic hydrodynamic (RHD) equations. First, the PCP property (i.e. preserving the positivity of the rest-mass density and the pressure and the bound of the velocity) is proved for the first-order accurate Lagrangian scheme with the HLLC Riemann solver and forward Euler time discretization. The key is that the intermediate states in the HLLC Riemann solver are shown to be admissible or PCP when the HLLC wave speeds are estimated suitably. Then, the higher-order accurate schemes are proposed by using the high-order accurate strong stability preserving (SSP) time discretizations and the scaling PCP limiter as well as the WENO reconstruction. Finally, several one- and two-dimensional numerical experiments are conducted to demonstrate the accuracy and the effectiveness of the PCP Lagrangian schemes in solving the special RHD problems involving strong discontinuities, or large Lorentz factor, or low rest-mass density or low pressure etc.

Keywords: Lagrangian scheme; physical-constraints-preserving scheme; special relativistic hydrodynamics; high order accuracy.

1 Introduction

The paper is concerned with the study of the physical-constraints-preserving (PCP) Lagrangian schemes (which preserve the positivity of the rest-mass density and pressure, and the bound of the fluid velocity) for one- and two-dimensional special relativistic hydrodynamic (RHD) equations. The $d$ dimensional governing equations of the special RHDs is a system of first-order quasi-linear partial differential equations, see e.g. [16], and in the laboratory frame, it can be written in the divergence form

[TABLE]

with the conservative vector $U$ and the flux \mbox{\boldmath\smallF}_{\ell}, defined respectively by

[TABLE]

where $D$, \mbox{\boldmath\smallm}=(m_{1},\cdots,m_{d}), $E$ and $p$ are the mass, momentum and total energy relative to the laboratory frame and the gas pressure, respectively, \mbox{\boldmath\smalle}_{\ell} is the row vector denoting the $\ell$th row of the unit matrix of size $d$. In the Lagrangian framework, (1.1) can be expressed as the integral form

[TABLE]

where $\varOmega(t)$ is the mass volume with the boundary $\partial\varOmega(t)$ and $n$ denotes the unit outward normal of $\partial\varOmega(t)$. The conservative vector $U$ can be explicitly related to the primitive variables \mbox{\boldmath\smallV}=(\rho,\mbox{\boldmath\smallu},p)^{T} by

[TABLE]

Here $\rho$ denotes the proper rest-mass density, and the Lorentz factor \gamma=(1-|\mbox{\boldmath\smallu}|^{2})^{-1/2} and the specific enthalpy $h=1+e+\frac{p}{\rho}$ when units are normalized so that the speed of light is $c=1$. The symbol $e$ is the specific internal energy and the fluid velocity \mbox{\boldmath\smallu}=(u_{1},\cdots,u_{d}). The system (1.1) should be closed by an additional equation of state (EOS). Our discussion in this paper will be restricted to the perfect gas, whose EOS is formulated as follows

[TABLE]

with the adiabatic index $\Gamma\in(1,2]$. Such restriction on $\Gamma$ is reasonable under the compressibility assumptions and the adiabatic index is taken as 5/3 for the mildly relativistic case and 4/3 for the ultra-relativistic case. For the EOS (1.5), the sound speed $c_{s}$ is given by

[TABLE]

which satisfies the following inequality

[TABLE]

It is an important result and will be used in proving the PCP property of the HLLC Riemann solver.

The inverse of (1.4) cannot be explicitly given and involves solving the nonlinear equation, e.g. the pressure equation

[TABLE]

for the EOS (1.5), where the Lorentz factor $\gamma$ has been expressed as \gamma=\big{(}1-|\mbox{\boldmath\smallm}|^{2}/(E+p)^{2}\big{)}^{-1/2}. Once $p$ is obtained by solving the equation (1.8), the rest-mass density $\rho$, the specific enthalpy $h$, and the velocity $u$ can be orderly calculated as follows

[TABLE]

The system (1.1) takes into account the relativistic description for the dynamics of the fluid (gas) at nearly speed of light. The relativistic fluid flow appears in investigating numerous astrophysical phenomena from stellar to galactic scales, e.g. coalescing neutron stars, core collapse supernovae, formation of black holes, active galactic nuclei, superluminal jets and gamma-ray bursts, etc. Due to the relativistic effect, especially the appearance of the Lorentz factor, the nonlinearity of the system (1.1) becomes much stronger than the non-relativistic case so that its analytic treatment is extremely difficult and challenging even though some literature studied the exact solution to the RHD equations in the special case, for instance the one-dimensional Riemann problem or the isentropic problem [19, 27, 17]. A primary and powerful approach to improve our understanding of the physical mechanisms in the RHDs is the numerical simulation. In comparison with the non-relativistic case, the numerical difficulties are coming from strongly nonlinear coupling between the RHD equations (1.1), which leads to no explicit expression of the primitive variable vector \mbox{\boldmath\smallV}=(\rho,\mbox{\boldmath\smallu},p)^{T} and the flux \mbox{\boldmath\smallF}_{\ell} in terms of $U$, and some physical constraints such as $\rho>0,p>0$, $E\geq D$ and |\mbox{\boldmath\smallu}|<c=1, etc.

The pioneering work in this field may trace back to the finite difference code via the artificial viscosity for the spherically symmetric general RHD equations in the Lagrangian coordinates [23, 24]. The first attempt to solve RHD equations numerically in the Eulerian coordinates was made in [35] by using an explicit finite difference method with monotonic transport and artificial viscosity technique. After those, the numerical study of the RHD equations has attracted more attention, and various exact or approximate Riemann solvers have been successively proposed, for instance the HLL (Harten-Lax-van Leer) method [31], the flux corrected transport method [8], the two-shock solvers [1, 5], the Roe solver [9], the upwind scheme [10], the kinetic schemes [46, 15], the flux-splitting method [7], the HLLC (HLL-Contact) scheme [25], and so forth. In addition, some higher-order accurate shock-capturing schemes were also developed for solving the RHD equations, such as ENO (essentially non-oscillatory) and weighted ENO (WENO) methods [6, 49, 33], finite volume local evolution Galerkin method [44], piecewise parabolic methods [20, 26], adaptive mesh refinement method [50], discontinuous Galerkin (DG) method [30], direct Eulerian GRP (generalized Riemann problem) schemes [47, 48, 45, 43], adaptive moving mesh methods [12, 13], space-time conservation element and solution element method [28] and Runge-Kutta DG methods with WENO limiter [54] and so on. The readers are also referred to the review articles [21, 11, 22] and references therein.

Unfortunately, in general, it cannot be proved that the above schemes are PCP, in other words, their solutions belong to the physically admissible set {\mathcal{G}}=\{\mbox{\boldmath\smallU}=(D,\mbox{\boldmath\smallm},E)^{T}~{}\big{|} \rho>0,~{}p>0,~{}|\mbox{\boldmath\smallu}|<c=1\}. Although they have been used to simulate some RHD flows successfully, they are still confronted with enormous risks of the calculation failure in solving the RHD problems with large Lorentz factor or low density or pressure, or strong discontinuity. If the negative density or pressure, or the larger velocity than the speed of light are numerically obtained, the eigenvalues of the Jacobian matrix or the Lorentz factor may become imaginary, leading directly to the ill-posedness of the discrete problem. In practice, such nonphysical numerical solutions are usually simply replaced with a “close” and “physical” one by performing recalculation with more diffusive schemes and smaller CFL number until the numerical solutions become physical, see e.g. [14, 50]. Obviously, to some extent such operations are not scientifically reasonable, so it is necessary and significant to develop high-order accurate physical-constraints-preserving (PCP) numerical schemes, whose solutions can satisfy those intrinsic physical constraints or properties: the positivity of the rest-mass density and the pressure and the bound of the velocity. Recently, there exist some works on the PCP numerical schemes for the RHD equations, such as the finite difference WENO schemes [39], non-central DG methods with bound-preserving limiter [29], and central DG methods [40]. Moreover, the admissible states and physical-constraints-preserving numerical schemes were first well-studied for the special relativistic magnetohydrodynamics [41, 42], and the provably physical-constraint-preserving methods were successfully designed for the general RHDs [36]. Motivated by [41, 42], positivity-preserving numerical schemes was analyzed for the non-relativistic ideal magnetohydrodynamics in [37, 38]. It is worthy of mentioning that all of the aforementioned schemes with the PCP property are formulated in the Eulerian framework.

The aim of this paper is to study the PCP Lagrangian schemes with the HLLC solver for the one- and two-dimensional RHD equations (1.1) with $d=1,2$. The HLLC approximate Riemann solver has a simpler form in the Lagrangian framework. First, we prove that the intermediate states in the HLLC Riemann solver are admissible or PCP (that is, the rest-mass density and pressure are positive and the fluid velocity magnitude is less than the speed of light) when the HLLC wave speeds are estimated suitably. Next, we derive the first-order PCP Lagrangian scheme with the HLLC solver and then propose the higher-order accurate PCP Lagrangian scheme via the WENO reconstruction and the scaling PCP limiter as well as the strong stability preserving (SSP) high order time discretization.

The paper is organized as follows. Section 2 gives the 1D PCP Lagrangian schemes. The PCP properties of intermediate states in the HLLC Riemann solver are proved, and then the first- and high-order PCP Lagrangian schemes are proposed. Section 3 presents the 2D PCP Lagrangian schemes. It needs more techniques to derive the PCP properties of the HLLC Riemann solver for the $x_{i}$-split 2D RHD equations, because an important difference from the 1D RHD equations comes from a purely multi-dimensional relativistic feature that the flow regions across the shock or rarefaction wave in the split 2D RHD equations are nonlinearly coupled through the Lorentz factor which is also built in terms of the tangential velocities. Some numerical experiments are conducted in Section 4 to demonstrate the performance such as the conservation, PCP property, accuracy of our Lagrangian schemes. Finally some remarkable conclusions are summarized in Section 5.

2 One-dimensional PCP Lagrangian schemes

This section considers the Lagrangian finite volume schemes for the special relativistic hydrodynamics (RHD) equations (1.1) or (1.3) with $d=1$, here the notations \mbox{\boldmath\smallF}_{1}, $u_{1}$ and $x_{1}$ are replaced with $F$, $u$ and $x$, respectively, and \mbox{\boldmath\small\mathcal{F}}_{\mbox{\boldmath\smalln}}(\mbox{\boldmath\smallU})=\left(0,p,pu\right)^{T}=:\mbox{\boldmath\small\mathcal{F}}(\mbox{\boldmath\smallU}).

Assume that the time interval $\{t>0\}$ is discretized as: $t_{n+1}=t_{n}+\Delta t^{n}$, $n=0,1,2,\cdots$, where $\Delta t^{n}$ is the time step size at $t=t_{n}$ and will be determined by the CFL type condition. The (dynamic) computational domain at $t=t_{n}$ is divided into $N$ cells: $I_{i}^{n}=[x^{n}_{i-\frac{1}{2}},x^{n}_{i+\frac{1}{2}}]$ with the sizes $\Delta x_{i}^{n}=x^{n}_{i+\frac{1}{2}}-x^{n}_{i-\frac{1}{2}}$ for $i=1,\cdots,N$. The mesh node $x^{n}_{i+\frac{1}{2}}$ is moved with the fluid velocity $u^{n}_{i+\frac{1}{2}}$, that is,

[TABLE]

For the RHD system (1.1) or (1.3) with $d=1$, the Lagrangian finite volume scheme with the Euler forward time discretization

[TABLE]

where \overline{\mbox{\boldmath\smallU}}_{i}^{n} and \overline{\mbox{\boldmath\smallU}}_{i}^{n+1} are approximate cell average values of the conservative vectors at $t_{n}$ and $t_{n+1}$ over the cells $I_{i}^{n}$ and $I_{i}^{n+1}$, respectively, the numerical flux \widehat{\mbox{\boldmath\small\mathcal{F}}} is a function of two variables and satisfies the consistency \widehat{\mbox{\boldmath\small\mathcal{F}}}(\mbox{\boldmath\smallU},\mbox{\boldmath\smallU})=\left(0,p,pu\right)^{T}, and \mbox{\boldmath\smallU}_{i+\frac{1}{2}}^{\mp} are the left- and right-limited approximations of $U$ at $x_{i+\frac{1}{2}}$, respectively, and obtained by using the initial reconstruction technique and \{\overline{\mbox{\boldmath\smallU}}_{i}^{n}\}.

In this paper, the numerical flux \widehat{\mbox{\boldmath\small\mathcal{F}}}(\mbox{\boldmath\smallU}^{-},\mbox{\boldmath\smallU}^{+}) is computed by means of the HLLC Riemann solver of (1.1) with $d=1$ [25], but in a coordinate system that is moving with the intermediate wave (i.e. contact discontinuity) speed $s^{*}$ (see below). Let us consider the Riemann problem

[TABLE]

where $\xi=x-s^{*}t$ and

[TABLE]

Figure 2.1 shows corresponding wave structure of the exact or HLLC approximate solution of the Riemann problem (2.2). There are three wave families associated with the eigenvalues of \partial\widetilde{\mbox{\boldmath\smallF}}/\partial\mbox{\boldmath\smallU}, and the contact discontinuity is always in accordance with the vertical axis $\xi=0$. Two constant intermediate states between the left and right acoustic waves are denoted by \mbox{\boldmath\smallU}^{\ast,-} and \mbox{\boldmath\smallU}^{\ast,+}, respectively, corresponding fluxes are denoted by \widetilde{\mbox{\boldmath\smallF}}^{\ast,\pm}. The symbols $s^{-}$ and $s^{+}$ denote the smallest and largest speeds of the acoustic waves, respectively, and are specifically set as

[TABLE]

where s_{\max}(\mbox{\boldmath\smallU}) and s_{\min}(\mbox{\boldmath\smallU}) are the maximum and minimum eigenvalues of the Jacobian matrix \partial{\mbox{\boldmath\smallF}}/\partial{\mbox{\boldmath\smallU}}, and can be explicitly given by

[TABLE]

Integrating the conservation law in (2.2) over the control volume $[\xi_{L},0]\times[t_{n},t_{n+1}]$, the left portion of Figure 2.1, gives the flux \widetilde{\mbox{\boldmath\smallF}}^{\ast,-} along the $t$-axis as follows

[TABLE]

Similarly, performing the same operation on the control volume $[0,\xi_{R}]\times[t_{n},t_{n+1}]$ gives

[TABLE]

By means of the the definitions of the intermediate states \mbox{\boldmath\smallU}^{\ast,\pm} [34]

[TABLE]

one can rewrite (2.6) and (2.7) as

[TABLE]

which are the Rankine-Hugoniot conditions for the right and left waves in Figure 2.1. If assuming \mbox{\boldmath\smallU}^{\ast,\pm}=(D^{\ast,\pm},m^{\ast,\pm},E^{\ast,\pm})^{T} and \widetilde{\mbox{\boldmath\smallF}}^{\ast,\pm}=\big{(}D^{\ast,\pm}(u^{\ast,\pm}-s^{\ast}),m^{\ast,\pm}(u^{\ast,\pm}-s^{\ast})+p^{\ast,\pm},E^{\ast,\pm}(u^{\ast,\pm}-s^{\ast})+p^{\ast,\pm}u^{\ast,\pm}\big{)}^{T}, (2.8) can explicitly give

[TABLE]

From those, one obtains

[TABLE]

where

[TABLE]

Because the system (2.9) is underdetermined, one can impose the conditions on the continuity of the pressure and velocity across the contact discontinuity, i.e.

[TABLE]

Combing (2.10) with (2.12) gives

[TABLE]

Eliminating $p^{*}$ gives a quadratic equation in terms of $s^{\ast}$ as follows

[TABLE]

where

[TABLE]

From (2.14), one has

[TABLE]

since the wave speed $s^{\ast}$ should be less than the speed of light $c=1$ [25]. Then (2.13) further gives

[TABLE]

From (2.9), one can obtain \mbox{\boldmath\smallU}^{\ast,\pm} and then substitute \mbox{\boldmath\smallU}^{\ast,\pm} into (2.8) to give the HLLC flux (approximating \widetilde{\mbox{\boldmath\smallF}}(\mbox{\boldmath\smallU})) in the Lagrangian framework \widehat{\mbox{\boldmath\small\mathcal{F}}}=\widetilde{\mbox{\boldmath\smallF}}^{\ast,+}=\widetilde{\mbox{\boldmath\smallF}}^{\ast,-}. It is worth noticing that $p^{\ast}$ in (2.16) and $u^{\ast}$ in (2.15) are different from the pressure and the velocity calculated from the resulting \mbox{\boldmath\smallU}^{\ast,\pm} via (1.8) and (1.9), so in general p^{\ast}\neq p(\mbox{\boldmath\smallU}^{\ast,\pm}), u^{\ast}\neq u(\mbox{\boldmath\smallU}^{\ast,\pm}) and \widetilde{\mbox{\boldmath\smallF}}^{\ast,\pm}\neq\widetilde{\mbox{\boldmath\smallF}}(\mbox{\boldmath\smallU}^{\ast,\pm}).

For the above HLLC solver, the following conclusion holds.

Lemma 2.1**.**

For the given \mbox{\boldmath\smallU}^{\pm}=(D^{\pm},m^{\pm},E^{\pm})^{T}, if the wave speed $s^{\pm}$ are estimated in (2.4) and $s^{\ast}$ is the speed of the contact discontinuity, then $A^{\pm}$ and $B^{\pm}$ defined in (2.11) satisfy the following inequalities

1. (i)

$A^{-}<0,~{}~{}A^{+}>0$; 2. (ii)

$A^{-}-B^{-}<0,~{}~{}~{}A^{+}+B^{+}>0$; 3. (iii)

$s^{-}A^{-}-B^{-}>0,~{}~{}s^{+}A^{+}-B^{+}>0$; 4. (iv)

$s^{+}A^{-}-B^{-}<0,~{}~{}s^{-}A^{+}-B^{+}<0$; 5. (v)

$s^{-}<u^{\pm}<s^{+},~{}~{}s^{-}<s^{\ast}<s^{+}$; 6. (vi)

$4(A^{\pm})^{2}-(s^{\pm}A^{\pm}+B^{\pm})^{2}>0$; 7. (vii)

$(A^{\pm})^{2}-(B^{\pm})^{2}-(D^{\pm})^{2}(s^{\pm}-u^{\pm})^{2}>0$.

The proof is presented in Appendix A.

2.1 First-order accurate scheme

For the first-order accurate Lagrangian scheme, \mbox{\boldmath\smallU}_{i\pm\frac{1}{2}}^{\mp} are calculated from the reconstructed piecewise constant function, namely, \mbox{\boldmath\smallU}_{i\pm\frac{1}{2}}^{\mp} in the scheme (2.1) are defined by

[TABLE]

Those form some local Riemann problems at the endpoints of each cell naturally, see the diagrammatic sketch in Figure 2.2. The endpoints of the cell will be evolved by

[TABLE]

Similar to the Godunov scheme, under a suitable CFL condition (see the coming Theorem 2.3), the waves in two neighboring local Riemann problems (e.g. centered at the points $x_{i-\frac{1}{2}}^{n}$ and $x_{i+\frac{1}{2}}^{n}$) do not interact with each other within a time step, so \overline{\mbox{\boldmath\smallU}}^{n+1}_{i} in the scheme (2.1) can be equivalently derived by exactly integrating the approximate Riemann solutions over the cell $[x^{n+1}_{i-\frac{1}{2}},x^{n+1}_{i+\frac{1}{2}}]$, i.e.

[TABLE]

where R_{h}(x/t,\overline{\mbox{\boldmath\smallU}}^{n}_{j-1},\overline{\mbox{\boldmath\smallU}}^{n}_{j}) is the approximate Riemann solution related to the states \overline{\mbox{\boldmath\smallU}}^{n}_{j-1} and \overline{\mbox{\boldmath\smallU}}^{n}_{j}, $j=i,i+1$. For the above HLLC Riemann solver, the integration of R_{h}(x/t,\overline{\mbox{\boldmath\smallU}}^{n}_{i-1},\overline{\mbox{\boldmath\smallU}}^{n}_{i}) will be equal to (x_{1}-x_{i-\frac{1}{2}}^{n+1})\mbox{\boldmath\smallU}^{\ast,+} with \mbox{\boldmath\smallU}^{\ast,+} is computed from \overline{\mbox{\boldmath\smallU}}^{n}_{i-1} and \overline{\mbox{\boldmath\smallU}}^{n}_{i}, while the integration of R(x/t,\overline{\mbox{\boldmath\smallU}}^{n}_{i},\overline{\mbox{\boldmath\smallU}}^{n}_{i+1}) will be (x_{i+\frac{1}{2}}^{n+1}-x_{2})\mbox{\boldmath\smallU}^{\ast,-} with \mbox{\boldmath\smallU}^{\ast,-} derived from \overline{\mbox{\boldmath\smallU}}^{n}_{i} and \overline{\mbox{\boldmath\smallU}}^{n}_{i+1}. Thus the first-order HLLC scheme is equivalent to

[TABLE]

Therefore in order to prove that the first-order HLLC scheme (2.20) is PCP, one just needs to show the intermediate states \mbox{\boldmath\smallU}^{\ast,-},\mbox{\boldmath\smallU}^{\ast,+}\in\mathcal{G} due to the convexity of $\mathcal{G}$.

Theorem 2.2**.**

For \mbox{\boldmath\smallU}^{\pm}\in\mathcal{G} and the wave speeds $s^{\pm}$ defined in (2.4), the intermediate states \mbox{\boldmath\smallU}^{\ast,\pm} obtained by the HLLC Riemann solver are PCP, that is to say, \mbox{\boldmath\smallU}^{\ast,\pm}\in\mathcal{G}.

Proof.

Following [39], one has to prove $D^{\ast,\pm}>0$, $E^{\ast,\pm}>0$, and $(E^{\ast,\pm})^{2}-(D^{\ast,\pm})^{2}-(m^{\ast,\pm})^{2}>0$.

(i) Due to (iv) of Lemma 2.1 and the first equality in (2.9), it is easy to get $D^{\ast,\pm}>0$.

(ii) For the left or right state, from (2.9) and (2.10), one has

[TABLE]

where $A^{\pm}$ and $B^{\pm}$ are defined in (2.11). Using (i) of Lemma 2.1 gives

[TABLE]

Combining (2.21), (2.22) with (i) and (v) in Lemma 2.1 further yields $\mathcal{A}^{-}<0$ and $\mathcal{A}^{+}>0$, which imply

[TABLE]

(iii) Define

[TABLE]

Using (2.10) gives

[TABLE]

where

[TABLE]

The conclusion (vi) in Lemma 2.1 can tells us $\mathcal{K}^{\pm}>0$ and then $\mathcal{B}^{\pm}>0$, which imply

[TABLE]

The proof is completed.

Based the above discussion, one can draw the following conclusion.

Theorem 2.3**.**

If \{\overline{\mbox{\boldmath\smallU}}_{i}^{n}\in\mathcal{G},\forall i=1,\cdots,N\} and the wave speeds $s^{\pm}$ estimated by (2.4), then \overline{\mbox{\boldmath\smallU}}_{i}^{n+1} obtained by the first-order Lagrangian scheme (2.1) with (2.17) and HLLC solver belong to the admissible state set $\mathcal{G}$ under the following time step restriction

[TABLE]

where the CFL number $\lambda\leq\frac{1}{2}$.

2.2 High-order accurate scheme

This section will develop the one-dimensional high-order accurate PCP Lagrangian finite volume scheme with the HLLC solver. It consists of three parts: the high-order accurate initial reconstruction, the scaling PCP limiter, and the high-order accurate time discretization.

For the known cell-average values \{\overline{\mbox{\boldmath\smallU}}_{i}^{n}\} of the solutions of the RHD equations (1.1) or (1.3) with $d=1$, by the aid of the local characteristic decomposition [54], the WENO reconstruction technique is applied to get the high-order approximations of the point values \mbox{\boldmath\smallU}_{i-\frac{1}{2}}^{+} and \mbox{\boldmath\smallU}_{i+\frac{1}{2}}^{-}, and then they can be used to give the HLLC flux \widetilde{\mbox{\boldmath\smallF}}^{\ast,\pm} and intermediate states \mbox{\boldmath\smallU}^{\ast,\pm}. For a $(K+1)$th-order accurate finite volume WENO scheme, as soon as the point values \mbox{\boldmath\smallU}_{i-\frac{1}{2}}^{+} and \mbox{\boldmath\smallU}_{i+\frac{1}{2}}^{-} are obtained by the WENO reconstruction at $t=t_{n}$, one can also give a polynomial vector \mbox{\boldmath\smallU}_{i}^{n}(x) of degree $K$ in principle such that \mbox{\boldmath\smallU}_{i}^{n}(x_{i-\frac{1}{2}})=\mbox{\boldmath\smallU}_{i-\frac{1}{2}}^{+},\mbox{\boldmath\smallU}_{i}^{n}(x_{i+\frac{1}{2}})=\mbox{\boldmath\smallU}_{i+\frac{1}{2}}^{-}, its cell average value over $I_{i}^{n}$ is equal to \overline{\mbox{\boldmath\smallU}}_{i}^{n}, and \mbox{\boldmath\smallU}^{n}_{i}(x) is a $(K+1)$th-order accurate approximation to \mbox{\boldmath\smallU}(x,t_{n}) in $I_{i}^{n}$. Such a polynomial vector can be gotten by using the Hermite type reconstruction technique [51] and satisfy exactly the $L$-point Gauss-Lobatto quadrature rule with $2L-3\geq K$, i.e.

[TABLE]

which gives a split of \overline{\mbox{\boldmath\smallU}}_{i}^{n}, where $\omega_{1},\cdots,\omega_{L}$ are the quadrature weights satisfying $\omega_{1}=\omega_{L}>0$. Practically, it does not need to explicitly obtain such a polynomial vector since the mean value theorem tell us that there exists some $x^{\ast}\in I_{i}^{n}$ such that \mbox{\boldmath\smallU}_{i}^{n}(x^{\ast})=\frac{1}{1-2\omega_{1}}\sum\limits_{\alpha=2}^{L-1}\omega_{\alpha}\mbox{\boldmath\smallU}_{i}^{n}(x_{i}^{\alpha})=:\mbox{\boldmath\smallU}^{\ast\ast}_{i} [53]. At this time, the split (2.25) can be simply replaced with

[TABLE]

and \mbox{\boldmath\smallU}^{\ast\ast}_{i} can be calculated by

[TABLE]

By adding and subtracting the term \Delta t^{n}\widehat{\mbox{\boldmath\small\mathcal{F}}}(\mbox{\boldmath\smallU}_{i-\frac{1}{2}}^{+},\mbox{\boldmath\smallU}_{i+\frac{1}{2}}^{-}) and using the split (2.26), the scheme (2.1) with high-order WENO reconstruction can be reformulated as follows

[TABLE]

where

[TABLE]

Because $\mathcal{H}_{1}$ and $\mathcal{H}_{L}$ do exactly mimic the first-order scheme (2.1) with (2.17) and the HLLC solver, one can know $\mathcal{H}_{1},\mathcal{H}_{L}\in\mathcal{G}$ if the two boundary values \mbox{\boldmath\smallU}_{i-\frac{1}{2}}^{+},\mbox{\boldmath\smallU}_{i+\frac{1}{2}}^{-}\in\mathcal{G}, the wave speeds $s^{\pm}$ estimated by (2.4), and the time stepsize $\Delta t^{n}$ satisfies the restriction

[TABLE]

with $\lambda\leq\frac{1}{2}$. Therefore, besides the time step restriction (2.30), the sufficient condition for \overline{\mbox{\boldmath\smallU}}_{i}^{n+1}\in\mathcal{G} is

[TABLE]

which can be ensured by using a scaling limiter, presented in the next subsection. Hence the aforementioned results can be summarized as follows.

Theorem 2.4**.**

If \mbox{\boldmath\smallU}_{i-\frac{1}{2}}^{+},\mbox{\boldmath\smallU}_{i+\frac{1}{2}}^{-},\mbox{\boldmath\smallU}_{i}^{\ast\ast}\in\mathcal{G} and the wave speeds $s^{\pm}$ estimated by (2.4), then \overline{\mbox{\boldmath\smallU}}_{i}^{n+1} obtained by the high-order Lagrangian scheme (2.1) with the WENO reconstruction and the HLLC solver belongs to the admissible state set $\mathcal{G}$ under the time stepsize restriction (2.30) with $\lambda\leq\frac{1}{2}$.

2.2.1 Scaling PCP limiter

This section uses the scaling PCP limiter, which has been used in [3, 52, 40], to limit \mbox{\boldmath\smallU}_{i-\frac{1}{2}}^{+},\mbox{\boldmath\smallU}_{i+\frac{1}{2}}^{-},\mbox{\boldmath\smallU}_{i}^{\ast\ast} such that the limited values \widetilde{\mbox{\boldmath\smallU}}_{i-\frac{1}{2}}^{+},\widetilde{\mbox{\boldmath\smallU}}_{i+\frac{1}{2}}^{-},\widetilde{\mbox{\boldmath\smallU}}_{i}^{\ast\ast} belong to $\mathcal{G}$ when \overline{\mbox{\boldmath\smallU}}_{i}^{n}\in\mathcal{G}. For the sake of brevity, the superscript $n$ will be omitted in this section and a small parameter $\varepsilon$ is taken as $10^{-13}$. Such limiting procedure can be implemented as follows.

First, let us enforce the positivity of the mass density. For each cell $I_{i}$, define the limiter

[TABLE]

and then limit the mass density as follows

[TABLE]

Define \widehat{\mbox{\boldmath\smallU}}_{i\pm\frac{1}{2}}^{\mp}=(\widehat{D}_{i\pm\frac{1}{2}}^{\mp},m_{i\pm\frac{1}{2}}^{\mp},E_{i\pm\frac{1}{2}}^{\mp})^{\mathrm{T}} and \widehat{\mbox{\boldmath\smallU}}_{i}^{\ast\ast}=(\widehat{D}_{i}^{\ast\ast},m_{i}^{\ast\ast},E_{i}^{\ast\ast})^{\mathrm{T}}.

Then, enforce the positivity of q(\mbox{\boldmath\smallU})=E-\sqrt{D^{2}+m^{2}}. For each cell $I_{i}$, define the limiter

[TABLE]

and then limit the conservative vectors as follows

[TABLE]

It is easy to check that \widetilde{\mbox{\boldmath\smallU}}_{i-\frac{1}{2}}^{+},\widetilde{\mbox{\boldmath\smallU}}_{i+\frac{1}{2}}^{-},\widetilde{\mbox{\boldmath\smallU}}_{i}^{\ast\ast}\in\mathcal{G}. Moreover, the above scaling PCP limiter does not destroy the original high order accuracy in smooth regions, see the detailed discussion in [52].

2.2.2 High-order time discretization

To get a globally high-order accurate scheme in time and space, we can further employ the strong stability preserving (SSP) Runge-Kutta (RK) method to replace the explicit Euler time discretization in (2.1) and (2.18). Similar to [3], for instance, to obtain a third-order accurate scheme in time, a third-stage SSP, explicit RK method may be used for the time discretization as follows.

Stage 1:

[TABLE]

Stage 2:

[TABLE]

Stage 3:

[TABLE]

where \mathcal{L}(\overline{\mbox{\boldmath\smallU}};i)=\widehat{\mbox{\boldmath\small\mathcal{F}}}(\mbox{\boldmath\smallU}_{i+\frac{1}{2}}^{-},\mbox{\boldmath\smallU}_{i+\frac{1}{2}}^{+})-\widehat{\mbox{\boldmath\small\mathcal{F}}}(\mbox{\boldmath\smallU}_{i-\frac{1}{2}}^{-},\mbox{\boldmath\smallU}_{i-\frac{1}{2}}^{+}).

The scaling PCP limiter described in Section 2.2.1 needs to be performed at each stage of the above RK method to limit the value of {\mbox{\boldmath\smallU}}_{i-\frac{1}{2}}^{+},{\mbox{\boldmath\smallU}}_{i+\frac{1}{2}}^{-},{\mbox{\boldmath\smallU}}_{i}^{\ast\ast}. Because each stage of the above SSP RK method is a convex combination of the forward Euler time discretization and $\mathcal{G}$ is convex, so is the above SSP RK method when the forward Euler method is conservative, stable and PCP.

3 Two-dimensional PCP Lagrangian schemes

This section considers the Lagrangian finite volume schemes for the special relativistic hydrodynamics (RHD) equations (1.1) or (1.3) with $d=2$. Here the notations $(m_{1},m_{2})$, $(u_{1},u_{2})$ and $(x_{1},x_{2})$ are replaced with \mbox{\boldmath\smallm}=(m_{x},m_{y}), \mbox{\boldmath\smallu}=(u_{x},u_{y}) and \mbox{\boldmath\smallx}=(x,y), respectively.

Assume that the time interval $\{t>0\}$ is divided into: $t_{n+1}=t_{n}+\Delta t^{n}$, $n=0,1,2,\cdots$, where the time step size $\Delta t^{n}$ will be determined by the CFL type condition. The (dynamic) computational domain $\Omega$ at $t=t_{n}$ is partitioned into $N_{x}\times N_{y}$ quadrilateral cells: $I_{ij}^{n}$ with the boundary $\partial I_{ij}^{n}$ and four vertices $\{(x_{i\pm\frac{1}{2}}^{n},y_{j\pm\frac{1}{2}}^{n}),(x_{i\pm\frac{1}{2}}^{n},y_{j\mp\frac{1}{2}}^{n})\}$, $i=1,\cdots,N_{x},j=1,\cdots,N_{y}$; and then the conservative Lagrangian finite volume scheme with the forward Euler time discretization for the governing equations (1.1) with $d=2$ can be given as follows

[TABLE]

where \overline{\mbox{\boldmath\smallU}}^{n}_{ij}=(\overline{D}_{ij}^{n},\overline{\mbox{\boldmath\smallm}}_{ij}^{n},\overline{E}_{ij}^{n})^{T} is the cell average approximation of $U$ at $t_{n}$ over the cell $I^{n}_{ij}$, $A^{n}_{ij}=\int_{I^{n}_{ij}}dxdy$ denotes the area of $I^{n}_{ij}$, \mbox{\boldmath\smallU}^{-}_{ij}=(D^{-}_{ij},\mbox{\boldmath\smallm}^{-}_{ij},E^{-}_{ij})^{T} and \mbox{\boldmath\smallU}^{+}_{ij}=(D^{+}_{ij},\mbox{\boldmath\smallm}^{+}_{ij},E^{+}_{ij})^{T} are the reconstructed limits of $U$ on the boundary $\partial I_{ij}^{n}$ from the inside and outside of the cell $I^{n}_{ij}$, respectively; $\partial I^{n,m}_{ij}$ is the $m$th cell edge of $\partial I_{ij}^{n}$, and $|l_{ij}^{m}|$ and \mbox{\boldmath\smalln}^{m}_{ij} are its length and unit outward normal vector from the inside to the outside of $I^{n}_{ij}$, respectively. Here \widehat{\mbox{\boldmath\small\mathcal{F}}} denotes the numerical flux, satisfying the consistency condition \widehat{\mbox{\boldmath\small\mathcal{F}}}_{\mbox{\boldmath\smalln}}(\mbox{\boldmath\smallU},\mbox{\boldmath\smallU})=(0,p\mbox{\boldmath\smalln}^{T},p\mbox{\boldmath\smallu}\cdot\mbox{\boldmath\smalln}^{T})^{T} and the conservation property \widehat{\mbox{\boldmath\small\mathcal{F}}}_{\mbox{\boldmath\smalln}}(\mbox{\boldmath\smallU}^{-},\mbox{\boldmath\smallU}^{+})=-\widehat{\mbox{\boldmath\small\mathcal{F}}}_{-\mbox{\boldmath\smalln}}(\mbox{\boldmath\smallU}^{+},\mbox{\boldmath\smallU}^{-}). The flux \widehat{\mbox{\boldmath\small\mathcal{F}}}_{\mbox{\boldmath\smalln}_{ij}^{m}} will be obtained by solving the local 1D Riemann problem along the vector normal to the edge $\partial I^{n,m}_{ij}$ with the HLLC Riemann solver. Different from the one-dimensional numerical flux \widehat{\mbox{\boldmath\small\mathcal{F}}}(\mbox{\boldmath\smallU}^{-},\mbox{\boldmath\smallU}^{+}) presented in Section 2, the flux \widehat{\mbox{\boldmath\small\mathcal{F}}}_{\mbox{\boldmath\smalln}_{ij}^{m}}(\mbox{\boldmath\smallU}^{-},\mbox{\boldmath\smallU}^{+}) has to be obtained by the HLLC Riemann solver of the split 2D system for (1.1) with $d=2$ in a moving coordinate system. The readers are referred to [48] for the emphasis on the differences between the system (1.1) with $d=1$ and the split 2D system for (1.1) with $d=2$. In view of that the current HLLC Riemann solver is essentially the same as that in Section 2 except for the special attention to the nonlinear coupling from the tangential velocity, some details of the current HLLC Riemann solver will be omitted here.

Let \mbox{\boldmath\smallU}^{\mp}=(D^{\mp},\mbox{\boldmath\smallm}^{\mp},E^{\mp})^{T} be the left and right states of the local Riemann problem of the split 2D system for (1.1) with $d=2$ in the direction $n$ as shown in Figure 3.1. At present, the wave structure of the exact or HLLC approximate solution of the local Riemann problem is the same as that in Figure 2.1. Two constant intermediate states between the left and right acoustic waves are denoted by \mbox{\boldmath\smallU}^{\ast,\pm}=(D^{\ast,\pm},\mbox{\boldmath\smallm}^{\ast,\pm},E^{\ast,\pm})^{T}, respectively, and the fluxes are denoted by \widetilde{\mbox{\boldmath\smallF}}^{\ast,\pm}_{\mbox{\boldmath\smalln}}, then for the left or right wave, the Rankine-Hugoniot conditions \widetilde{\mbox{\boldmath\smallF}}^{\ast,\pm}_{\mbox{\boldmath\smalln}}=\widetilde{\mbox{\boldmath\smallF}}_{\mbox{\boldmath\smalln}}(\mbox{\boldmath\smallU}^{\pm})-(s^{\pm}-s^{\ast})(\mbox{\boldmath\smallU}^{\pm}-\mbox{\boldmath\smallU}^{\ast,\pm}) can be similarly given in the component form

[TABLE]

where $m_{n}^{\pm}$ and $m_{\tau}^{\mp}$ are the normal and tangential components of \mbox{\boldmath\smallm}^{\mp}, i.e. m_{n}=\mbox{\boldmath\smallm}^{T}\cdot\mbox{\boldmath\smalln} and m_{\tau}=\mbox{\boldmath\smallm}^{T}\cdot\bm{\tau}, $u_{n}=m_{n}/(E+p)$ is the normal component of the velocity vector, and the wave speeds $s^{\pm}$ are estimated as follows

[TABLE]

with

[TABLE]

and $\sigma_{s}=c_{s}^{2}/[\gamma^{2}(1-c_{s}^{2})]$.

Combining the last equation of (3.2) with the second one gives the expression of $s^{\ast}$ in terms of $p^{\ast}$ as follows

[TABLE]

where

[TABLE]

Similar to (2.14), solving the system (3.5) gives

[TABLE]

where

[TABLE]

Based on those, the HLLC intermediate states \mbox{\boldmath\smallU}^{\ast,\pm} can be gotten from (3.2) and then the HLLC fluxes \widetilde{\mbox{\boldmath\smallF}}^{\ast,\pm}_{\mbox{\boldmath\smalln}} (approximating \widetilde{\mbox{\boldmath\smallF}}_{\mbox{\boldmath\smalln}}(\mbox{\boldmath\smallU})) and the numerical flux in the Lagrangian framework \widehat{\mbox{\boldmath\small\mathcal{F}}}_{\mbox{\boldmath\smalln}}=\widetilde{\mbox{\boldmath\smallF}}^{\ast,+}_{\mbox{\boldmath\smalln}}=\widetilde{\mbox{\boldmath\smallF}}^{\ast,-}_{\mbox{\boldmath\smalln}}.

Similar to Lemma 2.1, the following conclusion can be drawn.

Lemma 3.1**.**

For \mbox{\boldmath\smallU}^{\pm}=(D^{\pm},\mbox{\boldmath\smallm}^{\pm},E^{\pm})^{T} and the wave speeds $s^{\pm}$ estimated in (3.3), if $s^{\ast}$ is the speed of the contact discontinuity and $A^{\pm},B^{\pm}$ are defined by (3.6), then one has

1. (i)

$A^{-}<0,~{}~{}A^{+}>0$; 2. (ii)

$A^{-}-B^{-}<0,~{}~{}A^{+}+B_{+}>0$; 3. (iii)

$s^{-}A^{-}-B^{-}>0,~{}~{}s^{+}A^{+}-B^{+}>0$; 4. (iv)

$s^{+}A^{-}-B^{-}<0,~{}~{}s^{-}A^{+}-B^{+}<0$; 5. (v)

$s^{-}<u_{n}^{\pm}<s^{+},~{}~{}s^{-}<s^{\ast}<s^{+}$; 6. (vi)

$4(A^{\pm})^{2}-(s^{\pm}A^{\pm}+B^{\pm})^{2}>0$; 7. (vii)

$(A^{\pm})^{2}-(B^{\pm})^{2}-(D^{\pm})^{2}(s^{\pm}-u_{n}^{\pm})^{2}-(m_{\tau}^{\pm})^{2}(s^{\pm}-u_{n}^{\pm})^{2}>0$.

The proof is given in Appendix B. Based on this result, one can further get the following theorem, which corresponds to Theorem 2.2.

Theorem 3.2**.**

For \mbox{\boldmath\smallU}^{\pm}\in\mathcal{G}, and the wave speeds $s^{\pm}$ defined in (3.3), the intermediate states \mbox{\boldmath\smallU}^{\ast,\pm} obtained by the above HLLC Riemann solver are PCP, that is to say, \mbox{\boldmath\smallU}^{\ast,\pm}\in\mathcal{G}.

Proof.

Similarly, one has to prove $D^{\ast,\pm}>0$, $E^{\ast,\pm}>0$, and (E^{\ast,\pm})^{2}-(D^{\ast,\pm})^{2}-|\mbox{\boldmath\smallm}^{\ast,\pm}|^{2}>0.

(i) Due to the first equation in (3.2) and (v) in Lemma 3.1, it is easy to know that $D^{\ast,\pm}>0$.

(ii) According to (3.2) and (3.5), for the left or the right state one has

[TABLE]

Using the conclusion (i) in Lemma 3.1 gives

[TABLE]

Combining (3.10), (3.11) with Lemma 3.1 can easily show $\mathcal{A}^{-}<0$ and $\mathcal{A}^{+}>0$, which implies $E^{\ast,-}>0$ and $E^{\ast,+}>0$.

(iii) Define \mathcal{B}^{\pm}=\big{[}(E^{\ast,\pm})^{2}-(D^{\ast,\pm})^{2}-|\mbox{\boldmath\smallm}^{\ast,\pm}|^{2}\big{]}(s^{\pm}-s^{\ast})^{2}(A^{\pm}+s^{\pm}p^{\ast})^{2}, and then use (3.2) and (3.5) to give

[TABLE]

where

[TABLE]

Finally, using (vii) in Lemma 3.1 gives $\mathcal{K}^{\pm}>0$. Thus one has $\mathcal{B}^{\pm}>0$, which implies (E^{\ast,\pm})^{2}-(D^{\ast,\pm})^{2}-|\mbox{\boldmath\smallm}^{\ast,\pm}|^{2}>0. The proof is completed.

3.1 First-order accurate scheme

For the first-order accurate Lagrangian scheme, \mbox{\boldmath\smallU}_{m,ij}^{-} and \mbox{\boldmath\smallU}_{m,ij}^{+} in (3.1) are taken as the cell average values of the conservative vector $U$ in the cell $I_{ij}^{n}$ and its neighboring cell sharing the cell edge $\partial I_{ij}^{n,m}$, $m=1,2,3,4$. The scheme (1.3) becomes

[TABLE]

where \overline{\mbox{\boldmath\smallU}}_{m}^{\text{ext}(I^{n}_{ij})} is the cell average of $U$ over the neighboring cell of $I^{n}_{ij}$ sharing the edge $\partial I^{n,m}_{ij}$ with $I^{n}_{ij}$. The vertices of the cell will be evolved by

[TABLE]

where the velocity \mbox{\boldmath\smallu}_{i+\frac{1}{2},j+\frac{1}{2}}^{n} is calculated as follows:

[TABLE]

where \mbox{\boldmath\smallu}_{i+\frac{1}{2},\ell}^{\ast} and \mbox{\boldmath\smallu}_{\ell^{\prime},j+\frac{1}{2}}^{\ast}, $\ell=j,j+1$, $\ell^{\prime}=i,i+1$, are fluid velocities at midpoints of four cell edges with a common vertex (e.g. the point P in Figure 3.1), respectively. Here take the calculation of \mbox{\boldmath\smallu}_{i+\frac{1}{2},j}^{\ast} as an example. The velocity \mbox{\boldmath\smallu}_{i+\frac{1}{2},j}^{\ast} is gotten by using the local rotation transformation of $(u_{n},u_{\tau})_{i+\frac{1}{2},j}^{\ast}$, where $(u_{n})_{i+\frac{1}{2},j}^{\ast}=s_{i+\frac{1}{2},j}^{\ast}$, $(u_{\tau})_{i+\frac{1}{2},j}^{\ast}=\frac{1}{2}\left((u_{\tau})_{i,j}^{n}+(u_{\tau})_{i+1,j}^{n}\right)$, and $s_{i+\frac{1}{2},j}^{\ast}$ is the speed of contact discontinuity in the HLLC solver, see (3.7).

Because the flux \mbox{\boldmath\small\mathcal{F}}_{\mbox{\boldmath\smalln}}(\mbox{\boldmath\smallU}) in (1.3) can be written as follows

[TABLE]

and the geometrical relation \sum\limits_{m=1}^{4}\mbox{\boldmath\smalln}_{ij}^{m}|l_{ij}^{m}|=\mbox{\boldmath\small0} holds, utilizing the flux consistency gives

[TABLE]

Adding the identity \mbox{\boldmath\small0}=\Delta t^{n}\sum\limits_{m=1}^{4}\widehat{\mbox{\boldmath\small\mathcal{F}}}_{\mbox{\boldmath\smalln}_{ij}^{m}}(\overline{\mbox{\boldmath\smallU}}_{ij}^{n},\overline{\mbox{\boldmath\smallU}}_{ij}^{n})|l_{ij}^{m}| into the scheme (3.12) gives

[TABLE]

where

[TABLE]

It means that the scheme (3.12) can be expressed as a combination of $\mathcal{H}_{m}$ with positive coefficients, while $\mathcal{H}_{m}$ has the same form as the first order scheme (2.1) with (2.17) and HLLC solver. Thus the analysis in Section 2 can be applied to get the sufficient condition for that $\mathcal{H}_{m}$ is in the admissible state set $\mathcal{G}$. Summarizing those draws the following conclusion.

Theorem 3.3**.**

For the first order finite volume Lagrangian scheme (3.12), if \{\overline{\mbox{\boldmath\smallU}}_{ij}^{n}\in\mathcal{G},\forall i=1,\cdots,N_{x};j=1,\cdots,N_{y}\}, then we have \overline{\mbox{\boldmath\smallU}}_{ij}^{n+1}\in\mathcal{G} for all $i=1,\cdots,N_{x};j=1,\cdots,N_{y}$ under the wave speed estimates (3.3) and the following time step restriction

[TABLE]

where the CFL number $\lambda\leq\frac{1}{2}$.

3.2 High-order accurate scheme

In view of that the Lagrangian methods on the the quadrilateral grid with straight edge can be at most second order accurate anyway [2], this section focuses on developing the two-dimensional second-order accurate PCP Lagrangian finite volume schemes with the HLLC solver, based on the initial reconstruction, the scaling PCP limiter, and the second order Runge-Kutta time discretization.

The following discusses the initial reconstruction and the scaling PCP limiter. At this time, the two-dimensional Lagrangian finite volume scheme (3.1) should be replaced with

[TABLE]

which is derived by approximating the line integral in (3.1) via the $L$-point Gauss-Lobatto quadrature, where \{(\mbox{\boldmath\smallU}_{m,ij}^{\alpha})^{-},(\mbox{\boldmath\smallU}_{m,ij}^{\alpha})^{+}\} are the left and right limit values from the inside and outside of the cell $I^{n}_{ij}$, respectively, approximating the conservative variable $U$ at \mbox{\boldmath\smallx}_{m,ij}^{\alpha}, $\alpha=1,\dots,L$, which are the $L$ Gauss-Lobatto quadrature points mapped onto ${\partial I_{ij}^{n,m}}$. Those limit values are obtained by using the high-order WENO reconstruction used in [4]. To be specific, for the edge $\partial I^{n,1}_{ij}$ and the five known cell-average values \{\overline{\mbox{\boldmath\smallU}}_{i,j}^{n},\overline{\mbox{\boldmath\smallU}}_{i-1,j}^{n},\overline{\mbox{\boldmath\smallU}}_{i+1,j}^{n},\overline{\mbox{\boldmath\smallU}}_{i,j-1}^{n},\overline{\mbox{\boldmath\smallU}}_{i,j+1}^{n}\} of the solutions of the RHD equations (1.1) or (1.3) with $d=2$, we first transform those cell average values \{\overline{\mbox{\boldmath\smallU}}_{i,j}^{n}\} into $\{\overline{\mathbf{U}}_{i,j}^{n}\}$ by using the local rotational transformation from the $(x,y)$ coordinates to the local $(\xi,\eta)$ coordinates, where $\xi$ and $\eta$ are in the $\bm{n}^{n,1}_{ij}$ and $\bm{\tau}^{\tau,1}_{ij}$ directions, respectively, see Figure 3.2, and then calculate the values of corresponding characteristic variables \mbox{\boldmath\smallW}_{i,j}^{n}=\mbox{\boldmath\smallL}\overline{\mathbf{U}}_{i,j}^{n},\mbox{\boldmath\smallW}_{i-1,j}^{n}=\mbox{\boldmath\smallL}\overline{\mathbf{U}}_{i-1,j}^{n},\mbox{\boldmath\smallW}_{i+1,j}^{n}=\mbox{\boldmath\smallL}\overline{\mathbf{U}}_{i+1,j}^{n},\mbox{\boldmath\smallW}_{i,j+1}^{n}=\mbox{\boldmath\smallL}\overline{\mathbf{U}}_{i,j+1}^{n},\mbox{\boldmath\smallW}_{i,j-1}^{n}=\mbox{\boldmath\smallL}\overline{\mathbf{U}}_{i,j-1}^{n}, where \mbox{\boldmath\smallL}=\mbox{\boldmath\smallL}(\overline{\mbox{\boldmath\smallU}}_{i,j}^{n}) is the left eigen matrix of \partial\mbox{\boldmath\smallF}_{\mbox{\boldmath\smalln}}/\partial\mbox{\boldmath\smallU}. Let us consider the following four stencils $S_{1}=\{I^{n}_{i,j},I^{n}_{i+1,j},I^{n}_{i,j+1}\}$, $S_{2}=\{I^{n}_{i,j},I^{n}_{i-1,j},I^{n}_{i,j+1}\}$, $S_{3}=\{I^{n}_{i,j},I^{n}_{i+1,j},I^{n}_{i,j-1}\}$, and $S_{4}=\{I^{n}_{i,j},I^{n}_{i-1,j},I^{n}_{i,j-1}\}$.

For the stencil $S_{q}$, $q=1,2,3,4$, one can obtain a linear polynomial

[TABLE]

where $(\xi_{ij},\eta_{ij})$ is the barycenter of $I^{n}_{ij}$ and the coefficients \{\mbox{\boldmath\smalla}_{q},\mbox{\boldmath\smallb}_{q},\mbox{\boldmath\smallc}_{q}\} are determined by preserving the cell average values, e.g. for $q=1$

[TABLE]

Using those can give the final linear polynomial \widetilde{\mbox{\boldmath\smallW}}(\mbox{\boldmath\smallx})=\widetilde{\mbox{\boldmath\smallW}}(\xi,\eta)=\mbox{\boldmath\smalla}(\xi-\xi_{ij})+\mbox{\boldmath\smallb}(\eta-\eta_{ij})+\mbox{\boldmath\smallc} with the coefficients \{\mbox{\boldmath\smalla},\mbox{\boldmath\smallb},\mbox{\boldmath\smallc}\} determined by

[TABLE]

where the weights $\{{\omega}_{q}\}$ are defined by

[TABLE]

Using the polynomial \widetilde{\mbox{\boldmath\smallW}}(\mbox{\boldmath\smallx}) calculates its values at the Gauss-Lobatto point and then gives (\mathbf{U}_{1,ij}^{\alpha})^{-}=\mbox{\boldmath\smallL}^{-1}\widetilde{\mbox{\boldmath\smallW}}({\mbox{\boldmath\smallx}}^{\alpha}_{1,ij}). Finally, the values (\mbox{\boldmath\smallU}_{1,ij}^{\alpha})^{-} can be obtained by using the inverse rotation transformation. It is worth noting that as soon as the limit value $(\mathbf{U}_{m,ij}^{\alpha})^{-}$ and $(\mathbf{U}_{m,ij}^{\alpha})^{+}$ can give the HLLC flux and HLLC intermediate states corresponding to the local 1D Riemann problem located at the point \mbox{\boldmath\smallx}_{m,ij}^{\alpha}, $\alpha=1,\dots,L$.

To give a decomposition similar to (2.26), we use a coordinate transformation \mbox{\boldmath\smallx}=\mbox{\boldmath\smallx}_{ij}(\hat{x},\hat{y}) to transform the quadrilateral cell $I_{ij}^{n}$ in the $(x,y)$ plane to the unit square $\hat{I}_{0}=[-\frac{1}{2},\frac{1}{2}]\times[-\frac{1}{2},\frac{1}{2}]$ in the $(\hat{x},\hat{y})$ plane, see Figure 3.3, then define the set of the 2D Gauss-Lobatto quadrature points in the cell $I_{ij}^{n}$ by

[TABLE]

which are derived by inversely transformed the Gauss-Lobatto quadrature points in the unit square $\hat{I}_{0}$. Because only the second-order accurate scheme is considered here, one may apply the tensor product Simpson quadrature rule, in which the quadrature points consist of the cell vertices, the mid-points of each edge and the cell center, see Figure 3.3, i.e. $L=3$, $\omega_{1}=\omega_{3}=\frac{1}{6}$, and $\omega_{2}=\frac{2}{3}$. Based on those, the term \overline{\mbox{\boldmath\smallU}}^{n}_{ij}A^{n}_{ij} can be decomposed into

[TABLE]

for the four edges $\partial I^{n,m}_{ij}$, $m=1,2,3,4$, where $J_{ij}^{\alpha,\beta}=J_{ij}(\hat{x}_{\alpha},\hat{y}_{\beta})$, $J_{ij}(\hat{x},\hat{y})=\left|\frac{\partial(x_{ij}(\hat{x},\hat{y}),y_{ij}(\hat{x},\hat{y}))}{\partial(\hat{x},\hat{y})}\right|$ is the Jacobian for the coordinate transformation, and

[TABLE]

From (3.18), one has

[TABLE]

where $\widetilde{\omega}_{\alpha,\beta}=\omega_{\alpha}\omega_{\beta}J_{ij}^{\alpha,\beta}$.

Moreover, since the point values \mbox{\boldmath\smallU}_{1,ij}^{\alpha,1},\mbox{\boldmath\smallU}_{2,ij}^{3,\alpha},\mbox{\boldmath\smallU}_{3,ij}^{\alpha,3},\mbox{\boldmath\smallU}_{4,ij}^{1,\alpha} can be obtained from the WENO reconstruction, one can directly compute \mbox{\boldmath\smallU}_{m,ij}^{\ast\ast} from (3.18)

[TABLE]

Consequently, by adding and subtracting the term \Delta t^{n}\sum\limits_{m=2}^{4}\sum\limits_{\alpha=1}^{3}\omega_{\alpha}\widehat{\mbox{\boldmath\small\mathcal{F}}}_{\mbox{\boldmath\smalln}_{ij}^{m}}\big{(}(\mbox{\boldmath\smallU}_{1,ij}^{\alpha})^{-},(\mbox{\boldmath\smallU}_{m,ij}^{\alpha})^{-}\big{)}|l_{ij}^{m}| and using (3.20) and the fact that $\omega_{1}=\omega_{3}$, (\mbox{\boldmath\smallU}_{1,ij}^{\alpha})^{-}=\mbox{\boldmath\smallU}_{1,ij}^{\alpha,1}, (\mbox{\boldmath\smallU}_{2,ij}^{\alpha})^{-}=\mbox{\boldmath\smallU}_{2,ij}^{3,\alpha}, (\mbox{\boldmath\smallU}_{3,ij}^{\alpha})^{-}=\mbox{\boldmath\smallU}_{3,ij}^{\alpha,3}, and (\mbox{\boldmath\smallU}_{4,ij}^{\alpha})^{-}=\mbox{\boldmath\smallU}_{4,ij}^{1,\alpha}, the scheme (3.15) becomes

[TABLE]

where

[TABLE]

It is obvious that the equation of $\mathcal{H}_{1}^{\alpha}$ has the same type as the 2D first-order scheme (3.12), $\alpha=1,2,3$, while the equation of $\mathcal{H}_{m}^{\alpha}$ is similar to the 1D first-order scheme (2.1) with (2.17), $\alpha=1,2,3,m=2,3,4$. Meanwhile, \overline{\mbox{\boldmath\smallU}}^{n+1} is a convex combination of \mbox{\boldmath\smallU}_{m,ij}^{\ast\ast} and $\mathcal{H}_{m}^{\alpha}$, $\alpha=1,2,3,m=1,\cdots,4$. Thus, if those terms are PCP, then \overline{\mbox{\boldmath\smallU}}^{n+1}\in\mathcal{G} due to the convexity of the admissible state set $\mathcal{G}$.

Theorem 3.4**.**

If \overline{\mbox{\boldmath\smallU}}_{ij}^{n},\mbox{\boldmath\smallU}_{m,ij}^{\ast\ast}\in\mathcal{G} for all $m,i,j$, and the HLLC wave speeds are estimated in (3.3), then the high-order finite volume Lagrangian scheme (3.15) is PCP, i.e. \overline{\mbox{\boldmath\smallU}}_{ij}^{n+1}\in\mathcal{G} for all $i=1,\cdots,N_{x};j=1,\cdots,N_{y}$, under the following time stepsize restriction

[TABLE]

where the CFL number $\lambda\leq\frac{1}{2}$, $J_{ij}=\min\limits_{\alpha=1,\cdots,3}\{J_{ij}^{\alpha,1},J_{ij}^{3,\alpha},J_{ij}^{\alpha,3},J_{ij}^{1,\alpha}\}$, and

[TABLE]

Before ending this section, we discuss how to to limit \mbox{\boldmath\smallU}_{m,ij}^{\ast\ast} defined in (3.21) and (\mbox{\boldmath\smallU}_{m,ij}^{\alpha})^{-} reconstructed by the WENO technique such that the limited values \widetilde{\mbox{\boldmath\smallU}}_{m,ij}^{\ast\ast} and (\widetilde{\mbox{\boldmath\smallU}}_{m,ij}^{\alpha})^{-} belong to $\mathcal{G}$ when \overline{\mbox{\boldmath\smallU}}_{ij}^{n}\in\mathcal{G}. For the sake of brevity, the superscript $n$ will be omitted in this section and a small parameter $\varepsilon$ is taken as $10^{-13}$. Similar to the one-dimensional case, the scaling PCP limiter can be implemented as follows.

First, enforce the positivity of the mass density. For each cell $I_{ij}$, define

[TABLE]

and limit

[TABLE]

Define (\widehat{\mbox{\boldmath\smallU}}_{m,ij}^{\alpha})^{-}=((\widehat{D}_{m,ij}^{\alpha})^{-},(\mbox{\boldmath\smallm}_{m,ij}^{\alpha})^{-},(E_{m,ij}^{\alpha})^{-})^{T}.

Next, enforce the positivity of the term q(\mbox{\boldmath\smallU})=E-\sqrt{D^{2}+|\mbox{\boldmath\smallm}|^{2}}. For each cell $I_{ij}$, compute

[TABLE]

and then limit the point values

[TABLE]

It is easy to show that all those limited values are in the admissible state set $\mathcal{G}$ when \overline{\mbox{\boldmath\smallU}}_{ij}^{n}\in\mathcal{G}.

In order to get a 2D Lagrangian scheme of second order accuracy both in space and time, we replace the forward Euler time discretization with the second order Runge-Kutta time discretization in the scheme (3.15), which can be implemented as follows:

Stage 1:

[TABLE]

Stage 2:

[TABLE]

where \mathcal{L}(\overline{\mbox{\boldmath\smallU}};i,j)=\sum\limits_{m=1}^{4}\sum\limits_{\alpha=1}^{3}\omega_{\alpha}\widehat{\mbox{\boldmath\small\mathcal{F}}}_{\mbox{\boldmath\smalln}_{ij}^{m}}\big{(}(\mbox{\boldmath\smallU}_{m,ij}^{\alpha})^{-},(\mbox{\boldmath\smallU}_{m,ij}^{\alpha})^{+}\big{)}|l_{ij}^{m}|. Here, the velocity of the vertex \mbox{\boldmath\smallu}^{n}_{i+\frac{1}{2},j+\frac{1}{2}} can be computed by

[TABLE]

where \mbox{\boldmath\smallu}^{\ast,k}_{i+\frac{1}{2},j+\frac{1}{2}},k=1,\cdots,4 are fluid velocities at nodes of four cell edges sharing the common vertex \mbox{\boldmath\smallx}_{i+\frac{1}{2},j+\frac{1}{2}}^{n} (e.g. the point P in Figure 3.1), respectively. Here take the calculation of \mbox{\boldmath\smallu}^{\ast,1}_{i+\frac{1}{2},j+\frac{1}{2}} as an example. The velocity \mbox{\boldmath\smallu}_{i+\frac{1}{2},j+\frac{1}{2}}^{\ast,1} is gotten by using the local rotation transformation of $(u_{n},u_{\tau})_{i+\frac{1}{2},j+\frac{1}{2}}^{\ast,1}$, where $(u_{n})_{i+\frac{1}{2},j+\frac{1}{2}}^{\ast,1}=s_{i+\frac{1}{2},j+\frac{1}{2}}^{\ast,1}$, $(u_{\tau})_{i+\frac{1}{2},j+\frac{1}{2}}^{\ast,1}=\frac{1}{2}\left((u_{\tau}^{n})_{i+\frac{1}{2},j+\frac{1}{2}}^{-,1}+(u_{\tau}^{n})_{i+\frac{1}{2},j+\frac{1}{2}}^{+,1}\right)$, and $s_{i+\frac{1}{2},j+\frac{1}{2}}^{\ast,1}$ is the speed of contact discontinuity in the HLLC solver. And the computation of \mbox{\boldmath\smallu}^{(1)}_{i+\frac{1}{2},j+\frac{1}{2}} can be done by the similar way.

4 Numerical results

This section conducts some numerical experiments on several ultra-relativistic RHD problems with large Lorentz factor, or strong discontinuities, or low rest-mass density or pressure etc. to verify the accuracy, robustness and effectiveness of the studied PCP Lagrangian schemes. It is worth stressing that those ultra-relativistic RHD problems seriously challenge the numerical schemes. Unless otherwise stated, all the computations are restricted to the equation of state (1.5) with the adiabatic index $\Gamma=1.4$, and the time step size $\Delta t$ of the 1D (resp. 2D) first-order schemes is determined by (2.24) (resp. (3.14)), while it is decided by (2.30) (resp. (3.24)) for the 1D (resp. 2D) high order schemes, where the value of the parameter $\lambda$ is uniformly taken as $1/2$.

4.1 1D case

Example 4.1** (Accuracy test).**

It is to test the accuracy and PCP property of our schemes. The initial density and pressure are given by

[TABLE]

where $\rho_{\text{ref}}=10^{-7},K=10^{-1},L=0.3,\alpha=1$, and $\Gamma=5/3$. The initial velocity $u$ is specified by assuming that the Riemann invariant

[TABLE]

is constant. It describes an isentropic pulse moving in a smooth domain, similar to one in [50]. The computational domain is taken as $[-0.35,1]$, and the exact solution can be obtained by the method of characteristics.

The errors $\varepsilon_{1}$, $\varepsilon_{2}$ and $\varepsilon_{\infty}$ and corresponding orders of convergence at $t=0.02$ obtained by the first- and third-order Lagrangian schemes are shown in Tables 4.1 and 4.2, where the errors are defined by

[TABLE]

here $\mathbf{U}$ and $\mathbf{U}_{h}$ are the exact and numerical solutions at $t$, respectively. Table 4.2 also lists the proportions of the PCP limited cells at all time levels, denoted by $\Theta_{N}$. It is shown that the PCP limiter has been performed in the higher-order accurate schemes because of the low pressure, and the higher-order PCP Lagrangian schemes can achieve the theoretical accuracy.

Example 4.2** (Blast wave interaction).**

This is an initial-boundary-value problem for the 1D RHD equations and has been studied in [20, 39, 47]. The same initial setup is considered here. The computational domain is $[0,1]$ with outflow boundary conditions, and the initial condition is

[TABLE]

It is a severe test because the interaction is happened in a very narrow region and there meet the low density and pressure, and large velocity in the domain. If the PCP limiter is not employed, then the calculation of the second- or third-order schemes result in failure as soon as the computed pressure or density becomes negative. Figure 4.1 shows the numerical results at $t=0.43$ obtained by using the second- and third-order schemes. It can be found that the solutions within the interval $[0.5,0.53]$ obtained by using the third-order scheme are in good agreement with the exact solution, and the discontinuities are exactly and well captured on a coarse mesh with $400$ cells.

4.2 2D case

Example 4.3** (Accuracy test).**

Similar to [32], a 2D relativistic isentropic vortex problem is constructed here to test the accuracy of our Lagrangian schemes. First, in a coordinate system $S$ with the spacetime coordinates $(t,x,y)$, a steady, relativistic isentropic vortex solution $(\rho,u_{x},u_{y},p)$ of the 2D RHD equations is obtained as follows

[TABLE]

where the vortex strength is $\epsilon=10.0828$ such that the lowest density and lowest pressure are $7.8\times 10^{-15}$ and $1.78\times 10^{-20}$, respectively. Next, assume that a coordinate system $S^{\prime}$ with the spacetime coordinates $(t^{\prime},x^{\prime},y^{\prime})$ is in motion relative to the coordinate system $S$ with a constant velocity of magnitude $w$ along the $(1,1)$ direction, from the perspective of an observer stationary in $S$. Then the relation between the two coordinate systems is given by the Lorentz transformation as

[TABLE]

and the transformation between the velocities is

[TABLE]

Using those transformations can give a time-dependent solution $(\rho^{\prime},u^{\prime},v^{\prime},p^{\prime})$ in the coordinate system $S^{\prime}$

[TABLE]

The vortex in the coordinate system $S^{\prime}$ moves with a constant speed of magnitude $w$ in $(-1,-1)$ direction. Unlike the non-relativistic case, the relativistic circular vortex is contracted in $(1,1)$ direction due to the Lorentz contraction, thus it becomes elliptic in the coordinate system $S^{\prime}$.

Our numerical simulation is performed in an initial square $\Omega(0)=[-5,5]^{2}$ with $w=0.5$ and periodic boundary conditions, and the output time is $t=1$. For such problem, the second order Lagrangian scheme without the PCP limiter will break down because the numerical solution cannot be guaranteed to be in the admissible state set during the computation. Tables 4.3 and 4.4 list the errors $\varepsilon_{1}$, $\varepsilon_{2}$ and $\varepsilon_{\infty}$ and orders of convergence. Moreover, Table 4.4 also gives the proportions of the PCP limited cells at all time levels, denoted by $\Theta_{N}$. Figure 4.2 plots the deformed mesh and shifted elliptic vortex in the primitive variables at $t=1$ obtained by the second-order PCP Lagrangian scheme. It is shown that the first- and second-order PCP Lagrangian schemes can achieve the expected theoretical accuracy, and the PCP limiter has been performed in the second-order accurate schemes.

Example 4.4** (Blast problem on the Cartesian mesh).**

Consider a blast problem in a square domain $[0,1]^{2}$ with reflective boundary conditions. The initial data are specified as follows

[TABLE]

with $r=\sqrt{x^{2}+y^{2}}$. The domain is initially divided into a uniform Cartesian mesh with $60\times 60$ cells. If the second-order Lagrangian scheme is not PCP, then the negative density may be numerically obtained so that the calculation will result in failure.

Figure 4.3 plots the mesh and the density contour at $t=0.4$ obtained by using the second-order PCP Lagrangian scheme, and the density and pressure along the line $y=x$ obtained by the first- and second-order PCP Lagrangian schemes. The solid line in Figure 4.3(c)-(d) is the reference solution obtained by using a second order TVD Eulerian scheme with Lax-Friedrichs flux in the cylindrical coordinate and with $10000$ cells, and clearly shows that the solution consists of a “left-moving” rarefaction wave, a “right-moving” contact discontinuity and a “right-moving” shock wave. It is seen that our PCP Lagrangian schemes capture the contact discontinuity and the shock wave with high resolution, and the second-order scheme gives a better result than the first-order.

Example 4.5** (Blast problem on the polar mesh).**

It is to solve the problem in Example 4.4 on an equal-angled polar mesh with $100\times 10$ cells. Figure 4.4(a)-(b) shows the mesh and density at $t=0.4$ obtained by using the second-order PCP Lagrangian scheme, while Figure 4.4(c)-(d) plots the density and pressure with respect to the radial radius obtained by the first- and second-order schemes. The results show that the contact discontinuity can be exactly resolved and the second-order PCP Lagrangian scheme gives a better result than the first-order. Moreover, the results in Figure 4.4 are obviously better than those in Figure 4.3.

Example 4.6** (Strong blast problem on the polar mesh).**

Consider a strong blast problem with the initial data

[TABLE]

It is similar to Example 4.5 so that the flow pattern is similar to that in Example 4.5. However, it is very challenging for the Lagrangian scheme because there exists the extremely low pressure and the very narrow region between the contact discontinuity and the shock wave. If the scheme is not PCP, then its calculation may result in failure.

Figure 4.5(a)-(c) the pressure, density and the close-up of the density obtained by using the first and the second order schemes with $100\times 10$ cells, where the solid line is the reference solution obtained by a second-order TVD Eulerian scheme with Lax-Friedrichs flux in the cylindrical coordinate with $10000$ cells. Figure 4.5(d) shows the mesh with $100\times 10$ cells at $t=0.4$ obtained by using the second-order PCP Lagrangian scheme. It can be seen that our PCP schemes can capture the narrow region between the contact discontinuity and the shock wave well, and the second order scheme is better than the first-order.

Example 4.7** (Implosion problem).**

It is an implosion problem, similar to the Noh problem in the non-relativistic case. The computational domain in the polar coordinates $(r,\theta)$ is chosen as $[0,1]\times[0,\pi/2]$ and divided into an initial equal-angled polar mesh. The initial density and pressure of the fluid in the domain is $1$ and $10^{-12}$ respectively, and an inward radial velocity is set as $0.9$. A strong shock wave is generated by bringing the cold gas to rest at the origin, and will converge to the origin. Because the low pressure will appear in the solution, the negative pressure may be easily produced by using the non-PCP scheme.

Figure 4.6(a)-(b) shows the mesh with $100\times 10$ cells and the density contour at $t=0.6$ obtained by the second-order PCP Lagrangian scheme, while Figure 4.6(c)-(d) plots the radial density and pressure obtained by using the first- and second-order schemes, where the solid line is the reference solution obtained by using a first-order Eulerian scheme with Lax-Friedrichs flux in the cylindrical coordinate and with $10000$ cells. It is shown that the PCP scheme can be successfully simulate such extreme relativistic fluid flow, and the resolution of the second-order scheme is obviously better than the first-order even though there exists a wall heating phenomenon.

Example 4.8** (ICF-like test).**

The last problem is about a ICF-like test, which is similar to Example 6.7 in [18]. The initial target is defined in the polar coordinates $(r,\theta)$ by $[0,1]\times[0,\pi/2]$. There are two kinds of materials in the target with the same adiabatic index $\Gamma=5/3$. The target is divided into the internal part whose radius is $0.9$ and the external shell. The initial condition is specified as follows

[TABLE]

and on the external boundary of the shell, an additional pressure is added as

[TABLE]

The period of the perturbation of $\widetilde{p}$ is $\pi/3$.

Figures 4.7-4.8 show the mesh, density $\rho$, radial velocity $u_{r}$, and the pressure $p$ at $t=1$ obtained by using the first- and second-order PCP Lagrangian schemes with $100\times 30$ cells, respectively. We can see that the outmost mesh becomes sunken periodically due to the perturbation in the additional pressure, and the mesh points gather near the interface between the internal part and the external shell. Obviously, the second-order PCP Lagrangian scheme gives a better resolution near the interface than the first-order.

5 Conclusions

This paper studied the physical-constraints-preserving (PCP) Lagrangian finite volume schemes for one- and two-dimensional special relativistic hydrodynamic (RHD) equations. Our attention was paid to the Lagrangian schemes with the HLLC Riemann solver. First, we proved that the intermediate states in the HLLC Riemann solver were admissible or PCP (that is, the rest-mass density and pressure are positive and the fluid velocity magnitude is less than the speed of light) when the HLLC wave speeds were estimated suitably. It was worth noting that such PCP property has been observed in [25], but no rigorously mathematical proof was given there since it was confront with considerable difficulty and pretty challenging. Then we showed that the first-order accurate Lagrangian scheme with the HLLC Riemann solver and forward Euler time discretization was PCP and developed the higher-order accurate PCP Lagrangian schemes by using the high-order accurate strong stability preserving (SSP) time discretizations, the WENO reconstruction procedure and the scaling PCP limiter. Finally, several one- and two-dimensional numerical experiments were conducted to demonstrate the accuracy and the effectiveness of the PCP Lagrangian schemes in solving the special RHD problems involving large Lorentz factor, or low rest-mass density or low pressure or strong discontinuities etc.

Acknowledgements

This work was partially supported by the Special Project on High-performance Computing under the National Key R&D Program (No. 2016YFB0200603), Science Challenge Project (No. JCKY2016212A502), and the National Natural Science Foundation of China (Nos. 91330205, 91630310, 11421101, 11871111).

Appendix A Proof of Lemma 2.1

This appendix proves the properties (i)-(vii) presented in Lemma 2.1 orderly. For the sake of simplicity, denote s_{\min}^{\pm}=s_{\min}(\mbox{\boldmath\smallU}^{\pm}) and s_{\max}^{\pm}=s_{\max}(\mbox{\boldmath\smallU}^{\pm}).

(i) Since $A^{\pm}$ defined in (2.11) can be rewritten as

[TABLE]

where the superscript $\pm$ has been omitted, for left and right states state \mbox{\boldmath\smallU}^{\mp} one has

[TABLE]

and

[TABLE]

(ii) Due to the definition of $A^{\pm}$ and $B^{\pm}$ in (2.11), it can have

[TABLE]

and

[TABLE]

(iii) For the states \mbox{\boldmath\smallU}^{\pm}, it is easy to verify

[TABLE]

(iv) Because $A^{-}<0$, one has

[TABLE]

Similarly, the fact $A^{+}>0$ gives

[TABLE]

(v) Because

[TABLE]

for the left state \mbox{\boldmath\smallU}^{-} and

[TABLE]

for the right state \mbox{\boldmath\smallU}^{+}, it holds

[TABLE]

The remaindering is to prove show the inequality

[TABLE]

Because (2.13) gives

[TABLE]

and the wave speed is less than the speed of light $c=1$, two terms $s^{\ast}A^{-}-B^{-}$ and $s^{\ast}A^{+}-B^{+}$ should have the same sign. It means that

[TABLE]

which gives

[TABLE]

On the other hand, the properties (i) and (iii) give

[TABLE]

and the properties (i) and (iv) lead to

[TABLE]

so one has

[TABLE]

Combing (A.2) with (A.3) complete the proof of (A.1).

(vi) The left hand side of the inequality (vi) can be recast into

[TABLE]

with

[TABLE]

Thus, one has to prove $f_{1}^{\pm}f_{2}^{\pm}<0$ in order to draw the conclusion (vi).

For the left state \mbox{\boldmath\smallU}^{-}, because $s^{-}\in(-1,s_{\min}^{-}]$ with $s_{\min}^{-}=\frac{u^{-}-c_{s}^{-}}{1-c_{s}^{-}u^{-}}<1$, one has

[TABLE]

and

[TABLE]

where

[TABLE]

Similarly, for the right state \mbox{\boldmath\smallU}^{+}, the fact that $s^{+}\in[s_{\max}^{+},1)$ with $s_{\max}^{+}=\frac{u^{+}+c_{s}^{+}}{1+c_{s}^{+}u^{+}}>-1$ means

[TABLE]

and

[TABLE]

where

[TABLE]

In a word, the conclusion (vi) holds.

(vii) The left hand side of the inequality (vii) corresponds to the value of the following quadratic function $f^{\pm}(s)$ at $s=s^{\pm}$

[TABLE]

Because

[TABLE]

the parabolas $f^{\pm}(s)$ are convex and symmetric with the lines

[TABLE]

For the left state \mbox{\boldmath\smallU}^{-}, one has

[TABLE]

with

[TABLE]

Hence $f^{-}(s)$ is monotonically decreasing with $s\in(-1,s_{\min}^{-}]$. It means that $f^{-}(s)\geq f^{-}(s_{\min}^{-})$ for any $s\in(-1,s_{\min}^{-}]$ and thus $f^{-}(s^{-})\geq f^{-}(s_{\min}^{-})$ since $s^{-}=\min(s_{\min}^{-},s_{\min}^{+})$. To get the inequality (vii) for the left state \mbox{\boldmath\smallU}^{-}, the remaining is to prove $f^{-}(s_{\min}^{-})>0$. In fact, it is true because

[TABLE]

Similarly, for the right state \mbox{\boldmath\smallU}^{+}, $f^{+}(s)$ is monotonically increasing with $s\in[s_{\max}^{+},1)$ and then $f^{+}(s)\geq f^{+}(s_{\max}^{+})$ for any $s\in[s_{\max}^{+},1)$ since

[TABLE]

with

[TABLE]

To get the inequality (vii) for the left state \mbox{\boldmath\smallU}^{+}, the remaining is to prove $f^{+}(s_{\max}^{+})>0$. It is true since

[TABLE]

It completes the proof.

Appendix B Proof of Lemma 3.1

This appendix discusses the proof of Lemma 3.1. Here only the properties (vi) and (vii) are proved in detail, since the proof of the properties (i)-(v) is similar to that of Lemma 2.1. For the sake of simplicity, denote s_{\min}^{\pm}=s_{\min}(\mbox{\boldmath\smallU}^{\pm}), s_{\max}^{\pm}=s_{\max}(\mbox{\boldmath\smallU}^{\pm}), and define

[TABLE]

Due to (3.4), $s_{\min}^{\pm}$ and $s_{\max}^{\pm}$ can be expressed as

[TABLE]

Moreover, it can be verified that

[TABLE]

(vi) The left hand side of (vi) can be decomposed into the two parts

[TABLE]

where

[TABLE]

Thus the remaining is to prove that $f_{1}^{\pm}f_{2}^{\pm}<0$.

For the left state \mbox{\boldmath\smallU}^{-}, because $s^{-}\in(-1,s_{\min}^{-}]$ with s_{\min}^{-}=\frac{u_{n}^{-}(1-(c_{s}^{-})^{2})-c_{s}^{-}\alpha^{-}}{1-(c_{s}^{-})^{2}|\mbox{\boldmath\smallu}^{-}|^{2}}<u_{n}^{-}, one has

[TABLE]

and

[TABLE]

where

[TABLE]

Similarly, for the right state \mbox{\boldmath\smallU}^{+}, because $s^{+}\in[s_{\max}^{+},1)$ with s_{\max}^{+}=\frac{u_{n}^{+}\big{(}1-(c_{s}^{+})^{2}\big{)}+c_{s}^{+}\alpha^{+}}{1-(c_{s}^{+})^{2}|\mbox{\boldmath\smallu}^{+}|^{2}}>u_{n}^{+}, one has

[TABLE]

and

[TABLE]

where

[TABLE]

Therefore, the inequality (vi) is proved.

(vii) Let us consider the quadratic functions in terms of $s$

[TABLE]

which are symmetric with the lines

[TABLE]

Moreover, $f^{\pm}(s)$ are convex because

[TABLE]

and

[TABLE]

For the left state \mbox{\boldmath\smallU}^{-}, one has

[TABLE]

which implies that $f^{-}(s)$ is monotonically decreasing with $s\in(-1,s_{\min}^{-}]$, where

[TABLE]

Thus $f^{-}(s)\geq f^{-}(s_{\min}^{-})$ for $s\in(-1,s_{\min}^{-}]$, so that $f^{-}(s^{-})\geq f^{-}(s_{\min}^{-})$. The remaining is to prove $f^{-}(s_{\min}^{-})>0$. In fact, one has

[TABLE]

where

[TABLE]

and

[TABLE]

Hence it is true that $f^{-}(s_{\min}^{-})>0$ and thus the inequality (vii) for \mbox{\boldmath\smallU}^{-} is proved.

For the right state \mbox{\boldmath\smallU}^{+}, since

[TABLE]

with

[TABLE]

$f^{+}(s)$ is monotonically increasing with $s\in[s^{+},1)\subset[s_{\max}^{+},1)$ and then we have $f^{+}(s^{+})>f^{+}(s_{\max}^{+})$. Hence to prove the inequality (vii) for \mbox{\boldmath\smallU}^{+}, it suffices to prove $f^{+}(s_{\max}^{+})>0$. In fact, one has

[TABLE]

where

[TABLE]

and

[TABLE]

Hence it is true that $f^{+}(s_{\max}^{+})>0$ and thus the inequality (vii) for \mbox{\boldmath\smallU}^{+} holds. The proof is completed.