Second-order accurate genuine BGK schemes for the ultra-relativistic flow simulations

TL;DR

This paper introduces second-order accurate genuine BGK schemes based on the analytical solution of the Anderson-Witting model for simulating ultra-relativistic flows, including viscous effects, with improved accuracy and resolution.

Contribution

It presents the first derivation of genuine BGK schemes from the Anderson-Witting model for ultra-relativistic flows, incorporating particle collisions explicitly.

Findings

Accurate and stable for ultra-relativistic inviscid and viscous flows.

Higher resolution at contact discontinuities compared to existing schemes.

Validated through multiple 1D and 2D numerical experiments.

Abstract

This paper presents second-order accurate genuine BGK (Bhatnagar-Gross-Krook) schemes in the framework of finite volume method for the ultra-relativistic flows. Different from the existing kinetic flux-vector splitting (KFVS) or BGK-type schemes for the ultra-relativistic Euler equations, the present genuine BGK schemes are derived from the analytical solution of the Anderson-Witting model, which is given for the first time and includes the "genuine" particle collisions in the gas transport process. The BGK schemes for the ultra-relativistic viscous flows are also developed and two examples of ultra-relativistic viscous flow are designed. Several 1D and 2D numerical experiments are conducted to demonstrate that the proposed BGK schemes not only are accurate and stable in simulating ultra-relativistic inviscid and viscous flows, but also have higher resolution at the contact discontinuity than the KFVS or BGK-type schemes.

Tables

Table 1. Table 5.1: Example 5.1 : Numerical errors of n 𝑛 n in l 1 , l 2 superscript 𝑙 1 superscript 𝑙 2 l^{1},l^{2} -norms and convergence rates at t = 0.2 𝑡 0.2 t=0.2 with or without limiter.

$N$	With limiter				Without limiter
	$l^1$ error	$l^1$ order	$l^2$ error	$l^2$ order	$l^1$ error	$l^1$ order	$l^2$ error	$l^2$ order
25	1.6793e-03	–	2.5667e-03	–	6.0337e-04	-	6.7007e-04	-
50	4.9516e-04	1.7619	8.2151e-04	1.6436	1.5275e-04	1.9819	1.6965e-04	1.9818
100	1.3012e-04	1.9281	2.6823e-04	1.6148	3.8305e-05	1.9956	4.2559e-05	1.9950
200	3.4917e-05	1.8978	8.5622e-05	1.6474	9.5628e-06	2.0020	1.0621e-05	2.0025
400	8.2820e-06	2.0759	2.6141e-05	1.7117	2.3904e-06	2.0002	2.6550e-06	2.0001

Table 2. Table 5.2: Example 5.7 : Numerical errors at t = 0.1 𝑡 0.1 t=0.1 in l 1 , l 2 superscript 𝑙 1 superscript 𝑙 2 l^{1},l^{2} -norms and convergence rates with or without limiter.

N	With limiter				Without limiter
	$l^1$ error	$l^1$ order	$l^2$ error	$l^2$ order	$l^1$ error	$l^1$ order	$l^2$ error	$l^2$ order
25	1.7316e-03	-	2.5820e-03	-	6.1369e-04	-	6.8214e-04	-
50	5.3784e-04	1.6869	9.1457e-04	1.4974	1.5610e-04	1.9751	1.7340e-04	1.9759
100	1.4248e-04	1.9164	2.8992e-04	1.6574	3.9584e-05	1.9795	4.3962e-05	1.9798
200	3.8119e-05	1.9022	9.5759e-05	1.5982	9.8942e-06	2.0003	1.0989e-05	2.0002
400	1.0923e-05	1.8031	3.1779e-05	1.5914	2.4837e-06	1.9941	2.7586e-06	1.9941

Table 3. Table 5.3: Example 5.11 : Numerical errors of n 𝑛 n in l 1 superscript 𝑙 1 l^{1} and l 2 superscript 𝑙 2 l^{2} -norm and convergence rates at t = 1.2 𝑡 1.2 t=1.2 .

$N$	$l^1$ error	$l^1$ order	$l^2$ error	$l^2$ order
10	8.9214e-03	–	1.2028e-02	–
20	1.9291e-03	2.2094	2.4837e-03	2.2759
40	4.9766e-04	1.9546	6.2325e-04	1.9946
80	1.3682e-04	1.8629	1.6760e-04	1.8948

Equations

$$p_{\alpha }=g_{\alpha \beta }{p}^{\beta },p^{\alpha }=g^{\alpha \beta }{p}_{\beta },$$

$$p^{\alpha }\frac{\partial f}{\partial x^{\alpha }}=Q(f,f),$$

$$p^{\alpha }\frac{\partial f}{\partial x^{\alpha }}=-\frac{U_{\alpha }{p}^{\alpha }}{\tau c^2}(f-g),$$

$$U_{\beta }{T}^{\alpha \beta }=\epsilon g^{\alpha \beta }{U}_{\alpha },$$

$$g={\frac{nc^{3}}{8\pi k^{3}T^{3}}}\exp\left(-\frac{U_{\alpha}p^{\alpha}}{{kT}}\right)={\frac{nc^{3}}{8\pi k^{3}T^{3}}}\exp\left(-\frac{|\mbox{\boldmath\small$p$}|}{{kT}}\left(U_{0}-\sum_{i=1}^{3}U_{i}\frac{p^{i}}{|\mbox{\boldmath\small$p$}|}\right)\right),$$

$$\int_{\mathbb{R}^{3}}\frac{U_{\alpha}p^{\alpha}}{\tau c^{2}}(f-g)\Psi\frac{d^{3}\mbox{\boldmath\small$p$}}{p^{0}}=0,\ \ \mbox{\boldmath\small$\Psi$}=(1,p^{i},p^{0})^{T},$$

$${\partial }_{\alpha }{N}^{\alpha }=0,{\partial }_{\beta }{T}^{\alpha \beta }=0,$$

$$\displaystyle N^{\alpha}=c\int_{\mathbb{R}^{3}}p^{\alpha}f\frac{d^{3}\mbox{\boldmath\small$p$}}{p^{0}},\ \ T^{\alpha\beta}=c\int_{\mathbb{R}^{3}}p^{\alpha}p^{\beta}f\frac{d^{3}\mbox{\boldmath\small$p$}}{p^{0}}.$$

$$N^{\alpha }$$

$$T^{\alpha \beta }$$

$${\Delta }^{\alpha \beta }=g^{\alpha \beta }-\frac{1}{c^2}{U}^{\alpha }{U}^{\beta },$$

$$n=$$

$$n^{\alpha }=$$

$$\epsilon =$$

$${\pi }^{\alpha \beta }=$$

$$p+\varPi=-\frac{1}{3}\Delta_{\alpha\beta}T^{\alpha\beta}=\frac{1}{3c}\int_{\mathbb{R}^{3}}(E^{2}-m^{2}c^{4})f\frac{d^{3}\mbox{\boldmath\small$p$}}{p^{0}}.$$

$${\Delta }_{\mu \nu }^{\alpha \beta }=\frac{1}{2}({\Delta }_{\mu }^{\alpha }{\Delta }_{\nu }^{\beta }+{\Delta }_{\mu }^{\beta }{\Delta }_{\nu }^{\alpha }-\frac{2}{3}{\Delta }_{\mu \nu }{\Delta }^{\alpha \beta }).$$

$$(U_{\alpha})=(\gamma c,-\gamma\mbox{\boldmath\small$u$}),\ \ (U^{\alpha})=(\gamma c,\gamma\mbox{\boldmath\small$u$}),$$

$$n=c^{-2}{U}_{\alpha }{N}^{\alpha }>0.$$

$$\epsilon =nm{c}^2(G(\zeta )-{\zeta }^{-1}),$$

$$K_{\nu }(\zeta ):={\int }_0^{\infty }\cosh (\nu \vartheta )\exp (-\zeta \cosh \vartheta )d\vartheta ,\nu \ge 0.$$

$$\epsilon =3knT.$$

$$\displaystyle\mbox{\boldmath\small$X$}^{T}T^{\alpha\beta}\mbox{\boldmath\small$X$}$$

$$\displaystyle=c\int_{\mathbb{R}^{3}}(x_{\alpha}p^{\alpha})^{2}f\frac{d^{3}\mbox{\boldmath\small$p$}}{p^{0}}>0,$$

$$0<U_{\alpha }{T}^{\alpha \beta }{U}_{\beta }=\epsilon U_{\alpha }{g}^{\alpha \beta }{U}_{\beta },$$

$$|\mbox{\boldmath\small$u$}|=cU^{-1}_{0}\sqrt{U^{2}_{1}+U^{2}_{2}+U^{2}_{3}}<c.$$

$$\displaystyle\mbox{\boldmath\small$U$}\cdot\mbox{\boldmath\small$p$}$$

$$n=\frac{1}{c^{2}}U_{\alpha}N^{\alpha}=\frac{1}{c}\int_{\mathbb{R}^{3}}Ef\frac{d^{3}\mbox{\boldmath\small$p$}}{p^{0}}>0.$$

$$n$$

$$\epsilon$$

$$p$$

