BGK and Fokker-Planck models of the Boltzmann equation for gases with discrete levels of vibrational energy

TL;DR

This paper introduces BGK and Fokker-Planck models for the Boltzmann equation tailored to diatomic gases with discrete vibrational energy levels, crucial for high-temperature flow simulations such as atmospheric re-entry.

Contribution

It develops and analyzes new kinetic models that incorporate discrete vibrational energy modes, ensuring conservation, entropy properties, and deriving Navier-Stokes asymptotics.

Findings

Models satisfy conservation laws

Models adhere to H-theorem (entropy increase)

Derived Navier-Stokes asymptotics for these models

Abstract

We propose two models of the Boltzmann equation (BGK and Fokker-Planck models) for rarefied flows of diatomic gases in vibrational non-equilibrium. These models take into account the discrete repartition of vibration energy modes, which is required for high temperature flows, like for atmospheric re-entry problems. We prove that these models satisfy conservation and entropy properties (H-theorem), and we derive their corresponding compressible Navier-Stokes asymptotics.

Equations

$$M_{vib}[f](v,\epsilon ,i)=\frac{\rho }{{\sqrt{2\pi RT}}^3}\frac{1-e^{-T_0/T}}{RT}\exp \left(-\frac{\frac{1}{2}|u-v{|}^2+\epsilon +iR{T}_0}{RT}\right).$$

$$\rho ={⟨f⟩}_{v,\epsilon ,i},\rho u={⟨vf⟩}_{v,\epsilon ,i},\rho e={⟨\left(\frac{1}{2}|v-u{|}^2+\epsilon +iR{T}_0\right)f⟩}_{v,\epsilon ,i},$$

$${⟨M_{vib}[f]⟩}_{v,\epsilon ,i}=\rho ,{⟨v{M}_{vib}[f]⟩}_{v,\epsilon ,i}=\rho u.$$

$$\rho e_{tr}(T)={⟨(\frac{1}{2}(v-u{)}^2)M_{vib}[f]⟩}_{v,\epsilon ,i}=\frac{3}{2}\rho RT,$$

$$\rho e_{rot}(T)={⟨\epsilon M_{vib}[f]⟩}_{v,\epsilon ,i}=\rho RT,$$

$$\rho e_{vib}(T)={⟨(iR{T}_0)M_{vib}[f]⟩}_{v,\epsilon ,i}=\rho \frac{R{T}_0}{e^{T_0/T}-1}=\frac{\delta (T)}{2}\rho RT,$$

$$\delta (T)=\frac{2T_0/T}{e^{T_0/T}-1},$$

$${⟨(\frac{1}{2}(v-u{)}^2+\epsilon +iR{T}_0)M_{vib}[f]⟩}_{v,\epsilon ,i}=\rho e,$$

$$e=\frac{5+\delta (T)}{2}RT.$$

$$F(t,x,v,\epsilon )=\sum _{i}f(t,x,v,\epsilon ,i),\text{ and }G(t,x,v,\epsilon )=\sum _{i}iR{T}_{0}f(t,x,v,\epsilon ,i).$$

$$\rho ={⟨F⟩}_{v,\epsilon },\rho u={⟨vF⟩}_{v,\epsilon },\rho e={⟨(\frac{1}{2}(v-u{)}^2+\epsilon )F⟩}_{v,\epsilon }+{⟨G⟩}_{v,\epsilon }.$$

$$\begin{array}{rl}& \mathcal{H}(F,G)={⟨H(F,G)⟩}_{v,\epsilon },\text{ where }\\ & H(F,G)=\underset{f>0}{\inf }\left\{\sum _i{f}_i\log f_i\text{such that}\sum _i{f}_i=F,\sum _{i}iR{T}_0{f}_i=G\right\}.\end{array}$$

$$\mathcal{H}(F,G)={⟨F\log (F)+F\log \left(\frac{R{T}_{0}F}{R{T}_{0}F+G}\right)+\frac{G}{R{T}_0}\log \left(\frac{G}{R{T}_{0}F+G}\right)⟩}_{v,\epsilon }.$$

$$\mathcal{L}(f,\alpha ,\beta )=\sum _i{f}_i\log f_i-\alpha \left(\sum _i{f}_i-F\right)-\beta \left(\sum _{i}iR{T}_0{f}_i-G\right),$$

$$F=\sum _i{f}_i=\frac{A}{1-\exp \left(-BR{T}_0\right)},$$

$$G=\sum _{i}iR{T}_0{f}_i=\frac{AR{T}_0\exp \left(-BR{T}_0\right)}{{\left(1-\exp \left(-BR{T}_0\right)\right)}^2},$$

$$A=\frac{R{T}_0{F}^2}{R{T}_{0}F+G},B=\frac{1}{R{T}_0}\log \left(1+\frac{R{T}_{0}F}{G}\right).$$

$$H(F,G)=F\log (F)+F\log \left(\frac{R{T}_{0}F}{R{T}_{0}F+G}\right)+\frac{G}{R{T}_0}\log \left(\frac{G}{R{T}_{0}F+G}\right).$$

