Lattice Boltzmann model for three-dimensional decaying homogeneous isotropic turbulence

Physics Letters A 373 (2009) 1368–1373

Contents lists available at ScienceDirect

Physics Letters A www.elsevier.com/locate/pla

Lattice Boltzmann model for three-dimensional decaying homogeneous isotropic turbulence Hui Xu a , Wenquan Tao a,∗ , Yan Zhang b a b

School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an 710049, PR China Faculty of Science, Xi’an Jiaotong University, Xi’an 710049, PR China

a r t i c l e

i n f o

Article history: Received 9 March 2008 Received in revised form 9 September 2008 Accepted 18 January 2009 Available online 4 February 2009 Communicated by F. Porcelli PACS: 47.27.Gs 47.27.E-

a b s t r a c t We implement a lattice Boltzmann method (LBM) for decaying homogeneous isotropic turbulence based on an analogous Galerkin filter and focus on the fundamental statistical isotropic property. This regularized method is constructed based on orthogonal Hermite polynomial space. For decaying homogeneous isotropic turbulence, this regularized method can simulate the isotropic property very well. Numerical studies demonstrate that the novel regularized LBM is a promising approximation of turbulent fluid flows, which paves the way for coupling various turbulent models with LBM. © 2009 Elsevier B.V. All rights reserved.

Keywords: Lattice Boltzmann models Decaying homogeneous isotropic turbulence Taylor microscale Isotropic ratios Skewness Flatness

1. Introduction In the past decades, lattice Boltzmann method has become a powerful and efficient tool to solve various problems of fluid flows and drawn much attention of researchers in the computational fluid dynamics community [1]. Recently, LBM has been utilized to simulate turbulent fluid flow [1–4], and the results offer us some fundamental understanding for the turbulence [5]. Traditionally, Navier–Stokes equations (NSE) have been utilized to analyze and predict turbulent fluid flows. It is a well-accepted understanding that apart from some deterministic parameters, turbulent flow is also characterized by its statistical behavior [6,7]. The Boltzmann equation describes the phenomena of fluid flows and turbulence in the statistical physics framework, which possesses some inherent consistency with the physical processes of turbulent flows. Therefore LBM is of great scientific significance for study of turbulence in that it can redound to our better understanding of the meaning of turbulence or filtered turbulence. In the development of turbulence computational model and method valid over a wide range of flow parameters LBM can play an important role [5,6].

*

Corresponding author. Tel.: +86 29 82663232 802; fax: +86 29 82669106. E-mail address: [email protected] (W. Tao).

0375-9601/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.physleta.2009.01.058

In this Letter, we concentrate our attention on constructing a scheme to describe statistical turbulence on the basis of the orthogonal Hermite polynomial space, and reconsider the Reynolds stress tensor by the kinetic description. Generally speaking, based on the lattice Boltzmann scheme, Naver–Stokes equations can be recovered by the second-order orthogonal Hermite polynomial expansion of equilibrium distribution for the isothermal system. It is well known that for a higher-order orthonormal Hermite expansion, we can gain a representation of hydrodynamics beyond the Navier–Stokes equations [8]. For most incompressible turbulent flow systems represented by Navier–Stokes equations, however, it is unnecessary to turn to higher-order orthonormal Hermite expressions. Thus we will still keep the second-order orthonormal Hermite expressions, while examine the restrictions posed by the present treatment widely adopted in references and search for its improvement. When we recover Navier–Stokes equation by kinetic theory, there exists a restriction which makes us implementing computations in a finite second-order orthonormal Hermite polynomial subspace. However, for most engineering problems the fluid flow system is far from equilibrium, hence this restriction is often broken down and some unstable factors are introduced. Turbulent fluid flow is a typical example in which the nonlinear actions make system far from equilibrium. Under this circumstance, the single particle distribution function does not automatically lie

H. Xu et al. / Physics Letters A 373 (2009) 1368–1373