$$\frac{\partial\mbox{\boldmath\small$W$}}{\partial t}+\sum^{3}_{k=1}\frac{\partial\mbox{\boldmath\small$F$}^{k}(\mbox{\boldmath\small$W$})}{\partial x^{k}}=0,$$

$$\mbox{\boldmath\small$W$}=\left(N^{0},T^{0i},T^{00}\right)^{T}=\left(nU^{0},nhU^{0}U^{i},nhU^{0}U^{0}-p\right)^{T},$$

$$\mbox{\boldmath\small$F^{k}(W)$}=\left(N^{k},T^{ki},T^{k0}\right)^{T}=\left(nU^{k},nhU^{k}U^{i}+p\delta^{ik},nhU^{k}U^{0}\right)^{T}.$$

Full text

1 Introduction

Relativistic hydrodynamics (RHD) arise in astrophysics, nuclear physics, plasma physics and other fields. In many radiation hydrodynamics problems of astrophysical interest, the fluid moves at extremely high velocities near the speed of light, and relativistic effects become important. Examples of such flows are supernova explosions, the cosmic expansion, and solar flares.

The relativistic hydrodynamical equations are highly nonlinear, making the analytic treatment of practical problems extremely difficult. The numerical simulation is the primary and powerful way to study and understand the relativistic hydrodynamics. This work will mainly focus on the numerical methods for the special RHDs, where there is no strong gravitational field involved. The pioneering numerical work may date back to the finite difference code via artificial viscosity for the spherically symmetric general RHD equations in the Lagrangian coordinate [30, 31] and the finite difference method with the artificial viscosity technique for the multi-dimensional RHD equations in the Eulerian coordinate [48]. Since 1990s, the numerical study of the RHDs began to attract considerable attention, and various modern shock-capturing methods with an exact or approximate Riemann solver have been developed for the RHD equations. Some examples are the local characteristic approach [25], the two-shock approximation solvers [5, 8], the Roe solver [13], the flux corrected transport method [12], the flux-splitting method based on the spectral decomposition [11], the piecewise parabolic method [26, 33], the HLL (Harten-Lax-van Leer) method [42], the HLLC (Harten-Lax-van Leer-Contact) method [32] and the Steger-Warming flux vector splitting method [59]. The analytical solution of the Riemann problem in relativistic hydrodynamics was studied in [28]. Some other higher-order accurate methods have also been well studied in the literature, e.g. the ENO (essentially non-oscillatory) and weighted ENO methods [10, 9, 47], the discontinuous Galerkin (DG) method [40], the adaptive moving mesh methods [15, 16], the Runge-Kutta DG methods with WENO limiter [60, 61, 62], the direct Eulerian GRP schemes [56, 57, 52], and the local evolution Galerkin method [49]. Recently some physical-constraints-preserving (PCP) schemes were developed for the special RHD equations. They are the high-order accurate PCP finite difference weighted essentially non-oscillatory (WENO) schemes and discontinuous Galerkin (DG) methods proposed in [50, 51, 39]. The readers are also referred to the early review articles [27, 14] as well as references therein.

The gas-kinetic schemes present a gas evolution process from a kinetic scale to a hydrodynamic scale, where both inviscid and viscous fluxes are recovered from moments of a single time-dependent gas distribution function [34]. The development of gas-kinetic schemes, such as the kinetic flux vector splitting (KFVS) and Bhatnagar-Gross-Krook (BGK) schemes, has attracted much attention and significant progress has been made in the non-relativistic hydrodynamics. They utilize the well-known connection that the macroscopic governing equations are the moments of the Boltzmann equation whenever the distribution function is at equilibrium. The KFVS schemes are constructed by applying upwind technique directly to the collisionless Boltzmann equation, see e.g. [36, 24, 7, 17, 41, 35, 46, 45, 44]. Due to the lack of collision in the numerical flux calculations, the KFVS schemes smear the solutions, especially the contact discontinuity. To overcome this problem, the BGK schemes are constructed by taking into account the particle collisions in the whole gas evolution process within a time step, see e.g. [22, 54, 23]. Moreover, due to their specific derivation, they are also able to present the accurate Navier-Stokes solution in the smooth flow regime and have favorable shock capturing capability in the shock region. The kinetic beam scheme was first proposed for the relativistic gas dynamics in [55]. After that, the kinetic schemes for the ultra-relativistic Euler equations were developed in [19, 20, 21]. The BGK-type schemes [53, 46] were extended to the ultra-relativistic Euler equations in [18, 38] in order to reduce the numerical dissipation. Those kinetic schemes resulted directly from the moments of the relativistic J$\ddot{\text{u}}$ttner equilibrium distribution without including the “genuine” particle collisions in the gas transport process.

This paper will develop second-order genuine BGK schemes for the ultra-relativistic inviscid and viscous flow simulations. It is organized as follows. Section 2 introduces the special relativistic Boltzmann equation and discusses how to recover some macroscopic quantities from the kinetic theory. Section 3 presents the ultra-relativistic hydrodynamical equations through the Chapman-Enskog expansion. Section 4 develops second-order accurate genuine BGK schemes for the 1D and 2D ultra-relativistic Euler equations and 2D ultra-relativistic Navier-Stokes equations. Section 5 gives several numerical experiments to demonstrate accuracy, robustness and effectiveness of the proposed schemes in simulating inviscid and viscous ultra-relativistic fluid flows. Section 6 concludes the paper.

2 Preliminaries and notations

In the special relativistic kinetic theory of gases [6], a microscopic gas particle is characterized by the four-dimensional space-time coordinates (x^{\alpha})=(x^{0},\mbox{\boldmath\smallx}) and four-momentum vectors (p^{\alpha})=(p^{0},\mbox{\boldmath\smallp}), where $x^{0}=ct$, $c$ denotes the speed of light in vacuum, $t$ and $x$ are the time and 3D spatial coordinates, respectively, and the Greek index $\alpha$ runs from $0,1,2,3$. Besides the contravariant notation (e.g. $p^{\alpha}$), the covariant notation such as $p_{\alpha}$ will also be used in the following, while both notations $p^{\alpha}$ and $p_{\alpha}$ are related by

[TABLE]

where the Einstein summation convention over repeated indices has been used, $(g^{\alpha\beta})$ is the Minkowski space-time metric tensor and chosen as $(g^{\alpha\beta})=\text{diag}\{1,-1,-1,-1\}$, while $(g_{\alpha\beta})$ denotes the inverse of $(g^{\alpha\beta})$.

For a free relativistic particle, the relativistic energy-momentum relation (aka “on-shell” or “mass-shell” condition) E^{2}-|\mbox{\boldmath\smallp}|^{2}c^{2}=m^{2}c^{4} holds, where $m$ denotes the mass of each structure-less particle which is assumed to be the same for all particles. The “mass-shell” condition can be rewritten as $p^{\alpha}p_{\alpha}=m^{2}c^{2}$ if putting p^{0}=c^{-1}E=\sqrt{|\mbox{\boldmath\smallp}|^{2}+m^{2}c^{2}}, which becomes p^{0}=|\mbox{\boldmath\smallp}| in the ultra-relativistic limit, i.e. $m\to 0$.

Similar to the non-relativistic case, the relativistic Boltzmann equation describes the evolution of one-particle distribution function f(\mbox{\boldmath\smallx},t,\mbox{\boldmath\smallp}) in the phase space spanned by the space-time coordinates $x^{\alpha}$ and momentum $p^{\alpha}$ of particles. It reads

[TABLE]

where $Q(f,f)$ denotes the collision term and depends on the product of distribution functions of two particles at collision. In the literature, there exist several simple collision models. The Anderson-Witting model [4]

[TABLE]

is similar to the BGK model in the non-relativistic kinetic theory and will be considered in this paper, where $\tau$ is the relaxation time, in the Landau-Lifshitz frame, the hydrodynamic four-velocities $U_{\alpha}$ are defined by

[TABLE]

which implies that $(\varepsilon,U_{\alpha})$ is a generalized characteristic pair of $(T^{\alpha\beta},g^{\alpha\beta})$, $\varepsilon$ and $T^{\alpha\beta}$ are the energy density and energy-momentum tensor, respectively, and g=g(\mbox{\boldmath\smallx},t,\mbox{\boldmath\smallp}) denotes the distribution function at the local thermodynamic equilibrium, the so-called J$\ddot{\text{u}}$ttner equilibrium (or relativistic Maxwellian) distribution. In the ultra-relativistic case, it becomes [19]

[TABLE]

where $n$ and $T$ denote the number density and thermodynamic temperature, respectively, and $k$ is the Boltzmann’s constant. The Anderson-Witting model (2.2) can tend to the BGK model in the non-relativistic limit and the collision term $-\frac{U_{\alpha}p^{\alpha}}{\tau c^{2}}(f-g)$ satisfies the following identities

[TABLE]

which imply the conservation of particle number, momentum and energy

[TABLE]

where the particle four-flow $N^{\alpha}$ and the energy-momentum tensor $T^{\alpha\beta}$ are related to the distribution $f$ by

[TABLE]

In the Landau-Lifshitz decomposition, both $N^{\alpha}$ and $T^{\alpha\beta}$ are rewritten as follows

[TABLE]

