Modelling and numerical study of the polyatomic bitemperature Euler system

1.   Introduction

This work is devoted to the modelling and the numerical approximation of the nonconservative polyatomic bitemperature Euler system in the context of plasma physics. Physically, this model describes the interaction of one species of ions and one species of electrons, under the quasi-neutrality assumption.

The aim of this paper is more precisely to provide a construction and an approximation of the polyatomic Euler bitemperature system. This model is constituted of conservative equations for the mass and the momentum and two nonconservative equations for ionic and electronic energies. It is a variant of the system considered in [23]. The nonconservativity is due to source-terms but also to the presence of products of the velocity by pressure gradients. In Inertial Confinement Fusion applications, high temperature solutions involve shocks, for which those nonconservative products have to be determined. This can be done by defining paths, as proposed by Dal Maso, Le Floc'h and Murat ([24]). However, it is shown in [1] that the numerical adaptation of this theory given in [36] is delicate.

In order to compute physically consistent shocks, one can use an underlying kinetic system. This approach was usefully adopted in [3] for the monoatomic case. The model is constituted by a kinetic system coupled with Ampère and Poisson equations. Moreover, this construction leads to a kinetic numerical scheme starting from a semi-discretization in space and time of the kinetic model. A DVM approach of this model has also been considered in [20] together with an asymptotic preserving discretisation toward the bitemperature Euler system. However relaxation schemes and discrete BGK schemes have been developed in the polytropic situation for a general $\gamma$ law including the polyatomic case in [3], [6]. It seems important to consider situations where each species has its own $\gamma$ constant. Note that for smooth solutions, global in time existence has been proved under the condition that these constants are distinct, (see [5]).

We construct in this paper an extension of the model developed in [3] to a polyatomic setting with a continuous energy variable. This model is based on the use of several attractors like in [15] and is shown to satisfy an H theorem. We refer to [10], [13], [14], [19], [11], [34] for other BGK models for polyatomic gas mixtures.

In the present paper, the derivation is based again on an hydrodynamic limit performed from an underlying polyatomic kinetic model. In the present case, the unknowns of the kinetic equations are the distribution functions $f(t,x,v,I)$ depending on time $t$, space $x$, velocity $v$ and of one-dimensional internal energy parameter $I>0$. This energy parameter collects vibrational and rotational energies. Kinetic models with continuous energy variables have been introduced in [16] where the motivation was to develop a Monte-Carlo method. In [17], the authors derived a mesoscopic model of Boltzmann type associated to the previous microscopic model. This collision operator satisfies fundamental properties (H theorem, $\dotsc$). The generalization to mixtures has been performed in [26]. In [9], a compressible Navier-Stokes model has been derived for mixture of monoatomic and polyatomic gases. For some applications of such models we refer to [35], [38], [29]. In [35] the authors analyse the shock wave structure of some polyatomic gases. So by using an ESBGK model ([2], [21]) they show the presence of a double layer structure that is specific to the polyatomic setting. We mention that in [27], [30], the authors develop polyatomic models with a discrete internal energy variable.

As observed in [37], [3], it is possible to find an underlying kinetic model where the force term is related to the nonconservative terms. One advantage of the kinetic model, is its conservative form. In the present paper, the kinetic model describing the interspecies interaction is a two component polyatomic BGK model based on one continuous internal energy variable coupled with Ampère and Poisson equations. Hence starting from a standard semi-discretization of this model, the hydrodynamic limit leads to a numerical scheme for the bitemperature model. It must be noted that the obtained scheme is the same as the one obtained in section 3.2 of [3]. In this article however the scheme was obtained by a very different method involving models developed in [7] which are not based on a physically realistic kinetic approximation of the equations. The novelty here is that our polyatomic BGK model gives a physical justification to it in the general case where $\gamma_e$ and $\gamma_i$ may be distinct. In [4] diffusive terms have been constructed from a Chapman-Enskog expansion up to order $1$ leading to a Navier-Stokes type model. This last model generalizes the model that is studied in [22]. In [31], the authors perform a Chapman-Enskog expansion by introducing a small parameter representing the ratio between electronic and ionic molecular masses. They obtain an hyperbolic system with a parabolic regularisation on the electrons. The structure of shock wave solutions for this last model is considered in [39].