Table 1 Gauss–Hermite quadrature formulae in two and three dimensions for E 319,5 . The subscript FS denotes a fully symmetric set of points. Quadrature

ξα

E 319,5

(0, 0, 0)

ω(ξα )

P

√ (√3, 0 √, 0)FS ( 3, 3, 0)FS

1

1/3

6 12

1/18 1/36

in the finite second-order orthonormal Hermite polynomial subspace [9]. According to the idea of Galerkin method, this will lead to some instabilities and errors of computation. In this Letter, for incompressible decaying homogeneous isotropic turbulence (DHIT), we project the related distributions into a corresponding orthonormal Hermite polynomial subspace and implement our computations. This will be called Galerkin regularization. In fact, some regularization procedures have been adopted to regularize the distribution for some fluid flow problems [9,10], but these procedures differ from the Galerkin regularization and have not been utilized to make researches on the fundamental quantities (e.g. isotropy ratios) of three-dimensional (3D) turbulent flows. In the theoretical studies of turbulence, DHIT is a typical example often used to validate the models of turbulence [2,11]. LBM possesses the natural parallel characteristics, and hence it is an effective tool for computational fluid flows and furthermore for the turbulence. It is essential to implement the present (original) LBM for the fundamental isotropic researches of DHIT, investigate its shortcomings (if any) and to improve it. From the numerical results in [12], the original LBM seems not to simulate the isotropic quantities very well for DHIT. So, the study for improving the effectiveness of LBM is very significative for turbulence investigation. In the numerical study of this Letter we propose a new method for modeling turbulence by an improved LBM which will be called regularized LBM hereafter, and the computations performed in this Letter demonstrate that the regularized LBM is very useful and can be used to capture turbulence effectively. 2. Basic background of LBM Considering a general fluid flow in D-dimensional space (D = 2, 3), we take following lattice BGK formulation (LBGK) of fluid dynamics by the discrete distribution [9] f i (x + δt ξi , t + δt ) − f i (x, t ) = Ωi ,

(1)

Ωi = −

(2)

1

τ

 0

fi − fi

i = 0, . . . , d,

+ Fi,

where τ is the relaxation time, f i0 is the truncated Hermite expansion of the Maxwell–Boltzmann distribution at ξi and F i is the contribution from the body force. {ξi } are the abscissas of a sufficiently accurate Gauss–Hermite quadrature listed in Table 1. According to the Chapman–Enskog analysis, we can recover Navier–Stokes equation when the second-order moment terms are retained. And for the isothermal system, the corresponding equilibrium term exhibits the following form [8]

 f i0 = ω(ξi )ρ 1 + ξi · u +

1 2

(ξi · u )2 − u 2



 ,

(3)

where ω(ξi ) is the weighting factor which satisfies the following equation d 

ω(ξi ) = 1.

(4)

Hermite polynomials are a set of complete orthonormal basis of the Hilbert space with the following inner product [13]



 f , g =

ω f g dξ,

(5)

∞ where f and g are the bases of the space H∞ = {H(m) }m =0 and ( m) ω is the weight function. H is the mth Hermite polynomial as follows

H(0) = 1, (1)

Hi

(6a)

= ξi ,

(6b)

(2)

Hi j = ξi ξ j − δi j .

(6c)

The elements in the space H∞ satisfy the following orthonormal relation



ωH(i m) H(jn) dξ = δmn δij ,

(7)

where δij = 1 if i = {i 1 , . . . , im } is a permutation of j = { j 1 , . . . , jn }, and otherwise δij = 0. For any function f , we can approximately represent it by the space H N N  1 (n) (n) a H , n!

f ≈ ω(ξ )

(8)

n=0

where the parameter a(n) can be calculated by



f H(n) dξ.

a(n) =

(9)

