Macroscopic Lattice Boltzmann Method (MacLAB)
MMacroscopic Lattice Boltzmann Method (MacLAB)
Jian Guo ZhouDepartment of Computing and MathematicsManchester Metropolitan UniversityManchester, M1 5GD, [email protected]
Abstract
The birth of the lattice Boltzmann method (LBM) fulfils a dream that simple arithmeticcalculations can simulate complex fluid flows without solving complicated partial differentialflow equations. Its power and potential of resolving more and more challenging physicalproblems have been and are being demonstrated in science and engineering covering a widerange of disciplines such as physics, chemistry, biology, material science and image analysis.The method is a highly simplified model for fluid flows using a few limited fictitious particlesthat move one grid at a constant time interval and collide each other at a grid point onuniform lattices, which are the two routine steps for implementation of the method tosimulate fluid flows. The former represents fluids movement and the latter provides fluidviscosity for diffusion effect. As such, a real complex particle dynamics is approximated asa regular particle model using three parameters of lattice size, particle speed and collisionoperator. A fundamental question is “Are the two steps integral to the method or can thethree parameters be reduced to one for a minimal lattice Boltzmann method?”. Here, I showthat the collision step can be removed and the standard LBM can be reformulated into asimple macroscopic lattice Boltzmann method (MacLAB). This model relies on macroscopicphysical variables only and is completely defined by one basic parameter of lattice size δx ,bringing the LBM into a precise “Lattice” Boltzmann method. The viscous effect on flowsis naturally embedded through the particle speed, making it an ideal automatic simulatorfor fluid flows. Three additional advantages compared to the existing LBMs are that (i)physical variables can directly be retained as the boundary conditions; (ii) computationalmemory are much less required and (iii) the model is unconditional stable. The findings havebeen demonstrated and confirmed with numerical tests including flows that are independentof and dependent on fluid viscosity, 2D and 3D cavity flows, and an unsteady Taylor-Greenvortex flow. This provides an efficient and powerful model for resolving physical problemsin various disciplines of science and engineering. The LBM is characterised by its simplicity, parallel processing, and easy treatment of bound-ary conditions [1]. The first fully discrete model for fluid flows on a square lattice was proposedby Hardy et al. [2] in 1976. Ten years later, Frisch et al. [3] for the first time obtained a correct1 a r X i v : . [ c s . C E ] J a n attice gas automata (LGA) for Navier-Stokes equations using six-velocity hexagonal lattice. TheLGA comprises two steps: streaming and collision. The two steps are represented by lattice sizeand a collision operator on a uniform lattice. In physics, the former and the latter simulate thephenomena of fluid movement and diffusion, respectively, which determine the basic feature ofa LGA. Often, simulations generated using a LGA are very noisy due to its Boolean variablewith one for the presence and zero for the absence of particles [4, 5]. Also, the numerical proce-dure involves calculations of particle probability, which reduces the efficiency of the model. Toovercome these, the lattice Boltzmann method was proposed [6] and its basic difference fromthe LGA is that the Boolean variable is replaced with a particle distribution function. Such ap-proach eliminates the statistical noise in a LGA and retains all the advantages of locality in thekinetic form of a LGA [1]. McNamara and Zanetti [6] first used the lattice Boltzmann methodas an alternative to the LGA in 1988. As the collision operator takes a complex matrix form,this prevents the LBM from becoming a competing computational method. A breakthroughprogress has been made by Higuera and Jim´enez [7] who linearized the collision term aroundits local equilibrium state. This greatly simplifies the collision operator. Noble et al. [8] usedthis idea to express the collision operator as Ω αβ ( f eqβ − f β ), in which f β is the particle distribu-tion function; f eqβ is the local equilibrium distribution function; and Ω αβ is a collision matrix.Later, several researchers [9, 10] suggested a simple linearized form for the collision matrix byusing a single time relaxation towards the local equilibrium distribution, Ω αβ = − δ αβ /τ , whichis the Bhatnagar-Gross-Krook [11] collision operator. In the operator, δ αβ is the Kronecker deltafunction taking one when α = β , or otherwise zero, and τ is called the single relaxation timetaking a constant and is related to fluid viscosity. This leads to the single relaxation time latticeBoltzman method (SRT LBM) that has been most efficient and used so far, f α ( x j + e αj δt, t + δt ) = f α ( x j , t ) + 1 τ [ f eqα ( x j , t ) − f α ( x j , t )] , (1)where x j is a lattice coordinate along j -axis in Cartesian coordinate system, e.g., j = x, y in the two dimensional space; t is time; e αj is the j th component of the particle velocity vec-tor e α in α − link of the lattice and defined by time step δt and lattice size δx , e.g., e α =(0 , , ( e, , (0 , e ) , ( − e, , (0 , − e ) , ( e, e ) , ( − e, e ) , ( − e, − e ) , ( e, − e ) when α = 0 − e is the particlespeed and defined as e = δx/δt ; and f eqα is the local equilibrium distribution function given by f eqα = w α ρ (cid:18) e αi u i e + 92 e αi e αj u i u j e − u i u i e (cid:19) , (2)in which ρ is the fluid density and w α is a weighting factor depending on lattice pattern, e.g., w α = 4 / α = 0, w α = 1 / α = 1 − w α = 1 /
36 when α = 5 − ρ ( x j , t ) = (cid:88) α f α ( x j , t ) , u i ( x j , t ) = 1 ρ ( x j , t ) (cid:88) α e αi f α ( x j , t ) . (3)Since then, the study on the lattice Boltzmann method and applications of the method havereceived extensive attentions, making it become a very powerful modelling tool in many areassuch as thermodynamics [12], aerodynamics [13], multiphase flows [14], turbulent flows [15],2emodynamics [16], biomechanics [17], image analysis [18], biology [19], environmental science[20].In applications, it is found that the SRT LBM suffers from a numerical instability. To remedythis, the multiple-relaxation-time (MRT) collision operator was introduced in 1992 [21, 22]. Thisimproves the stability but it reduces the efficiency. To accelerate simulation, a two-relaxation-time collision operator was developed in 2008 [23], which has almost the same efficiency as theSRT LBM. As research progresses, it is noticed that MRT or TRT still suffers from numericalinstabilities when fluid flows with very small viscosity are simulated. After realising that thisis caused by an insufficient degree of Galilean invariance in the collision step, Geier et al. [24]proposed a cascaded lattice Boltzmann method by relaxaing the particle distribution function toits local equilibrium state in the central moment space, making the LBM stable for simulatingflows with small viscosity close to zero. In 2015, Geier et al. [25] further improved their centralmoment LBM using the cumulant in collision operator, which is called the cumulant latticeBoltzmann method (CLBM). Despite such enhancements, these schemes are more complicatedand computational efficiency is reduced as manipulation of matrix is involved. In addition, theyshare the same drawback as an existing LBM in that the boundary conditions for a physicalvariable such as velocity cannot be implemented without being converted to particle distributionfunctions, which further reduces efficiency and accuracy of the methods. Recently, Chen et al.[26] developed a simplified lattice Boltzmann method without evolution of particle distributionfunction, which successfully removes this drawback and enables a direct use of a physical variableas boundary conditions. However, the method involves the two steps of predictor and corrector,which is more complicated than the SRT LBM. Nevertheless, through all these research, theLBM has greatly been improved and developed to a point where it has become a very efficientand flexible alternative numerical method in computational physics. Its potential power is muchbeyond the original scope, being explored and demonstrated in various disciplines of science andengineering with time [27, 28, 29, 30].The above literature review highlights that the major research to improve the method havebeen carried out on the collision operator for resolving the well-known instability problem in theSRT LBM since 1992, leading to four representative variants of the method, MRT, TRT, Centralmoment and Cumulant LBMs. Even so, choosing suitable parameters or values for the collisionoperators is not clear and they cannot be tuned without using trial and error during simulations,which becomes more complicated in a non-single relaxation time scheme due to its complexity.This unnecessarily wastes time and computing resource, and it may become awkward to simulatea large-scale flow system, preventing the LBM from becoming an automatic simulator for anyscale flows when a super-fast computer such as a quantum computer becomes available oneday. In principle, “Everything should be made as simple as possible, but not simpler” (AlbertEinstein’s quote). Since the central problem comes from the collision operator, the problem maybe resolved forever if the collision operator can be removed. Due to the fact that the function ofthe collision is to relax the distribution function to its local equilibrium state, one may removethe collision operator and use the local equilibrium distribution function to replace the collisionby setting τ = 1 in Eq. (1). Following this idea, after mathematical manipulations, I obtain thefollowing simple macroscopic lattice Boltzmann model (MacLAB), ρ ( x j , t ) = (cid:88) f eqα ( x j − e αj δt, t − δt ) , (4)3 i ( x j , t ) = 1 ρ ( x j , t ) (cid:88) e αj f eqα ( x j − e αj δt, t − δt ) , (5)to determine physical density and velocity directly from the local equilibrium distribution func-tion without calculating the distribution function using Eq. (1) unlike existing LBMs (The the-oretical detail can be found in Methods). Apparently, the model involves the local equilibriumdistribution function only. However, doing so brings a new problem of how to consider the fluidviscosity in the absence of collision operator. This can be overcome from the findng throughthe recovery of the Navier-Stokes equations from Eqs. (4) and (5) that if the particle speed e isdetermined using e = 6 ν/δx, (6)and is employed in Eq. (2) instead of using e = δx/δt to calculate the local equilirium distributionfunction f eqα , the flow viscosity is naturally taken into account in the model. In this case, oncea lattice size δx is chosen, the model is ready to simulate a flow with a viscosity ν because( x j − e αj δt ) stands for a neighbouring lattice point; f eqα at time of ( t − δt ) represents its knownquantity at the current time; and the particle speed e is determined from Eq. (6) for use incomputation of f eqα . In addition, the time step δt is no longer an independent parameter butis calculated as δt = δx/e , which is used in simulations of unsteady flows. Consequently, onlythe lattice size δx is required in the MacLAB for simulation of fluid flows, bringing the LBMinto a precise “Lattice” Boltzmann method. This enables the model to become an automaticsimulator for any scale flows without tuning other simulation parameters, making it possibleand easy to model a large flow system when a super-fast computer such as a quantum computerbecomes available in the future. The model is unconditionally stable as it shares the same validcondition as that for f eqα , or the Mack number M = U c /e is much smaller than 1, in which U c is a characteristic flow speed. The Mack number can also be expressed as a lattice Reynoldsnumber of R le = U c δx/ν via Eq. (6). In practical simulations, it is found that the model is stableif R le = U m δx/ν < U m is the maximum flow speed and is used as the characteristicflow speed. The main features are that there is no collision operator and only macroscopicphysical variables such as density and velocity are required, which are directly used as boundaryconditions with a minimum memory requirement. The implementation of the model startingfrom the initial density and velocity is to (i) choose the lattice size δx and determine the particlespeed e from Eq. (6), (ii) calculate f eqα from Eq. (2) using density and velocity, (iii) update thedensity and velocity using Eqs. (4) and (5), (iv) apply the boundary conditions, and (v) repeatStep (ii) until a solution is reached. The only limitation of the described model is that, for verysmall viscosity or high speed flow, the chosen lattice size after satisfying R le < L iscalculated as N L = L/δx and N L is the lattice points); if the total lattice points is too big suchthat the demanding computations is beyond the current power of a computer, the simulationcannot be carried out. Such difficulties may be solved or relaxed through parallel computingusing computer techniques such as GPU processors and multiple servers, and will largely orcompletely removed using quantum computing when a quantum computer becomes available.In order to demonstrate the validation of the described model, I have carried out five numericalsimulations using D2Q9 and D3Q19 lattices for 2D and 3D flows, respectively. For D3Q19, theparticle velocity vector is e α = ( e αx , e αy , e αz ) = (0 , , , ( e, , , ( − e, , , (0 , e, , (0 , − e, , (0 , , e ) , (0 , , − e ) , ( e, e, , ( − e, − e, , ( − e, e, , ( e, − e, , (0 , e, e ) , (0 , − e, − e ) , (0 , − e, e ) , (0 , e, − e ) , ( e, , e ) , ( − e, , − e ) , ( − e, , e ) , ( e, , − e ), and the weighting factor w α for Eq. (2) is4 α = 1 / α = 0, w α = 1 /
18 when α = 1 − w α = 1 /
36 when α = 7 −
18. The SIunits are used with ρ = 1 in the numerical simulations. The first test is a Couette flow throughtwo parallel plates without a pressure gradient. The distance between the plates is h = 1. Thetop plate moves at velocity of u x = u = 0 . x stands for the streamwise direction and y for the vertical direction, the analyticalsolution is u x ( y ) = u h y, (7)which is the same test as that used by Chen et al. [26]. This is a very interested case as the steadyflow is independent of flow viscosity according to the theory (7). I use δx = 0 .
02 and 20 × x and y directions for three simulations of flows with three kinematic viscosities of ν = 0 . , ν = 0 .
001 and ν = 0 . ∂p/∂x = − . δx/ ( eρ ) ∂p/∂x [31].Both plates are fixed with zero velocities at top and bottom boundaries at which no calculationsare needed. The flow is affected by viscosity and the analytical solution is u ( y ) = u h y + 12 ρν ∂p∂x ( y − hy ) . (8)I simulate this flow using three viscosities of ν = 0 . , ν = 0 . , and ν = 0 . × u x = u and u y = 0 with u = 1; the other three sides are fixed, or no slip boundarycondition is applied, i.e., u x = 0 and u y = 0. The Reynolds number R e = u /ν = 1000. I use δx = 0 . ×
400 lattices in the simulation, which is carried out on the inside of thecavity excluding the four sides where velocities are retained as boundary conditions. After thesteady solution is obtained, the flow pattern in velocity vectors is shown in Fig 3, which closelyagrees with the well-known study by Ghia et al. [32]. The results are further compared againsttheir numerical solution for velocity profiles of u x and u y along y and x directions through thegeometric centre of the cavity in Figs. 4 and 5, respectively, demonstraing very good agreements.The fourth test is a 2D Taylor-Green vortex. This is an unsteady flow driven by decayingvortexes for which there is an exact solution of the incompressible Navier-Stokes equations andit is often applied to validation of a numerical method for solution to the incompressible Navier-Stokes equations. The initial conditions are u x ( x, y,
0) = − u cos( x ) sin( y ) and u y ( x, y,
0) = u sin( x ) cos( y ). The analytical solution are u x ( x, y, t ) = − u cos( x ) sin( y ) exp( − νt ) and u y ( x, y, t ) = u sin( x ) cos( y ) exp( − νt ). The time for an unsteady flow from initial state isaccumulated by its increase with time step δt . I use δx = 0 .
157 or 40 ×
40 lattices for square5omain of 2 π × π with kinematic viscosity of ν = 0 . u = 0 .
05, which gives the Reynoldsnumber of R e = 2 πu /ν = 10. The periodic boundary conditions are used. The simulation isrun for the total time of 30 seconds. The velocity field is plotted in Fig. 6, showing correctflow pattern. The velocity profiles for u x at x = π and u y at x = π/ y -direction aredepicted and compared with the analytical solutions in Fig. 7, showing excellent agreemeentsand confirming the accuracy of the method for an unsteady flow.The final test is a 3D cavity flow. This is again a well-known complex flow involving 3Dvortices within a simple cube with the dimensions of 1 × × x , spanwisedirection y and vertical direction z . No-slip boundary conditions, i.e, u x = 0 and u y = 0 and u z = 0, are applied to five fixed sides except for the top lid, where u x = u , u y = 0 and u z = 0with u = 1 are specified. The Reynolds number is R e = u /ν = 400. δx = 0 .
004 or total latticesof 250 × ×
250 are used and the simulation is undertaken only within the cube excludingthe boundaries where the velocities are retained. After the steady solution is reached, the flowcharacterises are displayed through the two dimensional planar projections of the velocity vectorfield on the x - z , y - z and x - y centroidal planes of the cube in Figs. 8, 9 and 10, respectively,demonstrating good agreed flow patterns with those by Wong and Baker [33]. In addition, thedistribution of the velocity component u x on the vertical plane centerline is widely used as a3D lid-driven cavity benchmark test. I compare this velocity component against the results byWong and Baker [33] and also by Jiang et al. [34] in Fig. 11, showning good agreements.In conclusion, the results demonstrate that the MacLAB is able to simulate fluid flows usingonly lattice size, bringing the LBM into a precise Lattice Boltzmann method. This takes theresearch on the method into a new era when future work may focus on improving on accuracy ofor formulating a new local equilibrium distribution function. The particle speed is determinedthrough the viscosity and lattice size and the time step δt is calculated as δt = δx/e . Themodel is unconditional stable as long as the valid condition for the local equilibrium distributionfunction holds. All these make the method an automatic simulator for fluid flows. The methodis straightforward to be extended for resolving other physical problems in different disciplines. Figures. y ( m ) u x (m/s)MacLABAnalytical Figure 1: Couette flow through two parallel plates without a pressure gradient. The distancebetween the plates is h = 1. The top plate moves at velocity of 0 . δx = 0 .
02 is used for three simulations of flowswith three kinematic viscosities of ν = 0 . , ν = 0 .
001 and ν = 0 . y ( m ) u x (m/s) n n n MacLABAnalytical
Figure 2: Couette flow through two parallel plates with a pressure gradient of ∂p/∂x = − . h = 1. Both plates are fixed with zero velocities at top andbottom boundaries where no calculations are needed. The steady numerical results are dependenton flow viscosity as confirmed in the simulations using the three viscosities of ν = 0 . , ν =0 . , and ν = 0 . y ( m ) x (m) Figure 3: 2D cavity flow within 1 × R e = 1000. The top lid moves at velocity of u x = 1 and u y = 0 and the other three sides are fixed, or no slip boundary condition is applied.After the steady solution is obtained, the flow pattern in velocity vectors shows a primary vortexand two secondary vortices. y ( m ) u x (m/s) MacLABGhia et al. Figure 4: 2D cavity flow within 1 × R e = 1000. The top lid moves at velocity of u x = 1 and u y = 0 and the other three sides are fixed, or no slip boundary condition is applied.After the steady solution is obtained, the comparison of velocity u x profile along y directionthrough the geometric centre of the cavity with the numerical solution by Ghia et al. [32].9 u y ( m / s ) x (m)MacLABGhia et al. Figure 5: 2D cavity flow within 1 × R e = 1000. The top lid moves at velocity of u x = 1 and u y = 0 and the other three sides are fixed, or no slip boundary condition is applied.After the steady solution is obtained, the comparison of velocity u y profile along x directionthrough the geometric centre of the cavity with the numerical solution by Ghia et al. [32]. p p p p Figure 6: Taylor-Green vortex within 2 π × π domain for R e = 10. The initial conditions are u x ( x, y,
0) = − u cos( x ) sin( y ) and u y ( x, y,
0) = u sin( x ) cos( y ) with u = 0 .
05. The periodicboundary conditions are used. Here shown is the flow pattern in velocity vectors at t = 30seconds, remaining the same vortex pattern as that at initial state.10 p )MacLAB u x /u MacLAB u y /u Analytical
Figure 7: Taylor-Green vortex within 2 π × π domain for R e = 10. The initial conditions are u x ( x, y,
0) = − u cos( x ) sin( y ) and u y ( x, y,
0) = u sin( x ) cos( y ) with u = 0 .
05. The periodicboundary conditions are used. Here shown is the comparisons of the relative velocity profilesfor u x /u at x = π and u y /u at x = π/ u x ( x, y, t ) = − u cos( x ) sin( y ) exp( − νt ) and u y ( x, y, t ) = u sin( x ) cos( y ) exp( − νt ) at t = 30 seconds. z ( m ) x (m) Figure 8: 3D cavity flow within 1 × × R e = 400. The top lid moves at velocity of u x = 1, u y = 0 and u z = 0 and the other five sides are fixed, or no slip boundary condition isapplied. After the solution is reached, the flow pattern in vectors in x − z centroidal plane showsthe primary and secondary vortices. 11 z ( m ) y (m) Figure 9: 3D cavity flow within 1 × × R e = 400. The top lid moves at velocity of u x = 1, u y = 0 and u z = 0 and the other five sides are fixed, or no slip boundary condition isapplied. After the solution is reached, the flow pattern in vectors in y − z centroidal plane showsone pair of strong secondary vortices at bottom and one pair of weak secondary vortices at top. y ( m ) x (m) Figure 10: 3D cavity flow within 1 × × R e = 400. The top lid moves at velocity of u x = 1, u y = 0 and u z = 0 and the other five sides are fixed, or no slip boundary condition isapplied. After the solution is reached, the flow pattern in vectors in x − y centroidal plane showsa pair of third vortices close to inflow boundary.12 z ( m ) u x (m/s)MacLABJiang et al.Wong & Baker Figure 11: 3D cavity flow within 1 × × R e = 400. The top lid moves at velocityof u x = 1, u y = 0 and u z = 0 and the other five sides are fixed, or no slip boundary conditionis applied. After the solution is reached, the comparisons of the distribution of the velocitycomponent u x on the vertical plane centerline with the restuls by Wong and Baker [33] andJiang et al. [34]. 13 ETHODS
I present the detail of the derivation for the present model. Setting τ = 1 in Eq. (1) leads to f α ( x j + e αj δt, t + δt ) = f eqα ( x j , t ) , (9)which can be rewritten as f α ( x j , t ) = f eqα ( x j − e αj δt, t − δt ) . (10)Taking (cid:80) Eq. (10) and (cid:80) e αi Eq. (10) yields (cid:88) f α ( x j , t ) = (cid:88) f eqα ( x j − e αj δt, t − δt ) , (11)and (cid:88) e αi f α ( x j , t ) = (cid:88) e αi f eqα ( x j − e αj δt, t − δt ) , (12)respectively. In the lattice Boltzmann method, the density and velocity are determined usingthe distribution function as ρ ( x j , t ) = (cid:88) α f α ( x j , t ) , u i ( x j , t ) = 1 ρ (cid:88) α e αi f α ( x j , t ) . (13)Combining Eq. (13) with Eqs. (11) and (12) results in the current MacLAB, Eqs. (4) and (5).Since the local equilibrium distribution function f eqα has the features of (cid:88) α f eqα ( x j , t ) = ρ ( x j , t ) , ρ (cid:88) α e αi f eqα ( x j , t ) = u i ( x j , t ) , (14)with reference to Eq. (13) the following relationships, (cid:88) α f α ( x j , t ) = (cid:88) α f eqα ( x j , t ) , (cid:88) α e αi f α ( x j , t ) = (cid:88) α e αi f eqα ( x j , t ) , (15)hold, which are the conditions that retain the conservation of the mass and momentum in thelattice Boltzmann method.Next, I prove that the continuity and Navier-Stokes equations can be recovered from Eqs. (4)and (5). Rewriting Eq. (1) as f α ( x j , t ) = f α ( x j − e αj δt, t − δt ) + 1 τ [ f eqα ( x j − e αj δt, t − δt ) − f α ( x j − e αj δt, t − δt )] . (16)Apparently, when τ = 1, the above equation becomes Eq. (10) that leads to Eqs. (4) and (5);hence Eq. (16) is a general equation and is used in the following derivation. Applying a Taylorexpansion to the two terms on the right hand side of Eq. (16) in time and space at point ( x , t )yields f α ( x j − e αj δt, t − δt ) = f α ( x j , t ) − δt (cid:32) ∂∂t + e αj ∂∂x j (cid:33) f α + 12 δt (cid:32) ∂∂t + e αj ∂∂x j (cid:33) f α + O ( δt ) , (17)and f eqα ( x j − e αj δt, t − δt ) = f eqα ( x j , t ) − δt (cid:32) ∂∂t + e αj ∂∂x j (cid:33) f eqα + 12 δt (cid:32) ∂∂t + e αj ∂∂x j (cid:33) f eqα + O ( δt ) . (18)14ccording to the Chapman-Enskog analysis, f α can be expanded around f (0) α f α = f (0) α + f (1) α δt + f (2) α δt + O ( δt ) . (19)After substituting Eqs. (17), (18) and (19) into Eq. (16), equating the coefficients results in forthe order ( δt ) f (0) α = f eqα , (20)for the order ( δt ) (cid:32) ∂∂t + e αj ∂∂x j (cid:33) f (0) α , = − τ f (1) α , (21)and for the order ( δt ) (cid:32) ∂∂t + e αj ∂∂x j (cid:33) f (1) α − (cid:32) ∂∂t + e αj ∂∂x j (cid:33) f (0) α = − τ f (2) α + 1 τ (cid:32) ∂∂t + e αj ∂∂x j (cid:33) f (1) α . (22)Substitution of Eq. (21) into the above equation gives (cid:32) ∂∂t + e αj ∂∂x j (cid:33) f (1) α − (cid:32) ∂∂t + e αj ∂∂x j (cid:33) ( − τ f (1) α ) = − τ f (2) α + 1 τ (cid:32) ∂∂t + e αj ∂∂x j (cid:33) f (1) α , (23)which is rearranged as (cid:18) − τ (cid:19) (cid:32) ∂∂t + e αj ∂∂x j (cid:33) f (1) α = − τ f (2) α . (24)From Eq. (21) + Eq. (24) × δt , I have (cid:32) ∂∂t + e αj ∂∂x j (cid:33) f (0) α + δt (cid:18) − τ (cid:19) (cid:32) ∂∂t + e αj ∂∂x j (cid:33) f (1) α = − τ ( f (1) α + δtf (2) α ) . (25)Now taking (cid:80) Eq. (25) leads to (cid:88) (cid:32) ∂∂t + e αj ∂∂x j (cid:33) f (0) α = 0 (26)as (cid:88) α f (1) α = (cid:88) α f (2) α = (cid:88) α e αi f (1) α = (cid:88) α e αi f (2) α = 0 (27)due to the condition of conservation of mass and momentum Eq. (14). Evaluating the terms inthe above equation using Eq. (2) produces the exact continuity equation, ∂ρ∂t + ∂ ( ρu j ) ∂x j = 0 . (28)Multipling Eq. (25) by e αi provides (cid:32) ∂∂t + e αj ∂∂x j (cid:33) e αi f (0) α + δt (cid:18) − τ (cid:19) (cid:32) ∂∂t + e αj ∂∂x j (cid:33) e αi f (1) α = − τ ( e αi f (1) α + δte αi f (2) α ) . (29)15aking (cid:80) Eq. (29) leads to (cid:88) (cid:32) ∂∂t + e αj ∂∂x j (cid:33) e αi f (0) α + δt (cid:18) − τ (cid:19) (cid:88) (cid:32) ∂∂t + e αj ∂∂x j (cid:33) e αi f (1) α = 0 (30)under the same condition (27) as that in the derivation of Eq. (26). Evaluating the terms inthe above equation using Eq. (2) produces the exact momentum equation, the Navier-Stokesequation at second-order accuracy on condition that the Mack number M = U c /e << ∂ ( ρu i ) ∂t + ∂ ( ρu i u j ) ∂x j = − ∂p∂x i + ν ∂ ( ρu i ) ∂x j , (31)where pressure p is defined as p = 13 ρe (32)and the kinmatic viscosity is ν = 16 (2 τ − eδx. (33)As τ takes a constant, use of τ = 1 will recovers the continuity and the Navier-Stoker equations atthe second-order accurate as the above derivation shows. In this case, Eq. (33) becomes Eq. (6),which determines the particle speed e . ADDITIONAL INFORMATIOND2Q9 and D3Q19 Lattice Structures . The D2Q9 uniform square and D3Q19 cubic latticesare depicted in Figs. 12a and 12b, respectively.
123 4 567 8 (a) D2Q9 square lattice. (b) D3Q19 cubic lattice.
Figure 12: Square and cubic Lattices for 2D and 3D flows.
References [1] S. Chen and G. D. Doolen. Lattice Boltzmann method for fluid flows.
Annual Review ofFluid Mechanics , 30:329–364, 1998. 162] J. Hardy, O. de Pazzis, and Y. Pomeau. Molecular dynamics of a classical lattice gas:Transport properties and time correlation functions.
Phys. Rev. A , 13:1949–1961, 1976.[3] U. Frisch, B. Hasslacher, and Y. Pomeau. Lattice-gas automata for the Navier-Stokes equa-tion.
Physical Review Letters , 56:1505–1508, 1986.[4] B. Chopard and M. Droz.
Cellular Automata Modeling of Physical Systems . CambridgeUniversity Press, UK, 1998.[5] J. P. Rivet and J. P. Boon.
Lattice Gas Hydrodynamics . Cambridge University Press, UK,2001.[6] G. R. McNamara and G. Zanetti. Use of the Boltzmann equation to simulate lattice-gasautomata.
Phys. Rev. Lett. , 61:2332–2335, 1988.[7] F. Higuera and J. Jim´enez. Boltzmann approach to lattice gas simulations.
Europhys lett. ,9:663–668, 1989.[8] D. R. Noble, S. Chen, J. G. Georgiadis, and R. O. Buckius. A consistent hydrodynamicboundary condition for the lattice Boltzmann method.
Physics of Fluids , 7:203–209, 1995.[9] Y. H. Qian.
Lattice Gas and lattice kinetic theory applied to the Navier-Stokes equations .PhD thesis, Universit´e Pierre et Marie Curie, Paris, 1990.[10] S. Chen, H. D. Chen, D. Martinez, and W. Matthaeus. Lattice Boltzmann model forsimulation of magnetohydrodynamics.
Phys. Rev. Lett. , 67:3776–3779, 1991.[11] P. L. Bhatnagar, E. P. Gross, and M. Krook. A model for collision processes in gases. i: smallamplitude processes in charged and neutral one-component system.
Phys. Rev. , 94:511–525,1954.[12] Ahad Zarghami and Harry E. A. Van den Akker. Thermohydrodynamics of an evaporatingdroplet studied using a multiphase lattice Boltzmann method.
Physical Review E , 95:043310,2017.[13] H. H. Mohamad and M. Masoud. Continuous and discrete adjoint approach based on latticeBoltzmann method in aerodynamic optimization part i: Mathematical derivation of adjointlattice Boltzmann equations.
Advances in Applied Mathematics and Mechanics , 6(5):570–589, 2014.[14] Xiaowen Shan and Hudong Chen. Lattice boltzmann model for simulating flows with mul-tiple phases and components.
Phys. Rev. E , 47:1815–1819, Mar 1993.[15] Hudong Chen, Satheesh Kandasamy, Steven Orszag, Rick Shock, Sauro Succi, and VictorYakhot. Extended Boltzmann kinetic equation for turbulent flows.
Science , 301:633–636,2003.[16] D. R. Golbert, P. J. Blanco, A. Clausse, and R. A. Feij´oo. Tuning a lattice-Boltzmannmodel for applications in computational hemodynamics.
Medical Engineering & Physics ,34:339–349, 2012. 1717] Sana Javed, Ayesha Sohail, Khadija Maqbool, Saad Ihsan Butt, and Qasim Ali Chaudhry.The lattice Boltzmann method and computational analysis of bone dynamics-i.
ComplexAdaptove Systems Modeling , 5:1–14, 2017.[18] Junhui Chen, Zhenhua Chai, Baochang Shi, and Wenhuan Zhang. Lattice Boltzmannmethod for filtering and contour detection of the natural images.
Computers and Math-ematics with Applications , 68:257–268, 2014.[19] M. Finck, D. H¨anel, and I. Wlokas. Simulation of nasal flowby lattice Boltzmann methods.
Computers in Biology and Medicine , 37:739–749, 2007.[20] J. G. Zhou, P. M. Haygarth, P. J. A. Withers, C. J. A. Macleod, P. D. Falloon, K. J. Beven,M. C. Ockenden, K. J. Forber, M. J. Hollaway, R. Evans, A. L. Collins, K. M. Hiscock,C. Wearing, R. Kahana, and M. L. Villamizar Velez. Lattice Boltzmann method for thefractional advection-diffusion equation.
Physical Review E , 93:043310, 2016.[21] D. d’Humi`eres. Generalized lattice Boltzmann equations. in rarefied gas dynamics. In B. D.Shizgal and D. P.Weaver, editors,
Rarefied Gas Dynamics: Theory and Simulations, Progressin Astronautics and Aeronautics , volume 159, pages 450–458. 1992.[22] Pierre Lallemand and Li-Shi Luo. Theory of the lattice Boltzmann method: Dispersion,dissipation, isotropy, galilean invariance, and stability.
Physical Review E , 61:6546–6562,2000.[23] Irina Ginzburg, Frederik Verhaeghe, and Dominique d’Humi`eres3. Two-relaxation-time lat-tice Boltzmann scheme: About parametrization,velocity, pressure andmixed boundary con-ditions.
Communications in Computational Physics , 3(2):427–478, 2008.[24] Martin Geier, Andreas Greiner, and Jan G. Korvink. Cascaded digital lattice Boltzmannautomata for high reynolds number flow.
Physical Review E , 73:066705, 2006.[25] Martin Geier, Martin Sch¨onherr, Andrea Pasquali, and Manfred Krafczyk. The cumulantlattice Boltzmann equation in three dimensions: Theory and validation.
Computers andMathematics with Applications , 70:507–547, 2015.[26] Z. Chen, C. Shu, Y. Wang, L. M. Yang, and D. Tan. A simplified lattice Boltzmann methodwithout evolution of distribution function.
Advances in Applied Mathematics and Mechanics ,9:1–22, 2017.[27] Sauro Succi.
The Lattice Boltzmann Equation for Fluid Dynamics and Beyond . OxfordUniversity Press, 2001.[28] Dieter Wolf-Gladrow.
Lattice-Gas Cellular Automata and Lattice Boltzmann Models .Springer Verlag, 2000.[29] Z. L. Guo and C. Shu.
Lattice Boltzmann Method and Its Applications in Engineering .World Scientific Publishing, 2013.[30] T. Kr¨uger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, and E. M. Viggen.
TheLattice Boltzmann Method: Principles and Practice . Springer Verlag, 2017.1831] J. G. Zhou.
Lattice Boltzmann Methods for Shallow Water Flows . Springer-Verlag, Berlin,2004.[32] U. Ghia, K.N. Ghia, and C.T. Shin. High-Re solutions for incompressible flow using theNavier-Stokes equations and a multigrid method.
Journal of Computational Physics , 48:387–411, 1982.[33] K. L. Wong and A. J. Baker. A 3d incompressible navierstokes velocityvorticity weak form2nite element algorithm.
International Journal for Numerical Methods in Fluids , 38:99–123,2002.[34] B. N. Jiang, T. L. Lin, and L. A. Povinelli. Large scale computation of incompressible viscousflows by least-squares finite element method.