A decoupled scheme based on the Hermite expansion to construct lattice Boltzmann models for the compressibleNavier-Stokes equations with arbitrary specific heat ratio

# A decoupled scheme based on the Hermite expansion to construct lattice Boltzmann models for the compressible Navier-Stokes equations with arbitrary specific heat ratio

Kainan Hu [    Hongwu Zhang Industrial Gas Turbine Laboratory, Institute of Engineering Thermophysics,
Shaojuan Geng Industrial Gas Turbine Laboratory, Institute of Engineering Thermophysics,
Chinese Academy of Sciences, Beijing, China
July 14, 2019
###### Abstract

A decoupled scheme based on the Hermite expansion to construct lattice Boltzmann models for the compressible Navier-Stokes equations with arbitrary specific heat ratio is proposed. The local equilibrium distribution function including the rotational velocity of particle is decoupled into two parts, i.e. the local equilibrium distribution function of the translational velocity of particle and that of the rotational velocity of particle. From these two local equilibrium functions, two lattice Boltzmann models are derived via the Hermite expansion, namely one is in relation to the translational velocity and the other is connected with the rotational velocity. Accordingly, the distribution function is also decoupled. After this, the evolution equation is decoupled into the evolution equation of the translational velocity and that of the rotational velocity. The two evolution equations evolve separately. The lattice Boltzmann models used in the scheme proposed by this work are constructed via the Hermite expansion, so it is easy to construct new schemes of higher-order accuracy. To validate the proposed scheme, a one dimensional shock tube simulation is performed. The numerical results agree with the analytical solutions very well.

PACS numbers

47.11.-j, 47.10.-g, 47.40.-x

###### pacs:
47.11.-j, 47.10.-g, 47.40.-x
preprint: APS/123-QED

Also at ]University of Chinese Academy Sciences,Beijing ,China

## I Introduction

The lattice Boltzmann method(LBM) has been successfully applied to isothermal fluidsChen et al. (1992); Chen and Doolen (1998). However, when it is applied to thermal fluids, the LBM encounters some difficulties. One of them is that the specific heat ratio in the macroscopic equations derived from the Bhatnager-Gross-Krook(BGK) equation via the Chapman-Enskog expansion is fixed, in other words, the specific heat ratio is not realistic. Several lattice Boltzmann(LB) schemes with flexible specific heat ratio have been proposedShi et al. (2001); Kataoka and Tsutahara (2004); Watari (2007); Tsutahara et al. (2008). These LB schemes are derived in a similar way. The discrete velocities and the local equilibrium distribution function are determined by a set of constraints which makes sure the macroscopic equations match the thermohydrodynamic equations with certain accuracy. Since 2006, a new way to construct LB models has been developedPhilippi et al. (2006, 2015); Shan et al. (2006); Mattila et al. (2014); Shim (2013a, b); Ansumali et al. (2003); Chikatamarla and Karlin (2006, 2009); Yudistiawan et al. (2010). Contrary to the previous way, the new way derives the discrete velocities and the equilibrium distribution function via the Hermite quadrature and the Hermite expansion. The LB models constructed by the new way are more stable than the LB models constructed by the previous way and it is easy to construct LB models of higher order. In this work, we apply the new way to constructing LB schemes for the compressible Navier-Stokes equations with flexible specific heat ratio. The local equilibrium distribution function including the rotational velocity of particle is decoupled into two parts — one is in relation to the translational velocity and the other is connected with the rotational velocity. The distribution function is also decoupled into two parts accordingly. Two LB models are derived via the Hermite expansion. One is for the distribution function of the translational velocity and the other is for that of the rotational velocity. After this, we decouple the evolution equation into the evolution equation of the translational velocity and that of the rotational velocity. The two evolution equations evolve separately. The decoupled scheme given above is validated by a shock tube simulation. The results of simulation agree with the analytical solutions very well.

## Ii Decoupling the local equilibrium distribution function including the rotational velocity of particle