where $\Delta^{\alpha\beta}$ is defined by

[TABLE]

satisfying $\Delta^{\alpha\beta}U_{\beta}=0$, the number density $n$, particle-diffusion current $n^{\alpha}$, energy density $\varepsilon$, and shear-stress tensor $\pi^{\alpha\beta}$ can be calculated by

[TABLE]

and the sum of thermodynamic pressure $p$ and bulk viscous pressure $\varPi$ is

[TABLE]

Here $E=U_{\alpha}p^{\alpha}$, $p^{<\alpha>}=\Delta^{{\alpha}}_{\gamma}p^{\gamma}$, $p^{<\alpha\beta>}=\Delta^{\alpha\beta}_{\gamma\delta}p^{\gamma}p^{{\delta}}$, and

[TABLE]

Remark 2.1

The quantities $n^{\alpha},\varPi$, and $\pi^{\alpha\beta}$ become zero at the local thermodynamic equilibrium $f=g$.

The following gives a general recovery procedure of the admissible primitive variables $n$, $u$, and $T$ from the nonnegative distribution f(\mbox{\boldmath\smallx},t,\mbox{\boldmath\smallp}), where $u$ is the macroscopic velocity in the $(x^{i})$ space. Such recovery procedure will be useful in our BGK scheme.

Theorem 2.1

For any nonnegative distribution f(\mbox{\boldmath\smallx},t,\mbox{\boldmath\smallp}) which is not always be zero, the number density $n$, velocity $u$ and temperature $T$ can be uniquely obtained as follows:

$T^{\alpha\beta}$* is positive definite and $(T^{\alpha\beta},g^{\alpha\beta})$ has only one positive generalized eigenvalue, i.e. the energy density $\varepsilon$, and $U_{\alpha}$ is corresponding generalized eigenvector satisfying $U_{0}=\sqrt{U^{2}_{1}+U^{2}_{2}+U^{2}_{3}+c^{2}}$. Thus, the macroscopic velocity $u$ can be calculated by \mbox{\boldmath\smallu}=-c(U^{-1}_{0}U_{1},U^{-1}_{0}U_{2},U^{-1}_{0}U_{3})^{T}, satisfying |\mbox{\boldmath\smallu}|<c and*

[TABLE]

where \gamma=(1-c^{-2}|\mbox{\boldmath\smallu}|^{2})^{-\frac{1}{2}} denotes the Lorentz factor. 2. 2.

The number density $n$ is calculated by

[TABLE] 3. 3.

The temperature $T$ solves the nonlinear algebraic equation

[TABLE]

where $\zeta=\frac{mc^{2}}{kT}$, $G(\zeta)=\frac{K_{3}(\zeta)}{K_{2}(\zeta)}$, and $K_{{\nu}}(\zeta)$ is modified Bessel function of the second kind, defined by

[TABLE]

In the ultra-relativistic case, $K_{2}(\zeta)$ and $K_{3}(\zeta)$ reduce to $\frac{2}{\zeta^{2}}$ and $\frac{8}{\zeta^{3}}$, respectively, so that one has $G(\zeta)=\frac{4}{\zeta}$, and then

[TABLE]

• Proof

Since the nonnegative distribution f(\mbox{\boldmath\smallx},t,\mbox{\boldmath\smallp}) is not identically zero, using the relation (2.7) gives

[TABLE]

for any nonzero vector \mbox{\boldmath\smallX}=(x_{0},x_{1},x_{2},x_{3})^{T}\in\mathbb{R}^{4}. Thus, the matrix $T^{\alpha\beta}$ is positive definite.

Thanks to $g^{\alpha\beta}=\text{diag}\{1,-1,-1,-1\}$ and (2.3), the matrix-pair $(T^{\alpha\beta},g^{\alpha\beta})$ has an unique positive generalized eigenvalue $\varepsilon$, satisfying

[TABLE]

which implies $U^{2}_{0}>U^{2}_{1}+U^{2}_{2}+U^{2}_{3}$. Thus, one can obtain $U_{0}=\sqrt{U^{2}_{1}+U^{2}_{2}+U^{2}_{3}+c^{2}}$ via multiplying $(U_{\alpha})$ by a scaling constant $c(U^{2}_{0}-U^{2}_{1}-U^{2}_{2}-U^{2}_{3})^{-1/2}$. As a result, the macroscopic velocity $u$ can be calculated by \mbox{\boldmath\smallu}=-c(U^{-1}_{0}U_{1},U^{-1}_{0}U_{2},U^{-1}_{0}U_{3})^{T}, satisfying

[TABLE] 2. 2.

For ${U}_{i}\in\mathbb{R}$ and ${p}^{i}\in\mathbb{R}$, $i=1,2,3$, using the Cauchy-Schwarz inequality gives

[TABLE]

which implies $E=U_{\alpha}p^{\alpha}>0$. Thus one has

[TABLE] 3. 3.

It is obvious that the positive temperature $T$ can be obtained from (2.20).

3 Ultra-relativistic hydrodynamic equations

This section gives the ultra-relativistic hydrodynamic equations, which can be derived from the Anderson-Witting model by using the Chapman-Enskog expansion. For the sake of convenience, units in which the speed of light and the Boltzmann’s constant are equal to one will be used here and hereafter.

3.1 Euler equations

In the ultra-relativistic limit, the macroscopic variables $n,\varepsilon,p$ are related to $g$ by

[TABLE]

If taking the zero order Chapman-Enskog expansion $f=g$ and using the conclusion in Remark 2.1, the ultra-relativistic Euler equations are derived as follows

[TABLE]

where

[TABLE]

and

[TABLE]

Here $i=1,2,3$ and $h=4T$ denotes the specific enthalpy. For the given conservative vector $W$, one can get the primitive variables $n,U^{k}$ and $p$ by [18]

[TABLE]

3.2 Navier-Stokes equations

If taking the first order Chapman-Enskog expansion

[TABLE]

with

[TABLE]

where $\nabla^{\alpha}=\Delta^{\alpha\beta}\partial_{\beta}$ and $\nabla^{<\alpha}U^{\beta>}=\Delta^{\alpha\beta}_{\gamma\delta}\nabla^{\gamma}U^{\delta}$, then (2.12), (2.14) and (2.15) give

[TABLE]

where $\lambda=\frac{4}{3T}p\tau$ and $\mu=\frac{4}{5}p\tau$. Based on those, the ultra-relativistic Navier-Stokes equations see [6] can be obtained as follows

[TABLE]

where

[TABLE]

and

[TABLE]

It shows that one cannot recover the values of primitive variables n,{\mbox{\boldmath\smallu}} and ${T}$ only from the given conservative vector $W$. In practice, the values of n,{\mbox{\boldmath\smallu}} and ${T}$ have to be recovered from the given $W$ and $F^{k}(W)$ or $N^{\alpha}$ and $T^{\alpha\beta}$ by using Theorem 2.1.

4 Numerical schemes

This section develops second-order accurate genuine BGK schemes for the 1D and 2D ultra-relativistic Euler and Navier-Stokes equations. The BGK schemes are derived from the analytical solution of the Anderson-Witting model (2.2), which is given for the first time and includes the “genuine” particle collisions in the gas transport process.

4.1 1D Euler equations

Consider the 1D ultra-relativistic Euler equations with \mbox{\boldmath\smallu}=(u,0,0)^{T} as

[TABLE]

where

[TABLE]

It is strictly hyperbolic because there are three real and distinct eigenvalues of the Jacobian matrix A(\mbox{\boldmath\smallW})=\partial\mbox{\boldmath\smallF}/\partial\mbox{\boldmath\smallW} [50]

[TABLE]

where $c_{s}=1/\sqrt{3}$ is the speed of sound.

Divide the spatial domain into a uniform mesh with the step size $\Delta x$ and the $j$th cell $I_{j}=(x_{j-\frac{1}{2}},x_{j+\frac{1}{2}})$, where $x_{j+\frac{1}{2}}=\frac{1}{2}(x_{j}+x_{j+1})$ and $x_{j}=j\Delta x$, $j\in\mathbb{Z}$. The time interval $[0,T]$ is also divided into a (non-uniform) mesh $\{t_{n+1}=t_{n}+\Delta t_{n},t_{0}=0,n\geqslant 0\}$, where the step size $\Delta t_{n}$ is determined by

[TABLE]

the constant $C$ denotes the CFL number, and $\bar{\varrho}_{j}$ denotes a suitable approximation of the spectral radius of A(\mbox{\boldmath\smallW}) within the cell $I_{j}$. For the given approximate cell-average values \{\mbox{\boldmath\small\bar{W}}^{n}_{j}\}, i.e.

[TABLE]

reconstruct a piecewise linear function as follows

[TABLE]

where \mbox{\boldmath\smallW}^{n,x}_{j} is the approximate slope in the cell $I_{j}$ obtained by using some slope limiter and $\chi_{j}(x)$ denotes the characteristic function of $I_{j}$.

In the 1D case, the Anderson-Witting model (2.2) reduces to

[TABLE]

whose analytical solution is given by

[TABLE]