For the standard LBM scheme, the equilibrium distribution (3) is represented by the Maxwell–Boltzmann equilibrium expansion in the space H2 for isothermal systems. 3. Galerkin regularization of distribution functions From Section 2, we realize that all of the distribution f α in Eq. (1) should be calculated in the space H2 , that is to say, the continuous Boltzmann equation with BGK collision term should be projected into the space H2 and then the computations are performed. When the calculation is implemented, because of the effects of nonlinear terms, the distribution f α will beyond the space H2 [9] and some unstable factors are resulted for fluid flows of low viscosity. When the flow system is not far from equilibrium, such error is small and ignorable. However, when it is far from equilibrium the induced error cannot be ignored [9]. In order to implement the regularization method for turbulence modeling, we introduce the following approximation ˜f α for f α

˜f i = ω(ξi )

2  1 (n) (n) a H (ξi ), n!

(10)

n=0

where a(n) is defined as follows [8,9,13] a(n) =

d 

f i H(n) (ξi ).

(11)

i =0

According to the Galerkin interpretation, ˜f i is the filtering result of neq f i in the space H2 . The nonequilibrium part f i of distribution f i with respect to truncated Maxwellian can also be filtered and calculated by the following formula in the space H2 [9],

˜f neq = ω(ξi ) H i

(2) (ξ

2

i =0

In fact, the above formula is described by the truncation of the Hermite polynomials up to second-order kinetic moments. The

1369

i)

d 

neq

fn

ξn ξn .

(12)

n=0

The post-streaming distribution can be split into two parts neq

fi = fi

+ f i(0) .

(13)

1370

H. Xu et al. / Physics Letters A 373 (2009) 1368–1373

It is to be noted that in the space H2 , f i0 = ˜f i0 . Thus Eq. (13) can be redefined as follows:

˜f i = ˜f neq + f (0) . i i

(14)

This implies that Eq. (10) can be calculated via Eq. (12) very easily. It is known that in the filtering process of the nonequilibrium terms only the nonconservation fields should be changed while mass and momentum remain invariant. According to the above discussion, we have d 

fi =

d 

i =0 d 

(15)

ξi ˜f i (u = u˜ ).

(16)

i =0

But the following formula is not admitted d 

ξi α ξi β f i =

d 

i0

Because f i0 is included in the space H2 ( f i0 = ˜f i0 ), we project Eq. (1) into the space H2 and gain the following LBM scheme without regard to the body force

˜f i (x + δt ξi , t + δt ) − ˜f i (x, t ) = Ω˜ i ,

τ

(18)

i = 0, . . . , d.

(19)

In the modeling of turbulence the strain rate must be calculated. For LBM, the strain rate related to the nonequilibrium part can be calculated by two forms 3

neq

2ρτ 3

S˜ α β = −

2ρτ

Πα β ,

neq

Πα β =



neq

ξi α ξi β f i

,

(20)

ξi α ξi β ˜f i

,

(21)

i

neq Π˜ α β ,

neq Π˜ α β =



neq

i

where S α β and S˜ α β dedicate S αβ =

1 2

(∂α u β + ∂β u α ),

S˜ α β =

1 2

(∂α u˜ β + ∂β u˜ α ).

(22)

In order to construct turbulence model for the current LBM model (18), we introduce the following filter or average bracket (23)

˜ i, F˜ i (x + δt ξi , t + δt ) − F˜ i (x, t ) = Ω

where

1





F˜ i − ˜f i0 ,

 f i0 

(24) (25)

 |u˜ |2 , 2

 ε = 2ν | S˜ α β |2 .

κ=

(30) (31)

(32)

we can rewrite Eq. (27) as follows





f i0 = F˜ i0 +

1 2

ω(ξi )ρ Γiα β σα β .

(33)

According to the above discussion, we can rewrite the R.H.S. of Eq. (24) by the following collision term for the case without the body force

where



1 

 1 1 F˜ i − F˜ i0 , ∗ = 1− ∗ τ τ τ

R.H.S. =

ω(ξi )ρ Γiα β σα β

2 F˜ i

neq

,

(34)

τ ∗ is the effective relaxation time.