$$D_{1}H(F,G)=1+\log \left(\frac{R{T}_0{F}^2}{R{T}_{0}F+G}\right),D_{2}H(F,G)=\frac{1}{R{T}_0}\log \left(\frac{G}{R{T}_{0}F+G}\right).$$

$$D_{11}H(F,G)=\frac{2}{F}-\frac{R{T}_0}{R{T}_{0}F+G},$$

$$D_{21}H(F,G)=D_{12}H(F,G),$$

$$\begin{array}{rl}& F{D}_{11}H(F,G)+G{D}_{21}H(F,G)=1,\\ & F{D}_{12}H(F,G)+G{D}_{22}H(F,G)=0.\end{array}$$

$$\mathcal{S}=\left\{(F_1,G_1)\text{ such that}{⟨F_1⟩}_{v,\epsilon }=\rho ,{⟨v{F}_1⟩}_{v,\epsilon }=\rho u,{⟨(\frac{1}{2}|v{|}^2+\epsilon )F_1+G_1⟩}_{v,\epsilon }=\rho e\right\}.$$

$$M_{vib}[F,G]=\frac{\rho }{{\sqrt{2\pi RT}}^3}\exp \left(-\frac{|v-u{|}^2}{2RT}\right)\frac{1}{RT}\exp \left(-\frac{\epsilon }{RT}\right)$$

$$\begin{array}{rl}& D_{1}H(F_1,G_1)(M_{vib}[F,G]-F_1)+D_{2}H(F_1,G_1)(e_{vib}(T)M_{vib}[F,G]-G_1)\\ & \le H(M_{vib}[F,G],e_{vib}(T)M_{vib}[F,G])-H(F_1,G_1)\le 0.\end{array}$$

$$\begin{array}{rl}\mathcal{L}(F_1,G_1,\alpha ,\beta ,\gamma )=& {⟨H(F_1,G_1)⟩}_{v,\epsilon }-\alpha ({⟨F_1⟩}_{v,\epsilon }-\rho )\\ & -\beta \cdot ({⟨v{F}_1⟩}_{v,\epsilon }-\rho u)-\gamma \left({⟨(\frac{1}{2}|v{|}^2+\epsilon )F_1+G_1⟩}_{v,\epsilon }-\rho e\right)\end{array}$$

$$D_{1}H(F_1,G_1)=\alpha +\beta \cdot v+\gamma \frac{1}{2}|v{|}^2,$$

$$D_{2}H(F_1,G_1)=\gamma ,$$

$${⟨F_1⟩}_{v,\epsilon }-\rho =0,$$

$${⟨v{F}_1⟩}_{v,\epsilon }-\rho u=0,$$

$${⟨(\frac{1}{2}|v{|}^2+\epsilon )F_1+G_1⟩}_{v,\epsilon }-\rho e=0,$$

$$\begin{array}{rl}& {\partial }_{F_1}\mathcal{L}(F_1,G_1,\alpha ,\beta ,\gamma )(\delta F_1)\\ & =⟨(D_{1}H(F_1,G_1)-(\alpha +\beta \cdot v+\gamma (\frac{1}{2}|v{|}^2+\epsilon )))\delta F_1{⟩}_{v,\epsilon },\end{array}$$

$$F_1$$

$$G_1$$

$${\partial }_{t}f+v\cdot \nabla f=\frac{1}{\tau }(M_{vib}[f]-f).$$

Full text

Keywords: Fokker-Planck model, BGK model, H-theorem, Rarefied Gas Dynamics, vibrational molecules

1 Introduction

Numerical simulation of atmospheric reentry flows requires to solve the Boltzmann equation of Rarefied Gas Dynamics. The standard method to do so is the Direct Simulation Monte Carlo (DSMC) method [1, 2], which is a particle stochastic method. However, it is sometimes interesting to have alternative numerical methods, like, for instance, methods based on a direct discretization of the Boltzmann equation (see [3]). This is hardly possible for the full Boltzmann equation (except for monatomic gases, see [4]), since this is still much too computationally expensive for real gases. But BGK like model equations [5] are very well suited for such deterministic codes: indeed, their complexity can be reduced by the well known reduced distribution technique [6], which leads to intermediate models between the full Boltzmann equation and moment models [7]. The Fokker-Planck model [8] is another model Boltzmann equation that can give very efficient stochastic particle methods, see [9].

These model equations have already been extended to polyatomic gases, so that they can take into account the internal energy of rotation of gas molecules. They contains correction terms that lead to correct transport coefficients: the ESBGK or Shakhov’s models [10, 11, 12], and the cubic Fokker-Planck and ES-Fokker-Planck [9, 13, 14, 15].

For high temperature flows, like in space reentry problems, the vibrational energy of molecules is activated, and has a significant influence on energy transfers in the gas flow. It is therefore interesting to extend the model equations to take this vibrational modes into account. Several extended BGK models have been recently proposed to do so, for instance [16, 17, 18, 19], and a recent Fokker-Planck model has been proposed earlier in [13].