where $v_{1}=p^{1}/p^{0}$ is the velocity of particle in $x$ direction, $x^{\prime}=x-v_{1}(t-t^{\prime})$ and $x^{\prime\prime}=x-v_{1}(t-t^{\prime\prime})$ are the particle trajectories, and $f_{0}$ is the initial particle velocity distribution function, i.e. f(x,0,\mbox{\boldmath\smallp})=f_{0}(x,\mbox{\boldmath\smallp}).

Taking the moments of (4.6) and integrating them over the space-time cell $I_{j}\times[t_{n},t_{n+1})$ yield

[TABLE]

where d\varXi=\frac{d^{3}\mbox{\boldmath\smallp}}{|\mbox{\boldmath\smallp}|}. Using the conservation constraints (2.5) gives

[TABLE]

which is the starting point of our 1D second-order accurate BGK scheme. If replacing the distribution f(x_{j\pm\frac{1}{2}},t,\mbox{\boldmath\smallp}) in (4.9) with an approximate distribution \hat{f}(x_{j\pm\frac{1}{2}},t,\mbox{\boldmath\smallp}), then one gets the following finite volume scheme

[TABLE]

where the numerical flux \hat{\mbox{\boldmath\smallF}}^{n}_{j+\frac{1}{2}} is given by

[TABLE]

with

[TABLE]

here $v_{1}=p^{1}/p^{0}$, $x^{\prime}=x_{j+\frac{1}{2}}-v_{1}(t-t^{\prime})$, $x^{\prime\prime}=x_{j+\frac{1}{2}}-v_{1}(t-t^{\prime\prime})$, f_{h,0}(x_{j+\frac{1}{2}}-v_{1}(t-t_{n}),\mbox{\boldmath\smallp})\approx f_{0}(x_{j+\frac{1}{2}}-v_{1}(t-t_{n}),\mbox{\boldmath\smallp}) and g_{h}(x^{\prime},t^{\prime},\mbox{\boldmath\smallp})\approx g(x^{\prime},t^{\prime},\mbox{\boldmath\smallp}). It is worth noting that it is very expensive to get $U_{\alpha}(x^{\prime\prime},t^{\prime\prime})$ and $U_{\alpha}(x^{\prime},t^{\prime})$ at the right hand side of (4.1). In practice, $U_{\alpha}(x^{\prime\prime},t^{\prime\prime})$ and $U_{\alpha}(x^{\prime},t^{\prime})$ in the first term may be approximated as $U_{\alpha,j+\frac{1}{2}}^{n}$ while $U_{\alpha}(x^{\prime},t^{\prime})$ in the second term may be simplified as $U_{\alpha,j+\frac{1}{2},L}^{n}$ or $U_{\alpha,j+\frac{1}{2},R}^{n}$ depending on the sign of $v_{1}$ and will be given in Section 4.1.1.

The remaining tasks are to derive the approximate initial velocity distribution function f_{h,0}(x_{j+\frac{1}{2}}-v_{1}(t-t_{n}),\mbox{\boldmath\smallp}) and equilibrium velocity distribution function g_{h}(x^{\prime},t^{\prime},\mbox{\boldmath\smallp}).

4.1.1 Equilibrium distribution $g_{0}$ at the point $(x_{j+\frac{1}{2}},t_{n})$

At the cell interface ${x=x_{j+\frac{1}{2}}}$, (4.5) gives the following left and right limiting values

[TABLE]

Using (2.4), \mbox{\boldmath\smallW}^{n}_{j+\frac{1}{2},L} and \mbox{\boldmath\smallW}^{n}_{j+\frac{1}{2},R} gives the J$\ddot{\text{u}}$ttner distributions at the left and right of cell interface $x=x_{j+\frac{1}{2}}$ as follows

[TABLE]

and the particle four-flow $N^{\alpha}$ and the energy-momentum tensor $T^{\alpha\beta}$ at the point $(x_{j+\frac{1}{2}},t_{n})$

[TABLE]

Using those and Theorem 2.1 calculates the macroscopic quantities $n^{n}_{j+\frac{1}{2}},T^{n}_{j+\frac{1}{2}}$, and $U_{{\alpha,j+\frac{1}{2}}}^{n}$, and then gives the J$\ddot{\text{u}}$ttner distribution function at the point $(x_{j+\frac{1}{2}},t_{n})$ as follows

[TABLE]

which will be used to derive the equilibrium velocity distribution g_{h}(x,t,\mbox{\boldmath\smallp}), see Section 4.1.3.

4.1.2 Initial distribution function f_{h,0}(x,\mbox{\boldmath\smallp})

Assuming that f(x,t,\mbox{\boldmath\smallp}) and g(x,t,\mbox{\boldmath\smallp}) are sufficiently smooth and borrowing the idea in the Chapman-Enskog expansion, f(x,t,\mbox{\boldmath\smallp}) is supposed to be expanded as follows

[TABLE]

with

[TABLE]

The conservation constraints (2.5) give the constraints on $A$ and $a$

[TABLE]

Setting $t=t_{n}$ and using (4.16) and the Taylor series expansion of f(x,t,\mbox{\boldmath\smallp}) with respect to $x$ from both sides of the cell interface $x=x_{j+\frac{1}{2}}$ give the following approximate initial non-equilibrium distribution function

[TABLE]

where $\tilde{x}=x-x_{j+\frac{1}{2}}$, $g_{L}$ and $g_{R}$ are given in (4.14), $(a_{L},A_{L})$ and $(a_{R},A_{R})$ are considered as the left and right limits of $(a,A)$ at the cell interface $x=x_{j+\frac{1}{2}}$ respectively. The slopes $a_{L}$ and $a_{R}$ come from the spatial derivative of J$\ddot{\text{u}}$ttner distribution and have unique correspondences with the slopes of the conservative variables $W$ by

[TABLE]

where

[TABLE]

Those correspondences form the linear system for the unknow \mbox{\boldmath\smalla}_{\omega}:=({a_{\omega,1},a_{\omega,2},a_{\omega,3}})^{T}

[TABLE]

where the coefficient matrix $M_{0}^{\omega}$ is given by

[TABLE]

Using the conservation constraints (4.18) and $a_{\omega}$ gives the linear system for $A_{\omega}$ as follows

[TABLE]

which can be cast into the following form

[TABLE]

with

[TABLE]

The rest is to calculate all elements of $M_{0}$ and $M_{1}$, whose superscript $L$ or $R$ has been omitted for the sake of convenience. In the ultra-relativistic limit, those can be exactly gotten. Because p^{0}=|\mbox{\boldmath\smallp}|, the triple integrals in $M_{0}$ and $M_{1}$ can be simplified by using polar coordinate transformation

[TABLE]

which implies d\Xi=|\mbox{\boldmath\smallp}|d|\mbox{\boldmath\smallp}|d\xi d\varphi. In fact, the above transformation can convert the triple integrals in the matrices $M_{0}$ and $M_{1}$ into a single integral with respect to |\mbox{\boldmath\smallp}| and a double integral with respect to $\xi$ and $\varphi$. On the other hand, in the 1D case, the integrands do not depend on the variable $\varphi$, so the double integral can further reduce to a single integral with respect to $\xi$ which can be exactly calculated. Those lead to

[TABLE]

and

[TABLE]

where

[TABLE]

4.1.3 Equilibrium velocity distribution g_{h}(x,t,\mbox{\boldmath\smallp})

Using \mbox{\boldmath\smallW}_{0}:=\mbox{\boldmath\smallW}_{j+\frac{1}{2}}^{n} derived in Section 4.1.1 and the approximate cell average values \mbox{\boldmath\small\bar{W}}_{j+1} and \mbox{\boldmath\small\bar{W}}_{j} reconstructs a cell-vertex based linear polynomial around the cell interface $x=x_{j+\frac{1}{2}}$ as follows

[TABLE]

where \mbox{\boldmath\smallW}_{0}^{x}=\frac{1}{\Delta x}(\mbox{\boldmath\small\bar{W}}_{j+1}-\mbox{\boldmath\small\bar{W}}_{j}). Again the Taylor series expansion of $g$ at the cell interface $x=x_{j+\frac{1}{2}}$ gives

[TABLE]

where $(a_{0},A_{0})$ are the values of $(a,A)$ at the point $(x_{j+\frac{1}{2}},t_{n})$. Similarly, the slope $a_{0}$ comes from the spatial derivative of J$\ddot{\text{u}}$ttner distribution and has a unique correspondence with the slope of the conservative variables $W$ by

[TABLE]

and then the conservation constraints and $a_{0}$ gives the following linear system

[TABLE]

Those can be rewritten as

[TABLE]

where \mbox{\boldmath\smalla}_{0}=(a_{0,1},a_{0,2},a_{0,3})^{T}, \mbox{\boldmath\smallA}_{0}=(A_{0,1},A_{0,2},A_{0,3})^{T}, and $M_{0}^{0}$ and $M_{1}^{0}$ can be calculated by (4.23) and (4.24) with $n,T$ and $U^{\alpha}$ instead of $n^{n}_{j+\frac{1}{2}},T^{n}_{j+\frac{1}{2}}$ and $U^{n,\alpha}_{j+\frac{1}{2}}$. Those systems can be solved by using the subroutine for (4.20) and (4.21).

