ACTA MECHANICA SINICA (English Series), Vol.13, No.3, August 1997 The Chinese Society of Theoretical and Applied Mechanics Chinese Journal of Mechanics Press, Beijing, China Allerton Press, INC., New York, U.S.A.
ISSN 0567-7718
A LATTICE BOLTZMANN MODEL FOR C O M P R E S S I B L E P E R F E C T GAS*
Hu Shouxin ( $ ] ~ ) Yan Guangwu ( ~ ' ~ ) Shi Weiping ( ~ )
(Mechanics Section, Department of Mathematics, Jilin University, Changchun 130023, China)
ABSTRACT:
A new lattice Boltzmann model for compressible perfect gas is proposed. The numerical example shows that it can be used to simulate shock wave and contact discontinuity. The results are comparable with those obtained by traditional methods. The ratio of specific heats 7 may be chosen according to the requirement of problems.
K E Y W O R D S : lattice Boltzmatm method, compressible, perfect gas, general case of the ratio of specific heats, two-speed three-energy-level model
1 INTRODUCTION Ten years ago, the Lattice Gas Automata (LGA) model was presented by Frisch, Hasslacher and Pomeau [1] to simulate the two-dimensional fluid flows. It has attracted much attention. Not very long since then, its rapid developments were obstructed by two fundamental difficult problems, that is, the density dependence of convection coefcient g(p) and the nonphysical velocity dependence of pressure. In 1992, these difficulties were removed in the lattice Boltzmann Method (LBM) developed by Qian and d'Humi~res et alfl], Chen HD and Chen SY et al.[3], Benzi and Succi et al. [4]. The LBM provides a 7-bit isothermal, incompressible model of fluid problems obeying the Navier-Stokes equations. Although it has showed strong spatial gradient phenomena in some examples of the wave propagation with initial large perturbation [5], generally speaking, it can only be used in incompressible problems with small velocity limit. Recently, many studies of the LBM has been concentrated on the compressible flows. Alexander and Chen et al.[6] formulated an isothermal LBM model which can include shocks. It has a selectable sound speed. This feature allows one to simulate compressible fluid flows with high Mach numbers. In paper [7], Qian and Orszag studied the nonlinear deviation of the LBGK model in compressible regimes, and presented a numerical simulation of a shock profile. Qian and Orszag Is] also simulated a regime of weak compressibility at high Reynolds numbers. The LBM scheme was used to study the Kolmogorov flows. Ancona [9] developed a class of "fully-Lagrangian" methods, and provided new perspective on the relationship Received 27 May 1996, revised 14 April 1997 * The project supported by the National Natural Science Foundationof China
Vo1.13, No.3
Hu Shouxin et al.: A Lattice Boltzmann Model for Compressible Gas
219
between LBM and finite difference techniques. Two simpler areas of application of this scheme were Burger's equations and one dimensional gas dynamics. However, the LBM model for computing strong discontinuity in gas dynamics has not been carried out yet. In Section 2 of this paper, we describe a new LBM model for two-dimensional compressible perfect gas. In the sense of first order of accuracy the Euler equations can be recovered from the model, and the ratio of specific heats "r may be chosen as a parameter. The problem of parameter choice is discussed in Section 3, while in Section 4, a numerical example, the Sod's problem of one dimensional shock tube, is given. It shows the formation of shock wave and contact discontinuity. 2 LATTICE BOLTZMANN MODEL FOR COMPRESSIBLE PERFECT GAS
The model presented in this paper is based on the FHP-III 7-bit model and the available LBM model, but we assume that the particles moving along every link are separated into two kinds, type A (a = 1 , . . - , 6) and type B ( a = 7 , . . . , 12), which are on two energy levels eA and eB(eA > 6B > 0), respectively, and the rest particles (a = 0) are on another energy level 6D(> 0). So it is actually a 13-bit model. The single particle distributions in the "shooting-in" and "shooting out" state at site r and time t are respectively denoted by fa (r, t) and f'a (r, t). We define the mass, momentum and total energy per site as follows
p=_,So
pu, = ~-~f~e~i
Oz
(1)
( i = 1,2)
(2) (3)
1
2
where e~ is the velocity vector of the moving particles along the link in the direction a, leal = c (a = 1,.--, 12), and E is the internal energy per unit mass. The evolution of the system from time t to t + A t is still divided into two steps: collision and streaming, and described by the following BGK-type Lattice Boltzmann Equation (LBE)
l(l~(r,t)-f~q(r,t)) f ~ ( r + e~At, t + At) = f ~ ( r , t ) = f ~ ( r , t ) -- T
I
(a=0,1,...
12)(4)
where r is the single-relaxation time, f~q is the local equilibrium distribution. We assume that f~q has the same expression as that in the available 7-bit LBM [2~4]
feq = Aop + All--levi + A2p?ziuje~iec~ j + A3pu 2 feq
=
(a = 1 , . . . , 6)
] (5)
Bop + Blpuieai + B2puiuieaiea j
+
S3pu 2
(OZ= 7 , . - . , 12) /
f~q = Dop + D3pu 2
We adopt the idea proposed by Yan [1~ that besides the conservation conditions of mass, momentum and energy:
220
ACTA MECHANICA SINICA (English Series)
E f ~ q~-'~-p
E
ct
1997
(6) (7)
2 +pE
feaqeai --" I~i
Z ~q6~
o~
1
(8)
the equilibrium distribution must satisfy the following flux conditions (the conditions of moment of higher order):
y f~ eq e,~ie~j = puiuj + P~ij
ot
(9)
ZJ:q6aeai-~-(lpu2+pE+p)ui
~t
(10)
The pressure p will be obtained from the equations of state and energy of perfect gas: P = (7 - 1)pE (11)
and the temperature T is in proportion to E. Substituting Eq.(5) into Eqs.(6),,~(10),we easily obtain the system of linear algebraic equations which determine the coefficientsAi,Bi,Di:
Ao + ]30 = (1 - Do)/b
Ao + Bo = (7 - 1)DE~ bcz
(12)
(13)
(14) (15) (16) (17) (18) (19)
- -
eAAo + eBBo = (E - eDDo)/b A1 q- B1 = D/bc 2 eAAI + eBB1 = D(u z + 27E)/2bc 2 A2 + B2 = D(D + 2)/2bc4 A3 q- B3 = - D / 2 b e 2
bcu b(A3 + B3) + D3 + ~ - ( A 2 + B2) = 0 beu 1 b(eAA3 + esB3) + "~'(eAA2 "beBB2) = ~
coD3
(20)
Here, D (= 2) is the space dimension, b (= 6 for hexagonal lattice) is the link number per site. We have D a2 Do = 1 - ( 7 - 1)Dc~ = 1 c2 (21) from Eqs.(12),~(13), and Da = - 1 / c 2 from Eqs.(17)~-(19). Here a = ~ / ~ speed. By introducing a complemental equation
eAA2 + eBB2 =
is the sound
D A-~'~
(22)
Vol.13, No.3
Hu Shouxin et al.: A Lattice Boltzmann Model for Compressible Gas
221
all pairs Ai, Bi can be solved easily from a system of two linear equations in two unknowns. Here, A is a chosen parameter, called separating factor. For example, by solving Eqs.(15) and (16) we get
A I = bc2 k 2 q- 7 E -
gB)/(SA ~ - eS)
(23) (241
D u2 BI= -~(---~--TE+eA)/(eA--e.)
It is worth while to point out that the coefficients Ai, Bi and Di may be functions of macroscopic state variables. Equation (4) is a finite difference scheme of f4- Its Taylor expansion (divided by At) is
014 014 Ot F e 4 i ~ +
A t ( i)2 f a . , 0214 0214 ~ 1 9 "-+ 0 ( A t 2) = - - v A t ( f 4 - fa eq) (25) 2 \ at 2 -ffze41"~--~.-'toxio~ e4'e43oziOxj)
.
Multiplying the above equation by 1, e4i, ea separately, summing over all a and noting Eqs.(6),~(8), we have Op Opui
+ ~
+ R~ = o
(26)
Opu~ O~r~j+ , 0--7- + ~-~j P ~ = 0
(27/ (28)
O[1 2 ~ OQi , ~spu + pE) + ~ + R3 = 0
where
~,~ = ~ Io,4,e4~,
o~
Q' = Z i4e~
4
(29)
According to Chapman-Enskog theory, we assume
fa
~-
feq
"Jr-
kfa(1) + k2f(2) +
"'"
(30)
where k is the Knudsen number and is a small parameter with the scaling condition k ,v A t = As/c, with As as the lattice unit. Putting (30) into (29), we have 7rij = ~ f~qe4ie~j + 0(At) = puiu~ + p&~ + 0(At)
Q, = 2
4
+ o(.,) :
+
+ .)., +
Substituting the above equations into (26)~(28), we obtain
g +-~-+
Op
Optq
R1 = 0
(31) (32)
+pu,)
Opu___!+ Opuiuj + P~ij + R2 = 0 Ot Oxj O(lpu2+pE)+ 0 '1 2 0
(33)
222
ACTA MECHANICA SINICA (English Series)
1997
The first two terms of the LHS of Eqs.(31)~(33) are the terms of Euler equations of gas dynamics, and Ri's are the truncation errors, Ri = O(At). In other words, the macroscopic Euler equations axe recovered from the difference scheme (4) in the sense of first-order accuracy. The same is true when using the technique of multiscale expansion.
3 CHOICE OF PARAMETERS
There are six parameters, CA, eB, eD, )~, C and % to be chosen in this model. The ratio of specific heats 7 may be chosen as well, but usually it is given beforehand according to the physical condition. The methods we used for choosing parameters are (1) The considerations based upon physical concepts In equilibrium states, the mass and internal energy on every energy level, as the total mass and internal energy, should be positive at every site, and the momentum on level A, B should not be opposite to the total momentum. We propose a little different but simpler criterion to meet these requirements, that is, the coefficients Ai, Bi, Di should not have opposite sign for all i. From Eqs.(13), (15), (17) and (18) we see, this criterion is equivalent to the conditions that Ai > O, Bi > O, Di > 0 (i = 0,1,2) and A3 < 0, B3 < 0, D3 < 0. The condition
zq= (1-~ leads to
C2 > It 2 A- D---a2 (35)
This implies the condition Do > 0. Actually, 1/c = A t ~ A s is the mesh ratio, and the condition (35) is similar to the CFL condition. The conditions A1 > 0 and B1 > 0 lead to (see Eqs.(23), (24))
u2
eA > - ~ + 7 E > eB
(36)
and from the conditions A0 > 0, B0 > 0, we get E - (1 - Do)eB > Doeo > E - (1 - Do)eA
(37)
(36) and (37) are the conditions for parameters eA, eB, eD. From the conditions A2 > 0, B2 > 0 and Aa < 0, B3 < 0, we get two conditions for parameter ~, (D + 2)eA (D + 2)es
2C2
>
>
2C2
(38)
+
In addition, Eq.(9) leads to
(eD +
12 1
> .X >
+
(39)
1~
feqr
1
2 +p
(40)
Vo1.13, No.3
Hu Shouxin et al.: A Lattice Boltzmann Model for Compressible Gas
223
Subtracting Eq.(40) from Eq.(8), we have
12
1
Let eA > eB > c2/2, ~D > 0, the internal energy pE will be positive if all f~q are positive. Combining with Eq.(36), we suggest an upper bound condition of speed c. u~ + 2 7 E = u 2 + (2) The requirement of numerical stability The detailed analysis on the residual terms P~ of Eqs.(31).~(33) gives a complicated expression of numerical viscosity terms, which verify the reasonableness of the above criterion based on physical concepts (presented in another paper), but the expression is different from the viscosity terms in the equations of the macroscopic fluid mechanics. So the LBM simulation of gas dynamics with physical viscosity is still left to be solved. The common factor k(r - 1/2) in the numerical viscosity terms should be positive to guarantee the positivity of numerical viscosity, namely, we must let r > 1/2, or
0 < 1/~" < 2
2 a2>c2 "y-1
(42)
Using the above conditions to choose proper parameters, we still have to do some numerical experiments. Our experiments show that, the numerical results are not sensitive to parameters eA, eB, ZD and A, they are relatively satisfied even if those conditions involved do not hold; however the CFL-like restriction (35) is important, its invalidation will generate numerial instability. The parameter ~" is an important factor which influence the width of discontinuities, a larger ~- yields wider discontinuities, if T - 1/2 is too small, probably, the oscillation behind discontinuities will occur, so as to give rise to numerical instability.
4 NUMERICAL RESULTS
A test problem, the Sod's problem In] with initial data on the left and right side
(pL,UL,PL) = (1,0,1) -0.5 < x < 0 0 < x < 0.5
(PR,uR,pR) =
(0.125,0,0.1)
is calculated by using this model. The numerical results are compared with the theoretical ones (Fig.l). It shows the formation of shock wave and contact discontinuities of density and temperature. Their widths are approximately those obtained by traditional numerical methods. The position of shock wave coincides with the theoretical predication. The lively process of shock evolution can be seen in a computer screen. The quantitative comparisons between numerical and theoretical results on some "platforms" in Fig.1 are given in Tab.1.
224
P
ACTA MECHANICA SINICA (English Series)
1
0.8 0.6 0.4 0.2 -0.6 -014 -012
1997
~X
P
1
0.8 0.6 0.4 0.2
o'.2 '0.4 - - 0.6
X
-0.6-0.4-0.2 0 0.2 0.4 0.6
(b)
E
X
(a)
1
ll
0.8 0.6 0.4
0.2 j 0 ~_-d
I
F
m
3 2.8 % ~ ~
2.6 2.4; 2.2 2 1.8 -0.6
-0.6 -0,4 -01.2
012
01.40.6
X
-0:4 -012 ;
(d)
012 0:4 0.6
X
(c)
Fig.1 Comparisons between numerical and theoretical results
E x a c t solution (line) and simulation results (circles) of p, p, u, a n d E are shown. Lattice size: 200 x 2. O u t p u t at 120 t i m e steps. P a r a m e t e r s : 7 = 1.4; c = 3; ~ = 1.75; 1 / r = 1.45; ~A = 18; ~V = 5.4; ~D = 0.9.
Table 1 Comparisons between numerical
and theoretical results
on some "plateforms"
P(==o.2) P(==o.x) P(==o.3) u(==o.=) E(==o.1) E(==o.3)
This m o d e l Theoretical results 0.300 1 0.3031 0.424 0 0.426 3 0.265 5 0.265 5 0.928 8 0.927 5 1.787 5 1.777 6 2.853 3 2.853 5
5 CONCLUDING
REMARKS
Some key points in our paper are as follows: (1) It is necessary to define tota]~ energy and internal energy for the recovery of energy equation. In many papers, the defined total energy is the total kinetic energy of particles 12 ET = ~ fac2/2, which corresponds to ~A = gB = ~ c ,gD = 0. This definition coincides with the physical concepts of perfect gas, but in one or two speed model it brings us two difficult problems: (i) it leads to 7 = 2 (so-called Ideal Case [9], see Eq.(41)), (ii) the
energy conservation can be derived from momentum flux conditions (9). To solve these problems, many researchers use multi-speed model (e.g. papers [12~15]) or introduce
Vo1.13, No.3
Hu Shouxin et al.: A Lattice Boltzmann Model for Compressible Gas
225
the concept of energy level (e.g. paper [16]). We take advantage of the merits of both, the present model is not only of multi-speed but also of multi-energy level. As a result, all equations of perfect gas are successfully included in the lattice Boltzmann model, and the ratio of specific heats appears as a chosen parameter (so-called General Case). The problems remained are those of accuracy and numerical stability. The other advantage is that the pressure p (or internal energy E ) in this model is a statistical quantity independent of p and pu~, so that, a wide range of sound speed (c8 = V ~ / P ) is allowable. (2) We adopt an idea that the local equilibrium distribution must satisfy conservation conditions and flux conditions of mass, momentum and energy. It makes the Euier equations, especially the energy equation, be recovered easily. This step is important in the simulation of gas dynamics. (3) In our model, the particle speed c should be chosen appropriately to meet the requirement of numerical stability. From the viewpoint of mechanics, the larger c allows the larger range of macroscopic speed u. This model preserves the main advantages of the available LBM model: noise-free, simplicity, and high parallelism, etc. The drawback of this model is too many parameters to be chosen, but, it also give us more possibilities to choose. Acknowledgements: The authors are indebted to Professor Li Ronghua, Huang Mingyou and Jin Xizhou for helpful discussions. The second author wishes to express sincere thanks to Professor Chen Yaosong for his instructions. REFERENCES 1 Frisch U, Hasslacher B, Pomeau Y. Lattice-gas automata for the Navier-Stokes equation. Phys Rev Lett, 1986, 56:1505~1508 2 Qian YH, d'Humi~res D, Lallemand P. Lattice BGK models for Navier-Stokes equation. Euro. phys Left, 1992, 17(6): 479~484 3 Chen HD, Chen SY and Matthaeus WH. Recovery of the Navier-Stokes equations using a lattice-gas Boltzmann method. Phys Rev A, 1992, 45:5339,,~5342 4 Benzi R, Succi S and Vergassola M. The lattice Boltzmann equation: theory and applications. Physics Reports, 1992, 222(3): 147~196 5 Yan GW,Shao X, Hu SX. Investigations of waves and experiments of its viscosity with LBM. Acta Scicntiarum Naturalium Univcrsitatis Jilinensis,1996, 118(4): 23~,,27 (in Chinese) 6 Alexander F J, Chen HD, Cheu SY et al. Lattice Boltzmann model for compressible fluids. Phys Rev A, 1992, 46:1967~1970 7 Qian YH and Orszag SA. Lattice BGK model for the Navier-Stokes equation: Nonlinear deviation in compressible regimes. Europhys Left, 1993, 21(3): 255~259 8 Qian YH and Orszag SA. Numerical simulation of weakly compressible Kolmogorov flow with kinetic model. 1994, preprint 9 Ancona MG. FuUy-Lagrangian and Lattice-Boltzmann methods for solving systems of conservation equation. J Comput Phys, 1994, 115:107~120 10 Yan Guangwu. Lattice Boltzmann equation for fluid dynamics. Ph D thesis, Mechanics Dep, Peking University, 1994
226
ACTA MECHANICA SINICA (English Series)
1997
11 Sod GA. A survey of several finite difference methods for systems of non linear hyperbolic conservation laws. J Comput Phys, 197'8, 27:1~31 12 Alexander F J, Chen SY, Sterling JD. Lattice Boltzmann thermohydrodynamics. Physical Rev, E, 1993, 4(4):2249~2252 13 McNamara GR, Garcia AL, Alder BJ. Stabilization of thermal lattice Boltzmann model. J Star Phys, 1995, 27(1/2): 395~408 14 Chen Y, Ohashi H, Akiyama M. Heat transfer in lattice BGK modeled fluid. J Star Phys, 1995,
s1(1/2):
n~s5
15 Shi WP, Hu SX, Yah GW. Lattice Boltzmann equation of multi-speed for compressible fluids. Acta Scien~iarum Naturalium Universitatis Jilinensis, 1996, 116(2): 1~12 (in Chinese) 16 Bernardin D, Sero-Guillaume OE, Sun CH. Multispecies 2D lattice g~s with energy levels: Diffusive properties. Physica, 1991, 47D: 169,~188