All these models assume a continuous vibrational energy repartition. However, while transitional and rotational energies in air can be considered as continuous for temperature larger than $1$K and $10$K, respectively, vibrational energy can be considered as continuous only for much larger temperatures ($2000$K for oxygen and $3300$K for nitrogen). For flows up to $3000$K around reentry vehicles, the discrete levels of vibrational energy must be used [20]. It seems that that the only BGK model that allows for this discrete repartition is the model of Morse [21].

In this paper, we consider a simpler version of this Morse BGK model for vibrating gases that allows for a discrete vibrational energy. We show that the complexity of this model can be reduced with the reduced distribution technique so that the discrete vibrational energy is eliminated. What is new here is that this construction allows us to prove that the corresponding reduced model satisfies the H-theorem. Moreover, the model is shown to give macroscopic Euler and Navier-Stokes equations in the dense regime, with temperature dependent heat capacities, as expected. This means that the reduced model is a good candidate for its implementation in a deterministic simulation code. Note that with this reduction, only higher order moments with respect to the vibration energy variable are lost: the macroscopic quantities of interest like pressure, temperature, and heat flux, are the same as in the non-reduced model. Moreover, since the reduced variable is not the velocity, this reduction does not require any assumption or special geometries.

An equivalent reduced Fokker-Planck model is also proposed, that has the same properties. However, this model is not based on a non-reduced model, since we are not able so far to define diffusion process for the discrete vibrational energy. Up to our knowledge, this is the first time such a Fokker-Planck model for vibration energy is proposed.

Our paper is organized as follows. In section 2, we present the kinetic description of a gas with vibrating molecules, and we discuss the mathematical properties of the reduced distributions that will be used for our models. Our BGK and Fokker-Planck models are presented in sections 3 and 4, respectively. In section 5, the hydrodynamic limits of our models, obtained by a Chapman-Enskog procedure, are discussed. Section 6 gives some perspectives of this work.

2 Kinetic description of a vibrating diatomic gas

2.1 Distribution function and local equilibrium

We consider a diatomic gas. We define $f(t,x,v,\varepsilon,i)$ the mass density of molecules with position $x$, velocity $v$, internal energy $\varepsilon$, and in the $i$th vibrational energy level, such that the corresponding vibrational energy is $iRT_{0}$, where $T_{0}=h\nu/k$ is a characteristic vibrational temperature of the molecule ($h$ and $k$ are the Planck and Boltzmann constant, while $\nu$ is the fundamental vibrational frequency of the molecule).

The corresponding local equilibrium distribution is defined by (see [1])

[TABLE]

Here, $\rho$ is the mass density of the gas, $T$ its temperature of equilibrium and $u$ its mean velocity, defined below.

2.2 Moments and entropy

The macroscopic quantities are defined by moments of $f$ as follows:

[TABLE]

where we use the notation $\left\langle\phi\right\rangle_{v,\varepsilon,i}=\sum_{i=0}^{\infty}\iint\phi(t,x,v,\varepsilon,i)\,dvd\varepsilon$ for any function $\phi$.

With standard Gaussian integrals and summation formula, it is easy to find that the moments of the equilibrium $M_{vib}[f]$ satisfy:

[TABLE]

At equilibrium, we can define the following energies of translation, rotation, and vibration

[TABLE]

where the number of degrees of freedom of vibrations is

[TABLE]

which is a non integer and temperature dependent number, while the number of degrees of freedom is $3$ for translation and $2$ for rotation.

The temperature $T$ is defined so that $M_{vib}[f]$ has the same energy as $f$:

[TABLE]

which gives the non linear implicit definition of $T$:

[TABLE]

Since the function $T\rightarrow e$ is monotonic, $T$ is uniquely defined by (7). Moreover, note that $\delta(T)$ is necessarily between 0 and 2, which means that vibrations add at most two degrees of freedom.

Finally, the entropy ${\cal H}(f)$ of $f$ is defined by ${\cal H}(f)=\left\langle f\log f\right\rangle_{v,\varepsilon,i}.$

2.3 Reduced distributions

For computational efficiency, it is interesting to define marginal, or reduced, distributions $F$ and $G$ by

[TABLE]

The macroscopic variables defined by $f$ can be obtained through $F$ and $G$ only, as it is shown in the following proposition by integrating with respect to $v$ and $\varepsilon$ and using the definition (2) of the moments.

Proposition 2.1** (Moments of the reduced distributions).**

The macroscopic variables $\rho$, $u$, and $e$, of $f$, defined by (2), satisfy

[TABLE]

where we use the notation $\left\langle\psi\right\rangle_{v,\varepsilon}=\iint\psi(t,x,v,\varepsilon)\,dvd\varepsilon$ for any function $\psi$.

This reduction procedure can be extended to the entropy functional as follows. First, to simplify the following relations, we use the notation $f_{i}(v,\varepsilon)$ for $f(v,\varepsilon,i)$. Then, we define the reduced entropy by