Up to now, all parameters in the initial gas distribution function $f_{h,0}$ and the equilibrium state $g_{h}$ have been determined. Substituting (4.19) and (4.26) into (4.1) gives our distribution function $\hat{f}$ at a cell interface $x=x_{j+\frac{1}{2}}$ as follows

[TABLE]

where $H[x]$ is the Heaviside function defined by

[TABLE]

Finally, substituting (4.27) into the integral (4.11) yields the numerical flux \hat{\mbox{\boldmath\smallF}}^{n}_{j+\frac{1}{2}}.

4.2 2D Euler equations

This section extends the above BGK scheme to the 2D ultra-relativistic Euler equations

[TABLE]

where

[TABLE]

with $h=4T$, $p=nT$, and \mbox{\boldmath\smallu}=(u_{1},u_{2},0)^{T}. Four real eigenvalues of the Jacobian matrix A_{1}(\mbox{\boldmath\smallW})=\partial\mbox{\boldmath\smallF}/\partial\mbox{\boldmath\smallW} and A_{2}(\mbox{\boldmath\smallW})=\partial\mbox{\boldmath\smallG}/\partial\mbox{\boldmath\smallW} can be given as follows

[TABLE]

where $k=1,2$, and $c_{s}=\frac{1}{\sqrt{3}}$ is the speed of sound.

Divide the spatial domain $\Omega$ into a rectangular mesh with the cell $I_{i,j}=\{(x,y)|x_{i-\frac{1}{2}}<x<x_{i+\frac{1}{2}},y_{j-\frac{1}{2}}<y<y_{j+\frac{1}{2}}\}$, where $x_{i+\frac{1}{2}}=\frac{1}{2}(x_{i}+x_{i+1}),y_{j+\frac{1}{2}}=\frac{1}{2}(y_{j}+y_{j+1})$, $x_{i}=i\Delta x,y_{j}=j\Delta y$, and $i,j\in\mathbb{Z}$. The time interval $[0,T]$ is also partitioned into a (non-uniform) mesh ${t_{n+1}=t_{n}+\Delta t_{n},t_{0}=0,n\geqslant 0}$, where the time step size $\Delta t_{n}$ is determined by

[TABLE]

the constant $C$ denotes the CFL number, and $\bar{\varrho}^{k}_{i,j}$ denotes the approximation of the spectral radius of A_{k}(\mbox{\boldmath\smallW}) over the cell $I_{i,j}$, $k=1,2$.

The 2D Anderson-Witting model becomes

[TABLE]

whose analytical solution can be given by

[TABLE]

where $v_{1}=p^{1}/p^{0}$ and $v_{2}=p^{2}/p^{0}$ are the particle velocities in $x$ and $y$ directions respectively, $\{x^{\prime}=x-v_{1}(t-t^{\prime}),y^{\prime}=y-v_{2}(t-t^{\prime})\}$ and $\{x^{\prime\prime}=x-v_{1}(t-t^{\prime\prime}),y^{\prime\prime}=y-v_{2}(t-t^{\prime\prime})\}$ are the particle trajectories, and f_{0}(x,y,\mbox{\boldmath\smallp}) is the initial particle velocity distribution function, i.e. f(x,y,0,\mbox{\boldmath\smallp})=f_{0}(x,y,\mbox{\boldmath\smallp}).

Taking the moments of (4.31) and integrating them over $I_{i,j}\times[t_{n},t_{n+1})$ yield the 2D finite volume scheme

[TABLE]

where \mbox{\boldmath\small\bar{W}}^{n}_{i,j} is the cell average approximation of conservative vector \mbox{\boldmath\smallW}(x,y,t) over the cell $I_{i,j}$ at time $t_{n}$, i.e.

[TABLE]

and

[TABLE]

where \hat{f}(x_{i+\frac{1}{2}},y_{j},t,\mbox{\boldmath\smallp})\approx f(x_{i+\frac{1}{2}},y_{j},t,\mbox{\boldmath\smallp}) and \hat{f}(x_{i},y_{j+\frac{1}{2}},t,\mbox{\boldmath\smallp})\approx f(x_{i},y_{j+\frac{1}{2}},t,\mbox{\boldmath\smallp}). Because the derivation of \hat{f}(x_{i},y_{j+\frac{1}{2}},t,\mbox{\boldmath\smallp}) is very similar to \hat{f}(x_{i+\frac{1}{2}},y_{j},t,\mbox{\boldmath\smallp}), we will mainly derive \hat{f}(x_{i+\frac{1}{2}},y_{j},t,\mbox{\boldmath\smallp}) with the help of (4.32) as follows

[TABLE]

where $\tilde{t}=t-t_{n}$, $x^{\prime}=x_{i+\frac{1}{2}}-v_{1}(t-t^{\prime}),y^{\prime}=y_{j}-v_{2}(t-t^{\prime})$ and $x^{\prime\prime}=x_{i+\frac{1}{2}}-v_{1}(t-t^{\prime\prime}),y^{\prime\prime}=y_{j}-v_{2}(t-t^{\prime\prime})$, and f_{h,0}(x_{i+\frac{1}{2}},y_{j}-v_{1}\tilde{t},\mbox{\boldmath\smallp}) and g_{h}(x^{\prime},y^{\prime},t^{\prime},\mbox{\boldmath\smallp}) are (approximate) initial distribution function and equilibrium velocity distribution function, respectively, which will be presented in the following. Similarly, in order to avoid expensive cost in getting $U_{\alpha}(x^{\prime\prime},y^{\prime\prime},t^{\prime\prime})$ or $U_{\alpha}(x^{\prime},y^{\prime},t^{\prime})$ along the particle trajectory, $U_{\alpha}(x^{\prime\prime},y^{\prime\prime},t^{\prime\prime})$ and $U_{\alpha}(x^{\prime},y^{\prime},t^{\prime})$ in (4.36) may be taken as a constant $U_{\alpha,{i+\frac{1}{2},j}}^{n}$, and $U_{\alpha}(x^{\prime},y^{\prime},t^{\prime})$ in the second term may be replaced with $U_{\alpha,{i+\frac{1}{2},j},L}^{n}$ or $U_{\alpha,{i+\frac{1}{2},j},R}^{n}$ which is given in Section 4.2.1.

4.2.1 Equilibrium distribution $g_{0}$ at the point $(x_{i+\frac{1}{2}},y_{j},t_{n})$

Using the cell average values \{\mbox{\boldmath\small\bar{W}}^{n}_{i,j}\} reconstructs a piecewise linear function

[TABLE]

where \mbox{\boldmath\smallW}^{n}_{i,j}(x,y):={\mbox{\boldmath\small\bar{W}}}^{n}_{i,j}+\mbox{\boldmath\smallW}^{n,x}_{i,j}(x-x_{i})+\mbox{\boldmath\smallW}^{n,y}_{i,j}(y-y_{j}), \mbox{\boldmath\smallW}^{n,x}_{i,j} and \mbox{\boldmath\smallW}^{n,y}_{i,j} are the $x$- and $y$-slopes in the cell $I_{i,j}$, respectively, and $\chi_{i,j}(x,y)$ is the characteristic function of the cell $I_{i,j}$. At the point $(x_{i+\frac{1}{2}},y_{j})$, the left and right limiting values of \mbox{\boldmath\smallW}_{h}(x,y,t_{n}) are given by

[TABLE]

Similar to the 1D case, with the help of \mbox{\boldmath\smallW}^{n}_{i+\frac{1}{2},j,L}, \mbox{\boldmath\smallW}^{n}_{i+\frac{1}{2},j,R} and J$\ddot{\text{u}}$ttner distribution (2.4), one can get $g_{L}$ and $g_{R}$ at $(x_{i+\frac{1}{2}},y_{j},{t_{n}})$. Then the particle four-flow $N^{\alpha}$ and the energy-momentum tensor $T^{\alpha\beta}$ at $(x_{j+\frac{1}{2}},y_{j},t_{n})$ can be defined by

[TABLE]

Using those and Theorem 2.1, the macroscopic quantities $n^{n}_{i+\frac{1}{2},j},T^{n}_{i+\frac{1}{2},j}$ and $U_{{\alpha,i+\frac{1}{2},j}}^{n}$ can be calculated and then the J$\ddot{\text{u}}$ttner distribution function $g_{0}$ at $(x_{i+\frac{1}{2}},y_{j},t_{n})$ is obtained.

Similarly, in the $y$-direction, \mbox{\boldmath\smallW}^{n}_{i,j+\frac{1}{2},L} and \mbox{\boldmath\smallW}^{n}_{i,j+\frac{1}{2},R} can also be given by (4.37) so that one has corresponding left and right equilibrium distributions $\tilde{g}_{L}$ and $\tilde{g}_{R}$. The particle four-flow $N^{\alpha}$ and the energy-momentum tensor $T^{\alpha\beta}$ at $(x_{i},y_{j+\frac{1}{2}},t_{n})$ are defined by

[TABLE]

which give $n^{n}_{i,j+\frac{1}{2}},T^{n}_{i,j+\frac{1}{2}}$, $U_{{\alpha,i,j+\frac{1}{2}}}^{n}$ and $g_{0}$ at $(x_{i},y_{j+\frac{1}{2}},t_{n})$.

