High Order Central Schemes Applied to Relativistic Multi-component Flow Models

The dynamics of inviscid multi-component relativistic fluids may be modeled by the relativistic Euler equations, augmented by one (or more) additional species equation(s). We use the high-resolution staggered central schemes to solve these equations. The equilibrium states for each component are coupled in space and time to have a common temperature and velocity. The current schemes can handle strong shocks and the oscillations near the interfaces are negligible, which usually happens in the multi-component flows. The schemes also guarantee the exact mass conservation for each component, the exact conservation of total momentum, and energy in the whole particle system. The central schemes are robust, reliable, compact and easy to implement. Several one-and two-dimensional numerical test cases are included in this paper, which validate the application of these schemes to relativistic multi-component flows.


Introduction
In recent years, relativistic gas dynamics plays an important role in areas of astrophysics, high energy particle beams, high energy nuclear collisions, and free-electron laser technology.The equations that describe the relativistic gas dynamics are highly nonlinear.For the practical problems, it is difficult to solve these equations analytically, therefore numerical solutions are perused.Several numerical methods for solving relativistic gas dy-namics have been reported.All these methods are mostly developed out of the existing reliable methods for solving the Euler equations of non-relativistic or Newtonian gas dynamics.
The first attempt to solve the equations of relativistic gas dynamics (RGD) was made by Wilson [1] [2] using an Eulerian explicit finite difference code with monotonic transport.The code relies on artificial viscosity technique [3] to handle shock wave.Despite of its popularity it turned out to be unable to accurately describe the extremely relativistic flows [4].In mid eighties, Norman and Winkler [5] proposed a reformulation of the difference equations with artificial viscosity consistent with relativistic dynamics of non-perfect fluids.Dean et al. [6] used flux correcting algorithms for RGD equations in the context of heavy ion collisions.
The development of numerical methods for the non-relativistic multi-component flows have attracted much attention in the past years, for example Fedkiw et al. [21] [22], Karni [23]- [25], Karni and Quirk [26], Marquina and Mulet [27].Moreover, Xu [28] used BGK-based gas-kinetic schemes to solve multi-component flows, while Lian and Xu [29] used the same scheme in order to solve the multi-component flows with chemical reactions.
This paper is an extension of the relativistic Euler equations to multi-component flows.We use the high-resolution non-oscillatory central schemes of Nessyahu and Tadmor [30] as well as Jiang and Tadmor [31] to solve these Euler equations.In this study, we consider only two-component flow, however, an extension to further components will result in addition of a continuity equation for the corresponding species.The central schemes are predictor-corrector methods which consist of two steps: starting with given cell averages, we first predict point values which are based on the non-oscillatory piecewise-linear reconstructions from the cell averages.At the second corrector step, we use staggered averaging, together with the predicted mid values, to realize the evolution of these averages.This results in a second-order, non-oscillatory central scheme.
The second order accuracy of these schemes is based on MUSCL-type reconstruction.Like upwind schemes, the reconstructed piecewise-polynomials used by the central schemes, also make use of non-linear limiters which guarantee the overall non-oscillatory nature of the approximate solution.But unlike the upwind schemes, central schemes do not require the intricate and time-consuming (approximate) Riemann solvers which are essential for the high-resolution upwind schemes.This advantage is especially important in the multi-dimensional case where there is no exact Riemann solver.Moreover, the central schemes are "genuinely multi-dimensional" in the sense that it does not necessitate dimensional splitting.Apart from these, central schemes do not produce spurious oscillations, such as carbuncle phenomena and odd-even decoupling which usually happens in the Godunov upwind schemes.The reason of this advantage is the presence of sufficient numerical dissipation in the central schemes.This is also the reason here in the multi-component flow that we see almost negligible or no oscillations at the gases interface [32].
The organization of this paper is as follows.In Section 2, we derive the three-dimensional Euler equations for the relativistic multi-component flows.We then discuss how to obtain the primitive variables from the conserved variables.In Section 3, we write the one-dimensional relativistic Euler equations for the dynamics of a mixture of two gases.Starting from the first order central scheme, we explain the high-resolution second order central schemes to solve these Euler equations, see [30] and references therein.In Section 4, we explain the scheme for the two-dimensional relativistic multi-component flows.We again start from the first order central schemes and then extend it to second order, see [31] and references therein.In Section 5, we present numerical test cases which include propagation of one-dimensional (1D) relativistic blast waves, collision, cylindrical explosion, interaction of an air shock with helium bubble and explosion in a square box.Finally, Conclusion is given in Section 6.