4. Numerical studies In this section, we implement the improved LBM simulation scheme presented above for the isotropy ratios of decaying homogeneous isotropic turbulence (DHIT). For DHIT our focus is to investigate the basic quantities by our regularization method. In the simulation, it is assumed that the initial incompressible velocity field u 0 satisfies the following 3D energy spectrum [12,14] E (k, 0) =



q20 kσ

σ

exp − 2 A kσp 2



k kp

2  ,

(35)

where q (t ) =



u 2i



∞ =2

E (k, t ) dk.

(36)



2

The value of q2 is such that the Mach number [12] M=

q

 1.

cs

(37)

And A possesses the following expression [12,14]

∞



kσ exp −

A=

has the following form



0 ˜f = ω(ξi )ρ 1 + ξi · u˜  + 1 (ξi α ξi β − δα β )u˜ α u˜ β  , i

σ 2

 k2 .

(38)

0

(26)

where the repeat indexes denote the Einstein summation. We can rewrite the above formula as follows



κ and ε are turbulence kinetic energy and dissipation rate

0

Now, making the bracket acting on Eq. (18), we get the following coarse-graining LBM scheme

τ

(29)

1

2

· : ˜f i → F˜ i .

Ω˜ i = −

where

κ2 ˜

S , ε αβ

Γi α β = ξi α ξi β − δα β , (17)

i0

S αβ = −

2

By introducing a second-order tensor

ξi α ξi β ˜f i (α = β).

 1 Ω˜ i = − ˜f i − f i0 ,

(28)

In the above equation, σα β = u˜ α u˜ β  and u˜ denotes the fluctuation velocity. For example, we could choose the following Reynolds tress tensor with κ − ε turbulence models 3

˜f i (ρ = ρ˜ ),

d 

i =0

u˜ α u˜ β  = u˜ α u˜ β  + u˜ α u˜ β .

σα β = − κ δα β + cμ

i =0

ξi f i =

where σα β is the Reynolds stress tensor. In fact, similar to the Reynolds average of turbulence modeling we have the following equality under the average meaning

     ˜f 0 = ω(ξi )ρ 1 + ξi · u˜  + 1 (ξi α ξi β − δα β ) u˜ α u˜ β  + σα β , i 2

(27)

In the numerical study, the initial velocity field u 0 in the physical space is easy to be obtained by the inverse Fourier transform related with the energy spectrum [14]. Here, a pressure field consistent with the initial velocity is required. Different methods for obtaining the pressure field are proposed in Refs. [12,14–16]. The iterative solution procedure for the pressure Possion’s equation is adopted [14]. The initial energy spectrum adopted is in the range of wave number [kmin , kmax ] = [2π / L , N π / L ] [12,17], where L is

H. Xu et al. / Physics Letters A 373 (2009) 1368–1373

the characteristic length scale and N determines the mesh resolution. Here, the simulation is implemented in a cubic box [0, 2π ]3 . According to the energy spectrum, the integral length scale is estimated by [18]

∞

3π

L (t ) =

k−1 E (k, t ) dk

∞

0

4

0

E (k, t ) dk

5

∞

λ(t ) =

(40)

.

Numerically, the Taylor microscale is calculated by

λ=

u 2i 1/2 (∂ u i /∂ xi )2 1/2

(41)

(the repeated index i does not imply summation). The Taylor microscale Reynolds number is defined by Reλ =

uλ

ν

kp

σ

λ

ν

u (0)

L (0)kmin

ηkmax

Reλ (0)

8

2

0.2282

0.000010358

0.01

0.33

1.0

220

Table 3 Basic quantities of DHIT.

1/2

E (k, t )

∞ 0 2 0 k E (k, t ) dk

Table 2 Main parameters for the simulation of DHIT.

(39)

,

and the Taylor microscale is expressed by



1371

t

I1

I2

I3

S

K 1 (0)

K 2 (0)

K 3 (0)

1 2 3 4 5 6 7 8

0.97 0.97 0.98 0.97 0.99 0.97 0.98 0.96

