Macroscopic Lattice Boltzmann Method

: The lattice Boltzmann method (LBM) is a highly simplified model for fluid flows using a few limited fictitious particles. It has been developed into a very efficient and flexible alternative numerical method in computational physics, demonstrating its great power and potential for resolving more and more challenging physical problems in science and engineering covering a wide range of disciplines such as physics, chemistry, biology, material science and image analysis. The LBM is implemented through the two routine steps of streaming and collision using the three parameters of the lattice size, particle speed and collision operator. A fundamental question is if the two steps are integral to the method or if the three parameters can be reduced to one for a minimal lattice Boltzmann method. In this paper, it is shown that the collision step can be removed and the standard LBM can be reformulated into a simple macroscopic lattice Boltzmann method (MacLAB). This model relies on macroscopic physical variables only and is completely defined by one basic parameter of the lattice size δ x , bringing the LBM into a precise “lattice” Boltzmann method. The viscous effect on flows is naturally embedded through the particle speed, making it an ideal automatic simulator for fluid flows. Three additional advantages compared to the existing LBMs are that: (i) physical variables can directly be retained as the boundary conditions; (ii) much less computational memory is required; and (iii) the model is unconditionally stable. The findings are demonstrated and confirmed with numerical tests including flows that are independent of and dependent on fluid viscosity, 2D and 3D cavity flows and an unsteady Taylor–Green vortex flow. This provides an efficient and powerful model for resolving physical problems in various disciplines of science and engineering.


Introduction
The birth of the lattice Boltzmann method (LBM) fulfils a dream that simple arithmetic calculations can simulate complex fluid flows without solving complicated partial differential flow equations.The method is a highly simplified discrete model for fluid flows using a few limited fictitious particles that move one grid at a constant time interval and collide with each other at a grid point on uniform lattices, which are the two routine steps for the implementation of the method to simulate fluid flows.As such, a real complex particle dynamics is approximated as a regular particle model using three parameters of the lattice size, particle speed and collision operator.The LBM is characterised by its simplicity, parallel processing and easy treatment of boundary conditions [1].The first fully discrete model for fluid flows on a square lattice was proposed by Hardy et al. [2] in 1976.Ten years later, Frisch et al. [3] for the first time obtained a correct lattice gas automata (LGA) for Navier-Stokes equations using a six-velocity hexagonal lattice.The LGA is comprised of two steps: streaming and collision.The two steps are represented by the lattice size, particle speed and a collision operator on a uniform lattice.In physics, the former and the latter simulate the phenomena of fluid movement and diffusion, respectively, which determine the basic feature of an LGA.Often, simulations generated using an LGA are very noisy due to its Boolean variable with one for the presence and zero for the absence of particles at a lattice point [4,5].Furthermore, the numerical procedure involves the calculations of particle probability, which reduces the efficiency of the model.To overcome these, the lattice Boltzmann method was proposed [6], and its basic difference from the LGA is that the Boolean variable is replaced with a particle distribution function.Such an approach eliminates the statistical noise in an LGA and retains all the advantages of locality in the kinetic form of an LGA [1].McNamara and Zanetti [6] first used the lattice Boltzmann method as an alternative to the LGA in 1988.As the collision operator takes a complex matrix form, this prevents the LBM from becoming a competing numerical method.A breakthrough in progress was made by Higuera and Jiménez [7] through linearisation of the collision term around its local equilibrium state.This greatly simplifies the collision operator.Noble et al. [8] used this idea to express the collision operator as Ω αβ ( f eq β − f β ), in which f β is the particle distribution 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 by using a single time relaxation towards the local equilibrium distribution, Ω αβ = −δ αβ /τ, which is popularly known as the Bhatnagar-Gross-Krook (BGK) collision operator [11], where δ αβ is the Kronecker delta function taking one when α = β, or otherwise zero, and τ is called the single relaxation time that takes a value.The correct choice of τ together with the lattice size and time step can achieve a desired fluid viscosity.This leads to the single relaxation time lattice Boltzmann method (SRT LBM), which has been the most efficient and used so far, where x j is a lattice coordinate along the j-axis in the Cartesian coordinate system, i.e., j = x, y in the two-dimensional space; t is time; e αj is the j th component of the particle velocity vector e α in α-link of the lattice and defined by time step δt and lattice size δx, e.g., e α = (0, 0), (e, 0), (0, e), (−e, 0), (0, −e), (e, e), (−e, e), (−e, −e), (e, −e) when α = 0 − 8 for nine particles moving in the two-dimensional uniform square lattice (D2Q9), in which e is the particle speed and defined as e = δx/δt; and f eq α is the local equilibrium distribution function given by: e αi e αj u i u j e 4 − 3 2 in which ρ is the fluid density and w α is a weighting factor depending on the lattice pattern, e.g., w α = 4/9 when α = 0, w α = 1/9 when α = 1 − 4 and w α = 1/36 when α = 5 − 8 on D2Q9.After the distribution function is calculated from the lattice Boltzmann Equation (1), the macroscopic psychical variables, density and velocity are simply updated as: Since then, the study of the lattice Boltzmann method and the applications of the method have received extensive attentions, making it become a very powerful modelling tool in many areas such as thermodynamics [12], aerodynamics [13], multiphase flows [14], turbulent flows [15], hemodynamics [16], biomechanics [17], image analysis [18], biology [19] and environmental science [20].
In applications, it is found that the SRT LBM suffers from a numerical instability.To remedy this, the multiple relaxation time (MRT) collision operator was introduced in 1992 [21,22].This improves the stability, but it reduces the efficiency.To accelerate the simulation, a two relaxation time collision operator was developed in 2008 [23], which has almost the same efficiency as the SRT LBM.As research progresses, it is noticed that MRT or TRT still suffers from numerical instabilities when fluid flows with very small viscosity are simulated.After realising that this is caused by an insufficient degree of Galilean invariance in the collision step, Geier et al. [24] proposed a cascaded lattice Boltzmann method by relaxing the particle distribution function to its local equilibrium state in the central moment space, making the LBM stable for simulating flows with small viscosity close to zero.In 2015, Geier et al. [25] further improved their central moment LBM using the cumulant in collision operator, which is called the cumulant lattice Boltzmann method (CLBM).Despite such enhancements, these schemes are more complicated and less efficient due to a manipulation of the matrices.In addition, these schemes share the same drawback as the existing LBM in that the boundary conditions for a physical variable such as velocity cannot be implemented without being converted to particle distribution functions, which further reduces the efficiency and accuracy of the methods.Recently, Chen et al. [26] developed a simplified lattice Boltzmann method without the evolution of the particle distribution function, which successfully removes this drawback and enables the direct use of a physical variable as the boundary conditions.However, the method involves the two steps of the predictor and corrector, which is more complicated than the SRT LBM.Nevertheless, through all this research, the LBM has greatly been improved and developed to a point where it has become a very efficient and flexible alternative numerical method in computational physics.Its potential power is much beyond the original scope, being explored and demonstrated in various disciplines of science and engineering with time [27][28][29][30].