Multi-Component Relativistic Euler Equations
For simplicity, we assume a model of mixture of two gases.An extension to more components is analogous.Let ρ ρ ρ = + denotes the total density of a mixture, with 1 ρ and 2 ρ as the rest mass densities of the first and second components, respectively.Also let 1 Y and 2 Y be the mass fractions of the first and second components.
We assume that both components are in thermal equilibrium and are perfect gases with specific heats at constant volume C and ratios of specific heats 1 γ , 2 γ .By stan- dard thermodynamic arguments, the ratio of specific heats γ of the mixture of gases is Using the Einstein summation convention the equations describing the motion of a two-component relativistic fluid are given by the six conservation laws ; ; where ( ) , 0, ,3 , µ ν =  and µ denote the covariant derivative with respect to coordinate x µ .Furthermore, u µ is the four-velocity of the mixture, and T µν is the stress-energy tensor, which for a perfect fluid can be written as .

T hu u pg
Here g µν is metric tensor 1, 0, 1, 1, 2, 3, 0, , ρ is the mixture density, p the fluid average pressure, and h the specific enthalpy of the fluid mixture defined by 1 , where ε is the specific internal energy.Note that we use natural units (i.e., the speed of light 1 c = ) through- out this study.In Minskowski space time and Cartesian coordinates ( ) , , , t x x x , the conservation equations given in Equation (2) can be written as with the conserved variables w and fluxes i f given as The six conserved quantities , , , , D D S S S and τ are the rest-mass densities of the two components, the three components of momentum density, and the energy density (measured relative to the rest mass density), respectively.They are all measured in laboratory frame, and are related to quantities in the local rest frame of the fluid (primitive variables) through the relations where i v are the components of three-velocity of fluid 0 , 1, 2, 3, The system of equations in Equation ( 5) along with definitions in Equations ( 6)-( 8) is closed by mean of an ideal equation of state (EOS) as given below ( ) .
We denote by s c the sound speed, defined by 2 , where s is the specific entropy, which is conserved along fluid lines.For EOS under consideration the speed of sound can be written as .
The Mach number of the flow is due to Königl [33] .
For any given initial macroscopic variables in space and time, ( ) ( ) The common values of density ρ , velocity i v , and pressure p can be obtained from the conservation re- quirements, .
From the above equations p, ρ and i v can be obtained by first solving an implicit function of pressure whose zero represents the pressure [7] [34].We have to find the root of the equation ) ( ) The monotonicity of ( ) [ ] is obtained from Equation (7) by taking in to account that (in our units) 1 v ≤ .Knowing p, Equation (15) then directly gives v and density.Similar to Aloy et al. [34], we obtained the solution ( ) 0 p η = by means of New- ton-Rahphson iteration in which the derivative of η , i.e. η′ , is approximated by where s c * is the speed of sound given by ( ) This approximation tends to the exact derivative when the solution is approached.On the other hand, it easily allows one to extend the present algorithm to general equation of state [34].

One-Dimensional Multi-Component Flows
Here, we are looking for spatially one-dimensional solutions of the two-component relativistic Euler equations.We only consider the solutions which depend on t and 1 x x = and satisfy ( ) , , , 0, 0 ( ) , p p t x = .The three dimensional system in Equation ( 2) then reduces to with the conserved variables w and fluxes f given as , , where where 0 2 1 u u = + .For any given initial macroscopic variables in space and time, ( ) ( ) the common values of density ρ , velocity v, and pressure p can be obtained from the conservation require- ments, From the above equations p, ρ , and v can be obtained by following the same procedure as given in relations ( 8)- (11).