The plan of the paper is the following. Section $2$ deals with the presentation of the different models that are used in this paper. In particular, we detail the eigenstructure of the bitemperature system. In section $3$, the fluid model is obtained starting from the polyatomic model. In section $4$ exact solutions for the model are computed and in section $5$, the numerical scheme is developed. The last part is devoted to some numerical tests in 1D and 2D configurations.

2.   The mathematical models

2.1.   The bitemperature Euler system

2.1.1.   The fluid model

The bitemperature Euler system reads

where $\rho > 0$ represents the total density of the plasma, $u$ is the average velocity of the plasma, $\mathcal{E}_e$ and $\mathcal{E}_i$ are the total electronic and ionic energies. $T_e$ and $T_i$ represent the electronic and ionic temperatures. The coefficient $\nu_{ei}$ is nonnegative. The physically realistic formulas which are given in [33] make the source term very stiff.

One has

where $\rho_e = n_e m_e$, $\rho_i = n_i m_i$ are electronic and ionic densities, electronic and ionic concentrations $n_e$ and $n_i$ being assumed to be linked by $Z = n_e / n_i \ge 1$. $Z$ is considered as constant. This situation corresponds to the quasi-neutral regime. $m_e$ and $m_i$ represent the electronic and ionic masses. The mass fractions

are then constant and $c_e$ and $c_i$ are given by

The total energies are linked to the internal electronic and ionic energies by

The electronic and ionic pressures and temperatures are related by $p_e = n_e k_B T_e$ and $p_i = n_i k_B T_i$. The electronic and ionic internal energies are then given by

where $\gamma_e, \gamma_i$ $\in$ $[1,3]$, and $k_B$ is the Boltzmann constant.

In the following we denote $\mathcal U = (\rho,\rho u, \mathcal{E}_e, \mathcal{E}_i)$, $U_ \beta = (\rho_ \beta,\rho_ \beta u , \mathcal{E}_ \beta)$.

2.1.2.   Hyperbolicity and characteristic fields

We complete here the calculations presented in [5]. The system is rewritten by using the variables $\mathcal{V} = (\rho,u,\varepsilon_e,\varepsilon_i)$:

The matrix of the system (5) writes

We get $4$ eigenvalues

where

The value of $a$ given by (7) corresponds to the global sound velocity which yields for the classical Euler system $\sqrt{\gamma p /\rho}$. The eigenvectors associated to the eigenvalue $u$ are equal to

This system is then hyperbolic diagonalisable. The fields 2 and 3 are linearly degenerate. Consider now the fields $1$ and $4$. The eigenvectors are

and

Hence the fields $1$ and $4$ are genuinely nonlinear.

2.2.   The kinetic model

In this section, we generalize the BGK model considered in [3] to a polyatomic setting with a continuous energy variable. We firstly precise the notations and next introduce the BGK model that is proved to satisfy the right conservation properties and an H theorem.

2.2.1.   Notations

Kinetic models for a mixture of two polyatomic gases are described by two distribution functions $f_{\beta}$ of species $\beta$ depending on time $t \in \mathbb{R}^+$, space $x \in \mathbb{R}^3$, velocity $v \in \mathbb{R}^3$ and on internal energy $I \in \mathbb{R}_+$.

Hydrodynamic quantities of species $\beta$ are defined for $\alpha_{\beta } \geq 0$ by

and

Velocities and temperatures of the mixture $u$ and $T$ are defined by

The parameters $\alpha_e$ and $\alpha_i$ are related to $\gamma_e$ and $\gamma_i$ by the formula ([26])

For example in the diatomic case, we have $\alpha_e = \alpha_i = 0$. The values $\alpha_e = \alpha_i = -1$ correspond to $\gamma_e = \gamma_i = 5/3$ which is a monoatomic mixture.

Define the entropy of the mixture by

and the entropy flux by

2.2.2.   A polyatomic BGK model

In this section we consider the following kinetic model for $\beta \in \{e;i\}$,

with

where

