NONCONSERVATIVE HYPERBOLIC SYSTEMS IN FLUID MECHANICS

This paper is devoted to the numerical approximation of nonconservative hyperbolic systems. More precisely, we consider the bitemperature Euler system and we propose two methods of discretization. The first one is a kinetic approach based on an underlying kinetic model. The second one deals with a Suliciu approach when magnetic fields are taken into account. Résumé. Cet article est consacré à l’approximation numérique de systèmes hyperboliques non conservatifs. On considère plus précisément le système d’Euler bitempérature pour lequel on a développé deux méthodes de discrétisation. La première approche est basée sur l’approximation d’un modèle cinétique sous-jacent tandis que la seconde correspond à une méthode de type Suliciu, lorsque les champs magnétiques sont pris en compte.


Introduction
This paper is devoted to the numerical approximation of some nonconservative hyperbolic systems arising in plasma physics. Classically, discontinuous solutions of such systems are not well-defined and several notions of weak solutions can exist. More precisely, there is a lack of Rankine-Hugoniot condition in order to define jump relations at a discontinuity. At the theoretical point of view, a general theory based on family of paths has been proposed in [20]. However the generalisation of this theory at the numerical point of view ( [23]) is unsatisfactory. Indeed even if the path is known, the numerical scheme is not able to give the right solution ( [2]). Moreover, the numerical solution depends on the viscosity of the scheme. Hence in the presence of shocks, different numerical schemes give different plateaux. In [19], the authors assume that the model is isentropic for the electrons and change the system into a conservative system. In ( [15]), the methodolgy of ( [23]) has been improved by controling the numerical viscosity in order to recover the right path.
In the present paper, we present two alternative methods that have been presented in ( [4]). The first one consists in considering an underlying kinetic model that is able to recover the bitemperature model from an hydrodynamic limit. Moreover the time and space discretisation of the kinetic model is classical because no nonconservative products appear. This procedure is also adapted to the discrete BGK models of the Aregba-Natalini formalism ( [5]). The second technique is based on the Suliciu method ( [24]) for a transverse magnetic configuration. The Suliciu relation system basically consists in linearizing the flux fonction, and more precisely the underlying pressure term, which encompasses the genuine nonlinearities of the model, while keeping unchanged the convective part of the system and the linearly degenerate characteristic field. This procedure enables the derivation of an extended relaxation system which is consistent with both the original system and its entropy inequality in the limit of a vanishing relaxation parameter. By extended, we mean that it contains additional variables associated with the linearization process. The relaxation models are hyperbolic with only linearly degenerate characteristic fields, so that the associated Riemann problem is easy to solve. This makes the whole procedure particularly adapted to construct efficient Godunov-type methods based on approximate Riemann solvers. Moreover the method can be designed to satisfy a discrete entropy inequality.
The structure of the paper is the following. In the second part, we present a unified approach of the kinetic methods in the spirit of [9]. We are able here to prove the compatibility of those models with the entropy properties of their hydrodynamic limits. In the third part, we present a relaxation scheme for the bitemperature Euler system taking into account magnetic fields and generalizing the results of ( [4]). We obtain in particular, an entropic scheme.

Kinetic models for a nonconservative Euler system
In this work we aim to define a kinetic framework for the nonconservative Euler system modeling a plasma with a constant ionization Z = ne ni , where n e and n i represent the electronic and ionic concentrations. This situation corresponds to the quasi-neutral regime.

Notations
The notations are the following. We use the subscripts e for the electrons, i for the ions. c e , c i are the mass fractions whereas m e and m i are the electronic and ionic molecular masses. The electronic and ionic densities ρ e and ρ i are defined by ρ e = ρc e = m e n e , ρ i = ρc i = m i n i .
So, by using c e + c i = 1 and the quasineutrality assumption, c e and c i are constant. Therefore, the model will be described with only one equation of conservation for the density. The electronic and ionic velocities u e and u i will be assumed at thermodynamical equilibrium in the model. So u e = u i = u, where u represents the velocity of the mixture. Ionic and electronic energies E e and E i defined by: will be described by two different dynamics. Hence we have two pressure laws and two temperatures: Finally, the system that is studied reads as a nonconservative system with source-terms as follows: Nonconservativity is due to the presence of source-terms and of products u∂ x p α . If γ e = γ i then (ρ, ρu, E e + E i ) satisfies the classical 3 × 3 Euler system for which solutions are known. Otherwise p e + p i cannot be written as a function of ρu and E e + E i . But even if γ e = γ i one has to find E e and E i and hence the solving of equations with nonconservative products cannot be avoided.
The system (1) is hyperbolic diagonalizable with eigenvalues: The fields 1 and 4 are genuinely nonlinear, while the fields 2 and 3 are linearly degenerate. In order to define admissible solutions, we have been able to exhibit a dissipative entropy for (1): denoting the Euler entropy for each species α = e, α = i, for U = (ρ, ρu, E e , E i ) we set If U is a smooth solution of (1) then The problem of defining weak admissible solutions for nonconservative systems has been addressed by several authors, (see [20], [7] and references therein). One can also refer to the numerical studies ( [19], [6], [14]). In the present work we use the kinetic approach to define the admissibility of solutions. Hence a weak solution is admissible if it can be obtained as a limit of an entropy-compatible underlying kinetic model.