We begin with the local equilibrium distribution function. The origin that the specific heat ratio is fixed is that gases are supposed to be monatomic, so there is only the translational free degree and the rotational free degree is limited. To describe diatomic gases or polyatomic gases, the rotational velocity of particle should be introduced. The local equilibrium distribution function including the rotational velocity of particle isWatari (2007); Xu et al. (2014)

 feq(ξ,η)= ρ1(2πRgT)D21(2nπRgT)12 ×exp[−(ξ−u)22RgT−η22nRgT], (1)

where is the density, is the absolute temperature, is the macroscopic velocity, is the translational velocity of particle, is the rotational velocity of particle, is the free degree of the rotational velocity of particle, is the dimension and is the universal gas constant.

The dimensionless local equilibrium distribution function is

 (2)

where , , , , , , ,and .

Omitting the tildes on , , , , , we simplify Formula (2)

 (3)

The dimensionless local equilibrium distribution function can be decoupled into the local equilibrium distribution function of and that of

 feq(ξ,η)=geq(ξ)heq(η), (4)

where

 geq(ξ)=ρ(2πθ)D2exp(−|ξ−u|22θ),
 heq(η)=1(2πθ)12exp(−η22θ).

Taking the moment integrals of , we obtain

 ρ= ∫geq(ξ)dξ (5a) ρu= ∫geq(ξ)ξdξ (5b) ρ(et+12u2)= ∫geq(ξ)12ξ2dξ (5c)

where is the translational internal energy.

Taking the moment integrals of , we get

 1= ∫heq(η)dη (6a) er= ∫heq(η)n2η2dη (6b)

where is the rotational internal energy.

Taking the moment integrals of the local equilibrium distribution function, we obtain

 ρ= ∫∫feq(ξ,η)dξdη (7a) ρu= ∫∫feq(ξ,η)ξdξdη (7b) ρ(E+12u2)= ∫∫feq(ξ,η)(12ξ2+n2η2)dξdη (7c)

where is the internal energy. It is the sum of the translational energy and the rotational energy .

According to the kinetic theory, we get

 ρ= ∫∫f(ξ,η)dξdη (8a) ρu= ∫∫f(ξ,η)ξdξdη (8b) ρ(E+12u2)= ∫∫f(ξ,η)(12ξ2+n2η2)dξdη (8c)

Now we assume the translational velocity of particle is independent of the rotational velocity of particle, so the distribution function can be decoupled into and

 f(ξ,η)=g(ξ)h(η), (9)

where we define

 f(ξ,η)=ρb(ξ, η)=ρb1(ξ)h(η), (10a) g(ξ)= ∫f(ξ,η)dη, (10b) h(η)= ∫b(ξ,η)dξ. (10c)

Section(III) and Appendix will discuss the reasonableness of this assumption.

It should be noticed that is a joint probability distribution, and are marginal probability distributions, is the product of and a marginal probability distribution.

According to Formula(8) and (10), the moments of the distribution function are

 ∫g(ξ)dξ= ∫∫f(ξ,η)dξdη=ρ, (11a) ∫g(ξ)ξdξ= ∫∫f(ξ,η)ξdξdη=ρu, (11b) ∫g(ξ)12ξ2dξ= ∫∫f(ξ,η)12ξ2dξdη=ρ(et+12u2). (11c)

Similar to Formula(11), the moments of the distribution function are

 ∫h(η)dη= ∫∫b(ξ,η)dξdη=1 (12a) ∫h(η)n2η2dη= ∫∫b(ξ,η)n2η2dξdη=er (12b)

## Iii Decoupling the evolution equation

We have decoupled the equilibrium distribution function and the distribution function in Section II

 feq(ξ,η)= geq(ξ)heq(η), f(ξ,η)= g(ξ)h(η).

After decoupling and , we can decouple the evolution equation of . The evolution equation of can be expressed as

 ∂f(ξ,η)∂t+ξ⋅∇f(ξ,η)=−1τ[f(ξ,η)−feq(ξ,η)]. (13)

where is the relaxation time.