0.96 0.98 0.97 0.96 0.98 0.97 0.98 0.97

1.02 1.01 1.03 1.02 1.01 1.03 1.02 1.03

0.51 0.48 0.49 0.50 0.51 0.51 0.50 0.49

2.97 3.00 3.02 3.02 2.94 3.01 2.94 3.00

2.96 2.98 2.92 2.95 2.99 3.02 2.97 3.02

2.99 2.99 2.97 3.03 3.00 2.98 2.92 3.05

(42)

,

where u is the average fluctuation velocity defined by

 1/ 2

 ∞ 2

u (t ) =

E (k, t ) dk

3

(43)

.

0

For any direct numerical simulation of turbulence flow, adequate resolution of the scale is of primary importance. In isotropic decay the resolution of the small scales is judged by the shape of the spectrum at high wave numbers [19]. It requires that the numerical cut-off wave number be at or beyond kmax η ≈ 1 [19], where η = (ν 3 / )1/4 is the Kolmogorov length scale, and  is the dissipation rate of energy is denoted by

∞ k2 E (k, t ) dk.

 (t ) = 2ν

(44)

0

Assuming isotropy, Eq. (44) is estimated by



 = 15ν

∂ ui ∂ xi

Fig. 1. Temporal Taylor micro-scale.

2  (45)

(the repeated index i does not imply summation). If η is chosen, the kinematic viscosity can be fixed by the definition of the Kolmogorov scale. The isotropy ratios are defined by the following formula [12]





I i = 3 u 2i /q2

(46)

(the repeated index i does not imply summation). It is known that for isotropic flows the skewness can be directly related to the production of dissipation [19]. The skewness is also a measure of the nonlinearity of the Navier–Stokes equations and is a characteristic of the turbulence flows. The skewness of the velocity derivatives is defined by Si = −

∂ u i /∂ xi 3 . (∂ u i /∂ xi 2 )3/2

(47)

The averaging skewness is defined by S = ( S 1 + S 2 + S 3 )/3. The flatness factor of ∂

(48) n

(∂ u i /∂ xni )4  . (∂ n u i /∂ xni )2 2

u i /∂ xni

is defined by

n

K i (n) =

(49)

The main parameters adopted in our simulation are given in Table 2. The initial lattice size is 2563 . The time t is normalized by 2π /u (0) [20].

In Table 3, the predicted values of basic quantities are given. It is clear that the skewness S of the velocity derivative is close to 0.5. The up-to-date results of the velocity derivative skewness are shown experimentally by Burattini et al. [21]. The current simulated results agree with the experimental results very well. The deviation of the current simulation is about ±1% and the deviation of literature [21] is about ±5%. The isotropic ratios of the velocity is given in Table 3. The isotropy ratio values of I 1 , I 2 and I 3 are within about ±2% from 1. The flatness factor values of the velocity are within about ±0.5% from 3. It is known that the isotropic ratio values should be equal to 1 and the flatness values of the velocity should be equal to 3 for DHIT [22]. From Table 2, it is known that the initial Taylor micro-scale λ = 0.2282. In Fig. 1, the values of Taylor micro-scale λ are given at different time. From the numerical results, it is known that λ declines at the beginning of simulation. Along the evolution of time, λ starts to grow after some decline. According to the theory of DHIT, Taylor micro-scale λ2 ∝ t [18]. This property is easy to be observed in Fig. 1. In Fig. 2, the three-dimensional kinetic energy spectrums is given in log–log axis at t = 1 and t = 5. In these figures, the slope values of straight lines in log–log axis are equal to −5/3. It is clear that the power laws of κ −5/3 are found in the current simulation as predicted by Kolmogorov [22]. In Fig. 3, the iso-vorticity surface shows the fine scale worm like vortical structures presented in the fully developed turbulence fluid flows. Meanwhile, the con-

1372

H. Xu et al. / Physics Letters A 373 (2009) 1368–1373

(a)

(a)

(b)

(b)