The following will derive the initial distribution function f_{h,0}(x,y,\mbox{\boldmath\smallp}) and equilibrium distribution g_{h}(x,y,t,\mbox{\boldmath\smallp}), separately.

4.2.2 Initial distribution function f_{h,0}(x,y,\mbox{\boldmath\smallp})

Borrowing the idea in the Chapman-Enskog expansion, f(x,y,t,\mbox{\boldmath\smallp}) is supposed to be of the form

[TABLE]

The conservation constraints (2.5) imply the constraints on $A,a$ and $b$

[TABLE]

Using the Taylor series expansion of $f$ at the cell interface $(x_{i+\frac{1}{2}},y_{j})$ gives

[TABLE]

where $\tilde{x}=x-x_{i+\frac{1}{2}},\tilde{y}=y-y_{j}$, and $(a_{\omega},b_{\omega},A_{\omega})$, $\omega=L,R$, are of the form

[TABLE]

The slopes $a_{\omega}$ and $b_{\omega}$ come from the spatial derivative of J$\ddot{\text{u}}$ttner distribution and have unique correspondences with the slopes of the conservative variables $W$ by the following linear systems for $a_{\omega}$ and $b_{\omega}$

[TABLE]

Those linear systems can also be expressed as follows

[TABLE]

where the coefficient matrix is defined by

[TABLE]

Substituting $a_{\omega}$ and $b_{\omega}$ into the conservation constraints (4.40) gives the linear systems for $A_{\omega}$ as follows

[TABLE]

which can be rewritten as

[TABLE]

where

[TABLE]

All elements of the matrices $M_{0}^{\omega}$, $M_{1}^{\omega}$ and $M_{2}^{\omega}$ can also be explicitly presented by using the coordinate transformation (4.22). If omitting the superscripts $L$ and $R$, then the matrices $M_{0}$, $M_{1}$, and $M_{2}$ are

[TABLE]

[TABLE]

and

[TABLE]

where $w^{1}=\xi,w^{2}=\sqrt{1-\xi^{2}}\sin\varphi,w^{3}=\sqrt{1-\xi^{2}}\cos\varphi$, and

[TABLE]

4.2.3 Equilibrium velocity distribution g_{h}(x,y,t,\mbox{\boldmath\smallp})

Using \mbox{\boldmath\smallW}_{0}:=\mbox{\boldmath\smallW}_{i+\frac{1}{2},j}^{n} derived in Section 4.2.1 and the cell averages \mbox{\boldmath\small\bar{W}}_{i+1,j} and \mbox{\boldmath\small\bar{W}}_{i,j} reconstructs a linear polynomial

[TABLE]

where \mbox{\boldmath\smallW}_{0}^{x}=\frac{1}{\Delta x}(\mbox{\boldmath\small\bar{W}}_{i+1,j}-\mbox{\boldmath\small\bar{W}}_{i,j}) and \mbox{\boldmath\smallW}_{0}^{y}=\frac{1}{2\Delta y}(\mbox{\boldmath\smallW}^{n}_{i+\frac{1}{2},j+1}-\mbox{\boldmath\smallW}^{n}_{i+\frac{1}{2},j-1}). Again using the Taylor series expansion of $g$ at the cell interface $(x_{i+\frac{1}{2}},y_{j})$ gives

[TABLE]

where $(a_{0},b_{0},A_{0})$ are the values of $(a,b,A)$ at the point $(x_{i+\frac{1}{2}},y_{j},t_{n})$. Similarly, the linear systems for $a_{0},b_{0}$ and $A_{0}$ can be derived as follows

[TABLE]

or

[TABLE]

where the elements of $M_{0}^{0},M_{1}^{0}$ and $M_{2}^{0}$ are given by (4.44), (4.45), and (4.46) with $n,T,U^{\alpha}$ instead of $n^{n}_{i+\frac{1}{2},j},T^{n}_{i+\frac{1}{2},j}$ and $U^{n,\alpha}_{i+\frac{1}{2},j}$.

Up to now, the initial gas distribution function $f_{h,0}$ and the equilibrium state $g_{h}$ have been given. Substituting (4.41) and (4.48) into (4.36) gives

[TABLE]

where $\tilde{t}=t-t_{n}$. Combining this \hat{f}(x_{i+\frac{1}{2}},y_{j},t,\mbox{\boldmath\smallp}) with (4.34) can get the numerical flux \hat{\mbox{\boldmath\smallF}}^{n}_{i+\frac{1}{2},j}. The numerical flux \hat{\mbox{\boldmath\smallG}}^{n}_{i,j+\frac{1}{2}} can be obtained in the same procedure.

4.3 2D Navier-Stokes equations

Because the previous simple expansion (4.16) or (4.39) cannot give the Navier-Stokes equations (3.11)-(3.13), one has to use the complicate Chapman-Enskog expansion (3.8)-(3.9) to design the genuine BGK schemes for the Navier-Stokes equations. On the other hand, for the Navier-Stokes equations, calculating the macroscopic quantities $n,U^{\alpha}$, and $p$ needs the value of the fluxes \mbox{\boldmath\smallF}^{k} besides $W$. More specially, one has to first calculate the energy-momentum tensor $T^{\alpha\beta}$ and particle four-flow $N^{\alpha}$ from the kinetic level and then use Theorem 2.1 to calculate $n,U^{\alpha}$, and $p$. It shows that there exists a very big difference between the genuine BGK schemes for the Euler and Navier-Stokes equations.

In order to obtain $T^{\alpha\beta}$ and $N^{\alpha}$ at $t=t_{n+1}$ from the kinetic level, multiplying (4.31) by $p^{k}/p^{0}$ gives

[TABLE]

Taking the moments of (4.31) and (4.50) and integrating them over the space-time domain $I_{i,j}\times[t_{n},t_{n+1})$ , respectively, yield

[TABLE]

where

[TABLE]

Our task is to get the approximate distributions $\hat{f}(x_{i+\frac{1}{2}},y_{j},t)$ and $\hat{f}(x_{i},y_{j+\frac{1}{2}},t)$ for the numerical fluxes and $\hat{f}(x_{i},y_{j},t)$ and $g(x_{i},y_{j},t)$ for the source terms. The following will focus on the derivation of $\hat{f}(x_{i+\frac{1}{2}},y_{j},t)$ with the help of the analytical solution (4.32) of the 2D Anderson-Witting model.

4.3.1 Initial distribution function f_{h,0}(x,y,t,\mbox{\boldmath\smallp})

This section derives the initial distribution function $f_{h,0}$ for $\hat{f}(x_{i+\frac{1}{2}},y_{j},t)$. The Chapman-Enskog expansion (3.8)-(3.9) is rewritten as follows

[TABLE]

where $A^{ce}=A_{\beta}^{ce}p^{\beta}+A^{ce}_{4}$, $a^{ce}=a^{ce}_{\beta}p^{\beta}+a^{ce}_{4}$, $b^{ce}=b^{ce}_{\beta}p^{\beta}+b^{ce}_{4}$, $c^{ce}=c^{ce}_{\beta}p^{\beta}+c^{ce}_{4}$, and

[TABLE]

It is observed from those expressions of $A^{ce},a^{ce},b^{ce}$, and $c^{ce}$ that one has to compute the time derivatives, which are not required in the Euler case. Those time derivatives are approximately computed by using the following second-order extrapolation method: for any smooth function $h(t)$, the first order derivative at $t=t_{n}$ is numerically obtained by

[TABLE]

Using the Chapman-Enskog expansion (4.53) and the Taylor series expansion in terms of $x$ gives the initial velocity distribution

[TABLE]

where $\tilde{x}=x-x_{i+\frac{1}{2}},\tilde{y}=y-y_{j}$, $g_{L}$ and $g_{R}$ denote the left and right J$\ddot{\text{u}}$ttner distributions at $x_{i+\frac{1}{2}}$ with $y=y_{j},t=t_{n}$, the Taylor expansion coefficients $(a_{L},b_{L})$ and $(a_{R},b_{R})$ are calculated by using the same procedure as in the Euler case, while the Chapman-Enskog expansion coefficients $a_{L}^{ce},a_{R}^{ce},b_{L}^{ce},b_{R}^{ce},c_{L}^{ce},c_{R}^{ce}$ and $A_{L}^{ce},A_{R}^{ce}$ are calculated by (4.54).

4.3.2 Equilibrium distribution functions g_{h}(x,y,t,\mbox{\boldmath\smallp})

In order to obtain the equilibrium distribution functions g_{h}(x,y,t,\mbox{\boldmath\smallp}) for $\hat{f}(x_{i+\frac{1}{2}},y_{j},t)$, the particle four-flow $N^{\alpha}$ and the energy-momentum tensor $T^{\alpha\beta}$ at $(x_{i+\frac{1}{2}},y_{j})$ and $t=t^{n}$ are defined by

[TABLE]