One-Dimensional Central Schemes
Let us begin by introducing the well-known first order Lax-Friedrichs (LxF) scheme for one-dimensional conservation laws.This first order scheme is then extended to a second order central scheme, see [30].We consider a piecewise-constant initial data, ( ) ( ) x x + of integration.This leads to the LxF scheme where t x λ ∆ = ∆ .The piecewise constant cells in each step are staggered with respect to those in the previous step.

Extension to Higher Order
Starting with a piecewise-constant solution in time and space, ( ) , one reconstruct a piecewise linear (MUSCL-type) approximation in space, namely ( ) where x i w abbreviates a first-order discrete slopes, see Figure 1.A possible computation of these slopes, which results in an overall non-oscillatory scheme (consult [30]), is given by family of discrete derivatives parameterized wit 1 2 θ ≤ ≤ , i.e., for any grid function { } i w we set , , , , . 2 Here, ∆ denotes the central differencing, This interpolate, is then evolved exactly in time and projected on the staggered cell-averages on the next time step, Consider the balance law over the control volume , , Here, The averaging of the linear data in Equation ( 26) at So far everything is exact.Moreover, the Courant-Friedrichs-Levy (CFL) condition guarantees that ( ) ) , are smooth functions of τ ; hence they can be integrated approximately by the midpoint rule at the expense of an ( ) Putting Equation (31) in Equation ( 30), we finally get ( ) ( ) By Taylor expansion and the conservation laws in Equation ( 19), we have This may serve as our approximate mid-values 1 2 n i w + within the permissible second-order accuracy requirement.Here, ( ) stands for an approximate numerical derivatives of the flux ( ) ( ) .
The fluxes

( )
x i f w are computed by applying the min-mod limiter to each of the component of f, i.e., , , , , . 2 Here, ∆ denotes the central differencing, ( ) ( ) , and MM denotes the min-mod nonlinear limiter given by Equation ( 28).This component wise approach is one of the main advantages offered by central schemes over corresponding characteristic decompositions required by upwind schemes [30] [31].It is important to emphasize that while using the central type LxF solver, we integrate over the entire Riemann fan, which consists of both the left and right going waves.On the one hand, this enables us to ignore the detailed knowledge about the exact (or approximate) generalized Riemann solver.On the other hand, this enables us to accurately compute the numerical flux, whose values are extracted from the smooth interface of two non-interacting Riemann problems.In summary, this family of central differencing scheme takes the easily implemented predictor-corrector form, ( ) ( ) ( )

Two-Dimensional Multi-Component Flows
Here we are looking for spatially two-dimensional solutions of the multi-component Euler equations.We only consider the solutions which depend on t, , , t x y , , , ,0 v v t x v t y = and ( ) , , p p t x y = .The three-dimensional system in Equation ( 5) then reduces to ( ) ( ) 0.
f w g w w t x y The conserved variables w and fluxes f, g are given by where Here, 0 2 2 1 2 1 u u u = + + .For any given initial macroscopic variables in space and time, ( ) ( ) The common values of density ρ , velocity i v , and pressure p can be obtained from the conservation re- quirements, .
From the above equations p, ρ and i v can be obtained by following the same procedure as given in Equa- tions ( 7)- (11).

Two-Dimensional Central Schemes
To approximate Equation (36), we begin with a piecewise constant solution w , the approximate cell-average at n t t = , associated with the cell , , , i j x y χ is a characteristic function of the cell , i j C .The arguments applied to the one-dimensional case can be easily extended to the higher dimensions.In the following, we will abbreviate to denote the normalized integral, i.e., normalized over its length, area, etc.Also let { } As given in Figure 2, the first integral has contribution from the four cells Simplifying the above balance law we finally get the following LxF scheme [31]