[TABLE]

In other words, for a given couple of reduced distributions $(F,G)$, we define the (non reduced) distribution that minimizes the marginal entropy $\sum_{i}f_{i}\log f_{i}$ among all the distributions that have the same marginal distributions $F$ and $G$. Then the reduced entropy is the integral with respect to $v$ and $\varepsilon$ of the corresponding marginal entropy.

Now it is possible to represent this reduced entropy as a function of $F$ and $G$ only, as it is shown in the following proposition.

Proposition 2.2** (Entropy).**

The reduced entropy $\mathcal{H}(F,G)$ defined by (9) is

[TABLE]

Proof.

The set $\left\{f>0\text{ such that }\sum_{i}f_{i}=F,\quad\sum_{i}iRT_{0}f_{i}=G\right\}$ is clearly convex, so that we can use a Lagrangian multiplier approach by finding if there exists a minimum of the function $\mathcal{L}$ defined through :

[TABLE]

where $\alpha$ and $\beta$ are real numbers and $\sum_{i}f_{i}\log f_{i}$ is a convex function of $f$ . The functional $\mathcal{L}$ is convex but no longer strictly convex. A minimum of $\mathcal{H}(F,G)$ necessarily satisfies $\frac{\partial\mathcal{L}}{\partial f}=0$, and it is easy to deduce that $f$ can be written $f_{i}(v,\varepsilon)=A\exp\left(-iBRT_{0}\right)$, where $A:=A(v,\varepsilon)$ and $B:=B(v,\varepsilon)$ are functions that are still to be determined.

The linear constraints give:

[TABLE]

where we have used the property $iRT_{0}f_{i}=-\frac{df_{i}}{dB}$ that comes from $f_{i}=A\exp\left(-iBRT_{0}\right)$. Solving this linear system gives

[TABLE]

so that

[TABLE]

A final integration with respect to $v$ and $\varepsilon$ gives the final result. ∎

The following proposition gives useful differential properties of the reduced entropy functional.

Proposition 2.3** (Properties of $H$).**

The partial derivatives of $H$ computed at $(F,G)$ are:

[TABLE] 2. 2.

We denote by $\displaystyle\mathbb{H}=\left(\begin{smallmatrix}D_{11}H(F,G)&\quad D_{12}H(F,G)\\ D_{12}H(F,G)&\quad D_{22}H(F,G)\end{smallmatrix}\right)$ the Hessian matrix of $H$. The second order derivatives are:

[TABLE]

and we have

[TABLE] 3. 3.

The function $(F,G)\mapsto H(F,G)$ is convex.

Proof.

Points 1 and 2 are given by direct calculations. For point 3, note that the determinant of the Hessian matrix $\mathbb{H}$, which is $\det\mathbb{H}=\frac{1}{G(RT_{0}F+G)}$ is positive. Moreover, its trace is positive too, so that the Hessian matrix is positive definite, and hence the function $H$ is convex. ∎

Now, we want to use this reduced entropy to define the corresponding reduced equilibrium. This is done by computing the minimum of the reduced entropy among all the reduced distributions $(F_{1},G_{1})$ that have the same moments as $(F,G)$, as it is stated in the following proposition.

Proposition 2.4** (Reduced equilibrium).**

Let $(F,G)$ be a couple of reduced distributions and $\rho$, $\rho u$, and $\rho e$ its moments as defined by (8). Let ${\cal S}$ be the convex set defined by

[TABLE]

The minimum of $\mathcal{H}$ on ${\cal S}$ is obtained for the couple $(M_{vib}[F,G],e_{vib}(T)M_{vib}[F,G])$ with

[TABLE]

where $e_{vib}(T)$ is the equilibrium vibrational energy defined by (5) and $\rho,u,T$ depend on $F$ and $G$ through the definition of the moments. 2. 2.

For every $(F_{1},G_{1})$ in $\mathcal{S}$, we have

[TABLE]

Proof.

First, the set $\mathcal{S}$ is clearly convex, and it is non empty, since it is easy to see that $(M_{vib}],e_{vib}(T)M_{vib})$ realizes the moments $\rho$, $\rho u$, and $\rho e$, and hence belongs to $\mathcal{S}$. Now, we define the following Lagrangian

[TABLE]

for every positive $(F_{1},G_{1})$, $\alpha\in\mathbb{R}$, $\beta\in\mathbb{R}^{3}$, $\gamma\in\mathbb{R}$. The reduced entropy can reach a minimum of ${\mathcal{S}}$ when $\mathcal{L}$ has its first derivatives equal to zero: it is a minimum if it is unique . Such a point, denoted by $(F_{1},G_{1},\alpha,\beta,\gamma)$ for the moment, is characterised by the fact that the partial derivatives of ${\mathcal{L}}$ vanish at $(F_{1},G_{1},\alpha,\beta,\gamma)$. This gives the following relations (using the cancellation of the $\mathcal{L}$ derivatives in $F_{1},G_{1},\alpha,\beta,\gamma$ respectively)