where $f_{L}$ and $f_{R}$ are the left and right limits of $f_{h,0}$ with $y=y_{j}$ at $x=x_{i+\frac{1}{2}}$. Using those definitions and Theorem 2.1, the macroscopic quantities $n^{n}_{i+\frac{1}{2},j},T^{n}_{i+\frac{1}{2},j}$ and $U_{\alpha,i+\frac{1}{2},j}^{n}$ can be obtained, and then one gets the J$\ddot{\text{u}}$ttner distribution function $g_{0}$ at $(x_{i+\frac{1}{2}},y_{j},t_{n})$. Similar to Section 4.2.3, we reconstruct a cell-vertex based linear polynomial and do the first-order Taylor series expansion of $g$ at the cell interface $(x_{i+\frac{1}{2}},y_{j})$, see (4.48). However, it is different from the Euler case that $A_{0}$ is obtained by

[TABLE]

where \mbox{\boldmath\smallW}^{t}_{0} is calculated by using the second-order extrapolation (4.55). After those, substituting $f_{h,0}$ and $g_{h}$ into (4.36) gets $\hat{f}(x_{i+\frac{1}{2}},y_{j},t)$. The distribution $\hat{f}(x_{i},y_{j+\frac{1}{2}},t)$ can be similarly obtained.

4.3.3 Derivation of source terms \mbox{\boldmath\smallS}_{{1},i,j} and \mbox{\boldmath\smallS}_{{2},i,j}

The rest is to calculate $\hat{f}(x_{i},y_{j},t)$ and $g(x_{i},y_{j},t)$ for the source terms \mbox{\boldmath\smallS}_{{1},i,j} and \mbox{\boldmath\smallS}_{{2},i,j}. The procedure is the same as the above except for taking the first-order Taylor series expansion at the cell-center $(x_{i},y_{j})$. To be more specific, $g$ and $f_{0}$ in the analytical solution (4.36) of 2D Anderson-Witting model are replaced with

[TABLE]

and

[TABLE]

where $(a_{0},b_{0},A_{0})$ are the Taylor expansion coefficients at $(x_{i},y_{j},t_{n})$ calculated by the same procedure as that for $\hat{f}(x_{i+\frac{1}{2}},y_{j},t)$, $\tilde{x}=x-x_{i}$, $\tilde{y}=y-y_{j}$, $g_{0}$ denotes the J$\ddot{\text{u}}$ttner distribution at $(x_{i},y_{j},t_{n})$, $a_{0}^{ce},b_{0}^{ce},c_{0}^{ce}$ and $A_{0}^{ce}$ are the Chapman-Enskog expansion coefficients at $(x_{i},y_{j},t_{n})$. It is worth noting that since $f_{h,0}$ is continuous at $(x_{i},y_{j})$, there is no need to consider whether the left or right states should be taken here. The subroutine for the coefficients in (4.56) can be used to get those in f_{h,0}({x,y},\mbox{\boldmath\smallp}).

In order to define the equilibrium state $g(x_{i},y_{j},t)$ in the source term, firstly we need to figure out the corresponding macroscopic quantities such as $N^{\alpha}$ and $T^{\alpha\beta}$ which can be obtained by taking the moments of $\hat{f}(x_{i},y_{j},t)$. Using the Theorem 2.1, the macroscopic quantities such as $n,T$ and $u$ can be obtained. Thus the J$\ddot{\text{u}}$ttner distribution function at cell center $(x_{i},y_{j})$ is derived according to the definition.

Until now, all distributions are derived and the second-order accurate genuine BGK scheme (4.51) is developed for the 2D ultra-relativistic Navier-Stokes equations.

5 Numerical experiments

This section will solve several 1D and 2D problems on the ultra-relativistic fluid flow to demonstrate the accuracy and effectiveness of the present genuine BGK schemes, which will be compared to the second-order accurate BGK-type and KFVS schemes [1, 37]. The collision time $\tau$ is taken as

[TABLE]

with $\tau_{m}=\frac{5\mu}{4p}$ for the viscous flow and $\tau_{m}=C_{1}\Delta t^{\alpha}_{n}$ for the inviscid flow, $C_{1}$, $C_{2}$ and $\alpha$ are three constants, $P_{L},P_{R}$ are the left and right limits of the pressure at the cell interface, respectively. Unless specifically stated, this section takes $C_{1}=0.001,C_{2}=1.5$ and $\alpha=1$, the time step-size $\Delta t_{n}$ is determined by the CFL condition (4.4) or (4.30) with the CFL number of 0.4, and the characteristic variables are reconstructed with the van Leer limiter.

5.1 1D Euler case

Example 5.1** (Accuracy test)**

To check the accuracy of our BGK method, we first solve a smooth problem which describes a sine wave propagating periodically in the domain $\Omega=[0,1]$. The initial conditions are taken as

[TABLE]

and corresponding exact solutions are given by

[TABLE]

The computational domain $\Omega$ is divided into $N$ uniform cells and the periodic boundary conditions are specified at $x=0,1$.

Table 5.1 gives the $l^{1}$- and $l^{2}$-errors at $t=0.2$ and corresponding convergence rates for the BGK scheme with $\alpha=2$ and $C_{1}=C_{2}=1$. The results show that a second-order rate of convergence can be obtained for our BGK scheme although the van Leer limiter loses slight accuracy.

Example 5.2** (Riemann problem I)**

This is a Riemann problem with the following initial data

[TABLE]

The initial discontinuity will evolve as a left-moving shock wave, a right-moving contact discontinuity, and a right-moving shock wave. Fig. 5.1 displays the numerical results at $t=0.5$ and their close-ups obtained by using our BGK scheme (“$\circ$”), the BGK-type scheme (“$\times$”), and the KFVS scheme (“+”) with 400 uniform cells in the domain $[0,1]$, where the solid lines denote the exact solutions. It can be seen that our BGK scheme resolves the contact discontinuity better than the second-order accurate BGK-type and KFVS schemes, and they can well capture such wave configuration.

Example 5.3** (Riemann problem II)**

The initial conditions of the second Riemann problem are

[TABLE]

Fig. 5.2 shows the numerical solutions at $t=0.5$ obtained by using our BGK scheme (“$\circ$”), the BGK-type scheme (“$\times$”), and the KFVS scheme (“+”) with 400 uniform cells within the domain $[0,1]$, where the solid line denotes the exact solution. It is seen that the solutions consist of a left-moving rarefaction wave, a contact discontinuity, and a right-moving shock wave, the computed solutions well accord with the exact solutions, and the rarefaction and shock waves are well resolved. Moreover, our BGK scheme exhibits better resolution of the contact discontinuity than the BGK-type and KFVS schemes.

Example 5.4** (Riemann problem III)**

The initial data are

[TABLE]

The initial discontinuity will evolve as a left-moving rarefaction wave, a stationary contact discontinuity, and a right-moving rarefaction wave.

Fig. 5.3 plots the numerical results at $t=0.5$ obtained by using our BGK scheme (“$\circ$”), the BGK-type scheme (“$\times$”), and the KFVS scheme (“+”) with 400 uniform cells in the domain $[0,1]$, where the solid line denotes the exact solution. It is seen that there is a undershoot near the contact discontinuity in the number density which usually happens in the non-relativistic cases.

Example 5.5** (Perturbed shock tube problem)**

The initial data are

[TABLE]

where $n_{r}=0.125-0.0875\sin(50(x-0.5))$. It is a perturbed shock tube problem, which has widely been used to test the ability of the shock-capturing schemes in resolving small-scale flow features in the non-relativistic flow.

Fig. 5.4 plots the numerical results at $t=0.5$ in the computational domain $\Omega=[0,1]$ obtained by using our BGK scheme (“$\circ$”), the BGK-type scheme (“$\times$”), and the KFVS scheme (“+”) with 400 uniform cells. Those are compared with the reference solution (the solid line) obtained by using the KFVS scheme with a finer mesh of 10000 uniform cells. It is seen that the shock wave is moving into a sinusoidal density field, some complex but smooth structures are generated at the left hand side of the shock wave when the shock wave interacts with the sine wave, and our BGK scheme is obviously better than the BGK-type and KFVS schemes in resolving those complex structures. Since the continuity equation in the Euler equations decouples from other equations for the pressure and velocity, one does not see the effect of perturbation in the pressure [18].

Example 5.6** (Collision of blast waves)**

It is about the collision of blast waves and simulated to evaluate the performance of the genuine BGK scheme and the BGK-type and KFVS schemes for the flow with strong discontinuities. The initial data are taken as follows

[TABLE]

Reflecting boundary conditions are specified at the two ends of the unit interval $[0,1]$.

Fig. 5.5 plots the numerical results at $t=0.75$ obtained by using our BGK scheme (“$\circ$”), the BGK-type scheme (“$\times$”), and the KFVS scheme (“+”) with 700 uniform cells within the domain $[0,1]$. It is found that the solutions at $t=0.75$ are bounded by two shock waves and those schemes can well resolve those shock waves. However, the genuine BGK scheme exhibits better resolution of the contact discontinuity than the BGK-type and KFVS schemes.

5.2 2D Euler case

Example 5.7** (Accuracy test)**

To check the accuracy of our BGK scheme, we solve a smooth problem which describes a sine wave propagating periodically in the domain $\Omega=[0,1]\times[0,1]$ at an angle $\alpha={45}^{\circ}$ with the $x$-axis. The initial conditions are taken as follows

[TABLE]

so that the exact solution can be given by