Macroscopic Lattice Boltzmann Model
The above literature review highlights that the major research to improve the method has been carried out on the collision operator for resolving the well-known instability problem in the SRT LBM since 1992, leading to four representative variants of the method, MRT, TRT, central moment and cumulant LBMs.Even so, choosing suitable parameters or values for the collision operators 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 simulate a large-scale flow system, preventing the LBM from becoming an automatic simulator for any scale flows when a super-fast computer such as a quantum computer becomes available one day.As the central problem comes from the collision operator, such a problem may be resolved forever if the collision operator can be removed.Since the function of the collision is to relax the distribution function towards its local equilibrium state, it may be removed and replaced with the local equilibrium distribution function by setting τ = 1 in Equation (1).Following this idea, after mathematical manipulations, we obtain the following simple macroscopic lattice Boltzmann model (MacLAB), to determine the physical density and velocity directly from the local equilibrium distribution function without calculating the distribution function using Equation (1) unlike existing LBMs (the theoretical derivation is given in Section 3).Apparently, the model involves the local equilibrium distribution function only.However, doing so brings a new problem of how to consider the fluid viscosity in the absence of the collision step.We overcome this from the finding, through the recovery of the Navier-Stokes equations from Equations ( 4) and ( 5), that if the particle speed e is determined using: which is employed in Equation (2) instead of using e = δx/δt to calculate the local equilibrium distribution function f eq α , the flow viscosity is naturally taken into account in the model.In the MacLAB, once a 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 known quantity at the current time; and the particle speed e is determined from Equation ( 6) for use in the computation of f eq α .In addition, the time step δt is no longer an independent parameter and needs to be determined as δt = δx/e, which is used in the simulations of unsteady flows.Consequently, only the lattice size δx is required in the MacLAB for the simulation of fluid flows, bringing the LBM into a precise "lattice" Boltzmann method.This enables the model to become an automatic simulator for any scale flows without tuning other simulation parameters, making it possible and easy to model a large flow system when a super-fast computer such as a quantum computer becomes available in the future.The model is unconditionally stable as it shares the same valid condition as that for f eq α , or the Mach number M = U c /e is much smaller than one, in which U c is a characteristic flow speed.The Mach number can also be expressed as a lattice Reynolds number of R le = U c δx/ν via Equation (6).In practical simulations, it is found that the model is stable if R le = U m δx/ν < 1 where U m is the maximum flow speed and is used as the characteristic flow speed.The main features are that there is no collision operator and only macroscopic physical variables such as the density and velocity are required, which are directly retained as boundary conditions with a minimum memory requirement.The implementation of the model starting from the initial density and velocity is to (i) choose the lattice size δx and determine the particle speed e from Equation ( 6), (ii) calculate f eq α from Equation (2) using the density and velocity, (iii) update the density and velocity using Equations ( 4) and ( 5), (iv) apply the boundary conditions, and (v) repeat Step (ii) until a solution is reached.The only limitation of the described model is that, for very small viscosity or high speed flow, the chosen lattice size after satisfying R le < 1 may turn out to generate very large lattice points (lattice points, e.g., for one dimension with a length of L are calculated as N L = L/δx, and N L is the lattice points); if the total lattice points are too big such that the demanding computations are beyond the current power of a computer, the simulation cannot be carried out.Such difficulties may be solved or relaxed through parallel computing using computer techniques such as GPU processors and multiple servers and will be largely or completely removed using quantum computing when a quantum computer becomes available.