[TABLE]

where $D_{1}H$ and $D_{2}H$ are defined in (12). For instance first relation comes from the fact that the derivative with respect to $F_{1}$ satisfies for every $\delta F_{1}$

[TABLE]

It is true for all $\delta F_{1}$ leading to the relation 15.

Now Combining equations (15) and (16), one gets that there exists four real numbers $A$, $B$, $C$, $D$ and one vector $E\in\mathbb{R}^{3}$, independent of $v$ and $\varepsilon$, such that:

[TABLE]

where $B$ and $C$ are necessarily non positive to ensure the integrability of $F_{1}$ and $G_{1}$. It is then standard to use equations (17)–(19) to get $F_{1}=M_{vib}(F,G)$ and $G_{1}=e_{vib}(T)M_{vib}(F,G)$.

Finally point 2 is a direct consequence of the convexity of $H$ and of the minimization property.

∎

3 A BGK model with vibrations

With the local equilibrium $M_{vib}[f]$ defined in (1), it is easy to derive the following BGK model:

[TABLE]

The macroscopic parameters $\rho$, $u$, and $T$ are defined through the moments $\rho$, $\rho u$ and $\rho e$ of $f$ (see (2)).

Like in the BGK model for monoatomic gases, it will be shown that the relaxation time of this BGK model is $\tau=\mu/p$, where $p=\rho RT$ is the pressure and $\mu$ the viscosity, that can be temperature dependent.

Now we have the following properties.

Property 3.1**.**

• •

Conservation: for BGK model (20) the mass, momentum and total energy are conserved:

[TABLE]

• •

H-theorem: for the entropy ${\cal H}(f)=\left\langle f\log f\right\rangle_{v,\varepsilon,i}$, we have

[TABLE]

The proof relies on standard arguments (definition of $M_{vib}[f]$ and convexity of $x\log x$) and is left to the reader.

3.1 A reduced BGK model with vibrations

For computational reasons, it is interesting to reduce the complexity of model (20) by using the usual reduced distribution technique [22]. We define the reduced distributions

[TABLE]

and by summation of (20) on $i$ we get the following closed system of two reduced equations:

[TABLE]

where the reduced Maxwellian is

[TABLE]

and the macroscopic quantities are defined by

[TABLE]

and $T$ is still defined by (7) which implies that $T$ depends both on $F$ and $G$: to avoid the heavy notation $T[F,G]$, it will still be denoted by $T$ in the following.

Note that this model can easily be reduced once again to eliminate the rotational energy variable. This gives a reduced system of three BGK equations, with three distributions.

It is interesting to compare our new model to the work of [23] and [19]: in these recent papers, the authors also proposed, independently, BGK and ES-BGK models for temperature dependent $\delta(T)$, like in the case of vibrational energy. However, they are not based on an underlying discrete vibrational energy partition, and the authors are not able to prove any H-theorem. Only a local entropy dissipation can be proved. The advantage of our new approach is that the reduced model, which is continuous in energy too, inherits the entropy property from the non-reduced model, and hence a H-theorem, as it is shown below.

3.2 Properties of the reduced model

System (21–22) naturally satisfies local conservation laws of mass, momentum, and energy. Moreover, the H-theorem holds with the reduced entropy $H(F,G)$ as defined in (9). Indeed, we have the

Proposition 3.1**.**

The reduced BGK system (21–22) satisfies the H-theorem

[TABLE]

where ${\cal H}(F,G)$ is the reduced entropy defined in (9).

Proof.

By differentiation we get

[TABLE]

where we have used (21–22) to replace the transport terms by relaxation ones, and point 5 of proposition 2.4 to obtain the inequality. ∎

4 A Fokker-Planck model with vibrations

It is difficult to derive a Fokker-Planck model for the distribution function $f$ with discrete energy levels. We find it easier to directly derive a reduced model, by analogy with the reduced BGK model (21–22) and by using our previous work [15] on a Fokker-Planck model for polyatomic gases. We remind that the original Fokker-Planck model for monoatomic gas can be derived from the Boltzmann collision operator under the assumption of small velocity changes through collisions and additional equilibrium assumptions (see [8]). In practice, the agreement of this model with the Boltzmann equation is observed even when the gas is far from equilibrium (see [9], for instance).

4.1 A reduced Fokker-Planck model with vibrations

First, we remind the Fokker-Planck model for a diatomic gas (without vibrations) obtained in [15]:

[TABLE]

where $f=f(t,x,v,\varepsilon)$ and the collision operator is

[TABLE]

where the macroscopic values are

[TABLE]

The internal energy $\varepsilon$ can be eliminated by the reduction technique (integration w.r.t $d\varepsilon$ and $\varepsilon d\varepsilon$) to get

[TABLE]

with the collision operators

[TABLE]

