- Email: [email protected]
Supersonic flow simulations by a three-dimensional multispeed thermal model of the finite difference lattice Boltzmann method
ARTICLE IN PRESS
Physica A 364 (2006) 129–144 www.elsevier.com/locate/physa
Supersonic flow simulations by a three-dimensional multispeed thermal model of the finite difference lattice Boltzmann method Minoru Wataria,, Michihisa Tsutaharab a
Wind Tunnel Technology Center, Japan Aerospace Exploration Agency, 7-44-1 Jindaiji Higashi, Chofu, Tokyo 182-8522, Japan b Graduate School of Science and Technology, Kobe University, 1-1 Rokkodai, Nada, Kobe 657-8501, Japan Received 15 December 2004; received in revised form 10 May 2005 Available online 4 October 2005
Abstract A three-dimensional multispeed thermal model of the finite difference lattice Boltzmann method (FDLBM) is proposed. In the FDLBM we can select particle velocities independently from the lattice configuration. Particle velocities of the proposed model consist of vectors pointing to the vertex from the center of a dodecahedron and icosahedron. The model stably performs simulations in a wide range, from subsonic to supersonic flows. To verify the model, simulations of Couette flow, normal shock wave and supersonic nozzle are performed. All results agree with the theoretical predictions. r 2005 Elsevier B.V. All rights reserved. Keywords: Finite difference; Lattice Boltzmann method; Three-dimensional; Couette flow; Supersonic flow; Laval nozzle; Normal shock
1. Introduction The lattice Boltzmann method (LBM) has developed into a promising numerical tool and demonstrates its unique capability in the fluid dynamics [1]. The multispeed thermal model, where several speeds of particle velocities are used, has been proposed to correctly represent heat characteristics and compressibility [2,3]. Ref. [4] has proposed solving the LBM using the finite difference scheme (FDLBM). The FDLBM has three advantages over the LBM. First, the FDLBM can secure the numerical stability by selecting an appropriate time increment. Secondly, it affords a generalized coordinate system that fits with the shape of the boundary. Thirdly, in the LBM, particle velocities are restricted to those that exactly link the lattice nodes in unit time. On the other hand, in the FDLBM as we do not need to consider that constraint, we can select particle velocities independently from the lattice configuration, which enables us to construct a model stable across a wide speed range. Two-dimensional (2-D) and three-dimensional (3-D) multispeed thermal models have been proposed. As far as we know, there have been a few 3-D multispeed thermal models [5–7]. Some of them show limited accuracy when they are used for numerical simulations.
Corresponding author.
E-mail address: [email protected] (M. Watari). 0378-4371/$ - see front matter r 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.physa.2005.06.103
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
130
The 3-D FDLBM model to be proposed in this paper utilizes the vectors pointing to the vertex from the center of a dodecahedron and icosahedron and is derived in a similar way as the 2-D model of our previous study [8]. The speed range that the model stably performs is wide, therefore, a supersonic nozzle simulation, presented in this paper, becomes possible. 2. Finite difference lattice Boltzmann method Below is a general description of the 3-D FDLBM thermal model. The evolution of the distribution function f kl for the particle velocity ckl is governed by the following equation: qf kl qf 1 þ ckla kl ¼ ðf kl f ð0Þ kl Þ, f qt qra
(1)
where the subscript k indicates a group of the particle velocities whose speed is ck and l indicates the particle’s direction. The subscript a indicates an x; y, or z component. The variable t is time, ra is the spatial coordinate, f ð0Þ kl is the local equilibrium distribution function, and f is the relaxation parameter. Macroscopic quantities: the density r, the velocity ua , and the internal energy e, are defined as XX f kl , (2) r¼ k
rua ¼
l
XX k
f kl ckla ,
(3)
l
XX c2 u2 r eþ f kl k . ¼ 2 2 k l
(4)
The local equilibrium distribution function is defined as a polynomial form with regard to the local velocity, which is obtained from the Maxwellian distribution. If appropriate sets of particle velocities are given and the relevant local equilibrium distribution function is defined, the above formulation is shown, by applying the Chapman Enskog expansion, to be equivalent to the following fluid equations, the Navier Stokes equations: qr q þ ðrua Þ ¼ 0. qt qra
(5)
q q q qub qua 2 qug ðrua Þ þ ðPdab þ rua ub Þ m þ dab ¼ 0, qt qrb qrb qra qrb 3 qrg
(6)
q u2 q u2 P q qe qub qua 2 qug r eþ rua e þ þ k0 þ mub þ dab ¼ 0, þ qt qra qra qra 2 2 r qra qrb 3 qrg
(7)
where the pressure P, the viscosity coefficient m, and the heat conductivity k0 have the following relations: P ¼ 23re,
(8)
m ¼ 23ref,
(9)
k0 ¼ 10 9 ref.
(10)
Same expressions: Eqs. (1)–(10), are obtained for nondimensional form if variables are divided as follows, using the reference density r0 , the reference length L, and the reference temperature T 0 (where
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
131
R is gas constant). r; f kl ; f ð0Þ kl ckla ; ua e ra t; f P m; k0
by r0 , pffiffiffiffiffiffiffiffiffiffi by RT 0 , by RT 0 , by L, L by pffiffiffiffiffiffiffiffiffiffi , RT 0 by r0 RT 0 , pffiffiffiffiffiffiffi by r0 L Rt0 .
ð11Þ
The nondimensional form is used for further discussion. The ratio of specific heats g, and the nondimensional ffi pffiffiffiffiffiffiffiffiffi temperature T and speed of sound c, which are normalized by T 0 and RT 0 respectively, are expressed as g ¼ 53,
(12)
2 T ¼ e, 3 pffiffiffiffiffiffi c ¼ gT .
(13) (14)
3. Particle velocities and local equilibrium distribution function To recover the Navier Stokes equations via the Chapman Enskog expansion, the local equilibrium distribution function has to satisfy the following velocity moments (see Appendix): X X ð0Þ f kl ¼ r, (15) k
l
XX k
(17)
2 ¼ r eðua dbg þ ub dga þ ug dab Þ þ ua ub ug , 3
(18)
c2k 10 u2 7 u2 eþ ckla cklb ¼ r e dab þ ua ub e þ . 9 3 2 3 2
(19)
f ð0Þ kl ckla cklb cklg
l
XX k
2 edab þ ua ub , 3
f ð0Þ kl ckla cklb ¼ r
l
XX k
(16)
l
XX k
f ð0Þ kl ckla ¼ rua ,
l
f ð0Þ kl
Eqs. (18) and (19) are requirements for dissipative characteristics. Therefore, the requirements play an important role for a flow where the viscosity is not neglected. As the right-hand side of Eq. (19) contains up to the fourth order of flow velocity u, the local equilibrium distribution function should retain up to the fourth order terms of flow velocity; consequently, the local equilibrium distribution function becomes containing the fourth rank tensor: 3u2 9u4 3 3u2 9 3u2 ð0Þ f kl ¼ rF k 1 1 þ þ cklx ux þ 2 1 cklx cklZ ux uZ 2e 8e 4e 32e2 4e 4e 9 27 þ cklx cklZ cklz ux uZ uz þ cklx cklZ cklz cklw ux uZ uz uw , ð20Þ 16e3 128e4 where F k is the weighting coefficient for the particle velocity ckl .
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
132
The left-hand side of Eq. (18) contains the third rank tensor plus the fourth rank tensor in the f ð0Þ kl . Therefore, it totally contains up to the seventh rank tensor. As a result, the tensors up to the seventh rank should be isotropic. There are two ways in selecting the particle velocity sets to satisfy the tensor requirements. One is to choose different sets of directions lðkÞ for different k and require the isotropy for moments suitably weighted [7]. Another is to find a basic group of velocities that satisfies the seventh rank isotropy and use several speeds of the basic group, as we did for the 2-D model [8]. By adjusting the particle speeds, a wider stable range is expected. Therefore in this paper, we seek the latter. As for being isotropic, the odd rank tensors should vanish and the even rank tensors should be the sum of all possible products of the Kronecker delta [9]. The second, fourth, and sixth rank isotropic tensors of m particle velocities are generally expressed as m X
ckla cklb ¼
l¼1 m X
ckla cklb cklg cklx ¼
l¼1 m X
m 2 c dab , 3 k
(21)
m 4 ð4Þ c D , 15 k abgx
ckla cklb cklg cklx cklZ cklz ¼
l¼1
m 6 ð6Þ c D , 105 k abgxZz
(22)
(23)
ð6Þ where Dð4Þ abgx and DabgxZz are
Dð4Þ abgx ¼ dab dgx þ dag dbx þ dax dbg ,
(24)
ð4Þ ð4Þ ð4Þ ð4Þ ð4Þ Dð6Þ abgxZz ¼ dab DgxZz þ dag DbxZz þ dax DbgZz þ daZ Dbgxz þ daz DbgxZ .
(25)
In the FDLBM we can adopt particle velocities independently from the lattice configuration. Vectors pointing to the vertex from the center of a dodecahedron and icosahedron, shown in Fig. 1, have a good isotropy [9]. Their unit vectors cla are summarized in Table 1. They are isotropic with regard to the odd rank. For the even rank, they are isotropic up to the fourth rank, however, anisotropic for the sixth rank. But if we combine the two groups into one group, the sixth rank tensor becomes almost isotropic, with about 1% or 2% errors, as shown in Table 2. Therefore, these 32 vectors are adopted as the basic group of particle velocities. Moving particle velocities ckla are defined as ckla ¼ ck cla .
(26)
The local equilibrium distribution function (20) is applied to conditions (15)–(19). Considering that the odd rank tensors vanish and the even rank tensors are expressed by Eqs. (21)–(23), the following equations to
Fig. 1. Dodecahedron has 20 vectors pointing to the vertex form the center. Icosahedron has 12 vectors.
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
133
Table 1
pffiffi pffiffi 2 Unit vectors pointing to the vertex from the center of a dodecahedron and icosahedron, where j ¼ 1þ2 5, l ¼ p1ffiffi3, and t ¼ pffiffiffiffiffiffiffiffiffi pffiffi 5þ 5
l : direction
Unit vector (clx ; cly ; clz )
Dodecahedron
l l l l
lð1; 1; 1Þ lð0; j1 ; jÞ lðj; 0; j1 Þ lðj1 ; j; 0Þ
Icosahedron
l ¼ 1–4 l ¼ 5–8 l ¼ 9–12
¼ 1–8 ¼ 9–12 ¼ 13–16 ¼ 17–20
tð0; j; 1Þ tð1; 0; jÞ tðj; 1; 0Þ
Table 2 Typical sixth rank tensors of the unit vectors of a dodecahedron, icosahedron, and the combined group Sixth rank tensor
Dodeca. m ¼ 20
Icosa. m ¼ 12
Dodeca.+Icosa. m ¼ 32
m P
clx clx clx clx clx clx
2:9630 (2.8571)
1:6000 (1.7143)
4:5630 (4.5714)
cly cly clx clx clx clx
0:3529 (0.5714)
0:5789 (0.3429)
0:9318 (0.9143)
cly cly cly cly clx clx
0:6842 (0.5714)
0:2211 (0.3429)
0:9053 (0.9143)
clz clz cly cly clx clx
0:2963 (0.1905)
0:0000 (0.1143)
0:2963 (0.3048)
l¼1 m P l¼1 m P l¼1 m P l¼1
Values in parentheses are those for isotropic tensors calculated by Eq. (23).
determine coefficients F k are derived: XX F k ¼ 1, k
X
F k c2k ¼
1 e, 16
(28)
F k c4k ¼
5 2 e, 24
(29)
F k c6k ¼
35 3 e, 36
(30)
F k c8k ¼
35 4 e. 6
(31)
k
X k
X k
X k
(27)
l
Five speeds of particle velocities are necessary to satisfy the above equations. We assume a rest particle ðc0 ¼ 0:0Þ and four groups of the 32 moving particles, whose speeds are c1 ; c2 ; c3 ; and c4 . Eqs. (27)–(31) are readily solved to give the following: F 0 ¼ 1 32ðF 1 þ F 2 þ F 3 þ F 4 Þ,
(32)
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
134
F1 ¼
F2 ¼
F3 ¼
F4 ¼
1 c21 ðc21 c22 Þðc21 c23 Þðc21 c24 Þ 35 4 35 2 5 2 2 c22 c23 c24 2 2 3 2 2 2 2 2 e ðc2 þ c3 þ c4 Þe þ ðc2 c3 þ c3 c4 þ c4 c2 Þe e , 6 36 24 16
c22 ðc22
1 c24 Þðc22 c21 Þ
c23 Þðc22
35 4 35 2 5 c2 c2 c2 e ðc3 þ c24 þ c21 Þe3 þ ðc23 c24 þ c24 c21 þ c21 c23 Þe2 3 4 1 e . 6 36 24 16
c23 ðc23
ð33Þ
ð34Þ
1 c21 Þðc23 c22 Þ
c24 Þðc23
35 4 35 2 5 c2 c2 c2 e ðc4 þ c21 þ c22 Þe3 þ ðc24 c21 þ c21 c22 þ c22 c24 Þe2 4 1 2 e . 6 36 24 16 1 c24 ðc24 c21 Þðc24 c22 Þðc24 c23 Þ 35 4 35 2 5 2 2 c21 c22 c23 2 2 3 2 2 2 2 2 e ðc1 þ c2 þ c3 Þe þ ðc1 c2 þ c2 c3 þ c3 c1 Þe e . 6 36 24 16
ð35Þ
ð36Þ
As far as the simulation being stably conducted, the values c1 ; c2 ; c3 ; and c4 do not affect the accuracy itself. We can utilize this freedom to obtain a stable range as wide as possible. We investigated the stable range by numerical experiments whether a small density disturbance grows in a uniform flow of speed U and internal energy e, for various combinations of the particle speeds. The stable range depends in a complicated manner on the selection of the particle speeds, the relaxation parameter, and the grid size. The combinations around the following speeds gave the best stable range as far as in the investigation: ðc0 ; c1 ; c2 ; c3 ; c4 Þ ¼ ð0:0; 1:0; 2:0; 3:0; 4:0Þ.
(37)
An example of the stable range for ðc0 ; c1 ; c2 ; c3 ; c4 Þ ¼ ð0:0; 1:0; 2:0; 3:0; 4:0Þ, f ¼ 0:002, and the grid size ðDx; DtÞ ¼ ð0:1; 0:001Þ is shown in Fig. 2. A theoretical stability study for the true optimization is being planned. 4. Verification of the proposed model The proposed thermal model was evaluated by numerical simulations. The evolution Eq. (1) was solved by applying the explicit scheme to time and the second upwind scheme to space. The distribution function of next step f new kl at node ðI; J; KÞ is qf kl qf kl qf kl 1 þ c þ c ¼ f c (38) f new Dt ðf kl f ð0Þ klx kly klz kl kl kl ÞDt, f qx qy qz 8 3f kl;I 4f kl;I1 þ f kl;I2 > > < qf kl 2Dx ¼ 3f 4f > qx kl;I kl;Iþ1 þ f kl;Iþ2 > : 2Dx 8 3f kl;J 4f kl;J1 þ f kl;J2 > > > < 2Dy qf kl ¼ 3f 4f > qy kl;J kl;Jþ1 þ f kl;Jþ2 > > : 2Dy
if cklx X0; if cklx o0:
if ckly X0; if ckly o0:
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
135
5 Mach=3
Mach=2.5 Mach=2 UNSTABLE
4
Mach=1.5 Speed U
3
Mach=1
2 STABLE
Mach=0.5
1
0
1
2 3 Internal energy e
4
5
Fig. 2. Stable range obtained by numerical experiments whether a small disturbance grows in a uniform flow for various U and e, and r ¼ 1. The conditions are ðc0 ; c1 ; c2 ; c3 ; c4 Þ ¼ ð0:0; 1:0; 2:0; 3:0; 4:0Þ, f ¼ 0:002, and ðDx; DtÞ ¼ ð0:1; 0:001Þ.
8 3f 4f kl;K1 þ f kl;K2 > > kl;K qf kl < 2Dz ¼ 3f 4f > qz kl;Kþ1 þ f kl;Kþ2 > : kl;K 2Dz
if cklz X0; if cklz o0;
where the node suffix: I; I 1; ::: are indicated only on the focused coordinate. 4.1. Couette flow simulation Dissipative characteristics were verified by Couette flow. The upper wall, which is H apart from the lower wall and has internal energy e2 , starts to move at speed U. The lower wall has e1 and is at rest. Theoretical internal energy distribution e along vertical axis y in a steady state is given as y m y y 1 e ¼ e1 þ ðe2 e1 Þ þ 0 U 2 . (39) H 2k H H As the ratio m=2k0 is constant (¼ 0:3) from Eqs. (9) and (10), the distribution does not depend on the relaxation parameter f. Fig. 3 shows the simulation result of the proposed model for various relaxation parameters. The result agrees with the theoretical prediction. The same simulations were conducted using the models of Refs. [5–7]. The result for the Ref. [6] model is shown in Fig. 4. The result shows dependence on f, which contradicts the theoretical prediction. We think that this contradiction is because the model retains up to the third order of flow velocity in the equilibrium distribution function and ensures the isotropy only up to the fifth rank: the model does not satisfy constraints (18) and (19) in Section 3. Actually, the 3-D model proposed by us in Ref. [7], which was constructed to overcome the Ref. [6] model’s deficiency by increasing polynomial terms in the equilibrium distribution
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
136
Internal energy subtracted by linear distribution
0.0010
0.0008
0.0006
0.0004
φ =0.05 φ =0.1
0.0002
φ=0.2 φ =0.4 Theory
0.0
0.5 Vertical position y/ H
1.0
Fig. 3. Couette flow simulation using the model proposed in this paper. Internal energy at the steady state for U ¼ 0:1 and e1 ¼ e2 ¼ 1:0. The internal energy subtracted by linear distribution, the last term in Eq. (39), is shown. The relaxation parameter is changed: f ¼ 0:05; 0:1; 0:2; 0:4. The results for all f overlap each other. The number of grids is 161, ðDy; DtÞ ¼ ð0:125; 0:0125Þ.
function and adding several sets of particle speeds to satisfy Eqs. (18) and (19), gives a correct result for the Couette simulation. The 3-D model of Ref. [5] did not give a rational solution: the temperature in the fluid showed a lower value than that of the walls. It seems that the author made a mistake in tensor expressions when deriving the 3-D model (see Ref. [7]). 4.2. Normal shock simulation Compressibility was confirmed by normal shock wave. A uniform supersonic flow: M 1 , r1 , e1 , P1 ð¼ 23 r1 e1 Þ is given at the upstream boundary. The pressure at the downstream boundary P2 is gradually increased. When the pressure P2 becomes a certain value, a stationary flow discontinuity or normal shock wave is generated. The ratios of flow parameters between upstream and downstream of the shock for various upstream Mach numbers M 1 were compared with the theoretical value: Rankine Hugoniot relations in a perfect gas. The result, shown in Fig. 5, agrees with the theory. Regarding the shock structure and the shock thickness, analytical solutions exist only when the shock is weak [10]. In the weak shock theory, the pressure along the x axis, normal to the shock, is expressed as the following, where the position x0 is taken at P ¼ 12 ðP2 þ P1 Þ: x x 1 1 0 P ¼ ðP2 þ P1 Þ þ ðP2 P1 Þ tanh . 2 2 d
(40)
The shock thickness d for monatomic perfect gas with BGK assumption is d¼
8aV 4cm ¼ ðP2 P1 Þðq2 V =qP2 Þs ðP2 P1 Þðg þ 1Þ
2 2 þ ðg 1Þ , 3
where V is the specific volume (¼ 1=r) and a is the sonic absorption parameter.
(41)
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
137
Internal energy subtracted by linear distribution
0.0020
0.0015
0.0010
φ = 0.05
0.0005
φ = 0.1 φ=0.2 φ = 0.4 Theory 0
0.5 Vertical position y/ H
1
Fig. 4. Couette flow simulation using the 3-D model of Ref. [6]. Internal energy at the steady state for U ¼ 0:1 and e1 ¼ e2 ¼ 1:0. The internal energy subtracted by linear distribution is shown. The relaxation parameter is changed: f ¼ 0:05; 0:1; 0:2; 0:4. The result shows dependence on f, which contradicts the theoretical prediction. The grid conditions are same as Fig. 3.
The shock structure and the shock thickness obtained in the simulations for M 1 ¼ 1:02 are shown in Fig. 6 and in Fig. 7, respectively. Both results agree with the theory. 4.3. Supersonic nozzle simulation Finally, to demonstrate the ability to handle a generalized coordinate system, a Laval nozzle that accelerates the flow from Mach 0.2 to Mach 2.0 was simulated. The generalized coordinate system fitted with the nozzle shape is shown in Fig. 8. The generalized coordinate ði; jÞ has the following relations with the cartesian ðx; yÞ. The increment Dx is constant. The increment Dy is a function of x, or i. x ¼ Dx i, y ¼ Dy j, 2 1125 1 ð0:036i þ 0:2Þ2 Dy ¼ 1þ Dx. 5776 0:036i þ 0:2 3 The evolution equation expressed in the generalized coordinate system is qf kl cklx qf kl 1 qDy qf kl 1 ckly cklx j þ þ ¼ ðf kl f ð0Þ kl Þ. Dy qx f qt Dx qi qj
ð42Þ
(43)
The explicit and the second upwind scheme is also applied to this generalized coordinate formulation. Upwind/downwind is judged by the signs of cklx in the i direction and ckly cklx ½qDy=qxj in the j direction.
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
138
3.5 Theory Pressure
Ratio of flow parameters
3.0
P2 / P1
Density ρ2 / ρ1 Int. energy e2 / e1 M2
2.5
2.0
1.5
1.0
0.5 1.0
1.2
1.4
1.6
Upstream Mach M1 Fig. 5. Normal shock simulation. Ratios of flow parameters between upstream and downstream of the shock for f ¼ 0:005 versus various upstream Mach numbers. The number of grids is 101. ðDx; DtÞ ¼ ð0:01; 0:001Þ for M 1 ¼ 1:01 to 1.3, and ðDx; DtÞ ¼ ð0:005; 0:0002Þ for M 1 ¼ 1:5.
1.60 Pressure (weak shock theory) Pressure Density Int. energy Mach
1.05
1.55
1.00
1.50
Density
Pressure, Int. energy, Mach
1.10
0.95
0
1
2
3 Position x
4
5
6
1.45
Fig. 6. Normal shock simulation. Shock structure for M 1 ¼ 1:02 and f ¼ 0:01. The number of grids is 101. ðDx; DtÞ ¼ ð0:04; 0:002Þ.
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
139
1.0 Weak shock theory Simulation
Normal shock thickness
0.8
0.6
0.4
0.2
0.000
0.005 0.010 Relaxation parameter
0.015
Fig. 7. Normal shock simulation. Thickness of the shock for M 1 ¼ 1:02. The relaxation parameter is changed: f ¼ 0:0005; 0:001; 0:002; 0:005; 0:01. The number of grids is 101. The grid sizes is changed from ðDx; DtÞ ¼ ð0:01; 0:0005Þ for f ¼ 0:0005 to ðDx; DtÞ ¼ ð0:04; 0:002Þ for f ¼ 0:01.
y
1 j=10
0 j=0
i=2 -1 j=-10 i=1 i=0 0
i=50
1
2
x
3
4
5
Fig. 8. Supersonic nozzle simulation. Generalized coordinate system fitted with nozzle shape is applied. The grids are 51 21. Dx is constant ð¼ 0:1Þ, Dy is a function of x defined in Eq. (42), and Dt ¼ 0:002.
Analytical solution based on one-dimensional (1-D) isentropic variation predicts the relation of the cross section area A, the internal energy, the density, and the pressure with the Mach number: A 1 2 g 1 2 gþ1=2ðg1Þ 1þ M ¼ , (44) A M g þ 1 2 e ¼ e
g1 g1=g r P 2 g 1 2 1 1þ M ¼ ¼ , r P gþ1 2
(45)
where the subscript indicates the values at the throat. To ensure isentropic variation, sufficiently low viscosity: f ¼ 0:002 was used in the simulation. The upper and lower edges are slip and adiabatic boundaries. Initially, a uniform fluid of r1 ; e1 ; P1 ð¼ 23 r1 e1 Þ in the nozzle is at rest. The density and the internal energy, therefore the pressure at the entrance boundary: r2 ; e2 ; P2 ð¼ 23 r2 e2 Þ, are gradually increased, while keeping the pressure at the exit boundary the initial value
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
140
P1 by adjusting the internal energy. The velocities at both boundaries and the density at the exit are extapolated from the inner fluid region. When the pressure at the entrance boundary becomes sufficiently high, a supersonic flow is exerted in the expansion part of the nozzle. The simulation result of Mach contour lines in the nozzle is shown in Fig. 9. The Mach number, the internal energy, the density, and the pressure are averaged at each section. The ratios of these flow parameters to the values at the throat, versus the averaged Mach number, are shown in Fig. 10. Agreement with the theory of isentropic acceleration is confirmed, though a slight undevelopment of inflative acceleration is observed at the end of the nozzle.
y
1
Mach 2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2
0
-1 0
1
2
3
x
4
5
Fig. 9. Supersonic nozzle simulation. Mach contours: M ¼ 0:2; 0:3; . . . ; 2:0 in the nozzle. Initially uniform fluid of ðr1 ; e1 ; P1 Þ ¼ ð0:286; 0:651; 0:124Þ is at rest in the nozzle. The pressure at the entrance is increased to ðr2 ; e2 ; P2 Þ ¼ ð1:0; 1:5; 1:0Þ.
3.0 Theory
Ratio of flow parameters
2.5
2.0
Area
A / A*
Int. energy
e / e*
Density
ρ / ρ*
Pressure
P / P*
1.5
1.0
0.5
0.0 0.2
0.5
1.0 Mach
1.5
2.0
Fig. 10. Supersonic nozzle simulation. Ratios of the averaged flow parameters: A=A ; e=e ; r=r , and P=P to the averaged Mach number.
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
141
5. Conclusions A thermally correct 3-D multispeed FDLBM model was proposed. To derive thermally correct fluid equations, the local equilibrium distribution function should contain up to the fourth order terms of flow velocity, and tensors made of particle velocities should be isotropic up to the seventh rank. Vectors pointing to the vertex from the center of a dodecahedron and icosahedron have a good isotropy. If these 32 vectors are combined as a group, the above tensor isotropy requirement can be realized. Particle velocities of the model consist of a rest particle and four groups of these vectors. The model was verified by simulations of Couette flow, normal shock wave, and supersonic nozzle. The versatility of the model was demonstrated, from subsonic to supersonic flows, by showing good agreements with the theoretical predictions. A theoretical stability study for the true optimization of the particle speeds is being planned. Appendix A. Chapman Enskog expansion The Chapman Enskog expansion for the FDLBM formulation is presented briefly to show that the formulation is equivalent to Navier Stokes equations if the local equilibrium distribution function satisfies the velocity moments’ condition [2,11]. The FDLBM formulation is qf kl qf 1 þ ckla kl ¼ ðf kl f ð0Þ kl Þ, f qt qra XX f kl , r¼ k
rua ¼
XX
f kl ckla ,
(A.3)
l
XX c2 u2 f kl k . r eþ ¼ 2 2 k l The velocity moments’ condition is X X ð0Þ f kl ¼ r,
l
XX k
(A.6)
c2k u2 ¼r eþ , 2 2
f ð0Þ kl ckla cklb f ð0Þ kl
2 ¼ r edab þ ua ub , 3
c2k 5 u2 ckla ¼ rua e þ , 3 2 2
l
(A.7)
(A.8)
(A.9)
f ð0Þ kl ckla cklb cklg
2 ¼ r eðua dbg þ ub dga þ ug dab Þ þ ua ub ug , 3
(A.10)
c2k ckla cklb f ð0Þ kl
10 u2 7 u2 eþ ¼r e dab þ ua ub e þ . 9 3 3 2
(A.11)
l
XX k
f ð0Þ kl
l
XX k
f ð0Þ kl ckla ¼ rua ,
l
XX k
(A.5)
l
XX k
(A.4)
l
XX k
(A.2)
l
k
k
(A.1)
2
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
142
Actually, since Eqs. (A.7) and (A.9) can be derived by contraction from (A.8) and (A.10), respectively, independent requirements are Eqs. (A.5), (A.6), (A.8), (A.10), and (A.11). We assume that the relaxation parameter f is small and the distribution function f kl is a deviation from the local equilibrium f ð0Þ kl . The distribution function is expanded around the local equilibrium: ð1Þ 2 ð2Þ f kl ¼ f ð0Þ kl þ ff kl þ f f kl þ .
Eq. (A.1) becomes " # " # ð0Þ ð1Þ ð1Þ qf ð0Þ qf qf qf ð1Þ ð2Þ kl kl þ ckla kl þ f kl þ f þ ckla kl þ f kl þ ¼ 0. qt qra qt qra
(A.12)
(A.13)
Using the relations (A.2)–(A.7), the following velocity moments of the perturbation distribution f ðnÞ kl ðn ¼ 1; 2; . . .Þ can be assumed to be zero: X X ðnÞ f kl ¼ 0, (A.14) k
l
XX k
(A.15)
l
XX k
f ðnÞ kl ckla ¼ 0,
l
f ðnÞ kl
c2k ¼ 0. 2
(A.16)
For the zero-th approximation: qf ð0Þ qf ð0Þ kl þ ckla kl þ f ð1Þ (A.17) kl ¼ 0, qt qra P P P P c2 P P taking the velocity moments: k l 1, k l ckla , k l 2k , and considering Eqs. (A.14)–(A.16), we obtain the following equations: q X X ð0Þ q X X ð0Þ f kl þ f ckla ¼ 0, qt k l qra k l kl
(A.18)
q X X ð0Þ q X X ð0Þ f kl ckla þ f ckla cklb ¼ 0, qt k l qrb k l kl
(A.19)
q X X ð0Þ c2k q X X ð0Þ c2k þ ckla ¼ 0. f kl f qt k l 2 qra k l kl 2
(A.20)
Substituting the velocity moments’ condition (A.5)–(A.9) into Eqs. (A.18)–(A.20), we obtain the Euler (compressible and inviscid fluid) equations: qr q þ ðrua Þ ¼ 0, qt qra
(A.21)
q q 2 ðrua Þ þ redab þ rua ub ¼ 0, qt qrb 3
(A.22)
q u2 q u2 2 r eþ ¼ 0. rua e þ þ e þ qt qra 2 2 3
(A.23)
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
143
Next, the first order OðfÞ term in Eq. (A.13) is included. Taking the velocity moments: P P P P c2k P P k l 1; k l ckla ; k l 2 , we obtain the following equations: X X X X q q f ð0Þ þ f ð0Þ ckla ¼ 0, (A.24) qt k l kl qra k l kl q X X ð0Þ q X X ð0Þ q X X ð1Þ f kl ckla þ f kl ckla cklb þ f f ckla cklb ¼ 0, qt k l qrb k l qrb k l kl
(A.25)
q X X ð0Þ c2k q X X ð0Þ c2k q X X ð1Þ c2k þ ckla þ f ckla ¼ 0. f kl f kl f qt k l qra k l kl 2 2 qra k l 2
(A.26)
ð0Þ The f ð1Þ kl is replaced by the f kl expression, using zero-th approximation (A.17): q X X ð0Þ q X X ð0Þ f kl þ f ckla ¼ 0, qt k l qra k l kl
q X X ð0Þ q X X ð0Þ f kl ckla þ f ckla cklb qt k l qrb k l kl q q X X ð0Þ q X X ð0Þ f f kl ckla cklb þ f ckla cklb cklg qrb qt k l qrg k l kl q X X ð0Þ c2k q X X ð0Þ c2k þ ckla f kl f qt k l 2 qra k l kl 2 q q X X ð0Þ c2k q X X ð0Þ c2k ckla þ ckla cklb f f kl f qra qt k l qrb k l kl 2 2
(A.27)
! ¼ 0,
ðA:28Þ
! ¼ 0.
ðA:29Þ
Substituting Eqs. (A.5)–(A.11) into the above equations, we obtain equations expressed by macroscopic quantities. You may be perplexed regarding how to handle time derivatives arising from the terms P P P P ð0Þ 2 q=qt k l f ð0Þ l f kl ½ck =2ckla . These time derivatives are replaced by the following spatial kl ckla cklb and q=qt k derivatives derived from Euler equations: qr qua qr ¼ r ua , (A.30) qt qra qra qua 2 qe 2e qr qua ¼ ub , 3 qra 3r qra qt qrb
(A.31)
qe qe 2e qua ¼ ua , qt qra 3 qra
(A.32)
then, we obtain the Navier Stokes equations: qr q þ ðrua Þ ¼ 0, qt qra q q 2 q 2 qub qua 2 qug ðrua Þ þ redab þ rua ub ref þ dab ¼ 0, qt qrb 3 qrb 3 qra qrb 3 qrg q u2 q u2 2 r eþ rua e þ þ e þ qt qra 2 2 3 q 10 qe 2 qub qua 2 qug ref þ refub þ dab ¼ 0. qra 9 qra 3 qra qrb 3 qrg
(A.33) (A.34)
ðA:35Þ
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
144
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11]
S. Chen, G. Doolen, Annu. Rev. Fluid Mech. 30 (1998) 329. G. McNamara, B. Alder, Physica A 194 (1993) 218. F.J. Alexander, S. Chen, J.D. Sterling, Phys. Rev. E 47 (1993) R2249. N. Cao, S. Chen, S. Jin, D. Martinez, Phys. Rev. E 55 (1997) R21. Y. Chen, H. Ohashi, M. Akiyama, Phys. Rev. E 50 (1994) 2776. N. Takada, Y. Yamakoshi, M. Tsutahara, Numerical analysis of fluid motions by three-dimensional thermal lattice Boltzmann model, JSME B, vol. 64, 1998, pp. 4 (in Japanese). M. Watari, M. Tsutahara, Phys. Rev. E 70 (2004) 16703. M. Watari, M. Tsutahara, Phys. Rev. E 67 (2003) 36306. S. Wolfram, J. Statist. Phys. 45 (1986) 471. L.D. Landau, E.M. Lifshitz, Fluid mechanics, Oxford, 1987. S. Chapman, T.G. Cowling, The Mathematical Theory of Non-uniform Gases, Cambridge University Press, London, 1970.
Physica A 364 (2006) 129–144 www.elsevier.com/locate/physa
Supersonic flow simulations by a three-dimensional multispeed thermal model of the finite difference lattice Boltzmann method Minoru Wataria,, Michihisa Tsutaharab a
Wind Tunnel Technology Center, Japan Aerospace Exploration Agency, 7-44-1 Jindaiji Higashi, Chofu, Tokyo 182-8522, Japan b Graduate School of Science and Technology, Kobe University, 1-1 Rokkodai, Nada, Kobe 657-8501, Japan Received 15 December 2004; received in revised form 10 May 2005 Available online 4 October 2005
Abstract A three-dimensional multispeed thermal model of the finite difference lattice Boltzmann method (FDLBM) is proposed. In the FDLBM we can select particle velocities independently from the lattice configuration. Particle velocities of the proposed model consist of vectors pointing to the vertex from the center of a dodecahedron and icosahedron. The model stably performs simulations in a wide range, from subsonic to supersonic flows. To verify the model, simulations of Couette flow, normal shock wave and supersonic nozzle are performed. All results agree with the theoretical predictions. r 2005 Elsevier B.V. All rights reserved. Keywords: Finite difference; Lattice Boltzmann method; Three-dimensional; Couette flow; Supersonic flow; Laval nozzle; Normal shock
1. Introduction The lattice Boltzmann method (LBM) has developed into a promising numerical tool and demonstrates its unique capability in the fluid dynamics [1]. The multispeed thermal model, where several speeds of particle velocities are used, has been proposed to correctly represent heat characteristics and compressibility [2,3]. Ref. [4] has proposed solving the LBM using the finite difference scheme (FDLBM). The FDLBM has three advantages over the LBM. First, the FDLBM can secure the numerical stability by selecting an appropriate time increment. Secondly, it affords a generalized coordinate system that fits with the shape of the boundary. Thirdly, in the LBM, particle velocities are restricted to those that exactly link the lattice nodes in unit time. On the other hand, in the FDLBM as we do not need to consider that constraint, we can select particle velocities independently from the lattice configuration, which enables us to construct a model stable across a wide speed range. Two-dimensional (2-D) and three-dimensional (3-D) multispeed thermal models have been proposed. As far as we know, there have been a few 3-D multispeed thermal models [5–7]. Some of them show limited accuracy when they are used for numerical simulations.
Corresponding author.
E-mail address: [email protected] (M. Watari). 0378-4371/$ - see front matter r 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.physa.2005.06.103
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
130
The 3-D FDLBM model to be proposed in this paper utilizes the vectors pointing to the vertex from the center of a dodecahedron and icosahedron and is derived in a similar way as the 2-D model of our previous study [8]. The speed range that the model stably performs is wide, therefore, a supersonic nozzle simulation, presented in this paper, becomes possible. 2. Finite difference lattice Boltzmann method Below is a general description of the 3-D FDLBM thermal model. The evolution of the distribution function f kl for the particle velocity ckl is governed by the following equation: qf kl qf 1 þ ckla kl ¼ ðf kl f ð0Þ kl Þ, f qt qra
(1)
where the subscript k indicates a group of the particle velocities whose speed is ck and l indicates the particle’s direction. The subscript a indicates an x; y, or z component. The variable t is time, ra is the spatial coordinate, f ð0Þ kl is the local equilibrium distribution function, and f is the relaxation parameter. Macroscopic quantities: the density r, the velocity ua , and the internal energy e, are defined as XX f kl , (2) r¼ k
rua ¼
l
XX k
f kl ckla ,
(3)
l
XX c2 u2 r eþ f kl k . ¼ 2 2 k l
(4)
The local equilibrium distribution function is defined as a polynomial form with regard to the local velocity, which is obtained from the Maxwellian distribution. If appropriate sets of particle velocities are given and the relevant local equilibrium distribution function is defined, the above formulation is shown, by applying the Chapman Enskog expansion, to be equivalent to the following fluid equations, the Navier Stokes equations: qr q þ ðrua Þ ¼ 0. qt qra
(5)
q q q qub qua 2 qug ðrua Þ þ ðPdab þ rua ub Þ m þ dab ¼ 0, qt qrb qrb qra qrb 3 qrg
(6)
q u2 q u2 P q qe qub qua 2 qug r eþ rua e þ þ k0 þ mub þ dab ¼ 0, þ qt qra qra qra 2 2 r qra qrb 3 qrg
(7)
where the pressure P, the viscosity coefficient m, and the heat conductivity k0 have the following relations: P ¼ 23re,
(8)
m ¼ 23ref,
(9)
k0 ¼ 10 9 ref.
(10)
Same expressions: Eqs. (1)–(10), are obtained for nondimensional form if variables are divided as follows, using the reference density r0 , the reference length L, and the reference temperature T 0 (where
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
131
R is gas constant). r; f kl ; f ð0Þ kl ckla ; ua e ra t; f P m; k0
by r0 , pffiffiffiffiffiffiffiffiffiffi by RT 0 , by RT 0 , by L, L by pffiffiffiffiffiffiffiffiffiffi , RT 0 by r0 RT 0 , pffiffiffiffiffiffiffi by r0 L Rt0 .
ð11Þ
The nondimensional form is used for further discussion. The ratio of specific heats g, and the nondimensional ffi pffiffiffiffiffiffiffiffiffi temperature T and speed of sound c, which are normalized by T 0 and RT 0 respectively, are expressed as g ¼ 53,
(12)
2 T ¼ e, 3 pffiffiffiffiffiffi c ¼ gT .
(13) (14)
3. Particle velocities and local equilibrium distribution function To recover the Navier Stokes equations via the Chapman Enskog expansion, the local equilibrium distribution function has to satisfy the following velocity moments (see Appendix): X X ð0Þ f kl ¼ r, (15) k
l
XX k
(17)
2 ¼ r eðua dbg þ ub dga þ ug dab Þ þ ua ub ug , 3
(18)
c2k 10 u2 7 u2 eþ ckla cklb ¼ r e dab þ ua ub e þ . 9 3 2 3 2
(19)
f ð0Þ kl ckla cklb cklg
l
XX k
2 edab þ ua ub , 3
f ð0Þ kl ckla cklb ¼ r
l
XX k
(16)
l
XX k
f ð0Þ kl ckla ¼ rua ,
l
f ð0Þ kl
Eqs. (18) and (19) are requirements for dissipative characteristics. Therefore, the requirements play an important role for a flow where the viscosity is not neglected. As the right-hand side of Eq. (19) contains up to the fourth order of flow velocity u, the local equilibrium distribution function should retain up to the fourth order terms of flow velocity; consequently, the local equilibrium distribution function becomes containing the fourth rank tensor: 3u2 9u4 3 3u2 9 3u2 ð0Þ f kl ¼ rF k 1 1 þ þ cklx ux þ 2 1 cklx cklZ ux uZ 2e 8e 4e 32e2 4e 4e 9 27 þ cklx cklZ cklz ux uZ uz þ cklx cklZ cklz cklw ux uZ uz uw , ð20Þ 16e3 128e4 where F k is the weighting coefficient for the particle velocity ckl .
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
132
The left-hand side of Eq. (18) contains the third rank tensor plus the fourth rank tensor in the f ð0Þ kl . Therefore, it totally contains up to the seventh rank tensor. As a result, the tensors up to the seventh rank should be isotropic. There are two ways in selecting the particle velocity sets to satisfy the tensor requirements. One is to choose different sets of directions lðkÞ for different k and require the isotropy for moments suitably weighted [7]. Another is to find a basic group of velocities that satisfies the seventh rank isotropy and use several speeds of the basic group, as we did for the 2-D model [8]. By adjusting the particle speeds, a wider stable range is expected. Therefore in this paper, we seek the latter. As for being isotropic, the odd rank tensors should vanish and the even rank tensors should be the sum of all possible products of the Kronecker delta [9]. The second, fourth, and sixth rank isotropic tensors of m particle velocities are generally expressed as m X
ckla cklb ¼
l¼1 m X
ckla cklb cklg cklx ¼
l¼1 m X
m 2 c dab , 3 k
(21)
m 4 ð4Þ c D , 15 k abgx
ckla cklb cklg cklx cklZ cklz ¼
l¼1
m 6 ð6Þ c D , 105 k abgxZz
(22)
(23)
ð6Þ where Dð4Þ abgx and DabgxZz are
Dð4Þ abgx ¼ dab dgx þ dag dbx þ dax dbg ,
(24)
ð4Þ ð4Þ ð4Þ ð4Þ ð4Þ Dð6Þ abgxZz ¼ dab DgxZz þ dag DbxZz þ dax DbgZz þ daZ Dbgxz þ daz DbgxZ .
(25)
In the FDLBM we can adopt particle velocities independently from the lattice configuration. Vectors pointing to the vertex from the center of a dodecahedron and icosahedron, shown in Fig. 1, have a good isotropy [9]. Their unit vectors cla are summarized in Table 1. They are isotropic with regard to the odd rank. For the even rank, they are isotropic up to the fourth rank, however, anisotropic for the sixth rank. But if we combine the two groups into one group, the sixth rank tensor becomes almost isotropic, with about 1% or 2% errors, as shown in Table 2. Therefore, these 32 vectors are adopted as the basic group of particle velocities. Moving particle velocities ckla are defined as ckla ¼ ck cla .
(26)
The local equilibrium distribution function (20) is applied to conditions (15)–(19). Considering that the odd rank tensors vanish and the even rank tensors are expressed by Eqs. (21)–(23), the following equations to
Fig. 1. Dodecahedron has 20 vectors pointing to the vertex form the center. Icosahedron has 12 vectors.
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
133
Table 1
pffiffi pffiffi 2 Unit vectors pointing to the vertex from the center of a dodecahedron and icosahedron, where j ¼ 1þ2 5, l ¼ p1ffiffi3, and t ¼ pffiffiffiffiffiffiffiffiffi pffiffi 5þ 5
l : direction
Unit vector (clx ; cly ; clz )
Dodecahedron
l l l l
lð1; 1; 1Þ lð0; j1 ; jÞ lðj; 0; j1 Þ lðj1 ; j; 0Þ
Icosahedron
l ¼ 1–4 l ¼ 5–8 l ¼ 9–12
¼ 1–8 ¼ 9–12 ¼ 13–16 ¼ 17–20
tð0; j; 1Þ tð1; 0; jÞ tðj; 1; 0Þ
Table 2 Typical sixth rank tensors of the unit vectors of a dodecahedron, icosahedron, and the combined group Sixth rank tensor
Dodeca. m ¼ 20
Icosa. m ¼ 12
Dodeca.+Icosa. m ¼ 32
m P
clx clx clx clx clx clx
2:9630 (2.8571)
1:6000 (1.7143)
4:5630 (4.5714)
cly cly clx clx clx clx
0:3529 (0.5714)
0:5789 (0.3429)
0:9318 (0.9143)
cly cly cly cly clx clx
0:6842 (0.5714)
0:2211 (0.3429)
0:9053 (0.9143)
clz clz cly cly clx clx
0:2963 (0.1905)
0:0000 (0.1143)
0:2963 (0.3048)
l¼1 m P l¼1 m P l¼1 m P l¼1
Values in parentheses are those for isotropic tensors calculated by Eq. (23).
determine coefficients F k are derived: XX F k ¼ 1, k
X
F k c2k ¼
1 e, 16
(28)
F k c4k ¼
5 2 e, 24
(29)
F k c6k ¼
35 3 e, 36
(30)
F k c8k ¼
35 4 e. 6
(31)
k
X k
X k
X k
(27)
l
Five speeds of particle velocities are necessary to satisfy the above equations. We assume a rest particle ðc0 ¼ 0:0Þ and four groups of the 32 moving particles, whose speeds are c1 ; c2 ; c3 ; and c4 . Eqs. (27)–(31) are readily solved to give the following: F 0 ¼ 1 32ðF 1 þ F 2 þ F 3 þ F 4 Þ,
(32)
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
134
F1 ¼
F2 ¼
F3 ¼
F4 ¼
1 c21 ðc21 c22 Þðc21 c23 Þðc21 c24 Þ 35 4 35 2 5 2 2 c22 c23 c24 2 2 3 2 2 2 2 2 e ðc2 þ c3 þ c4 Þe þ ðc2 c3 þ c3 c4 þ c4 c2 Þe e , 6 36 24 16
c22 ðc22
1 c24 Þðc22 c21 Þ
c23 Þðc22
35 4 35 2 5 c2 c2 c2 e ðc3 þ c24 þ c21 Þe3 þ ðc23 c24 þ c24 c21 þ c21 c23 Þe2 3 4 1 e . 6 36 24 16
c23 ðc23
ð33Þ
ð34Þ
1 c21 Þðc23 c22 Þ
c24 Þðc23
35 4 35 2 5 c2 c2 c2 e ðc4 þ c21 þ c22 Þe3 þ ðc24 c21 þ c21 c22 þ c22 c24 Þe2 4 1 2 e . 6 36 24 16 1 c24 ðc24 c21 Þðc24 c22 Þðc24 c23 Þ 35 4 35 2 5 2 2 c21 c22 c23 2 2 3 2 2 2 2 2 e ðc1 þ c2 þ c3 Þe þ ðc1 c2 þ c2 c3 þ c3 c1 Þe e . 6 36 24 16
ð35Þ
ð36Þ
As far as the simulation being stably conducted, the values c1 ; c2 ; c3 ; and c4 do not affect the accuracy itself. We can utilize this freedom to obtain a stable range as wide as possible. We investigated the stable range by numerical experiments whether a small density disturbance grows in a uniform flow of speed U and internal energy e, for various combinations of the particle speeds. The stable range depends in a complicated manner on the selection of the particle speeds, the relaxation parameter, and the grid size. The combinations around the following speeds gave the best stable range as far as in the investigation: ðc0 ; c1 ; c2 ; c3 ; c4 Þ ¼ ð0:0; 1:0; 2:0; 3:0; 4:0Þ.
(37)
An example of the stable range for ðc0 ; c1 ; c2 ; c3 ; c4 Þ ¼ ð0:0; 1:0; 2:0; 3:0; 4:0Þ, f ¼ 0:002, and the grid size ðDx; DtÞ ¼ ð0:1; 0:001Þ is shown in Fig. 2. A theoretical stability study for the true optimization is being planned. 4. Verification of the proposed model The proposed thermal model was evaluated by numerical simulations. The evolution Eq. (1) was solved by applying the explicit scheme to time and the second upwind scheme to space. The distribution function of next step f new kl at node ðI; J; KÞ is qf kl qf kl qf kl 1 þ c þ c ¼ f c (38) f new Dt ðf kl f ð0Þ klx kly klz kl kl kl ÞDt, f qx qy qz 8 3f kl;I 4f kl;I1 þ f kl;I2 > > < qf kl 2Dx ¼ 3f 4f > qx kl;I kl;Iþ1 þ f kl;Iþ2 > : 2Dx 8 3f kl;J 4f kl;J1 þ f kl;J2 > > > < 2Dy qf kl ¼ 3f 4f > qy kl;J kl;Jþ1 þ f kl;Jþ2 > > : 2Dy
if cklx X0; if cklx o0:
if ckly X0; if ckly o0:
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
135
5 Mach=3
Mach=2.5 Mach=2 UNSTABLE
4
Mach=1.5 Speed U
3
Mach=1
2 STABLE
Mach=0.5
1
0
1
2 3 Internal energy e
4
5
Fig. 2. Stable range obtained by numerical experiments whether a small disturbance grows in a uniform flow for various U and e, and r ¼ 1. The conditions are ðc0 ; c1 ; c2 ; c3 ; c4 Þ ¼ ð0:0; 1:0; 2:0; 3:0; 4:0Þ, f ¼ 0:002, and ðDx; DtÞ ¼ ð0:1; 0:001Þ.
8 3f 4f kl;K1 þ f kl;K2 > > kl;K qf kl < 2Dz ¼ 3f 4f > qz kl;Kþ1 þ f kl;Kþ2 > : kl;K 2Dz
if cklz X0; if cklz o0;
where the node suffix: I; I 1; ::: are indicated only on the focused coordinate. 4.1. Couette flow simulation Dissipative characteristics were verified by Couette flow. The upper wall, which is H apart from the lower wall and has internal energy e2 , starts to move at speed U. The lower wall has e1 and is at rest. Theoretical internal energy distribution e along vertical axis y in a steady state is given as y m y y 1 e ¼ e1 þ ðe2 e1 Þ þ 0 U 2 . (39) H 2k H H As the ratio m=2k0 is constant (¼ 0:3) from Eqs. (9) and (10), the distribution does not depend on the relaxation parameter f. Fig. 3 shows the simulation result of the proposed model for various relaxation parameters. The result agrees with the theoretical prediction. The same simulations were conducted using the models of Refs. [5–7]. The result for the Ref. [6] model is shown in Fig. 4. The result shows dependence on f, which contradicts the theoretical prediction. We think that this contradiction is because the model retains up to the third order of flow velocity in the equilibrium distribution function and ensures the isotropy only up to the fifth rank: the model does not satisfy constraints (18) and (19) in Section 3. Actually, the 3-D model proposed by us in Ref. [7], which was constructed to overcome the Ref. [6] model’s deficiency by increasing polynomial terms in the equilibrium distribution
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
136
Internal energy subtracted by linear distribution
0.0010
0.0008
0.0006
0.0004
φ =0.05 φ =0.1
0.0002
φ=0.2 φ =0.4 Theory
0.0
0.5 Vertical position y/ H
1.0
Fig. 3. Couette flow simulation using the model proposed in this paper. Internal energy at the steady state for U ¼ 0:1 and e1 ¼ e2 ¼ 1:0. The internal energy subtracted by linear distribution, the last term in Eq. (39), is shown. The relaxation parameter is changed: f ¼ 0:05; 0:1; 0:2; 0:4. The results for all f overlap each other. The number of grids is 161, ðDy; DtÞ ¼ ð0:125; 0:0125Þ.
function and adding several sets of particle speeds to satisfy Eqs. (18) and (19), gives a correct result for the Couette simulation. The 3-D model of Ref. [5] did not give a rational solution: the temperature in the fluid showed a lower value than that of the walls. It seems that the author made a mistake in tensor expressions when deriving the 3-D model (see Ref. [7]). 4.2. Normal shock simulation Compressibility was confirmed by normal shock wave. A uniform supersonic flow: M 1 , r1 , e1 , P1 ð¼ 23 r1 e1 Þ is given at the upstream boundary. The pressure at the downstream boundary P2 is gradually increased. When the pressure P2 becomes a certain value, a stationary flow discontinuity or normal shock wave is generated. The ratios of flow parameters between upstream and downstream of the shock for various upstream Mach numbers M 1 were compared with the theoretical value: Rankine Hugoniot relations in a perfect gas. The result, shown in Fig. 5, agrees with the theory. Regarding the shock structure and the shock thickness, analytical solutions exist only when the shock is weak [10]. In the weak shock theory, the pressure along the x axis, normal to the shock, is expressed as the following, where the position x0 is taken at P ¼ 12 ðP2 þ P1 Þ: x x 1 1 0 P ¼ ðP2 þ P1 Þ þ ðP2 P1 Þ tanh . 2 2 d
(40)
The shock thickness d for monatomic perfect gas with BGK assumption is d¼
8aV 4cm ¼ ðP2 P1 Þðq2 V =qP2 Þs ðP2 P1 Þðg þ 1Þ
2 2 þ ðg 1Þ , 3
where V is the specific volume (¼ 1=r) and a is the sonic absorption parameter.
(41)
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
137
Internal energy subtracted by linear distribution
0.0020
0.0015
0.0010
φ = 0.05
0.0005
φ = 0.1 φ=0.2 φ = 0.4 Theory 0
0.5 Vertical position y/ H
1
Fig. 4. Couette flow simulation using the 3-D model of Ref. [6]. Internal energy at the steady state for U ¼ 0:1 and e1 ¼ e2 ¼ 1:0. The internal energy subtracted by linear distribution is shown. The relaxation parameter is changed: f ¼ 0:05; 0:1; 0:2; 0:4. The result shows dependence on f, which contradicts the theoretical prediction. The grid conditions are same as Fig. 3.
The shock structure and the shock thickness obtained in the simulations for M 1 ¼ 1:02 are shown in Fig. 6 and in Fig. 7, respectively. Both results agree with the theory. 4.3. Supersonic nozzle simulation Finally, to demonstrate the ability to handle a generalized coordinate system, a Laval nozzle that accelerates the flow from Mach 0.2 to Mach 2.0 was simulated. The generalized coordinate system fitted with the nozzle shape is shown in Fig. 8. The generalized coordinate ði; jÞ has the following relations with the cartesian ðx; yÞ. The increment Dx is constant. The increment Dy is a function of x, or i. x ¼ Dx i, y ¼ Dy j, 2 1125 1 ð0:036i þ 0:2Þ2 Dy ¼ 1þ Dx. 5776 0:036i þ 0:2 3 The evolution equation expressed in the generalized coordinate system is qf kl cklx qf kl 1 qDy qf kl 1 ckly cklx j þ þ ¼ ðf kl f ð0Þ kl Þ. Dy qx f qt Dx qi qj
ð42Þ
(43)
The explicit and the second upwind scheme is also applied to this generalized coordinate formulation. Upwind/downwind is judged by the signs of cklx in the i direction and ckly cklx ½qDy=qxj in the j direction.
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
138
3.5 Theory Pressure
Ratio of flow parameters
3.0
P2 / P1
Density ρ2 / ρ1 Int. energy e2 / e1 M2
2.5
2.0
1.5
1.0
0.5 1.0
1.2
1.4
1.6
Upstream Mach M1 Fig. 5. Normal shock simulation. Ratios of flow parameters between upstream and downstream of the shock for f ¼ 0:005 versus various upstream Mach numbers. The number of grids is 101. ðDx; DtÞ ¼ ð0:01; 0:001Þ for M 1 ¼ 1:01 to 1.3, and ðDx; DtÞ ¼ ð0:005; 0:0002Þ for M 1 ¼ 1:5.
1.60 Pressure (weak shock theory) Pressure Density Int. energy Mach
1.05
1.55
1.00
1.50
Density
Pressure, Int. energy, Mach
1.10
0.95
0
1
2
3 Position x
4
5
6
1.45
Fig. 6. Normal shock simulation. Shock structure for M 1 ¼ 1:02 and f ¼ 0:01. The number of grids is 101. ðDx; DtÞ ¼ ð0:04; 0:002Þ.
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
139
1.0 Weak shock theory Simulation
Normal shock thickness
0.8
0.6
0.4
0.2
0.000
0.005 0.010 Relaxation parameter
0.015
Fig. 7. Normal shock simulation. Thickness of the shock for M 1 ¼ 1:02. The relaxation parameter is changed: f ¼ 0:0005; 0:001; 0:002; 0:005; 0:01. The number of grids is 101. The grid sizes is changed from ðDx; DtÞ ¼ ð0:01; 0:0005Þ for f ¼ 0:0005 to ðDx; DtÞ ¼ ð0:04; 0:002Þ for f ¼ 0:01.
y
1 j=10
0 j=0
i=2 -1 j=-10 i=1 i=0 0
i=50
1
2
x
3
4
5
Fig. 8. Supersonic nozzle simulation. Generalized coordinate system fitted with nozzle shape is applied. The grids are 51 21. Dx is constant ð¼ 0:1Þ, Dy is a function of x defined in Eq. (42), and Dt ¼ 0:002.
Analytical solution based on one-dimensional (1-D) isentropic variation predicts the relation of the cross section area A, the internal energy, the density, and the pressure with the Mach number: A 1 2 g 1 2 gþ1=2ðg1Þ 1þ M ¼ , (44) A M g þ 1 2 e ¼ e
g1 g1=g r P 2 g 1 2 1 1þ M ¼ ¼ , r P gþ1 2
(45)
where the subscript indicates the values at the throat. To ensure isentropic variation, sufficiently low viscosity: f ¼ 0:002 was used in the simulation. The upper and lower edges are slip and adiabatic boundaries. Initially, a uniform fluid of r1 ; e1 ; P1 ð¼ 23 r1 e1 Þ in the nozzle is at rest. The density and the internal energy, therefore the pressure at the entrance boundary: r2 ; e2 ; P2 ð¼ 23 r2 e2 Þ, are gradually increased, while keeping the pressure at the exit boundary the initial value
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
140
P1 by adjusting the internal energy. The velocities at both boundaries and the density at the exit are extapolated from the inner fluid region. When the pressure at the entrance boundary becomes sufficiently high, a supersonic flow is exerted in the expansion part of the nozzle. The simulation result of Mach contour lines in the nozzle is shown in Fig. 9. The Mach number, the internal energy, the density, and the pressure are averaged at each section. The ratios of these flow parameters to the values at the throat, versus the averaged Mach number, are shown in Fig. 10. Agreement with the theory of isentropic acceleration is confirmed, though a slight undevelopment of inflative acceleration is observed at the end of the nozzle.
y
1
Mach 2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2
0
-1 0
1
2
3
x
4
5
Fig. 9. Supersonic nozzle simulation. Mach contours: M ¼ 0:2; 0:3; . . . ; 2:0 in the nozzle. Initially uniform fluid of ðr1 ; e1 ; P1 Þ ¼ ð0:286; 0:651; 0:124Þ is at rest in the nozzle. The pressure at the entrance is increased to ðr2 ; e2 ; P2 Þ ¼ ð1:0; 1:5; 1:0Þ.
3.0 Theory
Ratio of flow parameters
2.5
2.0
Area
A / A*
Int. energy
e / e*
Density
ρ / ρ*
Pressure
P / P*
1.5
1.0
0.5
0.0 0.2
0.5
1.0 Mach
1.5
2.0
Fig. 10. Supersonic nozzle simulation. Ratios of the averaged flow parameters: A=A ; e=e ; r=r , and P=P to the averaged Mach number.
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
141
5. Conclusions A thermally correct 3-D multispeed FDLBM model was proposed. To derive thermally correct fluid equations, the local equilibrium distribution function should contain up to the fourth order terms of flow velocity, and tensors made of particle velocities should be isotropic up to the seventh rank. Vectors pointing to the vertex from the center of a dodecahedron and icosahedron have a good isotropy. If these 32 vectors are combined as a group, the above tensor isotropy requirement can be realized. Particle velocities of the model consist of a rest particle and four groups of these vectors. The model was verified by simulations of Couette flow, normal shock wave, and supersonic nozzle. The versatility of the model was demonstrated, from subsonic to supersonic flows, by showing good agreements with the theoretical predictions. A theoretical stability study for the true optimization of the particle speeds is being planned. Appendix A. Chapman Enskog expansion The Chapman Enskog expansion for the FDLBM formulation is presented briefly to show that the formulation is equivalent to Navier Stokes equations if the local equilibrium distribution function satisfies the velocity moments’ condition [2,11]. The FDLBM formulation is qf kl qf 1 þ ckla kl ¼ ðf kl f ð0Þ kl Þ, f qt qra XX f kl , r¼ k
rua ¼
XX
f kl ckla ,
(A.3)
l
XX c2 u2 f kl k . r eþ ¼ 2 2 k l The velocity moments’ condition is X X ð0Þ f kl ¼ r,
l
XX k
(A.6)
c2k u2 ¼r eþ , 2 2
f ð0Þ kl ckla cklb f ð0Þ kl
2 ¼ r edab þ ua ub , 3
c2k 5 u2 ckla ¼ rua e þ , 3 2 2
l
(A.7)
(A.8)
(A.9)
f ð0Þ kl ckla cklb cklg
2 ¼ r eðua dbg þ ub dga þ ug dab Þ þ ua ub ug , 3
(A.10)
c2k ckla cklb f ð0Þ kl
10 u2 7 u2 eþ ¼r e dab þ ua ub e þ . 9 3 3 2
(A.11)
l
XX k
f ð0Þ kl
l
XX k
f ð0Þ kl ckla ¼ rua ,
l
XX k
(A.5)
l
XX k
(A.4)
l
XX k
(A.2)
l
k
k
(A.1)
2
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
142
Actually, since Eqs. (A.7) and (A.9) can be derived by contraction from (A.8) and (A.10), respectively, independent requirements are Eqs. (A.5), (A.6), (A.8), (A.10), and (A.11). We assume that the relaxation parameter f is small and the distribution function f kl is a deviation from the local equilibrium f ð0Þ kl . The distribution function is expanded around the local equilibrium: ð1Þ 2 ð2Þ f kl ¼ f ð0Þ kl þ ff kl þ f f kl þ .
Eq. (A.1) becomes " # " # ð0Þ ð1Þ ð1Þ qf ð0Þ qf qf qf ð1Þ ð2Þ kl kl þ ckla kl þ f kl þ f þ ckla kl þ f kl þ ¼ 0. qt qra qt qra
(A.12)
(A.13)
Using the relations (A.2)–(A.7), the following velocity moments of the perturbation distribution f ðnÞ kl ðn ¼ 1; 2; . . .Þ can be assumed to be zero: X X ðnÞ f kl ¼ 0, (A.14) k
l
XX k
(A.15)
l
XX k
f ðnÞ kl ckla ¼ 0,
l
f ðnÞ kl
c2k ¼ 0. 2
(A.16)
For the zero-th approximation: qf ð0Þ qf ð0Þ kl þ ckla kl þ f ð1Þ (A.17) kl ¼ 0, qt qra P P P P c2 P P taking the velocity moments: k l 1, k l ckla , k l 2k , and considering Eqs. (A.14)–(A.16), we obtain the following equations: q X X ð0Þ q X X ð0Þ f kl þ f ckla ¼ 0, qt k l qra k l kl
(A.18)
q X X ð0Þ q X X ð0Þ f kl ckla þ f ckla cklb ¼ 0, qt k l qrb k l kl
(A.19)
q X X ð0Þ c2k q X X ð0Þ c2k þ ckla ¼ 0. f kl f qt k l 2 qra k l kl 2
(A.20)
Substituting the velocity moments’ condition (A.5)–(A.9) into Eqs. (A.18)–(A.20), we obtain the Euler (compressible and inviscid fluid) equations: qr q þ ðrua Þ ¼ 0, qt qra
(A.21)
q q 2 ðrua Þ þ redab þ rua ub ¼ 0, qt qrb 3
(A.22)
q u2 q u2 2 r eþ ¼ 0. rua e þ þ e þ qt qra 2 2 3
(A.23)
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
143
Next, the first order OðfÞ term in Eq. (A.13) is included. Taking the velocity moments: P P P P c2k P P k l 1; k l ckla ; k l 2 , we obtain the following equations: X X X X q q f ð0Þ þ f ð0Þ ckla ¼ 0, (A.24) qt k l kl qra k l kl q X X ð0Þ q X X ð0Þ q X X ð1Þ f kl ckla þ f kl ckla cklb þ f f ckla cklb ¼ 0, qt k l qrb k l qrb k l kl
(A.25)
q X X ð0Þ c2k q X X ð0Þ c2k q X X ð1Þ c2k þ ckla þ f ckla ¼ 0. f kl f kl f qt k l qra k l kl 2 2 qra k l 2
(A.26)
ð0Þ The f ð1Þ kl is replaced by the f kl expression, using zero-th approximation (A.17): q X X ð0Þ q X X ð0Þ f kl þ f ckla ¼ 0, qt k l qra k l kl
q X X ð0Þ q X X ð0Þ f kl ckla þ f ckla cklb qt k l qrb k l kl q q X X ð0Þ q X X ð0Þ f f kl ckla cklb þ f ckla cklb cklg qrb qt k l qrg k l kl q X X ð0Þ c2k q X X ð0Þ c2k þ ckla f kl f qt k l 2 qra k l kl 2 q q X X ð0Þ c2k q X X ð0Þ c2k ckla þ ckla cklb f f kl f qra qt k l qrb k l kl 2 2
(A.27)
! ¼ 0,
ðA:28Þ
! ¼ 0.
ðA:29Þ
Substituting Eqs. (A.5)–(A.11) into the above equations, we obtain equations expressed by macroscopic quantities. You may be perplexed regarding how to handle time derivatives arising from the terms P P P P ð0Þ 2 q=qt k l f ð0Þ l f kl ½ck =2ckla . These time derivatives are replaced by the following spatial kl ckla cklb and q=qt k derivatives derived from Euler equations: qr qua qr ¼ r ua , (A.30) qt qra qra qua 2 qe 2e qr qua ¼ ub , 3 qra 3r qra qt qrb
(A.31)
qe qe 2e qua ¼ ua , qt qra 3 qra
(A.32)
then, we obtain the Navier Stokes equations: qr q þ ðrua Þ ¼ 0, qt qra q q 2 q 2 qub qua 2 qug ðrua Þ þ redab þ rua ub ref þ dab ¼ 0, qt qrb 3 qrb 3 qra qrb 3 qrg q u2 q u2 2 r eþ rua e þ þ e þ qt qra 2 2 3 q 10 qe 2 qub qua 2 qug ref þ refub þ dab ¼ 0. qra 9 qra 3 qra qrb 3 qrg
(A.33) (A.34)
ðA:35Þ
ARTICLE IN PRESS M. Watari, M. Tsutahara / Physica A 364 (2006) 129–144
144
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11]
S. Chen, G. Doolen, Annu. Rev. Fluid Mech. 30 (1998) 329. G. McNamara, B. Alder, Physica A 194 (1993) 218. F.J. Alexander, S. Chen, J.D. Sterling, Phys. Rev. E 47 (1993) R2249. N. Cao, S. Chen, S. Jin, D. Martinez, Phys. Rev. E 55 (1997) R21. Y. Chen, H. Ohashi, M. Akiyama, Phys. Rev. E 50 (1994) 2776. N. Takada, Y. Yamakoshi, M. Tsutahara, Numerical analysis of fluid motions by three-dimensional thermal lattice Boltzmann model, JSME B, vol. 64, 1998, pp. 4 (in Japanese). M. Watari, M. Tsutahara, Phys. Rev. E 70 (2004) 16703. M. Watari, M. Tsutahara, Phys. Rev. E 67 (2003) 36306. S. Wolfram, J. Statist. Phys. 45 (1986) 471. L.D. Landau, E.M. Lifshitz, Fluid mechanics, Oxford, 1987. S. Chapman, T.G. Cowling, The Mathematical Theory of Non-uniform Gases, Cambridge University Press, London, 1970.