Integrating Formula(13) on , and substituting Formula(10b) into it, we obtain the evolution equation of

 ∂g(ξ)∂t+ξ⋅∇g(ξ)=−1τ[g(ξ)−g(ξ)eq]. (14)

Substituting Formula(4) and Formula(9) into Formula(13), we obtain

 ∂g(ξ)h(η)∂t+ξ⋅∇g(ξ)h(η)=−1τ[g(ξ)h(η)−geq(ξ)heq(η)]. (15)

Expanding Formula(15) and simplifying it, we obtain

 h(η)[∂g(ξ)∂t+ ξ⋅∇g(ξ)]+g(ξ)[∂h(η)∂t+ξ⋅∇h(η)] =−1τ[g(ξ)h(η)−geq(ξ)heq(η)]. (16)

Substituting Formula(14) into (III), integrating on and simplifying it, we obtain the evolution equation of

 ∂h(η)∂t+u⋅∇h(η)=−1τ[h(η)−h(η)eq]. (17)

Discretizing the evolution equation of and in the discrete velocity space, we obtain the discrete evolution equations of and

 ∂gi∂t+ξi⋅∇gi=−1τ(gi−geqi), (18a) ∂hj∂t+u⋅∇hj=−1τ(hj−heqj), (18b)

Where and are the discrete form of and . Formula(18a) is the evolution equation of the discrete translational velocity and formula(18b) is the evolution equation of the discrete rotational velocity . It should be noticed Formula(18a) is independent of and formula(18b) is indirect connected to Formula(18a) via the macroscopic velocity .

From the two evolution equations of and i.e. Formula(18a) and (18b), the Navier-Stokes equations with flexible specific heat ratio via the Chapman-Enskog expansion can be derived,

 ∂∂tρ +∇⋅ρu=0, (19a) ∂∂tρ u+∇⋅(ρuu+Pδ) = ∇⋅μ[(∇u+u∇)−2D+n∇⋅uδ], (19b) ∂∂tρ (E+12u2)+∇⋅ρu(E+12u2+Pρ) = ∇⋅μu(∇u+u∇−2D∇⋅uδ)+∇⋅κ∇E (19c)

where is the pressure, is the dynamic viscosity coefficient, is the heat conductivity, and the specific heat ratio is defined as

 γ=D+n+2D+n. (20)

The derivation shows that it is reasonable to assume the distribution function can be decoupled into and .

The appendix will give the derivation in details.

## Iv LB models for the translational velocity and the rotational velocity

In this section, we derive LB models from and respectively via the Hermite expansion. The process of deriving LB models via the Hermite expansion has been discussed intensively by X.ShanShan et al. (2006, 2010), C.PhilippiPhilippi et al. (2006, 2015); Mattila et al. (2014) and JW.ShimShim (2013a, b). In this work, we only discuss two-dimensional fluids. The case of three dimension is similar. Employing the Hermite expansion, we construct a two-dimensional LB model of fourth-order accuracy, i.e. D2Q37, from . The discrete particle velocities and the weights of D2Q37 are showed in Table(1). The discrete equilibrium distribution function of D2Q37 is

 geqi(ξ)=ωiρ4∑k=01k!a(k)⋅H(k), (21)

where

 a(0)⋅H(0)= 1, a(1)⋅H(1)= ξ⋅u, a(2)⋅H(2)= (ξ⋅u)2+(θ−1)(θ2−D)−u2, a(3)⋅H(3)= (ξ⋅u)[(ξ⋅u)2−3u2 +3(θ−1)(u2−D−2)], a(4)⋅H(4)= (ξ⋅u)4−6(ξ⋅u)2u2+3u4 +6(θ−1)[(ξ⋅u)2(u2−D−4) +(D+2−u2)ξ2] +3(θ−1)2[u4−2(D+2)u2+D(D+2)],

and .

It should be noticed that the LB model given above is different from the LB model given byPhilippi et al. (2006).

In a similar way, a one-dimensional LB model of fourth-order accuarcy can be derived from . Here, we adopt the D1Q7 model proposed by JW.ShimShim (2013b).The discrete velocities and the weights are shown in Table(II).