[TABLE]

The computational domain $\Omega$ is divided into $N\times N$ uniform cells and the periodic boundary conditions are specified.

Table 5.2 gives the $l^{1}$- and $l^{2}$- errors at $t=0.1$ and corresponding convergence rates for the BGK scheme with $\alpha=2$ and $C_{1}=C_{2}=1$. The results show that the 2D BGK scheme is second-order accurate and the van Leer limiter affects the accuracy.

To verify the capability of our genuine BGK scheme in capturing the complex 2D relativistic wave configurations, we will solve three inviscid problems: explosion in a box, cylindrical explosion, and ultra-relativistic jet problems.

Example 5.8** (Implosion in a box)**

This example considers a 2D Riemann problem inside a squared domain $[0,2]\times[0,2]$ with reflecting walls. A square with side length of 0.5 embedded in the center of the outer box of side length of 2. The number density is 4 and the pressure is 10 inside the small box while both the density and the pressure are 1 outside of the small box. The fluid velocities are zero everywhere.

Figs. 5.6 and 5.7 give the contours of the density, pressure and velocities at time $t=3$ and $12$ obtained by our BGK scheme on the uniform mesh of $400\times 400$ cells, respectively. The results show that the genuine BGK scheme captures the complex wave interaction well. Fig. 5.8 gives a comparison of the numerical densities along the line $y=1$ calculated by using the genuine BGK scheme (“$\circ$”), BGK-type scheme (“$\times$”), and KFVS scheme (“+”) respectively. Obviously, the genuine BGK scheme resolves the complex wave structure better than the BGK-type and KFVS schemes.

Example 5.9** (Cylindrical explosion problem)**

Initially, there is a high-density, high-pressure circle with a radius of 0.2 embedded in a low density, low pressure medium within a squared domain $[0,1]\times[0,1]$. Inside the circle, the number density is 2 and the pressure is 10, while outside the circle the number density and pressure are 1 and 0.3, respectively. The velocities are zero everywhere.

Fig. 5.9 displays the the contour plots at $t=0.2$ obtained by using the BGK scheme on the mesh of $200\times 200$ uniform cells. The results show that a circular shock wave and a circular discontinuity travel away from the center, and a circular rarefaction wave propagates toward the center of the circle. Fig. 5.10 gives a comparison of the number density and pressure along the line $y=0.5$ obtained by the BGK, BGK-type, and KFVS schemes, respectively. The symbols “$\circ$” , “$\times$” and “+” denote the solutions obtained by using the BGK, BGK-type and KFVS schemes. It can be observed that all of them give closer results. However, the BGK scheme resolves the discontinuities better than the KFVS.

Example 5.10** (Ultra-relativistic jet)**

The dynamics of relativistic jet relevant in astrophysics has been widely studied by numerical methods in the literature [3, 58, 29]. This test simulates a relativistic jet with the computational region $[0,12]\times[-3.5,3.5]$ and $\alpha=C_{1}=C_{2}=1$. The initial states for the relativistic jet beam are

[TABLE]

where the subscripts $b$ and $m$ correspond to the beam and medium, respectively.

The initial relativistic jet is injected through a unit wide nozzle located at the middle of left boundary while a reflecting boundary is used outside of the nozzle. Outflow boundary conditions with zero gradients of variables are imposed at the other part of the domain boundary. Fig. 5.11 shows the numerical results at $t=5,6,7,8$ obtained by our BGK scheme on the mesh of $600\times 350$ uniform cells. The average speed of the jet head is 0.91 which matches the theoretical estimate 0.87 in [29].

5.3 Navier-Stokes case

This section designs two examples of viscous flow to test the genuine BGK scheme (4.51) for the ultra-relativistic Navier-Stokes equations. Because the extrapolation (4.55) requires the numerical solutions at $t=t_{n-1}$ and $t_{n-2}$, the “initial” data at first several time levels have to be specified for the BGK scheme in advance. In the following examples, the macroscopic variables at $t=t_{0}+0.5\Delta t_{0}$ and $t_{0}+\Delta t_{0}$ are first obtained by using the initial data, time partial derivatives at $t=t_{0}$, and BGK scheme proposed in Section 4.3, where the first order partial derivatives in time are derived by using the exact solutions. Then, the time partial derivatives at $t=t_{0}+\Delta t_{0}$ for the macroscopic variables are calculated by using the extrapolation (4.55), and the solutions are further evolved in time by the BGK scheme with the extrapolation (4.55).

Example 5.11** (longitudinally boost-invariant system)**

For ease of numerical implementation, this test focuses on the longitudinally boost-invariant systems. They are conveniently described in curvilinear coordinates $x_{m}=(\tilde{t},y,z,\eta)$, where $\tilde{t}=\sqrt{t^{2}-x^{2}}$ is the longitudinal proper time, $\eta=\frac{1}{2}\ln\left(\frac{t+x}{t-x}\right)$ is the space-time rapidity and $(y,z)$ are the usual Cartesian coordinates in the plane transverse to the beam direction $x$. The systems are realized by assuming a specific “scaling” velocity profile $u_{1}=x/t$ along the beam direction, and the initial conditions are independent on the longitudinal reference frame (boost invariance), that is to say, they do not depend on $\eta$. The readers are referred to [43] for more details.

Our computations consider the boost-invariant longitudinal expansion without transverse flow, so that the relativistic Navier-Stokes equations read

[TABLE]

Since $u_{1}=\frac{x}{t}$, $U^{0}=t/\tilde{t}$ and $U^{1}=x/\tilde{t}$, it holds that $\partial_{\alpha}U^{\alpha}=1/\tilde{t}$. Thus the equation for $n$ becomes

[TABLE]

The analytical solutions can be given by

[TABLE]

where $C_{1}=p_{0}(t_{0}^{2}-x_{0}^{2})^{{\frac{2}{3}}}-\frac{4\mu}{3}(t_{0}^{2}-x_{0}^{2})^{{\frac{1}{6}}}$ and $C_{2}=n_{0}\sqrt{t_{0}^{2}-x_{0}^{2}}$. We take $x_{0}=0,t_{0}=1,p_{0}=1,n_{0}=1$, $\mu=0.0005$, and $\Omega=[-\frac{t_{0}}{2},\frac{t_{0}}{2}]$. Moreover, the time partial derivatives of $n,{u_{1}},p$ at $t={t_{0}}$ are given by the exact solution.

Fig. 5.12 shows the number density, velocity and pressure at $t=1.2$ obtained by our 1D BGK scheme with 20 cells (“${\triangle}$”) and 40 cells (“$\circ$”), respectively. The results show that the numerical results predicted by our BGK scheme fit the exact solutions very well. Table 5.3 lists the $l^{1}$- and $l^{2}$-errors at $t=1.2$ and corresponding convergence rates for our BGK scheme. Those data show that a second-order rate of convergence can be obtained by our BGK scheme.

Example 5.12** (Heat conduction)**

This test considers the problem of heat conduction between two parallel plates, which are assumed to be infinite and separated by a distance $H$. Moreover, both plates are always stationary. The temperatures of the lower and upper plates are given by $T_{0}$ and $T_{1}$, respectively. The viscosity $\mu$ is a constant.

Based on the above assumptions, the Navier-Stokes equations can be simplified as

[TABLE]

whose analytic solution is gotten as follows

[TABLE]

Our computation takes $H=1,p=0.8,u_{1}=0.2,u_{2}=0,\mu=5\times 10^{-3},T_{0}=0.1,T_{1}=1.0002T_{0}$, and $0.5(T_{0}+T_{1})$ as the initial value for the temperature $T$ in the entire domain. Moreover, the initial time partial derivatives are given by $n_{t}(x,0)=0,v_{1t}(x,0)=0,v_{2t}(x,0)=0$ and $p_{t}(x,0)=0$. Because $u_{1}\neq 0$, the the 2D BGK scheme should be used for numerical simulation.

The left figure in Fig. 5.13 plots the numerical temperature (“$\circ$”) obtained by the 2D BGK scheme in comparison with the steady-state analytic solution (solid line) given by (5.6). It is seen that the numerical solution is well comparable with the analytic. The right figure in Fig. 5.13 shows convergence of the temperature to the steady state measured in the $l^{1}$-error between the numerical and analytic solutions.

6 Conclusions

The paper developed second-order accurate genuine BGK schemes in the framework of finite volume method for the 1D and 2D ultra-relativistic flows. Different from the existing KFVS or BGK-type schemes for the ultra-relativistic Euler equations the present genuine BGK schemes were derived from the analytical solution of the Anderson-Witting model, which was given for the first time and included the “genuine” particle collisions in the gas transport process. The genuine BGK schemes were also developed for the ultra-relativistic viscous flows and two ultra-relativistic viscous examples were designed. Several 1D and 2D numerical experiments were conducted to demonstrate that the proposed BGK schemes were accurate and stable in simulating ultra-relativistic inviscid and viscous flows, and had higher resolution at the contact discontinuity than the KFVS or BGK-type schemes. The present BGK schemes could be easily extended to the 3D Cartesian grid for the ultra-relativistic flows and it was interesting to develop the genuine BGK schemes for the special and general relativistic flows.

Acknowledgements

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