Recovery of the Navier-Stokes Equations
Setting τ = 1 in Equation (1) leads to: which can be rewritten as: Taking ∑ Equation ( 8) and ∑ e αi Equation ( 8) yields: ∑ e αi f α (x j , t) = ∑ e αi f eq α (x j − e αj δt, t − δt), (10) respectively.In the lattice Boltzmann method, the density and velocity are determined using the distribution function as: Combining Equation (11) with Equations ( 9) and ( 10) results in the current MacLAB, Equations ( 4) and (5).Since the local equilibrium distribution function f eq α has the features of: with reference to Equation ( 11), the following relationships, hold, which are the conditions that retain the conservation of the mass and momentum in the lattice Boltzmann method.
Next, we prove that the continuity and Navier-Stokes equations can be recovered from Equations ( 4) and (5).Rewriting Equation (1) as: Apparently, when τ = 1, the above equation becomes Equation (8), which leads to Equations ( 4) and ( 5); hence, Equation ( 14) is a general equation and is used in the following derivation.Applying a Taylor expansion to the two terms on the right-hand side of Equation ( 14) in time and space at point (x, t) yields: and: According to the Chapman-Enskog analysis, f α can be expanded around f (0) α : After substituting Equations ( 15)-( 17) into Equation ( 14), equating the coefficients of δt with various orders results for the order (δt) 0 in: for the order (δt) 1 in: and for the order (δt) 2 in: α . ( Substitution of Equation ( 19) into the above equation gives: which is rearranged as: From Equation (19) + Equation ( 22) ×δt, we have: Now, taking ∑Equation (23) leads to: as: due to the condition of the conservation of mass and momentum Equation (13).Evaluating the terms in the above equation using Equation (2) produces the exact continuity equation, Multiplying Equation ( 23) by e αi provides: Taking ∑ Equation ( 27) leads to: under the same condition (25) as that in the derivation of Equation (24).Evaluating the terms in the above equation using Equation (2) produces the exact momentum equation, the Navier-Stokes equation at second-order accuracy on the condition that the Mach number M = U c /e << 1, where pressure p is defined as: and the kinematic viscosity is: As τ takes a constant, the use of τ = 1 will recover the continuity and the Navier-Stokes equations at the second-order accuracy, as the above derivation shows.In this case, Equation ( 31) becomes Equation ( 6), which determines the particle speed e in the proposed model.

Couette Flow
The first test is a Couette flow through two parallel plates without a pressure gradient.The distance between the plates is h = 1.The top plate moves at a velocity of u x = u 0 = 0.1 in the streamwise direction, and the bottom plate is fixed.If x stands for the streamwise direction and y for the vertical direction, the analytical solution is: which is the same test as that used by Chen et al. [26].This is a very interesting case as the steady flow is independent of the flow viscosity according to the theory (32).We use δx = 0.02 and 20 × 50 lattices in the x and y directions for three simulations of flows with three kinematic viscosities of ν 1 = 0.01, ν 2 = 0.001 and ν 3 = 0.0006, respectively.The periodic boundary conditions are applied at the inflow and outflow boundaries.After the steady solutions are reached, the results are indeed independent of the viscosities, and one of those is shown in Figure 1, demonstrating excellent agreement with the analytical solution.

Couette Flow with a Pressure Gradient
The second test is the same flow as that in Test 1 except that a pressure gradient of ∂p/∂x = −0.0001 is specified, which is added to the right-hand side of Equation ( 5) as +δx/(eρ)∂p/∂x [31].Both plates are fixed with zero velocities at the top and bottom boundaries at which no calculations are needed.The flow is affected by viscosity, and the analytical solution is: The flow is simulated using three viscosities of ν 1 = 0.003, ν 2 = 0.001 and ν 3 = 0.0006.The numerical results are plotted in Figure 2, showing the effect of viscosity on the flow, in excellent agreement with the analytical solutions.This confirms the unique feature that the current model can simulate viscous flow correctly due to the use of Equation ( 6), although no explicit effect of viscosity on flows is taken into account.

2D Cavity Flow
The third test is a 2D cavity flow, which is a well-known complex flow within a simple geometry.The domain is a 1 × 1 square.The boundary conditions are that the top lid moves at a velocity of u x = u 0 and u y = 0 with u 0 = 1; the other three sides are fixed, or no slip boundary condition is applied, i.e., u x = 0 and u y = 0.The Reynolds number R e = u 0 /ν = 1000.We use δx = 0.0025 or 400 × 400 lattices in the simulation, which is carried out on the inside of the cavity excluding the four sides where the velocities are retained as boundary conditions.After the steady solution is obtained, the flow pattern in the velocity vectors is shown in Figure 3, which closely agrees with the well-known study by Ghia et al. [32].The results are further compared against their numerical solutions for velocity profiles of u x and u y along the y and x directions through the geometric centre of the cavity in Figures 4 and 5, respectively, demonstrating very good agreements.After the steady solution is obtained, the comparison of the velocity u x profile along the y direction through the geometric centre of the cavity with the numerical solution by Ghia et al. [32], showing good agreement.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 the velocity u y profile along the x direction through the geometric centre of the cavity with the numerical solution by Ghia et al. [32], showing good agreement.

2D Taylor-Green Vortex
The fourth test is a 2D Taylor-Green vortex.This is an unsteady flow driven by decaying vortexes for which there is an exact solution of the incompressible Navier-Stokes equations, and it is often applied to the validation of a numerical method for the solution to the incompressible Navier-Stokes equations.The initial conditions are u x (x, y, 0) = −u 0 cos(x) sin(y) and u y (x, y, 0) = u 0 sin(x) cos(y).The analytical solutions are u x (x, y, t) = −u 0 cos(x) sin(y) exp(−2νt) and u y (x, y, t) = u 0 sin(x) cos(y) exp(−2νt).The time for an unsteady flow from the initial state is accumulated by its increase with time step δt.We use δx = 0.157 or 40 × 40 lattices for a square domain of 2π × 2π with a kinematic viscosity of ν = 0.0314 and u 0 = 0.05, which gives the Reynolds number of R e = 2πu 0 /ν = 10.The periodic boundary conditions are used.The simulation is run for the total time of 30 s.The velocity field is plotted in Figure 6, showing the correct flow pattern.The velocity profiles for u x at x = π and u y at x = π/2 along the y direction are depicted and compared with the analytical solutions in Figure 7, showing excellent agreements and confirming the accuracy of the method for an unsteady flow.Taylor-Green vortex within the 2π × 2π domain for R e = 10.The initial conditions are u x (x, y, 0) = −u 0 cos(x) sin(y) and u y (x, y, 0) = u 0 sin(x) cos(y) with u 0 = 0.05.The periodic boundary conditions are used.Shown here is the comparisons of the relative velocity profiles for u x /u 0 at x = π and u y /u 0 at x = π/2, matching the analytical solutions of u x (x, y, t) = −u 0 cos(x) sin(y) exp(−2νt) and u y (x, y, t) = u 0 sin(x) cos(y) exp(−2νt) at t = 30 s.

3D Cavity Flow
The final test is a 3D cavity flow.This is again a well-known complex flow involving 3D vortices within a simple cube with dimensions of 1 × 1 × 1 in streamwise direction x, spanwise direction y and vertical direction z.No-slip boundary conditions, i.e, u x = 0, u y = 0 and u z = 0, are applied to five fixed sides except for the top lid, where u x = u 0 , u y = 0 and u z = 0 with u 0 = 1 are specified.The Reynolds number is R e = u 0 /ν = 400.δx = 0.004 or total lattices of 250 × 250 × 250 are used, and the simulation is undertaken only within the cube, excluding the boundaries where the velocities are retained.After the steady solution is reached, the flow characterises are displayed through the two-dimensional planar projections of the velocity vector field on the x-z, y-z and x-y centroidal planes of the cube in Figures 8-10, respectively, demonstrating flow patterns in good agreement with those by Wong and Baker [33].In addition, the distribution of the velocity component u x on the vertical plane centerline is widely used as a 3D lid-driven cavity benchmark test.We compare this velocity component against the results by Wong and Baker [33] and also by Jiang et al. [34] in Figure 11, showing good agreements.After the steady solution is reached, the comparisons of the distribution of the velocity component u x on the vertical plane centerline with the results by Wong and Baker [33] and Jiang et al. [34].

Conclusions
A novel macroscopic lattice Boltzmann method (MacLAB) for the Navier-Stokes equations is proposed and validated.It can accurately simulate fluid flows using only the lattice size, bringing the LBM into a precise lattice Boltzmann method and revolutionising the standard LBM theory.This takes the research on the method into a new era in which future work may focus on improving the accuracy of or formulating a new local equilibrium distribution function.The particle speed is determined through the viscosity and lattice size, and the time step δt is calculated as δt = δx/e.The model is unconditional stable as long as the valid condition for the local equilibrium distribution function holds.When a super-fast computer such as a quantum computer becomes available in the future, this will enable the MacLAB to take any fine lattice to simulate most challenging fluid flows such as turbulent flows without relying on a turbulent flow model.All these make the method an automatic simulator for fluid flows in all relevant problems.The method is straightforward to extend for resolving other physical problems in different disciplines such as chemical and environmental engineering.

Figure 1 .Figure 2 .
Figure 1.Couette flow through two parallel plates without a pressure gradient.The distance between the plates is h = 1.The top plate moves at a velocity of 0.1 in the streamwise direction, and the bottom plate is fixed, where no calculations are required at the top or bottom boundaries.The periodic boundary conditions are applied at the inflow and outflow boundaries.δx = 0.02 is used for three simulations of flows with three kinematic viscosities of ν 1 = 0.01, ν 2 = 0.001 and ν 3 = 0.0006, respectively.All the steady numerical results are almost identical and are independent of flow viscosity, as shown here in the comparison of one numerical result with the analytical solution.MacLAB, macroscopic lattice Boltzmann method.

Figure 3 .
Figure 3. 2D cavity flow within a 1 × 1 square for R e = 1000.The top lid moves at a 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 vortex and two secondary vortices.

Figure 4 .
Figure 4. 2D cavity flow within a 1 × 1 square for R e = 1000.The top lid moves at a 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 the velocity u x profile along the y direction through the geometric centre of the cavity with the numerical solution by Ghia et al.[32], showing good agreement.

Figure 5 .
Figure5.2D cavity flow within a 1 × 1 square for 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 the velocity u y profile along the x direction through the geometric centre of the cavity with the numerical solution by Ghia et al.[32], showing good agreement.

Figure 6 .
Figure 6.Taylor-Green vortex within the 2π × 2π domain for R e = 10.The initial conditions are u x (x, y, 0) = −u 0 cos(x) sin(y) and u y (x, y, 0) = u 0 sin(x) cos(y) with u 0 = 0.05.The periodic boundary conditions are used.Shown here is the flow pattern in velocity vectors at t = 30 s, the same vortex pattern remaining as that at the initial state.
Figure 7.Taylor-Green vortex within the 2π × 2π domain for R e = 10.The initial conditions are u x (x, y, 0) = −u 0 cos(x) sin(y) and u y (x, y, 0) = u 0 sin(x) cos(y) with u 0 = 0.05.The periodic boundary conditions are used.Shown here is the comparisons of the relative velocity profiles for u x /u 0 at x = π and u y /u 0 at x = π/2, matching the analytical solutions of u x (x, y, t) = −u 0 cos(x) sin(y) exp(−2νt) and u y (x, y, t) = u 0 sin(x) cos(y) exp(−2νt) at t = 30 s.

Figure 8 .
Figure8.3D cavity flow within a 1 × 1 × 1 cube for R e = 400.The top lid moves at a velocity of u x = 1, u y = 0 and u z = 0, and the other five sides are fixed, or no slip boundary condition is applied.After the steady solution is reached, the flow pattern in vectors in the x − zcentroidal plane show the primary and secondary vortices.

Figure 9 .Figure 10 .Figure 11 .
Figure9.3D cavity flow within a 1 × 1 × 1 cube for R e = 400.The top lid moves at a velocity of u x = 1, u y 0 and u z = 0, and the other five sides are fixed, or no slip boundary condition is applied.After the steady solution is reached, the flow pattern in vectors in the y − z centroidal plane shows one pair of strong secondary vortices at the bottom and one pair of weak secondary vortices at the top.