A Lifting Relation from Macroscopic Variables to Mesoscopic Variables in Lattice Boltzmann Method: Derivation, Numerical Assessments and Coupling Computations Validation
AA Lifting Relation from Macroscopic Variablesto Mesoscopic Variables in Lattice BoltzmannMethod: Derivation, Numerical Assessmentsand Coupling Computations Validation
Hui Xu a , b , Huibao Luan a , Yaling He a , Wenquan Tao a , ∗ a Key Laboratory of Thermo-Fluid Science & Engineering, Xi’an JiaotongUniversity, Xi’an, Shaanxi, 710049, China. b Institut Jean le Rond d’Alembert, UMR CNRS 7190, Universit´e Pierre et MarieCurie - Paris 6, 4 Place Jussieu case 162 Tour 55-65, 75252 Paris Cedex 05,France ∗ Corresponding author. Tel:+86-29-82669106. Fax:+86-29-82669106.
Email addresses: [email protected] or [email protected] (Hui Xu), [email protected] (Wenquan Tao ).
Preprint submitted to Elsevier 17 September 2018 a r X i v : . [ phy s i c s . f l u - dyn ] O c t bstract In this paper, analytic relations between the macroscopic variables and themesoscopic variables are derived for lattice Boltzmann methods (LBM). Theanalytic relations are achieved by two different methods for the exchange fromvelocity fields of finite-type methods to the single particle distribution func-tions of LBM. The numerical errors of reconstructing the single particle distri-bution functions and the non-equilibrium distribution function by macroscopicfields are investigated. Results show that their accuracy is better than the ex-isting ones. The proposed reconstruction operator has been used to implementthe coupling computations of LBM and macro-numerical methods of FVM.The lid-driven cavity flow is chosen to carry out the coupling computationsbased on the numerical strategies of domain decomposition methods (DDM).The numerical results show that the proposed lifting relations are accurateand robust.Key words: LBM, Navier-Stokes equations, non-equilibrium distribution func-tions, multi-scale perturbation expansion, coupling computation, FVM2
Introduction
In the past decades, LBM has been widely used to simulate fluid flow problems[1,2], including complex turbulent fluid flows [3,4] and multiscale modeling[5,6].This method is based on the Boltzmann kinetic equation which is used to de-scribe a number of interacting populations of particles. As described in [7],“The LBE could potentially play a twofold function-as a telescope for theatomistic scale and a microscope for the macroscopic scale”. In [8] dense flu-ids flow past and through a carbon nano tube (CNT) was studied by a hybridmodel coupling LBM and MDS. The authors pointed out that replacing thefinite volume solver by a LBM aims to take advantage of the mesoscopicmodeling inherent in LB simulations. Thus LBM is a mesoscopic method innature is a widely-accepted understanding in the literature. The macroscopicparameters such as fluid density, velocity and pressure can be obtained viasome averages of the mesoscopic variable which conform the basic conservationlaws of mass and momentum [2]. In practical applications of LBM to simulatea macroscopic problem, a crucial problem is confronted, that is, a reasonableinitial meso-field must be specified to start the evolution process. The first ini-tializing method was proposed in [9] in 1993. Recently, several methods havebeen proposed to improve the accuracy of numerical results and reduce theinitial layers (oscillation layers) [10,11]. Such oscillations have a numerical ori-gin and are due to the artificial compressibility of LBM. Here, “ initial layer ”refers to such a computational stage within which the macroscopic parametersare oscillating. When the initial data is not well-prepared, there is an initiallayer during which the solution adapts itself to match the profile dictated bythe environment. For the LBM, the existence of the initial layers is a commonphenomenon [10]. In this paper, we will derive the lifting relations betweenthe macroscopic variables and the mesoscopic variables in LBM by two ways.3ccording to the authors’ knowledge, the proposed lifting relations in thispaper are different from those in the existing literature [9,10,11,12,13,14,15].The proposed relations will offer us some new views about the reconstructionof nonequilibrium distribution functions in LBM.Challenging multiscale phenomena or processes are widely existed in materialscience, chemical engineering process, energy and power engineering, and otherengineering fields. Generally speaking, for a multiscale problem, we often mustuse different methods to numerically model the processes at different geomet-ric sub-regions and exchange solution information at interface [16,17,18,19].Such coupling computations are widely adopted in the present-day multiscalesimulation. As indicated above LBM is a kind of mesoscopic methods, whichis a candidate to implement the meso-macro or micro-meso coupling compu-tations in engineering applications [7]. So, the proposed method not only canbe used to obtain a better initial field for LBM, but also can be adopted inthe multi-scale computation. For example in [7] the possibility of couplingLBM with molecular dynamics simulation (MDS) was investigated and foundthat with proper time and geometric scales the two numerical methods canbe coupled. And in [8] such coupling simulation was conducted. In the ex-isting literatures the coupling of finite difference method (FDM, which is amacrosopic method) with LBM was adopted in [19,20,21], but the proposedcoupling method is similar to a multigrid method and a simple regulariza-tion formula is used in their computations. The regularization formula in [19]only considers the first-order approximation of the single particle distribu-tion function and the coupling formula in [20] is only used to deal with theone-dimensional reaction-diffusion system. In [8] the coupling between LBMand MDS was implemented by exchange of velocity and velocity gradient atthe interface region. In this paper, the proposed meso-macro (or micro-meso)4oupling is expected to be used for domain decomposition methods, in whichLBM and macro-type numerical method (or micro-type numerical method andLBM) are adopted in different sub-domain and information is exchanged at theinterface. We believe that our proposed relation is more useful method for en-gineering multiscale computations. In addition, the proposed coupling methodcan also be used to carry out the multigrid computations and equation-freemultiscale (EFM) computations [22]. It is well-known that LBM is very power-full for the parallel computing on a low cost [23,24]. So, the proposed relationcan be used in the parallel simulations for multiscale simulations of complexfluid flows based on the refinement strategies.To the authors’ understanding the glossary “lifting relation” means that macro-scopic variables in a lower degree-of-freedom (DoF) system are upscaled tomeso/microscopic variables in a higher DoF system. Generally, it is difficultto establish the one-to-one map from a lower DoF system to a higher DoFsystem, although the lower DoF system can be seemed to be an approximateor approaching form of a higher DoF system in some referred scales. This sit-uation happens when numerical results of different scales are coupled at thesame location. For example when MDS and continuum method are coupled,reference [25] indicated that it is straightforward to obtain the continuumquantities (such as velocity, pressure) from the particle description by averag-ing over the local region and over time, but the reverse problem, generatingmeso/microscopic particle configuration from known macroscopic quantities isnon-trivial and must necessarily be non-unique. The glossary “lifting relation”in the title of this paper is proposed based on the concept of the DoF of thegoverning equations.In this paper, we will give two methods to establish the relations betweenvariables of the Navier-Stokes equations and variables of LBM. Numerical5ests demonstrate that the proposed methods of computing non-equilibriumdistribution functions are effective and accurate.The rest of the paper is organized as follow. In section 2, the details of multi-scale derivation of non-equilibrium distribution functions is given. In section 3,the non-equilibrium distribution functions are obtained by Boltzmann-BGKequations. In section 4, the performances of the proposed relations to recon-struct non-equilibrium distribution functions are demonstrated by numericaltests. Finally, some conclusions are given.
In this section, we will review LBM and the corresponding macroscopic equa-tion. Based on this review, we will derive a relation for lifting macroscopicvariables to microscopic variables by multiscale approach.
We now introduce the lattice Boltzmann-BGK model as a solver for theweakly-compressible Navier-Stokes equations. LBM is built up from the latticegas cellular automata models [2]. The numerical scheme of LBM is establishedbased on a finite discrete-velocity model of the Boltzmann-BGK equation andcan be expressed as follows f i (x + δt c i , t + δt ) − f (x , t ) = Ω i , (1)where f i represents the single-particle distribution function along the directionc i ( i = 0 , . . . , n ), c i is the element of the discrete velocity set V = { c , . . . , c n } .Ω i denotes the collision operator which is non-dimensional. The macroscopic6ariables, the density ρ and the velocity u, are defined locally by the distribu-tion functions as follows ρ (x , t ) = n (cid:88) i =0 f i (x , t ) = n (cid:88) i =0 f eq i (x , t ) , (2)u(x , t ) = 1 ρ (cid:88) c i ∈V c i f i (x , t ) = 1 ρ (cid:88) c i ∈V c i f (eq) i (x , t ) . (3)For the standard LBM, the collision operator is defined by the so-called BGKcollision Ω BGK i = − τ lbm [ f i (x , t ) − f (eq) i (x , t )] . (4)For the convenience of comparison, from here, we use the similar notations in[26]. The local equilibrium distribution f (eq) i is defined by f (eq) i (x , t ) = f L ( eq ) i (x , t ) + f Q(eq) i (x , t ) , (5)where f L(eq) i (x , t ) and f Q(eq) i (x , t ) denote the linear part and the quadratic partof the equilibrium distribution, respectively. The linear part is given by f L(eq) i (x , t ) = ω i ρ (1 + 1 c s c i · u(x , t )) , (6)and the quadratic part is expressed by f Q(eq) i (x , t ) = ω i c s ρ (u(x , t )u(x , t )) : Σ i , (7)where c s is the lattice sound speed of the model, ω i denotes the weight andΣ i is a second-order tensor defined byΣ iαβ = c iα c iβ − c s δ αβ . (8)The tensor product definition between two first order tensors a and b is givenas follows ( ab ) αβ = a α b β , (9)7nd the corresponding second-order tensor :-product between A and B is givenby A : B = d (cid:88) α,β =1 A αβ B αβ , (10)where d denotes the spatial dimension.In this paper, we mainly focus on the standard LBM. By the Chapman-Enskogexpansion, under the small M a number restriction (
M a ≤ . ∂ t ρ + ∂ α ( ρu α ) + O ( δt ) = 0 , (11) ∂ t ( ρu α ) + ∂ β ( ρu α u β ) = − ∂ α p + ν∂ β ( ρ ( ∂ α u β + ∂ β u α )) + O ( δt ) + O ( δtu ) , (12)where p is defined by p = c s ρ. It is clear that the recovered Navier-Stokes equations are weakly compressible[2,27,28]. So, the density is coupled with the pressure field in LBM. In Eq.(12), the second term of R.H.S can be rewritten as ν∂ β ( ρ ( ∂ α u β + ∂ β u α )) = νρ ( ∂ β ∂ β u α ) + ν ( ∂ β ρ )( ∂ α u β + ∂ β u α ) + νρ∂ α ∂ β u β . (13)And the corresponding third-order term O ( δtu ) is given by O ( δtu ) = − σ∂ β ∂ γ ( ρu α u β u γ ) . (14)The fluid viscosity ν is defined by ν = c s ( τ lbm −
12 ) δt, (15)and σ is given by σ = νc s . (16)In Eq. (13), the third term of R.H.S will vanish for a divergence-free field.But the second term will not vanish, if the density ρ is nonhomogeneous in8he spatial domain. The Navier-Stokes equations are recovered by LBM underthe low Mach condition. Physically, LBM is a weakly compressible model forsolving Navier-Stokes equations.At this point, we describe two situations where the lifting relation is useful. Thefirst situation is using the lifting relation to get a good initial field of the densitydistribution function from specified velocity and pressure fields. As indicatedabove the recovered Navier-Stokes equations are weakly compressible, hencepressure field is coupled with the density field by the equation of state ( p = c s ρ ). In engineering computations, the weakly-compressible flow is often usedas an approximation of the incompressible flow. For the lifting function, theconsideration should be made from the weakly compressible side. The non-homogeneous character of the initial density is very significant for an initialroutine of LBM in the proposed lift relation. This significance can be observedfrom the follow-up derivations. For the initial processes, if the initial pressurefield is given, the lifting relation can be used to obtain the initial distributionfunctions consistent with the recovered Navier-Stokes equations. In anotherdevelopment when we couple LBM with other macroscopic solver of Navier-Stokes equations, we need to pass the macroscopic variables (pressure andvelocity fields) to an approximate single particle distribution functions or thenon-equilibrium distribution functions. At this time, a macroscopic equationrelating to the given velocity and pressure to the particle distribution functionof LBM become very useful. The major goal of the present paper is to derivesuch a lift relation, or a reconstruction operator as depicted in [9].For the convenience of deriving such an equation, some changes are made forthe form of Eq. (12). We first rewrite Eq. (12) as ∂ t ( ρu α ) + ∂ β ( ρu α u β ) = ρ ( ∂ t u α + u β ∂ β u α ) + u α ( ∂ t ρ + ∂ β ( ρu β )) . (17)9f the initial velocity field is divergence-free, we have ∂ t ( ρu α ) + ∂ β ( ρu α u β ) = ρ ( ∂ t u α + ∂ β ( u α u β )) + u α ( ∂ t ρ + ∂ β ( ρu β )) . (18)The neglecting of the term ρu α ∂ β u β is a widely accepted approximation. Ac-cording to Eq. (11), we have ∂ t ( ρu α ) + ∂ β ( ρu α u β ) = ρ ( ∂ t u α + ∂ β ( u α u β )) . (19)Now, combining Eq. (12), Eq. (13) and Eq. (17), we gain ∂ t u α + u β ∂ β u α = − ∂ α pρ + ν ( ∂ β ∂ β u α + ∂ α ∂ β u β ) + ν ∂ β ρρ ( ∂ α u β + ∂ β u α ) . (20) The coupled macro-micro/mesoscale simulation is a rapidly developing areaof research that deals with processes covering several order of geometries. Forsuch numerical approach, one needs to construct an initial condition u ( x,
0) forthe meso/microscopic simulator, which is corresponding to the initial macro-scopic variable U ( x, u ( x,
0) represents the meso/microscopic statevariables and U ( x,
0) stands for macroscopic state variables. As indicatedabove this procedure is called lif ting [22] or reconstruction [29] step. Thelifting (reconstruction) operator µ is defined by u ( x,
0) = µ ( U ( x, . (21)The lifting procedure leads to a one-to-many mapping. After the initializationof the meso/microscopic variables by the reconstruction operator µ , they willbe evolved by the meso/microscopic simulator. In this paper, LBM is adoptedas the mesoscopic simulator. As indicated in [18,20] the macroscopic state10ariables are easy to be achieved. To transfer the micro/meso-scale parametersinto macro parameters we need some restriction [22] or compression [29]operators. Conceptually, this operator M is defined by U ( x, t ) = M ( u ( x, t )) . (22)For LBM, the operator M is implemented by Eq. (2) and Eq. (3). Our at-tention will put on the development of the reconstruction operator µ by themulti-scale analysis. As discussed above the reconstruction operator in multi-scale computation is corresponding to the lifting relation in an initial problem.In the following we will derive the operator from the initial problem aspect.To obtain an appropriate initial field, we turn to a simple multiscale pertur-bation expansion. We separate the time scale into two different time scales, t = (cid:15)t (diffusive time-scale) and t = (cid:15) t (convective time-scale). The timederivative ∂ t is expanded using a small parameter (cid:15) , which normally is pro-portional to the small Knudsen number ( Kn < .
1) [27], ∂ t = (cid:15)∂ t + (cid:15) ∂ t + O ( (cid:15) ) . (23)Similarly, introducing space scale x = (cid:15)x , the corresponding spatial derivativeis not expanded beyond the first-order term [27] ∂ α = (cid:15)∂ α + O ( (cid:15) ) . (24)The single-particle distribution function is expanded as follows [27] f i (x , t ) = f (0) i (x , t ) + (cid:15)f (1) i (x , t ) + (cid:15) f (2) i (x , t ) + . . . . (25)By the Taylor expansion, from Eq. (1), we get δt ( ∂ t + c iα ∂ α ) f i (x , t ) + δt ( ∂ t + c iα ∂ α ) f i (x , t ) + O ( δt ) = Ω i . (26)11ombining Eq.(23)-Eq.(25) with Eq.(26), we obtain f (0) i (x , t ) = f (eq) i (x , t ) (27)and (cid:15)f (1) i (x , t ) + (cid:15) f (2) i (x , t ) = − τ lbm [( (cid:15)∂ t + (cid:15) ∂ t + (cid:15)c iα ∂ α ) δt +( (cid:15) ∂ t c iα ∂ α + (cid:15) ∂ t c iα ∂ α + (cid:15) c iα c iβ ∂ α ∂ β ) δt ]( f (0) i (x , t ) + (cid:15)f (1) i (x , t )) + O ( δt ) . (28)For first order of (cid:15) , we get f (1) i (x , t ) = − τ lbm δt ( ∂ t + c iα ∂ α (x , t )) f (eq) i + O ( δt ) . (29)According to Eq. (2)-Eq. (3), we have following equations in the first-orderscale of (cid:15) [30] ∂ t ρ + ∂ α ( ρu α ) + O ( δt ) = 0 , (30) ∂ t ( ρu α ) + ∂ β ( ρu α u β + c s ρδ αβ ) + O ( δt ) = 0 , (31)Then, Eq. (31) can be rewritten as ρ∂ t ( u α ) + ρu β ∂ β ( ρu α + c s ρδ αβ ) + O ( δt ) = 0 . (32)By matching small scales, from Eq. (28), we can get up to the second orderequations of the small parameter (cid:15) : f (2) i = − τ lbm δt∂ t f (0) i − δt ( τ lbm −
12 )( ∂ t + c iβ ∂ β ) f (0) i + O ( δt ) . (33)Then, we can get [30] ∂ t ρ + O ( δt ) = 0 , (34) ∂ t ( ρu α ) = ν∂ β ( ρ ( ∂ α u β + ∂ β u α )) + O ( δt + δtu ) (35)Furthermore, from Eq. (29), we have f (1) i (x , t ) = − τ lbm δt ( ∂ t + c iα ∂ α )( f L(eq) i (x , t ) + f Q(eq) i (x , t )) + O ( δt ) . (36)12n the derivation of Eq. (36), we introduce the following formulas accordingto the chain rule of derivatives [31] ∂ t f (eq) i (x , t ) = ∂ ρ f (eq) i (x , t ) ∂ t ρ + ∂ u β f ( eq ) i (x , t ) ∂ t u β , (37) ∂ α f (eq) i (x , t ) = ∂ ρ f (eq) i (x , t ) ∂ α ρ + ∂ u β f (eq) i (x , t ) ∂ α u β . (38)Now, the equilibrium function can be differentiated by the macroscopic vari-ables as follows [31] ∂ ρ f (eq) i (x , t ) = 1 ρ f (eq) i (x , t ) , (39) ∂ u β f (eq) i (x , t ) = ∂ u β f L(eq) i (x , t ) + ∂ u β f Q(eq) i (x , t ) . (40)According to Eq. (6) and Eq. (7), we have ∂ u β f L(eq) i = ω i ρ c s c iβ , (41) ∂ u β f Q(eq) i = ω i ρ c s (2 c iα c iβ u α − c s u β ) = ω i ρ ( c s c iα c iβ u α − c s u β ) . (42)So, we have ∂ u β f (eq) i = ω i ρ [ 1 c s ( c iβ − u β ) + 1 c s c iα c iβ u α ] . (43)Come here we can have following corollaries. Corollary 1
From Eq. (36), for the first-order approximation of (cid:15) , there ex-ists a lifting relation from the macroscopic variables to the microscopic variable f (1) i f (1) i = − τ lbm δt { ( c iα − u α ) ∂ α ρ∂ ρ f eqi + ( c iα − u α ) ∂ α u β ∂ u β f (eq) i − ρ∂ ρ f (eq) i ∂ α u α − ρ ∂ α p∂ u α f (eq) i } = − τ lbm δt { ( c iα − u α ) ρ ∂ α ρf (eq) i + ( c iα − u α ) ∂ α u β ω i ρ [ c s ( c iβ − u β ) + c s c iβ c iγ u γ ] − f (eq) i ∂ α u α − ρ ∂ β p ω i ρ [ c s ( c iβ − u β ) + c s c iβ c iγ u γ ] } . (44)13 orollary 2 From Eq. (33), for the second-order scale of (cid:15) , we have the fol-lowing approximation f (2) i ≈ − τ lbm δt∂ t f (eq) i , (45) where the second-order derivative of f (0) i is ignored. Hence, we can easily establish an approximation for f (2) i by the method anal-ogous to the approximation of f (1) i as follows ∂ t f (eq) i (x , t ) = ∂ ρ f (eq) i (x , t ) ∂ t ρ + ∂ u β f ( eq ) i (x , t ) ∂ t u β . (46)By Eq. (34), we have ∂ t f (eq) i (x , t ) = ∂ u β f (eq) i (x , t ) ∂ t u β = 1 ρ ∂ u β f (eq) i (x , t ) ∂ t ( ρu β ) . (47)From Eq. (35) and Eq. (43), it is easy to obtain ∂ t f (eq) i = νω i [ 1 c s ( c iβ − u β ) + 1 c s c iβ c iγ u γ ] ∂ α ( ρ ( ∂ β u α + ∂ α u β )) . (48)So, we have (cid:15) f (2) i ≈ − τ δtνω i [ 1 c s ( c iβ − u β ) + 1 c s c iβ c iγ u γ ] ∂ α ( ρ ( ∂ β u α + ∂ α u β )) . (49)By a simple derivation, we have ∂ α ( ρ ( ∂ β u α + ∂ α u β )) = ∂ α ρ ( ∂ β u α + ∂ α u β ) + ρ ( ∂ β ∂ α u α + ∂ α u β ) . (50)From Eqs. (49) ∼ (50), we have (cid:15) f (2) i ≈ − τ lbm δtνω i [ c s ( c iβ − u β ) + c s c iβ c iγ u γ ]( ∂ α ρ ( ∂ β u α + ∂ α u β )+ ρ ( ∂ β ∂ α u α + ∂ α u β )) (51)Therefore, we get the following approximation of the non-equilibrium distri-bution function from Eq. (25) f (neq) i ≈ (cid:15)f (1) i + (cid:15) f (2) i , (52)14hat is, f (neq) i (x , t ) ≈ − τ lbm δt { u T,iα ρ ∂ α ρf (eq) i + u T,iα ∂ α u β ω i ρ [ c s u T,iβ + c s c iβ c iγ u γ ] − f (eq) i ∂ α u α − ρ ∂ β p ω i ρ [ c s u T,iβ + c s c iβ c iγ u γ ]+ ω i ρ [ c s u T,iβ + c s c iβ c iγ u γ ]( ν ρ ∂ α ρ ( ∂ β u α + ∂ α u β ) + ν ( ∂ β ∂ α u α + ∂ α u β )) } , (53)where u T,iα = c iα − u α (u T = c i − u, peculiar velocity). Since the velocity fieldis divergence-free, we have f (neq dfv) i ≈ − τ lbm δt { u T,iα ρ ∂ α ρf (eq) i + u T,iα ∂ α u β ω i ρ [ c s u T,iβ + c s c iβ c iγ u γ ] −− ρ ∂ β p ω i ρ [ c s u T,iβ + c s c iβ c iγ u γ ]+ ω i ρ [ c s u T,iβ + c s c iβ c iγ u γ ]( ν ρ ∂ α ρ ( ∂ β u α + ∂ α u β ) + ν∂ α u β ) } . (54)Here, we also introduce an approximation of ∂ u β f (eq) i by ignoring the higher-order terms of u as adopted in [31] ∂ u β f (eq) i = ω i ρ [ 1 c s u T,iβ + 1 c s c iβ c iγ u γ ] ≈ u T,iβ c s f (eq) i . (55)Now, we have f (neq dfv) i ≈ − τ lbm δt c s f (eq) i u T,iβ ( u T,iα ∂ α u β + ν ρ ∂ α ρ ( ∂ β u α + ∂ α u β ) + ν∂ α u β ) . (56)Rewriting the above formula, we obtain f (neq dfv) i ≈ − τ lbm δtf (eq) i c s u T,iβ ( u T,iα ∂ α u β + ν∂ α u β + ν ρ ∂ α ρS αβ ) , (57)where S αβ = ∂ β u α + ∂ α u β .In all, we can get an approximation of the single-particle distribution function15or divergence-free velocity fields as follows f i ≈ f (eq) i { − c s τ lbm δtu T,iβ ( u T,iα ∂ α u β + ν∂ α u β + ν ρ ∂ α ρS αβ ) } , (58)By a similar deviation, we can get an approximation of the single-particledistribution function for weak-compressible velocity fields as follows: f i ≈ f (eq) i { − c s τ lbm δt [ u T,iβ ( u T,iα ∂ α u β + ν ( ∂ α u β + ∂ β ∂ α u α )+ ν ρ ∂ α ρS αβ ) − c s ∂ α u α ] } , (59)Now we compare our results with that published in literatures.1. T. Imamura et al [31] obtained the following formula f i ≈ f (eq) i + (cid:15)f (1) i = f (eq) i [1 − τ lbm δt ( 3 u T,iα u T,iβ c − δ αβ ) ∂ β u α ] (60)They only used f (1) i to approximate the single-partial distribution functions. Itis well-known that in order to recover the correct Navier-Stokes equations, f (2) i is needed. From this point of view, the approaching form (59) of the distribu-tion functions are more accurate than (60). If the divergence-free velocity fieldis considered, Eq. (58) is also superior to Eq. (60) because Eq. (58) containsthe information of f (2) i which is related with molecule viscosity and densitygradient. As for the lifting relation it is certainly essential to involve moleculeviscosity and density gradient [19,20].2. Skodors [9] gave the following formula (ignoring the term of O ( M a )) f (neq) , S i = − τ lbm δ t ω i (cid:34) c s c i c i : ∇ ( ρ u) − ∇ · ( ρ u) (cid:35) . (61)Guo and Zhao [32] further simplified Eq.(61) and obtained the following rela-tion f (neq) , G i = − τ lbm δ t ω i ρ c s c i c i : ∇ u . (62)16t is very clear that Eqs.(58) and (59) are totally different from Eqs.(61) and(62), respectively. The co-existence of Eqs. (58) ∼ (59) and Eqs.(61) ∼ (62) asthe lift correlation for the same situation may be regarded as the witness thatthe transformation from one-to-many must necessarily be not unique [25].Some comparisons will be performed in Sec. 4 between Eqs. (58) ∼ (59) andEqs.(61) ∼
62) for schemes of D2Q9 and D2Q17. It turns out that the accuracyof Eqs. (58) ∼ (59) derived in this paper is better than that of Eqs.(61) ∼ (62).The derivation procedures of Eqs. (58) ∼ (59) kept the information of the f (2)2 and other more details which are important to reduce the reconstruction rel-ative errors. The Boltzmann equation [33] describes the statistical distribution of particlesin a fluid. It is one of the most important equations of non-equilibrium statisti-cal mechanics, which deals with systems far from thermodynamic equilibrium[34]. The Boltzmann equation is described by ∂f (x , v , t ) ∂t + v · ∇ x f (x , v , t ) + 1 m F ( x ) ∇ v f (x , v , t ) = Ω( f (x , v , t )) . (63)The Boltzmann equation (63) is an equation for the time t evolution of thedistribution (properly a density) function f (x , v , t ) in one-particle phase space,where x = ( x , x , · · · , x d ) ∈ R d and v = ( v , v , · · · , v d ) ∈ R d (d denotesthe spatial dimension) are position and velocity, respectively. The equilibriumdistribution function f eq (x , v , t ) can be determined by f (eq) (x , v , t ) = n (x , t ) (cid:32) m πκT (x , t ) (cid:33) d/ exp (cid:34) − m κT (x , t ) (v − u(x , t )) (cid:35) (64)17ere, the quantities T (x , t ), n (x , t ) and u(x , t ) represent the local tempera-ture , the local particle-number distribution density and the local velocity [2,34],repectively. u T = v − u(x , t ) is the so called thermal velocity. m denotesthe single-particle mass which is set to be unity for convenience. In orderto simplify the complex collisional term, the following conserved relaxationtime approximation is used to describe the collision term through only onecharacteristic frequency[34] ∂f (x , v , t ) ∂t + v · ∇ f (x , v , t ) = − τ ( f (x , v , t ) − f (eq) (x , v , t )) , (65)where the external force term is not considered and ∇ denotes ∇ x . τ representsthe relaxation time.In order to solve Eq.(65), the velocity space is discretized [2] and we gain ∂f i (x , t ) ∂t + c i · ∇ f i (x , t ) = − τ ( f i (x , t ) − f (eq) i (x , t )) , (66)where w i denotes the integral weight factor, f i (x , t ) = w i f (x , c i , t ) and f (eq) i (x , t ) = w i f (eq) (x , c i , t ). Furthermore, along the characteristic line, the time-discretizationform of Eq.(66) can be expressed as [2,36] f i (x + c i δt, t + δt ) = f i (x , t ) − τ lbm ( f i (x , t ) − f (eq) i (x , t )) , i = 0 , , ..., N. (67)where f i is the probability distribution function (PDF) along the ith direction, f (eq) i is its corresponding equilibrium PDF, δt is the time step, c i is the particlevelocity in the ith direction, and N is the number of the discrete particlevelocities. Note: τ lbm = τ /δt which is a dimensionless relaxation time. Thelocal macro quantities are defined by Eqs. (2) and (3).At the low fluid flow velocity (or low Mach number), an approximate form ofthe equilibrium distribution function f (eq) i is described by the discrete equilib-rium distribution, Eqs. (5) ∼ (7). 18ow, we consider the conserved relaxation time approximation of the Boltz-mann equation (65). The right hand side of Eq. (65) represents the differencebetween the distribution function and a local Maxwell distribution. This dif-ference is termed non-equilibrium distribution defined by f (neq) (x , v , t ) = f (x , v , t ) − f (eq) (x , v , t ) . (68)Then, Eq. (65) can be rewritten as follows (cid:32) ∂∂t + v · ∇ (cid:33) f (neq) (x , v , t ) + (cid:32) ∂∂t + v · ∇ (cid:33) f (eq) (x , v , t ) = − τ f (neq) (x , v , t ) . (69)In the hydrodynamic region [33], the first term on the left-hand side of Eq.(69) can be neglected compared with the right-hand side [34]. Then, we obtain (cid:32) ∂∂t + v · ∇ (cid:33) f (eq) (x , v , t ) = − τ f (neq) (x , v , t ) . (70)In terms of the Maxwell equilibrium distribution and assuming a uniformtemperature of the system, we can obtain f (eq) (x , v ,t ) n (x ,t ) (cid:16) ∂∂t + v · ∇ (cid:17) n (x , t ) − f (eq) (x , v , t ) (cid:16) ∂∂t + v · ∇ (cid:17) u κT = − τ f (neq) (x , v , t ) , (71)where T = T (x , t ) = constant. In Eq. (71), the left-hand term can be rewrittenas follows f (eq) (x , v ,t ) n (x ,t ) (cid:16) ∂∂t + v · ∇ (cid:17) n (x , t ) = f (eq) (x , v ,t ) n (x ,t ) (cid:16) ∂∂t + u(x , t) · ∇ (cid:17) n (x , t )+ f (eq) (x , v ,t ) n (x ,t ) u T · ∇ n (x , t ) (72)In order to satisfy the mass conservation condition of the fluid flow system,the first term of the right-hand side in Eq. (72) should be equal to zero. Hence,19e have the following equation f (eq) (x , v , t ) (cid:16) ∂∂t + v · ∇ (cid:17) u c s − f (eq) (x , v ,t ) n (x ,t ) (u T · ∇ n (x , t ) − n (x , t ) ∇ · u(x , t )) = τ f (neq) (x , v , t ) , (73)where c s = √ κT . The term u is the thermal fluctuation energy, thus thenon-equilibrium is determined by the material derivative of this thermal fluc-tuation energy. The quantity (cid:16) ∂∂t + v · ∇ (cid:17) u can be determined by the dynam-ical equation corresponding to the micro dynamical system. Here, we rewrite (cid:16) ∂∂t + v · ∇ (cid:17) u as follows12 (cid:32) ∂∂t + v · ∇ (cid:33) u = − u T · (cid:32) ∂∂t + u(x , t ) · ∇ (cid:33) u(x , t ) − u T · (u T ·∇ )u(x , t ) . (74)Generally, the governing equation of the macroscopic physical quantity is rep-resented by DDt u(x , t ) = (cid:32) ∂∂t + u(x , t ) · ∇ (cid:33) u(x , t ) = F (x , u(x , t ) , t ) . (75)Normally, the macroscopic physical quantity u(x , t ) in the governing equationis known. So, F (x , u(x , t ) , t ) can be determined easily. For fluid flow problems,taking u(x , t ) as fluid velocity, then F (x , u(x , t ) , t ) can be estimated by fluidacceleration. The term u T · (u T · ∇ )u in Eq. (74) can be determined by u(x , t )and the spatial derivatives of u(x , t ).The lattice Boltzmann model is a special discrete form of the BGK latticeBolzmann equation with respect to temporal and spatial variables. For LBMthe equilibrium distribution, Eq. (5), is a polynomial-truncated approximationof the Maxwell distribution up to O ( | u | ), so Eq. (73) can be applied to LBM20irectly as follows f (eq) i (x , t ) (cid:16) ∂∂t + c i · ∇ (cid:17) u i, T c s − f (eq) i (x , v ,t ) n (x ,t ) (u i, T · ∇ n (x , t ) − n (x , t ) ∇ · u(x , t )) = τ f (neq) i (x , t ) , (76)where u i, T = c i − u(x , t ). Now, the non-equilibrium distribution function canbe denoted by f (neq) i (x , t ) = − τf (eq) i (x ,t ) c s (cid:104) u i, T · (cid:16) F (x , u(x , t ) , t ) + (u i, T · ∇ )u + c s n (x ,t ) ∇ n (x , t ) (cid:17) − c s ∇ · u(x , t ) (cid:105) . (77)The derivation of Eq.(77) is completed based on the rigorous inherent physicalconsistency in the hydrodynamic region and and the derivation is independenton the spatial dimension. Meanwhile, the Maxwell equilibrium distribution isregarded as the tool to implement the analysis.It is worth pointing out that for DnQb LBM, F (x , u(x , t ) , t ) can easily bedetermined from the recovered Naiver-Stokes equations, so the obtained non-equilibrium distribution function formulas (77) and (59) are identical. Thus,by using different derivation method we come to the same conclusion.In addition, according to Eqs. (58),(59) and (77), it can be seen that thenon-equilibrium distribution functions have the following form f (neq) i = f (eq) i λ i ( ρ, u) , (78)where λ i ( ρ, u) is a perturbative parameter with respect to ρ and u. The pa-rameter λ i ( ρ, u) in Eq. (78) needs to satisfy the following constraints (cid:88) i f (eq) i λ i ( ρ, u) = 0 , (cid:88) c i ∈V c i f (eq) i λ i ( ρ, u) = 0 (79)21 Numerical Tests
In this section, the non-equilibrium distribution function will be validated bynumerical tests. The numerical tests focus on validating the precision of thereconstruction operator and the correctness of the coupling computations. It’sworth noting that the word “multiscale simulation” used in this paper is re-ferred to the coupling between numerical methods of microscale (moleculardynamics simulation), mesoscale (LBM) and macroscale(say, FVM) adoptedin neighboring computational regions. And for such coupling the major con-cern is the transformation of solutions from macro(or meso)scales to meso(ormicro)scales at the interface. The focus of the following presentation is to val-idate the correctness of the proposed operators. Because of space limitationthe effect of the grid fineness on the numerical solution will not be conducted.Reference [35] can be referred. The effect of the mesh size on the accuracy ofthe reconstruction operator will be presented in Sec. 4.2.
In order to validate Formula (77), the D2Q9 [36] and D2Q17 [37] LBM areadopted to simulate 2D fluid flows. At low Mach number (
M a = u(x , t ) / c s (cid:28) F α (x , u(x , t ) , t ) = − ∂ α pρ + ν ( ∂ β ∂ β u α + ∂ α ∂ β u β ) + ν ∂ β ρρ ( ∂ α u β + ∂ β u α )(80)where ν = c s ( τ lbm −
12 ) δt. (81)The details of the macroscopic dynamic equation corresponding to D2Q17LBM are omitted (see [37]). Now, the non-equilibrium distribution in Eq.(77)22an be determined directly by the right-hand side of Eq.(20). For any giveninitial velocity and density fields, each term in the right-hand side of Eq.(20)can be calculated. In order to validate the precision of the proposed method,the following two basic quantities are defined (cid:98) f i (x , t ) = f (eq) i (x , t ) + (cid:98) f (neq) i (x , t ) , (82) f (neq) i (x , t ) = f i (x , t ) − f (eq) i (x , t ) (83)where (cid:98) f (neq) i (x , t ) is called reconstructed non-equilibrium distribution functionand is calculated by Eq.(77) and (cid:98) f i (x , t ) is the reconstructed single-particledistribution function. f (neq) i (x , t ) and f i (x , t ) denote the real non-equilibriumdistribution function and the real single-particle distribution function, respec-tively. Here, we give two kinds of relative error definitions: single particledistribution function reconstruction error, single particle non-equilibrium dis-tribution function reconstruction errorE( f i , (cid:98) f i ) = (cid:118)(cid:117)(cid:117)(cid:116) N um × ( n + 1) (cid:88) x (cid:88) i | (cid:98) f i (x , t ) − f i (x , t ) | f i (x , t ) , (84)E( f (neq) i , (cid:98) f (neq) i ) = (cid:118)(cid:117)(cid:117)(cid:116) N um × ( n + 1) (cid:88) x (cid:88) i | (cid:98) f (neq) i (x , t ) − f (neq) i (x , t ) | f (neq) i (x , t ) (85)where N um denotes the number of lattice nodes.In order to demonstrate the proposed method, a freely-decaying 2D turbulenceproblem will be simulated by the proposed method. This turbulence problemoften makes the local discrete single-particle distribution functions to be farfrom the local discrete equilibrium distribution functions, which yields a richvelocity structure. The freely-decaying 2D turbulence is implemented in aperiodic box Ω = [0 , π ] × [0 , π ]. A 2D random velocity field will be specifiedas the initial condition. The initial fields are initialized by random phase in23ourier spectral space and the initial spectrum is given by [38] E ( k,
0) = a s u k − p (cid:32) kk p (cid:33) (2 s +1) exp − (cid:18) s + 12 (cid:19) (cid:32) kk p (cid:33) (86)where s = 0 , , , · · · , and the normalization constant a s is given by a s = (2 s + 1) s +1 / s s ! . All the results presented below correspond to s = 3, k p = 16, u = { . , . } and ρ = 2 .
7. The lattice size is 512 × L is equalto 0.12953. The Reynolds number ( Re L = Lu /ν ) is equal to 111 . u = 0 . t = 1000 δt , it is clear that the recon-structed single-particle distribution functions and the non-equilibrium distri-bution functions coincide with the real single-particle distribution functionsand non-equilibrium distribution functions very well for D2Q9 and D2Q17 inFigs (1)-(2). The corresponding relative errors E( f i , (cid:98) f i ) are about 0 . . f (neq) i , (cid:98) f (neq) i ) are about 16 . . et al [31] is used tocalculate the single-particle non-equilibrium distribution functions, the rel-ative errors E( f (neq) i , (cid:98) f (neq) i ) are up to about 21 .
65% and 18 .
13% for D2Q9and D2Q17, respectively. We also adopted Eqs. (61) in [9] and (62)in [32]to do the same calculations. The relative errors E( f (neq) i , (cid:98) f (neq) i ) of the single-particle non-equilibrium distribution functions can be up to about 80% atmany lattice nodes. In Fig. 5, the numerical relation between f (neq) i and (cid:98) f (neq) i for the method in [32]. The mean relative error E( f (neq) i , (cid:98) f (neq) i ) is larger than43 .
74% for D2Q9. In the statistical procedure, we ignore the points with very24mall f (neq) i and (cid:98) f (neq) i ( f (neq) i , (cid:98) f (neq) i < − ) for the method in [32]. Here, wemust point out that when f (neq) i and (cid:98) f (neq) i are very small, the relative errorsE( f (neq) i , (cid:98) f (neq) i ) of the methods in [9,32] are very large. In such a circumstance,the relative error of the non-equilibrium distribution functions by Eq. (77) isalso a bit larger, but it still less than that computed by Eq. (60) [31] andmuch less than that computed by Eqs. (61) ∼ (62) of [9] and [32], respectively.Similar results can be observed for the case of u = 0 .
01 at t = 10000 δt forD2Q9 and D2Q17. For the simplicity of presentation, they are not providedhere.In addition we also found that when the single-particle distribution functionsand non-equilibrium distribution functions are reconstructed, the results fromD2Q17 model show a better accuracy than that of D2Q9 model. Meanwhile,from the both models, more accurate results can be gained when the Machnumber is reduced. Such results are very reasonable, and can be understoodas follows. First, D2Q17 model is more accurate to approach Maxwell dis-tribution function in discrete velocity spaces than D2Q9 model. Second, lowMach number will lead to a reduction of the truncated errors for approach-ing Maxwell distributions and a better recovering Navier-Stokes equation. It isproved [37] that D2Q17 model can eliminate the third-order term of statisticalvelocity in recovered Navier-Stokes equation.Finally, attention is turned to the comparison of vorticity by the real f i (x , t )and the reconstructed (cid:98) f i (x , t ) in Figs. 6 ∼
7, where the vorticity contour fig-ures are given for u = 0 . u = 0 .
01, respectively.In order to show thequantitative sense of the vorticity reconstruction error, we choose 100 and1000 time-series samples for u = 0 . u = 0 . L -relative departures of the reconstructed vorticity are 0 . ± . u = 0 . . ± . u = 0 . . ± . = 0 .
01) and 0 . ± . u = 0 . In order to validate the approach behaviors versus different grid sizes, we givethe convergence properties of the D2Q9 and D2Q17 models by different meshscales. The 2D Taylor-Green vortex problem is chosen as the intial fields u = − A cos( k x )sin( k y ) F ( t ) v = A k k sin( k x )cos( k y ) F ( t ) p = p − A (cid:34) cos(2 k x ) + k k cos(2 k y ) (cid:35) F ( t ) (87)where F ( t ) = exp [ − ν ( k + k ) t ], A = 0 . k = k = 4 and p = ρ c s . Thecomputational domain Ω = [0 , π ] and Re = 10000. The periodic boundaryconditions are applied in both directions. The initial distribution functionsare initialized by the reconstruction operator. The reconstruction L and L relative errors of the distribution functions are calculated at the time steps n = { , , , , } corresponding to the mesh resolutions h = { / , / , / , / , / } respectively. In Fig. 8, the relative errors aregiven in the log-log coordinates. From the results, it is clear that for theD2Q9 model and the D2Q17 model, they nearly have the same convergence26ates which are approximately equal to 2.6. However, the relative errors ofthe D2Q17 model are smaller than that of the D2Q9 model. That means thereconstruction precision can be improved when the number of the discretevelocity increases. This conclusion is consistent with the result in Sec. 4.1. In order to illustrate the feasibility of the recommended reconstruction opera-tor, the lid-driven cavity flow is simulated by the coupled LBM-FVM method.The computational domain is decomposed in two regions in which the LBMand FVM methods are used respectively (see Fig. 9-(a)). The coarseness andfineness of the grids can adjusted according to the zone spatial scale in eachregion. If the grid systems at the interface of overlap subregions are not iden-tical, space interpolation at the interface is required when transferring theinformation at the interface. In this paper, the identical mesh structures areused for FVM and LBM for convenience to avoid the spatial interpolation(see Fig. 9-(b)). In order to implement the coupling computations, the overlapSchwartz alternative procedure is used to handle the computations.Numerical simulations were carried out for cavity flow of Re = 100 ,
400 and1000 on a grid 200 × L = 1.The boundaries of the cavity are stationary walls, except the top-boundarywith a uniform tangential velocity ( u t,Re =100 = 3 . × − , u t,Re =400 = 1 . × − , u t,Re =1000 = 3 . × − ). Fig. 10 shows plots of the stream functionfor the Reynolds number considered. These plots give a clear picture of theoverall flow pattern and the effect of Reynolds number on the structure ofthe recirculating eddies in the cavity. The smoothness of the stream functiondistribution, especially around the overlap region confirms the correctness of27he information transfer at the interface. To further quantify these results,the velocity profiles along the vertical and horizontal centerlines of the cavityare shown in Fig. 11. The results are in close agreement with the benchmarksolution [39]. The smoothness and consistency of velocity distribution in theoverlap region is presented in Fig. 12 where a local, enlarged view of thevector plot in the overlap region is shown. Clearly, the vectors in the overlapregion are quiet consistent between the LBM results and the FVM results.Figs. 13 and 14 show the contours of horizontal and vertical velocity. It is seenthat these physical quantities are all smooth across the interface. Accordingto the authors’ numerical experience, the smoothness of vorticity contour isthe most difficult to obtain for such coupled computation, because vorticity ifthe derivative of velocity. The contours of vorticity distribution are shown inFig. 15. Over all, the smoothness on the overlap region are quite good, with aminor bumpiness of the left-hand vortex contours for the case of Re = 100.In all, by the proposed lifting relation, we can couple the mesoscopic LBMwith FVM to implement the domain decomposition coupling-computations.This paves the way for implementing multiscale computations based on LBMand macro-numerical methods of finite-family.It should be noted that we also tried the coupling computations based on thedistribution function f i (x , t ) reconstructed by Eq. (61) of [9] and (62) of [32].Unfortunately, all of our tries were unsuccessful and converged solutions couldnot be obtained. 28 onclusion In this paper, we derive the relation to lift the macroscopic variables to themicroscopic variables for LBM. Two methods of derivation are conducted andthey lead to the same result. Numerical tests demonstrate that the derivedlifting relation possesses good precision. The proposed lifting relation offers away to implement the multiscale-computations involving LBM more efficientlyand robustly.
Acknowledgment
This work was supported by the Key Projects National Natural Science Foun-dation of China (51136004) and the National Basic Research Program (973)(2007CB206902). We appreciate the referee’s valuable comments on our work.
References [1] R. Benzi, S. Succi and M. Vergassola, The lattice Boltzmann equation: theoryand applications, Phys. Rep. 222 (1992) 145-197.[2]
S. Succi , The Lattice Boltzmann Equation for Fluid Dynamics and Beyond.(Oxford Uniersity Press, Oxford, UK, 2001).[3]
S. A. Orszag, H. Chen, S. Succi and J. Latt , Turbulence effects on kineticequations J. Scie. Comput. 28(213) (2006) 459-466.[4]
H. Chen, S. Kandasamy, S. A. Orszag, R. Shock, S. Succi andV.Yakhot , Extended Boltzmann kinetic equation for turbulent flows Science.301 (2003) 633-636.
5] Q. J. Kang, D. X. Zhang and S. Y. Chen , Unified lattice boltzmann method forflow in multiscale porous media, Phys. Rev. E 66 (2002) 056307.[6] M. G. Fyta, S. Melchionna, E. Kaxiras and S. Succi, Multiscale couplingof molecular dynamics and hydrodynamics: application to DNA translocationthrough a nanopore, Multiscale Model. Sim. 5 (2006) 1156-1173.[7]
S. Succi, O. Filippova, G. Smithand Kaxiras E. , Applying the LatticeBoltzmann Equation to Multi-scale Fluid Problems, Computing in Science andEngineering, 3(6) (2001), 26-37.[8]
A. Dupuis, E. M. Kotsalis and P. Koumoutsakos , Coupling LatticeBoltzmann and Molecular Dynamics Models for Dense Fluids, PHYSICALREVIEW E, 75 (2007), 046704.[9]
P. A. Skordos , Initial and boudary conditions for the lattice Boltzmannmethod, Phys. Rev. E 48(6) 1993, 4823-4841.[10]
A. Caiazzo , Analysis of lattice Boltzmann initialization routines, J. Stat. Phys.121 (2005), 37-48.[11]
R. Mei, L.-S. Luo, P. Lallemand, and D. d’Humi ` e res , Consistent initialconditions for lattice Boltzmann simulations, Computers and Fluids 35 (8/9)(2006), 855-862.[12] A.A. Mohamad and S. Succi, A note on equilibrium boundary conditions inlattice Boltzmann fluid dynamic simulations, Eur. Phys. J. Special Topics 171(2009) 213-221.[13] M. Junk and Z.X. Yang, Convergence of lattice Boltzmann methods for Navier-Stokes flows in periodic and bounded domains, Numerische Mathematik, 112(1)(2009) 65-87.[14] M. Junk and Z.X. Yang, Outflow boundary conditions for the lattice Boltzmannmethod, Progress in Computational Fluid Dynamics, 8(1/4) (2008) 38-48.
15] M. Junk and Z.X. Yang, Asymptotic analysis of lattice Boltzmann boundaryconditions, J. Stat. Phys. 121 (2005) 3-35.[16]
S. T. O’Connel, P. A. Thompson , Molecular dynamics-continuum hybridcomputations: A tool for studying complex fluid flows, Physics Review E, 52 (6)(1995), R5792-R5795.[17]
F F. Abraham
Dynamically spanning the length scales from the quantum tothe continuum, International Journal of Modern Physics C, 11 (6) (2000), 1135-1148.[18]
J. Liu, S. Y. Chen, X. B. Nie and M. O. Robbins , A continuum atomisticsimulation of heat transfer in micro- and nano-flow, J. computational Physics, 227(2007), 279-291.[19]
P. Albuquerque, D. Alemani, B. Chopard and P. Leone , Couplinga Lattice Boltzmann and a Finite Difference Scheme, Computational Science,ICCS- 04, Kracow, June 6-9, 2004. LCNS 3039, Bubak, M.; Albada, G.D.v.; Sloot,P.M.A.; Dongarra, J. (Eds.) Springer Verlag, Berlin.[20]
P. Van Leemput, W. Vanroose and D. Roose , Numerical and analyticalspatial coupling of a lattice Boltzmann model and a partial differential equation.In Model Reduction and Coarse-Graining Approaches for Multiscale Phenomena,(A.N. Gorban, N. Kazantzis, I.G. Kevrekidis, H.C. Ottinger, C. Theodoropouloseds.), p. 423-441. Springer, 2006.[21] P. Van Leemput, W. Vanroose and D. Roose, Mesoscale analysis of the equation-free constrained runs initialization scheme, (SIAM) Multiscale Modeling andSimulation, 6 (4) (2007): 1234-1255.[22]
I. G. Kevrekidis, C. W. Gear, J. M. Hyman, P. G. Kevrekidis,O. Runborg and C. Theodoropoulos , Equation-free, coase-grainedmultiscale computation: Enabling microscopic simulators to perform system-levelanalysis, Communications in Mathematical Science, 1(4) (2003) 715-762. M. D. Mazzeo, P. V. Coveney,
HemeLB: A high performance parallellattice-Boltzmann code for large scale fluid flow in complex geometries, ComputerPhysics Communications, 178 (12) (2008): 894-914.[24] G. Amati, S. Succi, R. Piva, Massively Parallel Lattice-Boltzmann Simulationof Turbulent Channel Flow, International Journal of Modern Physics C, 8 (4)(1997): 869-877.[25]
X. B. Nie, S. Y. Chen, W. N. E and M. O. Robbins , A continuum andmolecular dynamics hybrid method for micro- and nano-fluid flow. J. Fluid Mech.500 (2004), 55-64.[26]
M. Junk, A. Klar, L. S. Luo
Asymptotic analysis of the lattice Boltzmannequation, Journal of Computational Physics, 210 (2005) 676-704.[27] S. Chen, G. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev.Fluid Mech. 161 (1998) 329.[28] D. Ricot, V. Maillard, C. Bailly, Numerical simulation of unsteady cavity flowusing Lattice Boltzmann Method, in: AIAA-Paper 2002-2532, 2002.[29]
W. N. E and B. Engquist , The heterogeneous multiscale methods,Communications in Mathematical Science, 1(1) (2003) 87-133.[30]
Z. L. Guo and C. G. Zheng , Theory and application of lattice BoltzmannMethod. (Science Press, Beijing, 2008)[31]
T. Imamura, K. Suzuki, T. Nakamura and M. Yoshida , Acceleration ofsteady-state lattice Boltzmann simulation on non-uniform mesh using local timestep method, J. Comput. Phys. 202 (2005) 645-663.[32]
Z. L. Guo and T. S. Zhao , Explicit finite-difference lattice Boltzmannmethod for curvilinear coordinates, Phys. Rev. E, 67 (2003) 066709.[33]
S.Chapman and T. G. Cowling , The mathematical theory of nonuniformgases, 3rd. ed Cambridge University Press, Cambridge, 1970. F. Schwabl , Statistical Mechanics, 2nd. ed, Springer-Verlag Berlin Heidelberg,(2006).[35]
H.B. Luan, H. Xu, L. Chen, D. L. Sun, Y. L. He and W. Q. Tao ,Evaluation of the coupling scheme of FVM and LBM for fluid flows aroundcomplex geometries, Int. J. Heat Mass Tran. 54 (2011), 1975-1985.[36]
Y. H. Qian, D. d’Humieres, and P. Lallemand , Lattice BGK Models forNavier-Stokes Equation, Europhys. Lett., 17(6) (1992), 479-484.[37]
Y.H. Qian and Y. Zhou , Complete Galilean-Invariant Lattice BGK Modelsfor the Navier-Stokes Equation, Europhys. Lett., 42 (1998) 359-364.[38]
J. R. Chasnov , On the decay of two-dimensional homogeneous turbulence,Phys. Fluids, 9(1) (1997) 171-180.[39]
U. Ghia , K. N. Ghia, and C. T. Shin , High-Re Solutions for IncompressibleFlow using the Navier-Stokes Equationuations and a Multigrid Method, J. Compt.Phys., 48 (1982) 387-411. .40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.800.400.450.500.550.600.650.700.750.80 Real fi vs Restructed fi Fitting line slope = 0.99758 R ec on s t r u c t e d f i Real fi (a)Linear regression between f i (x , t ) and (cid:98) f i (x , t ) -0.03 -0.02 -0.01 0.00 0.01 0.02 0.03 0.04-0.02-0.010.000.010.02 Real fineq vs Restructed fineq Fitting line slope = 0.83265 R ec on s t r u c t e d f i n e q Real fi neq (b)Linear regression between f (neq) i (x , t ) and (cid:98) f (neq) i (x , t ) Fig. 1. Linear regression (D2Q9, u = 0 . t = 1000 δt , i = 2): (a)Fit theline (cid:98) f i (x , t ) = af i (x , t ) + b ,where a = 0 . b = 0 . (cid:98) f (neq) i (x , t ) = af (neq) i (x , t ) + b ,where a = 0 . b = − . × − . Stan-dard deviation:(a) σ = 0 . σ = 9 . × − . .40 0.45 0.50 0.55 0.60 0.65 0.70 0.750.400.450.500.550.600.650.700.75 Real fi vs Restructed fi Fitting line slope = 0.99806 R ec on s t r u c t e d f i Real fi (a)Linear regression between f i (x , t ) and (cid:98) f i (x , t ) -0.02 -0.01 0.00 0.01 0.02-0.015-0.010-0.0050.0000.0050.0100.0150.020 Real fineq vs Restructed fineq Fitting line slope = 0.84218 R ec on s t r u c t e d f i n e q Real fi neq (b)Linear regression between f (neq) i (x , t ) and (cid:98) f (neq) i (x , t ) Fig. 2. Linear regression (D2Q17, u = 0 . t = 1000 δt , i = 2): (a)Fit theline (cid:98) f i (x , t ) = af i (x , t ) + b ,where a = 0 . b = 0 . (cid:98) f (neq) i (x , t ) = af (neq) i (x , t ) + b ,where a = 0 . b = − . × − . Stan-dard deviation:(a) σ = 0 . σ = 8 . × − . .57 0.58 0.59 0.60 0.610.570.580.590.600.61 R ec on s t r u c t e d f i Real fi
Real fi vs Restructed fi Fitting line slope = 0.99925 (a)(a)Linear regression between f i (x , t ) and (cid:98) f i (x , t ) -0.002 -0.001 0.000 0.001 0.002-0.0015-0.0010-0.00050.00000.00050.00100.00150.0020 Real fineq vs Restructed fineq Fitting line slope = 0.83655 R ec on s t r u c t e d f i n e q Real fi neq (b)Linear regression between f (neq) i (x , t ) Fig. 3. Linear regression (D2Q9, u = 0 . t = 10000 δt , i = 2): (a)Fit the line (cid:98) f i (x , t ) = af i (x , t ) + b ,where a = 0 . b = 4 . × − ;(b)Fit the line (cid:98) f (neq) i (x , t ) = af (neq) i (x , t ) + b ,where a = 0 . b = − . × − . Standarddeviation:(a) σ = 3 . × − ;(b) σ = 1 . × − . .58 0.59 0.60 0.61 0.620.580.590.600.610.62 Real fi vs Restructed fi Fitting line slope = 0.99963 R ec on s t r u c t e d f i Real fi (a)Linear regression between f i (x , t ) and (cid:98) f i (x , t ) -0.0016 -0.0008 0.0000 0.0008 0.0016-0.0015-0.0010-0.00050.00000.00050.00100.0015 Real fineq vs Restructed fineq Fitting line slope = 0.84764 R ec on s t r u c t e d f i n e q Real fi neq (b)Linear regression between f (neq) i (x , t ) and (cid:98) f (neq) i (x , t ) Fig. 4. Linear regression (D2Q17, u = 0 . t = 10000 δt , i = 2): (a)Fit the line (cid:98) f i (x , t ) = af i (x , t ) + b ,where a = 0 . b = 2 . × − ;(b)Fit the line (cid:98) f (neq) i (x , t ) = af (neq) i (x , t ) + b ,where a = 0 . b = − . × − . Standarddeviation:(a) σ = 3 . × − ;(b) σ = 9 . × − . ig. 5. Linear regression (D2Q9, u = 0 . t = 1000 δt , i = 2): Fit the line (cid:98) f i (x , t ) = af i (x , t ) + b ,where a = 0 . b = − . × − . f i (x , t )(b)Vorticity contour plot by the reconstructed (cid:98) f i (x , t ) Fig. 6. Vorticity contour plots (D2Q9, u = 0 . t = 1000 δt ): (a)Vorticity contourplot by the real f i (x , t ) ; (b)Vorticity contour plot by the reconstructed (cid:98) f i (x , t ) f i (x , t )(b)Vorticity contour plot by the reconstructed (cid:98) f i (x , t ) Fig. 7. Vorticity contour plots (D2Q9, u = 0 . t = 10000 δt ): (a)Vorticity contourplot by the real f i (x , t ) ; (b)Vorticity contour plot by the reconstructed (cid:98) f i (x , t ) −3 −2 −1 −5 −4 −3 −2 Mesh scale h L & L r e l a t i v e e rr o r s L relative errorL relative errorRef. line slope = 2.6 (a) the D2Q9 model −3 −2 −1 −6 −5 −4 −3 Mesh scale h L & L r e l a t i v e e rr o r s L relative errorL relative errorRef. line slope = 2.6 (b) the D2Q17 model Fig. 8. Convergence rates of the reconstruction operator for D2Q9 and D2Q17. × Fig. 9. Geometric structure and mesh partition: (a)Interface structure between tworegions of FVM and LBM; (b)Grid layout for a 2D lid-driven cavity (200 × Re = 100 (b) Re = 400 (c) Re = 1000 Fig. 10. Contour plots of streamline for different Reynolds numbers (a) Horizontal velocity profiles43b)Vertical velocity profiles