A Second-Order Extension in 2D
A two-dimensional extension of the second order central scheme was introduced in [31].As in one-dimensional case, this staggered scheme can be viewed as an extension to the first-order LxF Scheme.A piecewise-linear interpolant is reconstructed from the calculated cell-averages at time n t , ( ) , , , .w are discrete slopes in the x-and y-directions, respectively, which are reconstructed from the given cell averages.To guarantee second-order accuracy, these slopes should approximate the corresponding derivatives, , , , , , .
x n y n i j i j i j i j w x w t x y O x w y w t x y O y x x A possible computation of these slopes, which results in an overall non-oscillatory schemes is given by family of discrete derivatives parameterized with 1 2 θ ≤ ≤ , for example Here MM denotes the min-mod nonlinear limiter given by Equation ( 28).This guarantees that the corresponding piecewise-linear reconstruction in Equation ( 45), ( ) , , n w t x y , is co-monotone with the underlying piecewise-constant approximation, ( ) , , , . Similar to one-dimensional case, the construction of the central scheme proceeds with a second step of an exact evolution.The integration of Equation ( 36) over volume We begin by evaluating the cell average ( ) , , d d , we find an average of the reconstructed polynomial in (45), ( ) ( ) . 4 16 x x y y w t x y x y w w w x y x y w w w Continuing in a counter clockwise direction, we have . 4 16 w t x y x y w w w . 4 16 w t x y x y w w w . 4 16 By adding the last four integrals, we find that the exact staggered averages of the reconstructed solution at So far everything is exact.We now turn to approximating the four fluxes on the right of Equation (48), starting with the one along the east face (consult Figure 3), i.e.

∫ ∫
We use midpoint quadrature rule for second-order approximation of the temporal integral, ( ) , , d d , ,  , , d , and, for the reasons to be clarified below, we use the second-order rectangular quadrature rule for the spatial integration across the y-axis, yielding n n j n n t i i j i j t y J f w t x y y t f w f w In similar manner, we approximate the remaining fluxes, ( ) ( ) n n i n n t j i j i j t x I g w t x y x t g w g w n n j n n t i i j i j t y J f w t x y y t f w f w n n i n n t j i j i j t x I g w t x y x t g w g w The fluxes in Equations ( 54)-(57) use the midpoint values, , , , and we take advantage of utilizing these mid-values for the spatial integration by the rectangular rule.Namely, since these mid-values are secured at the smooth center of their cells, , i j C , bounded away from the jump discontinuities along the edges, we may use Taylor expansion, ( ) ( ) are one-dimensional discrete slopes of the fluxes in the x-and y-directions, of the type reconstructed in Equation (46).We find these slopes in the same way as done for the conservative filed variables using min-mod procedure.Inserting these values, together with the staggered averages computed in Equation ( 53), into Equation (48), we get new staggered averages at 1 , In summary, we end up with a simple two-step predictor-corrector scheme in Equations ( 58) and (59).Starting with the cell averages,

Numerical Test Cases
Here we present one-and two-dimensional numerical problems in order to validate the application of central schemes for the solution of multi-component flow problems.
Problem 1: Propagation of relativistic blast waves: , , , , 10.0, 0.0,13.33,1.4,1.0 if 0.5, , , , , 1.0, 0.0, 0.66 10 ,1.67,1.0 if 0.5, where the computational domain 0 1 x ≤ ≤ is subdivided into 400 mesh points.This test problem has been con- sidered by several authors in one-component case, for example, Hawley, Smarr and Wilson [35], Schneider et al. [16], and Martí and Müller where the computational domain is 0 1 x ≤ ≤ with 400 mesh points.This problem was first considered in sin- gle-component case by Norman and Winkle [36].The flow pattern is similar to that of problem 1, but more extreme as shown Figure 5.In case of 4000 mesh points the relativistic effects reduces the post-shock state to a thin dense shell with a width of only about 2% of the grid length at t = 0.35.The fluid in the shell moves with  (i.e., shell 3.35 Γ = ), while jump in density in the shell reaches a value of 8.17.Problem 3: A similar problem in non-relativistic case was solved by Quirk and Karni [26]  Here, the computational domain is 0 1 x ≤ ≤ .We compute the approximate solution of this problem with a grid of 400 mesh points at time 0.65 t = . The reference solution is obtained at 3000 mesh points.The results are shown in Figure 6.