Fig. 2. Three-dimensional kinetic energy spectrum: (a) at t = 1; (b) at t = 5.

(c)

Fig. 3. Iso-vorticity surface at t = 5:

√

ω2 = 0.244.

Fig. 4. Contour plot of vorticity on a plane at t = 5: (a) the plane at x = π ; (b) the plane at y = π ; (c) the plane at z = π .

H. Xu et al. / Physics Letters A 373 (2009) 1368–1373

tour plot slices of vorticity are given at the mid vertical surfaces of the different axes in Fig. 4. From the numerically simulated results, it is shown that the model in this Letter gives the well-predicted values of DHIT. These results demonstrate that the current model can be used to simulate the turbulence fluid flows very well. In this Letter, we implement the 3D simulation of the classical DHIT problem with the regularized LBM to investigate the basic quantities. The results show that the regularized LBM can hold the isotropy ratios closer to 1 than original LBM. The regularized LBM shows some reliable perspectives for simulating turbulence, such as κ − ε model. These results can pave the way to reconsider the LBM for turbulence problem, especially for DNS. Acknowledgements This work has been supported by Chinese NSF Grant No. 50636050. The suggestions of Chief Scientist Prof. Hudong Chen in Exa Corporation are kindly acknowledged by the first author. The first author thanks Dr. Weiqin Li for fruitful discussions on FFT. The authors thank the referees for their valuable comments and suggestions that helped to improve the results of the article. The simulation is performed in National High Performance Computing Center in Xi’an Jiaotong University.

1373

References [1] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Oxford Univ. Press, Oxford, UK, 2001. [2] S.A. Orszag, H. Chen, S. Succi, J. Latt, J. Sci. Comput. 28 (213) (2006) 459. [3] H. Chen, S.A. Orszag, I. Staroselsky, S. Succi, J. Fluid Mech. 517 (2004) 301. [4] H. Chen, S. Succi, S.A. Orszag, Phys. Rev. E 59 (1999) R2527. [5] H. Chen, S.A. Orszag, I. Staroselsky, S. Succi, J. Fluid Mech. 519 (2004) 301. [6] S.S. Girimaji, Phys. Rev. Lett. 99 (2007) 034501. [7] U. Frish, Turbulence: The legancy of A.N. Kolmogorov, Cambridge Univ. Press, 1995. [8] X. Shan, X. Yuan, H. Chen, J. Fluid Mech. 550 (2006) 413. [9] R. Zhang, X. Shan, H. Chen, Phys. Rev. E 74 (2006) 046703. [10] J. Lätt, B. Chopard, Math. Comput. Simul. 72 (2006) 165. [11] J. Lätt, B. Chopard, S. Succi, F. Toschi, Phys. A 362 (2006) 6. [12] P. Burattini, P. Lavoie, A. Agrawal, L. Djenidi, R.A. Antonia, Phys. Rev. E 73 (2006) 066304. [13] X. Shan, X. He, Phys. Rev. Lett. 80 (1998) 65. [14] P. Orlandi, Fluid Flow Phenomena: A Numerical Toolkit, Kluwer Academic, Dordrecht, 1999. [15] R. Mei, L. Luo, P. Lallemand, D. d’Humiéres, Comput. Fluids 35 (2006) 855. [16] H. Yu, S.S. Girimaji, L. Luo, Phys. Rev. E 71 (2005) 016708. [17] S. Pope, Turbulent Flows, Cambridge Univ. Press, 2000. [18] W.D. McComb, The Physics of Fluid Turbulence, Clarendon Press, Oxford, 1990. [19] N.N. Mansour, A.A. Wray, Phys. Fluids 6 (2) (1994) 808. [20] R.A. Antonia, P. Orlandi, J. Fluid Mech. 505 (2004) 123. [21] P. Burattini, P. Lavoie, R.A. Antonia, Exp. Fluids 45 (2008) 523. [22] P.A. Davidson, Turbulence: An Introduction for Scientists and Engineers, Oxford Univ. Press, 2004.