- Email: [email protected]
A moment conservation-based non-free parameter compressible lattice Boltzmann model and its application for flux evaluation at cell interface
Computers & Fluids 79 (2013) 190–199
Contents lists available at SciVerse ScienceDirect
Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d
A moment conservation-based non-free parameter compressible lattice Boltzmann model and its application for flux evaluation at cell interface L.M. Yang a, C. Shu b,⇑, J. Wu a a b
Department of Aerodynamics, College of Aerospace Engineering, Nanjing University of Aeronautics and Astronautics, Yudao Street, Nanjing 210016, Jiangsu, China Department of Mechanical Engineering, National University of Singapore, 10 Kent Ridge Crescent, Singapore 119260, Singapore
a r t i c l e
i n f o
Article history: Received 19 December 2012 Received in revised form 5 March 2013 Accepted 26 March 2013 Available online 9 April 2013 Keywords: LBM Conservation forms of moments Compressible flows Euler equations Non-free parameter D1Q4 model FV-LBM
a b s t r a c t Based on the idea of constructing equilibrium distribution functions directly from the conservation forms of moments, a platform for developing non-free parameter lattice Boltzmann models is presented in this work. It is found that the existing compressible lattice Boltzmann models such as D1Q4L2 and D1Q5L2 models of Qu et al. [28,29] and D1Q5 model of Kataoka and Tsutahara [22] can be derived by the platform. This paper goes further to determine the lattice velocities of a non-free parameter D1Q4 model by incorporating two additional higher order conservation forms of moments. Since the lattice velocities are determined physically rather than specified artificially, the non-free parameter D1Q4 model can be applied to simulate compressible flows with a wide range of Mach numbers. The developed non-free parameter D1Q4 model is then applied to a local Riemann problem at the cell interface to establish a new Riemann flux solver for the solution of Euler equations by the finite volume method (FVM). Some test problems, such as Sod shock tube, shock reflection, compressible flows around NACA0012 airfoil, hypersonic flows around a blunt body and double Mach reflection, are simulated to illustrate the capability of present solver. Numerical results show that the present non-free parameter D1Q4 model can provide accurate results with faster convergence rate. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction As an alternative computational fluid dynamics (CFD) approach, the lattice Boltzmann method (LBM) has received an increasing attention in recent years [1–3]. Compared with traditional CFD methods, LBM has many unique advantages: simple scheme, linear convective term, and easy parallelism. Due to these advantages, LBM has been successfully applied to simulate various flow problems. Hitherto, the applications of LBM include the multiphase flow [4–9], porous media flow [10,11], particle suspension flow [12], magneto hydrodynamics [13], and chemical reaction flow [14], etc. However, most of developed lattice Boltzmann models are limited for incompressible flow simulation. The main reason is that the Maxwellian distribution function, which is originated from Boltzmann equation, is difficult to manipulate due to its complication. Usually, the equilibrium distribution function used in LBM is just simply derived by applying the truncated Taylor series expansion to Maxwellian function in terms of Mach number, which inevitably limits the range of Mach number [15–19]. Recently, Hu et al. [20] and Yan et al. [21] developed two two-dimensional (2D)
⇑ Corresponding author. Tel.: +65 6516 6476; fax: +65 6779 1459. E-mail address: [email protected] (C. Shu). 0045-7930/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compfluid.2013.03.020
models. One is the 13-bit with three rest energy levels, and the other is the 17-bit with two rest energy levels. Both of them can recover Euler equations with the streaming–collision process. Although some compressible flow cases with shock waves and contact discontinuities were successfully simulated, the drawback of these two models is obvious. There are a number of free parameters in the models, which have to be tuned carefully to make simulations be stable. Kataoka and Tsutahara [22,23] also presented several models for the compressible Euler and Navier–Stokes (NS) equations. The traditional finite difference scheme (FDS) with the first-order temporal discretization and the second-order upwind spatial discretization of the kinetic equation were applied to solve the differential form of the lattice Boltzmann equation (LBE). The drawback of these models is that the numerical instability occurs when the local Mach number exceeds 1. Besides, Sun [24–26] and Sun and Hsu [27] proposed the adaptive LBM, in which the pattern of lattice velocities varies with the mean flow velocity and the internal energy. Their method has been successfully applied to simulate some subsonic flows and supersonic flows with strong shock waves. Also aiming at the simulation of compressible flows, Qu et al. [28,29] recently presented several models by applying a new deriving method. At first, the Maxwellian distribution function is replaced by a circular function and consequently the integral in the infinite domain of velocity space is changed to the line
L.M. Yang et al. / Computers & Fluids 79 (2013) 190–199
integral along the circle. After that, the continuous circular function is distributed to discrete points. In such a way, the conservation of moments is kept when the line integral is replaced by the weighted sum of functional values at discrete points. The functional value at the discrete velocity point is in fact the equilibrium distribution function. Inspired by this idea, Qu et al. [28,29] developed the D1Q4L2 model, D1Q5L2 model and D2Q13L2 model. By using the proposed models, some supersonic flows with strong shock waves were simulated successfully. However, since the equilibrium distribution function in these models is derived by distributing the circular function along the lattice velocity directions via the Lagrange interpolation, it would inevitably contain some interpolation errors, even though the circular function satisfies the conservation forms of moments. Besides, in these models, the lattice velocities are artificially specified. With different lattice velocities, the feasible range of Mach number may be quite different. Other interesting contributions in the development of compressible lattice Boltzmann models include the work of McNamara et al. [30] and Dellar [31]. In the work of McNamare et al. [30], a system of 21 conservation equations is used to determine the equilibrium distribution functions of thermal lattice Boltzmann model. Dellar [31] derived two one-dimensional models, unsplit seven velocity model and split 4 + 3 velocity model, from the solution of seven conservation equations. Similar to the work of Qu et al. [28,29], in the models of McNamare et al. [30] and Dellar [31], the lattice velocities are still needed to be specified artificially. Inspired by the idea of constructing lattice Boltzmann models directly from the conservation forms of moments, we propose a platform to determine both the equilibrium distribution functions and the lattice velocities without artificial specification for simulation of compressible flows [32]. In the LBM framework, we start from the conservation forms of moments, i.e., the macroscopic conservative variables and their moments are written as a weighted sum of equilibrium distribution functions. Then the equilibrium distribution functions are considered as unknowns and are solved from the conservation forms of moments. Besides, the lattice velocities can be determined explicitly by incorporating additional higher order conservation forms of moments [33,34]. As a platform, the proposed approach is convenient to develop new lattice Boltzmann models as long as a well-posed equation system of momentum conservation forms is provided. In fact, it was found that D1Q4L2 and D1Q5L2 models of Qu et al. [28,29] and D1Q5 model of Kataoka and Tsutahara [22] are special cases of the present approach. To simulate compressible flows with shock waves, the developed one-dimensional nonfree parameter lattice Boltzmann model is not applied globally in the whole domain as done by the model of Kataoka and Tsutahara [22]. Instead, it is applied to a local Riemann problem at the cell interface in the solution of Euler equations by the finite volume method. In other words, the developed new lattice Boltzmann model is actually applied locally to establish a Riemann flux solver. To validate the developed model, some test cases, such as Sod shock tube, shock reflection, compressible flows around NACA0012 airfoil, hypersonic flows around a blunt body and double Mach reflection, are simulated. As compared with the data in the literature, good agreements are achieved.
@ q @ qua þ ¼0 @t @xa @ qua @ qua ub þ pdab þ ¼0 @t @xb @ qE @ qðE þ pÞua þ ¼0 @t @xa
191
ð1Þ
where
p ¼ qRT E¼
ð2Þ
1 RT ua ua þ 2 c1
ð3Þ
t is the time and xa is the spatial coordinate. q, ua, p, T and E are the density, velocity in the a direction, pressure, temperature and total energy of the mean flow, respectively. R is the specific gas constant and c is the specific heat ratio. Note that the summation convention is applied to the subscript a and b. There are three conservative quantities and three convective fluxes in Eq. (1). They are the density q, momentum (density flux) qua, momentum flux quaub + pdab, energy qE, energy flux q(E + p)ua. To recover Euler equations by the lattice Boltzmann model, we require the relations between these conservative quantities and convective fluxes and the equilibrium density distribution functions. Let ni,a be the lattice velocity in the xa direction of the ith particle, and k be the particle potential energy introduced to control the specific heat ratio. According to the work of [22,24,28], the equilibrium density distribution functions gi must satisfy the flowing relations in order to recover the Euler equations,
q¼
I X gi i¼1
qua ¼
I X g i ni;a i¼1
qua ub þ pdab ¼
I X g i ni;a ni;b
ð4Þ
i¼1
I X 1 gi ni;a ni;a þ k 2 2 i¼1 I X q 1 ½ua ua þ ðb þ 2ÞRTub ¼ gi ni;a ni;a þ k ni;b 2 2 i¼1
q
ðua ua þ bRTÞ ¼
Here, I denotes the total number of discrete lattice velocities, b is a given constant representing the number of degrees of freedom of molecules. For monatomic gas and diatomic gas, b is taken as 3 and 5 respectively. Equation system (4) actually constitutes conservation forms of moments of the equilibrium distribution functions from the zeroth order to the third order. Thanks to the fourth relation which is in charge of the energy, we have,
1 RT ðua ua þ bRTÞ ¼ q ua ua þ 2 2 c1
q
ð5Þ
From Eq. (5), we can get the expression of
c¼
bþ2 b
ð6Þ
2.1. Conservation forms of moments for recovering Euler equations
For diatomic gas, c = 1.4. By substituting the first and the third relations into the fourth relation of Eq. (4) and using Eq. (6), we can finally obtain the relation between the particle potential energy and the potential energy of the mean flow as
In this work, only the inviscid compressible flow simulation is considered. The conventional Euler equations can be written as
D k ¼ 1 ðc 1Þ e 2
2. Methodology
ð7Þ
192
L.M. Yang et al. / Computers & Fluids 79 (2013) 190–199
where D is the abbreviation of the dimension, D = 1 means onedimensional. e is the potential energy of the mean flow, which is dep fined as e ¼ ðc1Þ q.
Generally, the existent approaches for developing lattice Boltzmann models can be divided into two categories: indirect method [20–29] and direct method [30,31]. The indirect method is usually based on a known function. For example, those methods presented in [20–23] are based on the truncated Maxwellian function, the adaptive lattice Boltzmann method proposed by Sun [24–26] and Sun and Hsu [27] are based upon a simple delta function, and models of Qu et al. [28,29] are constructed from the circular function. In contrast, the direct method [30,31] derives equilibrium distribution functions directly from the conservation forms of moments. In this work, we follow the idea of the direct method to derive equilibrium distribution functions and lattice velocities in the compressible LBM. It is noted that, different from the direct method of [30,31], the particle potential energy in this work is considered to be independent of lattice velocity direction as it is the result of internal motion (rotational velocity which is independent of translational velocity). The derivation procedure of present method is described as follows. For one-dimensional (1D) case, Eq. (4) contains five relations. The distribution of the discrete lattice velocities for 1D case is shown in Fig. 1. Then, the corresponding form of Eq. (4) can be written as
ð8Þ
þ 2kðg 1 þ g 2 þ g 3 þ g 4 þ g 5 Þ
q½u2 þ ðb þ 2ÞRTu ¼ g 2 d31 g 3 d31 þ g 4 d32 g 5 d32 þ 2kðg 2 d1 g 3 d1 þ g 4 d2 g 5 d2 Þ
ð9Þ
g 1 ¼ g 1 ðq; u; cÞ 2
g3 ¼
2
ðq g 1 Þd1 d2 qd2 u þ qd1 u2 þ qd1 c2 þ qu3 þ 3 quc2 2
2
2d1 ðd1 d2 Þ ðq
2 g 1 Þd1 d2
þq
2 d2 u
þ qd1 u2 þ qd1 c2 qu3 3 quc2 2
2
2d1 ðd1 d2 Þ −d2
− d1
0
d1
2
2
ðq g 1 Þd1 d2 qd1 u qd2 u2 qd2 c2 þ qu3 þ 3 quc2 2
2
2d2 ðd1 d2 Þ
ð10Þ
qðd21 d22 d21 u2 d21 c2 d22 u2 d22 c2 þ u4 þ 6u2 c2 þ c4 Þ 2 2
d1 d2
ð11Þ
then Eq. (10) gives
g1 ¼
g2 ¼
g3 ¼
qðd21 d22 d21 u2 d21 c2 d22 u2 d22 c2 þ u4 þ 6u2 c2 þ c4 Þ 2 2
d1 d2
qð3d1 uc2 d1 d22 u d22 u2 d22 c2 þ u4 þ 6u2 c2 þ c4 þ d1 u3 Þ 2
2
2
2d1 ðd1 d2 Þ
qð3d1 uc2 þ d1 d22 u d22 u2 d22 c2 þ u4 þ 6u2 c2 þ c4 d1 u3 Þ 2
2
2
2d1 ðd1 d2 Þ
g4 ¼
g5 ¼
qðd21 d2 u þ 3d2 uc2 d21 u2 d21 c2 þ u4 þ 6u2 c2 þ c4 þ d2 u3 Þ 2
2
2
2d2 ðd1 d2 Þ
qðd21 d2 u 3d2 uc2 d21 u2 d21 c2 þ u4 þ 6u2 c2 þ c4 d2 u3 Þ 2
2
2
2d2 ðd1 d2 Þ ð12Þ
Obviously, equilibrium distribution functions in Eq. (12) are exactly the same as those in D1Q5L2 model [28]. In the work of Kataoka and Tsutahara [22], the particle potential energy is defined to be varying with lattice velocity. In fact, in their model, the particle potential energy is taken as zero except for the static particle. For the 1D case, by applying our platform, the equilibrium distribution functions can be derived as
ðb 1Þqc2
g20 2
g2 ¼
2
ðq g 1 Þd1 d2 qd2 u þ qd1 u2 þ qd1 c2 þ qu3 þ ðb þ 2Þquc2 2
2
g3 ¼
2
ðq g 1 Þd1 d2 þ qd2 u þ qd1 u2 þ qd1 c2 qu3 ðb þ 2Þ quc2 2
2
2d1 ðd1 d2 Þ 2
2
ðq g 1 Þd1 d2 þ qd1 u qd2 u2 qd2 c2 qu3 ðb þ 2Þquc2 2
2
2d2 ðd1 d2 Þ 2
g5 ¼
2
2d1 ðd1 d2 Þ
2
ðq g 1 Þd1 d2 qd1 u qd2 u2 qd2 c2 þ qu3 þ ðb þ 2Þquc2 2
2
2d2 ðd1 d2 Þ ð13Þ
We can solve Eq. (9) easily by using the well-known software Matlab or Maple. The final solutions of Eq. (9) are
g2 ¼
g1 ¼
g4 ¼
q g1 ¼ g2 þ g3 þ g4 þ g5 qu ¼ g 2 d1 g 3 d1 þ g 4 d2 g 5 d2 qu2 þ qc2 ¼ g 2 d21 þ g 3 d21 þ g 4 d22 þ g 5 d22 qu3 þ 3qc2 u ¼ g 2 d31 g 3 d31 þ g 4 d32 g 5 d32
2
2d2 ðd1 d2 Þ
If we take g1 the same as the one in D1Q5L2 model of Qu et al. [28], that is
g1 ¼
where g1, g2, g3, g4, g5 are equilibrium density distribution functions associated with lattice velocities of 0, d1, d1, d2, d2 as shown in Fig. 1. Note that Eq. (8) contains five unknowns gi(i = 1–5) and there are only four independent relations since the fourth relation has been used to determine the particle potential energy k. As a consequence, we must regard one of gi as given in order to have a wellposed equation system. pffiffiffiffiffiffiffiffiffi If we set g1 as known and take pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c ¼ Dðc 1Þe (c ¼ p=q for 1D case), Eq. (8) can be rewritten as
2
ðq g 1 Þd1 d2 þ qd1 u qd2 u2 qd2 c2 qu3 3 quc2 2
g5 ¼
2.2. Development of non-free parameter D1Q4 model from conservation forms of moments
q ¼ g1 þ g2 þ g3 þ g4 þ g5 qu ¼ g 2 d1 g 3 d1 þ g 4 d2 g 5 d2 qu2 þ p ¼ g 2 d21 þ g 3 d21 þ g 4 d22 þ g 5 d22 qðu2 þ bRTÞ ¼ g 2 d21 þ g 3 d21 þ g 4 d22 þ g 5 d22
2
g4 ¼
d2
Fig. 1. Distribution of discrete lattice velocities for 1D model.
2 0
where g stands for the particle potential energy of the stationary particle. Once again, equilibrium distribution functions in Eq. (13) derived by our platform are exactly the same as those in the D1Q5 model of Kataoka and Tsutahara [22]. From the above analysis, it is natural to draw a conclusion that the models of Qu et al. [28] and Kataoka and Tsutahara [22] are only the special cases of current platform. As mentioned earlier, Eq. (9) only contains four independent relations. Consequently, a D1Q4 model can be easily developed by applying this method when only four lattice velocities are defined in the 1D case. This can be done by removing the static particle as shown in Fig. 1. Therefore, the general form of D1Q4 model is
193
L.M. Yang et al. / Computers & Fluids 79 (2013) 190–199
g1 ¼ g2 ¼ g3 ¼ g4 ¼
qðd1 d22 d22 u þ d1 u2 þ d1 c2 þ u3 þ 3uc2 Þ 2 2d1 ðd1
2 d2 Þ
qðd1 d22 þ d22 u þ d1 u2 þ d1 c2 u3 3uc2 Þ 2
2
2d1 ðd1 d2 Þ
qðd21 d2 þ d21 u d2 u2 d2 c2 u3 3uc2 Þ 2
ð14Þ
2
2d2 ðd1 d2 Þ
qðd21 d2 d21 u d2 u2 d2 c2 þ u3 þ 3uc2 Þ 2
2
2d2 ðd1 d2 Þ
Here, g1, g2, g3, g4 are associated with lattice velocities of d1, d1, d2, d2. It is indicated that the equilibrium distribution functions in Eq. (14) are exactly the same as those in the D1Q4L2 model of Qu et al. [28,29]. In the work of Qu et al. [28,29], the discrete lattice velocities d1 and d2 are artificially specified, an example of taking d1 as 1 and d2 as 2. However, numerical results show that when adopting the hybrid LBM with traditional finite volume method into flow simulation, this D1Q4L2 model may not work very well. To overcome this difficulty, some other forms of d1 and d2 have been proqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 posed, such as d1 = d2/4 and d2 ¼ 12 u2n þ ðcc1Þ (where un is the normal velocity at the interface). By using these improved forms of d1 and d2, the D1Q4L2 model has been successfully applied to simulate some one- and two-dimensional compressible flows. Nevertheless, the permitted free-stream Mach number of the mean flow is still not very high. Even worse, an unexpected spurious oscillation in the vicinity of the sonic line may be generated for some situations. In order to determinate the discrete lattice velocities d1 and d2 of Eq. (14) reasonably, more physical and mathematical constraints are needed to supplement Eq. (9). In this work, we add two higher order conservation forms of moments [33,34], i.e.,
qu4 þ 6qc2 u2 þ 3qc4 ¼ g 1 d41 þ g 2 d41 þ g 3 d42 þ g 4 d42 qu5 þ 10qc2 u3 þ 15qc4 u ¼ g 1 d51 g 2 d51 þ g 3 d52 g 4 d52 ð15Þ The first expression of Eq. (15) is the same as the constraint needed to recover Navier–Stokes equations, and the second relation is associated with Burnett correction to heat flux. Note that equation systems (9) and (15) provide six independent relations, which can be used to determine six unknowns, g1, g2, g3, g4 and d1, d2. We name this model as non-free parameter D1Q4 model, in which d1 and d2 can be derived as
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi d1 ¼ u2 þ 3c2 4u2 c2 þ 6c4 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi d2 ¼ u2 þ 3c2 þ 4u2 c2 þ 6c4
ð16Þ
Since g1, g2, g3, g4 still satisfy Eq. (9), the solution given by Eq. (14) can still be used in the non-free parameter D1Q4 model. 2.3. Flux evaluation by D1Q4 model in the FV solution of Euler equations In the developed non-free parameter D1Q4 model shown in Section 2.2, the lattice velocities and particle potential energy depend on the local fluid properties. When this model is applied globally in the whole domain as done by D1Q5 model of Kataoka and Tsutahara [22], the Euler equations may not be recovered accurately. To remove this drawback, the developed D1Q4 model is only applied to a local Riemann problem at the cell interface to evaluate inviscid flux in the solution of Euler equations by the finite volume method. According to the feature of local Riemann problem, the upstream and downstream fluid properties are constants. This means that along each streaming direction, the lattice
velocity is always a constant. So, locally, the non-free parameter D1Q4 model works in the same manner (constant lattice velocity) as D1Q5 model of Kataoka and Tsutahara [22]. In the solution of Euler equations by FVM, we need to evaluate the inviscid flux at the cell interface. To do this job by using the developed non-free parameter D1Q4 model, we need to apply it to a local Riemann problem where the equilibrium distribution functions at the interface are divided into two parts, located respectively on two sides of the interface according to the directions of lattice velocities. Then the equilibrium distribution functions across the interface can be determined easily. This process is illustrated in Fig. 2. For the one-dimensional case, Euler equations given from Eq. (1) can be written as
2
3
q qu
@ 6 4 @t q
½u2 þ bRT
2
2
7 @ 6 5þ 4 @x q 2
3
qu qu2 þ p ½u2 þ ðb þ 2ÞRTu
7 5¼0
ð17Þ
According to Eq. (8), the flux vector at the cell interface can be written as
2
N X g i di
3
7 6 7 6 i¼1 7 6 7 6 N X 7 6 6 7 6 7 ~ 7 qu2 þ p F ¼ 4 F2 5 ¼ 4 g i di di 5¼6 7 6 7 6 q 2 i¼1 ½u þ ðb þ 2ÞRTu F3 7 6 2 6X N 1 7 5 4 di g i 2 di di þ k 2
F1
3
2
3
qu
ð18Þ
i¼1
When multi-dimensional problems are considered, the above 1D model needs to be applied along the normal direction of interface. For example, for the 2D case, as shown in Fig. 3, we can use the normal velocity Un to replace u in the 1D model (18). Suppose that the unit normal vector at the interface is (nx, ny) for the 2D case. Then we have the following relationships,
U n ¼ nx u þ ny v ;
U s ¼ nx v ny u
u ¼ nx U n ny U s ;
ð19Þ
v ¼ nx U s þ ny Un
ð20Þ
where Un and Us are the normal and tangential velocity components at the interface, and u and v are the velocity components in the x and y directions respectively. The two-dimensional Euler equation can be written as
~ @~ @W F x @~ Fy þ þ ¼0 @t @x @y
ð21Þ
where
2 6 ~ ¼6 W 6 4 2 6 6 ~ Fx ¼ 6 4
q qu qv q 2
½u2 þ v 2 þ bRT qu
3 7 7 7 5
qu2 þ p quv q
½u2 þ v 2 þ ðb þ 2ÞRTu 2
3
2
7 6 7 ~ 6 7; F y ¼ 6 5 4
3
qv quv qv þ p 2
q
½u2 þ v 2 þ ðb þ 2ÞRTv 2
7 7 7 5
ð22Þ Integrating Eq. (21) over a control volume I gives Nf ~I dW 1X ~ F ni Si ¼ V I i¼1 dt
ð23Þ
~ I is the vector of conservative variables for the control cell where W I, VI and Nf are the volume and the number of faces of the control
194
L.M. Yang et al. / Computers & Fluids 79 (2013) 190–199
g 2L
g4L
g1L
g 4R
g 2R
R
L
R
L
g3L g1R
splitting
g 4R
g 2R
g1L
g3L
g3R
interface
interface
Fig. 2. Streaming process of D1Q4 model at the cell interface.
cell I, respectively, and Si is the area of the ith face of the control cell. The flux vector at the interface is given by
2 6 6 ~ F n ¼ nx~ F x þ ny~ Fy ¼ 6 4
q 2
qU n
3
qU n u þ pnx qU n v þ pny
7 7 7 5
ð24Þ
U n ½u2 þ v 2 þ ðb þ 2ÞRT
In the above equations, the tangential velocity Us is still unknown. One feasible way to calculate the terms of F1Us and F 1 U 2s is given by
F 1 U s ¼ F L1 U Ls þ F R1 U Rs
ð27Þ
F 1 U 2s ¼ F L1 ðU Ls Þ2 þ F R1 ðU Rs Þ2
ð28Þ
L
Using Eq. (20), Eq. (24) can also be written as
3 2 3 qU n qU n 7 6 7 6 6 ðqU n U n þ pÞnx qU n U s ny 7 6 ðqU n U n þ pÞnx 7 ~ 7¼6 7 F n ¼6 6 ðqU n U n þ pÞny þ qU n U s nx 7 6 ðqU n U n þ pÞny 7 5 4 5 4 2 2 2 q q U ½U þ U þ ðb þ 2ÞRT U ½U þ ðb þ 2ÞRT s n n 2 n 2 n 3 2 0 6 qU n U s ny 7 7 6 þ6 ð25Þ 7 4 qU n U s nx 5 2
q 2
U n U 2s
Note that in the last component of flux vector in Eq. (25), we have used the relation of u2 þ v 2 ¼ U 2n þ U 2s (velocity magnitude does not P change). If we define F 1 ¼ qU n ¼ Ni¼1 g i di , and apply one-dimensional results of Eq. (18) to the normal direction, the flux vector at the interface (Eq. (25)) can be finally written as
2
3
N X
g i di 7 6 7 6 i¼1 7 6 7 6 2 3 6X N 7 F1 6 g i di di nx þ F 1 U s ðny Þ 7 7 6F 7 6 7 6 i¼1 6 27 ~ 7 Fn ¼ 6 7 ¼ 6 7 N 4 F3 5 6 7 6 X 7 6 g d d n þ F U ðn Þ y 1 s x i i i 7 6 F4 7 6 i¼1 7 6 N 7 6X 4 2 5 1 1 dg dd þk þ F U i i 2
i i
2
1
ð26Þ
s
i¼1
R
where U s and U s are the tangential velocities at the two sides of the interface. F L1 and F R1 denote the mass fluxes calculated by the equilibrium distribution functions (only consider the components pointing to the interface) at the two sides of the interface. In particular, as the non-free parameter D1Q4 model is applied, we have L L R R F L1 ¼ g L1 d1 þ g L3 d2 and F R1 ¼ g R2 d1 g R4 d2 . Once the flux vector at the interface is obtained from Eq. (26), we can simply solve Eq. (23) by multi-stage Runge–Kutta method. Note that in the present simulation of compressible flows, the LBM is used to calculate the flux vector at the interface, while the FVM is used to discretize the Euler equations. The advantages of LBM and FVM are well combined in present work [35], i.e., efficient calculation of flux vector by LBM and geometric flexibility by FVM. The present work is different from the work of [36,37], where the discrete velocity Boltzmann equation (DVBE) is solved directly. DVBE is a set of partial differential equations. The number of unknowns in DVBE (same as the number of lattice velocities) is much larger than that in the Euler equations. In addition, the time step used for solving DVBE is usually very small due to severe stability condition. On the other hand, it is indicated that the present solver can be considered to fall within the broader class of kinetic flux vector splitting scheme (KFVS) [38–41]. KFVS is similar to the flux vector splitting scheme (FVS), and it constructs the Riemann solver by applying Boltzmann equation. In the conventional KFVS [38–41], the continuous Maxwellian function or its variants are usually used to evaluate the flux. In contrast, the present work evaluates the flux of Euler equations by applying the non-free parameter lattice Boltzmann model. 3. Numerical examples
v Uτ
Un
u
In order to validate the present non-free parameter D1Q4 (nonfree D1Q4) model and its flux evaluation for compressible flow problems, the Sod shock tube, shock reflection, flows around NACA0012 airfoil, flows around a blunt body and double Mach reflection are simulated. Unless otherwise stated, the Venkatakrishnan’s limiter [42,43] is applied to reconstruct conservative variables at the interface in this work. 3.1. Sod shock tube The first case is the Sod shock tube. This is a one-dimensional problem, in which the initial condition is given as
interface Fig. 3. Applying 1D model into 2D cases.
ðqL ; uL ; pL Þ ¼ ð1; 0; 1Þ;
0:5 < x < 0
ðqR ; uR ; pR Þ ¼ ð0:125; 0; 0:1Þ 0 < x < 0:5
ð29Þ
L.M. Yang et al. / Computers & Fluids 79 (2013) 190–199
In this case, we set q0 = 1, L0 = 1 (reference density and reference length). The mesh spacing is chosen as Dx = 1/250 and the time step is taken as Dt = 0.001. To reconstruct conservative variables at the interface, the van Leer’s MUSCL approach [44] is used. The results of non-free D1Q4 model are compared with those of D1Q4L2 model of Qu et al. [28,29] and D1Q5 model of Kataoka and Tsutahara [22], to show the performance and capability of present model for this one-dimensional case. At the same time, the theoretical data and numerical results of Roe scheme are also involved to make comparison. We follow the work of Qu et al. [28,29] to take d1 = 1, d2 = 2 in the D1Q4L2 model, and the work of Kataoka and Tsutahara [22] to take d1 = 1, d2 = 3, g0 = 2 in the D1Q5 model. Fig. 4 shows the computed density, velocity, pressure, and internal energy profiles obtained by the non-free D1Q4 model (dashed lines) at t = 0.22. Also displayed in the figure are the exact solutions (solid lines), the results of D1Q4L2 model of Qu et al. [28,29] (long dash lines), D1Q5 model of Kataoka and Tsutahara [22] (dash dot lines) and Roe scheme (dash double dot lines). It is observed from Fig. 4 that
195
the performances of three models are basically the same. However, in the vicinity of the expansion wave, the results of D1Q4L2 model and D1Q5 model are steeper than the present non-free D1Q4 model. This is not a drawback of the non-free D1Q4 model. Since the incorporated two higher order conservation forms of moments can mimic the real physics, it is easier for the non-free D1Q4 model to ensure the positive density and pressure in the case of strong shocks and discontinuities. Based on the comparison among the exact solution, the results of non-free D1Q4 model and Roe scheme, it is clear that the results of non-free D1Q4 model agree fairly well with those of Roe scheme and can essentially match with the exact solution.
3.2. Shock reflection In order to investigate the effect of lattice velocity values on the performance of lattice Boltzmann models, the shock reflection is simulated. The computational domain with a rectangle of length
Fig. 4. Comparison of density, pressure, velocity, and internal energy profiles for the 1D Sod shock tube.
196
L.M. Yang et al. / Computers & Fluids 79 (2013) 190–199
4 and height 1 is divided into 280 80 cells. The boundary conditions are: a reflecting condition along the bottom boundary, supersonic outflow along the right boundary. The Dirichlet conditions on the left and upper boundaries are given by
ðq; u; v ; pÞjð0;y;tÞ ¼ ð1:0; 2:9; 0:0; 1=1:4Þ ðq; u; v ; pÞjðx;1;tÞ ¼ ð1:69997; 2:61934; 0:50633; 1:52819Þ
ð30Þ
In this case, the lattice velocities of D1Q4L2 model [28,29] d1 and d2 are taken as several different values for comparison. Fig. 5 shows the pressure contours calculated by the non-free D1Q4 model. Fig. 6 shows the pressure profiles at y = 0.5 obtained by the non-free D1Q4 model and D1Q4L2 model with various values of d1 and d2. Also presented in this figure is the result of the adaptive lattice Boltzmann method [24]. Clearly, the results of non-free D1Q4 model and D1Q4L2 model with d1 = 1, d2 = 2 and d1 = 2, d2 = 4 agree fairly well with each other, and they are all steeper than the result reported by Sun [24]. However, the result computed by D1Q4L2 model with d1 = 3, d2 = 6 shows some oscillation in the vicinity of shocks. Additionally, numerical results also show that D1Q4L2 model becomes unstable when the values of lattice velocities exceed the range of d1 = [1, 3], d2 = [2, 6]. Therefore, we can conclude that the non-free D1Q4 model is easier to use for the user than D1Q4L2 model, since there is no any free parameter needed to be chosen artificially in D1Q4 model.
Fig. 6. Comparison of pressure profile for the shock refection.
3.3. Flows around NACA0012 airfoil The third case is the flow around NACA0012 airfoil. For this case, the free-stream Mach number is taken as 0.8 and the angle of attack is chosen as 1.25°. Unstructured grid with cell number of 10382 is used, and its partial view is shown in Fig. 7. When D1Q4L2 model of Qu et al. [28,29] is applied, the computation will diverge if d1 and d2 are still taken as constants throughout the whole domain such as those in Sod shock tube case and shock qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 reflection case. Therefore, d1 = d2/4 and d2 ¼ 12 u2n þ ðcc1Þ are chosen in this test case in order to avoid divergence. The simulation results by non-free D1Q4 model are shown in Figs. 8 and 9 (solid lines). Also presented in the figures are the results of D1Q4L2 model (dashed lines) and Roe scheme (dash dot lines). Obviously, the current results quantitatively agree well with those of D1Q4L2 model. As can be seen from Fig. 9, the non-free D1Q4 model has a faster convergence rate than that of both D1Q4L2 model and Roe scheme. This is an appealing feature of non-free parameter D1Q4 model. In the meantime, the comparison of the lift and drag coefficients (Cl and Cd) for this test case obtained by non-free D1Q4 model and D1Q4L2 model and Roe scheme is made in Table 1 with reference data of Stolcis and Johnston [45]. It can be observed that the results of non-free D1Q4 model and D1Q4L2 model are close each other, and they basically agree well with the reference data. 3.4. Flows around a blunt body To further compare current non-free D1Q4 model with D1Q4L2 model of Qu et al. [28,29], the simulation of compressible flows
Fig. 5. Pressure contours of the shock refection.
Fig. 7. Partial view of the grid around NACA0012 airfoil.
around a circular cylinder with free-stream Mach number of 8 is carried out. Same as the case of flows around NACA0012 airfoil, qin ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 d1 = d2/4 and d2 ¼ 12 u2n þ ðcc1Þ are chosen to avoid divergence for D1Q4L2 model. However, with these choices of d1 and d2, the maximum free-stream Mach number without divergence only reaches 3 for this problem. Alternatively, to improve the Mach number of simulation, d1 = abs(|un| |c|) and d2 = |un| + |c| are chosen when the local Mach number exceeds 1. In this case, a uniform mesh size of 160 40 in the cylindrical coordinate system is used. The computed density, pressure, and Mach number contours by non-free D1Q4 model are shown in Fig. 10 (left side). Also displayed in the figure are the results of D1Q4L2 model (right side). As shown in this figure, an unexpected spurious oscillation appears in the vicinity of the sonic line in the results of D1Q4L2 model. On the contrary, the simulation results by non-free D1Q4 model are smooth in the same region. Additionally, the numerical experiments also show that the present nonfree D1Q4 model can be applied to cases with very high freestream Mach number. The reason may be due to the fact that d1 and d2 in our model are physically chosen rather than artificially specified.
L.M. Yang et al. / Computers & Fluids 79 (2013) 190–199
197
Fig. 8. Pressure coefficient distribution on the NACA0012 airfoil surface.
Fig. 9. Convergence history of non-free D1Q4 model, D1Q4L2 model of Qu et al. [28,29] and Roe scheme for flows around NACA0012 airfoil.
Table 1 Comparison of lift (Cl) and drag (Cd) coefficients for NACA 0012 case. References
Cl
Cd
Stolcis and Johnston [45] Roe scheme D1Q4L2 Non-free D1Q4
0.3397 0.2836 0.3059 0.3041
0.0228 0.0215 0.0236 0.0237
3.5. Double Mach reflection In order to illustrate the capability of present model for dealing with strong shock wave problems, the classic double Mach refection with high pressure ratio (pressure ratio is 116.5) is simulated. This problem is initialized by generating a shock wave diagonally into the reflecting wall. In fact, simulating such a problem is
Fig. 10. Comparison of density, pressure, and Mach number contours for flows around a circular cylinder (non-free D1Q4 (left) and D1Q4L2 (right)) at M1 = 8.
L.M. Yang et al. / Computers & Fluids 79 (2013) 190–199
inflow
inflow
1
outflow
ux = 8.25sin (π / 3) u0
ux = 0
uy = −8.25cos(π / 3) u0
uy = 0
ρ = 8ρ0
ρ = 1.4ρ 0
p = 116.5 p0
outflow
198
p = p0
Ma = 10
π /3 0
inflow
1/6
reflecting
3
Fig. 11. Schematic diagram of the double Mach refection.
4. Conclusions
Fig. 12. Density and pressure contours of double Mach refection.
equivalent to simulating a horizontal shock which encounters a wedge of 30°, as shown in Fig. 11. Uniform mesh size of 360 100 is used. The computed density and pressure contours at t = 0.2 are shown in Fig. 12. These results are in good agreement with available results in the literatures [26,46]. The results of Woodward and Colella [46] are shown in Fig. 13. The complex features such as oblique shocks and triple points are well captured.
Following the idea of developing lattice Boltzmann models directly from the conservation forms of moments, a platform for deriving the non-free parameter model is presented in this paper. In the platform, the conservation forms of moments are considered to form an inverse problem. That is, equilibrium distribution functions associated with lattice velocities are considered as unknowns, which are solved from the conservation forms of moments. For the 1D case, it was found that the 5 conservation forms of moments only provide 4 independent equations, which can only solve for 4 unknowns. It was also found that if the equilibrium distribution function for the static particle (lattice velocity is zero) is taken the same as that in the work of Qu et al. [28,29] and Kataoka and Tsutahara [22], their D1Q5L2 and D1Q5 models can be fully recovered from the platform. To avoid using free parameters such as lattice velocities, this paper further incorporates two additional higher order conservation forms of moments and forms a system of 6 independent equations. The solution of this system leads to the non-free parameter D1Q4 model.
Fig. 13. Density and pressure contours of double Mach refection (computed with the PPMDE scheme [46]).
L.M. Yang et al. / Computers & Fluids 79 (2013) 190–199
For simulation of compressible flows with shock waves, Euler equations are solved by the finite volume method. At the cell interface, the flux is evaluated by the developed non-free parameter D1Q4 model. In other words, the proposed non-free parameter D1Q4 model is only applied locally at the cell interface rather than globally in the whole domain. Numerical experiments, including the Sod shock tube, shock reflection, flows around NACA0012 airfoil, flows around a blunt body and double Mach reflection, show that the compressible transonic flows, supersonic flows and hypersonic flows with strong shock waves can be well simulated by the developed non-free parameter D1Q4 model. In addition, it has a faster convergence rate than the D1Q4L2 model and the popularly-used Roe scheme. Various numerical experiments also show that the present nonfree parameter D1Q4 model can be applied to problems with very high Mach number. The good performance of non-free parameter D1Q4 model may be due to the fact that its lattice velocities are determined physically rather than specified artificially. References [1] Ginzburg I. Truncation errors, exact and heuristic stability analysis of tworelaxation-times lattice Boltzmann schemes for anisotropic advection– diffusion equation. Commun Comput Phys 2012;11:1439–502. [2] Zhang T, Shi BC, Chai ZH, Rong FM. Lattice BGK model for incompressible axisymmetric flows. Commun Comput Phys 2012;11:1569–90. [3] Singh S, Krithivasan S, Karlin IV, Succi S, Ansumali S. Energy conserving lattice Boltzmann models for incompressible flow simulations. Commun Comput Phys 2013;13:603–13. [4] Swift MR, Orlandini E, Osborn WR, Yeomans JM. Lattice Boltzmann simulations of liquid–gas and binary fluid systems. Phys Rev E 1996;54:5041–52. [5] He XY, Chen SY, Zhang RY. A lattice Boltzmann scheme for incompressible multiphase flow and its application in simulation of Rayleigh-Taylor instability. J Comput Phys 1999;152:642–63. [6] Luo LS. Theory of the lattice Boltzmann method: lattice Boltzmann models for nonideal gases. Phys Rev E 2000;62:4982–96. [7] Inamuro T, Ogata T, Tajima S, Konishi N. A lattice Boltzmann method for incompressible two-phase flows with large density differences. J Comput Phys 2004;198:628–44. [8] Huang HB, Krafczyk M, Lu XY. Forcing term in single-phase and Shan-Chentype multiphase lattice Boltzmann models. Phys Rev E 2011;84:046710. [9] Srivastava S, Perlekar P, Biferale L, Sbragaglia M, ten Thije Boonkkamp JHM, Toschi F. A study of fluid interfaces and moving contact lines using the lattice Boltzmann method. Commun Comput Phys 2013;13:725–40. [10] Guo ZL, Zhao TS. Lattice Boltzmann model for incompressible flows through porous media. Phys Rev E 2002;66:036304. [11] Tang GH, Tao WQ, He YL. Gas slippage effect on microscale porous flow using the lattice Boltzmann method. Phys Rev E 2005;72:056301. [12] Ladd AJC, Verberg R. Lattice-Boltzmann simulations of particle-fluid suspensions. J Stat Phys 2001;104:1191–251. [13] Chen SY, Chen HD, Martinez D, Matthaeus W. Lattice Boltzmann model for simulation of magnetohydrodynamics. Phys Rev Lett 1991;67:3776–9. [14] Chen S, Dawson SP, Doolen GD, Janecky DR, Lawniczak A. Lattice methods and their applications to reacting systems. Comput Chem Eng 1995;19:617–46. [15] Qian YH. Simulating thermohydrodynamics with lattice BGK models. J Sci Comput 1993;8:231–42. [16] Alexander FJ, Chen S, Sterling JD. Lattice Boltzmann thermohydrodynamics. Phys Rev E 1993;47:R2249–52. [17] Chen Y, Ohashi H, Akiyama M. Thermal lattice Bhatnagar–Gross–Krook model without nonlinear deviations in macrodynamic equations. Phys Rev E 1994;50:2776–83.
199
[18] Watari M, Tsutahara M. Two-dimensional thermal model of the finitedifference lattice Boltzmann method with high spatial isotropy. Phys Rev E 2003;67:036306. [19] Watari M, Tsutahara M. Possibility of constructing a multispeed Bhatnagar– Gross–Krook thermal model of the lattice Boltzmann method. Phys Rev E 2004;70:016703. [20] Hu SX, Yan GW, Shi WP. A lattice Boltzmann model for compressible perfect gas. Acta Mech Sin – PRC 1997;13:218–26. [21] Yan GW, Chen YS, Hu SX. Simple lattice Boltzmann model for simulating flows with shock wave. Phys Rev E 1999;59:454–9. [22] Kataoka T, Tsutahara M. Lattice Boltzmann method for the compressible Euler equations. Phys Rev E 2004;69:056702. [23] Kataoka T, Tsutahara M. Lattice Boltzmann method for the compressible Navier–Stokes equations with flexible specific-heat ratio. Phys Rev E 2004;69:R035701. [24] Sun CH. Lattice-Boltzmann models for high speed flows. Phys Rev E 1998;58:7283–7. [25] Sun CH. Adaptive lattice Boltzmann model for compressible flows: viscous and conductive properties. Phys Rev E 2000;61:2645–53. [26] Sun CH. Simulations of compressible flows with strong shocks by adaptive lattice Boltzmann model. J Comput Phys 2000;161:70–84. [27] Sun CH, Hsu AT. Three-dimensional lattice Boltzmann model for compressible flows. Phys Rev E 2003;68:016303. [28] Qu K, Shu C, Chew YT. Alternative method to construct equilibrium distribution functions in lattice-Boltzmann method simulation of inviscid compressible flows at high Mach number. Phys Rev E 2007;75:036706. [29] Qu K, Shu C, Chew YT. Simulation of shock-wave propagation with finite volume lattice Boltzmann method. Int J Mod Phys C 2007;18:447–54. [30] McNamara GR, Garcia AL, Alder BJ. Stabilization of thermal lattice Boltzmann models. J Stat Phys 1995;81:395–408. [31] Dellar PJ. Two routes from the Boltzmann equation to compressible flow of polyatomic gases. Prog Comput Fluid Dyn 2008;8:84–96. [32] Yang LM, Shu C, Wu J. Development and comparative studies of three non-free parameter lattice Boltzmann models for simulation of compressible flows. Adv Appl Math Mech 2012;4:454–72. [33] Xu K. Gas-kinetic schemes for unsteady compressible flow simulations. VKI for fluid dynamics lecture series; 1998. [34] Xu K. A gas-kinetic BGK scheme for the Navier–Stocks equations and its connection with artificial dissipation and Godunov method. J Comput Phys 2001;171:289–335. [35] Ji CZ, Shu C, Zhao N. A lattice Boltzmann method-based flux solver and its application to solve shock tube problem. Mod Phys Lett B 2009;23:313–6. [36] Chen HD. Volumetric formulation of the lattice Boltzmann method for fluid dynamics: basic concept. Phys Rev 1998;58:3955–63. [37] Ubertini S, Succi S, Bellal G. Lattice Boltzmann schemes without coordinates. Philos Trans R Soc Lond A 2004;362:1763–71. [38] Mandal JC, Deshpande SM. Kinetic flux vector spitting for Euler equations. Comput Fluids 1994;23:447–78. [39] Coquel F, Helluy P, Schneider J. Second-order entropy diminishing scheme for the Euler equations. Int J Numer Meth Fluids 2006;50:1029–61. [40] Chou SY, Baganoff D. Kinetic flux-vector splitting for the Navier–Stokes equations. J Comput Phys 1997;130:217–30. [41] Liu SH, Xu K. Entropy analysis of kinetic flux vector splitting schemes for the compressible Euler equations. Z Angew Math Phys 2001;52:62–78. [42] Venkatakrishnan V. On the accuracy of limiters and convergence to steady state solutions. AIAA paper 93-0880; 1993. [43] Venkatakrishnan V. Convergence to steady-state solutions of the Euler equations on unstructured grids with limiters. J Comput Phys 1995;118:120–30. [44] van Leer B. Towards the ultimate conservative difference scheme V. A second order sequel to Godunov’s method. J Comput Phys 1979;32:101–36. [45] Stolcis L, Johnston LJ. Solution of the Euler equations on unstructured grids for two-dimensional compressible flow. Aeronaut J 1990;94:181–95. [46] Woodward P, Colella P. The numerical simulation of two-dimensional fluid flow with strong shocks. J Comput Phys 1984;54:115–73.
Contents lists available at SciVerse ScienceDirect
Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d
A moment conservation-based non-free parameter compressible lattice Boltzmann model and its application for flux evaluation at cell interface L.M. Yang a, C. Shu b,⇑, J. Wu a a b
Department of Aerodynamics, College of Aerospace Engineering, Nanjing University of Aeronautics and Astronautics, Yudao Street, Nanjing 210016, Jiangsu, China Department of Mechanical Engineering, National University of Singapore, 10 Kent Ridge Crescent, Singapore 119260, Singapore
a r t i c l e
i n f o
Article history: Received 19 December 2012 Received in revised form 5 March 2013 Accepted 26 March 2013 Available online 9 April 2013 Keywords: LBM Conservation forms of moments Compressible flows Euler equations Non-free parameter D1Q4 model FV-LBM
a b s t r a c t Based on the idea of constructing equilibrium distribution functions directly from the conservation forms of moments, a platform for developing non-free parameter lattice Boltzmann models is presented in this work. It is found that the existing compressible lattice Boltzmann models such as D1Q4L2 and D1Q5L2 models of Qu et al. [28,29] and D1Q5 model of Kataoka and Tsutahara [22] can be derived by the platform. This paper goes further to determine the lattice velocities of a non-free parameter D1Q4 model by incorporating two additional higher order conservation forms of moments. Since the lattice velocities are determined physically rather than specified artificially, the non-free parameter D1Q4 model can be applied to simulate compressible flows with a wide range of Mach numbers. The developed non-free parameter D1Q4 model is then applied to a local Riemann problem at the cell interface to establish a new Riemann flux solver for the solution of Euler equations by the finite volume method (FVM). Some test problems, such as Sod shock tube, shock reflection, compressible flows around NACA0012 airfoil, hypersonic flows around a blunt body and double Mach reflection, are simulated to illustrate the capability of present solver. Numerical results show that the present non-free parameter D1Q4 model can provide accurate results with faster convergence rate. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction As an alternative computational fluid dynamics (CFD) approach, the lattice Boltzmann method (LBM) has received an increasing attention in recent years [1–3]. Compared with traditional CFD methods, LBM has many unique advantages: simple scheme, linear convective term, and easy parallelism. Due to these advantages, LBM has been successfully applied to simulate various flow problems. Hitherto, the applications of LBM include the multiphase flow [4–9], porous media flow [10,11], particle suspension flow [12], magneto hydrodynamics [13], and chemical reaction flow [14], etc. However, most of developed lattice Boltzmann models are limited for incompressible flow simulation. The main reason is that the Maxwellian distribution function, which is originated from Boltzmann equation, is difficult to manipulate due to its complication. Usually, the equilibrium distribution function used in LBM is just simply derived by applying the truncated Taylor series expansion to Maxwellian function in terms of Mach number, which inevitably limits the range of Mach number [15–19]. Recently, Hu et al. [20] and Yan et al. [21] developed two two-dimensional (2D)
⇑ Corresponding author. Tel.: +65 6516 6476; fax: +65 6779 1459. E-mail address: [email protected] (C. Shu). 0045-7930/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compfluid.2013.03.020
models. One is the 13-bit with three rest energy levels, and the other is the 17-bit with two rest energy levels. Both of them can recover Euler equations with the streaming–collision process. Although some compressible flow cases with shock waves and contact discontinuities were successfully simulated, the drawback of these two models is obvious. There are a number of free parameters in the models, which have to be tuned carefully to make simulations be stable. Kataoka and Tsutahara [22,23] also presented several models for the compressible Euler and Navier–Stokes (NS) equations. The traditional finite difference scheme (FDS) with the first-order temporal discretization and the second-order upwind spatial discretization of the kinetic equation were applied to solve the differential form of the lattice Boltzmann equation (LBE). The drawback of these models is that the numerical instability occurs when the local Mach number exceeds 1. Besides, Sun [24–26] and Sun and Hsu [27] proposed the adaptive LBM, in which the pattern of lattice velocities varies with the mean flow velocity and the internal energy. Their method has been successfully applied to simulate some subsonic flows and supersonic flows with strong shock waves. Also aiming at the simulation of compressible flows, Qu et al. [28,29] recently presented several models by applying a new deriving method. At first, the Maxwellian distribution function is replaced by a circular function and consequently the integral in the infinite domain of velocity space is changed to the line
L.M. Yang et al. / Computers & Fluids 79 (2013) 190–199
integral along the circle. After that, the continuous circular function is distributed to discrete points. In such a way, the conservation of moments is kept when the line integral is replaced by the weighted sum of functional values at discrete points. The functional value at the discrete velocity point is in fact the equilibrium distribution function. Inspired by this idea, Qu et al. [28,29] developed the D1Q4L2 model, D1Q5L2 model and D2Q13L2 model. By using the proposed models, some supersonic flows with strong shock waves were simulated successfully. However, since the equilibrium distribution function in these models is derived by distributing the circular function along the lattice velocity directions via the Lagrange interpolation, it would inevitably contain some interpolation errors, even though the circular function satisfies the conservation forms of moments. Besides, in these models, the lattice velocities are artificially specified. With different lattice velocities, the feasible range of Mach number may be quite different. Other interesting contributions in the development of compressible lattice Boltzmann models include the work of McNamara et al. [30] and Dellar [31]. In the work of McNamare et al. [30], a system of 21 conservation equations is used to determine the equilibrium distribution functions of thermal lattice Boltzmann model. Dellar [31] derived two one-dimensional models, unsplit seven velocity model and split 4 + 3 velocity model, from the solution of seven conservation equations. Similar to the work of Qu et al. [28,29], in the models of McNamare et al. [30] and Dellar [31], the lattice velocities are still needed to be specified artificially. Inspired by the idea of constructing lattice Boltzmann models directly from the conservation forms of moments, we propose a platform to determine both the equilibrium distribution functions and the lattice velocities without artificial specification for simulation of compressible flows [32]. In the LBM framework, we start from the conservation forms of moments, i.e., the macroscopic conservative variables and their moments are written as a weighted sum of equilibrium distribution functions. Then the equilibrium distribution functions are considered as unknowns and are solved from the conservation forms of moments. Besides, the lattice velocities can be determined explicitly by incorporating additional higher order conservation forms of moments [33,34]. As a platform, the proposed approach is convenient to develop new lattice Boltzmann models as long as a well-posed equation system of momentum conservation forms is provided. In fact, it was found that D1Q4L2 and D1Q5L2 models of Qu et al. [28,29] and D1Q5 model of Kataoka and Tsutahara [22] are special cases of the present approach. To simulate compressible flows with shock waves, the developed one-dimensional nonfree parameter lattice Boltzmann model is not applied globally in the whole domain as done by the model of Kataoka and Tsutahara [22]. Instead, it is applied to a local Riemann problem at the cell interface in the solution of Euler equations by the finite volume method. In other words, the developed new lattice Boltzmann model is actually applied locally to establish a Riemann flux solver. To validate the developed model, some test cases, such as Sod shock tube, shock reflection, compressible flows around NACA0012 airfoil, hypersonic flows around a blunt body and double Mach reflection, are simulated. As compared with the data in the literature, good agreements are achieved.
@ q @ qua þ ¼0 @t @xa @ qua @ qua ub þ pdab þ ¼0 @t @xb @ qE @ qðE þ pÞua þ ¼0 @t @xa
191
ð1Þ
where
p ¼ qRT E¼
ð2Þ
1 RT ua ua þ 2 c1
ð3Þ
t is the time and xa is the spatial coordinate. q, ua, p, T and E are the density, velocity in the a direction, pressure, temperature and total energy of the mean flow, respectively. R is the specific gas constant and c is the specific heat ratio. Note that the summation convention is applied to the subscript a and b. There are three conservative quantities and three convective fluxes in Eq. (1). They are the density q, momentum (density flux) qua, momentum flux quaub + pdab, energy qE, energy flux q(E + p)ua. To recover Euler equations by the lattice Boltzmann model, we require the relations between these conservative quantities and convective fluxes and the equilibrium density distribution functions. Let ni,a be the lattice velocity in the xa direction of the ith particle, and k be the particle potential energy introduced to control the specific heat ratio. According to the work of [22,24,28], the equilibrium density distribution functions gi must satisfy the flowing relations in order to recover the Euler equations,
q¼
I X gi i¼1
qua ¼
I X g i ni;a i¼1
qua ub þ pdab ¼
I X g i ni;a ni;b
ð4Þ
i¼1
I X 1 gi ni;a ni;a þ k 2 2 i¼1 I X q 1 ½ua ua þ ðb þ 2ÞRTub ¼ gi ni;a ni;a þ k ni;b 2 2 i¼1
q
ðua ua þ bRTÞ ¼
Here, I denotes the total number of discrete lattice velocities, b is a given constant representing the number of degrees of freedom of molecules. For monatomic gas and diatomic gas, b is taken as 3 and 5 respectively. Equation system (4) actually constitutes conservation forms of moments of the equilibrium distribution functions from the zeroth order to the third order. Thanks to the fourth relation which is in charge of the energy, we have,
1 RT ðua ua þ bRTÞ ¼ q ua ua þ 2 2 c1
q
ð5Þ
From Eq. (5), we can get the expression of
c¼
bþ2 b
ð6Þ
2.1. Conservation forms of moments for recovering Euler equations
For diatomic gas, c = 1.4. By substituting the first and the third relations into the fourth relation of Eq. (4) and using Eq. (6), we can finally obtain the relation between the particle potential energy and the potential energy of the mean flow as
In this work, only the inviscid compressible flow simulation is considered. The conventional Euler equations can be written as
D k ¼ 1 ðc 1Þ e 2
2. Methodology
ð7Þ
192
L.M. Yang et al. / Computers & Fluids 79 (2013) 190–199
where D is the abbreviation of the dimension, D = 1 means onedimensional. e is the potential energy of the mean flow, which is dep fined as e ¼ ðc1Þ q.
Generally, the existent approaches for developing lattice Boltzmann models can be divided into two categories: indirect method [20–29] and direct method [30,31]. The indirect method is usually based on a known function. For example, those methods presented in [20–23] are based on the truncated Maxwellian function, the adaptive lattice Boltzmann method proposed by Sun [24–26] and Sun and Hsu [27] are based upon a simple delta function, and models of Qu et al. [28,29] are constructed from the circular function. In contrast, the direct method [30,31] derives equilibrium distribution functions directly from the conservation forms of moments. In this work, we follow the idea of the direct method to derive equilibrium distribution functions and lattice velocities in the compressible LBM. It is noted that, different from the direct method of [30,31], the particle potential energy in this work is considered to be independent of lattice velocity direction as it is the result of internal motion (rotational velocity which is independent of translational velocity). The derivation procedure of present method is described as follows. For one-dimensional (1D) case, Eq. (4) contains five relations. The distribution of the discrete lattice velocities for 1D case is shown in Fig. 1. Then, the corresponding form of Eq. (4) can be written as
ð8Þ
þ 2kðg 1 þ g 2 þ g 3 þ g 4 þ g 5 Þ
q½u2 þ ðb þ 2ÞRTu ¼ g 2 d31 g 3 d31 þ g 4 d32 g 5 d32 þ 2kðg 2 d1 g 3 d1 þ g 4 d2 g 5 d2 Þ
ð9Þ
g 1 ¼ g 1 ðq; u; cÞ 2
g3 ¼
2
ðq g 1 Þd1 d2 qd2 u þ qd1 u2 þ qd1 c2 þ qu3 þ 3 quc2 2
2
2d1 ðd1 d2 Þ ðq
2 g 1 Þd1 d2
þq
2 d2 u
þ qd1 u2 þ qd1 c2 qu3 3 quc2 2
2
2d1 ðd1 d2 Þ −d2
− d1
0
d1
2
2
ðq g 1 Þd1 d2 qd1 u qd2 u2 qd2 c2 þ qu3 þ 3 quc2 2
2
2d2 ðd1 d2 Þ
ð10Þ
qðd21 d22 d21 u2 d21 c2 d22 u2 d22 c2 þ u4 þ 6u2 c2 þ c4 Þ 2 2
d1 d2
ð11Þ
then Eq. (10) gives
g1 ¼
g2 ¼
g3 ¼
qðd21 d22 d21 u2 d21 c2 d22 u2 d22 c2 þ u4 þ 6u2 c2 þ c4 Þ 2 2
d1 d2
qð3d1 uc2 d1 d22 u d22 u2 d22 c2 þ u4 þ 6u2 c2 þ c4 þ d1 u3 Þ 2
2
2
2d1 ðd1 d2 Þ
qð3d1 uc2 þ d1 d22 u d22 u2 d22 c2 þ u4 þ 6u2 c2 þ c4 d1 u3 Þ 2
2
2
2d1 ðd1 d2 Þ
g4 ¼
g5 ¼
qðd21 d2 u þ 3d2 uc2 d21 u2 d21 c2 þ u4 þ 6u2 c2 þ c4 þ d2 u3 Þ 2
2
2
2d2 ðd1 d2 Þ
qðd21 d2 u 3d2 uc2 d21 u2 d21 c2 þ u4 þ 6u2 c2 þ c4 d2 u3 Þ 2
2
2
2d2 ðd1 d2 Þ ð12Þ
Obviously, equilibrium distribution functions in Eq. (12) are exactly the same as those in D1Q5L2 model [28]. In the work of Kataoka and Tsutahara [22], the particle potential energy is defined to be varying with lattice velocity. In fact, in their model, the particle potential energy is taken as zero except for the static particle. For the 1D case, by applying our platform, the equilibrium distribution functions can be derived as
ðb 1Þqc2
g20 2
g2 ¼
2
ðq g 1 Þd1 d2 qd2 u þ qd1 u2 þ qd1 c2 þ qu3 þ ðb þ 2Þquc2 2
2
g3 ¼
2
ðq g 1 Þd1 d2 þ qd2 u þ qd1 u2 þ qd1 c2 qu3 ðb þ 2Þ quc2 2
2
2d1 ðd1 d2 Þ 2
2
ðq g 1 Þd1 d2 þ qd1 u qd2 u2 qd2 c2 qu3 ðb þ 2Þquc2 2
2
2d2 ðd1 d2 Þ 2
g5 ¼
2
2d1 ðd1 d2 Þ
2
ðq g 1 Þd1 d2 qd1 u qd2 u2 qd2 c2 þ qu3 þ ðb þ 2Þquc2 2
2
2d2 ðd1 d2 Þ ð13Þ
We can solve Eq. (9) easily by using the well-known software Matlab or Maple. The final solutions of Eq. (9) are
g2 ¼
g1 ¼
g4 ¼
q g1 ¼ g2 þ g3 þ g4 þ g5 qu ¼ g 2 d1 g 3 d1 þ g 4 d2 g 5 d2 qu2 þ qc2 ¼ g 2 d21 þ g 3 d21 þ g 4 d22 þ g 5 d22 qu3 þ 3qc2 u ¼ g 2 d31 g 3 d31 þ g 4 d32 g 5 d32
2
2d2 ðd1 d2 Þ
If we take g1 the same as the one in D1Q5L2 model of Qu et al. [28], that is
g1 ¼
where g1, g2, g3, g4, g5 are equilibrium density distribution functions associated with lattice velocities of 0, d1, d1, d2, d2 as shown in Fig. 1. Note that Eq. (8) contains five unknowns gi(i = 1–5) and there are only four independent relations since the fourth relation has been used to determine the particle potential energy k. As a consequence, we must regard one of gi as given in order to have a wellposed equation system. pffiffiffiffiffiffiffiffiffi If we set g1 as known and take pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c ¼ Dðc 1Þe (c ¼ p=q for 1D case), Eq. (8) can be rewritten as
2
ðq g 1 Þd1 d2 þ qd1 u qd2 u2 qd2 c2 qu3 3 quc2 2
g5 ¼
2.2. Development of non-free parameter D1Q4 model from conservation forms of moments
q ¼ g1 þ g2 þ g3 þ g4 þ g5 qu ¼ g 2 d1 g 3 d1 þ g 4 d2 g 5 d2 qu2 þ p ¼ g 2 d21 þ g 3 d21 þ g 4 d22 þ g 5 d22 qðu2 þ bRTÞ ¼ g 2 d21 þ g 3 d21 þ g 4 d22 þ g 5 d22
2
g4 ¼
d2
Fig. 1. Distribution of discrete lattice velocities for 1D model.
2 0
where g stands for the particle potential energy of the stationary particle. Once again, equilibrium distribution functions in Eq. (13) derived by our platform are exactly the same as those in the D1Q5 model of Kataoka and Tsutahara [22]. From the above analysis, it is natural to draw a conclusion that the models of Qu et al. [28] and Kataoka and Tsutahara [22] are only the special cases of current platform. As mentioned earlier, Eq. (9) only contains four independent relations. Consequently, a D1Q4 model can be easily developed by applying this method when only four lattice velocities are defined in the 1D case. This can be done by removing the static particle as shown in Fig. 1. Therefore, the general form of D1Q4 model is
193
L.M. Yang et al. / Computers & Fluids 79 (2013) 190–199
g1 ¼ g2 ¼ g3 ¼ g4 ¼
qðd1 d22 d22 u þ d1 u2 þ d1 c2 þ u3 þ 3uc2 Þ 2 2d1 ðd1
2 d2 Þ
qðd1 d22 þ d22 u þ d1 u2 þ d1 c2 u3 3uc2 Þ 2
2
2d1 ðd1 d2 Þ
qðd21 d2 þ d21 u d2 u2 d2 c2 u3 3uc2 Þ 2
ð14Þ
2
2d2 ðd1 d2 Þ
qðd21 d2 d21 u d2 u2 d2 c2 þ u3 þ 3uc2 Þ 2
2
2d2 ðd1 d2 Þ
Here, g1, g2, g3, g4 are associated with lattice velocities of d1, d1, d2, d2. It is indicated that the equilibrium distribution functions in Eq. (14) are exactly the same as those in the D1Q4L2 model of Qu et al. [28,29]. In the work of Qu et al. [28,29], the discrete lattice velocities d1 and d2 are artificially specified, an example of taking d1 as 1 and d2 as 2. However, numerical results show that when adopting the hybrid LBM with traditional finite volume method into flow simulation, this D1Q4L2 model may not work very well. To overcome this difficulty, some other forms of d1 and d2 have been proqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 posed, such as d1 = d2/4 and d2 ¼ 12 u2n þ ðcc1Þ (where un is the normal velocity at the interface). By using these improved forms of d1 and d2, the D1Q4L2 model has been successfully applied to simulate some one- and two-dimensional compressible flows. Nevertheless, the permitted free-stream Mach number of the mean flow is still not very high. Even worse, an unexpected spurious oscillation in the vicinity of the sonic line may be generated for some situations. In order to determinate the discrete lattice velocities d1 and d2 of Eq. (14) reasonably, more physical and mathematical constraints are needed to supplement Eq. (9). In this work, we add two higher order conservation forms of moments [33,34], i.e.,
qu4 þ 6qc2 u2 þ 3qc4 ¼ g 1 d41 þ g 2 d41 þ g 3 d42 þ g 4 d42 qu5 þ 10qc2 u3 þ 15qc4 u ¼ g 1 d51 g 2 d51 þ g 3 d52 g 4 d52 ð15Þ The first expression of Eq. (15) is the same as the constraint needed to recover Navier–Stokes equations, and the second relation is associated with Burnett correction to heat flux. Note that equation systems (9) and (15) provide six independent relations, which can be used to determine six unknowns, g1, g2, g3, g4 and d1, d2. We name this model as non-free parameter D1Q4 model, in which d1 and d2 can be derived as
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi d1 ¼ u2 þ 3c2 4u2 c2 þ 6c4 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi d2 ¼ u2 þ 3c2 þ 4u2 c2 þ 6c4
ð16Þ
Since g1, g2, g3, g4 still satisfy Eq. (9), the solution given by Eq. (14) can still be used in the non-free parameter D1Q4 model. 2.3. Flux evaluation by D1Q4 model in the FV solution of Euler equations In the developed non-free parameter D1Q4 model shown in Section 2.2, the lattice velocities and particle potential energy depend on the local fluid properties. When this model is applied globally in the whole domain as done by D1Q5 model of Kataoka and Tsutahara [22], the Euler equations may not be recovered accurately. To remove this drawback, the developed D1Q4 model is only applied to a local Riemann problem at the cell interface to evaluate inviscid flux in the solution of Euler equations by the finite volume method. According to the feature of local Riemann problem, the upstream and downstream fluid properties are constants. This means that along each streaming direction, the lattice
velocity is always a constant. So, locally, the non-free parameter D1Q4 model works in the same manner (constant lattice velocity) as D1Q5 model of Kataoka and Tsutahara [22]. In the solution of Euler equations by FVM, we need to evaluate the inviscid flux at the cell interface. To do this job by using the developed non-free parameter D1Q4 model, we need to apply it to a local Riemann problem where the equilibrium distribution functions at the interface are divided into two parts, located respectively on two sides of the interface according to the directions of lattice velocities. Then the equilibrium distribution functions across the interface can be determined easily. This process is illustrated in Fig. 2. For the one-dimensional case, Euler equations given from Eq. (1) can be written as
2
3
q qu
@ 6 4 @t q
½u2 þ bRT
2
2
7 @ 6 5þ 4 @x q 2
3
qu qu2 þ p ½u2 þ ðb þ 2ÞRTu
7 5¼0
ð17Þ
According to Eq. (8), the flux vector at the cell interface can be written as
2
N X g i di
3
7 6 7 6 i¼1 7 6 7 6 N X 7 6 6 7 6 7 ~ 7 qu2 þ p F ¼ 4 F2 5 ¼ 4 g i di di 5¼6 7 6 7 6 q 2 i¼1 ½u þ ðb þ 2ÞRTu F3 7 6 2 6X N 1 7 5 4 di g i 2 di di þ k 2
F1
3
2
3
qu
ð18Þ
i¼1
When multi-dimensional problems are considered, the above 1D model needs to be applied along the normal direction of interface. For example, for the 2D case, as shown in Fig. 3, we can use the normal velocity Un to replace u in the 1D model (18). Suppose that the unit normal vector at the interface is (nx, ny) for the 2D case. Then we have the following relationships,
U n ¼ nx u þ ny v ;
U s ¼ nx v ny u
u ¼ nx U n ny U s ;
ð19Þ
v ¼ nx U s þ ny Un
ð20Þ
where Un and Us are the normal and tangential velocity components at the interface, and u and v are the velocity components in the x and y directions respectively. The two-dimensional Euler equation can be written as
~ @~ @W F x @~ Fy þ þ ¼0 @t @x @y
ð21Þ
where
2 6 ~ ¼6 W 6 4 2 6 6 ~ Fx ¼ 6 4
q qu qv q 2
½u2 þ v 2 þ bRT qu
3 7 7 7 5
qu2 þ p quv q
½u2 þ v 2 þ ðb þ 2ÞRTu 2
3
2
7 6 7 ~ 6 7; F y ¼ 6 5 4
3
qv quv qv þ p 2
q
½u2 þ v 2 þ ðb þ 2ÞRTv 2
7 7 7 5
ð22Þ Integrating Eq. (21) over a control volume I gives Nf ~I dW 1X ~ F ni Si ¼ V I i¼1 dt
ð23Þ
~ I is the vector of conservative variables for the control cell where W I, VI and Nf are the volume and the number of faces of the control
194
L.M. Yang et al. / Computers & Fluids 79 (2013) 190–199
g 2L
g4L
g1L
g 4R
g 2R
R
L
R
L
g3L g1R
splitting
g 4R
g 2R
g1L
g3L
g3R
interface
interface
Fig. 2. Streaming process of D1Q4 model at the cell interface.
cell I, respectively, and Si is the area of the ith face of the control cell. The flux vector at the interface is given by
2 6 6 ~ F n ¼ nx~ F x þ ny~ Fy ¼ 6 4
q 2
qU n
3
qU n u þ pnx qU n v þ pny
7 7 7 5
ð24Þ
U n ½u2 þ v 2 þ ðb þ 2ÞRT
In the above equations, the tangential velocity Us is still unknown. One feasible way to calculate the terms of F1Us and F 1 U 2s is given by
F 1 U s ¼ F L1 U Ls þ F R1 U Rs
ð27Þ
F 1 U 2s ¼ F L1 ðU Ls Þ2 þ F R1 ðU Rs Þ2
ð28Þ
L
Using Eq. (20), Eq. (24) can also be written as
3 2 3 qU n qU n 7 6 7 6 6 ðqU n U n þ pÞnx qU n U s ny 7 6 ðqU n U n þ pÞnx 7 ~ 7¼6 7 F n ¼6 6 ðqU n U n þ pÞny þ qU n U s nx 7 6 ðqU n U n þ pÞny 7 5 4 5 4 2 2 2 q q U ½U þ U þ ðb þ 2ÞRT U ½U þ ðb þ 2ÞRT s n n 2 n 2 n 3 2 0 6 qU n U s ny 7 7 6 þ6 ð25Þ 7 4 qU n U s nx 5 2
q 2
U n U 2s
Note that in the last component of flux vector in Eq. (25), we have used the relation of u2 þ v 2 ¼ U 2n þ U 2s (velocity magnitude does not P change). If we define F 1 ¼ qU n ¼ Ni¼1 g i di , and apply one-dimensional results of Eq. (18) to the normal direction, the flux vector at the interface (Eq. (25)) can be finally written as
2
3
N X
g i di 7 6 7 6 i¼1 7 6 7 6 2 3 6X N 7 F1 6 g i di di nx þ F 1 U s ðny Þ 7 7 6F 7 6 7 6 i¼1 6 27 ~ 7 Fn ¼ 6 7 ¼ 6 7 N 4 F3 5 6 7 6 X 7 6 g d d n þ F U ðn Þ y 1 s x i i i 7 6 F4 7 6 i¼1 7 6 N 7 6X 4 2 5 1 1 dg dd þk þ F U i i 2
i i
2
1
ð26Þ
s
i¼1
R
where U s and U s are the tangential velocities at the two sides of the interface. F L1 and F R1 denote the mass fluxes calculated by the equilibrium distribution functions (only consider the components pointing to the interface) at the two sides of the interface. In particular, as the non-free parameter D1Q4 model is applied, we have L L R R F L1 ¼ g L1 d1 þ g L3 d2 and F R1 ¼ g R2 d1 g R4 d2 . Once the flux vector at the interface is obtained from Eq. (26), we can simply solve Eq. (23) by multi-stage Runge–Kutta method. Note that in the present simulation of compressible flows, the LBM is used to calculate the flux vector at the interface, while the FVM is used to discretize the Euler equations. The advantages of LBM and FVM are well combined in present work [35], i.e., efficient calculation of flux vector by LBM and geometric flexibility by FVM. The present work is different from the work of [36,37], where the discrete velocity Boltzmann equation (DVBE) is solved directly. DVBE is a set of partial differential equations. The number of unknowns in DVBE (same as the number of lattice velocities) is much larger than that in the Euler equations. In addition, the time step used for solving DVBE is usually very small due to severe stability condition. On the other hand, it is indicated that the present solver can be considered to fall within the broader class of kinetic flux vector splitting scheme (KFVS) [38–41]. KFVS is similar to the flux vector splitting scheme (FVS), and it constructs the Riemann solver by applying Boltzmann equation. In the conventional KFVS [38–41], the continuous Maxwellian function or its variants are usually used to evaluate the flux. In contrast, the present work evaluates the flux of Euler equations by applying the non-free parameter lattice Boltzmann model. 3. Numerical examples
v Uτ
Un
u
In order to validate the present non-free parameter D1Q4 (nonfree D1Q4) model and its flux evaluation for compressible flow problems, the Sod shock tube, shock reflection, flows around NACA0012 airfoil, flows around a blunt body and double Mach reflection are simulated. Unless otherwise stated, the Venkatakrishnan’s limiter [42,43] is applied to reconstruct conservative variables at the interface in this work. 3.1. Sod shock tube The first case is the Sod shock tube. This is a one-dimensional problem, in which the initial condition is given as
interface Fig. 3. Applying 1D model into 2D cases.
ðqL ; uL ; pL Þ ¼ ð1; 0; 1Þ;
0:5 < x < 0
ðqR ; uR ; pR Þ ¼ ð0:125; 0; 0:1Þ 0 < x < 0:5
ð29Þ
L.M. Yang et al. / Computers & Fluids 79 (2013) 190–199
In this case, we set q0 = 1, L0 = 1 (reference density and reference length). The mesh spacing is chosen as Dx = 1/250 and the time step is taken as Dt = 0.001. To reconstruct conservative variables at the interface, the van Leer’s MUSCL approach [44] is used. The results of non-free D1Q4 model are compared with those of D1Q4L2 model of Qu et al. [28,29] and D1Q5 model of Kataoka and Tsutahara [22], to show the performance and capability of present model for this one-dimensional case. At the same time, the theoretical data and numerical results of Roe scheme are also involved to make comparison. We follow the work of Qu et al. [28,29] to take d1 = 1, d2 = 2 in the D1Q4L2 model, and the work of Kataoka and Tsutahara [22] to take d1 = 1, d2 = 3, g0 = 2 in the D1Q5 model. Fig. 4 shows the computed density, velocity, pressure, and internal energy profiles obtained by the non-free D1Q4 model (dashed lines) at t = 0.22. Also displayed in the figure are the exact solutions (solid lines), the results of D1Q4L2 model of Qu et al. [28,29] (long dash lines), D1Q5 model of Kataoka and Tsutahara [22] (dash dot lines) and Roe scheme (dash double dot lines). It is observed from Fig. 4 that
195
the performances of three models are basically the same. However, in the vicinity of the expansion wave, the results of D1Q4L2 model and D1Q5 model are steeper than the present non-free D1Q4 model. This is not a drawback of the non-free D1Q4 model. Since the incorporated two higher order conservation forms of moments can mimic the real physics, it is easier for the non-free D1Q4 model to ensure the positive density and pressure in the case of strong shocks and discontinuities. Based on the comparison among the exact solution, the results of non-free D1Q4 model and Roe scheme, it is clear that the results of non-free D1Q4 model agree fairly well with those of Roe scheme and can essentially match with the exact solution.
3.2. Shock reflection In order to investigate the effect of lattice velocity values on the performance of lattice Boltzmann models, the shock reflection is simulated. The computational domain with a rectangle of length
Fig. 4. Comparison of density, pressure, velocity, and internal energy profiles for the 1D Sod shock tube.
196
L.M. Yang et al. / Computers & Fluids 79 (2013) 190–199
4 and height 1 is divided into 280 80 cells. The boundary conditions are: a reflecting condition along the bottom boundary, supersonic outflow along the right boundary. The Dirichlet conditions on the left and upper boundaries are given by
ðq; u; v ; pÞjð0;y;tÞ ¼ ð1:0; 2:9; 0:0; 1=1:4Þ ðq; u; v ; pÞjðx;1;tÞ ¼ ð1:69997; 2:61934; 0:50633; 1:52819Þ
ð30Þ
In this case, the lattice velocities of D1Q4L2 model [28,29] d1 and d2 are taken as several different values for comparison. Fig. 5 shows the pressure contours calculated by the non-free D1Q4 model. Fig. 6 shows the pressure profiles at y = 0.5 obtained by the non-free D1Q4 model and D1Q4L2 model with various values of d1 and d2. Also presented in this figure is the result of the adaptive lattice Boltzmann method [24]. Clearly, the results of non-free D1Q4 model and D1Q4L2 model with d1 = 1, d2 = 2 and d1 = 2, d2 = 4 agree fairly well with each other, and they are all steeper than the result reported by Sun [24]. However, the result computed by D1Q4L2 model with d1 = 3, d2 = 6 shows some oscillation in the vicinity of shocks. Additionally, numerical results also show that D1Q4L2 model becomes unstable when the values of lattice velocities exceed the range of d1 = [1, 3], d2 = [2, 6]. Therefore, we can conclude that the non-free D1Q4 model is easier to use for the user than D1Q4L2 model, since there is no any free parameter needed to be chosen artificially in D1Q4 model.
Fig. 6. Comparison of pressure profile for the shock refection.
3.3. Flows around NACA0012 airfoil The third case is the flow around NACA0012 airfoil. For this case, the free-stream Mach number is taken as 0.8 and the angle of attack is chosen as 1.25°. Unstructured grid with cell number of 10382 is used, and its partial view is shown in Fig. 7. When D1Q4L2 model of Qu et al. [28,29] is applied, the computation will diverge if d1 and d2 are still taken as constants throughout the whole domain such as those in Sod shock tube case and shock qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 reflection case. Therefore, d1 = d2/4 and d2 ¼ 12 u2n þ ðcc1Þ are chosen in this test case in order to avoid divergence. The simulation results by non-free D1Q4 model are shown in Figs. 8 and 9 (solid lines). Also presented in the figures are the results of D1Q4L2 model (dashed lines) and Roe scheme (dash dot lines). Obviously, the current results quantitatively agree well with those of D1Q4L2 model. As can be seen from Fig. 9, the non-free D1Q4 model has a faster convergence rate than that of both D1Q4L2 model and Roe scheme. This is an appealing feature of non-free parameter D1Q4 model. In the meantime, the comparison of the lift and drag coefficients (Cl and Cd) for this test case obtained by non-free D1Q4 model and D1Q4L2 model and Roe scheme is made in Table 1 with reference data of Stolcis and Johnston [45]. It can be observed that the results of non-free D1Q4 model and D1Q4L2 model are close each other, and they basically agree well with the reference data. 3.4. Flows around a blunt body To further compare current non-free D1Q4 model with D1Q4L2 model of Qu et al. [28,29], the simulation of compressible flows
Fig. 5. Pressure contours of the shock refection.
Fig. 7. Partial view of the grid around NACA0012 airfoil.
around a circular cylinder with free-stream Mach number of 8 is carried out. Same as the case of flows around NACA0012 airfoil, qin ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 d1 = d2/4 and d2 ¼ 12 u2n þ ðcc1Þ are chosen to avoid divergence for D1Q4L2 model. However, with these choices of d1 and d2, the maximum free-stream Mach number without divergence only reaches 3 for this problem. Alternatively, to improve the Mach number of simulation, d1 = abs(|un| |c|) and d2 = |un| + |c| are chosen when the local Mach number exceeds 1. In this case, a uniform mesh size of 160 40 in the cylindrical coordinate system is used. The computed density, pressure, and Mach number contours by non-free D1Q4 model are shown in Fig. 10 (left side). Also displayed in the figure are the results of D1Q4L2 model (right side). As shown in this figure, an unexpected spurious oscillation appears in the vicinity of the sonic line in the results of D1Q4L2 model. On the contrary, the simulation results by non-free D1Q4 model are smooth in the same region. Additionally, the numerical experiments also show that the present nonfree D1Q4 model can be applied to cases with very high freestream Mach number. The reason may be due to the fact that d1 and d2 in our model are physically chosen rather than artificially specified.
L.M. Yang et al. / Computers & Fluids 79 (2013) 190–199
197
Fig. 8. Pressure coefficient distribution on the NACA0012 airfoil surface.
Fig. 9. Convergence history of non-free D1Q4 model, D1Q4L2 model of Qu et al. [28,29] and Roe scheme for flows around NACA0012 airfoil.
Table 1 Comparison of lift (Cl) and drag (Cd) coefficients for NACA 0012 case. References
Cl
Cd
Stolcis and Johnston [45] Roe scheme D1Q4L2 Non-free D1Q4
0.3397 0.2836 0.3059 0.3041
0.0228 0.0215 0.0236 0.0237
3.5. Double Mach reflection In order to illustrate the capability of present model for dealing with strong shock wave problems, the classic double Mach refection with high pressure ratio (pressure ratio is 116.5) is simulated. This problem is initialized by generating a shock wave diagonally into the reflecting wall. In fact, simulating such a problem is
Fig. 10. Comparison of density, pressure, and Mach number contours for flows around a circular cylinder (non-free D1Q4 (left) and D1Q4L2 (right)) at M1 = 8.
L.M. Yang et al. / Computers & Fluids 79 (2013) 190–199
inflow
inflow
1
outflow
ux = 8.25sin (π / 3) u0
ux = 0
uy = −8.25cos(π / 3) u0
uy = 0
ρ = 8ρ0
ρ = 1.4ρ 0
p = 116.5 p0
outflow
198
p = p0
Ma = 10
π /3 0
inflow
1/6
reflecting
3
Fig. 11. Schematic diagram of the double Mach refection.
4. Conclusions
Fig. 12. Density and pressure contours of double Mach refection.
equivalent to simulating a horizontal shock which encounters a wedge of 30°, as shown in Fig. 11. Uniform mesh size of 360 100 is used. The computed density and pressure contours at t = 0.2 are shown in Fig. 12. These results are in good agreement with available results in the literatures [26,46]. The results of Woodward and Colella [46] are shown in Fig. 13. The complex features such as oblique shocks and triple points are well captured.
Following the idea of developing lattice Boltzmann models directly from the conservation forms of moments, a platform for deriving the non-free parameter model is presented in this paper. In the platform, the conservation forms of moments are considered to form an inverse problem. That is, equilibrium distribution functions associated with lattice velocities are considered as unknowns, which are solved from the conservation forms of moments. For the 1D case, it was found that the 5 conservation forms of moments only provide 4 independent equations, which can only solve for 4 unknowns. It was also found that if the equilibrium distribution function for the static particle (lattice velocity is zero) is taken the same as that in the work of Qu et al. [28,29] and Kataoka and Tsutahara [22], their D1Q5L2 and D1Q5 models can be fully recovered from the platform. To avoid using free parameters such as lattice velocities, this paper further incorporates two additional higher order conservation forms of moments and forms a system of 6 independent equations. The solution of this system leads to the non-free parameter D1Q4 model.
Fig. 13. Density and pressure contours of double Mach refection (computed with the PPMDE scheme [46]).
L.M. Yang et al. / Computers & Fluids 79 (2013) 190–199
For simulation of compressible flows with shock waves, Euler equations are solved by the finite volume method. At the cell interface, the flux is evaluated by the developed non-free parameter D1Q4 model. In other words, the proposed non-free parameter D1Q4 model is only applied locally at the cell interface rather than globally in the whole domain. Numerical experiments, including the Sod shock tube, shock reflection, flows around NACA0012 airfoil, flows around a blunt body and double Mach reflection, show that the compressible transonic flows, supersonic flows and hypersonic flows with strong shock waves can be well simulated by the developed non-free parameter D1Q4 model. In addition, it has a faster convergence rate than the D1Q4L2 model and the popularly-used Roe scheme. Various numerical experiments also show that the present nonfree parameter D1Q4 model can be applied to problems with very high Mach number. The good performance of non-free parameter D1Q4 model may be due to the fact that its lattice velocities are determined physically rather than specified artificially. References [1] Ginzburg I. Truncation errors, exact and heuristic stability analysis of tworelaxation-times lattice Boltzmann schemes for anisotropic advection– diffusion equation. Commun Comput Phys 2012;11:1439–502. [2] Zhang T, Shi BC, Chai ZH, Rong FM. Lattice BGK model for incompressible axisymmetric flows. Commun Comput Phys 2012;11:1569–90. [3] Singh S, Krithivasan S, Karlin IV, Succi S, Ansumali S. Energy conserving lattice Boltzmann models for incompressible flow simulations. Commun Comput Phys 2013;13:603–13. [4] Swift MR, Orlandini E, Osborn WR, Yeomans JM. Lattice Boltzmann simulations of liquid–gas and binary fluid systems. Phys Rev E 1996;54:5041–52. [5] He XY, Chen SY, Zhang RY. A lattice Boltzmann scheme for incompressible multiphase flow and its application in simulation of Rayleigh-Taylor instability. J Comput Phys 1999;152:642–63. [6] Luo LS. Theory of the lattice Boltzmann method: lattice Boltzmann models for nonideal gases. Phys Rev E 2000;62:4982–96. [7] Inamuro T, Ogata T, Tajima S, Konishi N. A lattice Boltzmann method for incompressible two-phase flows with large density differences. J Comput Phys 2004;198:628–44. [8] Huang HB, Krafczyk M, Lu XY. Forcing term in single-phase and Shan-Chentype multiphase lattice Boltzmann models. Phys Rev E 2011;84:046710. [9] Srivastava S, Perlekar P, Biferale L, Sbragaglia M, ten Thije Boonkkamp JHM, Toschi F. A study of fluid interfaces and moving contact lines using the lattice Boltzmann method. Commun Comput Phys 2013;13:725–40. [10] Guo ZL, Zhao TS. Lattice Boltzmann model for incompressible flows through porous media. Phys Rev E 2002;66:036304. [11] Tang GH, Tao WQ, He YL. Gas slippage effect on microscale porous flow using the lattice Boltzmann method. Phys Rev E 2005;72:056301. [12] Ladd AJC, Verberg R. Lattice-Boltzmann simulations of particle-fluid suspensions. J Stat Phys 2001;104:1191–251. [13] Chen SY, Chen HD, Martinez D, Matthaeus W. Lattice Boltzmann model for simulation of magnetohydrodynamics. Phys Rev Lett 1991;67:3776–9. [14] Chen S, Dawson SP, Doolen GD, Janecky DR, Lawniczak A. Lattice methods and their applications to reacting systems. Comput Chem Eng 1995;19:617–46. [15] Qian YH. Simulating thermohydrodynamics with lattice BGK models. J Sci Comput 1993;8:231–42. [16] Alexander FJ, Chen S, Sterling JD. Lattice Boltzmann thermohydrodynamics. Phys Rev E 1993;47:R2249–52. [17] Chen Y, Ohashi H, Akiyama M. Thermal lattice Bhatnagar–Gross–Krook model without nonlinear deviations in macrodynamic equations. Phys Rev E 1994;50:2776–83.
199
[18] Watari M, Tsutahara M. Two-dimensional thermal model of the finitedifference lattice Boltzmann method with high spatial isotropy. Phys Rev E 2003;67:036306. [19] Watari M, Tsutahara M. Possibility of constructing a multispeed Bhatnagar– Gross–Krook thermal model of the lattice Boltzmann method. Phys Rev E 2004;70:016703. [20] Hu SX, Yan GW, Shi WP. A lattice Boltzmann model for compressible perfect gas. Acta Mech Sin – PRC 1997;13:218–26. [21] Yan GW, Chen YS, Hu SX. Simple lattice Boltzmann model for simulating flows with shock wave. Phys Rev E 1999;59:454–9. [22] Kataoka T, Tsutahara M. Lattice Boltzmann method for the compressible Euler equations. Phys Rev E 2004;69:056702. [23] Kataoka T, Tsutahara M. Lattice Boltzmann method for the compressible Navier–Stokes equations with flexible specific-heat ratio. Phys Rev E 2004;69:R035701. [24] Sun CH. Lattice-Boltzmann models for high speed flows. Phys Rev E 1998;58:7283–7. [25] Sun CH. Adaptive lattice Boltzmann model for compressible flows: viscous and conductive properties. Phys Rev E 2000;61:2645–53. [26] Sun CH. Simulations of compressible flows with strong shocks by adaptive lattice Boltzmann model. J Comput Phys 2000;161:70–84. [27] Sun CH, Hsu AT. Three-dimensional lattice Boltzmann model for compressible flows. Phys Rev E 2003;68:016303. [28] Qu K, Shu C, Chew YT. Alternative method to construct equilibrium distribution functions in lattice-Boltzmann method simulation of inviscid compressible flows at high Mach number. Phys Rev E 2007;75:036706. [29] Qu K, Shu C, Chew YT. Simulation of shock-wave propagation with finite volume lattice Boltzmann method. Int J Mod Phys C 2007;18:447–54. [30] McNamara GR, Garcia AL, Alder BJ. Stabilization of thermal lattice Boltzmann models. J Stat Phys 1995;81:395–408. [31] Dellar PJ. Two routes from the Boltzmann equation to compressible flow of polyatomic gases. Prog Comput Fluid Dyn 2008;8:84–96. [32] Yang LM, Shu C, Wu J. Development and comparative studies of three non-free parameter lattice Boltzmann models for simulation of compressible flows. Adv Appl Math Mech 2012;4:454–72. [33] Xu K. Gas-kinetic schemes for unsteady compressible flow simulations. VKI for fluid dynamics lecture series; 1998. [34] Xu K. A gas-kinetic BGK scheme for the Navier–Stocks equations and its connection with artificial dissipation and Godunov method. J Comput Phys 2001;171:289–335. [35] Ji CZ, Shu C, Zhao N. A lattice Boltzmann method-based flux solver and its application to solve shock tube problem. Mod Phys Lett B 2009;23:313–6. [36] Chen HD. Volumetric formulation of the lattice Boltzmann method for fluid dynamics: basic concept. Phys Rev 1998;58:3955–63. [37] Ubertini S, Succi S, Bellal G. Lattice Boltzmann schemes without coordinates. Philos Trans R Soc Lond A 2004;362:1763–71. [38] Mandal JC, Deshpande SM. Kinetic flux vector spitting for Euler equations. Comput Fluids 1994;23:447–78. [39] Coquel F, Helluy P, Schneider J. Second-order entropy diminishing scheme for the Euler equations. Int J Numer Meth Fluids 2006;50:1029–61. [40] Chou SY, Baganoff D. Kinetic flux-vector splitting for the Navier–Stokes equations. J Comput Phys 1997;130:217–30. [41] Liu SH, Xu K. Entropy analysis of kinetic flux vector splitting schemes for the compressible Euler equations. Z Angew Math Phys 2001;52:62–78. [42] Venkatakrishnan V. On the accuracy of limiters and convergence to steady state solutions. AIAA paper 93-0880; 1993. [43] Venkatakrishnan V. Convergence to steady-state solutions of the Euler equations on unstructured grids with limiters. J Comput Phys 1995;118:120–30. [44] van Leer B. Towards the ultimate conservative difference scheme V. A second order sequel to Godunov’s method. J Comput Phys 1979;32:101–36. [45] Stolcis L, Johnston LJ. Solution of the Euler equations on unstructured grids for two-dimensional compressible flow. Aeronaut J 1990;94:181–95. [46] Woodward P, Colella P. The numerical simulation of two-dimensional fluid flow with strong shocks. J Comput Phys 1984;54:115–73.