As suggested in [32] and developed in [3] the definition of $u^{\#}$ and $T^{\#}$ as

and

gives the conservations of impulsion and total energy for the model (13, 14, 15). These definitions allow in particular to consider a model such that the molecular mass ratios of the species are different. This situation corresponds to the case $\tau_{ei} \not = \tau_{ie}$.

The model (13, 14, 15, 16, 17) is coupled to Ampère and Poisson equations through the electric field $E$ as

$j$ represents the plama current, $\overline{\rho}$ the total charge and $\varepsilon_{0}$ is the permittivity. $j$ and $\overline{\rho}$ are defined by

When one of the two components is monoatomic the model has to be slightly modified. In [8], the authors study for moments systems the link beween polyatomic and monoatomic models. The connexion between monoatomic and polyatomic models can be made by passing to the limit in some parameters. When the two components are monoatomic we refer to [3].

In the case of one monoatomic component and one polyatomic component the model has to be clearly written. Consider for example that electrons are monoatomic whereas ions remain polyatomic. In that case, the model is described with the distributions $f_{e}(t,x,v)$ and $f_{i}(t,x,v,I)$. In that case, the model reads

with

The hydrodynamic quantities of the monoatomic species are computed as

The total impulsion $\rho u$ and the total energy $\mathcal{E}$ for the system reads

In that case $u^{\#}$ is given by (16) and $T^{\#}$ becomes

As previously, $u^{\#}$ and $T^{\#}$ are defined in such a way that the conservation of impulsion and total energy are satisfied. $\mathcal{M}_{i}$ and $\overline{\mathcal{M}_{i}}(f_{e} , f_{i})$ are given by (14, 15).

In the present case the definitions of $j$ and $\overline{\rho}$ given in (20, 21) become

and the electric field is still given by (18, 19). For the sake of conciseness, we continue the presentation only for a two polyatomic species mixture. However, the following steps can be easily adapted when one of the components is monoatomic.

2.2.3.   Properties of the model

Proposition 1. The model (13, 14, 15, 16, 17) conservesthe mass per species, the total impulsion and the total energy.

The proof is straighforward and based as in [3] on the definition of the fictitious quantities (16, 17).

Theorem 2.1. The model (13, 14, 15, 16, 17) satisfiesan H theorem. The model satisfies the entropic inequality

The equality holds in the above equation if and only if there exists($n_e,n_i, \, u, \,T$) $\in$ $\mathbb{R}^2_+ \times \mathbb{R}^{3}\times \mathbb{R}_+$ such that

An important feature of the polyatomic model (13, 14, 15, 16, 17) is that it satisfies an entropy dissipation property that is compatible with the macroscopic one. The entropy dissipation property has already been obtained in [3] for the system (1), starting directly from the fluid system. In the present case, we are able to show that the entropy of the system (1) is compatible with the Boltzmann entropy, see subsection 3.3.

3.   Construction of the fluid model

3.1.   Scaling on the one dimensional BGK model

Suppose that the system (13, 14, 15, 16, 17) is space homogeneous in two directions. So we assume that the system is even in $v_2$ and $v_3$. Hence, the distribution function $f_{\beta}$ of species $\beta$ depends on time $t \in \mathbb{R}_+$, space $x \in \mathbb{R}$, velocity $v = (v_1 , v_2, v_3)$ and on the energy variable $I \in \mathbb{R}_+$. The model (13, 14, 15, 16, 17) can be rewritten in this case

where $\varepsilon$ is a nonnegative parameter proportional to the Knudsen number. In that case, Maxwellian distributions (14, 15) become

where $u^{\#}$ and $T^{\#}$ are defined in (16) and (17).

3.2.   Hydrodynamic limit

Proposition 2. The system (27, 28, 29, 16, 17) convergesformally to the nonconservative bitemperature Euler system where $E$ is given by generalized Ohm's law

and

Proof. Performing a Chapman-Enskog expansion, it comes that each component of the species is expanded as

so that

Moreover, Gauss equation in (27) implies that $n_i q_i = -n_e q_e + \mathcal{O}(\varepsilon^2)$. So, $n_e = Z n_i + \mathcal{O}(\varepsilon^2)$. Ampère equation then gives $u_e = u_i+\mathcal{O}(\varepsilon^2) = u$. Plugging (32) into the first equation of (27) leads to

for $\{ \beta; \, \delta \} \in \{e;i\} , \; \; \beta \not = \delta$.

Mass conservation equation is obtained by integrating w.r.t $v$ and $I$. The equation of the conservation of the impulsion is then recovered by taking the second moment of (35). Moreover, proceeding as in [3], we get Ohm's law (30). Next the energy equation on the electrons writes

Moreover a direct computation gives

So, according to the relation (10) defining $T$, we get

and $\nu_{ei}$ is given by (31).

3.3.   Entropy dissipation

As proved in [3], the system (1) owns a dissipative entropy-entropy flux pair

where

Here, $C$ is a nonegative constant. With the same method as in [3] we can prove

Theorem 3.1. Let $(\eta, Q)$ be defined by (36)-(37). $\eta$ is a strictly convex dissipative entropy for system (1) and $Q$ is the related entropy flux. More precisely, any smooth solution of the system satisfies the following equality:

If $\mathcal U$ is a weak solution of system (1) obtained as the hydrodynamic limit of the kinetic model, then it satisfies the following inequality

As in [3], the inequality (39) is obtained by using the proof of the H theorem (Theorem 2.1).

4.   Some exact solutions

We compute contact discontinuities and rarefaction waves for system (1) in the homogeneous case $\nu_{ei} = 0$. As $u$ is a double eigenvalue we have contact discontinuities of two kinds. Moreover the computation of the rarefaction waves is different from the one for the classical Euler system.

4.1.   Contact discontinuities

A contact discontinuity is a weak solution $\mathcal{U} = (\rho,\rho u, \mathcal{E}_e, \mathcal{E}_i)$ of (1) such that $u$ is a constant and

Here, $\rho_L$, $\rho_R$, $\mathcal{E}_{\beta,L}$, $\mathcal{E}_{\beta,R}$ are constant. In that case the homogeneous system related to (1) can be written under the conservative form:

Rankine-Hugoniot jump conditions are

One can realize those conditions by taking $[\rho] = 0$ or not. In the first case, nontrivial contact discontinuities are obtained with nonequal left and right values of the ionic and electronic pressures. This case is specific to the bitemperature model. The case $[\rho]\not = 0$ appears also in contact discontinuities for the $3 \times 3$ monotemperature Euler system.

For contact discontinuities there is no entropy dissipation:

4.2.   Rarefaction waves

A rarefaction wave is a selfsimilar continuous, piecewise $C^1$ solution of (1) with $\nu_{ei} = 0$. As we look for a smooth solution, we may use the variable $\mathcal{V} = (\rho,u,\varepsilon_e,\varepsilon_i)$. The rarefaction waves are given by solutions $\mathcal{V}(x,t) = \mathcal{W}\left(\frac{x}{t}\right)$ of the homogeneous system related to (5) that is, denoting $y = x/t$:

where $A$ is given by (6). Rarefaction waves are closely related to the integral curves of the eigenvectors of $A$, as soon as the fields are genuinely nonlinear, see [18] for example. Let us consider the eigenvalue $\lambda_+ = u+a$ with $a$ given by (7), with the eigenvector $r_4$ satisfying (8). We solve the ODS

We set $\mathcal{W}_{R} = \Phi(\xi_0)$ with $\xi_0>0$ and $\Psi(\xi) = \lambda_+(\Phi(\xi))$. $\Psi$ is an increasing monotone function. We set

We have $\mathcal{W} \in C^0( \mathbb{R})$ and (41). Hence we have to solve for $\xi>0$:

We find

As

we have also:

If $\gamma_i \not = \gamma_e$ we cannot parametrize the rarefaction curves by the pressure as one does for the monotemperature Euler system. Hence from (42), the rarefaction curves can be parametrized as follows

So

We retrieve the fact that $\rho_R>\rho_L$. We have also $p_ \beta > p_{ \beta,L}$ and the specific entropy by species is constant. It is defined by

For all $\xi$:

hence

5.   Numerical approximation

In this section we derive a numerical scheme starting from a semi-discretization of the kinetic model. In [3], this approach has been developed for a monoatomic gas mixture. We follow the lines of [3] and obtain the scheme for the polyatomic case. For the sake of completeness we give the details.

First we recall that for $P_ \beta$ defined by

one has

where $F_ \beta$ is the flux of $3 \times 3$ Euler equations with $\gamma_ \beta$ pressure law.

The spatial discretisation is defined by the step $\Delta x$ and the cells $C_j = ]x_{j- \frac{1}{2}},x_{j+ \frac{1}{2}}[$. We consider that $\Delta x$ is constant whereas the time step $\Delta t$: $t_0 = 0$, $t_{n+1} = t_{n}+ \Delta t$ can be variable.

We use a finite volume approach: for any unknown $V(x,t)$, we look for the approximation $V^n_j$ of the average of $V$ at time $t_n$ on the cell $C_j$. Suppose that $\mathcal U^n$ is known ($n \geq 0$).

First step. For $\beta = e,i$ we set

and $U^n_{ \beta,j} = (\rho^n_{ \beta,j},\rho^n_{ \beta,j}u^n_j, \mathcal{E}_{ \beta,j})$. $f_e^n(v,I)$, $f_i^n(v,I)$ are computed by projection on the maxwellian states:

Second step. We solve the transport equations $\partial_t f_\beta + v_1 \partial_x f_\beta = 0$ by using a numerical flux $h_{ \beta,j+ \frac{1}{2}}(v,I) = h_ \beta(f_{ \beta,j}(v,I),f_{ \beta,j+1}(v,I),v,I)$. Here we choose the HLL type flux

where $\lambda_1\leq 0 \leq \lambda_2$ are distinct real constants to be fixed. We set for all $(v,I)$

Remark 1. This scheme is stable for $\lambda_1 \leq v_1 \leq \lambda_2$ and appropriate CFL condition. As we are going to integrate this formula w.r.t. $v_1 \in \mathbb{R}$, we have no chance to prove stability at this level. The resulting macroscopic scheme owns however stability properties that can be proven by a different way, as we shall show at the end of this section.

We apply $P_e$ on $f_e$, $P_i$ on $f_i$, and obtain $U^{n+ \frac{1}{2}}_{e,j}$ and $U^{n+ \frac{1}{2}}_{i,j}$:

Denoting for $\beta = e,i$

we have by using (44), (46) and (47)

and

Hence $U^{n+ \frac{1}{2}}_{ \beta,j}$ is computed by the HLL scheme ($\lambda_1 \leq 0 \leq \lambda_2$).

Remark 2. We could have chosen the upwind scheme to approximate the transport equations, but it is more difficult to integrate the formulas w.r.t. $v_1$ in this case, because the numerical flux depends on the sign of $v_1$.

It is easy to prove the following result.

Proposition 3. For $\beta = e,i$

where $\rho$ is defined in (2).

Third step. We take into account the force terms and the source terms. For all $j \in \mathbb{Z}$, $\alpha, \beta \in \{e,i\}$ and $\beta \not = \alpha$, we define

and

One obtains the following equations for $\beta = e,i$:

Subsequently, it is necessary to ensure that the quasineutrality constraints are satisfied, which correspond to Maxwell-Gauss and Maxwell-Ampère equations in the limit $\varepsilon \rightarrow 0$:

By proposition 3 the first condition is satisfied and $\rho^{n+1}_j = \rho^{n+ \frac{1}{2}}_{j}.$ The second condition is equivalent to $u^{n+1}_{i,j} = u^{n+1}_{e,j} = u^{n+1}_j .$ As a consequence if $\mathcal U^{n+1}_j = (\rho^{n+1}_j , \rho^{n+1}_j u^{n+1}_j , \mathcal{E}^{n+1}_{e,j}, \mathcal{E}^{n+1}_{i,j})$ then $U^{n+1}_{e,j}$ and $U^{n+1}_{i,j}$ satisfy (45), so our notation is consistent. By applying these properties to (52) for $\beta = e, i$, one gets:

Hence, by multiplying the first equation by $c_i$ and the second equation by $c_e$, and then by substracting one to the other, one obtains, analoguously to the continuous case, the discrete generalized Ohm law:

where

Finally, defining $F_{j+ \frac{1}{2}}$ by

we get a scheme that is consistent with the Euler system (1):

Proposition 4. For any $n \geq 0$ if $\mathcal U^n = \{ \mathcal U^n_{j}\}_{j \in \mathbb{Z}}$ is the approximate solution of the system (1) at time $t_n$, we set

The kinetic flux $h_{ \beta}$ is defined by (47). The numerical fluxes $F_{ \beta,j+ \frac{1}{2}}$, $\delta_{j+ \frac{1}{2}}$ and $F_{j+ \frac{1}{2}}$ are then defined by (48), (53), (54). The approximate solution at time $t_{n+1}$ is defined by the implicit scheme

Remark 3. The scheme is implicit but easy to implement because the first two equations of (56) give $\rho^{n+1}_j$ and $u^{n+1}_j$ explicitly and the two last equations can be expressed as a linear $2\times 2$ system for the unknown $T_e$, $T_i$.

This numerical scheme is the same as the one obtained in section 3.2 of [3] and for which a second order two-dimensional extension is presented in [6], ($\lambda_1 \lambda_2 \leq 0$). However, in these two articles, the scheme was obtained by a very different method involving models developed in [7]. Those models are formally like discrete BGK equations but are not based on a physically realistic kinetic representation of the equations. The novelty here is that our polyatomic BGK model gives a physical justification to it, including the case $\gamma_e \not = \gamma_i$.

We recall that stability and entropy properties where proved in these references. In particular, a discrete entropy inequality holds under appropriate choices of $\lambda_1$ and $\lambda_2$ and a CFL condition using the sound velocity of each species. These theoritical conditions give rise however to too much numerical diffusion. We replace them by using the global sound velocity defined in (7) (see [6]):

6.   Numerical results

As pointed out in section 5, the presented scheme has already been tested in [3] and [6]. The aim of this section is to investigate more precisely the polyatomic case with two tests where $\gamma_e = \frac{5}{3}$ and $\gamma_i = \frac{7}{5}$. We make use of the second order extension in space and time with affine reconstruction and Heun scheme developed in [6].

6.1.   Double rarefaction

We solve the bitemperature Euler system with the following Riemann data, whose orders of magnitude are those encountered in physical situations:

We set $Z = 1$ and a final computation time $t = 4.0901 \times 10^{-7}\, sec$. The numerical simulations are performed on the interval $[0,1]$ with 2000 cells. The values of $\lambda_1$ and $\lambda_2$ are computed at every time-step with the condition (57). We set $\nu_{ei} = 0$ so that the solution consists of two rarefaction waves propagating to left and right, the contact discontinuity being trivial. In order to determine the analytical solution, denoting $(\rho,u,\varepsilon_e,\varepsilon_i)$ the intermediate state, one has to find $\xi>0$ such that:

with

Hence we find numerically $\xi >0$ such that

The numerical results are depicted in Figures 1 and 2. We compare the exact and computed results for density, velocity and temperatures. As already observed in [3] for $\gamma_e = \gamma_i = 5/3$, a peak of ionic temperature happens at $x = 1/2$. This peak is similar as the one observed classically for the monotemperature $3 \times 3$ Euler system.

6.2.   Two dimensional implosion

In this test case, we consider an implosion problem, introduced in [28] in the monoatomic case $\gamma_e = \gamma_i = \frac{5}{3}$. We compared our approach with the conservative one of this paper in [6]. Here we set $\gamma_e = \frac{5}{3}$ and $\gamma_i = \frac{7}{5}$ and keep the other parameters unchanged, that is: the physical domain is the square $[-1,1]\times[-1,1]$. Initial data for this Riemann problem is as follows: $\rho = 1$ kg.m$^{-3}$, $u = 0$ m.s$^{-1}$ and temperatures are given by:

The relaxation frequency $\nu_{ei}$ is chosen realistically, according to the formulae given by the NRL formulary [33]. Due to stiffness of the source term one has $T_e = T_i$ very quickly.

Thanks to symmetry properties of the problem, it is only necessary to solve it on the quarter domain $[0,1]\times[0,1]$, equipped with suitable boundary conditions.

On Figure 3, are given the isovalues of the total density and of the electronic temperature at time $t = 4.0901\times 10^{-7}$ sec. We can see that the qualitative behaviour is the same as in the monoatomic case, see Figure 7 of [6]. Here also one can observe a shock propagating towards the center of the domain, a contact discontinuity, and a rarefaction wave propagating to the exterior. However the values of densities and temperatures are different.

When the shock front reaches the center, a peak of density occurs. This peak occurs at time $t = 8.798 \times 10^{-7}$ sec. in the monoatomic case, while it occurs here at time $t = 9.2 \times 10^{-7}$ sec., see Figure 4. The temperatures are also maximal at the center when this peak occurs, as shown in Figure 5.