Note that the two velocity drift-diffusion terms in the two previous equations have exactly the same structure as the one in the non-reduced model (24). However, it is interesting to note that the energy drift-diffusion term of (24) gives, after reduction, a relaxation operator in the ${\cal G}$ equation. Moreover by reducing the model we lose some moments of initial distribution functions (higher moments in internal energy notably) but we are still able to capture energies and fluxes which are generally the main quantities of interest.

By analogy, now we propose the following reduced Fokker-Planck model for a diatomic gas with vibrations. Note that now, the model is still with variables $x$, $v$, and $\varepsilon$: only the discrete energy levels $i$ are eliminated. This model is

[TABLE]

with

[TABLE]

where the macroscopic values are defined as in (23) and  (7). Again, note that the temperature $T$ depends on $F$ and $G$.

Note that we do not derive this reduced Fokker-Planck model directly from a model with discrete vibrational energy as for the BGK model, since we are not able so far to define a discrete diffusion operator. As mentioned above, this model is obtained by analogy with the Fokker-Plank model proposed for polyatomic gases. Its derivation from reduction of a discrete in energy Fokker-Plank model will be studied in a future work.

4.2 Properties of the reduced model

Using direct calculations and dissipation properties as in [15] we can prove the following proposition.

Proposition 4.1**.**

The collision operator conserves the mass, momentum, and energy:

[TABLE]

the reduced entropy $\mathcal{H}(F,G)$ satisfies the H-theorem:

[TABLE]

and we have the equilibrium property

[TABLE]

Proof.

The conservation property is the consequence of direct integration of (27). The equilibrium property can be proved as follows. First, note that the Maxwellian $M_{vib}[F,G]$ can be written as

[TABLE]

with $\Omega=\begin{pmatrix}RT&0\\ 0&2\varepsilon RT\end{pmatrix}$. To shorten the notations, $M_{vib}[F,G]$ will be simply denoted by $M_{vib}$ below, and $e_{vib}(T)$ will be simply denoted by $e_{vib}$ as well. Then the collision operators can be written in the compact form

[TABLE]

Then an integration by part gives the following identity for $D_{F}(F,G)$:

[TABLE]

Consequently, if $D_{F}(F,G)=0$, since the integrand in the previous relation is a definite positive form, the gradient is necessarily zero, and hence $F=M_{vib}$. For the equilibrium property of $G$, the proof is a bit more complicated. First, we have

[TABLE]

Consequently, if $D_{G}(F,G)=0$, and since $F=M_{vib}$, we have

[TABLE]

which comes from (8) and $F=M_{vib}$. Therefore, we obtain

[TABLE]

and again this gives $G=e_{vib}M_{vib}$, which concludes the proof of the equilibrium property.

The proof of the H-theorem is much longer. First, by differentiation one gets that the quantity $\mathcal{D}(F,G)=\partial_{t}\mathcal{H}(F,G)+\nabla_{x}\cdot\left\langle vH(F,G)\right\rangle_{v,\varepsilon}$ satisfies:

[TABLE]

from (21–22). Then the proof is based on the convexity of $H(F,G)$: while for the BGK model we only used the first derivatives of $H$, we now use the positive-definiteness of the Hessian matrix of $H$. To do so we integrate by parts $\mathcal{D}(F,G)$ and multiply by $\tau$ so that:

[TABLE]

To use the positive definiteness of the Hessian matrix $\mathbb{H}$ of $H$, we introduce the following vectors:

[TABLE]

and we decompose the partial derivatives of $F$ and $G$ in factor of $D_{11}F$, $D_{22}F$, $D_{12}F$ as follows:

[TABLE]

This gives

[TABLE]

Now this expression can be considerably simplified by using property (13), and we get

[TABLE]

Then the first two terms are simplified by using an integration by parts and relations (8) and (7) to get

[TABLE]

The terms with the Hessian are clearly negative, since $\mathbb{H}$ is positive definite. Then we have

[TABLE]

Note that from (8) the first term can be written as

[TABLE]

and can be factorized with the second term to find

[TABLE]

We can now prove that the integrand of the right-hand side is non-positive. Indeed, assume for instance that the first factor is non-positive, that is to say $e_{vib}(T)F-G\leq 0$. By using $e_{vib}(T)=\frac{RT_{0}}{e^{T_{0}/T}-1}$ (see definition (5)), it is now very easy to prove the following relations

[TABLE]

that is to say the second factor of the integrand is non-negative.

Consequently, we have proved $\tau\mathcal{D}(F,G)\leq 0$, which concludes the proof. ∎

5 Hydrodynamic limits for reduced models

With a convenient nondimensionalization, the relaxation time $\tau$ of the reduced BGK model (21)–(22) and the Fokker-Planck model (25)-(26) is replaced by $\mathrm{Kn}\thinspace\tau$, where $\mathrm{Kn}\thinspace$ is the Knudsen number, which can be defined as a ratio between the mean free path and a macroscopic length scale. It is then possible to look for macroscopic models derived from BGK and Fokker-Planck reduced models, in the asymptotic limit of small Knudsen numbers.