Underlying kinetic (BGK) models for the conservative Euler system
We start from BGK models for the Euler monotemperature equations. Denoting U = (ρ, ρu, E) ∈ Ω ⊂ R 3 , F (U ) = (ρu, ρu 2 + p, u(E + p)), the Euler system is a system of conservation laws: We follow the framework proposed by F. Bouchut in [9]. We define a measure space (X, dξ), a real valued function a defined on X, a "maxwellian function" M from R 3 × X onto R k , and a "moment operator" P from X to L(R k , R 3 ) such that for all U ∈ Ω: Formally if lim ε→0 f ε = f , then f (x, t, ξ) = M (U (f )(x, t), ξ) and U (f ) is a solution of (2). In [9], conditions are given for the existence of microscopic entropies compatible with all the entropies of the macroscopic limit.

Example 1: a BGK model for a polyatomic gas with a continuous energy variable
Here and M (U, v, I) ∈ R, with α = α(γ).
As f ε (x, t, v, I) ∈ R, it is a rank one model, (see [3], [13]). This model is compatible with the physical entropy of Euler system.

Example 2: a discrete velocity BGK model
Here (see [5]). For any Euler entropy, the existence of related microscopic entropies is ensured under Liu's subcharacteristic condition, see [9]:
In order to approximate the nonconservative products we define a linear operator N such that In the case of example 1: In the case of example 2: The equations for f ε e and f ε i are coupled with the ones for the electric field E: The source-term B αβ is such that if ε → 0, then When ε tends to 0, we have formally Quasineutrality holds: ρ = ρc e = ρc i and c e , c i are constant. By taking the moments, it comes that Considering Hence U = (ρ, ρu, E e , E i ) is a solution of system (1).
Theorem 2.1. Suppose that there exists microscopic entropies for the kinetic model (3) related to the entropy η. Let U be a solution of the Euler bitemperature model (1) obtained by passing to the limit in (3). Then one has the entropy inequality We define such a solution U as an admissible solution.
In the case of examples 1 and 2 above, the microscopic entropies exist, see [4] for details.

The conservative case
In the monotemperature case the numerical flux for the transport equation is chosen: (U n j ) j being known we set By applying the moment operator P : In the case of example 1 one chooses and therefore In the case of example 2 (λ 1 ≤ 0 ≤ λ 2 ) we apply the upwind scheme to approximate f n j (i). We obtain the same scheme as for example 1. But in the present case, we are able to prove discrete entropy inequalities. Indeed, the set of velocities being bounded, we get the monotonicity properties of our scheme.

2.4.2.
An entropic numerical scheme in the bitemperature case U 0 = (U 0 j ) j∈Z being given, we set: • Step 1: projection onto equilibrium Suppose that n ≥ 0 being fixed, we have U n , U n e , U n i with ρ α = ρc α . We define Coupling with Maxwell-Ampère and Poisson equations: As q e = −e and q i = Ze, we get as in the continous case ρ n+1 e,j = c e ρ n+1 The numerical flux per species is defined as and is consistent with (2). We then compute E n+1 as in the continuous case: The nonconservative products are approximated by defining δ n Finally the numerical scheme reads as follows: Theorem 2.2. In the example 2 with the subcharacteristic condition, the following discrete entropy inequality holds:

Double rarefaction wave
We consider a double rarefaction test case, for ν ei = 0, where it is possible to compute an exact solution. Consider a Riemann problem with left data ρ − , u − , T e,− , T i,− and right data ρ + , u + , T e,+ , T i,+ defined by In this test case, the kinetic scheme is compared with the exact solution and a relaxation scheme based on a Suliciu method. The exact solution is composed of two rarefaction waves travelling in opposite direction. Fig.  1 compares the ionic and the electronic temperatures for the kinetic scheme, a Suliciu scheme and the exact solution.
2.6. Implosion test case with γ e = 9 7 , γ i = 7 5 This case has been studied in [21] with γ e = γ i = 5 3 . It is a Riemann problem whose left data are ρ − , u − , T e,− , T i,− and right data are ρ + , u + , T e,+ , T i,+ . They are given by:  cells and a standard second order extension in space and time. As γ i = γ e the density and the velocity depend on the value of the coupling coefficient ν ei . We consider two values : ν ei = 0, and the one given in the NRL plasma formulary [22]. Densities and velocities are displayed in Figure 2. In Figure 3, electronic and ionic temperatures are depicted. As expected at this computation time, T e = T i when ν ei = 0.

Bitemperature MHD model
The aim of this part is to generalise the results of ( [4]) by taking into account magnetic fields. In particular, we consider in the following the case of a transverse magnetic field.

Introduction
In this part we shall not be able to impose any condition (except the entropy inequalities) in order to define weak solutions to the system so that the solutions may differ from the ones computed with other numerical schemes. Nevertheless, we consider any possible solution that satisfies an entropy inequality. The novelty of this work is to derive a new efficient numerical method. More precisely we obtain an explicit finite volume scheme, which is firstly consistant with smooth solutions. Moreover, in case of discontinuous solutions, we are able to exhibit sufficient CFL conditions, under which the scheme is conservative with respect to the physical conserved quantities, preserves the positivity of densities and internal energies electrons and ions and is entropy satisfying.

Main results
The nonconservative system we deal with is obtained from a kinetic conservative model. More precisely, we firstly introduce an underlying conservative kinetic model coupled to Maxwell equations. Then the bitemperature Euler system with transverse magnetic field is established from this kinetic model by hydrodynamic limit. This point is an extension of section 2 and is not developped here. In this work, we derive of a finite volume method to approximate weak solutions. It is obtained by solving a relaxation system of Suliciu type, and is similar to HLLC type solvers. The solver is shown in particular to preserve positivity of density and internal energies. Moreover we use a local minimum entropy principle to prove discrete entropy inequalities, ensuring the robustness of the scheme.
In the following we introduce the nonconservative system we want to deal with and we present the numerical scheme we have used to approximate weak solutions of the system. Then we briefly point out the main properties of the scheme and we give some numerical results.

Bitemperature MHD system
If dependency is only one spatial variable x, the nonconservative two species Euler equations with transverse magnetic field are given by: Here u = (u 1 , u 2 ) is the average velocity of the plasma, B 3 is the vertical component of the magnetic field.
The source terms for exchanges between electron and ion species S ei and S ie are defined by where ν 1,ei ≥ 0 is the frequency exchange between temperatures, ν 2,ei , ν 2,ie ≥ 0 are frequencies due to the drift velocy u 2 ± ∂ x B 3 . We do not mention here the explicit formulas here because we want to focus on the homogeneous part of the system. Moreover the homogeneous system associated to (4) -(9) is endowed with an entropy inequality: where s α , α = e, i is the classical specific entropy Here C is a nonnegative constant.

Numerical approximation
The numerical approximation we use has two steps: • First step. We use an approximate Riemann solver for the homogeneous system. Let U n+1/2 be the obtained solution. • Second step. We take the temperatures interaction into account implicitly: the approximate solution of system at time t n+1 is defined by This system is linear and owns an explicit solution.
From here we will focus on the first step of the numerical approximation. We will explain how to derive an efficient approximate Riemann solvers for the homogeneous part of the system (4) - (9). Indeed a finite volume scheme for this homogeneous quasilinear system is classically built following Godunov's approach, by considering piecewise constant approximation of and invoking an approximate Riemann solver at the interface between two cells.
In order to get those approximate Riemann solvers, we use a standard relaxation approach, introduced in [8] for the gas dynamic equations. This approach has been developed in [11] for the MHD equations, in [10] for shallow elastic fluids, in [12] for shallow water MHD equations, in [4] for the bitemperature Euler system. An abstract general description can be found in [18], and related works are [16,17]. This technique enables to naturally handle the entropy inequality (11), and also preserves the positivity of density and internal energies.