The discrete equilibrium distribution function is

 heqj(η)=ωj4∑j=01j!a(j)⋅H(j), (22)

where

 a(0)⋅H(0)= 1, a(1)⋅H(1)= 0, a(2)⋅H(2)= (θ−1)(θ2−D), a(3)⋅H(3)= 0, a(4)⋅H(4)= 3(θ−1)2[η4−2(D+2)η2+D(D+2)],

and .

D2Q37 and D1Q7 are both models of fourth-order accuracy, so the scheme given above is of fourth-order accuracy. In this way, the higher-order of accuracy can be achieved easily.

We can also construct or adopt other LB models. But it should be noticed that the scaling factors of LB models derived form and should be equal or else interpolation is necessary.

## V Calculation procedure

In this section, we first discretize the evolution equations of and in time and space, then we give the computational algorithm.

### v.1 Discretized the evolution equation in space and time

Now we discretize the discrete evolution equations of and in time and space. The first order difference is employed for the time discretization and the convection term is performed by the third order upwind scheme. The discretized form of Formula(18a) is

 gi(x,t+Δt)=gi(x,t)−Δtξi⋅∇gi−Δtτ[gi(x,t)−geqi(x,t)], (23)

where is the time increment, and the convection term along the coordinate is

 ξix∂gi∂x=12ξi,x+|ξi,x|6Δx[gi(x−2Δx,y)−6gi(x−Δx,y)+3gi(x,y)+2gi(x+Δx,y)]+12ξi,x−|ξi,x|6Δx[−gi(x+2Δx,y)+6gi(x+Δx,y)−3gi(x,y)−2gi(x−Δx,y)],

and is the space increment. The convection term along the coordinate is similar. In a similar way, the discretized form of the discrete evolution equation of , i.e. Formula(18b), is

 hj(x,t+Δt)=hj(x,t)−Δtu⋅∇hj−Δtτ[hj(x,t)−heqj(x,t)], (24)

where the convection term is similar with that of Formula(18a)

 ux∂hj∂x=12ux+|ux|6Δx[hj(x−2Δx,y)−6hj(x−Δx,y)+3hj(x,y)+2hj(x+Δx,y)]+12ux−|ux|6Δx[−hj(x+2Δx,y)+6hj(x+Δx,y)−3hj(x,y)−2hj(x−Δx,y)].

The convection term along the coordinate is similar.

### v.2 Computational algorithm

The computational algorithm is as fellow :

(1) Update by Formula(23);

(2) Update by Formula(24);

(3) Calculate the density , the macroscopic velocity and the translational internal energy

 ρ= ∑igi (25a) ρu= ∑igiξi (25b) 12ρu2+ρet= ∑igi12ξ2i (25c)

where the translational internal energy is

 et=D2T. (26)

The pressure is defined as ;

(4) Calculate the rotational internal energy

 er=∑jn2hjη2j, (27)

and combine Formula(27) with(25c), then we obtain

 12ρu2+ρE=∑igi12ξ2i+ρ∑jhjn2η2j. (28)

As defined above, the internal energy is the sum of the translational internal energy and the rotational internal energy

 E=et+er=D+nDet. (29)

Substituting Formula(29) into (28) we obtain the internal energy

 E=∑igi12ξ2i+ρ∑jn2hjη2j−ρ12u2ρ. (30)

Substituting Formula(26) and (29) into (30) we obtain the absolute temperature

 T=2D+n∑igi12ξ2i+ρ∑jn2hjη2j−ρ12u2ρ. (31)

(5) Implement the boundary conditions.

## Vi Numerical validation

In this section, we apply the decoupled scheme given above to simulating a shock tube. The grid is . The initial condition of the left tube is and that of the right tube is . The specific heat ratio is , the rotational free degree is and the relaxation time is . All of these macroscopic variables are dimensionless. The periodic boundary condition is employed for the up and down boundaries and the open boundary condition is employed for the left and right boundaries.

