BGK and Fokker-Planck models of the Boltzmann equation for gases with discrete levels of vibrational energy
aa r X i v : . [ phy s i c s . c l a ss - ph ] M a r BGK and Fokker-Planck models of the Boltzmann equation for gases with discretelevels of vibrational energy
J. Mathiaud , , L. Mieussens CEA-CESTA15 avenue des sablières - CS 6000133116 Le Barp Cedex, France ([email protected]) Univ. Bordeaux, Bordeaux INP, CNRS, IMB, UMR 5251, F-33400 Talence, France. ([email protected])
Abstract
We propose two models of the Boltzmann equation (BGK and Fokker-Planck models) forrarefied flows of diatomic gases in vibrational non-equilibrium. These models take into accountthe 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 andentropy properties (H-theorem), and we derive their corresponding compressible Navier-Stokesasymptotics.
Keywords: Fokker-Planck model, BGK model, H-theorem, Rarefied Gas Dynamics, vibrationalmolecules
Numerical simulation of atmospheric reentry flows requires to solve the Boltzmann equation ofRarefied 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 interestingto have alternative numerical methods, like, for instance, methods based on a direct discretizationof the Boltzmann equation (see [3]). This is hardly possible for the full Boltzmann equation (exceptfor 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, theircomplexity can be reduced by the well known reduced distribution technique [6], which leads tointermediate 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 particlemethods, see [9].These model equations have already been extended to polyatomic gases, so that they can takeinto account the internal energy of rotation of gas molecules. They contains correction terms thatlead to correct transport coefficients: the ESBGK or Shakhov’s models [10, 11, 12], and the cubicFokker-Planck and ES-Fokker-Planck [9, 13, 14, 15].For high temperature flows, like in space reentry problems, the vibrational energy of moleculesis activated, and has a significant influence on energy transfers in the gas flow. It is thereforeinteresting to extend the model equations to take this vibrational modes into account. Several1xtended BGK models have been recently proposed to do so, for instance [16, 17, 18, 19], and arecent Fokker-Planck model has been proposed earlier in [13].All these models assume a continuous vibrational energy repartition. However, while transitionaland rotational energies in air can be considered as continuous for temperature larger than Kand K, respectively, vibrational energy can be considered as continuous only for much largertemperatures (
K for oxygen and
K for nitrogen). For flows up to
K around reentryvehicles, the discrete levels of vibrational energy must be used [20]. It seems that that the onlyBGK 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 thatallows for a discrete vibrational energy. We show that the complexity of this model can be reducedwith the reduced distribution technique so that the discrete vibrational energy is eliminated. Whatis new here is that this construction allows us to prove that the corresponding reduced modelsatisfies the H-theorem. Moreover, the model is shown to give macroscopic Euler and Navier-Stokesequations in the dense regime, with temperature dependent heat capacities, as expected. This meansthat the reduced model is a good candidate for its implementation in a deterministic simulationcode. Note that with this reduction, only higher order moments with respect to the vibration energyvariable 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 definediffusion process for the discrete vibrational energy. Up to our knowledge, this is the first time sucha 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 withvibrating molecules, and we discuss the mathematical properties of the reduced distributions thatwill 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-Enskogprocedure, are discussed. Section 6 gives some perspectives of this work.
We consider a diatomic gas. We define f ( t, x, v, ε, i ) the mass density of molecules with position x ,velocity v , internal energy ε , and in the i th vibrational energy level, such that the correspondingvibrational energy is iRT , where T = hν/k is a characteristic vibrational temperature of themolecule ( h and k are the Planck and Boltzmann constant, while ν is the fundamental vibrationalfrequency of the molecule).The corresponding local equilibrium distribution is defined by (see [1]) M vib [ f ]( v, ε, i ) = ρ √ πRT − e − T /T RT exp − | u − v | + ε + iRT RT ! . (1)Here, ρ 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: ρ = h f i v,ε,i , ρu = h vf i v,ε,i , ρe = (cid:28)(cid:18) | v − u | + ε + iRT (cid:19) f (cid:29) v,ε,i , (2)where we use the notation h φ i v,ε,i = P ∞ i =0 RR φ ( t, x, v, ε, i ) dvdε for any function φ .With standard Gaussian integrals and summation formula, it is easy to find that the momentsof the equilibrium M vib [ f ] satisfy: h M vib [ f ] i v,ε,i = ρ, h vM vib [ f ] i v,ε,i = ρu. At equilibrium, we can define the following energies of translation, rotation, and vibration ρe tr ( T ) = (cid:28) ( 12 ( v − u ) ) M vib [ f ] (cid:29) v,ε,i = 32 ρRT, (3) ρe rot ( T ) = h εM vib [ f ] i v,ε,i = ρRT, (4) ρe vib ( T ) = h ( iRT ) M vib [ f ] i v,ε,i = ρ RT e T /T − δ ( T )2 ρRT, (5)where the number of degrees of freedom of vibrations is δ ( T ) = 2 T /Te T /T − , (6)which is a non integer and temperature dependent number, while the number of degrees of freedomis for translation and for rotation.The temperature T is defined so that M vib [ f ] has the same energy as f : (cid:28) ( 12 ( v − u ) + ε + iRT ) M vib [ f ] (cid:29) v,ε,i = ρe, which gives the non linear implicit definition of T : e = 5 + δ ( T )2 RT. (7)Since the function T → e is monotonic, T is uniquely defined by (7). Moreover, note that δ ( T ) isnecessarily between 0 and 2, which means that vibrations add at most two degrees of freedom.Finally, the entropy H ( f ) of f is defined by H ( f ) = h f log f i v,ε,i . For computational efficiency, it is interesting to define marginal, or reduced, distributions F and G by F ( t, x, v, ε ) = X i f ( t, x, v, ε, i ) , and G ( t, x, v, ε ) = X i iRT f ( t, x, v, ε, i ) . The macroscopic variables defined by f can be obtained through F and G only, as it is shown inthe following proposition by integrating with respect to v and ε and using the definition (2) of themoments. 3 roposition 2.1 (Moments of the reduced distributions) . The macroscopic variables ρ , u , and e ,of f , defined by (2) , satisfy ρ = h F i v,ε , ρu = h vF i v,ε , ρe = (cid:28) ( 12 ( v − u ) + ε ) F (cid:29) v,ε + h G i v,ε . (8) where we use the notation h ψ i v,ε = RR ψ ( t, x, v, ε ) dvdε for any function ψ . This reduction procedure can be extended to the entropy functional as follows. First, to simplifythe following relations, we use the notation f i ( v, ε ) for f ( v, ε, i ) . Then, we define the reduced entropyby H ( F, G ) = h H ( F, G ) i v,ε , where H ( F, G ) = inf f> (X i f i log f i such that X i f i = F, X i iRT f i = G ) . (9)In other words, for a given couple of reduced distributions ( F, G ) , we define the (non reduced)distribution that minimizes the marginal entropy P i f i log f i among all the distributions that havethe same marginal distributions F and G . Then the reduced entropy is the integral with respect to v and ε 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 shownin the following proposition. Proposition 2.2 (Entropy) . The reduced entropy H ( F, G ) defined by (9) is H ( F, G ) = (cid:28) F log( F ) + F log (cid:18) RT FRT F + G (cid:19) + GRT log (cid:18) GRT F + G (cid:19)(cid:29) v,ε . (10) Proof.
The set { f > such that P i f i = F, P i iRT f i = G } is clearly convex, so that we can usea Lagrangian multiplier approach by finding if there exists a minimum of the function L definedthrough : L ( f, α, β ) = X i f i log f i − α X i f i − F ! − β X i iRT f i − G ! , where α and β are real numbers and P i f i log f i is a convex function of f . The functional L isconvex but no longer strictly convex. A minimum of H ( F, G ) necessarily satisfies ∂ L ∂f = 0 , andit is easy to deduce that f can be written f i ( v, ε ) = A exp ( − iBRT ) , where A := A ( v, ε ) and B := B ( v, ε ) are functions that are still to be determined.The linear constraints give: F = X i f i = A − exp ( − BRT ) ,G = X i iRT f i = ART exp ( − BRT )(1 − exp ( − BRT )) , iRT f i = − df i dB that comes from f i = A exp ( − iBRT ) . Solvingthis linear system gives A = RT F RT F + G , B = 1 RT log (cid:18) RT FG (cid:19) . so that H ( F, G ) = F log( F ) + F log (cid:18) RT FRT F + G (cid:19) + GRT log (cid:18) GRT F + G (cid:19) . (11)A final integration with respect to v and ε gives the final result.The following proposition gives useful differential properties of the reduced entropy functional. Proposition 2.3 (Properties of H ) .
1. The partial derivatives of H computed at ( F, G ) are: D H ( F, G ) = 1 + log (cid:18) RT F RT F + G (cid:19) , D H ( F, G ) = 1 RT log (cid:18) GRT F + G (cid:19) . (12)
2. We denote by H = (cid:16) D H ( F,G ) D H ( F,G ) D H ( F,G ) D H ( F,G ) (cid:17) the Hessian matrix of H . The second orderderivatives are: D H ( F, G ) = 2 F − RT RT F + G , D H ( F, G ) = − RT F + G ,D H ( F, G ) = D H ( F, G ) , D H ( F, G ) = FG ( RT F + G ) , and we have F D H ( F, G ) + GD H ( F, G ) = 1 ,F D H ( F, G ) + GD H ( F, G ) = 0 . (13)
3. The function ( F, G ) H ( F, G ) is convex.Proof. Points 1 and 2 are given by direct calculations. For point 3, note that the determinant ofthe Hessian matrix H , which is det H = G ( RT 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. Thisis done by computing the minimum of the reduced entropy among all the reduced distributions ( F , G ) 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 ρ , ρu ,and ρe its moments as defined by (8) . Let S be the convex set defined by S = ( ( F , G ) such that h F i v,ε = ρ, h vF i v,ε = ρu, (cid:28) ( 12 | v | + ε ) F + G (cid:29) v,ε = ρe ) . . The minimum of H on S is obtained for the couple ( M vib [ F, G ] , e vib ( T ) M vib [ F, G ]) with M vib [ F, G ] = ρ √ πRT exp (cid:18) − | v − u | RT (cid:19) RT exp (cid:16) − εRT (cid:17) (14) where e vib ( T ) is the equilibrium vibrational energy defined by (5) and ρ, u, T depend on F and G through the definition of the moments.2. For every ( F , G ) in S , we have D H ( F , G )( M vib [ F, G ] − F ) + D H ( F , G )( e vib ( T ) M vib [ F, G ] − G ) ≤ H ( M vib [ F, G ] , e vib ( T ) M vib [ F, G ]) − H ( F , G ) ≤ . Proof.
First, the set 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 ρ , ρu , and ρe , and hence belongs to S . Now, we define the following Lagrangian L ( F , G , α, β, γ ) = h H ( F , G ) i v,ε − α ( h F i v,ε − ρ ) − β · ( h vF i v,ε − ρu ) − γ (cid:28) ( 12 | v | + ε ) F + G (cid:29) v,ε − ρe ! for every positive ( F , G ) , α ∈ R , β ∈ R , γ ∈ R . The reduced entropy can reach a minimum of S when L has its first derivatives equal to zero: it is a minimum if it is unique . Such a point,denoted by ( F , G , α, β, γ ) for the moment, is characterised by the fact that the partial derivativesof L vanish at ( F , G , α, β, γ ) . This gives the following relations (using the cancellation of the L derivatives in F , G , α, β, γ respectively) D H ( F , G ) = α + β · v + γ | v | , (15) D H ( F , G ) = γ, (16) h F i v,ε − ρ = 0 , (17) h vF i v,ε − ρu = 0 , (18) (cid:28) ( 12 | v | + ε ) F + G (cid:29) v,ε − ρe = 0 , (19)where D H and D H are defined in (12). For instance first relation comes from the fact that thederivative with respect to F satisfies for every δF ∂ F L ( F , G , α, β, γ )( δF )= h ( D H ( F , G ) − ( α + β · v + γ ( 12 | v | + ε ))) δF i v,ε , It is true for all δF 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 ∈ R , independent of v and ε , such that: F = A exp (cid:0) E · v + B | v | + Cε (cid:1) ,G = DF , B and C are necessarily non positive to ensure the integrability of F and G . It is thenstandard to use equations (17)–(19) to get F = M vib ( F, G ) and G = e vib ( T ) M vib ( F, G ) .Finally point 2 is a direct consequence of the convexity of H and of the minimization property. With the local equilibrium M vib [ f ] defined in (1), it is easy to derive the following BGK model: ∂ t f + v · ∇ f = 1 τ ( M vib [ f ] − f ) . (20)The macroscopic parameters ρ , u , and T are defined through the moments ρ , ρu and ρe of f (see (2)).Like in the BGK model for monoatomic gases, it will be shown that the relaxation time of thisBGK model is τ = µ/p , where p = ρRT is the pressure and µ the viscosity, that can be temperaturedependent.Now we have the following properties. Property 3.1. • Conservation: for BGK model (20) the mass, momentum and total energyare conserved: ∂ t * v | v | f + v,ε,i + ∇ x · * v v | v | f + v,ε,i = 0 . • H-theorem: for the entropy H ( f ) = h f log f i v,ε,i , we have ∂ t H ( f ) + ∇ x · h vf log f i v,ε,i = 1 τ h ( M vib [ f ] − f ) log f i v,ε,i ≤ . The proof relies on standard arguments (definition of M vib [ f ] and convexity of x log x ) and isleft to the reader. For computational reasons, it is interesting to reduce the complexity of model (20) by using theusual reduced distribution technique [22]. We define the reduced distributions F = X i f ( t, x, v, ε, i ) , and G = X i iRT f ( t, x, v, ε, i ) , and by summation of (20) on i we get the following closed system of two reduced equations: ∂ t F + v · ∇ x F = 1 τ ( M vib [ F, G ] − F ) , (21) ∂ t G + v · ∇ x G = 1 τ ( δ ( T )2 RT M vib [ F, G ] − G ) , (22)where the reduced Maxwellian is M vib [ F, G ] = ρ √ πRT exp (cid:18) − | v − u | RT (cid:19) RT exp (cid:16) − εRT (cid:17) , ρ = h F i v,ε , ρu = h vF i v,ε , ρe = (cid:28) ( 12 ( v − u ) + ε ) F (cid:29) v,ε + h G i v,ε , (23)and T is still defined by (7) which implies that T depends both on F and G : to avoid the heavynotation 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 δ ( T ) , like in the case of vibrational energy. However, they are not based on an underlying discretevibrational energy partition, and the authors are not able to prove any H-theorem. Only a localentropy 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, andhence a H-theorem, as it is shown below. System (21–22) naturally satisfies local conservation laws of mass, momentum, and energy. More-over, the H-theorem holds with the reduced entropy H ( F, G ) as defined in (9). Indeed, we havethe Proposition 3.1.
The reduced BGK system (21–22) satisfies the H-theorem ∂ t H ( F, G ) + ∇ x · h vH ( F, G ) i v,ε ≤ , where H ( F, G ) is the reduced entropy defined in (9) .Proof. By differentiation we get ∂ t H ( F, G ) + ∇ x · h vH ( F, G ) i v,ε = h D H ( F, G )( ∂ t F + v ∇ x F ) + D H ( F, G )( ∂ t G + v ∇ x G ) i v,ε = 1 τ (cid:28) D H ( F, G )( M vib [ F, G ] − F ) + D H ( F, G )( δ ( T )2 RT M vib [ F, G ] − G ) (cid:29) v,ε ≤ where we have used (21–22) to replace the transport terms by relaxation ones, and point 5 ofproposition 2.4 to obtain the inequality. It is difficult to derive a Fokker-Planck model for the distribution function f with discrete energylevels. We find it easier to directly derive a reduced model, by analogy with the reduced BGKmodel (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 theBoltzmann collision operator under the assumption of small velocity changes through collisions andadditional equilibrium assumptions (see [8]). In practice, the agreement of this model with theBoltzmann equation is observed even when the gas is far from equilibrium (see [9], for instance).8 .1 A reduced Fokker-Planck model with vibrations First, we remind the Fokker-Planck model for a diatomic gas (without vibrations) obtained in [15]: ∂ t f + v · ∇ x f = D ( f ) , (24)where f = f ( t, x, v, ε ) and the collision operator is D ( f ) = 1 τ (cid:0) ∇ v · (cid:0) ( v − u ) f + RT ∇ v f (cid:1) + 2 ∂ ε ( εf + RT ε∂ ε f ) (cid:1) , where the macroscopic values are ρ = h f i v,ε , ρu = h f v i v,ε , ρe = (cid:28) f (cid:18)
12 ( v − u ) + ε (cid:19)(cid:29) v,ε = 52 ρRT. The internal energy ε can be eliminated by the reduction technique (integration w.r.t dε and εdε )to get ∂ t F + v · ∇ x F = D ( F , G ) ,∂ t G + v · ∇ x G = D ( F , G ) , with the collision operators D F ( F , G ) = 1 τ ∇ v · (cid:0) ( v − u ) F + RT ∇ v F (cid:1) ,D G ( F , G ) = 1 τ ∇ v · (cid:0) ( v − u ) G + RT ∇ v G (cid:1) + 2 τ ( RT F − G ) . Note that the two velocity drift-diffusion terms in the two previous equations have exactly the samestructure as the one in the non-reduced model (24). However, it is interesting to note that theenergy drift-diffusion term of (24) gives, after reduction, a relaxation operator in the G equation.Moreover by reducing the model we lose some moments of initial distribution functions (highermoments in internal energy notably) but we are still able to capture energies and fluxes which aregenerally the main quantities of interest.By analogy, now we propose the following reduced Fokker-Planck model for a diatomic gas withvibrations. Note that now, the model is still with variables x , v , and ε : only the discrete energylevels i are eliminated. This model is ∂ t F + v · ∇ x F = D F ( F, G ) , (25) ∂ t G + v · ∇ x G = D G ( F, G ) , (26)with D F ( F, G ) = 1 τ (cid:0) ∇ v · (cid:0) ( v − u ) F + RT ∇ v F (cid:1) + 2 ∂ ε ( εF + RT ε∂ ε F ) (cid:1) ,D G ( F, G ) = 1 τ (cid:0) ∇ v · (cid:0) ( v − u ) G + RT ∇ v G (cid:1) + 2 ∂ ε ( εG + RT ε∂ ε G ) (cid:1) + 2 τ ( e vib ( T ) F − G ) , (27)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 discretevibrational energy as for the BGK model, since we are not able so far to define a discrete diffusionoperator. As mentioned above, this model is obtained by analogy with the Fokker-Plank modelproposed for polyatomic gases. Its derivation from reduction of a discrete in energy Fokker-Plankmodel will be studied in a future work. 9 .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: h (1 , v ) D F ( F, G ) i v,ε = 0 and (cid:28) ( 12 | v | + ε ) D F ( F, G ) + D G ( F, G ) (cid:29) v,ε = 0 , the reduced entropy H ( F, G ) satisfies the H-theorem: ∂ t H ( F, G ) + ∇ x · h vH ( F, G ) i v,ε = D ( F, G ) ≤ , and we have the equilibrium property ( D F ( F, G ) = 0 and D G ( F, G ) = 0) ⇔ ( F = M vib [ F, G ] and G = e vib ( T ) M vib [ F, G ]) . Proof.
The conservation property is the consequence of direct integration of (27). The equilibriumproperty can be proved as follows. First, note that the Maxwellian M vib [ F, G ] can be written as M vib [ F, G ] = ρ (2 π ) / ( RT ) / exp − (cid:18) v − u ε (cid:19) T Ω − (cid:18) v − u ε (cid:19)! , with Ω = (cid:18) RT
00 2 εRT (cid:19) . 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 bewritten in the compact form D F ( F, G ) = 1 τ ∇ v,ε · (cid:18) Ω M vib ∇ v,ε FM vib (cid:19) ,D G ( F, G ) = 1 τ ∇ v,ε · (cid:18) Ω M vib ∇ v,ε GM vib (cid:19) + 2 τ ( e vib F − G ) . Then an integration by part gives the following identity for D F ( F, G ) : (cid:28) D F ( F, G ) FM vib (cid:29) v,ε = − τ *(cid:18) ∇ v,ε FM vib (cid:19) T Ω M vib ∇ v,ε FM vib + v,ε . Consequently, if D F ( F, G ) = 0 , since the integrand in the previous relation is a definite positiveform, the gradient is necessarily zero, and hence F = M vib . For the equilibrium property of G , theproof is a bit more complicated. First, we have (cid:28) D G ( F, G ) Ge vib M vib (cid:29) v,ε = − τ e vib *(cid:18) ∇ v,ε GM vib (cid:19) T Ω M vib ∇ v,ε GM vib + v,ε + (cid:28) τ ( e vib F − G ) Ge vib M vib (cid:29) v,ε . Consequently, if D G ( F, G ) = 0 , and since F = M vib , we have e vib *(cid:18) ∇ v,ε GM vib (cid:19) T Ω M vib ∇ v,ε GM vib + v,ε = 2 τ (cid:28) ( e vib M vib − G ) Ge vib M vib (cid:29) v,ε = − τ (cid:28) ( e vib M vib − G ) e vib M vib (cid:29) v,ε + 2 τ h e vib M vib − G i v,ε ≤ τ h e vib M vib − G i v,ε = 2 τ ( ρe vib − h G i v,ε ) = 0 , F = M vib . Therefore, we obtain e vib *(cid:18) ∇ v,ε GM vib (cid:19) T Ω M vib ∇ v,ε GM vib + v,ε ≤ , 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 D ( F, G ) = ∂ t H ( F, G ) + ∇ x · h vH ( F, G ) i v,ε satisfies: D ( F, G ) = h D H ( F, G )( ∂ t F + v · ∇ x F ) + D H ( F, G )( ∂ t G + v · ∇ x G ) i v,ε = h D H ( F, G ) D F ( F, G ) + D H ( F, G ) D G ( F, G ) i v,ε , (28)from (21–22). Then the proof is based on the convexity of H ( F, G ) : while for the BGK model weonly 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 D ( F, G ) and multiply by τ so that: τ D ( F, G ) = − X i =1 h ∂ v i ( F ) D H ( F, G ) ( F ( v i − u i ) + RT ∂ v i F ) i v,ε − h ∂ ε ( F ) D H ( F, G ) (
F ε + RT ε∂ ε F ) i v,ε − X i =1 h ∂ v i ( G ) D H ( F, G ) ( F ( v i − u i ) + RT ∂ v i F ) i v,ε − h ∂ ε ( G ) D H ( F, G ) (
F ε + RT ε∂ ε F ) i v,ε − X i =1 h ∂ v i ( F ) D H ( F, G ) ( G ( v i − u i ) + RT ∂ v i G ) i v,ε − h ∂ ε ( F ) D H ( F, G ) ( Gε + RT ε∂ ε G ) i v,ε − X i =1 h ∂ v i ( G ) D H ( F, G ) ( G ( v i − u i ) + RT ∂ v i G ) i v,ε − h ∂ ε ( G ) D H ( F, G ) ( Gε + RT ε∂ ε G ) i v,ε +2 (cid:28) ( e vib ( T ) F − G ) 1 RT log (cid:18) GRT F + G (cid:19)(cid:29) v,ε To use the positive definiteness of the Hessian matrix H of H , we introduce the following vectors: V i = ( F ( v i − u i ) + RT ∂ v i F, G ( v i − u i ) + RT ∂ v i G ) E = ( F ε + RT ε∂ ε F, Gε + RT ε∂ ε G ) , and we decompose the partial derivatives of F and G in factor of D F , D F , D F as follows: ( ∂ v i ( F ) , ∂ v i ( G )) = 1 RT V i − ( F v i − u i RT , G v i − u i RT )( ∂ ε ( F ) , ∂ ε ( G )) = 1 ε E − ( F RT , G RT ) . τ D ( F, G ) = X i =1 (cid:28)(cid:18) F v i − u i RT (cid:19) D H ( F, G ) ( F ( v i − u i ) + RT ∂ v i F ) (cid:29) v,ε +2 (cid:28)(cid:18) F RT (cid:19) D H ( F, G ) (
F ε + RT ε∂ ε F ) (cid:29) v,ε + X i =1 (cid:28)(cid:18) G v i − u i RT (cid:19) D H ( F, G ) ( F ( v i − u i ) + RT ∂ v i F ) (cid:29) v,ε +2 (cid:28)(cid:18) G RT (cid:19) D H ( F, G ) (
F ε + RT ε∂ ε F ) (cid:29) v,ε + X i =1 (cid:28)(cid:18) F v i − u i RT (cid:19) D H ( F, G ) ( G ( v i − u i ) + RT ∂ v i G ) (cid:29) v,ε +2 (cid:28)(cid:18) f RT (cid:19) D H ( F, G ) ( gε + RT ε∂ ε G ) (cid:29) v,ε + X i =1 (cid:28)(cid:18) G v i − u i RT (cid:19) D H ( F, G ) ( G ( v i − u i ) + RT ∂ v i G ) (cid:29) v,ε +2 (cid:28)(cid:18) G RT (cid:19) D H ( F, G ) ( Gε + RT ε∂ ε G ) (cid:29) v,ε − X i =1 (cid:10) V Ti H V i (cid:11) v,ε − (cid:10) E T H E (cid:11) v,ε +2 (cid:28) ( e vib ( T ) F − G ) 1 RT log (cid:18) GRT F + G (cid:19)(cid:29) v,ε Now this expression can be considerably simplified by using property (13), and we get τ D ( F, G ) = X i =1 (cid:28)(cid:18) v i − u i RT (cid:19) ( F ( v i − u i ) + RT ∂ v i F ) (cid:29) v,ε +2 (cid:28) RT ( F ε + RT ε∂ ε F ) (cid:29) v,ε − X i =1 V ti H V i − E t H E − (cid:28) ( e vib ( T ) F − G ) 1 RT log (cid:18) GRT F + G (cid:19)(cid:29) v,ε . Then the first two terms are simplified by using an integration by parts and relations (8) and (7)12o get τ D ( F, G ) = 2 RT ( ρe vib ( T ) − h G i v,ε ) − X i =1 V ti H V i − E t H E +2 (cid:28) ( e vib ( T ) F − G ) 1 RT log (cid:18) GRT F + G (cid:19)(cid:29) v,ε . The terms with the Hessian are clearly negative, since H is positive definite. Then we have τ D ( F, G ) ≤ RT ( ρe vib ( T ) − h G i v,ε )+2 (cid:28) ( e vib ( T ) F − G ) 1 RT log (cid:18) GRT F + G (cid:19)(cid:29) v,ε . Note that from (8) the first term can be written as RT ( ρe vib ( T ) − h G i v,ε ) = 2 RT h e vib ( T ) F − G i v,ε , and can be factorized with the second term to find τ D ( F, G ) ≤ (cid:28) ( e vib ( T ) F − G ) (cid:18) RT log (cid:18) GRT F + G (cid:19) + 1 RT (cid:19)(cid:29) v,ε . We can now prove that the integrand of the right-hand side is non-positive. Indeed, assume forinstance that the first factor is non-positive, that is to say e vib ( T ) F − G ≤ . By using e vib ( T ) = RT e T /T − (see definition (5)), it is now very easy to prove the following relations e vib ( T ) F − G ≤ ⇔ T log (cid:18) GRT F + G (cid:19) ≥ − RT that is to say the second factor of the integrand is non-negative.Consequently, we have proved τ D ( F, G ) ≤ , which concludes the proof. With a convenient nondimensionalization, the relaxation time τ of the reduced BGK model (21)–(22) and the Fokker-Planck model (25)-(26) is replaced by Kn τ , where Kn is the Knudsen number,which can be defined as a ratio between the mean free path and a macroscopic length scale. It isthen 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 modelis ∂ t F + v · ∇ x F = 1Kn τ ( M vib [ F, G ] − F ) , (29) ∂ t G + v · ∇ x G = 1Kn τ ( δ ( T )2 T M vib [ F, G ] − G ) , (30)13here M vib [ F, G ] can be defined by (14) with R = 1 . Similarly, the relations (3)–(7) between thetranslational, internal, and total energies and the temperature, have to be read with R = 1 innon-dimensional variables. We remind that T is still a function of F and G . The Fokker-Planckmodel is ∂ t F + v · ∇ x F = D F ( F, G ) , (31) ∂ t G + v · ∇ x G = D G ( F, G ) , (32)with D F ( F, G ) = 1Kn τ (cid:0) ∇ v · (cid:0) ( v − u ) F + T ∇ v F (cid:1) + 2 ∂ ε ( εF + T ε∂ ε F ) (cid:1) ,D G ( F, G ) = 1Kn τ (cid:0) ∇ v · (cid:0) ( v − u ) G + T ∇ v G (cid:1) + 2 ∂ ε ( εG + T ε∂ ε G ) (cid:1) + 2Kn τ ( e vib ( T ) F − G ) . (33) In this section, we prove the following proposition:
Proposition 5.1.
The mass, momentum, and energy densities ( ρ, ρu, E = ρu + ρe ) of the solu-tions of the reduced BGK and the Fokker-Planck models satisfy the equations ∂ t ρ + ∇ x · ρu = 0 ,∂ t ρu + ∇ x · ( ρu ⊗ u ) + ∇ p = O (Kn ) ,∂ t E + ∇ x · ( E + p ) u = O (Kn ) , (34) which are the Euler equations, up to O (Kn ) . The non-conservative form of these equations is ∂ t ρ + ∇ x · ρu = 0 ,ρ ( ∂ t u + ( u · ∇ x ) u ) + ∇ p = O (Kn ) ,∂ t T + u · ∇ x T + Tc v ( T ) ∇ x · u = O (Kn ) , (35) where c v ( T ) = ddT e ( T ) is the specific heat capacity at constant volume.Proof. The reduced BGK model (21)–(22) is multiplied by , v , and | v | + ε and integrated withrespect to v and ε , which gives the following conservation laws: ∂ t ρ + ∇ x · ρu = 0 ,∂ t ρu + ∇ x · ( ρu ⊗ u ) + ∇ x σ ( F ) = 0 ,∂ t E + ∇ x · Eu + ∇ x · σ ( F ) u + ∇ x · q ( F, G ) = 0 , (36)where σ ( F ) = h ( v − u ) ⊗ ( v − u ) F i v,ε is the stress tensor, and q ( F, G ) = (cid:10) ( v − u )( | v − u | + ε ) F (cid:11) v,ε + h ( v − u ) G i v,ε is the heat flux.When Kn is very small, if all the time and space derivatives of F and G are O (1) with respectto Kn , then (29)–(30) imply F = M vib [ F, G ] + O (Kn ) and G = e vib ( T ) M vib [ F, G ] + O (Kn ) . Thenit is easy to find that σ ( F ) = σ ( M vib [ F, G ]) + O (Kn ) = pI + O (Kn ) , where I is the unit tensor,and q ( F, G ) = q ( M vib [ F, G ] , e vib ( T ) M vib [ F, G ]) + O (Kn ) = O (Kn ) , since the heat flux is zero at14quilibrium, which gives the Euler equations (35). The same analysis can be applied for the reducedFokker-Planck model (31)–(33).Finally, the non conservative form is readily obtained from the conservative form. Note anotherformulation of the energy equation that will be useful below: ∂ t e vib ( T ) + u · ∇ x e vib ( T ) + T e ′ vib ( T ) c v ( T ) ∇ x · u = O (Kn ) , (37)where e ′ vib ( T ) = ddT e vib ( T ) . 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 (Kn ) , the Navier-Stokes equations ∂ t ρ + ∇ · ρu = 0 ,∂ t ρu + ∇ · ( ρu ⊗ u ) + ∇ p = −∇ · σ,∂ t E + ∇ · ( E + p ) u = −∇ · q − ∇ · ( σu ) , (38) where the shear stress tensor and the heat flux are given by σ = − µ (cid:0) ∇ u + ( ∇ u ) T − α ∇ · u (cid:1) , and q = − κ ∇ · T, (39) and where the following values of the viscosity and heat transfer coefficients (in dimensional vari-ables) are µ = τ p, and κ = µc p ( T ) for BGK ,µ = 12 τ p, and κ = 23 µc p ( T ) for Fokker-Planck , (40) while the volumic viscosity coefficient is α = c p ( T ) c v ( T ) − for both models, and c p ( T ) = ddT ( e ( T )+ p/ρ ) = c v ( T ) + R is the specific heat capacity at constant pressure. Moreover, the corresponding Prandtlnumber is Pr = µc p ( T ) κ = 1 for BGK , and for Fokker-Planck . (41)Note that both models do not provide a correct Prandtl number, which can lead to errors forthe computation of fluxes in numerical simulations. This is a usual problem with single parametermodels like BGK or Fokker-Planck: this can be corrected by a modification of the models likewith the ES-BGK or ES-FP approaches, as it has been done for polyatomic gases (see [11, 15] forinstance). The usual Chapman-Enskog method is applied as follows. We decompose F and G as F = M vib [ F, G ] + Kn F and G = e vib ( T ) M vib [ F, G ] + Kn G , which gives σ ( F ) = pI + Kn σ ( F ) , and q ( F, G ) = Kn q ( F , G ) . σ ( F ) and q ( F , G ) up to O (Kn ) . This is done by using the previousexpansions and (21) and (22) to get F = − τ ( ∂ t M vib [ F, G ] + v · ∇ x M vib [ F, G ]) + O (Kn ) ,G = − τ ( ∂ t e vib ( T ) M vib [ F, G ] + v · ∇ x e vib ( T ) M vib [ F, G ]) + O (Kn ) . This gives the following approximations σ ( F ) = − τ h ( v − u ) ⊗ ( v − u )( ∂ t M vib [ F, G ] + v · ∇ x M vib [ F, G ]) i v,ε + O (Kn ) , (42)and q ( F , G ) = − τ (cid:28) ( v − u )( 12 | v − u | + ε )( ∂ t M vib [ F, G ] + v · ∇ x M vib [ F, G ]) (cid:29) v,ε − τ h ( v − u )( ∂ t e vib ( T ) M vib [ F, G ] + v · ∇ x e vib ( T ) M vib [ F, G ]) i v,ε + O (Kn ) . (43)Now it is standard to write ∂ t M vib [ F, G ] and ∇ x M vib [ F, G ] as functions of derivatives of ρ , u ,and T , and then to use Euler equations (34) to write time derivatives as functions of the spacederivatives only. After some algebra, we get ∂ t ( M vib ( F, G )) + v · ∇ x ( M vib ( F, G )) = ρT M ( V ) e − J (cid:18) A · ∇ T √ T + B : ∇ u (cid:19) + O (Kn ) , (44)where V = v − u √ T , J = εT , M ( V ) = 1(2 π ) exp( − | V | A = (cid:18) | V | J − (cid:19) V, B = V ⊗ V − (cid:18) c v (cid:18) | V | + J (cid:19) + e ′ vib ( T ) c v ( T ) (cid:19) I. Then we introduce (44) into (42) to get σ ij ( F ) = − τ ρT (cid:10) V i V j B kl M e − J (cid:11) V,J ∂ x l u k + O (Kn ) , where we have used the change of variables ( v, ε ) ( V, J ) in the integral (the term with A vanishesdue to the parity of M ). Then standard Gaussian integrals (see appendix A) give σ ( F ) = − µ (cid:0) ∇ u + ( ∇ u ) T − α ∇ · u I (cid:1) + O (Kn ) , with µ = τ ρT and α = c p c v − , 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 ∂ t ( e vib M vib ( F, G )) + v · ∇ x ( e vib M vib ( F, G )) = ρT M ( V ) (cid:18) ˜ A · ∇ T √ T + ˜ B : ∇ u (cid:19) + O (Kn ) , (45)where ˜ A = (cid:18) | V | J −
72 +
T e ′ vib ( T ) e vib (cid:19) V, ˜ B = V ⊗ V − (cid:18) c v (cid:18) | V | + J (cid:19) + e ′ vib ( T ) c v ( T ) + T e ′ vib ( T ) c v ( T ) e vib (cid:19) I. q ( F , G ) as given in (43) can be reduced to q i ( F , G ) = − τ ρT (cid:28) | V | V i A j M e − J (cid:29) V,J + (cid:10) V i J A j M e − J (cid:11) V,J ! ∂ x j T − τ ρ D V i ˜ A j M e − J E V,J ∂ x j T. Using again Gaussian integrals , we get q ( F , G ) = − κ ∇ x T, where κ = µc p ( T ) with c p ( T ) = ddT ( e ( T ) + pρ ) = + e ′ vib ( T ) = 1 + c v ( T ) in a non-dimensional form. Here, we rather use the decomposition F = M vib (1 + Kn F ) and G = e vib M vib (1 + Kn G ) , whichgives σ ( F ) = pI + Kn σ ( M vib F ) and q ( F, G ) = Kn q ( M vib F , e vib M vib G ) , in which, for clarity, the dependence of M vib on F and G has been omitted, and the dependenceof e vib on T as well. Finding F and G is less simple than for the BGK model: however, thecomputations 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 D F ( F, G ) = 1 τ M vib L F ( F ) + O (Kn ) ,D G ( F, G ) = 1 τ e vib M vib L G ( F , G ) + O (Kn ) , where L F and L G are linear operators defined by L F ( F ) = 1 M vib (cid:16) ∇ v · ( T M vib ∇ v F ) + ∂ ε (2 T εM vib ∂ ε F ) (cid:17) ,L G ( F , G ) = 1 e vib M vib (cid:16) ∇ v · ( T e vib M vib ∇ v G ) + 2 ∂ ε ( T εe vib M vib ∂ ε G ) + 2( F − G ) (cid:17) . (46)Then the Fokker-Planck equations (31)-(32) suggest to look for an approximation of F and G up to O (Kn ) as solutions of ∂ t M vib + v · ∇ x M vib = 1 τ M vib ( F, G ) L F ( F ) ∂ t e vib M vib + v · ∇ x e vib M vib = 1 τ e vib M vib ( F, G ) L G ( F , G ) . By using (44)-(45), these relations are equivalent, up to another O (Kn ) approximation, to L F ( F ) = τ (cid:18) A · ∇ T √ T + B : ∇ u (cid:19) , and L G ( F , G ) = τ (cid:18) ˜ A · ∇ T √ T + ˜ B : ∇ u (cid:19) , (47)where A , B , ˜ A , and ˜ B are the same as for the BGK equation in the previous section.17ow, we rewrite L F ( F ) and L G ( F , G ) , defined in (46), by using the change of variables V = v − u √ T and G = εT to get L F ( F ) = − V · ∇ V F + ∇ V · ( ∇ V F ) + 2 ((1 − J ) ∂ J F + J ∂ JJ F ) ,L G ( F , G ) = L F ( G ) + 2( F − G ) . Then simple calculation of derivatives show that A , B , ˜ A , and ˜ B satisfy the following properties L F ( A ) = − A, L F ( B ) = − B,L G ( A, ˜ A ) = − A, L G ( B, ˜ B ) = − B. Therefore, we look for F and G as solution of (47) under the following form F = aA · ∇ T √ T + bB : ∇ u and G = ˜ a ˜ A · ∇ T √ T + ˜ b ˜ B : ∇ u, and we find ˜ a = a = − / and ˜ b = b = 1 / .Finally, using these relations into σ and q and using some Gaussian integrals (see appendix A)give σ ( M vib F ) = − µ (cid:0) ∇ u + ( ∇ u ) T − α ∇ · u I (cid:1) and q ( M vib F , e vib M vib G ) = − κ ∇ x T, where α = c p c v − , µ = τ ρT , and κ = µc p ( T ) , which is the announced result, in a non-dimensionalform. In this paper, we have proposed to different models (BGK and Fokker-Planck) of the Boltzmannequation that allow for vibrational energy discrete modes. These models satisfy the conservationand entropy property (H-theorem), and the vibration energy variable can be eliminated by theusual reduced distribution function. The low complexity of the reduced BGK model can make itattractive to be implemented in a deterministic code, while the Fokker-Planck model can be easilysimulated with a stochastic method.Of course, since these models are based on a single time relaxation, they cannot allow for multiplerelaxation times scales. This is not physically correct, since it is known that the relaxation times fortranslational, rotational, and vibrational energies are very different. However, standard procedurescan be used to extend our model, like the ellipsoidal-statistical approach, already used to correctthe Prandtl number of the BGK model [11] and Fokker-Plank models [14, 15].
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 ( V ) = π ) exp( − | V | ) . Wedenote by h φ i = R R φ ( V ) dV for any function φ . It is standard to derive the following integral18elations (see [24], for instance), written with the Einstein notation: h M i V = 1 , h V i V j M i V = δ ij , h V i M i V = 1 , h| V | M i V = 3 , h V i V j V k V l M i V = δ ij δ kl + δ ik δ jl + δ il δ jk , h V i V j M i V = 1 + 2 δ ij h V i V j | V | M i V = 5 δ ij , h| V | M i V = 15 , h V i V j | V | M i V = 35 δ ij , h| V | M i = 105 , while all the integrals of odd power of V are zero. Note that the first relation of each line impliesthe other relations of the same line: these relations are given here to improve the readability of thepaper. From the previous Gaussian integrals, it can be shown that for any × matrix C , we have h V i V j C kl V k V l M i V = C ij + C ji + C ii δ ij . Finally, we have also used the following relations: Z + ∞ J e − J dJ = Z + ∞ e − J dJ = 1 , and also + ∞ X i =0 e − iT /T = 11 − e − T /T and + ∞ X i =0 ie − iT /T = e − T /T (1 − e − T /T ) . References [1] G. A. Bird.
Molecular Gas Dynamics and the Direct Simulation of Gas Flows . Oxford Engi-neering Science Series, 2003.[2] Thomas E. Schwartzentruber Iain D. Boyd.
Nonequilibrium Gas Dynamics and MolecularSimulation . Cambridge Aerospace Series. Cambridge University Press, 2017.[3] G. Dimarco and L. Pareschi. Numerical methods for kinetic equations.
Acta Numerica ,23:369–520, 2014.[4] Luc Mieussens. A survey of deterministic solvers for rarefied flows (invited).
AIP ConferenceProceedings , 1628(1):943–951, 2014.[5] E.P. Gross, P.L. Bhatnagar, and M. Krook. A model for collision processes in gases.
Physicalreview , 94(3):511–525, 1954.[6] C. K. Chu. Kinetic-theoretic description of the formation of a shock wave.
Phys. Fluids , 8(1):12,1965.[7] H. Struchtrup.
Macroscopic Transport Equations for Rarefied Gas Flows Approximation Meth-ods in Kinetic Theory . Interaction of Mechanics and Mathematics. Springer, 2005.[8] C. Cercignani.
The Boltzmann Equation and Its Applications , volume 68. Springer-Verlag,Lectures Series in Mathematics, 1988. 199] M.H. Gorji, M. Torrilhon, and Patrick Jenny. Fokker-Planck model for computational studiesof monatomic rarefied gas flows.
Journal of fluid mechanics , 680:574–601, August 2011.[10] Jr. Lowell H. Holway. New statistical models for kinetic theory: Methods of construction.
Physics of Fluids , 9(9):1658–1673, 1966.[11] P. Andries, P. Le Tallec, J.-P. Perlat, and B. Perthame. The Gaussian-BGK model of boltzmannequation with small prandtl number.
Eur. J. Mech. B-Fluids , pages 813–830, 2000.[12] E. M. Shakhov. Generalization of the Krook relaxation kinetic equation.
Izv. Akad. NaukSSSR. Mekh. Zhidk. Gaza , 1(5):142–145, 1968.[13] M. Hossein Gorji and Patrick Jenny. A Fokker-Planck based kinetic model for diatomic rarefiedgas flows.
Physics of fluids , 25(6):062002–, June 2013.[14] J. Mathiaud and L. Mieussens. A Fokker–Planck model of the Boltzmann equation with correctPrandtl number.
Journal of Statistical Physics , 162(2):397–414, Jan 2016.[15] J. Mathiaud and L. Mieussens. A Fokker–Planck model of the Boltzmann equation withcorrect Prandtl number for polyatomic gases.
Journal of Statistical Physics , 168(5):1031–1055,Sep 2017.[16] Behnam Rahimi and Henning Struchtrup. Capturing non-equilibrium phenomena in rarefiedpolyatomic gases: A high-order macroscopic model.
Physics of Fluids , 26(5):052001, 2014.[17] Zhao Wang, Hong Yan, Qibing Li, and Kun Xu. Unified gas-kinetic scheme for diatomicmolecular flow with translational, rotational, and vibrational modes.
Journal of ComputationalPhysics , 350:237 – 259, 2017.[18] Takashi Arima, Tommaso Ruggeri, and Masaru Sugiyama. Rational extended thermodynamicsof a rarefied polyatomic gas with molecular relaxation processes.
Phys. Rev. E , 96:042143, Oct2017.[19] S. Kosuge, Hung-Wen Kuo, and Kazuo Aoki. A kinetic model for a polyatomic gas withtemperature-dependent specific heats and its application to shock-wave structure. submitted,2019.[20] J. D. Anderson.
Hypersonic and high-temperature gas dynamics second edition . AmericanInstitute of Aeronautics and Astronautics, 2006.[21] T. F. Morse. Kinetic model for gases with internal degrees of freedom.
Phys. Fluids , 7(159),1964.[22] A. B. Huang and D. L. Hartley. Nonlinear rarefied couette flow with heat transfer.
Phys.Fluids , 11(6):1321, 1968.[23] C. Baranger, G. Marois, J. Mathé, J. Mathiaud, and L. Mieussens. A BGK model for hightemperature rarefied gas flows.
Work in progress , 2018.[24] S. Chapman and T.G. Cowling.