Resolution of the homogeneous system
We introduce new variables π e , π i , the relaxed pressures, and a intended to parametrize the speed. The form of the relaxation system is as follows, The approximate Riemann solver can be defined as follows, starting from left and right values U l , U r at an interface.
• Solve the Riemann problem of the system (14)- (22) with initial data obtained by completing U l , U r by the equilibrium relations π e,L = p e,L ≡ (γ e − 1)ρ L ε e,L , π i,L = p i,L ≡ (γ i − 1)ρ L ε i,L , π e,R = p e,R ≡ (γ e − 1)ρ R ε e,R , and with suitable positive values of a l , a r that will be discussed further on, essentially in Section 3.4.4. This step is a projection procedure. It corresponds to add in the equations (20, 21) the terms − ρ τ (p e −π e ) and − ρ τ (p i − π i ) and let the parameter τ tends to 0. • Retain in the solution only the variables ρ, ρu 1 , ρu 2 , E e , E i , B 3 . The result is a vector called R(x/t, U l , U r ).
Intuitively, the solver is consistent because of the equations (14)- (19), that are consistent with the equations (4)- (9). The specific values used for a do not play any role in this consistency.

Godunov scheme
Following the Godunov approach, the numerical scheme can be defined by the approximate Riemann solver as follows. We consider a mesh of cells discrete times t n with t n+1 − t n = ∆t, and cell values U n i approximating the average of U over the cell i at time t n . We can then define an approximate solution U appr (t, x) for t n ≤ t < t n+1 and x ∈ R by where x i = (x i−1/2 + x i+1/2 )/2. This definition is coherent under a half CFL condition, formulated as The new values at time t n+1 are defined by Notice that it is only in this averaging procedure that the choice of the particular pseudo-conservative variable U as (13) is involved. We can follow the computations of [8, Section 2.3], the only difference being that the system is not conservative. We obtain the update formula where The variable ξ stands for x/t, and the pseudo-conservative flux is chosen as In (28), the fourth and fifth components could be chosen differently since the two energy equations in our system are not conservative. We can remark that the choice of F has no influence on the update formula (26).

Intermediate states
In the solution to the Riemann problem, the speeds corresponding to the previous eigenvalues will be denoted by: Thus we get a 3-wave solver with two intermediate states. The variables take the values "L" for x/t < Σ 1 , "L*" for Σ 1 < x/t < Σ 2 , "R*" for Σ 2 < x/t < Σ 3 , "R" for Σ 3 < x/t, see Figure 4. We skip the details but Figure 4. Intermediate states in the Riemann solution using the Riemann invariants of the relaxation system, we obtain the following intermediate states: For α = e, i we have Finally, using previous formulas one can compute the speeds The maximal propagation speed is then The corresponding CFL condition is Remark 3.1. Let us notice that c i π e − c e π i is equal to c i w 1,e − c e w 1,i . As a consequence, c i π e − c e π i is a Riemann invariant for both extreme eigenvalues. This means that this quantity remains constant through the related contact discontinuities, so that u∂ x (c i π e − c e π i ) = 0 there. For the central discontinuity, u is constant so that u∂ x (c i π e − c e π i ) = ∂ x (u (c i π e − c e π i )) this product is also well defined in the usual weak sense.

Positivity of density and internal energies
Here we provide some sufficient conditions on a L and a R in order to satisfy (35) and the realisability of the intermediate states, that is the positivity of ρ * L , ρ * R , ε * e,L , ε * i,L , ε * e,R and ε * i,R . First, we obtain that are sufficient conditions to preserve the positivity of ρ * L and ρ * R . Second, from a straightforward calculation using (30), ε * e,L , ε * i,L , ε * e,R and ε * i,R are positive if

Numerical tests
In this section we perform numerical approximations in order to evaluate both accuracy and robustness of the scheme.
Firstly, Fig.1 compares the results given by the Suliciu method for the system (1) with an analytical solution for two rarefaction waves.
Next Here we solve the homogeneous part of the system (4) - (9) . We consider two test cases. The first one corresponds to a classical shock tube Sod problem whereas the second one is a symmetric rarefaction waves problem. The left and right states of the Riemann problems are written in table 1 in terms of density, velocity, transverse magnetic field, electronic and ionic temperatures. The numerical solution for test 1 consists of, from left to right, a left rarefaction wave, a material contact and a right shock. The results are shown in Figure 5. We observe numerically the right structure of waves and the convergence of the scheme.
The numerical solution for test 2 consists of two rarefaction waves travelling in opposite directions. The results are displayed in Figure 6 and show no oscillations in the results. This illustrates the robustness of the scheme when dealing with vacuum data.