For convenience, these models are re-written below in non-dimensional form. The BGK model is

[TABLE]

where $M_{vib}[F,G]$ can be defined by (14) with $R=1$. Similarly, the relations (3)–(7) between the translational, internal, and total energies and the temperature, have to be read with $R=1$ in non-dimensional variables. We remind that $T$ is still a function of $F$ and $G$. The Fokker-Planck model is

[TABLE]

with

[TABLE]

5.1 Euler limit

In this section, we prove the following proposition:

Proposition 5.1**.**

The mass, momentum, and energy densities $(\rho,\rho u,E=\frac{1}{2}\rho u^{2}+\rho e)$ of the solutions of the reduced BGK and the Fokker-Planck models satisfy the equations

[TABLE]

which are the Euler equations, up to $O(\mathrm{Kn}\thinspace)$. The non-conservative form of these equations is

[TABLE]

where $c_{v}(T)=\frac{d}{dT}e(T)$ is the specific heat capacity at constant volume.

Proof.

The reduced BGK model (21)–(22) is multiplied by $1$, $v$, and $\frac{1}{2}|v|^{2}+\varepsilon$ and integrated with respect to $v$ and $\varepsilon$, which gives the following conservation laws:

[TABLE]

where $\sigma(F)=\left\langle(v-u)\otimes(v-u)F\right\rangle_{v,\varepsilon}$ is the stress tensor, and $q(F,G)=\left\langle(v-u)(\frac{1}{2}|v-u|^{2}+\varepsilon)F\right\rangle_{v,\varepsilon}+\left\langle(v-u)G\right\rangle_{v,\varepsilon}$ is the heat flux.

When $\mathrm{Kn}\thinspace$ is very small, if all the time and space derivatives of $F$ and $G$ are $O(1)$ with respect to $\mathrm{Kn}\thinspace$, then (29)–(30) imply $F=M_{vib}[F,G]+O(\mathrm{Kn}\thinspace)$ and $G=e_{vib}(T)M_{vib}[F,G]+O(\mathrm{Kn}\thinspace)$. Then it is easy to find that $\sigma(F)=\sigma(M_{vib}[F,G])+O(\mathrm{Kn}\thinspace)=pI+O(\mathrm{Kn}\thinspace)$ , where $I$ is the unit tensor, and $q(F,G)=q(M_{vib}[F,G],e_{vib}(T)M_{vib}[F,G])+O(\mathrm{Kn}\thinspace)=O(\mathrm{Kn}\thinspace)$, since the heat flux is zero at equilibrium, which gives the Euler equations (35). The same analysis can be applied for the reduced Fokker-Planck model (31)–(33).

Finally, the non conservative form is readily obtained from the conservative form. Note another formulation of the energy equation that will be useful below:

[TABLE]

where $e_{vib}^{\prime}(T)=\frac{d}{dT}e_{vib}(T)$. ∎

5.2 Compressible Navier-Stokes limit

In this section, we shall prove the following proposition:

Proposition 5.2**.**

The moments of the solution of the BGK and Fokker-Planck kinetic models (21)-(22) and (25)-(26) satisfy, up to $O(\mathrm{Kn}\thinspace^{2})$, the Navier-Stokes equations

[TABLE]

where the shear stress tensor and the heat flux are given by

[TABLE]

and where the following values of the viscosity and heat transfer coefficients (in dimensional variables) are

[TABLE]

while the volumic viscosity coefficient is $\alpha=\frac{c_{p}(T)}{c_{v}(T)}-1$ for both models, and $c_{p}(T)=\frac{d}{dT}(e(T)+p/\rho)=c_{v}(T)+R$ is the specific heat capacity at constant pressure. Moreover, the corresponding Prandtl number is

[TABLE]

Note that both models do not provide a correct Prandtl number, which can lead to errors for the computation of fluxes in numerical simulations. This is a usual problem with single parameter models like BGK or Fokker-Planck: this can be corrected by a modification of the models like with the ES-BGK or ES-FP approaches, as it has been done for polyatomic gases (see [11, 15] for instance).

5.2.1 Proof for the BGK model

The usual Chapman-Enskog method is applied as follows. We decompose $F$ and $G$ as $F=M_{vib}[F,G]+\mathrm{Kn}\thinspace F_{1}$ and $G=e_{vib}(T)M_{vib}[F,G]+\mathrm{Kn}\thinspace G_{1}$, which gives

[TABLE]

Then we have to approximate $\sigma(F_{1})$ and $q(F_{1},G_{1})$ up to $O(\mathrm{Kn}\thinspace)$. This is done by using the previous expansions and (21) and (22) to get

[TABLE]

This gives the following approximations

[TABLE]

and

[TABLE]

Now it is standard to write $\partial_{t}M_{vib}[F,G]$ and $\nabla_{x}M_{vib}[F,G]$ as functions of derivatives of $\rho$, $u$, and $T$, and then to use Euler equations (34) to write time derivatives as functions of the space derivatives only. After some algebra, we get

[TABLE]