− ×
The velocities are zero everywhere.The specific heat ratios for helium and air are 1.67 and 1.4, while specific heats at constant volume are equal to 1 for both air and helium.The solution consists of a circular shock wave propagating outwards from the origin, followed by a circular contact discontinuity propagating in the same direction, and a circular rarefaction wave travelling towards the origin.The results are shown in Figure 7 for 400 mesh points at 15 .0 = t .Problem 5: A 1.16 Ms = shock wave in air hits a Helium cylindrical bubble: In this example, we introduce a single planar shock, moving in the air, with a cylindrical bubble of Helium.A schematic description of computational set-up is shown in Figure 8, where reflection boundary conditions are used on the upper and lower boundaries, while out flow boundary conditions on are used on the left and right boundaries.The bubble is assumed to be in both thermal and mechanical equilibrium with the surrounding air.The non-dimentionalized initial data are      heat ratios for helium and air are 1.67 and 1.4, while specific heats at constant volume are equal to 2.42 and 0.72, respectively.The results are shown in Figure 10.

Conclusion
A high-resolution staggered central scheme was applied to solve the relativistic multicomponent flow equations.The equilibrium states for each component are coupled in space and time to have a common temperature and velocity.Several case studies of the one-and two-dimensional flows were considered to validate the schemes.It was found that the suggested schemes have capability to resolve strong shocks without spurious oscillations.The schemes also guarantee the exact mass conservation for each component and the exact conservation of total momentum and energy in the whole particle system.The central schemes were found to be robust, reliable, compact and easy to implement for such complicated model equations.
uniqueness of the solution.The lower bound of the phys- ically allowed domain, min p MM denotes the min-mod nonlinear li- miter

Figure 2 .
Figure 2. Floor plane of the staggered grid.

.
It has contribution from the four intersecting cells, , Starting with the intersecting cell , i j C at the corner (see Fig- ure 2),
use the first-order predictor in Equation (58) for the evolution of the mid- followed by the second-order corrector in Equation (59) for computation of the new cell averages, This results in a second-order accurate non-oscillatory central schemes.As in the one- dimensional case, no exact (approximate) Riemann solvers are involved.The non-oscillatory behavior of the scheme things on the reconstructed discrete slopes, x w ,

4 . Problem 2 :
[7] [12], etc.It involves the formation of an intermediate state bounded by a shock wave propagating to the right and transonic rarefaction wave propagating to the left.The fluid in the intermediate state moves at a mildly relativistic speed ( right.Flow particles accumulate in a dense shell behind the shock wave compressing the fluid by a factor of 5 and heating it up to values of internal energy much larger than the rest-mass energy.Hence the fluid is extremely relativistic in thermo-dynamical point of view, but mildly relativistic dynamically.The results are shown in Figure The initial data are:

Problem 4 :
Cylindrical Explosion Problem: Consider a square domain [ ] [ ] 0,1 0,1 × .The initial data are constant in two regions separated by a circle of radius 0.2 centered at ( ) 0.5, 0.5 .Inside the circle is helium with density 10.0 and pressure 13.33, while outside is air of density 1.0 and pressure equal to 6 0.066 10 .

Figure 10 .
Figure 10.Problem 6 (explosion in a box) results at different times. .
. It consists of a 1.14 Ms = shock tube filled with air, where shock wave moves to the left.In the pre-shock wave stage, a bubble of Helium is set.The initial data are as follow