Fig(1) gives the results of simulation employing the decoupled scheme given above at , i.e. time

 t=stepX×r=1801000×1.1969797752=0.1504.

The analytical solutions Sod (1978) at the same time are also given. The black lines show the simulation results and the gray lines show the analytical resolutions. It can be seen from Fig(1) that the simulation results agree with the analytical resolutions very well.

Tab(3) gives the relative error of density , pressure , absolute temperature and the velocity . The relative error is defined as

 Error=∑i|xnume,i−xanal,i|∑i|xanal,i|, (32)

where is the numerical solutions, is the analytical solutions and . Tab(3) shows that the maximum of relative error is . This relative error is acceptable.

## Vii Conclusion

This work proposes a way based on the Hermite expansion to construct LB schemes for the compressible Navier-Stokes equations with arbitrary specific heat ratio. The equilibrium distribution function , the distribution function and the evolution function are decoupled into two parts, namely one is in relation to the translational velocity and the other is connected with the rotational velocity . The two evolution equations evolve separately. The translational velocity is discretized in a two- or three- dimensional LB model and the rotational velocity is discretized in another one-dimensional LB model. The Hermite expansion is applied to deriving these two LB models. The correct flexible specific heat ratio is obtained and correct relation between the temperature and the internal energy is derived via the Chapman-Enskog expansion. The decoupled scheme is validated by a shock tube simulation. The simulation results agree with the analytical resolutions very well.

The LB models used in the decoupled scheme is same as the ones used in the schemes with fixed specific heat ratio. It is not necessary for the decoupled scheme to construct new LB models specially. The models with fixed specific heat ratio can applied to the decoupled scheme without any recommendation. This is different from the existing schemes which construct new models in order to adjust the specific heat ratio.

The decoupled scheme proposed by this work can make use of the models constructed via the Hermite expansion, so the process of constructing new schemes is simple and higher-order accuracy can be achieved easily. For the same reason, the decoupled scheme is more stable than the schemes constructed by the early way, which derives LB models via a try-error method.

*

## Appendix A Derivation of the Navier-Stokes equations from the evolution equations of g(ξ) and h(η) via the Chapman-Enskog expansion

In this appendix, we derive the Navier-Stokes equations with flexible specific heat ratio from the evolution equations of and via the Chapman-Enskog expansion.

Expanding the distribution functions and , the derivatives of the time and the space in terms of the Kundsen number we obtain

 ∇= ϵ∇1, (33a) ∂∂t=ϵ∂∂t1 +ϵ2∂∂t2, (33b) gi=g(0)i+ϵ g(1)i+ϵ2g(2)i, (33c) hj=h(0)j+ϵ h(1)j+ϵ2h(2)j. (33d)

Substituting Formula(33c) into the evolution equation of the translational velocity, i.e. Formula(18a) and comparing the order of we obtain

 g(0)i=g(eq)i, (34a) (∂∂t1+ ξi⋅∇1)g(0)i+1τg(1)i=0, (34b) ∂g(0)i∂t2+(∂∂t1 +ξi⋅∇1)g(1)i+1τg(2)i=0. (34c)

Considering the discrete form of Formula(11) in the discrete velocity space

 ∑igi= ∑igeqi=ρ, (35a) ∑igi= ∑igeqi=ρu, (35b) ∑igi12ξ2i= ∑igeqi12ξ2i=ρ(12u2+et), (35c)

we obtain

 ∑ig(n)i=0,∑ig(n)iξi=0,∑ig(n)i12ξ2i=0,n=1,2. (36)

Substituting Formula(33d) into the evolution equation of the rotational velocity,i.e.Formula(18b) and comparing the order of we obtain

 h(0)j=h(eq)j, (37a) (∂∂t1+ u⋅∇1)h(0)j+1τh(1)j=0, (37b) ∂h(0)j∂t2+(∂∂t1 +u⋅∇1)h(1)j+1τh(2)j=0. (37c)