where

[TABLE]

Then we introduce (44) into (42) to get

[TABLE]

where we have used the change of variables $(v,\varepsilon)\mapsto(V,J)$ in the integral (the term with $A$ vanishes due to the parity of $M_{0}$). Then standard Gaussian integrals (see appendix A) give

[TABLE]

with $\mu=\tau\rho T$ and $\alpha=\frac{c_{p}}{c_{v}}-1$, which is the announced result, in a non-dimensional form.

For the heat flux, we use the same technique. First for $e_{vib}(T)M_{vib}[F,G]$ we obtain

[TABLE]

where

[TABLE]

Then $q(F_{1},G_{1})$ as given in (43) can be reduced to

[TABLE]

Using again Gaussian integrals , we get

[TABLE]

where $\kappa=\mu c_{p}(T)$ with $c_{p}(T)=\frac{d}{dT}(e(T)+\frac{p}{\rho})=\frac{7}{2}+e_{vib}^{\prime}(T)=1+c_{v}(T)$ in a non-dimensional form.

5.2.2 Proof for the Fokker-Planck model

Here, we rather use the decomposition $F=M_{vib}(1+\mathrm{Kn}\thinspace F_{1})$ and $G=e_{vib}M_{vib}(1+\mathrm{Kn}\thinspace G_{1})$, which gives

[TABLE]

in which, for clarity, the dependence of $M_{vib}$ on $F$ and $G$ has been omitted, and the dependence of $e_{vib}$ on $T$ as well. Finding $F_{1}$ and $G_{1}$ is less simple than for the BGK model: however, the computations are very close to what is done in the standard monatomic Fokker-Planck model (see [14] for instance), so that we only give the main steps here (see appendix A for details).

First, the decomposition is injected into (33) to get

[TABLE]

where $L_{F}$ and $L_{G}$ are linear operators defined by

[TABLE]

Then the Fokker-Planck equations (31)-(32) suggest to look for an approximation of $F_{1}$ and $G_{1}$ up to $O(\mathrm{Kn}\thinspace)$ as solutions of

[TABLE]

By using (44)-(45), these relations are equivalent, up to another $O(\mathrm{Kn}\thinspace)$ approximation, to

[TABLE]

where $A$, $B$, $\tilde{A}$, and $\tilde{B}$ are the same as for the BGK equation in the previous section.

Now, we rewrite $L_{F}(F_{1})$ and $L_{G}(F_{1},G_{1})$, defined in (46), by using the change of variables $V=\frac{v-u}{\sqrt{T}}$ and $G=\frac{\varepsilon}{T}$ to get

[TABLE]

Then simple calculation of derivatives show that $A$, $B$, $\tilde{A}$, and $\tilde{B}$ satisfy the following properties

[TABLE]

Therefore, we look for $F_{1}$ and $G_{1}$ as solution of (47) under the following form

[TABLE]

and we find $\tilde{a}=a=-1/3$ and $\tilde{b}=b=1/2$.

Finally, using these relations into $\sigma$ and $q$ and using some Gaussian integrals (see appendix A) give

[TABLE]

where $\alpha=\frac{c_{p}}{c_{v}}-1$, $\mu=\frac{\tau}{2}\rho T$, and $\kappa=\frac{2}{3}\mu c_{p}(T)$, which is the announced result, in a non-dimensional form.

6 Conclusion

In this paper, we have proposed to different models (BGK and Fokker-Planck) of the Boltzmann equation that allow for vibrational energy discrete modes. These models satisfy the conservation and entropy property (H-theorem), and the vibration energy variable can be eliminated by the usual reduced distribution function. The low complexity of the reduced BGK model can make it attractive to be implemented in a deterministic code, while the Fokker-Planck model can be easily simulated with a stochastic method.

Of course, since these models are based on a single time relaxation, they cannot allow for multiple relaxation times scales. This is not physically correct, since it is known that the relaxation times for translational, rotational, and vibrational energies are very different. However, standard procedures can be used to extend our model, like the ellipsoidal-statistical approach, already used to correct the Prandtl number of the BGK model [11] and Fokker-Plank models [14, 15].

Appendix A Gaussian integrals and other summation formulas

In this section, we give some integrals and summation formula that are used in the paper.

First, we remind the definition of the absolute Maxwellian $M_{0}(V)=\frac{1}{(2\pi)^{\frac{3}{2}}}\exp(-\frac{|V|^{2}}{2})$. We denote by $\langle\phi\rangle=\int_{{\mathbb{R}}^{3}}\phi(V)\,dV$ for any function $\phi$. It is standard to derive the following integral relations (see [24], for instance), written with the Einstein notation:

[TABLE]

while all the integrals of odd power of $V$ are zero. Note that the first relation of each line implies the other relations of the same line: these relations are given here to improve the readability of the paper. From the previous Gaussian integrals, it can be shown that for any $3\times 3$ matrix $C$, we have

[TABLE]

Finally, we have also used the following relations:

[TABLE]

and also

[TABLE]