- Email: [email protected]
Optimal relaxation collisions for lattice Boltzmann methods
Computers and Mathematics with Applications 65 (2013) 172–185
Contents lists available at SciVerse ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
Optimal relaxation collisions for lattice Boltzmann methods Fuzhang Zhao ∗ APD Optima 266 Robinson Drive, Tustin, CA 92782, USA
article
info
Keywords: Functioned polylogarithm polynomial Lattice Boltzmann equations Optimal relaxation time Truncation error
abstract The optimal relaxation time of about 0.8090 has been proposed to balance the efficiency, stability, and accuracy at a given lattice size of numerical simulations with lattice Boltzmann methods. The optimal lattice size for a desired Reynolds number can be refined by reducing the Mach number for incompressible flows. The functioned polylogarithm polynomials are defined and used to express the lattice Boltzmann equations at different time scales and to analyze the impact of relaxation times and lattice sizes on truncation errors. Smaller truncation errors can be achieved when relaxation times are greater than 0.5 and less than 1.0. The steady-state lid-driven cavity flow was chosen to validate the code of lattice Boltzmann procedures. The applications of the optimal relaxation parameters numerically balance the stability, efficiency, and accuracy through Hartmann flow. The optimal relaxation time can also be used to select the initial lattice size for the channel flow over a square cylinder with a given Mach number. © 2011 Elsevier Ltd. All rights reserved.
1. Introduction The lattice Boltzmann method (LBM) has emerged as a promising alternative to the numerical scheme for high performance computational fluid dynamics (CFD) in the past twenty years. The numerical approaches based on the Boltzmann equation cover a wide range of the Knudsen number for practical applications. No effort for solving linear systems of equations is needed. The method is particularly competitive for fluid flow applications with complexities of geometries, boundaries, and interfacial dynamics [1]. The most commonly used collision operators are the single relaxation time (SRT) [2] and the multiple relaxation time (MRT) [3,4]. The LBMs can be applied to physical problems by means of nondimensionalization [5]. When both the physical system and the lattice system are equivalently scaled, three key dimensionless quantities: lattice size, relaxation time or collision frequency, and Mach number are yet to be determined for any CFD problems with a given Reynolds number. Different combinations of the key dimensionless quantities may numerically impact the stability, accuracy, and efficiency of the LBMs. The optimal relaxation time has been derived to balance the numerical properties. The optimal relaxation parameters are applicable for both the SRT and the MRT collision operators since no particular assumptions relative to collision operators have been made during their derivations. The lattice Boltzmann method, which can be derived from the kinetic theory of gases, evolved from the lattice-gas automata. The primary variable is a one-particle probability distribution function f (x, e, t ). The macroscopic flow variables such as density, velocity, pressure, temperature, and stress can be expressed as the velocity moments of the one-particle probability distribution function. The general transport equation reads
(∂t + e · ∇ x + a · ∇ e )f (x, e, t ) = Q , ∗
Tel.: +1 714 389 0228. E-mail address: [email protected].
0898-1221/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.camwa.2011.06.005
(1)
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
173
where t is the time, e is the particle velocity, x is the special location, a is the acceleration due to an external force, and Q is a collision operator. To complete the Boltzmann equation from (1), the collision operator has to be specified. Four assumptions are made: (a) only binary collisions are taken into consideration, (b) the velocity of a molecule is not correlated to its space location, (c) the effect of the external force on the collision cross section is neglected, and (d) wall effects are ignored [6]. Under these assumptions, Boltzmann related the collision operator to the one-particle probability distribution function. One of the most commonly used collision operators is the linearized SRT or BGK collision operator [2] 1 Q = − [f − f (eq) ],
(2)
τ
where τ is a relaxation time. With the Boltzmann H-theorem, the Maxwell–Boltzmann or equilibrium one-particle probability distribution function was given by f
(eq)
=
ρ
exp −
(2π RT )D/2
(e − u)2
2RT
,
(3)
where D, R, T , ρ , and u are the dimensions of space, gas constant, temperature, density, and velocities, respectively. To numerically solve for f given in (1)–(3), the velocity space e is simplified by a finite set of velocity vectors ei , which is chosen for isotropy. For the D2Q9 model, the discrete velocity vectors can be practically expressed as e0 = (0, 0), e1 = −e3 = (c , 0), e2 = −e4 = (0, c ), e5 = −e7 = (c , c ), and e6 = −e8 = (−c , c ). The lattice speed c = δ x/δ t is the ratio of a lattice spacing δ x to a time step δ t. The low Mach number expansion of (3) at a constant temperature is given by [7] fi
(eq)
9 3 3 = ρwi 1 + 2 (ei · u) + 4 (ei · u)2 − 2 (u · u) , c
2c
2c
(4)
where the weighting factors wi are specifically given by w0 = 4/9, w1 = w2 = w3 = w4 = 1/9, and w5 = w6 = w7 = w8 = 1/36. The leading moments of the truncated equilibrium distribution function (4) for incompressible flows can be evaluated as
ρ=
8
fi =
8
i=0
ρ uα =
8
fi
3 1 3
,
(5)
i=0
eiα fi =
i=1
1
(eq)
8
eiα fi
(eq)
,
(6)
i=1
c 2 ρδαβ + ρ uα uβ =
8
eiα eiβ fi
(eq)
,
(7)
i=1
c 2 ρ(δαβ uγ + δβγ uα + δγ α uβ ) =
8
eiα eiβ eiγ fi
(eq)
,
(8)
i=1
where eiα is the projection of ei on α -axis. In (6) through (8), the summation is taken from 1 to 8 because of a rest particle devised for better computational stability. For ideal gases, the equation of state relates density ρ to pressure p with the following equation p = ρ RT =
1 3
ρc2.
(9)
Eqs. (1) and (2), by introducing the collision frequency ω = 1/τ , can be discretized in space x and time t with the lattice spacing of δ x = ei δ t as fi (x + ei δ t , t + δ t ) − fi (x, t ) = −ω[fi (x, t ) − fi
(eq)
(x, t )].
(10)
This is the lattice Boltzmann equation with the BGK simplification. The lattice BGK model can be solved in two steps: first by collision and then by streaming. Therefore, the relaxation of local collision and advection of particle distribution functions evolve from one lattice to another by (eq) fi (x, t + δ t ) = (1 − ω)fi (x, t ) + ωfi (x, t ),
(11)
fi (x + ei δ t , t + δ t ) = fi (x, t + δ t ),
(12)
174
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
where fi denotes the post-collision state of the particle distribution function. Since the collision step is local and the streaming step is uniform, they can be implemented straightforwardly. The kinematic viscosity in the Navier–Stokes equations for the D2Q9, D3Q15, and D3Q19 models is
ν=
2τ − 1
c 2 δt . 6 It is worthwhile to notice that the kinematic viscosity for other models might be different.
(13)
2. Optimal relaxation parameters According to Buckingham’s rule, there are usually fewer dimensionless groups than there are physical quantities. The groups make up the real characteristic and meaningful variables of the problem. The two physical quantities of density and dynamic viscosity appear in the standard incompressible Navier–Stokes equations without the gravitation term and by introducing several nondimensional variables, the Reynolds number has been shown to be the only dimensionless group in the corresponding nondimensional form of the equations. Since all flows of the same Reynolds number are mathematically comparable, the Reynolds number will be used to derive the optimal relaxation time and the optimal lattice size. The Reynolds number is usually written as Re =
U0 L
ν
,
(14)
where U0 is the characteristic fluid velocity and L is the characteristic length that the flow is going through or around. Since the equilibrium one-particle probability distribution function (4) is isothermal, the Mach number can be given by
√ Ma =
U0
3U0
=
cs
c
,
(15)
√ where cs is a sound speed, and under an isothermal condition, we have a constant sound speed cs = c / 3. The collision frequency ω determines the fraction of particles that collide in a small volume during a time interval δ t, based on the collision interval theory [6]. The lattice size Nx implies the total number of particles in the x-coordinate. Choosing an appropriate collision frequency or a relaxation time, a lattice size, and a Mach number, and scaling nondimensional quantities of interest are crucial for the lattice Boltzmann methods to successfully apply to practical physical problems with the desired Reynolds number. As we have known, the relative error of the velocity is proportional to Ma2 , and increasing the collision frequency ω will decrease the stability and the efficiency. Thus, it is necessary to establish the relationship between the collision frequency ω and Mach number Ma. The balance among the stability, efficiency, and accuracy at a given lattice size can then be obtained. Substituting (13) into (14), using (15), and rearranging yields
ω 2
1
=
√
1+
Ma 2 3L Re δ x
.
(16)
√
Since the first dimensionless group (ω/2) is inversely related to the second dimensionless group (Ma/Re2 3L/δ x) in (16), increasing one dimensionless group will decrease the other dimensionless group. For simulations of incompressible flows with a given Re, we usually use a smaller Ma to reduce the compressibility and to enhance the accuracy at a given lattice size (L/δ x). Reducing Ma, however, will increase ω and as a result, the efficiency will be lower. As the collision frequency ω approaches 2.0, the stability will deteriorate. Considering both the stability and the efficiency, we wish to obtain a smaller collision frequency. As we have mentioned, based on (16), a smaller ω requires a greater Ma. Thus, there must be a balance among the efficiency, stability, and accuracy for LBM simulations at a given lattice size. For the balance, the two dimensionless groups in (16) must be equalized
ωo
√ =
Ma 2 3L
.
2 Re δ x Substituting (17) into (16), replacing ω with ωo , and rearranging gives
ωo2 + 2ωo − 4 = 0.
(17)
(18)
The optimal collision frequency, ωo , or the reciprocal of the optimal relaxation time, τo , can be obtained from the positive root of the quadratic equation (18)
ωo =
1
τo
√
5 − 1.
=
(19)
The accompanied equation useful for selecting the optimal lattice size, Nox , can also be obtained by rearranging (17) and using (19)
√ Nox = 1 +
L
δx
=1+
√
15 −
3 Re
12
Ma
,
(20)
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
175
and the representative optimal mean free pass is
√ λo =
√
6π ( 5 − 1)
δx 24 if the Knudsen number is defined as Kn =
λ L
=
Ma
Re
π 2
(21)
.
(22)
Based on the derivation given above, the optimal relaxation time applies to incompressible flows with all discretization models since the optimal relaxation time is independent of discretization model. The optimal lattice size, however, is dependent on discretization model because the kinematic viscosity may vary with the discretization model. For the kinematic viscosity with the form other than (13), a new equation for the optimal lattice size needs to be derived by following the procedures from the derivation of (16) to that of (20). For incompressible CFD applications with a desired Re, the Mach number must be controlled to reduce compressibility effects. Reducing Ma will increase accuracies for the observed quantities such as the velocity and the density since the relative errors for both the velocity and the density are proportional to Ma2 . For the desired Reynolds number, the selection of an optimal lattice size can be coordinated with the Mach number using (20). The ratio of the Mach number to the Reynolds number can also be expressed by the Knudsen number Kn in (22) for isothermal ideal gases. The optimal relaxation parameters can be applied to both the SRT and the MRT collision operators. In the case of the MRT collision operator, we may have the relationship of s8 = s9 = ωo = 1/τo for the D2Q9 model. The bulk viscosity has a relationship with s2 . Therefore, the bulk and shear viscosities can be independently controlled in the MRT collision operator. 3. Analytical analysis In order to study truncation errors, lattice Boltzmann equations (LBE) for different time scales will be derived based on the Chapman–Enskog expansion. A new polynomial which relates to the polylogarithm function will be defined and used. As a result, the incompressible Navier–Stokes equations along with truncation errors will be recovered. 3.1. Functioned polylogarithm polynomial Slightly different from the definition in [8], the functioned polylogarithm polynomial in τ of order m is defined as
Lm (τ ) ≡
1 (−1)m Li−m 1 − , τ −1 τ
(23)
where Li−m is the polylogarithm of a negative integer order and m is a non-negative integer. The functioned polylogarithm polynomial has many useful properties. The two essential properties have been derived. The first property is the following recursion formula d[Lm (τ )]
− τ Lm (τ ), with L0 (τ ) = 1. (24) dτ The second property of the functioned polylogarithm polynomial of consecutive order has the following inequality relationship Lm+1 (τ ) = τ (1 − τ )
|Lm (τ )| 6 |Lm+1 (τ )| m = 1, 2, . . . , 1 6 τ < ∞.
(25)
The first six orders of the polynomials have been obtained from (24)
L1 (τ ) = −τ ,
(26)
L2 (τ ) = τ (2τ − 1),
(27)
L3 (τ ) = −τ (6τ − 6τ + 1), 2
(28)
L4 (τ ) = τ (2τ − 1)(12τ − 12τ + 1), 2
(29)
L5 (τ ) = −τ (120τ − 240τ + 150τ − 30τ + 1),
(30)
L6 (τ ) = τ (2τ − 1)(360τ − 720τ + 420τ − 60τ + 1).
(31)
4
3
2
4
3
2
Their applicable maximum roots are exactly obtained as
√ τ3 =
1 2
+
3 6
√ ,
τ4 =
1 2
+
6 6
√
,
τ5 =
1 2
+
1 8
+
√
105
120
,
τ6 =
1 2
+
1 6
+
15 60
,
where the subscripts of τ , such as 3, 4, 5, 6, denote the order of the polynomials obtained from (24).
(32)
176
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
Fig. 1. Maximum and averaged applicable roots.
All the non-zero roots are positive and the maximum roots satisfy the following relationship 0.5 < τ3 < τ4 < τ5 < τ6 < · · · < τm < 1.0.
(33)
Since positive kinematic viscosities require 0.5 < τ < ∞, the values of τ1 = 0 and τ2 = 0.5 cannot be used in lattice Boltzmann calculations. The maximum roots and the averaged applicable roots for the first fifteen orders are plotted in Fig. 1. The applicable roots refer to all the roots of (23) that are greater than 0.5. Some details on the explicit higher order functioned polylogarithm polynomials are given in Appendix. The impact of using different relaxation times on the truncation errors will be demonstrated by the following Chapman–Enskog expansion analysis. 3.2. Lattice Boltzmann equations in different time scales To write the lattice Boltzmann equations in different time scales ϵ = δ t, the Chapman–Enskog expansion is usually applied [9]. With the following equations fi =
∞
ϵ n fi(n) ,
(34)
n =0
∞ ∂ ∂ = ϵn , ∂t ∂ tn n=0
(35)
∞ ϵn ∂ ∂ n fi (x + ei δ t , t + δ t ) = + ei · fi (x, t ), n! ∂ t ∂x n =0
(36)
and
we can expand the LBE (10) into a series of LBEs in different orders of the time scale ϵ . They are, by letting D0 = ∂t0 + ei · ∂x , expressed as fi fi fi fi fi
(0)
= fi(eq) ,
(1)
= L1 (τ )D0 fi(0) ,
(37) (38)
(0)
(2)
= L1 (τ )
∂ fi L2 (τ ) 2 (0) + D0 fi , ∂ t1 2!
(39)
(3)
= L1 (τ )
∂ fi(0) ∂ f (0) L3 (τ ) 3 (0) + L2 (τ )D0 i + D 0 fi , ∂ t2 ∂ t1 3!
(40)
(4)
= L1 (τ )
(0) ∂ fi(0) ∂ f (0) L3 (τ ) 2 ∂ fi L4 (τ ) 4 (0) + L2 (τ )D0 i + D0 + D0 fi + A4 , ∂ t3 ∂ t2 2! ∂ t1 4!
(41)
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
177
and for the equation of order ϵ m , one may generally have fi
(m)
= L1 (τ ) +
(0) ∂ fi(0) ∂ f (0) L3 (τ ) 2 ∂ fi + L2 (τ )D0 i + + ··· D0 ∂ tm−1 ∂ tm−2 2! ∂ tm−3
Lm−1 (τ )
(m − 2)!
−2 Dm 0
∂ fi(0) Lm (τ ) m (0) + D0 fi + Am , ∂ t1 m!
(42)
where A4 and Am are the terms including other higher order and cross derivatives for the orders ϵ 4 and ϵ m equations, respectively. 3.3. Truncation errors Since the Burnett and superBurnett equations are obtained from the third and fourth order truncations of the Chapman–Enskog expansion of the Boltzmann equation, the Boltzmann equation in general is not necessarily converged to the Navier–Stokes equations [10]. It is useful, therefore, to study the impact of different combinations of relaxation times and lattice sizes on truncation errors. Similar to the procedure in [11], the incompressible Navier–Stokes equations with truncation errors may be recovered by neglecting the divergent term and convection terms of O(Ma2 ) or higher order for momentum equations
∂ρ + ∇ · (ρ u) = RC , ∂t ∂(ρ u) + ∇ · (ρ uu) = −cs2 ∇ρ + ν∇ 2 (ρ u) + RαM , ∂t
(43) (44)
where the truncation errors RC and Rα M for continuity and momentum equations boil down to
RC = ϵ 2
+ϵ Rα M = ϵ
L2
C2,1 +
L1
m−1
2
L2 L1
L2 L1
L3 L1
C2,2 + ϵ 3
Cm−1,1 +
Mα 2,1 +
+ · · · + ϵ m−1
L3
L1 L2
L1
L3 L1
L2 L1
C3,1 +
L3
L1 Lm
Cm−1,2 + · · · +
Mα 2,2 + ϵ
3
L2 L1
C3,2 +
L1
Lm L1
C3,3 + · · ·
L1
Cm−1,m−1 ,
Mα 3,1 +
Mα m−1,1 + · · · +
L4
L3 L1
Mα 3,2 +
(45)
L4 L1
Mα 3,3
Mα m−1,m−1 ,
(46)
where Cj,k and Mα j,k for j ∈ [2, 3, . . . , m − 1] and k ∈ [1, 2, . . . , j], independent of a relaxation time, are the higher order derivatives of density, ρ , and momentum, ρ uα , respectively, with respect to different time scales. For the continuity truncation errors, we have C2,1 = −
∂[D0 (ρ)] , ∂ t1
C3,1 = −
∂[D0 (ρ)] ∂ 2ρ − , ∂ t2 2!∂ t12
C2,2 = −
D30 (ρ) 3!
,
C3,2 = −
(47)
∂[D20 (ρ)] , 2!∂ t1
C3,3 = −
D40 (ρ) 4!
,
(48)
and Cm−1,m−1 = −
Dm 0 (ρ) m!
.
(49)
For the momentum truncation errors, we get Mα 2,1 = −
∂[D0 (ρ uα )] , ∂ t1
Mα 3,1 = −
∂[D0 (ρ uα )] ∂ 2 (ρ uα ) − , ∂ t2 2!∂ t12
M α 2 ,2 = −
D30 (ρ uα ) 3!
,
Mα 3,2 = −
(50)
∂[D20 (ρ uα )] , 2!∂ t1
Mα 3,3 = −
D40 (ρ uα ) 4!
,
(51)
and Mα m−1,m−1 = −
Dm 0 (ρ uα ) m!
.
(52)
178
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
In order to reduce the high order truncation errors, we can use the averaged applicable roots at different orders shown in Fig. 1. There are many ways to use the roots of the functioned polylogarithm polynomial to construct a relaxation time. A reasonable guess is the average-weighted relaxation time, which takes the following form
τa =
τ¯3 + ϵ τ¯4 + ϵ 2 τ¯5 + · · · + ϵ m−3 τ¯m , 1 + ϵ + ϵ 2 + · · · + ϵ m−3
(53)
where τ¯m denotes the averaged applicable root in order m. The averaged roots at any order are in the range of τ3 ≈ 0.7887 and τ4 ≈ 0.9082. As the polynomial order m increases, the effects of the averaged applicable roots fade away. Thus, the average τa should be closer to the lower limit τ3 and the optimal relaxation time τo . In other words, the optimal relaxation time reduces the higher order truncation errors of the Navier–Stokes equations and therefore increases the accuracy of simulations. If higher order truncation errors were the only concern, using τ3 as a relaxation time and a refined mesh would be a good choice to reduce truncation errors, making the Boltzmann equation converge to the Navier–Stokes equations. The absolute values of the lowest order truncation errors for continuity and momentum equations take the following forms, respectively
L2 (τ ) ∂[D0 (ρ)] , |RC 2,1 | = ϵ L1 (τ ) ∂ t1 L2 (τ ) ∂[D0 (ρ uα )] . |RαM2,1 | = ϵ 2 L1 (τ ) ∂ t1 2
(54)
(55)
For incompressible flows, we can select a lattice size based on Ma and Re. Thus, the lattice spacing, ϵ , is independent of a relaxation time. The truncation errors for the continuity equation and the momentum equations share the same form in terms of ϵ and Lm (τ ). The functioned polylogarithm polynomials have the following relationship
L2 (0.5 < τ < 1) L2 (1 6 τ < ∞) < L (0.5 < τ < 1) L (1 6 τ < ∞) . 1 1
(56)
With (54), (55), and (56), we obtain
|RC 2,1 (0.5 < τ < 1, ϵ)| < |RC 2,1 (1 6 τ < ∞, ϵ)|,
(57)
|RαM2,1 (0.5 < τ < 1, ϵ)| < |RαM2,1 (1 6 τ < ∞, ϵ)|.
(58)
The inequalities (57) and (58) indicate that the truncation errors using relaxation times within the range of 0.5 < τ < 1.0 are smaller than those with other applicable relaxation times at a given lattice size. 4. Numerical analysis with applications 4.1. Steady-state lid-driven cavity flow In order to demonstrate the benefits of the optimal relaxation time and the corresponding optimal lattice size, the algorithm has been implemented with MATLAB. The two-dimensional lid-driven cavity flow has been widely adopted for a CFD code validation because plenty of benchmark data is available [12]. The computational domain is a square. The top fluid moving with the constant horizontal velocity of U0 drives the flow inside the cavity. All the other velocity boundary conditions are set to zero. In the LBM simulations, the MLS scheme [13,14] was adopted to treat all the velocity boundary conditions. For LB simulations, the criterion of a steady state is satisfied by the following Euclidean norms Eu =
∥u(x, t + δ t ) − u(x, t )∥ 6 Ea , ∥u(x, t + δ t )∥
(59)
where Ea is the allowable tolerance for errors. In the lid-driven cavity flow simulations, Ea = 10−7 is used. The isothermal LBM is a nearly incompressible scheme. When we use the scheme to model incompressible flows, the compressibility effect on a solution must be carefully monitored and controlled. The compressibility can be described by the coefficient of variation of a density
φ=
1
ρ¯
(ρ − ρ) ¯ 2 /N ,
(60)
where the average density is given by the following equation
ρ¯ =
N
ρ
.
(61)
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
179
Table 1 Input and output for two-dimensional lid-driven cavity flow.
ω
τ
U0
φ
Nt
1.800 1.600 1.400 1.236 1.000 0.800 0.600
0.5556 0.6250 0.7143 0.8091 1.0000 1.2500 1.6667
0.0144676 0.0325521 0.0558036 0.0804848 0.1302083 0.1953125 0.3038194
7.6126 × 10−5 3.1069 × 10−4 8.4017 × 10−4 1.6761 × 10−3 4.1936 × 10−3 9.1303 × 10−3 2.1379 × 10−2
– – 103,172 80,783 59,975 46,096 33,763
Fig. 2. Velocity comparisons of two-dimensional cavity flow with Re = 100.
The lid-driven cavity flow with Re = 100, N = 129 × 129 lattices, and different relaxation times were chosen for the simulations by the LBM. In order to keep the Reynolds number a constant, different velocities U0 corresponding to the selected relaxation times have been obtained. The computational efficiency is measured by the total number of time steps Nt when the steady-state condition (59) is satisfied. The input and output for the benchmark case are listed in Table 1. Generally, as the relaxation time increases, the compressibility increases but the accuracy decreases while the efficiency increases. Specifically, on the one hand, when ω = 1.8 or 1.6 is used, the LBM does not converge to the desired tolerance of Ea = 10−7 despite a relatively low compressibility of φ = 7.6126 × 10−5 and φ = 3.1069 × 10−4 , respectively. On the other hand, when ω = 0.6 is used, the result possesses a relatively high compressibility of φ = 2.1379 × 10−2 , even though it is as efficient as Nt = 33763. The results in Table 1 can also be explained in the collision interval theory. When the fraction of particles for collisions is overdetermined (ω > ωo ), it will take more time to converge and the change from f to f (0) is small and linear in a time interval. When it is underdetermined (ω < ωo ), despite less time for convergence, the change from f to f (0) is large and nonlinear in a time interval, causing an increase in compressibility. The optimal collision frequency ωo determines the balanced and representative fraction of particles in collisions during a time interval and the optimal lattice size Nox given in (20) helps us to select the total number of particles within a given domain. The velocity results from the current simulation along the vertical and horizontal lines through the geometric center for Re = 100 agree with the results obtained from the finite difference method with a multigrid solver given in [12] and they are depicted in Fig. 2. The lid-driven cavity flow, however, has no analytical solution. The accuracy can be estimated based on the Richardson extrapolation. The results for the lattice 513 × 513 has been extrapolated from those of lattices 257 × 257 and 129 × 129, and are used as an ‘‘exact’’ solution for the studied cavity flow by ue ≈ u(513×513) = [4u(257×257) − u(129×129) ]/3. The relative error of velocity can then be calculated by |u − ue |/ |ue | for the cavity flow. For accuracy estimation, three different lattices, 257 × 257, 129 × 129, and 65 × 65, with Re = 100 and ωo ≈ 1.236, were used in LBM simulations of the lid-driven cavity flow. The relative change of velocity defined in (59) as a function of the total time steps has been plotted in Fig. 3. As time continues, the percentage error drops steeply from 1.00 to 0.01 for all the three different meshes, changes slope, and converges to Ea = 10−7 at the total time steps of 26,493, 80,783, and 272,333 from meshes of 65 × 65, 129 × 129, and 257 × 257, respectively. The details depicted in the inset graph of Fig. 3 demonstrate an oscillating convergence.
180
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
Fig. 3. Iteration behaviors of two-dimensional cavity flow with Re = 100.
The horizontal velocities along the vertical line through the geometric center of the cavity were also selected for accuracy estimations. The Errors calculated are 0.00393, 0.00109, and 0.00027 for meshes 65 × 65 (Ma = 0.2788), 129 × 129 (Ma = 0.1394), and 257 × 257 (Ma = 0.0697), respectively. Thus, the relative velocity error is proven to be proportional to Ma2 . Also monitored and controlled, the corresponding compressibility effects are 6.12155 × 10−3 , 1.67611 × 10−3 , and 4.50300 × 10−4 . 4.2. Hartmann flow The magnetohydrodynamic (MHD) effects have many engineering applications such as the electromagnetic casting of metals, the confinement of plasmas, magnetic flowmeters, and so forth. Hartmann flow refers to a steady-state incompressible flow of an electrically conducting fluid through a pair of parallel plates under a transversely prescribed magnetic field. The force driving the fluid along the parallel plates is the pressure gradient. The magnetic fields induce currents in the moving conductive fluid, which create forces on the fluid, and also change the magnetic field itself. With the boundary conditions of Ux (±H ) = 0 and Bx (±H ) = 0 on the plates, the analytical solution of Hartmann flow takes the following form [15] Ux (y) Un Bx (y) Bn
coth(Ha)
=
1−
Ha 1
=
sinh(Ha y¯ ) sinh(Ha)
Ha
cosh(Ha y¯ ) cosh(Ha)
,
(62)
− y¯ ,
(63)
√
where Un = GH 2 /(ρν), Bn = GH 2 / ρνη, G = −dp/dx, y¯ = y/H, and the Hartmann number is given by B0 H Ha = √ ,
(64)
ρνη
where B0 is the transverse magnetic field, and η is the magnetic diffusivity. The detailed procedures of using the LBM to simulate MHD problems can be found in [15]. The equilibrium distribution functions for both the hydrodynamics and the magnetic field are adopted respectively as fi
′ (eq)
(eq)
gi
= fi
(eq)
+ wi
9
2c 4
1 2
|B| |ei | − (ei · B) 2
2
= wi [B + 3ei · (uB − Bu)],
2
,
(65) (66)
where the magnetic field vector B is given by B=
8 i=0
gi =
8 i =0
(eq)
gi
,
(67)
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
181
Fig. 4. Error as a function of relaxation time domain for Hartmann flow.
and the gi follows the vector Boltzmann equation (eq)
gi (x + ei δ t , t + δ t ) − gi (x, t ) = −ωm [gi (x, t ) − gi
(x, t )].
(68)
The magnetic diffusivity through the Chapman–Enskog expansion analysis of (68) can be recovered for the D2Q9 model as
η=
2τm − 1 6
c 2 δt ,
(69)
where τm is a magnetic relaxation time. By the same token, the magnetic optimal collision frequency or the inverse of the magnetic optimal relaxation time has been obtained
ωmo =
√
1
τmo
=
5 − 1.
(70)
For the constraint of ∇ · B = 0, the magnetic isothermal Mach number may be defined based on the Alfvén velocity as
√ 3 Mam = B0 √ .
(71)
ρ
Therefore, the optimal lattice size, Nmox , can be refined by reducing the magnetic Mach number (71) based on the following equation
√ Nmox = 1 +
5−1 24
2
w
Ha −1 . Mam
(72)
Under the conditions of ω = ωo = ωmo , the magnetic optimal lattice size (72) is equivalent to the hydrodynamic one (20). The lattice size of Nx × Ny = 9 × 101, the half distance between plates of H = 1, and the pressure gradient of G = 2.0833 × 10−4 were chosen for the Hartmann flow simulations. The initial conditions of both velocity, ux , and magnetic field, bx , were assumed to be zero. The boundary conditions were prescribed using the nonequilibrium extrapolation scheme [16]. In order to reduce the impact of boundary conditions on the accuracy of simulations, the second order extrapolation scheme was used for the nonequilibrium distribution functions. The pressure gradient was implemented by applying both inlet pressure and outlet pressure boundary conditions. The allowable tolerance to reach a steady state is also set to be 10−7 for both the velocity field and the magnetic field based on the Euclidean norm (59). The overall relative errors between analytical and numerical solutions for both velocity and magnetic field can be evaluated by
E=
(Ux − ux )2 + (Bx − bx )2 , (Ux )2 + (Bx )2
(73)
where ux and bx are the streamwise velocity and induction magnetic field calculated by LBM, respectively. The LBM simulations of Hartmann flow with Ha = 50 were initially conducted within the relaxation time domain of τ , τm ∈ [0.65, 1.15]. The relative error defined in (73) as a function of the relaxation time domain is plotted in Fig. 4. Better results have been achieved when equalized relaxation times, τ = τm , are used. The accuracy as a function of equalized relaxation times has then been plotted in Fig. 5. When the optimal relaxation times of τ = τm = 0.8090, rather than τ3 or
182
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
Fig. 5. Error as a function of equalized relaxation times for Hartmann flow.
τ4 , are used, the error of E = 0.005412 is the least among the selected relaxation times for the simulations. For the extreme cases of τ = τm = 0.55 and 2.00, the simulations failed to converge. For the Hartmann flow with Ha = 100, a similar study has been conducted. Although the minimum error of E = 0.0088039 was calculated at τ = τm = τ3 and the slightly greater error of E = 0.0088396 was obtained at τ = τm = τo , the efficiency of Nt = 1212 and the compressibility of φ = 2.2030 × 10−5 at τo is better than the efficiency of Nt = 1454 and the compressibility of φ = 2.2032 × 10−5 at τ3 , correspondingly. As we have known, the induced magnetic field strength tends to be zero as Ha → 0 or Ha → ∞, meaning that the extreme values of Bx must exist at a Hartmann number of Ha ∈ (0, ∞). The minimum value as a function of Ha can be obtained Bx,min Bn
=
1
Ha
sinh(Ha y¯ s ) sinh(Ha)
− y¯ s ,
(74)
where the stationary point as a function of Ha is given below
y¯ s =
1 Ha
ln
sinh(Ha)
Ha
+
sinh(Ha) Ha
2
−1 .
(75)
Based on (74) and (75), the minimum induced magnetic field has the value of Bx,min /Bn ≈ −0.102818 located at y¯ s ≈ 0.633185 for Ha ≈ 3.299915. For the Hartmann flow with Ha = 3.3, a similar study has been conducted with G = 4.1667 × 10−6 . Although the minimum error of E = 0.000345 was calculated at τ = τm = 0.7 and the slightly greater error of E = 0.000464 was obtained at τ = τm = τo , the efficiency of Nt = 71,888 and the compressibility of φ = 4.4049 × 10−7 at τo are better than the efficiency of Nt = 108,991 and the compressibility of φ = 4.4120 × 10−7 at τ = τm = 0.7, correspondingly. It was also noticed that the simulation failed to converge at τ = τm = 0.6. Considering the overall balance among the stability, efficiency, and accuracy at a given lattice size, the optimal relaxation time should be used for all the simulations. The numerical results for Ha = 3.3, 50.0, and 100.0 at the optimal relaxation time are shown in Fig. 6(a) for streamwise velocities, and in Fig. 6(b) for induced magnetic fields, in which some data points were skipped for clarity. The relevant data for Ha = 3.3, 50, and 100 is listed in Table 2. E (ux ) and E (bx ) are relative errors for velocity and magnetic field, respectively. The comparisons between analytical and numerical results of Hartmann flow show that the optimal relaxation time balances the stability, efficiency, and accuracy at a given lattice size. 4.3. Channel flow around a square cylinder A flow over a blunt object has many practical applications in buildings, bridges, and thermal mass flowmeters with resistance temperature detectors. The fluid flow under certain conditions past the object creates alternating low pressure vortices on the downstream, forming a vortex shedding. In order to capture this dynamic phenomenon with the lattice Boltzmann method, a channel with a size of 3 × 1 has been selected, in which a square cylinder with a blockage ratio of 1 to 5 is vertically symmetrically located and horizontally positioned from inlet to the front square wall with a distance of 0.6. To guarantee two-dimensional laminar flows, the values of Re = 100–150 is generally used. In order to further demonstrate the versatility of the optimal relaxation time, the lattice size of 1051 × 351, τo = 0.8090, Re = 100, and Ma = 0.2549 were used to
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
(a) Streamwise velocity.
183
(b) Induced magnetic field.
Fig. 6. Comparisons between analytical and numerical results with different Hartmann numbers at the optimal relaxation time.
Fig. 7. Velocities U and V for the channel flow over a block. Table 2 Relevant data for Hartmann flow. Ha
τ
E
E (ux )
E (b x )
φ
Nt
100.0 100.0 100.0 50.0 50.0 50.0 3.3 3.3 3.3
0.7887 0.8090 0.9083 0.7887 0.8090 0.9083 0.7000 0.8090 0.9083
0.008804 0.008839 0.014164 0.005475 0.005412 0.006478 0.000345 0.000464 0.000602
0.009569 0.009818 0.014086 0.005582 0.005665 0.006799 0.000262 0.000454 0.000614
0.005828 0.004647 0.014400 0.005118 0.004512 0.005324 0.000772 0.000551 0.000472
2.2032 × 10−5 2.2030 × 10−5 2.2027 × 10−5 2.2030 × 10−5 2.2028 × 10−5 2.2024 × 10−5 4.4120 × 10−7 4.4049 × 10−7 4.4043 × 10−7
1454 1212 884 2749 2556 1853 108,991 71,888 53,405
simulate the channel flow around a square cylinder. The initial conditions of zero velocity and constant density were applied. A parabolic stream velocity profile was prescribed at the channel inlet, and the pressure at the inlet was extrapolated. The solid wall was imposed with the no-slip boundary conditions. The outlet was treated with the equilibrium open boundary conditions [17]. The direct observation quantity of velocity at the centerline was selected and plotted at the time step of 100,000 and shown in Fig. 7. The dimensionless streamwise velocity, U, first gradually decreases from 1, then drastically drops to zero as the horizontal distance approaches the front wall of the square block. As the distance away from the back wall of the block increases, the velocity further drops below zero due to the low pressure zone created behind the square block, quickly climb
184
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
Fig. 8. Vorticity (top) and pressure (bottom) images for the channel flow over a square block with Re = 100.
Fig. 9. Streaklines around a square cylinder for Re = 100.
up, and then slowly inch up over 1 due to compressibility with oscillations. Each oscillation drop in velocity corresponds to a low pressure zone. The cross stream velocity, V , is zero through the block and also oscillates corresponding to the low pressure zones shown in Fig. 7. The indirect observation quantities of vorticity and pressure were used to describe the vortex shedding. The vorticity and pressure images are generated at the time step of 100,000 and shown in Fig. 8. The most effective way to visualize the transient flows with a vortex shedding is the streaklines. The streaklines of the channel flow is shown in Fig. 9. The gray lines represent the positive vortex with a counterclockwise rotation and the black lines depict the negative vortex with a clockwise rotation which nicely matches with the vorticity image shown in Fig. 8, and concisely captures the vortex shedding. When a relaxation time of τ = 0.5263 was used for the same case, however, the simulation with the corresponding lattice size of 91 × 31 from (16) failed to converge after 471 time steps. The simulation with the refined mesh of 436 × 146 corresponding to τ = 0.6250 also blew up after 2353 iterations. When the values of both Ma and Re are fixed in nearly incompressible flows, refining lattice mesh requires an increase in relaxation time. The optimal relaxation time can, therefore, be used to select an initial lattice size. If a further refinement is needed, we have to choose a relaxation time that is greater than the optimal one. 5. Conclusions The optimal relaxation time of 0.8090 and the optimal lattice size have been proposed to balance the numerical stability, efficiency, and accuracy at a given lattice size in incompressible flow simulations with lattice Boltzmann methods. The functioned polylogarithm polynomial has been defined and its two essential properties have been derived. The polynomials are useful in expressing the lattice Boltzmann equations in different time scales and the truncation errors as well. Truncation errors within the range of 0.5 < τ < 1.0 are smaller than those at other applicable relaxation times at a given lattice size. The applications of the optimal relaxation time, through the Hartmann flow numerical simulations, have been shown to be stable, efficient, and accurate. Although the optimal relaxation collisions have only been implemented and tested with the SRT collision operator in two-dimensional cases, they are equally applicable for the LBM with the MRT collision operator
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
185
in three-dimensional domains. A further study is needed for the simulations of compressible flows with the relaxation times of τo 6 τ < ∞. Acknowledgments The author would like to express his gratitude to Professor Li-Shi Luo for constructive discussions at ICMMES-2010, Dr. Chiun Wang for the encouragement and helpful discussions, and reviewers for improving the quality of this paper. The author is also indebted to his family Jianming and Jiesi Zhao. Without them, there would not have been such an enjoyable excursion into the lattice Boltzmann methods for computational fluid dynamic problems. Appendix. Explicit functioned polylogarithm polynomials The first six orders of explicit forms of the functioned polylogarithm polynomial in (23) have been given in (26) through (31) and the seventh order through the fifteenth order polynomials have also been derived. They are
L7 (τ ) = −5040τ 7 + 15120τ 6 − 16800τ 5 + 8400τ 4 − 1806τ 3 + 126τ 2 − τ
(A.1)
L8 (τ ) = 40320τ − 141120τ + 191520τ − 126000τ + 40824τ − 5796τ + 254τ − τ
(A.2)
L9 (τ ) = −362880τ + 1451520τ − 2328480τ + 1905120τ − 834120τ + 186480τ 4 − 18150τ 3 + 510τ 2 − τ
(A.3)
8
7
6
9
5
8
4
7
3
6
2
5
L10 (τ ) = 3628800τ − 16329600τ + 30240000τ − 29635200τ + 16435440τ − 5103000τ 5 + 818520τ 4 − 55980τ 3 + 1022τ 2 − τ 10
9
8
7
6
(A.4)
L11 (τ ) = −39916800τ 11 + 199584000τ 10 − 419126400τ 9 + 479001600τ 8 − 322494480τ 7 + 129230640τ 6 − 29607600τ 5 + 3498000τ 4 − 171006τ 3 + 2046τ 2 − τ
(A.5)
L12 (τ ) = 479001600τ − 2634508800τ + 6187104000τ − 8083152000τ + 6411968640τ − 3162075840τ 7 + 953029440τ 6 − 165528000τ 5 + 14676024τ 4 − 519156τ 3 + 4094τ 2 − τ
(A.6)
L13 (τ ) = −6227020800τ + 37362124800τ − 97037740800τ + 142702560000τ − 130456085760τ 9 + 76592355840τ 8 − 28805736960τ 7 + 6711344640τ 6 − 901020120τ 5 + 60780720τ 4 − 1569750τ 3 + 8190τ 2 − τ
(A.7)
L14 (τ ) = 87178291200τ 14 − 566658892800τ 13 + 1612798387200τ 12 − 2637143308800τ 11 + 2731586457600τ 10 − 1863435974400τ 9 + 843184742400τ 8 − 248619571200τ 7 + 45674188560τ 6 − 4809004200τ 5 + 249401880τ 4 − 4733820τ 3 + 16382τ 2 − τ
(A.8)
12
11
13
10
12
9
11
8
10
L15 (τ ) = −1307674368000τ + 9153720576000τ − 28332944640000τ + 50999300352000τ − 59056027430400τ 11 + 45950224320000τ 10 − 24359586451200τ 9 + 8734434508800τ 8 − 2060056318320τ 7 + 302899156560τ 6 − 25292030400τ 5 + 1016542800τ 4 − 14250606τ 3 + 32766τ 2 − τ . 15
14
13
12
(A.9)
References [1] S. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1998) 329–364. [2] P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems, Phys. Rev. 94 (1954) 511–525. [3] D. d’Humières, Generalized lattice-Boltzmann equations, in: B.D. Shizgal, D.P. Weave (Eds.), Rarefied Gas Dynamics: Theory and Simulations, in: Prog. Astronaut. Aeronaut., vol. 159, AIAA, Washington, DC, 1992, pp. 450–458. [4] P. Lallemand, L.-S. Luo, Theory of the lattice Boltzmann method: dispersion, dissipation, isotropy, Galilean invariance and stability, Phys. Rev. E 61 (2000) 6546–6562. [5] C.K. Aidun, J.R. Clausen, Lattice-Boltzmann method for complex flows, Annu. Rev. Fluid Mech. 42 (2010) 439–472. [6] K. Huang, Statistical Mechanics, John Wiley & Sons, Inc., New York, 1987. [7] Y. Qian, D. d’Humières, P. Lallemand, Lattice BGK models for Navier–Stokes equation, Europhys. Lett. 17 (1992) 479–484. [8] D.J. Holdych, D.R. Noble, J.G. Georgiadis, R.O. Buckius, Truncation error analysis of lattice Boltzmann methods, J. Comput. Phys. 193 (2004) 595–619. [9] S. Chapman, T.G. Cowling, The Mathematical Theory of Nonuniform Gases, Cambridge University Press, Cambridge, 1970. [10] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Oxford University Press Inc., New York, 2001. [11] X. He, L.-S. Luo, Lattice Boltzmann model for the incompressible Navier–Stokes equation, J. Stat. Phys. 88 (1997) 927–944. [12] U. Ghia, K.N. Ghia, C.T. Shin, High-Re solutions for incompressible flow using the Navier–Stokes equations and a multigrid method, J. Comput. Phys. 48 (1982) 387–411. [13] R. Mei, L.-S. Luo, W. Shyy, An accurate curved boundary treatment in the lattice Boltzmann method, J. Comput. Phys. 155 (1999) 307–330. [14] R. Mei, W. Shyy, D. Yu, L.-S. Luo, Lattice Boltzmann method for 3-D flows with curved boundary, J. Comput. Phys. 161 (2000) 680–699. [15] P.J. Dellar, Lattice kinetic schemes for magnetohydrodynamics, J. Comput. Phys. 179 (2002) 95–126. [16] Z.-L. Guo, C.-G. Zheng, B.-C. Shi, Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method, Chinese Phys. 11 (2002) 366–374. [17] S. Izquierdo, P. Martnez-Lera, N. Fueyo, Analysis of open boundary effects in unsteady lattice Boltzmann simulations, Comput. Math. Appl. 58 (2009) 914–921.
Contents lists available at SciVerse ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
Optimal relaxation collisions for lattice Boltzmann methods Fuzhang Zhao ∗ APD Optima 266 Robinson Drive, Tustin, CA 92782, USA
article
info
Keywords: Functioned polylogarithm polynomial Lattice Boltzmann equations Optimal relaxation time Truncation error
abstract The optimal relaxation time of about 0.8090 has been proposed to balance the efficiency, stability, and accuracy at a given lattice size of numerical simulations with lattice Boltzmann methods. The optimal lattice size for a desired Reynolds number can be refined by reducing the Mach number for incompressible flows. The functioned polylogarithm polynomials are defined and used to express the lattice Boltzmann equations at different time scales and to analyze the impact of relaxation times and lattice sizes on truncation errors. Smaller truncation errors can be achieved when relaxation times are greater than 0.5 and less than 1.0. The steady-state lid-driven cavity flow was chosen to validate the code of lattice Boltzmann procedures. The applications of the optimal relaxation parameters numerically balance the stability, efficiency, and accuracy through Hartmann flow. The optimal relaxation time can also be used to select the initial lattice size for the channel flow over a square cylinder with a given Mach number. © 2011 Elsevier Ltd. All rights reserved.
1. Introduction The lattice Boltzmann method (LBM) has emerged as a promising alternative to the numerical scheme for high performance computational fluid dynamics (CFD) in the past twenty years. The numerical approaches based on the Boltzmann equation cover a wide range of the Knudsen number for practical applications. No effort for solving linear systems of equations is needed. The method is particularly competitive for fluid flow applications with complexities of geometries, boundaries, and interfacial dynamics [1]. The most commonly used collision operators are the single relaxation time (SRT) [2] and the multiple relaxation time (MRT) [3,4]. The LBMs can be applied to physical problems by means of nondimensionalization [5]. When both the physical system and the lattice system are equivalently scaled, three key dimensionless quantities: lattice size, relaxation time or collision frequency, and Mach number are yet to be determined for any CFD problems with a given Reynolds number. Different combinations of the key dimensionless quantities may numerically impact the stability, accuracy, and efficiency of the LBMs. The optimal relaxation time has been derived to balance the numerical properties. The optimal relaxation parameters are applicable for both the SRT and the MRT collision operators since no particular assumptions relative to collision operators have been made during their derivations. The lattice Boltzmann method, which can be derived from the kinetic theory of gases, evolved from the lattice-gas automata. The primary variable is a one-particle probability distribution function f (x, e, t ). The macroscopic flow variables such as density, velocity, pressure, temperature, and stress can be expressed as the velocity moments of the one-particle probability distribution function. The general transport equation reads
(∂t + e · ∇ x + a · ∇ e )f (x, e, t ) = Q , ∗
Tel.: +1 714 389 0228. E-mail address: [email protected].
0898-1221/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.camwa.2011.06.005
(1)
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
173
where t is the time, e is the particle velocity, x is the special location, a is the acceleration due to an external force, and Q is a collision operator. To complete the Boltzmann equation from (1), the collision operator has to be specified. Four assumptions are made: (a) only binary collisions are taken into consideration, (b) the velocity of a molecule is not correlated to its space location, (c) the effect of the external force on the collision cross section is neglected, and (d) wall effects are ignored [6]. Under these assumptions, Boltzmann related the collision operator to the one-particle probability distribution function. One of the most commonly used collision operators is the linearized SRT or BGK collision operator [2] 1 Q = − [f − f (eq) ],
(2)
τ
where τ is a relaxation time. With the Boltzmann H-theorem, the Maxwell–Boltzmann or equilibrium one-particle probability distribution function was given by f
(eq)
=
ρ
exp −
(2π RT )D/2
(e − u)2
2RT
,
(3)
where D, R, T , ρ , and u are the dimensions of space, gas constant, temperature, density, and velocities, respectively. To numerically solve for f given in (1)–(3), the velocity space e is simplified by a finite set of velocity vectors ei , which is chosen for isotropy. For the D2Q9 model, the discrete velocity vectors can be practically expressed as e0 = (0, 0), e1 = −e3 = (c , 0), e2 = −e4 = (0, c ), e5 = −e7 = (c , c ), and e6 = −e8 = (−c , c ). The lattice speed c = δ x/δ t is the ratio of a lattice spacing δ x to a time step δ t. The low Mach number expansion of (3) at a constant temperature is given by [7] fi
(eq)
9 3 3 = ρwi 1 + 2 (ei · u) + 4 (ei · u)2 − 2 (u · u) , c
2c
2c
(4)
where the weighting factors wi are specifically given by w0 = 4/9, w1 = w2 = w3 = w4 = 1/9, and w5 = w6 = w7 = w8 = 1/36. The leading moments of the truncated equilibrium distribution function (4) for incompressible flows can be evaluated as
ρ=
8
fi =
8
i=0
ρ uα =
8
fi
3 1 3
,
(5)
i=0
eiα fi =
i=1
1
(eq)
8
eiα fi
(eq)
,
(6)
i=1
c 2 ρδαβ + ρ uα uβ =
8
eiα eiβ fi
(eq)
,
(7)
i=1
c 2 ρ(δαβ uγ + δβγ uα + δγ α uβ ) =
8
eiα eiβ eiγ fi
(eq)
,
(8)
i=1
where eiα is the projection of ei on α -axis. In (6) through (8), the summation is taken from 1 to 8 because of a rest particle devised for better computational stability. For ideal gases, the equation of state relates density ρ to pressure p with the following equation p = ρ RT =
1 3
ρc2.
(9)
Eqs. (1) and (2), by introducing the collision frequency ω = 1/τ , can be discretized in space x and time t with the lattice spacing of δ x = ei δ t as fi (x + ei δ t , t + δ t ) − fi (x, t ) = −ω[fi (x, t ) − fi
(eq)
(x, t )].
(10)
This is the lattice Boltzmann equation with the BGK simplification. The lattice BGK model can be solved in two steps: first by collision and then by streaming. Therefore, the relaxation of local collision and advection of particle distribution functions evolve from one lattice to another by (eq) fi (x, t + δ t ) = (1 − ω)fi (x, t ) + ωfi (x, t ),
(11)
fi (x + ei δ t , t + δ t ) = fi (x, t + δ t ),
(12)
174
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
where fi denotes the post-collision state of the particle distribution function. Since the collision step is local and the streaming step is uniform, they can be implemented straightforwardly. The kinematic viscosity in the Navier–Stokes equations for the D2Q9, D3Q15, and D3Q19 models is
ν=
2τ − 1
c 2 δt . 6 It is worthwhile to notice that the kinematic viscosity for other models might be different.
(13)
2. Optimal relaxation parameters According to Buckingham’s rule, there are usually fewer dimensionless groups than there are physical quantities. The groups make up the real characteristic and meaningful variables of the problem. The two physical quantities of density and dynamic viscosity appear in the standard incompressible Navier–Stokes equations without the gravitation term and by introducing several nondimensional variables, the Reynolds number has been shown to be the only dimensionless group in the corresponding nondimensional form of the equations. Since all flows of the same Reynolds number are mathematically comparable, the Reynolds number will be used to derive the optimal relaxation time and the optimal lattice size. The Reynolds number is usually written as Re =
U0 L
ν
,
(14)
where U0 is the characteristic fluid velocity and L is the characteristic length that the flow is going through or around. Since the equilibrium one-particle probability distribution function (4) is isothermal, the Mach number can be given by
√ Ma =
U0
3U0
=
cs
c
,
(15)
√ where cs is a sound speed, and under an isothermal condition, we have a constant sound speed cs = c / 3. The collision frequency ω determines the fraction of particles that collide in a small volume during a time interval δ t, based on the collision interval theory [6]. The lattice size Nx implies the total number of particles in the x-coordinate. Choosing an appropriate collision frequency or a relaxation time, a lattice size, and a Mach number, and scaling nondimensional quantities of interest are crucial for the lattice Boltzmann methods to successfully apply to practical physical problems with the desired Reynolds number. As we have known, the relative error of the velocity is proportional to Ma2 , and increasing the collision frequency ω will decrease the stability and the efficiency. Thus, it is necessary to establish the relationship between the collision frequency ω and Mach number Ma. The balance among the stability, efficiency, and accuracy at a given lattice size can then be obtained. Substituting (13) into (14), using (15), and rearranging yields
ω 2
1
=
√
1+
Ma 2 3L Re δ x
.
(16)
√
Since the first dimensionless group (ω/2) is inversely related to the second dimensionless group (Ma/Re2 3L/δ x) in (16), increasing one dimensionless group will decrease the other dimensionless group. For simulations of incompressible flows with a given Re, we usually use a smaller Ma to reduce the compressibility and to enhance the accuracy at a given lattice size (L/δ x). Reducing Ma, however, will increase ω and as a result, the efficiency will be lower. As the collision frequency ω approaches 2.0, the stability will deteriorate. Considering both the stability and the efficiency, we wish to obtain a smaller collision frequency. As we have mentioned, based on (16), a smaller ω requires a greater Ma. Thus, there must be a balance among the efficiency, stability, and accuracy for LBM simulations at a given lattice size. For the balance, the two dimensionless groups in (16) must be equalized
ωo
√ =
Ma 2 3L
.
2 Re δ x Substituting (17) into (16), replacing ω with ωo , and rearranging gives
ωo2 + 2ωo − 4 = 0.
(17)
(18)
The optimal collision frequency, ωo , or the reciprocal of the optimal relaxation time, τo , can be obtained from the positive root of the quadratic equation (18)
ωo =
1
τo
√
5 − 1.
=
(19)
The accompanied equation useful for selecting the optimal lattice size, Nox , can also be obtained by rearranging (17) and using (19)
√ Nox = 1 +
L
δx
=1+
√
15 −
3 Re
12
Ma
,
(20)
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
175
and the representative optimal mean free pass is
√ λo =
√
6π ( 5 − 1)
δx 24 if the Knudsen number is defined as Kn =
λ L
=
Ma
Re
π 2
(21)
.
(22)
Based on the derivation given above, the optimal relaxation time applies to incompressible flows with all discretization models since the optimal relaxation time is independent of discretization model. The optimal lattice size, however, is dependent on discretization model because the kinematic viscosity may vary with the discretization model. For the kinematic viscosity with the form other than (13), a new equation for the optimal lattice size needs to be derived by following the procedures from the derivation of (16) to that of (20). For incompressible CFD applications with a desired Re, the Mach number must be controlled to reduce compressibility effects. Reducing Ma will increase accuracies for the observed quantities such as the velocity and the density since the relative errors for both the velocity and the density are proportional to Ma2 . For the desired Reynolds number, the selection of an optimal lattice size can be coordinated with the Mach number using (20). The ratio of the Mach number to the Reynolds number can also be expressed by the Knudsen number Kn in (22) for isothermal ideal gases. The optimal relaxation parameters can be applied to both the SRT and the MRT collision operators. In the case of the MRT collision operator, we may have the relationship of s8 = s9 = ωo = 1/τo for the D2Q9 model. The bulk viscosity has a relationship with s2 . Therefore, the bulk and shear viscosities can be independently controlled in the MRT collision operator. 3. Analytical analysis In order to study truncation errors, lattice Boltzmann equations (LBE) for different time scales will be derived based on the Chapman–Enskog expansion. A new polynomial which relates to the polylogarithm function will be defined and used. As a result, the incompressible Navier–Stokes equations along with truncation errors will be recovered. 3.1. Functioned polylogarithm polynomial Slightly different from the definition in [8], the functioned polylogarithm polynomial in τ of order m is defined as
Lm (τ ) ≡
1 (−1)m Li−m 1 − , τ −1 τ
(23)
where Li−m is the polylogarithm of a negative integer order and m is a non-negative integer. The functioned polylogarithm polynomial has many useful properties. The two essential properties have been derived. The first property is the following recursion formula d[Lm (τ )]
− τ Lm (τ ), with L0 (τ ) = 1. (24) dτ The second property of the functioned polylogarithm polynomial of consecutive order has the following inequality relationship Lm+1 (τ ) = τ (1 − τ )
|Lm (τ )| 6 |Lm+1 (τ )| m = 1, 2, . . . , 1 6 τ < ∞.
(25)
The first six orders of the polynomials have been obtained from (24)
L1 (τ ) = −τ ,
(26)
L2 (τ ) = τ (2τ − 1),
(27)
L3 (τ ) = −τ (6τ − 6τ + 1), 2
(28)
L4 (τ ) = τ (2τ − 1)(12τ − 12τ + 1), 2
(29)
L5 (τ ) = −τ (120τ − 240τ + 150τ − 30τ + 1),
(30)
L6 (τ ) = τ (2τ − 1)(360τ − 720τ + 420τ − 60τ + 1).
(31)
4
3
2
4
3
2
Their applicable maximum roots are exactly obtained as
√ τ3 =
1 2
+
3 6
√ ,
τ4 =
1 2
+
6 6
√
,
τ5 =
1 2
+
1 8
+
√
105
120
,
τ6 =
1 2
+
1 6
+
15 60
,
where the subscripts of τ , such as 3, 4, 5, 6, denote the order of the polynomials obtained from (24).
(32)
176
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
Fig. 1. Maximum and averaged applicable roots.
All the non-zero roots are positive and the maximum roots satisfy the following relationship 0.5 < τ3 < τ4 < τ5 < τ6 < · · · < τm < 1.0.
(33)
Since positive kinematic viscosities require 0.5 < τ < ∞, the values of τ1 = 0 and τ2 = 0.5 cannot be used in lattice Boltzmann calculations. The maximum roots and the averaged applicable roots for the first fifteen orders are plotted in Fig. 1. The applicable roots refer to all the roots of (23) that are greater than 0.5. Some details on the explicit higher order functioned polylogarithm polynomials are given in Appendix. The impact of using different relaxation times on the truncation errors will be demonstrated by the following Chapman–Enskog expansion analysis. 3.2. Lattice Boltzmann equations in different time scales To write the lattice Boltzmann equations in different time scales ϵ = δ t, the Chapman–Enskog expansion is usually applied [9]. With the following equations fi =
∞
ϵ n fi(n) ,
(34)
n =0
∞ ∂ ∂ = ϵn , ∂t ∂ tn n=0
(35)
∞ ϵn ∂ ∂ n fi (x + ei δ t , t + δ t ) = + ei · fi (x, t ), n! ∂ t ∂x n =0
(36)
and
we can expand the LBE (10) into a series of LBEs in different orders of the time scale ϵ . They are, by letting D0 = ∂t0 + ei · ∂x , expressed as fi fi fi fi fi
(0)
= fi(eq) ,
(1)
= L1 (τ )D0 fi(0) ,
(37) (38)
(0)
(2)
= L1 (τ )
∂ fi L2 (τ ) 2 (0) + D0 fi , ∂ t1 2!
(39)
(3)
= L1 (τ )
∂ fi(0) ∂ f (0) L3 (τ ) 3 (0) + L2 (τ )D0 i + D 0 fi , ∂ t2 ∂ t1 3!
(40)
(4)
= L1 (τ )
(0) ∂ fi(0) ∂ f (0) L3 (τ ) 2 ∂ fi L4 (τ ) 4 (0) + L2 (τ )D0 i + D0 + D0 fi + A4 , ∂ t3 ∂ t2 2! ∂ t1 4!
(41)
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
177
and for the equation of order ϵ m , one may generally have fi
(m)
= L1 (τ ) +
(0) ∂ fi(0) ∂ f (0) L3 (τ ) 2 ∂ fi + L2 (τ )D0 i + + ··· D0 ∂ tm−1 ∂ tm−2 2! ∂ tm−3
Lm−1 (τ )
(m − 2)!
−2 Dm 0
∂ fi(0) Lm (τ ) m (0) + D0 fi + Am , ∂ t1 m!
(42)
where A4 and Am are the terms including other higher order and cross derivatives for the orders ϵ 4 and ϵ m equations, respectively. 3.3. Truncation errors Since the Burnett and superBurnett equations are obtained from the third and fourth order truncations of the Chapman–Enskog expansion of the Boltzmann equation, the Boltzmann equation in general is not necessarily converged to the Navier–Stokes equations [10]. It is useful, therefore, to study the impact of different combinations of relaxation times and lattice sizes on truncation errors. Similar to the procedure in [11], the incompressible Navier–Stokes equations with truncation errors may be recovered by neglecting the divergent term and convection terms of O(Ma2 ) or higher order for momentum equations
∂ρ + ∇ · (ρ u) = RC , ∂t ∂(ρ u) + ∇ · (ρ uu) = −cs2 ∇ρ + ν∇ 2 (ρ u) + RαM , ∂t
(43) (44)
where the truncation errors RC and Rα M for continuity and momentum equations boil down to
RC = ϵ 2
+ϵ Rα M = ϵ
L2
C2,1 +
L1
m−1
2
L2 L1
L2 L1
L3 L1
C2,2 + ϵ 3
Cm−1,1 +
Mα 2,1 +
+ · · · + ϵ m−1
L3
L1 L2
L1
L3 L1
L2 L1
C3,1 +
L3
L1 Lm
Cm−1,2 + · · · +
Mα 2,2 + ϵ
3
L2 L1
C3,2 +
L1
Lm L1
C3,3 + · · ·
L1
Cm−1,m−1 ,
Mα 3,1 +
Mα m−1,1 + · · · +
L4
L3 L1
Mα 3,2 +
(45)
L4 L1
Mα 3,3
Mα m−1,m−1 ,
(46)
where Cj,k and Mα j,k for j ∈ [2, 3, . . . , m − 1] and k ∈ [1, 2, . . . , j], independent of a relaxation time, are the higher order derivatives of density, ρ , and momentum, ρ uα , respectively, with respect to different time scales. For the continuity truncation errors, we have C2,1 = −
∂[D0 (ρ)] , ∂ t1
C3,1 = −
∂[D0 (ρ)] ∂ 2ρ − , ∂ t2 2!∂ t12
C2,2 = −
D30 (ρ) 3!
,
C3,2 = −
(47)
∂[D20 (ρ)] , 2!∂ t1
C3,3 = −
D40 (ρ) 4!
,
(48)
and Cm−1,m−1 = −
Dm 0 (ρ) m!
.
(49)
For the momentum truncation errors, we get Mα 2,1 = −
∂[D0 (ρ uα )] , ∂ t1
Mα 3,1 = −
∂[D0 (ρ uα )] ∂ 2 (ρ uα ) − , ∂ t2 2!∂ t12
M α 2 ,2 = −
D30 (ρ uα ) 3!
,
Mα 3,2 = −
(50)
∂[D20 (ρ uα )] , 2!∂ t1
Mα 3,3 = −
D40 (ρ uα ) 4!
,
(51)
and Mα m−1,m−1 = −
Dm 0 (ρ uα ) m!
.
(52)
178
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
In order to reduce the high order truncation errors, we can use the averaged applicable roots at different orders shown in Fig. 1. There are many ways to use the roots of the functioned polylogarithm polynomial to construct a relaxation time. A reasonable guess is the average-weighted relaxation time, which takes the following form
τa =
τ¯3 + ϵ τ¯4 + ϵ 2 τ¯5 + · · · + ϵ m−3 τ¯m , 1 + ϵ + ϵ 2 + · · · + ϵ m−3
(53)
where τ¯m denotes the averaged applicable root in order m. The averaged roots at any order are in the range of τ3 ≈ 0.7887 and τ4 ≈ 0.9082. As the polynomial order m increases, the effects of the averaged applicable roots fade away. Thus, the average τa should be closer to the lower limit τ3 and the optimal relaxation time τo . In other words, the optimal relaxation time reduces the higher order truncation errors of the Navier–Stokes equations and therefore increases the accuracy of simulations. If higher order truncation errors were the only concern, using τ3 as a relaxation time and a refined mesh would be a good choice to reduce truncation errors, making the Boltzmann equation converge to the Navier–Stokes equations. The absolute values of the lowest order truncation errors for continuity and momentum equations take the following forms, respectively
L2 (τ ) ∂[D0 (ρ)] , |RC 2,1 | = ϵ L1 (τ ) ∂ t1 L2 (τ ) ∂[D0 (ρ uα )] . |RαM2,1 | = ϵ 2 L1 (τ ) ∂ t1 2
(54)
(55)
For incompressible flows, we can select a lattice size based on Ma and Re. Thus, the lattice spacing, ϵ , is independent of a relaxation time. The truncation errors for the continuity equation and the momentum equations share the same form in terms of ϵ and Lm (τ ). The functioned polylogarithm polynomials have the following relationship
L2 (0.5 < τ < 1) L2 (1 6 τ < ∞) < L (0.5 < τ < 1) L (1 6 τ < ∞) . 1 1
(56)
With (54), (55), and (56), we obtain
|RC 2,1 (0.5 < τ < 1, ϵ)| < |RC 2,1 (1 6 τ < ∞, ϵ)|,
(57)
|RαM2,1 (0.5 < τ < 1, ϵ)| < |RαM2,1 (1 6 τ < ∞, ϵ)|.
(58)
The inequalities (57) and (58) indicate that the truncation errors using relaxation times within the range of 0.5 < τ < 1.0 are smaller than those with other applicable relaxation times at a given lattice size. 4. Numerical analysis with applications 4.1. Steady-state lid-driven cavity flow In order to demonstrate the benefits of the optimal relaxation time and the corresponding optimal lattice size, the algorithm has been implemented with MATLAB. The two-dimensional lid-driven cavity flow has been widely adopted for a CFD code validation because plenty of benchmark data is available [12]. The computational domain is a square. The top fluid moving with the constant horizontal velocity of U0 drives the flow inside the cavity. All the other velocity boundary conditions are set to zero. In the LBM simulations, the MLS scheme [13,14] was adopted to treat all the velocity boundary conditions. For LB simulations, the criterion of a steady state is satisfied by the following Euclidean norms Eu =
∥u(x, t + δ t ) − u(x, t )∥ 6 Ea , ∥u(x, t + δ t )∥
(59)
where Ea is the allowable tolerance for errors. In the lid-driven cavity flow simulations, Ea = 10−7 is used. The isothermal LBM is a nearly incompressible scheme. When we use the scheme to model incompressible flows, the compressibility effect on a solution must be carefully monitored and controlled. The compressibility can be described by the coefficient of variation of a density
φ=
1
ρ¯
(ρ − ρ) ¯ 2 /N ,
(60)
where the average density is given by the following equation
ρ¯ =
N
ρ
.
(61)
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
179
Table 1 Input and output for two-dimensional lid-driven cavity flow.
ω
τ
U0
φ
Nt
1.800 1.600 1.400 1.236 1.000 0.800 0.600
0.5556 0.6250 0.7143 0.8091 1.0000 1.2500 1.6667
0.0144676 0.0325521 0.0558036 0.0804848 0.1302083 0.1953125 0.3038194
7.6126 × 10−5 3.1069 × 10−4 8.4017 × 10−4 1.6761 × 10−3 4.1936 × 10−3 9.1303 × 10−3 2.1379 × 10−2
– – 103,172 80,783 59,975 46,096 33,763
Fig. 2. Velocity comparisons of two-dimensional cavity flow with Re = 100.
The lid-driven cavity flow with Re = 100, N = 129 × 129 lattices, and different relaxation times were chosen for the simulations by the LBM. In order to keep the Reynolds number a constant, different velocities U0 corresponding to the selected relaxation times have been obtained. The computational efficiency is measured by the total number of time steps Nt when the steady-state condition (59) is satisfied. The input and output for the benchmark case are listed in Table 1. Generally, as the relaxation time increases, the compressibility increases but the accuracy decreases while the efficiency increases. Specifically, on the one hand, when ω = 1.8 or 1.6 is used, the LBM does not converge to the desired tolerance of Ea = 10−7 despite a relatively low compressibility of φ = 7.6126 × 10−5 and φ = 3.1069 × 10−4 , respectively. On the other hand, when ω = 0.6 is used, the result possesses a relatively high compressibility of φ = 2.1379 × 10−2 , even though it is as efficient as Nt = 33763. The results in Table 1 can also be explained in the collision interval theory. When the fraction of particles for collisions is overdetermined (ω > ωo ), it will take more time to converge and the change from f to f (0) is small and linear in a time interval. When it is underdetermined (ω < ωo ), despite less time for convergence, the change from f to f (0) is large and nonlinear in a time interval, causing an increase in compressibility. The optimal collision frequency ωo determines the balanced and representative fraction of particles in collisions during a time interval and the optimal lattice size Nox given in (20) helps us to select the total number of particles within a given domain. The velocity results from the current simulation along the vertical and horizontal lines through the geometric center for Re = 100 agree with the results obtained from the finite difference method with a multigrid solver given in [12] and they are depicted in Fig. 2. The lid-driven cavity flow, however, has no analytical solution. The accuracy can be estimated based on the Richardson extrapolation. The results for the lattice 513 × 513 has been extrapolated from those of lattices 257 × 257 and 129 × 129, and are used as an ‘‘exact’’ solution for the studied cavity flow by ue ≈ u(513×513) = [4u(257×257) − u(129×129) ]/3. The relative error of velocity can then be calculated by |u − ue |/ |ue | for the cavity flow. For accuracy estimation, three different lattices, 257 × 257, 129 × 129, and 65 × 65, with Re = 100 and ωo ≈ 1.236, were used in LBM simulations of the lid-driven cavity flow. The relative change of velocity defined in (59) as a function of the total time steps has been plotted in Fig. 3. As time continues, the percentage error drops steeply from 1.00 to 0.01 for all the three different meshes, changes slope, and converges to Ea = 10−7 at the total time steps of 26,493, 80,783, and 272,333 from meshes of 65 × 65, 129 × 129, and 257 × 257, respectively. The details depicted in the inset graph of Fig. 3 demonstrate an oscillating convergence.
180
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
Fig. 3. Iteration behaviors of two-dimensional cavity flow with Re = 100.
The horizontal velocities along the vertical line through the geometric center of the cavity were also selected for accuracy estimations. The Errors calculated are 0.00393, 0.00109, and 0.00027 for meshes 65 × 65 (Ma = 0.2788), 129 × 129 (Ma = 0.1394), and 257 × 257 (Ma = 0.0697), respectively. Thus, the relative velocity error is proven to be proportional to Ma2 . Also monitored and controlled, the corresponding compressibility effects are 6.12155 × 10−3 , 1.67611 × 10−3 , and 4.50300 × 10−4 . 4.2. Hartmann flow The magnetohydrodynamic (MHD) effects have many engineering applications such as the electromagnetic casting of metals, the confinement of plasmas, magnetic flowmeters, and so forth. Hartmann flow refers to a steady-state incompressible flow of an electrically conducting fluid through a pair of parallel plates under a transversely prescribed magnetic field. The force driving the fluid along the parallel plates is the pressure gradient. The magnetic fields induce currents in the moving conductive fluid, which create forces on the fluid, and also change the magnetic field itself. With the boundary conditions of Ux (±H ) = 0 and Bx (±H ) = 0 on the plates, the analytical solution of Hartmann flow takes the following form [15] Ux (y) Un Bx (y) Bn
coth(Ha)
=
1−
Ha 1
=
sinh(Ha y¯ ) sinh(Ha)
Ha
cosh(Ha y¯ ) cosh(Ha)
,
(62)
− y¯ ,
(63)
√
where Un = GH 2 /(ρν), Bn = GH 2 / ρνη, G = −dp/dx, y¯ = y/H, and the Hartmann number is given by B0 H Ha = √ ,
(64)
ρνη
where B0 is the transverse magnetic field, and η is the magnetic diffusivity. The detailed procedures of using the LBM to simulate MHD problems can be found in [15]. The equilibrium distribution functions for both the hydrodynamics and the magnetic field are adopted respectively as fi
′ (eq)
(eq)
gi
= fi
(eq)
+ wi
9
2c 4
1 2
|B| |ei | − (ei · B) 2
2
= wi [B + 3ei · (uB − Bu)],
2
,
(65) (66)
where the magnetic field vector B is given by B=
8 i=0
gi =
8 i =0
(eq)
gi
,
(67)
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
181
Fig. 4. Error as a function of relaxation time domain for Hartmann flow.
and the gi follows the vector Boltzmann equation (eq)
gi (x + ei δ t , t + δ t ) − gi (x, t ) = −ωm [gi (x, t ) − gi
(x, t )].
(68)
The magnetic diffusivity through the Chapman–Enskog expansion analysis of (68) can be recovered for the D2Q9 model as
η=
2τm − 1 6
c 2 δt ,
(69)
where τm is a magnetic relaxation time. By the same token, the magnetic optimal collision frequency or the inverse of the magnetic optimal relaxation time has been obtained
ωmo =
√
1
τmo
=
5 − 1.
(70)
For the constraint of ∇ · B = 0, the magnetic isothermal Mach number may be defined based on the Alfvén velocity as
√ 3 Mam = B0 √ .
(71)
ρ
Therefore, the optimal lattice size, Nmox , can be refined by reducing the magnetic Mach number (71) based on the following equation
√ Nmox = 1 +
5−1 24
2
w
Ha −1 . Mam
(72)
Under the conditions of ω = ωo = ωmo , the magnetic optimal lattice size (72) is equivalent to the hydrodynamic one (20). The lattice size of Nx × Ny = 9 × 101, the half distance between plates of H = 1, and the pressure gradient of G = 2.0833 × 10−4 were chosen for the Hartmann flow simulations. The initial conditions of both velocity, ux , and magnetic field, bx , were assumed to be zero. The boundary conditions were prescribed using the nonequilibrium extrapolation scheme [16]. In order to reduce the impact of boundary conditions on the accuracy of simulations, the second order extrapolation scheme was used for the nonequilibrium distribution functions. The pressure gradient was implemented by applying both inlet pressure and outlet pressure boundary conditions. The allowable tolerance to reach a steady state is also set to be 10−7 for both the velocity field and the magnetic field based on the Euclidean norm (59). The overall relative errors between analytical and numerical solutions for both velocity and magnetic field can be evaluated by
E=
(Ux − ux )2 + (Bx − bx )2 , (Ux )2 + (Bx )2
(73)
where ux and bx are the streamwise velocity and induction magnetic field calculated by LBM, respectively. The LBM simulations of Hartmann flow with Ha = 50 were initially conducted within the relaxation time domain of τ , τm ∈ [0.65, 1.15]. The relative error defined in (73) as a function of the relaxation time domain is plotted in Fig. 4. Better results have been achieved when equalized relaxation times, τ = τm , are used. The accuracy as a function of equalized relaxation times has then been plotted in Fig. 5. When the optimal relaxation times of τ = τm = 0.8090, rather than τ3 or
182
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
Fig. 5. Error as a function of equalized relaxation times for Hartmann flow.
τ4 , are used, the error of E = 0.005412 is the least among the selected relaxation times for the simulations. For the extreme cases of τ = τm = 0.55 and 2.00, the simulations failed to converge. For the Hartmann flow with Ha = 100, a similar study has been conducted. Although the minimum error of E = 0.0088039 was calculated at τ = τm = τ3 and the slightly greater error of E = 0.0088396 was obtained at τ = τm = τo , the efficiency of Nt = 1212 and the compressibility of φ = 2.2030 × 10−5 at τo is better than the efficiency of Nt = 1454 and the compressibility of φ = 2.2032 × 10−5 at τ3 , correspondingly. As we have known, the induced magnetic field strength tends to be zero as Ha → 0 or Ha → ∞, meaning that the extreme values of Bx must exist at a Hartmann number of Ha ∈ (0, ∞). The minimum value as a function of Ha can be obtained Bx,min Bn
=
1
Ha
sinh(Ha y¯ s ) sinh(Ha)
− y¯ s ,
(74)
where the stationary point as a function of Ha is given below
y¯ s =
1 Ha
ln
sinh(Ha)
Ha
+
sinh(Ha) Ha
2
−1 .
(75)
Based on (74) and (75), the minimum induced magnetic field has the value of Bx,min /Bn ≈ −0.102818 located at y¯ s ≈ 0.633185 for Ha ≈ 3.299915. For the Hartmann flow with Ha = 3.3, a similar study has been conducted with G = 4.1667 × 10−6 . Although the minimum error of E = 0.000345 was calculated at τ = τm = 0.7 and the slightly greater error of E = 0.000464 was obtained at τ = τm = τo , the efficiency of Nt = 71,888 and the compressibility of φ = 4.4049 × 10−7 at τo are better than the efficiency of Nt = 108,991 and the compressibility of φ = 4.4120 × 10−7 at τ = τm = 0.7, correspondingly. It was also noticed that the simulation failed to converge at τ = τm = 0.6. Considering the overall balance among the stability, efficiency, and accuracy at a given lattice size, the optimal relaxation time should be used for all the simulations. The numerical results for Ha = 3.3, 50.0, and 100.0 at the optimal relaxation time are shown in Fig. 6(a) for streamwise velocities, and in Fig. 6(b) for induced magnetic fields, in which some data points were skipped for clarity. The relevant data for Ha = 3.3, 50, and 100 is listed in Table 2. E (ux ) and E (bx ) are relative errors for velocity and magnetic field, respectively. The comparisons between analytical and numerical results of Hartmann flow show that the optimal relaxation time balances the stability, efficiency, and accuracy at a given lattice size. 4.3. Channel flow around a square cylinder A flow over a blunt object has many practical applications in buildings, bridges, and thermal mass flowmeters with resistance temperature detectors. The fluid flow under certain conditions past the object creates alternating low pressure vortices on the downstream, forming a vortex shedding. In order to capture this dynamic phenomenon with the lattice Boltzmann method, a channel with a size of 3 × 1 has been selected, in which a square cylinder with a blockage ratio of 1 to 5 is vertically symmetrically located and horizontally positioned from inlet to the front square wall with a distance of 0.6. To guarantee two-dimensional laminar flows, the values of Re = 100–150 is generally used. In order to further demonstrate the versatility of the optimal relaxation time, the lattice size of 1051 × 351, τo = 0.8090, Re = 100, and Ma = 0.2549 were used to
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
(a) Streamwise velocity.
183
(b) Induced magnetic field.
Fig. 6. Comparisons between analytical and numerical results with different Hartmann numbers at the optimal relaxation time.
Fig. 7. Velocities U and V for the channel flow over a block. Table 2 Relevant data for Hartmann flow. Ha
τ
E
E (ux )
E (b x )
φ
Nt
100.0 100.0 100.0 50.0 50.0 50.0 3.3 3.3 3.3
0.7887 0.8090 0.9083 0.7887 0.8090 0.9083 0.7000 0.8090 0.9083
0.008804 0.008839 0.014164 0.005475 0.005412 0.006478 0.000345 0.000464 0.000602
0.009569 0.009818 0.014086 0.005582 0.005665 0.006799 0.000262 0.000454 0.000614
0.005828 0.004647 0.014400 0.005118 0.004512 0.005324 0.000772 0.000551 0.000472
2.2032 × 10−5 2.2030 × 10−5 2.2027 × 10−5 2.2030 × 10−5 2.2028 × 10−5 2.2024 × 10−5 4.4120 × 10−7 4.4049 × 10−7 4.4043 × 10−7
1454 1212 884 2749 2556 1853 108,991 71,888 53,405
simulate the channel flow around a square cylinder. The initial conditions of zero velocity and constant density were applied. A parabolic stream velocity profile was prescribed at the channel inlet, and the pressure at the inlet was extrapolated. The solid wall was imposed with the no-slip boundary conditions. The outlet was treated with the equilibrium open boundary conditions [17]. The direct observation quantity of velocity at the centerline was selected and plotted at the time step of 100,000 and shown in Fig. 7. The dimensionless streamwise velocity, U, first gradually decreases from 1, then drastically drops to zero as the horizontal distance approaches the front wall of the square block. As the distance away from the back wall of the block increases, the velocity further drops below zero due to the low pressure zone created behind the square block, quickly climb
184
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
Fig. 8. Vorticity (top) and pressure (bottom) images for the channel flow over a square block with Re = 100.
Fig. 9. Streaklines around a square cylinder for Re = 100.
up, and then slowly inch up over 1 due to compressibility with oscillations. Each oscillation drop in velocity corresponds to a low pressure zone. The cross stream velocity, V , is zero through the block and also oscillates corresponding to the low pressure zones shown in Fig. 7. The indirect observation quantities of vorticity and pressure were used to describe the vortex shedding. The vorticity and pressure images are generated at the time step of 100,000 and shown in Fig. 8. The most effective way to visualize the transient flows with a vortex shedding is the streaklines. The streaklines of the channel flow is shown in Fig. 9. The gray lines represent the positive vortex with a counterclockwise rotation and the black lines depict the negative vortex with a clockwise rotation which nicely matches with the vorticity image shown in Fig. 8, and concisely captures the vortex shedding. When a relaxation time of τ = 0.5263 was used for the same case, however, the simulation with the corresponding lattice size of 91 × 31 from (16) failed to converge after 471 time steps. The simulation with the refined mesh of 436 × 146 corresponding to τ = 0.6250 also blew up after 2353 iterations. When the values of both Ma and Re are fixed in nearly incompressible flows, refining lattice mesh requires an increase in relaxation time. The optimal relaxation time can, therefore, be used to select an initial lattice size. If a further refinement is needed, we have to choose a relaxation time that is greater than the optimal one. 5. Conclusions The optimal relaxation time of 0.8090 and the optimal lattice size have been proposed to balance the numerical stability, efficiency, and accuracy at a given lattice size in incompressible flow simulations with lattice Boltzmann methods. The functioned polylogarithm polynomial has been defined and its two essential properties have been derived. The polynomials are useful in expressing the lattice Boltzmann equations in different time scales and the truncation errors as well. Truncation errors within the range of 0.5 < τ < 1.0 are smaller than those at other applicable relaxation times at a given lattice size. The applications of the optimal relaxation time, through the Hartmann flow numerical simulations, have been shown to be stable, efficient, and accurate. Although the optimal relaxation collisions have only been implemented and tested with the SRT collision operator in two-dimensional cases, they are equally applicable for the LBM with the MRT collision operator
F. Zhao / Computers and Mathematics with Applications 65 (2013) 172–185
185
in three-dimensional domains. A further study is needed for the simulations of compressible flows with the relaxation times of τo 6 τ < ∞. Acknowledgments The author would like to express his gratitude to Professor Li-Shi Luo for constructive discussions at ICMMES-2010, Dr. Chiun Wang for the encouragement and helpful discussions, and reviewers for improving the quality of this paper. The author is also indebted to his family Jianming and Jiesi Zhao. Without them, there would not have been such an enjoyable excursion into the lattice Boltzmann methods for computational fluid dynamic problems. Appendix. Explicit functioned polylogarithm polynomials The first six orders of explicit forms of the functioned polylogarithm polynomial in (23) have been given in (26) through (31) and the seventh order through the fifteenth order polynomials have also been derived. They are
L7 (τ ) = −5040τ 7 + 15120τ 6 − 16800τ 5 + 8400τ 4 − 1806τ 3 + 126τ 2 − τ
(A.1)
L8 (τ ) = 40320τ − 141120τ + 191520τ − 126000τ + 40824τ − 5796τ + 254τ − τ
(A.2)
L9 (τ ) = −362880τ + 1451520τ − 2328480τ + 1905120τ − 834120τ + 186480τ 4 − 18150τ 3 + 510τ 2 − τ
(A.3)
8
7
6
9
5
8
4
7
3
6
2
5
L10 (τ ) = 3628800τ − 16329600τ + 30240000τ − 29635200τ + 16435440τ − 5103000τ 5 + 818520τ 4 − 55980τ 3 + 1022τ 2 − τ 10
9
8
7
6
(A.4)
L11 (τ ) = −39916800τ 11 + 199584000τ 10 − 419126400τ 9 + 479001600τ 8 − 322494480τ 7 + 129230640τ 6 − 29607600τ 5 + 3498000τ 4 − 171006τ 3 + 2046τ 2 − τ
(A.5)
L12 (τ ) = 479001600τ − 2634508800τ + 6187104000τ − 8083152000τ + 6411968640τ − 3162075840τ 7 + 953029440τ 6 − 165528000τ 5 + 14676024τ 4 − 519156τ 3 + 4094τ 2 − τ
(A.6)
L13 (τ ) = −6227020800τ + 37362124800τ − 97037740800τ + 142702560000τ − 130456085760τ 9 + 76592355840τ 8 − 28805736960τ 7 + 6711344640τ 6 − 901020120τ 5 + 60780720τ 4 − 1569750τ 3 + 8190τ 2 − τ
(A.7)
L14 (τ ) = 87178291200τ 14 − 566658892800τ 13 + 1612798387200τ 12 − 2637143308800τ 11 + 2731586457600τ 10 − 1863435974400τ 9 + 843184742400τ 8 − 248619571200τ 7 + 45674188560τ 6 − 4809004200τ 5 + 249401880τ 4 − 4733820τ 3 + 16382τ 2 − τ
(A.8)
12
11
13
10
12
9
11
8
10
L15 (τ ) = −1307674368000τ + 9153720576000τ − 28332944640000τ + 50999300352000τ − 59056027430400τ 11 + 45950224320000τ 10 − 24359586451200τ 9 + 8734434508800τ 8 − 2060056318320τ 7 + 302899156560τ 6 − 25292030400τ 5 + 1016542800τ 4 − 14250606τ 3 + 32766τ 2 − τ . 15
14
13
12
(A.9)
References [1] S. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1998) 329–364. [2] P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems, Phys. Rev. 94 (1954) 511–525. [3] D. d’Humières, Generalized lattice-Boltzmann equations, in: B.D. Shizgal, D.P. Weave (Eds.), Rarefied Gas Dynamics: Theory and Simulations, in: Prog. Astronaut. Aeronaut., vol. 159, AIAA, Washington, DC, 1992, pp. 450–458. [4] P. Lallemand, L.-S. Luo, Theory of the lattice Boltzmann method: dispersion, dissipation, isotropy, Galilean invariance and stability, Phys. Rev. E 61 (2000) 6546–6562. [5] C.K. Aidun, J.R. Clausen, Lattice-Boltzmann method for complex flows, Annu. Rev. Fluid Mech. 42 (2010) 439–472. [6] K. Huang, Statistical Mechanics, John Wiley & Sons, Inc., New York, 1987. [7] Y. Qian, D. d’Humières, P. Lallemand, Lattice BGK models for Navier–Stokes equation, Europhys. Lett. 17 (1992) 479–484. [8] D.J. Holdych, D.R. Noble, J.G. Georgiadis, R.O. Buckius, Truncation error analysis of lattice Boltzmann methods, J. Comput. Phys. 193 (2004) 595–619. [9] S. Chapman, T.G. Cowling, The Mathematical Theory of Nonuniform Gases, Cambridge University Press, Cambridge, 1970. [10] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Oxford University Press Inc., New York, 2001. [11] X. He, L.-S. Luo, Lattice Boltzmann model for the incompressible Navier–Stokes equation, J. Stat. Phys. 88 (1997) 927–944. [12] U. Ghia, K.N. Ghia, C.T. Shin, High-Re solutions for incompressible flow using the Navier–Stokes equations and a multigrid method, J. Comput. Phys. 48 (1982) 387–411. [13] R. Mei, L.-S. Luo, W. Shyy, An accurate curved boundary treatment in the lattice Boltzmann method, J. Comput. Phys. 155 (1999) 307–330. [14] R. Mei, W. Shyy, D. Yu, L.-S. Luo, Lattice Boltzmann method for 3-D flows with curved boundary, J. Comput. Phys. 161 (2000) 680–699. [15] P.J. Dellar, Lattice kinetic schemes for magnetohydrodynamics, J. Comput. Phys. 179 (2002) 95–126. [16] Z.-L. Guo, C.-G. Zheng, B.-C. Shi, Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method, Chinese Phys. 11 (2002) 366–374. [17] S. Izquierdo, P. Martnez-Lera, N. Fueyo, Analysis of open boundary effects in unsteady lattice Boltzmann simulations, Comput. Math. Appl. 58 (2009) 914–921.