Considering the discrete form of Formula(12) in the discrete velocity space

 ∑jhj= ∑jheqj=1, (38a) ∑jhjn2η2j= ∑jheqjn2η2j=er, (38b)

we obtain

 ∑jh(n)j=0,∑jh(n)jn2η2j=0,n=1,2. (39)

Some velocity moments of and will be used in the derivation of the Navier-Stokes equations and we list them as fellow

 ∑igeqi= ρ, (40a) ∑igeqiξi= ρu, (40b) ∑igeqiξiξi= ρuu+Pδ, (40c) ∑igeqiξiξiξi= ρuuu+Puδ, (40d) ∑igeqi12ξ2i= ρ(12u2+et), (40e) ∑igeqi12ξ2iξi= ρ(12u2+et)u, (40f) ∑igeqi12ξ2iξiξi= P(2Det+12u2+et)δ, +(2P+12ρu2+ρet)uu, (40g)

where . Here, the Grad notes is usedGrad (1949). Two velocity moments of will be used in the following parts

 ∑jheqj =1, (41a) ∑jheqjn2η2 =er. (41b)

### a.1 Derivation of the continuity equation

Taking the zeroth order moment of Formula(34b), we obtain

 ∑i(∂∂t1+ξi⋅∇1)g(0)i+1τ∑ig(1)i=0 (42)

Substituting Formula(36) and (40) into (42), the continuity equation of the first order is obtained

 ∂∂t1ρ+∇1⋅ρu=0. (43)

Taking the zeroth order moment of Formula(34c), we obtain

 ∂∑ig(0)i∂t2+∑i(∂∂t1+ξi⋅∇1)g(1)i+1τ∑ig(2)i=0. (44)

Substituting Formula(36) and (40) into (44), then summing on we obtain the continuity equation of the second order

 ∂ρ∂t2=0. (45)

Making use of Formula(33b) and combining the continuity equation of the first and second order, i.e. Formula(43) and (45), the continuity equation is obtained

 ∂∂tρ+∇⋅ρu=0. (46)

### a.2 Derivation of the momentum conservation equation

Taking the first order moment of Formula(34b)

 ∑i(∂∂t1+ξi⋅∇1)g(0)iξi+1τ∑ig(1)iξi=0, (47)

and inserting Formula(36) and (40), we get the conservation momentum equation of the first order

 ∂∂t1ρu+∇1⋅(ρuu+Pδ)=0. (48)

Taking the first order moment of Formula(34c)

 ∂∑ig(0)iξi∂t2+∑i(∂∂t1+ξi⋅∇1)g(1)iξi+1τ∑ig(2)iξi=0, (49)

substituting Formula(34b) into (49) and simplifying, we obtain

 ∂∑ig(0)iξi∂t2−τ∇1⋅ (∂∂t1∑iξiξig(0)i +∇1⋅∑iξiξiξig(1)i)=0. (50)

Inserting the moments of i.e. Formula(40), into Formula(A.2) we obtain

 ∂ρu∂t2−τ∇1⋅[∂∂t1(ρuu+Pδ)+∇1⋅∑iξiξiξig(1)i]=0, (51)

where is a difficult point to simplify. To simplify , we should obtain the energy conservation equation of the first order firstly.

Multiplying to Formula(34b) and summing on we obtain

 ∑i(∂∂t1+ξi⋅∇1)g(0)i12ξ2i+1τ∑ig(1)i12ξ2i=0. (52)

Substituting the moments of i.e. Formula(40), into Formula(52), we obtain the translational internal energy conversation equation of the first order

 ∂∂t1(12ρu2+ρet)+∇1⋅(12ρu2+ρet+P)u=0. (53)

Multiplying to Formula(37b), summing on and inserting the moments of ,i.e. Formula(41), we obtain the rotational internal energy conversation equation of the first order in nonconservation form

 ∂∂t1er+u⋅∇1er=0. (54)

Multiplying to Formula(54) and multiplying to Formula(43), then adding up we obtain

 ρ(∂∂t1er+u⋅∇1er)+er(∂∂t1