Numerical extraction of a macroscopic pde and a lifting operator from a lattice Boltzmann model
NNumerical extraction of a macroscopic PDE and a liftingoperator from a Lattice Boltzmann model
Ynte Vanderhoydonc ∗ Wim Vanroose † Abstract
Lifting operators play an important role in starting a lattice Boltzmann model from a giveninitial density. The density, a macroscopic variable, needs to be mapped to the distributionfunctions, mesoscopic variables, of the lattice Boltzmann model. Several methods proposedas lifting operators have been tested and discussed in the literature. The most famous meth-ods are an analytically found lifting operator, like the Chapman-Enskog expansion, and anumerical method, like the Constrained Runs algorithm, to arrive at an implicit expressionfor the unknown distribution functions with the help of the density. This paper proposes alifting operator that alleviates several drawbacks of these existing methods. In particular,we focus on the computational expense and the analytical work that needs to be done. Theproposed lifting operator, a numerical Chapman-Enskog expansion, obtains the coefficients ofthe Chapman-Enskog expansion numerically. Another important feature of the use of liftingoperators is found in hybrid models. There the lattice Boltzmann model is spatially coupledwith a model based on a more macroscopic description, for example an advection-diffusion-reaction equation. In one part of the domain, the lattice Boltzmann model is used, while inanother part, the more macroscopic model. Such a hybrid coupling results in missing data atthe interfaces between the different models. A lifting operator is then an important tool sincethe lattice Boltzmann model is typically described by more variables than a model based ona macroscopic partial differential equation.
Keywords: Lifting operator, missing data, lattice Boltzmann models, macroscopic partial dif-ferential equations, hybrid models, Chapman-Enskog expansion, Constrained Runs, numericalChapman-Enskog expansion.
A lifting operator is, in a multiscale method, an important tool that maps macroscopic variables tomicroscopic/mesoscopic variables. In kinetic models, for example, a lifting operator will map loworder moments, like the density ρ ( x , t ) that counts the number of particles in a point x ∈ D x ⊂ R n , n ∈ N to a distribution function f ( x , v , t ) that counts the number of particles in a point ( x , v )in phase space where x ∈ D x ⊂ R n and the velocity v ∈ D v ⊂ R n . For practical applications oneuses n ∈ { , , } and time t ≥ ρ ( x , t ) can be represented by an advection-diffusion-reactionequation. While the microscopic/mesoscopic level is typically described by a Boltzmann equationthat evolves the distribution function f ( x , v , t ). A lattice Boltzmann model (LBM) is a specialdiscretization of the Boltzmann equation that is, for example, used to simulate complex fluidsystems. ∗ Dept. Mathematics and Computer Science, Universiteit Antwerpen, Middelheimlaan 1, 2020 Antwerpen, Bel-gium ( [email protected] ). † Dept. Mathematics and Computer Science, Universiteit Antwerpen, Middelheimlaan 1, 2020 Antwerpen, Bel-gium ( [email protected] ). a r X i v : . [ c s . C E ] F e b ome examples of complex flows for which LBMs are used are flows in complicated geometries,multiphase and turbulent flows. Applications can be found in [3, 25] and a more recent review isgiven by Aidun et al. [1]. The application of the lattice Boltzmann model to multiscale physics influids is discussed in [26]. Banda et al. [2] present a high order relaxation system for the multiscalelattice Boltzmann equation to obtain the incompressible Navier-Stokes limit.Kevrekidis et al. introduced a lifting operator to couple different scales in a dynamical systemin the equation-free framework [17]. This allows to model the dynamics at the macroscopic levelby using short bursts of the microscopic simulation.A problem of lattice Boltzmann methods is the determination of initial conditions, usuallygiven by macroscopic variables. During initialization and spatial coupling in a hybrid model, aone-to-many map needs to be created, known as a lifting operator. In this article we discuss ahybrid LBM and PDE model that uses a lattice Boltzmann model in one part of the domain whileanother part of the domain is described by a macroscopic partial differential equation. Thesedifferent levels of description create missing data at the interfaces that can be resolved with alifting operator.Hybrid approaches have been formulated for various flow problems. A Lennard-Jones particledynamics is coupled with a compressible Navier-Stokes in [12]. A Boltzmann, respectively a latticeBoltzmann, model is coupled to the Navier-Stokes equations in [19] and [18]. A LBM is also coupledwith a Navier-Stokes in Peano, an adaptive mesh refinement framework with spacetree grids [20].Furthermore, the Boltzmann equation is coupled to the Euler equations in [4].Coupled models also play an important role in the simulation of materials. A review ofatomistic-to-continuum coupling is found in [22]. A more detailed overview for coupling methodsin hybrid models can be found in [13] and [11]. In [13] adaptive mesh and algorithm refinementis used in parts of the domain where a continuum description is replaced by a particle descrip-tion. Coupled molecular dynamics and lattice Boltzmann models based on Schwarz’s alternatingmethod is presented in [11]. Dimarco et al. [10] merge deterministic methods for the equilibriumpart with particle methods for the nonequilibrium part and present results for the Boltzmannequation with Bhatnagar-Gross-Krook (BGK) approximation.The coupling — that will be discussed in this article — of LBMs with reaction-diffusion PDEsis earlier considered in [29, 28, 33].We propose a general lifting operator that maps densities to distribution functions. It is illus-trated for a LBM, but we believe that it is applicable to general discretizations of the Boltzmannequation and it can map more moments to the corresponding distribution functions.The new method will be compared to the Chapman-Enskog expansion [6], a well known an-alytical method for the initialization of a lattice Boltzmann model, and the Constrained Runs(CR) algorithm [30], a numerical lifting operator. The Chapman-Enskog expansion writes thedistribution functions as an analytical series of the density. The Constrained Runs algorithm isbased on the attraction of the dynamics toward the slow manifold and expresses, in an implicitway, the unknown distribution functions with the help of the density in successive grid points. Anumerical comparison of these methods is given in [33] for hybrid models that spatially couple adiffusion PDE model and a LBM. Although the Constrained Runs algorithm is very accurate, itsmajor drawback is the computational expense. It achieves a high accuracy for the coupling of alattice Boltzmann model with a diffusion-reaction PDE [29]. However, lifting in the CR-algorithmrequires many additional LBM steps. This computational cost is too expensive to be useful in morecomplex problems in higher dimensions. These drawbacks became clear in [33] when comparingthe different methods numerically. The intention of this paper is to alleviate them.In this paper we propose a numerical Chapman-Enskog expansion that seriously reduces thecomputational cost of the lifting. It combines the idea of the Constrained Runs algorithm with theChapman-Enskog expansion and does not need an analytical derivation as the Chapman-Enskogexpansion. This lifting operator is calculated before the simulation and finds the coefficients ofthe Chapman-Enskog expansion numerically. Once these coefficients are found the applicationof the lifting operator is just a stencil computation as cheap as the analytical Chapman-Enskogexpansion. As a spin-off it also extracts the macroscopic PDE from the lattice Boltzmann model.This allows us to construct the hybrid model without deriving the macroscopic PDE analytically.2he numerical results show that the new lifting operator can also reach a high accuracy. Although,we illustrate and benchmark the new method on academic model problems, we believe that it isapplicable to other discretizations of the Boltzmann equation. Furthermore, for the clarity of thepresentation we have kept the boundary between the LBM and PDE domain fixed. In a realapplication this boundary might be moved adaptively, triggered by an error estimate similar as inadaptive mesh refinement.This work is organized as follows. In Section 2 the model problem is defined. It focuses, inparticular, on a hybrid model that consists of a LBM in one part of the domain and a macroscopicequivalent PDE in another part. Section 3 gives an overview of existing lifting techniques thatare used in the literature. The Chapman-Enskog expansion and the Constrained Runs algorithmare respectively considered in Sections 3.1 and 3.2. These methods are discussed in Section 3.3.We tend to remove these drawbacks by considering a numerical Chapman-Enskog expansion inSection 4. Section 5 contains the numerical results. In Section 5.1 the proposed lifting operator istested in a setting of restriction and lifting. The application of the lifting operator to the hybridLBM and PDE model is considered in Sections 5.2 and 5.3. We conclude and give an outlook inSection 6. Kinetic models make use of the Boltzmann equation [25] that describes the evolution of a distribu-tion function f ( x , v , t ) (function space C R ( D )) that counts the number of particles or individualsin point x ∈ D x ⊂ R n , n ∈ N , with a velocity v ∈ D v ⊂ R n , at time t ≥
0. The equation is ∂∂t f ( x , v , t ) + v ∂∂ x f ( x , v , t ) + F ( x , t ) ∂∂ v f ( x , v , t ) = Ω . (1)This is an evolution law in phase space where F ( x , t ) is the external force and Ω an integraloperator that models the reorganization of the velocity distribution due to collisions or otherinteractions.The collision operator can be approximated by a simpler Bhatnagar-Gross-Krook (BGK) modelΩ = ω ( f eq ( x , v , t ) − f ( x , v , t )) [15] in which the equilibrium distribution f eq ( x , v , t ) is given by theMaxwell-Boltzmann distribution [25]. The BGK approximation represents a relaxation towardsequilibrium with an associated time scale τ = 1 /ω .At the moment it is still computationally expensive to simulate or analyze a Boltzmann modelnumerically. The development of efficient numerical methods for models based on the Boltzmannequation is therefore an active research field. A possible increase in efficiency can be obtained byconstructing a hybrid model. The kinetic model is then replaced with a macroscopic descriptionin the regions of the spatial domain where this is justified, for example, away from reaction fronts.These macroscopic models are cheaper to simulate.Possible macroscopic models for fluid dynamics are described by the Navier-Stokes and Eulerequations. One can derive both the Navier-Stokes as the Euler equations from the Boltzmannequation [34]. A lifting operator then transforms the variables of the PDE to the distributionfunction of the Boltzmann equation at the boundaries between the domains.In this paper we study a simple model that allows a detailed study of the lifting operator bothfor the initialization of the distribution function and in a hybrid context. The model uses a latticeBoltzmann model with an equilibrium distribution function that only depends on the density.While the macroscopic PDE is an advection-diffusion-reaction equation for the density only. Thissimple model allows a detailed analysis yet it is general enough to expect that the results can beextended to more realistic Boltzmann models. The correspondence of the lattice Boltzmann withthe Boltzmann equation is discussed in [34, 24]. A lattice Boltzmann model (LBM) [25, 34] is a special discretization of Eq. (1). It describes theevolution of one-particle distribution functions f i ( x , t ) = f ( x , v i , t ) discretized in space x , time3 and velocity v i . The velocities are taken from a discrete set defined by the geometry of thegrid. The functions are represented as f i : X × T → R with X × T the space-time grid withspace steps ∆ x i in the direction of velocity v i , time step ∆ t and T = { , ∆ t, t, . . . } ∩ [0 , T ].Representation DdQq used for the description of LBMs stands for d dimensions and q velocitydirections. D1Q3, for example, considers in a one-dimensional spatial domain only three valuesfor the velocity v i = c i ∆ x/ ∆ t with c i = i , i ∈ {− , , } the dimensionless grid velocities.The remaining of this section contains the description of the lattice Boltzmann equation in onedimension but can easily be generalized to more dimensions.The lattice Boltzmann equation (LBE) describing the evolution of the distribution functions(with BGK approximation and no external force in Eq. (1)) is f i ( x + c i ∆ x, t + ∆ t ) = (1 − ω ) f i ( x, t ) + ωf eqi ( x, t ) . (2)The equilibrium distributions are given by f eqi ( x, t ) = ρ ( x, t ), i ∈ {− , , } [27] in which theparticle density ρ ( x, t ) is defined as the zeroth order moment of the distribution functions ρ ( x, t ) = (cid:80) i ∈{− , , } f i ( x, t ). These equilibrium distributions correspond to a local diffusive equilibrium.The focus of this paper concerns the initialization of a lattice Boltzmann model. Starting theLBM scheme from a given initial density includes some arbitrariness. The distribution functions f i ( x,
0) at time t = 0 need to be constructed from a given density ρ ( x, This section contains descriptions that represent macroscopic equivalent PDEs specific for LBMs.Partial differential equations model, at a macroscopic scale, the evolution of the moments of theparticle distribution functions like density ρ ( x, t ) = (cid:80) i f i ( x, t ), momentum φ ( x, t ) = (cid:80) i v i f i ( x, t )or energy ξ ( x, t ) = (cid:80) i v i f i ( x, t ).The transition between the distribution functions and the moments is straightforward sincethe matrix M below is invertible. ρφξ = − f f f − = M f f f − . (3)When we look at these functions in a point x at time t , they can be represented either as( f , f , f − ) T ∈ R or as ( ρ, φ, ξ ) T ∈ R . If we focus on the complete discretization in space,with n spatial grid points, the function spaces are R × n .It can be shown that the diffusion PDE and the LBM are macroscopic equivalent [28] whenconsidering D1Q3 and ∂ρ∂t = D ∂ ρ∂x , D = 2 − ω ω ∆ x ∆ t , f eqi ( x, t ) = 13 ρ ( x, t ) , i ∈ {− , , } . (4)This can be checked by using a Chapman-Enskog expansion. Here f i ( x, t ), i ∈ {− , , } is writtenas a series, each term containing higher order derivatives of ρ ( x, t ) [6, 28, 5]. f i = f [0] i + f [1] i ∆ x + f [2] i ∆ x + f [3] i ∆ x + . . . , (5)where f [0] i = f eqi = ρ , f [1] i = − i ω ∂ρ∂x , f [2] i = − ω ( ω − i − ∂ ρ∂x . The macroscopic diffusion PDE (4) is obtained by summing the series of the Chapman-Enskogexpansion in (5) over the velocities. This considers purely diffusive effects.For advection-diffusion problems with uniform velocity field a on D1Q3 the equilibrium distri-bution functions are given by [27] f eqi ( x, t ) = 13 (cid:18) icac s + a c s (cid:19) , i ∈ {− , , } , c s = 23 c , c = ∆ x ∆ t , (6)4ith the equivalent macroscopic description ∂ρ∂t + a ∂ρ∂x = D ∂ ρ∂x , D = 2 − ω ω ∆ x ∆ t . (7)Similar results can be obtained for higher dimensional problems. In two spatial dimensions,represented by x and y , with equal space steps ∆ x = ∆ y , the macroscopic equivalence is given by ∂ρ∂t = D ∂ ρ∂x + D ∂ ρ∂y , D = 2 − ω ω ∆ x ∆ t , f eqi ( x , t ) = w i ρ ( x , t ) , i ∈ { , . . . , } , x = (cid:18) xy (cid:19) , w = 13 , w = . . . = w = 16 , (8)for D2Q5 [27] and ∂ρ∂t + a x ∂ρ∂x + a y ∂ρ∂y = D ∂ ρ∂x + D ∂ ρ∂y , D = 2 − ω ω ∆ x ∆ t ,f eqi ( x , t ) = w i ρ (cid:18) v i · a c s + ( v i · a ) c s − a c s (cid:19) , i ∈ { , . . . , } , a = ( a x , a y ) , c s = 13 ∆ x ∆ t , x = (cid:18) xy (cid:19) , w = 49 , w i = 19 , i ∈ { , . . . , } , w i = 136 , i ∈ { , . . . , } , (9)for D2Q9 [27].With such analytical expressions from the Chapman-Enskog expansion available, a lifting op-erator can be constructed. Indeed, f i ( x, t ) is then written as a series in function of the givendensity ρ ( x, t ). This allows us to construct the distribution functions from a given initial densitynecessary to initialize the LBM.We will consider the lattice Boltzmann model as the ‘exact’ model and these PDEs of thedensity as the macroscopic approximations. From this, we can construct a hybrid model problemas outlined in Section 2.3. This section deals with the construction of the hybrid model problem. Consider the one-dimensionalproblem D1Q3 but bear in mind that a similar construction of a hybrid model can be done in higherdimensions. In particular, we couple a lattice Boltzmann model with the macroscopic equivalentPDE. The resulting hybrid domain for D1Q3 is presented in Figure 1. A one-dimensional domain[ a, b ] is considered that couples the PDE (4) on [ a, l [ with the LBE (2) on [ l, b ]. Furthermore, weassume periodic boundary conditions. Note that in both subdomains the same grid spacings areused in space (∆ x ) and in time (∆ t ) and the boundary l remains fixed for all times. Using adifferent space-time grid is of particular interest for future work but the same spacings are usedto highlight the coupling error. Similarly, the boundary may be moved adaptively as in adaptivemesh refinement. In an actual physical problem it will be important to use the hybrid model thatgives a Boltzmann description where shock waves, contact discontinuities or sharp gradients occur.Since these move in time it might be useful to work with a moving interface method as in [8] and[9]. However, this is not the focus of the current paper.On the domain [ a, l [, we discretize the PDE (4) with cell centered central differences in spaceand forward Euler time discretization. The grid points x j with j ∈ { , , . . . , p } cover this domainand for these points it holds that ρ ( x j , t + ∆ t ) = ρ ( x j , t ) + D ∆ t ∆ x ( ρ ( x j − , t ) − ρ ( x j , t ) + ρ ( x j +1 , t )) . (10)For the grid points x j with j ∈ { p + 1 , . . . , n − } in [ l, b ] the LBE (2) holds.The full domain has an initial condition ρ ( x j , ∀ j ∈ { , , . . . , n − } . The lifting operator isrequired to formulate the initial conditions of the LBM domain f i ( x j , ∀ j ∈ { p + 1 , . . . , n − } .5 l bx − = x n − x x . . . . . .x p x p +1 x n − x n = x ρ ( x , t ) ρ ( x p , t ) f − f f PDE domain LBM domain
Figure 1: The domain [ a, b ] in the hybrid model is split into [ a, l [ on which we solve the PDEmodel and [ l, b ] on which we solve the LBM. The solid points ( • ) represent the grid for the density ρ of the discrete PDE, the circles ( ◦ ) represent the LBM variables ( f , f , f − ) T . The periodicboundary conditions and the coupling are implemented with ghostcells which are drawn by dashedcircles. The density in the ghostcells of the PDE domain, in x − and x p +1 , are found by taking (cid:80) i f i in x n − and x p +1 , respectively. However, the ghostcells for the LBM domain, in x p and x n ,require a lifting operator that lifts ρ to ( f , f , f − ) T in these points.The periodic boundary conditions lead to the following boundary conditions for the PDEdomain: ∀ t : ρ ( x − , t ) = (cid:80) i f i ( x n − , t ) and ρ ( x p +1 , t ) = (cid:80) i f i ( x p +1 , t ).The aim is to construct the boundary conditions of the LBM domain in such a way that ∀ t > ∀ j ∈ { , , . . . , n − } the macroscopic density defined as ρ ( x j , t ) = (cid:40) ρ ( x j , t ) if j ∈ { , . . . , p } , (cid:80) i f i ( x j , t ) if j ∈ { p + 1 , . . . , n − } , (11)behaves as the density of a LBM solved on the full domain.To formulate these boundary conditions, a lifting operator is required that maps the density ρ ( x, t ) in the ghost points x and x p , the unknown of the PDE, to the distribution functions f i ( x, t ), i ∈ {− , , } of the LBM.This can be generalized by considering a higher dimensional spatial domain. The remainingderivations in this paper focus on one dimension although they can also be generalized to moredimensions. This section gives an overview of existing lifting operators that map densities to distributionfunctions. An analytical expansion that expresses the distribution functions as a series of thedensity and its spatial derivatives is given in Section 3.1, while a numerical method is presentedin Section 3.2. Section 3.3 discusses these methods which results in the motivation to propose alifting operator based on the combined ideas of the analytical and numerical method.
The Chapman-Enskog expansion [6, 28, 5], already discussed in Section 2.3, can be used as alifting operator. It constructs a mapping from ρ ( x, t ) to f i ( x, t ), i ∈ {− , , } (D1Q3). An alternative numerical procedure is the Constrained Runs algorithm discussed in this section. Itis well known that in phase space the dynamics are quickly attracted toward a slow manifold [14].For the problem we are studying the dynamics on the slow manifold can be parameterized by thedensity ρ ( x, t ). The distribution functions are then of the form { f ( ρ ( x, t )) , f ( ρ ( x, t )) , f − ( ρ ( x, t )) } .6ecause of Eq. (3), it is equivalent to determine { f ( ρ ) , f ( ρ ) , f − ( ρ ) } or { ρ, φ ( ρ ) , ξ ( ρ ) } . Themissing distribution functions { f , f , f − } can be found by determining φ and ξ for a given ρ suchthat φ , ξ and ρ lie on the slow manifold. This is the basic idea of the Constrained Runs (CR)algorithm that was proposed by Gear et al. [14] for stiff singularly perturbed ordinary differentialequations (ODEs) to map macroscopic initial data to missing microscopic variables. It uses thenumerical simulator to find the missing data such that the evolution is close to the slow manifold.This Constrained Runs algorithm can be applied to lattice Boltzmann models [30]. The stateof the lattice Boltzmann model can be split into u = ( ρ ) and v = (cid:18) φξ (cid:19) , where u ∈ R n and v ∈ R n for a LBM with n spatial grid points. The density ρ is known so u is given, while v is unknown since φ and ξ are missing. Denote the known initial conditions as u (0) = u ∈ R n .The idea is now to initialize v such that the evolution of v under the LBM is smooth of order m . The smoothness condition is defined by d m +1 v ( t ) dt m +1 (cid:12)(cid:12)(cid:12)(cid:12) t =0 = 0 , which is approximated by ∆ m +1 v ( t ) ≈ ∆ t m +1 d m +1 v ( t )d t m +1 , (12)where ∆ m is the well-known forward finite difference stencil on v ( t ), v ( t + ∆ t ), . . . . For m = 0, theconverged v satisfies the smoothness condition, up to a certain tolerance, and it is an approximationto the point of intersection with the slow manifold. This is schematically represented in Figure2. This iteration is always stable and the point of intersection is found to first order accuracycompared to the Chapman-Enskog expansion for the LBM with BGK collisions for one-dimensionalreaction-diffusion problems [30]. For m ≥
1, multiple LBM steps are necessary to estimate thederivative. For m = 1, this is often interpreted as a backward linear extrapolation in time [32].The number of LBM steps used in the backward extrapolation determines the accuracy of thescheme. Higher order schemes increase the accuracy but they can become unstable. In [32] thisinstability is circumvented by formulating the point of intersection as a fixed point v k +1 = C m ( u , v k ) , (13)where C m denotes one step of the CR-algorithm and m is related to the order of the time derivativethat is set to zero in the backward extrapolation in time. In general Eq. (13) is nonlinear and thefixed point can be found by a Newton-Krylov iteration. However, this requires many additionalLBM evaluations to construct the Jacobian. Similarly, matrix-free methods like GMRES stillrequire many matrix-vector products since the spectrum is unfavorable for fast convergence [31].In [33] the CR-algorithm is combined with Newton’s method by performing local updates at theghost points of the hybrid model to reduce the size of the Jacobian. The methods discussed in Sections 3.1 and 3.2 are well known methods to construct a liftingoperator for LBMs. However, each of these methods has some drawbacks. As noted earlier, adrawback of the use of the Chapman-Enskog expansion (Section 3.1) is the necessity to constructthe expressions analytically. Therefore its use is limited to a few examples where the expansion isknown. However, its computational cost is limited to the calculation of the numerical approxima-tion of the derivatives. The Chapman-Enskog expansion then becomes a stencil operator and thecost of the application grows linearly with the number of points where lifting is required.7igure 2: Sketch of the first few steps of the Constrained Runs algorithm for the LBM with a con-stant backward extrapolation in time. The solid line shows the evolution of { ρ ( x, t ) , φ ( x, t ) , ξ ( x, t ) } along the slow manifold for a given grid point x . For a given ρ we search for the intersectionof the plane with the slow manifold. We start iterating with ρ , the known density, and initialguesses φ and ξ for the missing moments ( v ). After each step of the LBM, the density is resetto its initial value ρ but the moments evolve during the LBM time simulation. This results in v k , the k -th iterate of the CR-algorithm. This algorithm finds an approximation for the missingvalues φ and ξ on the slow manifold.The Constrained Runs scheme (Section 3.2) can be used to approximate these expressionsnumerically. However, the lifting method can become computationally expensive since it requiresmany evaluations of the underlying lattice Boltzmann model to construct the Jacobian matrix.Even with the matrix-free methods and local updates discussed at the end of Section 3.2 it stillremains computationally expensive to use in practice, especially in higher dimensional problems.As an advantage of the Constrained Runs algorithm, we should note that the lifting errorcan be smaller than the modeling error, the difference in density between the LBM and its PDEapproximation, by using the CR-algorithm in the hybrid model discussed in Section 2.3 [29].Section 5.4 contains a comparison of the computational cost of these existing methods in the senseof hybrid models.The focus of this paper is to obtain an alternative lifting operator that reduces the computa-tional cost but holds the advantage of achieving the modeling error. In this section, we construct a lifting operator that alleviates the computational expense of theCR-algorithm. It combines the ideas of Constrained Runs and the Chapman-Enskog expansion.Instead of using Constrained Runs to find for each grid point the missing moments φ and ξ of thedistribution functions, we use Constrained Runs to find the unknown coefficients of the Chapman-Enskog expansion. This has several advantages that will be discussed at the end of Section 4.5.The derivations in this section are again based on one-dimensional problems but can easily begeneralized to more dimensions. 8 .1 Distribution functions as a series of the density This section shows that the solution f i ( x, t ), i ∈ {− , , } of a LBM with an infinite domain andparameters ∆ x , ∆ t and ω can be written as a series of ρ ( x, t ), the macroscopic density. We ini-tially characterize distribution functions f i as smooth functions that are sufficiently differentiablefunctions in time and space which implies that the same holds for the density, a sum of thesedistribution functions. The smoothness condition will be specified below. This condition can bejustified when the lattice spacing ∆ x is much bigger than the mean free path [16, 7]. Then thedistribution functions can be written as f i ( x, t ) = f eqi ( x, t ) + α i ∂ρ∂x + β i ∂ ρ∂x + δ i ∂ ρ∂x + (cid:15) i ∂ ρ∂x + . . . + γ i ∂ρ∂t + ζ i ∂ ρ∂t + . . . + η i ∂ ρ∂x∂t + . . . , (14)where α = α α α − ∈ R , β = β β β − ∈ R , . . . , (15)are fixed constants that only depend on ω , ∆ x and ∆ t . The derivation of this expansion is outlinedin the remaining of this section.Since the functions f i ( x, t ) are infinitely differentiable, a Taylor expansion can be constructed.The distribution functions in point x + i ∆ x , i ∈ {− , , } at time t + ∆ t are given by f i ( x + i ∆ x, t + ∆ t ) = f i ( x, t ) + ∂f i ∂x i ∆ x + ∂ f i ∂x i ∆ x ∂f i ∂t ∆ t + ∂ f i ∂t ∆ t ∂ f i ∂x∂t i ∆ x ∆ t + . . . . Combined with the assumption that f i is a solution of the LBE (2) on an infinite domain, we endup with f i ( x, t ) = f eqi ( x, t ) − i ∆ xω ∂f i ∂x − i ∆ x ω ∂ f i ∂x − ∆ tω ∂f i ∂t − ∆ t ω ∂ f i ∂t − i ∆ x ∆ tω ∂ f i ∂x∂t − . . . . With the notation L i for the functional L i = − i ∆ xω ∂∂x − i ∆ x ω ∂ ∂x − ∆ tω ∂∂t − ∆ t ω ∂ ∂t − i ∆ x ∆ tω ∂ ∂x∂t − . . . , (16)we can rewrite the LBE into a set of three coupled PDEs for the distribution functions(1 − L i ) f i ( x, t ) = f eqi ( x, t ) , ∀ i ∈ {− , , } , (17)that holds for x ∈ ] − ∞ , ∞ [ and t ∈ [0 , ∞ [.The solution can be found by performing a Picard or fixed point iteration f ( n +1) i = L i f ( n ) i + f eqi , (18)with initial guess f ( − i = 0 that results in f (0) i = f eqi =: g i ,f (1) i = g i + L i f (0) i = g i + L i g i ,f (2) i = g i + L i g i + L i ( L i g i ) ,. . .f ( n ) i = n (cid:88) k =0 L ki g i , f ( n ) i the n -th iterate. This iteration converges if the error between subsequent iterations goesto zero.In contrast to traditional iterations, which require convergence for any initial guess, Eq. (18) isa fixed point iteration with initial guess zero and a smooth right hand side. It is only necessary toshow convergence for this particular case. To discuss this convergence we introduce the 2-norm, (cid:107) . (cid:107) , to show what happens between subsequent iterations. The absolute difference of subsequentiterations is given by (cid:107) f ( n +1) i − f ( n ) i (cid:107) = (cid:107)L n +1 i g i (cid:107) . This goes to zero if f i is smooth enough, implying smoothness on ρ and g i = f eqi ( x, t ) such thatlim n →∞ (cid:107)L n +1 i g i (cid:107) = 0. This smoothness condition depends on the parameters of the LBM, ∆ x ,∆ t and ω , and the derivatives of g i . For example, when g i can be described by a polynomial, wehave that there exists a k such that for all n > k applies that (cid:107)L n +1 i g i (cid:107) = 0.We end up with the series (14) that consists of the vectors of constants given in (15), thedensity and its derivatives. Once the constants are determined, the lifting operator — that isnecessary to initialize the LBM and to determine the ghost points in the hybrid model — can beconstructed. How these constants are found is discussed in Sections 4.2 and 4.3. Section 4.2 dealswith the analytical derivation while Section 4.3 is concerned with the numerical procedure. As asurplus, it allows us to find the corresponding macroscopic PDE as outlined in Section 4.4. With the help of expansion (14) it is possible to build a lifting operator that constructs thedistribution functions for a given density. The focus of this section lies in the determination of thevectors of constants (15), the coefficients of such a lifting operator (14). To simplify the discussionand notation we limit ourselves to a truncated series f ( x, t ) = f eq ( x, t ) + α ∂ρ∂x + β ∂ ρ∂x + γ ∂ρ∂t , (19)where α , β and γ are the vectors containing the constants. The method is easily generalized toinclude higher order terms which will be considered in Section 4.3.2.Using the fact that Eq. (19) is valid for every possible grid point, we can consider three gridpoints x j , x k and x l and set up a linear system for the nine unknowns, namely three vectors eachcontaining three constants. Where j , k and l are certain indices determined in a later stage of thepaper. ∂ρ ( x j ) ∂x ∂ ρ ( x j ) ∂x ∂ρ ( x j ) ∂t∂ρ ( x j ) ∂x ∂ ρ ( x j ) ∂x ∂ρ ( x j ) ∂t∂ρ ( x j ) ∂x ∂ ρ ( x j ) ∂x ∂ρ ( x j ) ∂t∂ρ ( x k ) ∂x ∂ ρ ( x k ) ∂x ∂ρ ( x k ) ∂t∂ρ ( x k ) ∂x ∂ ρ ( x k ) ∂x ∂ρ ( x k ) ∂t∂ρ ( x k ) ∂x ∂ ρ ( x k ) ∂x ∂ρ ( x k ) ∂t∂ρ ( x l ) ∂x ∂ ρ ( x l ) ∂x ∂ρ ( x l ) ∂t∂ρ ( x l ) ∂x ∂ ρ ( x l ) ∂x ∂ρ ( x l ) ∂t∂ρ ( x l ) ∂x ∂ ρ ( x l ) ∂x ∂ρ ( x l ) ∂t α α α − β β β − γ γ γ − = f ( x j , t ) − f eq ( x j , t ) f ( x j , t ) − f eq ( x j , t ) f − ( x j , t ) − f eq − ( x j , t ) f ( x k , t ) − f eq ( x k , t ) f ( x k , t ) − f eq ( x k , t ) f − ( x k , t ) − f eq − ( x k , t ) f ( x l , t ) − f eq ( x l , t ) f ( x l , t ) − f eq ( x l , t ) f − ( x l , t ) − f eq − ( x l , t ) . (20) For a given f i ( x, t ) where i ∈ {− , , } , the linear system (20) will give the coefficients α , β and γ . However, linear system (20) only delivers the correct coefficients if f i is smooth enoughsuch that f eqi satisfies the smoothness condition. This is the case when f i lies on the slow manifold.The Constrained Runs algorithm offers a way to reach the slow manifold in an iterative way.We combine the ideas of the CR-algorithm to reach the slow manifold and the Chapman-Enskog expansion to find the unknown constants on this slow manifold. The numerical procedureto do so is given in Section 4.3. 10f a PDE in closed form exists that describes the evolution of ρ in the form of ρ t + aρ x = Dρ xx ,then the linear system (20) will be singular. Indeed, the PDE will give a relation between ρ t , ρ x and ρ xx in each of the grid points x j , x k and x l . As a result, every element in the last threecolumns of the linear system (20) can be written as a linear combination of the first six columns.In practice, however, the PDE is only an approximation and the system will be close to singular.This clarifies why we do not solve the linear system in (20) and first focus on one that is notclose to singular in Section 4.3. Section 4.4 explains why such a PDE exists. From the previous discussion it is clear that the coefficients of the lifting operator can be extractedfrom a linear system once f approaches the slow manifold. Next, the extraction of the coefficientsis combined with the CR-algorithm to bring f close to the slow manifold. This discussion is limited, as in Section 4.2, to the first few terms of the expansion. The singularsystem can be avoided by taking a series that only contains spatial derivatives. Such a series canrepresent the same state since the time derivative ∂ t ρ is often related to the spatial derivativesthrough a macroscopic PDE. For example, suppose that a PDE of the form ρ t + aρ x = Dρ xx describes the behavior of ρ . It is then possible to eliminate ∂ t ρ from the expansion. The coefficientsare then ˜ α = α − γa and ˜ β = β + γD . The distribution functions are now series with only spatialderivatives. f ( x, t ) = f eq ( x, t ) + ˜ α ∂ρ∂x + ˜ β ∂ ρ∂x . (21) Remark 1
Rewrite ˜ α and ˜ β as α and β but bear in mind that it considers different coefficientsin Eq. (19) and Eq. (21) . Again, once the distribution functions are close to the slow manifold, we can extract thecoefficients α and β from the linear system ∂ρ ( x j ) ∂x ∂ ρ ( x j ) ∂x ∂ρ ( x j ) ∂x ∂ ρ ( x j ) ∂x ∂ρ ( x j ) ∂x ∂ ρ ( x j ) ∂x ∂ρ ( x k ) ∂x ∂ ρ ( x k ) ∂x ∂ρ ( x k ) ∂x ∂ ρ ( x k ) ∂x ∂ρ ( x k ) ∂x ∂ ρ ( x k ) ∂x α α α − β β β − = f ( x j , t ) − f eq ( x j , t ) f ( x j , t ) − f eq ( x j , t ) f − ( x j , t ) − f eq − ( x j , t ) f ( x k , t ) − f eq ( x k , t ) f ( x k , t ) − f eq ( x k , t ) f − ( x k , t ) − f eq − ( x k , t ) . (22)To reach the slow manifold we combine Constrained Runs with the extraction of the coefficients.Consider a numerical function h ( α, β ; ρ, m ) as described in Function 1. This function takes asinput α and β and as parameters a fixed density ρ and an integer m , the order of the smoothnesscondition. It first constructs, with this input, a state f with the help of series (21). This state isthen used to perform multiple LBM steps. For each of these steps we can find the correspondingmoments φ and ξ . On the moments, we can use the CR-algorithm (Section 3.2) to find new11oments that are closer to the slow manifold by considering the finite difference approximationsof the m -th order smoothness condition, d m +1 φdt m +1 = 0 and d m +1 ξdt m +1 = 0 . (23)These new moments result in new coefficients α and β , by applying the linear system Eq. (22) onthe distribution functions f i , corresponding to the new φ and ξ and the given ρ .The idea is now to determine α and β such that they are invariant under this numerical function h ( α, β ; ρ, m ). Indeed, if the initial and final state can be described by the same α and β then thelifted f is close to the slow manifold since it is a fixed point of the underlying CR-iteration.Instead of performing a regular fixed point iteration with h ( α, β ; ρ, m ), a Newton iterationis used that finds a := ( α, β ) ∈ R such that a = h ( a ; ρ, m ). This reduces the computationalcost significantly because the size of the Jacobian system with α and β is much smaller then theJacobian of the original Constrained Runs algorithm. The latter involves the moments in everygrid point and this becomes very large. Function 1 h ( α, β ; ρ, m ) Require:
Guess on coefficients α , β , given density ρ , order m to use in Eq. (23). Construct lifting operator f = f eq + α ∂ρ∂x + β ∂ ρ∂x (Eq. (21)). Compute corresponding moments φ and ξ by applying Eq. (3). Perform m + 1 LBM time steps to compute d m +1 φdt m +1 = 0 and d m +1 ξdt m +1 = 0 by using forward finitedifference formulas. This results in new moments φ and ξ . { find moments closer to slowmanifold } Revert back to distribution functions f by applying Eq. (3). Select grid points x j , x k and x l to construct linear system (22). Solve the system for α and β . return α , β .The numerical function h ( α, β ; ρ, m ) has a density ρ ( x,
0) as a parameter and the solution forthe coefficients is independent of its choice of ρ . The coefficients are determined by functional L i (16) that only depends on constants ω , the spatial grid size ∆ x and time step ∆ t . Since thecoefficients do not depend on time, we can choose an arbitrary density ρ ( x,
0) such that Eq. (22)is easily solvable.The choice of the grid points x j and x k in (22) should be such that the condition number ofthe matrix is optimal. In addition, the spatial derivatives that are considered needs to exist andshould not become zero during the LBM evolution since otherwise we would end up with singularlinear systems.Furthermore, the test domain used in the LBM inside the function h ( α, β ; ρ, m ) can be sig-nificantly smaller than the domain of the original LBM problem. A smaller test domain will notaffect the constant coefficients of the lifting operator. However, it should use the same ∆ x and ∆ t as the LBM of interest since the coefficients depend on the chosen spacings in space and time. Thechoice for the test domain, density and indices is further discussed in Section 5.1 for the consideredmodel problem. There are two ways to increase the accuracy. First, more terms in the expansion can be consideredsuch that more derivatives of the density are taken into account. Second, we can enforce a higherorder smoothness in the CR-algorithm. Both methods are outlined below.The proposed method can easily be extended by considering more terms with higher orderderivatives in the truncated series (21). For example, consider the expansion f = f eq + α ∂ρ∂x + β ∂ ρ∂x + δ ∂ ρ∂x + (cid:15) ∂ ρ∂x (24)12hat now requires the determination of more coefficients that are found by considering — inaddition to x j and x k — additional grid points x l and x m . This leads to a larger system ofunknowns but will give better results.Higher order smoothness can be enforced on the moments φ and ξ as in the CR-algorithm byconsidering a higher order m in Eq. (23). This requires more LBM steps and uses a higher orderfinite difference formula to estimate the derivatives in time.For further conclusions and results higher order derivatives and higher order smoothness aretaken into account. Next, we derive from Eq. (14), the macroscopic PDE by summing over the velocities. Using (cid:80) i f i ( x, t ) = ρ ( x, t ) = (cid:80) i f eqi ( x, t ) results in a macroscopic PDE for the density. (cid:32)(cid:88) i γ i (cid:33) ∂ρ∂t + (cid:32)(cid:88) i ζ i (cid:33) ∂ ρ∂t + . . . + (cid:32)(cid:88) i η i (cid:33) ∂ ρ∂x∂t + . . . = − (cid:32)(cid:88) i α i (cid:33) ∂ρ∂x − (cid:32)(cid:88) i β i (cid:33) ∂ ρ∂x − (cid:32)(cid:88) i δ i (cid:33) ∂ ρ∂x − (cid:32)(cid:88) i (cid:15) i (cid:33) ∂ ρ∂x − . . . . Series (14) derived in this setting leads to the classical Chapman-Enskog expansion [28, 5] thatwe obtained in Section 2.2. Indeed, Eq. (14) is written as f i ( x, t ) = g i + L i g i + L i g i + . . . , because of the application of the fixed point iteration. When the series is truncated after thesecond order spatial derivative and the first order time derivative, we end up with f i ( x, t ) = 13 ρ − i ∆ x ω ∂ρ∂x + (cid:18) i ∆ x ω − i ∆ x ω (cid:19) ∂ ρ∂x − ∆ t ω ∂ρ∂t . Summing f i ( x, t ) over i ∈ {− , , } we obtain the same macroscopic diffusion PDE (4). Substi-tuting this PDE in the series to remove the time derivative leads to the classical Chapman-Enskogexpansion in Eq. (5).This macroscopic PDE was used for the removal of the time derivative in Eq. (19) and itsreplacement with Eq. (21). Furthermore, it is important to note that the macroscopic PDE is notnecessarily of the reaction-diffusion prototype. Truncating (14) after more terms and taking morederivatives into account results in a better output but it will lead to the term ∂ ρ∂t which gives aless comfortable macroscopic PDE.This relation to the macroscopic PDE can now be integrated in the numerical Chapman-Enskogmethod. Once the fixed point described in Section 4.3.1 is found, we have α and β that lifts ρ to the distribution functions close to the slow manifold. By performing two more LBM steps, ∂ρ∂t can be calculated by using a forward finite difference formula. System (20) can be applied to findthe vectors of constants of this larger system that include the time derivative. There are now twopossibilities: either the resulting system is non-singular and it can be solved for the coefficientsand only an approximate PDE can be found as is considered above. Or it is too singular to besolved accurately but then the PDE can be extracted from the nullspace of the system.Let us first discuss the situation where the matrix in (20) is non-singular. The system canthen be solved for α , β and γ . The approximate PDE can be determined by summing over theobtained coefficients. ∂ρ∂t = − (cid:80) i α i (cid:80) i γ i ∂ρ∂x − (cid:80) i β i (cid:80) i γ i ∂ ρ∂x . This PDE is only approximate. Otherwise, if it would hold exactly, the system would be singularas expected. 13or a singular system, we know that one or more of the eigenvalues will be zero with a cor-responding null eigenvector. Focusing on the null eigenvector v = { v , v , . . . , v } , we know that Av = 0 with A the matrix in system (20). Using this, we obtain ∂ρ ( x j ) ∂x v + ∂ ρ ( x j ) ∂x v + ∂ρ ( x j ) ∂t v = 0 , from which we conclude that the resulting PDE looks like ∂ρ ( x j ) ∂t = − v v ∂ρ ( x j ) ∂x − v v ∂ ρ ( x j ) ∂x . Remark that same PDE will be found when considering the equation in grid points x k and x l instead of x j . The results of the previous sections are now combined in an algorithm that delivers a liftingoperator and an approximate macroscopic PDE. This can be used, for example, to constructa hybrid model. The pseudocode is presented in Algorithm 2 while the complete algorithm ispresented below. The algorithm starts by searching for the lifting operator on the basis of thespatial derivatives. Thereafter, it inserts time derivatives and calculates the coefficients of themacroscopic PDE.Start with an initial guess for { α, β, δ, (cid:15) } in Eq. (24). Apply Function 1 h ( α, β, δ, (cid:15) ; ρ, m ) for agiven ρ and a certain m for the order of smoothness. This results in coefficients { α, β, δ, (cid:15) } thatrepresent distribution functions closer to the slow manifold. The lifting operator is constructed atthis point.When these distribution functions are found based on the spatial derivatives only, we still needto determine the corresponding PDE by considering the null eigenvector or by a summation of thecoefficients as discussed in Section 4.4. By performing two extra LBM steps — to estimate thetime derivative with a forward finite difference formula — the coefficient γ belonging to the timederivative of the expansion below can be numerically calculated. f = f eq + α ∂ρ∂x + β ∂ ρ∂x + δ ∂ ρ∂x + (cid:15) ∂ ρ∂x + γ ∂ρ∂t . Since the PDE can be obtained from the numerically constructed distribution functions, the PDEobtained through the Chapman-Enskog expansion does not need to be obtained analytically.
Remark 2
Note that one can also consider f eq ( x, t ) = κρ ( x, t ) and determine the constants ofvector κ in a similar setting by using an extra grid point to obtain a larger system of unknowns. Remark 3
In this paper we have chosen to use the same ∆ x and ∆ t in the LBM as in the PDE.However, their stability properties may be different. Our specific LBM simulation is stable in the2-norm when ≤ ω ≤ [23]. However, it is not necessary that the macroscopic PDE, when itis discretized with the same ∆ x and ∆ t and forward Euler, is also stable. Indeed, when ω → for a fixed ∆ x and ∆ t the resulting diffusion coefficient D grows, see (4) , and this can lead to aninstability. The pseudocode of the numerical Chapman-Enskog expansion as a lifting operator is given inAlgorithm 2 together with the determination of the transport coefficients of the PDE to constructa hybrid model. The proposed algorithm has several advantages. In contrast to the Chapman-Enskog expansion no analytical work is required. Compared to the Constrained Runs algorithm itsignificantly reduces the number of unknowns in the lifting since it only needs to find the coefficients(vectors of constants) rather than the full state of the distribution functions. Furthermore, it canbe done off-line before the calculations. Indeed, once the coefficients are found they can be reusedevery time step to realize the lifting. As an extra surplus, the PDE can be determined to constructhybrid models. 14 lgorithm 2
Pseudocode numerical Chapman-Enskog expansion
Require:
Test domain that defines ρ ( x, a = { α, β, δ, (cid:15) } = zeros(12,1) and auser-defined tolerance tol, parameter m for higher order smoothness. repeat a k +1 = a k − (cid:0) J ( a k ) (cid:1) − h ( a k ; ρ, m ) with h defined in Function 1. until (cid:107) a k +1 − a k (cid:107) < tol.Result for coefficients { α, β, δ, (cid:15) } belonging to the spatial derivatives.Distribution functions f = f eq + α ∂ρ∂x + β ∂ ρ∂x + δ ∂ ρ∂x + (cid:15) ∂ ρ∂x closer to the slow manifold.Perform 2 more LBM steps to determine ∂ρ∂t numerically by a forward finite difference formula.Construct system (20) — including higher order spatial derivatives — to achieve coefficients { α, β, δ, (cid:15), γ } .Determine coefficients of PDE by using the nullspace of the system or by summation of thecoefficients (Section 4.4). return Lifting operator f = f eq + α ∂ρ∂x + β ∂ ρ∂x + δ ∂ ρ∂x + (cid:15) ∂ ρ∂x and macroscopic PDE. The new lifting operator is now illustrated in several examples. First we benchmark its accuracyagainst a reference solution that is reconstructed. This is done in Section 5.1. In Section 5.2 werecall the one-dimensional hybrid model of Figure 1. Two-dimensional problems are discussedin Section 5.3. The important comparison of the additional required LBM steps to perform thelifting in a hybrid model is presented in Section 5.4.
The proposed lifting operator can be tested against a reference distribution function f c . Thisreference solution is calculated by performing 1000 lattice Boltzmann steps starting from an initialstate that corresponds to the equilibrium distribution function of a given density ρ .The lifting operator can now be evaluated by restricting the reference distribution function f c to its density ρ = (cid:80) i f i ( x, t ) and lift it back to a distribution function f by using the proposedlifting operator. The resulted f will be compared with f c with the help of the 2-norm (cid:107) f − f c (cid:107) . Example 1
The considered model problem has the following parameters for a one-dimensionaldomain of length L . L = 10 , n = 200 , ∆ x = Ln = 0 . , ∆ t = 0 . , ρ ( x,
0) = e − ( x − L ) , ω = 0 . . For these parameters the classical Chapman-Enskog expansion predicts a diffusion coefficient D = 1 (Eq. (4) ). To reproduce the numerical results linked to Example 1 we include some extra information onhow to determine the indices of (22). As mentioned in Section 4.3.1 we can choose an arbitraryinitial density ρ ( x,
0) and a test domain to determine the constants of the lifting operator (21).For example, consider ρ ( x,
0) = x + 1 / x for unknowns α, β and ρ ( x,
0) = x + 1 / x + 1 / x forunknowns α, β, δ . x is defined by the spatial nodes of the test domain defined below. Furthermore,the test domain reduces the computational expense compared to the actual spatial domain.Consider the domain parameters of Example 1. Since we know that the constants are onlyaffected by these space and time steps, we should consider — together with ρ ( x,
0) — a testdomain with the same step sizes since the vectors of constants (15) are affected by these choices.The test domain is of length L test = 3 such that n test = 60 since ∆ x = 0 .
05. This number of gridpoints will make it possible to choose the indices x j , x k , . . . such that system (22) is not close tosingular. We can return now to the question which indices should be used in system (22). Focus15n the fact that we do not want an effect of wrongly chosen boundary conditions in the smallertest domain. The grid points should be taken far enough from the edges and in points such thatthe system (22) does not become singular. The indices can be, for example, 10, 20, 30, . . . spreadover the test domain of 60 grid points.Note that one can also focus on local updates around the considered grid points x j , x k , . . . .When the number of iterations needed in Newton’s method are known, one knows how many LBMsteps will be performed to find the new coefficients α and β . Then the size of the test domain canbe shrunk to a smaller domain around x j , x k , . . . . This is the same idea as used in [33] to performlocal updates for the CR-algorithm.To compare the proposed lifting operator with the existing ones discussed in Section 3, theresults for (cid:107) f − f c (cid:107) of the different lifting operators are included in Tables 1, 2 and 3.Table 1 contains results obtained with the analytical Chapman-Enskog expansion as a liftingoperator. The first column gives the order of the expansion and the second column shows (cid:107) f − f c (cid:107) .As expected a better accuracy is obtained with higher order expansions. For example, by takingthe third derivative of the density into account an error of 1.24e-5 is achieved.Table 1: The error (cid:107) f − f c (cid:107) is presented to test the exact Chapman-Enskog expansion as a liftingoperator. The reference solution is obtained after 1000 LBM steps in the model problem fromExample 1. Construction lifting operator: (cid:107) f − f c (cid:107) exact Chapman-Enskog expansion f = f eq f = f eq + α exact ∂ρ∂x f = f eq + α exact ∂ρ∂x + β exact ∂ρ ∂x f = f eq + α exact ∂ρ∂x + β exact ∂ρ ∂x + δ exact ∂ ρ∂x In Table 2 we use the Constrained Runs algorithm of various orders of accuracy to numericallylift the density to distribution functions. Different types of backward extrapolation are listed inthe first column of the table while the corresponding 2-norm (cid:107) f − f c (cid:107) is described in the secondcolumn. There we see that very accurate results can be found for the higher order versions. Notethat these methods find for each grid point the moments φ and ξ of the distribution functions.Together with ρ , the corresponding distribution functions are found by Eq. (3). Since this givesa local solution it can give accurate results. The last column contains (cid:107) f − f c (cid:107) when some extraadvection effect is included, which shows similar results as the pure diffusion problem.Table 2: The error (cid:107) f − f c (cid:107) with the Constrained Runs algorithm (combined with Newton’smethod) for various orders of accuracy as a lifting operator. The reference solution is obtainedfor the model problem in Example 1 by performing 1000 LBM time steps before restricting andlifting. The last column contains results when an extra advection effect of a = 0 .
66 is included —which changes the used equilibrium distribution functions as noted in Eq. (6).
Construction lifting operator: (cid:107) f − f c (cid:107) (cid:107) f − f c (cid:107) extrapolation CR-algorithm pure diffusion D = 1 plus advection a = 0 . Table 3 shows the results with the proposed numerical Chapman-Enskog expansion as a liftingoperator. We clearly see that taking more terms in the expansion, i.e. more derivatives of ρ ( x, t )in the lifting, leads to a better lifting operator. In the same table we show the results with higherorder smoothness conditions by using Eq. (23) with higher order m . As in the CR-algorithm,16igher order smoothness does result in a significant improvement. The accuracy increases to9.45e-11 when up to the sixth spatial derivative is taken into account. Including advection in thistable will also show similar results but these are not added.Table 3: The error (cid:107) f − f c (cid:107) with the numerical Chapman-Enskog expansion as a lifting opera-tor. The reference distribution function f c is obtained by performing 1000 LBM time steps withparameters listed in Example 1. Each of the lifting operators is calculated as a fixed point for thecoefficients. The first table shows the results with one LBM step before updating the moments,implying an update of the coefficients. There are results listed for increasing number of terms inthe expansion, implying an increasing number of considered coefficients. The second table showsthe same results where two LBM steps are used to estimate the smoothness. The third and fourthtable show the results with a quadratic, respectively cubic computation of the finite differenceapproximation in Eq. (23). The constant coefficients found via the numerical procedure are com-pared to those found exactly with the classical Chapman-Enskog expansion in columns 3, 4 and5. Construction lifting operator: (cid:107) f − f c (cid:107) (cid:107) α − α exact (cid:107) (cid:107) β − β exact (cid:107) (cid:107) δ − δ exact (cid:107) Numerical Chapman-Enskogconstant computation of fixed point f = f eq + α ∂ρ∂x f = f eq + α ∂ρ∂x + β ∂ρ ∂x f = f eq + α ∂ρ∂x + β ∂ρ ∂x + δ ∂ ρ∂x f = f eq + α ∂ρ∂x + β ∂ρ ∂x + δ ∂ ρ∂x + (cid:15) ∂ ρ∂x f = f eq + α ∂ρ∂x + β ∂ρ ∂x + δ ∂ ρ∂x + (cid:15) ∂ ρ∂x + θ ∂ ρ∂x f = f eq + α ∂ρ∂x + β ∂ρ ∂x + δ ∂ ρ∂x + (cid:15) ∂ ρ∂x + θ ∂ ρ∂x + ι ∂ ρ∂x (cid:107) f − f c (cid:107) (cid:107) α − α exact (cid:107) (cid:107) β − β exact (cid:107) (cid:107) δ − δ exact (cid:107) f = f eq + α ∂ρ∂x f = f eq + α ∂ρ∂x + β ∂ρ ∂x f = f eq + α ∂ρ∂x + β ∂ρ ∂x + δ ∂ ρ∂x f = f eq + α ∂ρ∂x + β ∂ρ ∂x + δ ∂ ρ∂x + (cid:15) ∂ ρ∂x f = f eq + α ∂ρ∂x + β ∂ρ ∂x + δ ∂ ρ∂x + (cid:15) ∂ ρ∂x + θ ∂ ρ∂x f = f eq + α ∂ρ∂x + β ∂ρ ∂x + δ ∂ ρ∂x + (cid:15) ∂ ρ∂x + θ ∂ ρ∂x + ι ∂ ρ∂x (cid:107) f − f c (cid:107) (cid:107) α − α exact (cid:107) (cid:107) β − β exact (cid:107) (cid:107) δ − δ exact (cid:107) f = f eq + α ∂ρ∂x f = f eq + α ∂ρ∂x + β ∂ρ ∂x f = f eq + α ∂ρ∂x + β ∂ρ ∂x + δ ∂ ρ∂x f = f eq + α ∂ρ∂x + β ∂ρ ∂x + δ ∂ ρ∂x + (cid:15) ∂ ρ∂x f = f eq + α ∂ρ∂x + β ∂ρ ∂x + δ ∂ ρ∂x + (cid:15) ∂ ρ∂x + θ ∂ ρ∂x f = f eq + α ∂ρ∂x + β ∂ρ ∂x + δ ∂ ρ∂x + (cid:15) ∂ ρ∂x + θ ∂ ρ∂x + ι ∂ ρ∂x (cid:107) f − f c (cid:107) (cid:107) α − α exact (cid:107) (cid:107) β − β exact (cid:107) (cid:107) δ − δ exact (cid:107) f = f eq + α ∂ρ∂x f = f eq + α ∂ρ∂x + β ∂ρ ∂x f = f eq + α ∂ρ∂x + β ∂ρ ∂x + δ ∂ ρ∂x f = f eq + α ∂ρ∂x + β ∂ρ ∂x + δ ∂ ρ∂x + (cid:15) ∂ ρ∂x f = f eq + α ∂ρ∂x + β ∂ρ ∂x + δ ∂ ρ∂x + (cid:15) ∂ ρ∂x + θ ∂ ρ∂x f = f eq + α ∂ρ∂x + β ∂ρ ∂x + δ ∂ ρ∂x + (cid:15) ∂ ρ∂x + θ ∂ ρ∂x + ι ∂ ρ∂x A comparison of Tables 1, 2 and 3 shows that the proposed numerical lifting operator leadsto better results than the analytically found Chapman-Enskog expansion. As can be seen is theConstrained Runs algorithm a good lifting method, but, as will be discussed in Section 5.4, thecomputational expense of this method brings down the beauty of it. Table 4 of Section 5.4 contains17 comparison of the number of additional LBM steps required for each of the lifting operators. Inthe CR-algorithm this additional cost can be attributed to the construction of the Jacobian matrix.These additional LBM steps make the method computationally very expensive. In two dimensionsthis method becomes prohibitive. This makes the numerical Chapman-Enskog lifting operator agood alternative to the CR-algorithm that gives a similar accuracy at a limited computationalcost.
To compare the results of the numerical Chapman-Enskog expansion with the earlier proposedlifting operators discussed in Section 3, Figures 3 and 4 show the results of the absolute difference | ρ hybrid − ρ LBM | for the exact Chapman-Enskog expansion and those obtained with the ConstrainedRuns algorithm. ρ hybrid is the density of the hybrid model and ρ LBM the density of a full LBM. ρ LBM is the reference solution to compare the hybrid solution with. It considers a LBM on thewhole spatial domain [ a, b ] with the parameters outlined in Example 1 and the domain representedin Figure 1. The lifting operators are used both to initialize the LBM and to find the ghost pointsof the LBM domain. Figure 3 shows the absolute differences by using as lifting operator theexact Chapman-Enskog expansion respectively up to zeroth (top left), first (top right), second(bottom left) and third order (bottom right). Figure 4 shows the absolute differences with thelifting operator based on the Constrained Runs algorithm in combination with Newton’s methodfor respectively a constant (top left), linear (top right), quadratic (bottom left) and cubic (bottomright) extrapolation in time. These results were obtained in [33] by considering local updates atthe ghost points of the LBM domain.When the numerical Chapman-Enskog expansion (up to the sixth spatial derivative) is used inour one-dimensional hybrid model problem, Figure 5 is obtained. Here, we have two possibilities.First, act as if we know the PDE (4) obtained from the exact Chapman-Enskog expansion. | ρ hybrid − ρ LBM | is given in the left Figure 5 for which the hybrid domain is shown in Figure 1 and the PDE isthe one given in Eq. (4). Second, use the PDE that is obtained from the proposed lifting operatorthrough summing the proposed lifting operator or considering the nullspace as explained in Section4.4. With this PDE, the result for | ρ hybrid − ρ LBM | is shown in the right Figure 5.As can be seen in Figure 5, a change in the PDE — by considering the PDE obtained throughthe numerical Chapman-Enskog expansion — results in an even smaller modeling error comparedto the one obtained via the classical Chapman-Enskog expansion.Changing the parameters of the model such that advection is included, is considered below.The figures show similar results with advection-term a = 0 .
66. The model problem remains theone from Example 1. The only difference is the change in the equilibrium distribution as shownin (6). Figure 6 contains the comparison results obtained through the CR-algorithm. Figure 7(left) shows | ρ hybrid − ρ LBM | with the numerical Chapman-Enskog expansion as a lifting operatorand the PDE obtained through the Chapman-Enskog expansion while Figure 7 (right) uses thePDE obtained from the proposed lifting operator. This section generalizes the previous one. Two spatial dimensions are considered. Two-dimensionalproblems can take different discrete sets of velocities into account. In Section 5.3.1 results for D2Q5are presented while Section 5.3.2 contains results for D2Q9.
The hybrid test domain for this section is represented in Figure 8 for D2Q5. Again, the domainis split into subdomains. One part of the domain is described by the LBM while another part isdescribed by a macroscopic PDE. Example 2 describes the parameters for the model problem inthis two-dimensional setting. 18
DE LBM00 . . | ρ h y b r i d − ρ L B M | j : x j , j ∈ { , . . . , } PDE LBM05 · − | ρ h y b r i d − ρ L B M | j : x j , j ∈ { , . . . , } PDE LBM01 · − | ρ h y b r i d − ρ L B M | j : x j , j ∈ { , . . . , } PDE LBM01 · − | ρ h y b r i d − ρ L B M | j : x j , j ∈ { , . . . , } Figure 3: | ρ hybrid − ρ LBM | after 200 time steps. The difference is also shown at earlier time slots,but shifted down for clarity. The lines represent time steps between one and 200. The topline corresponds to time step 191 while the bottom line represents time step 11. The lines inbetween correspond to jumps with 10 time steps from 11 to 21, . . . ,181 to 191. The domain thatis considered, is shown in Figure 1. The lifting operator is f i ( x, t ) = 1 / ρ ( x, t ) (top left), firstorder (top right), second order (bottom left) and third order Chapman-Enskog (bottom right)respectively. The lifting operator is used both to find the ghost points of the LBM domain andfor the creation of the initial state for the LBM region. The model problem parameters are listedin Example 1. 19 DE LBM00 . | ρ h y b r i d − ρ L B M | j : x j , j ∈ { , . . . , } PDE LBM05 · − | ρ h y b r i d − ρ L B M | j : x j , j ∈ { , . . . , } PDE LBM05 · − | ρ h y b r i d − ρ L B M | j : x j , j ∈ { , . . . , } PDE LBM05 · − | ρ h y b r i d − ρ L B M | j : x j , j ∈ { , . . . , } Figure 4: | ρ hybrid − ρ LBM | after 200 time steps where the lifting operator based on the CR-algorithmis used in combination with the method of Newton. We show results for constant (top left), linear(top right), quadratic (bottom left) and cubic (bottom right) backward extrapolation, respectively.The difference is also shown at earlier time slots, but shifted down for clarity. The domain that isconsidered, is shown in Figure 1. The model problem parameters are listed in Example 1. PDE LBM05 · − | ρ h y b r i d − ρ L B M | j : x j , j ∈ { , . . . , } PDE LBM05 · − | ρ h y b r i d − ρ L B M | j : x j , j ∈ { , . . . , } Figure 5: | ρ hybrid − ρ LBM | after 200 time steps in the model problem of Example 1. The domainthat is considered, is shown in Figure 1. To deal with the initial error and the error in the ghostpoints of the LBM domain the numerical Chapman-Enskog expansion (order spatial expansion 6)is used. The considered PDE in the hybrid domain is the analytically known PDE (4) in the leftfigure and the one that is obtained from the numerical Chapman-Enskog expansion in the rightfigure. 20 DE LBM00 . | ρ h y b r i d − ρ L B M | j : x j , j ∈ { , . . . , } PDE LBM01 · − · − | ρ h y b r i d − ρ L B M | j : x j , j ∈ { , . . . , } PDE LBM01 · − · − | ρ h y b r i d − ρ L B M | j : x j , j ∈ { , . . . , } PDE LBM01 · − · − | ρ h y b r i d − ρ L B M | j : x j , j ∈ { , . . . , } Figure 6: | ρ hybrid − ρ LBM | after 200 time steps where the lifting operator based on the CR-algorithmis used in combination with the method of Newton. We show results for constant (top left), linear(top right), quadratic (bottom left) and cubic (bottom right) backward extrapolation, respectively.The difference is also shown at earlier time slots, but shifted down for clarity. The domain thatis considered, is shown in Figure 1. Model parameters are listed in Example 1 with an extraadvection coefficient a = 0 . PDE LBM01 · − · − | ρ h y b r i d − ρ L B M | j : x j , j ∈ { , . . . , } PDE LBM01 · − · − | ρ h y b r i d − ρ L B M | j : x j , j ∈ { , . . . , } Figure 7: | ρ hybrid − ρ LBM | after 200 time steps in the model problem of Example 1 with advectioneffect a = 0 .
66. The domain that is considered, is shown in Figure 1. To deal with the initial errorand the error in the ghost points of the LBM domain the numerical Chapman-Enskog expansion(order spatial expansion 4) is used. The considered PDE in the hybrid domain is the analyticallyknown PDE (7) in the left figure and the one that is obtained from the numerical Chapman-Enskogexpansion in the right figure. 21 xample 2
The considered model problem has the following parameters for a two-dimensionaldomain — described by 5 possible velocity directions (D2Q5) — of length L × L (with n thenumber of grid points). L = 10 , n = 200 , ∆ x = ∆ y = Ln , ∆ t = 0 . , ω = 1 . . For these parameters the classical Chapman-Enskog expansion predicts a diffusion coefficient D = 1 (Eq. (8) ). The comparison of | ρ hybrid − ρ LBM | is represented in Figure 9 for Example 2. The different liftingoperators are used to obtain distribution functions from a given density. The used lifting operatorsare the equilibrium distribution function in the top left figure, the first order Chapman-Enskogexpansion in the top right, the second order Chapman-Enskog expansion in the middle left, thenumerical Chapman-Enskog expansion of spatial order expansion 4 in the middle right and thebottom — depending on the used PDE in the hybrid model. albx x p x p +1 x n − LBM domainPDE domain yxf i ( x j ,y k ,t ) i ∈ { , ,..., } for D2Q5 j ∈ { ,...,p } k ∈ { ,...,n − } ρ ( x j ,y k ,t ) j ∈ { p + 1 ,...,n − } k ∈ { ,...,n − } Figure 8: The two-dimensional spatial domain [ a, b ] × [ a, b ] ⊂ R in the hybrid model is split into[ a, l [ × [ a, b ] on which we solve the LBM and [ l, b ] × [ a, b ] on which we solve the PDE model. Thesolid points ( • ) represent the grid for the density ρ of the discrete PDE, the circles ( ◦ ) representthe LBM variables f i ( x, y, t ) , i ∈ { , . . . , } for D2Q5. The periodic boundaries and the couplingare implemented with ghostcells which are drawn by dashed circles. The density in the ghostcellsof the PDE domain, in ( x p , y k ), k ∈ { , . . . , n − } and ( x n , y k ), are found by taking (cid:80) i f i in( x p , y j ) and ( x , y k ), respectively. However, the ghostcells for the LBM domain, in ( x − , y k ) and( x p +1 , y k ), require a lifting operator that lifts ρ to the distribution functions in these points. This section takes more directions for the velocities into account. Example 3 contains the modelproblem parameters for D2Q9.
Example 3
The considered model problem has the following parameters for a two-dimensionaldomain — described by 9 possible velocity directions (D2Q9) — of length L × L (with n thenumber of grid points). L = 10 , n = 200 , ∆ x = ∆ y = Ln , ∆ t = 0 . , ω = 1 . . | ρ hybrid − ρ LBM | is plotted at time step 200. To deal with the initialerror and the error in the ghost points of the LBM we use: the equilibrium distribution function(top left), the first order Chapman-Enskog expansion (top right), the second order Chapman-Enskog expansion (middle left), the numerical Chapman-Enskog expansion (order expansion 4)where the PDE in the hybrid domain is the analytically known PDE given in (8) (middle right)and the numerical Chapman-Enskog expansion where the considered PDE in the hybrid domainis the one that is obtained from the numerical Chapman-Enskog expansion (bottom).23 or these parameters the classical Chapman-Enskog expansion predicts a diffusion coefficient D = 1 (Eq. (8) ). First consider no advection in the equilibrium distribution functions ( a = (0; 0)). The resultsfor | ρ hybrid − ρ LBM | are presented in Figure 10. The used lifting operators are the equilibriumdistribution (top left), the first order Chapman-Enskog expansion (top right), the second orderChapman-Enskog expansion (middle left), the numerical Chapman-Enskog expansion (order ex-pansion 4) where the PDE in the hybrid domain is the analytically known PDE given in (9)(middle right) and the numerical Chapman-Enskog expansion where the considered PDE in thehybrid domain is the one that is obtained from the numerical Chapman-Enskog expansion (bot-tom).When advection ( a = (1; 0 . | ρ hybrid − ρ LBM | when the numerical Chapman-Enskog expansion is used to lift density to distributionfunctions. This section compares the computational cost of the lifting operators in the one-dimensional testproblem (Section 5.2).The motivation for this paper is to bring down this cost. Especially the Constrained Runsalgorithm requires many additional LBM steps to lift the density in the ghost points of the LBMdomain. While numerical Chapman-Enskog only requires a single calculation with a fixed costthat can be done off-line before the simulation. This significantly reduces the cost of the lifting.A detailed analysis of the lifting cost in terms of additional LBM steps is listed in Table 4. Thetable is an extension of the results of [33] with results for the classical Chapman-Enskog expansion,Constrained Runs algorithm combined with Newton’s method and the numerical Chapman-Enskogexpansion.It can be seen that the total number of LBM steps for the CR-algorithm are listed per ghostpoint and per time step. The number for the numerical Chapman-Enskog expansion is the totalfor the entire domain and at all time steps since the calculations for the coefficients are doneoff-line.In two-dimensional problems the numerical Chapman-Enskog expansion still has a limitedcomputational cost. Only a few additional coefficients need to be determined associated with theextra spatial derivatives.Note that the computational cost of applying the numerical Chapman-Enskog lifting operatoris the same as applying the analytical Chapman-Enskog operator. For each grid point we need thederivatives of ρ , which can be calculated by finite differences using the densities at neighboringgrid points. This article proposes a numerical lifting operator for lattice Boltzmann models (LBMs) that mapsa given density to the corresponding distribution functions. This new lifting operator is based onthe Chapman-Enskog expansion that writes the missing distribution functions as analytical seriesof the density and its derivatives. The coefficients of this expansion are now determined througha numerical method, in contrast to the original expansion where they are found analytically. Thenumerical method is based on the Constrained Runs algorithm that relies on the attraction of thedynamics toward the slow manifold.A systematic numerical comparison of the accuracy and the computational cost between theanalytical Chapman-Enskog expansion, the Constrained Runs algorithm and the new lifting oper-ator is performed in this article. The cheapest way to lift is with the Chapman-Enskog expansion.However, the analytical expressions are not always available for the system of interest. An al-ternative numerical lifting operator is Constrained Runs (CR), but its computational cost grows24igure 10: | ρ hybrid − ρ LBM | after 200 time steps with model problem presented in Example 3without advection terms. To deal with the initial error and the error in the ghost points of theLBM the equilibrium distribution is used (top left), the first order Chapman-Enskog expansion(top right), the second order Chapman-Enskog expansion (middle left), the numerical Chapman-Enskog expansion (order expansion 4) where the PDE in the hybrid domain is the analyticallyknown PDE given in (8) (middle right) and the numerical Chapman-Enskog expansion where theconsidered PDE in the hybrid domain is the one that is obtained from the numerical Chapman-Enskog expansion (bottom). 25igure 11: | ρ hybrid − ρ LBM | after 200 time steps with model problem presented in Example 3 withadvection a = (1; 0 . ρ to the distribution function f . The table shows these additionalsteps to construct the Jacobian operator for the Newton iteration. The listed values for theCR-algorithm even are per ghost point and per time step. While those of the proposed liftingoperator are the total number for the entire domain and for all time steps together since thevectors of constants needs to be determined only once and can be reused throughout the rest ofthe simulation. Lifting operator Number of LBM steps Total numberiterations to perform of LBM stepsone iteration
Exact Chapman-Enskog / / 0
CR-algorithm per ghost pointtype of extrapolation in time per time step
Constant 3 19 × × × × Numerical Chapman-Enskog for entire domainwith 18 unknowns and all time steps { α, β, δ, (cid:15), θ, ι } Constant 3/(4) 19 57+2=59Linear 3/(4) 19 × × × f ( x, v, t ). And secondly, it can be done off-linebefore the calculations. Indeed, once the coefficients are found they can be reused every time stepand every grid point to realize the lifting, at no significant additional cost. A third advantage isthat the expansion gives, as a spin-off, the transport coefficients of the macroscopic PDE.The new lifting operator, the numerical Chapman-Enskog expansion, is then used in a hybriddomain that spatially couples a macroscopic partial differential equation (PDE) with a latticeBoltzmann model. This creates a missing data problem at the interfaces since the PDE model hastoo few variables to provide the LBM with the correct boundary conditions.The numerical Chapman-Enskog expansion deals with this mismatch in variables. It mapsthe variables of the PDE model to those of the LBM. We evaluate and compare various liftingoperators. In particular, we have focused on a simple LBM and PDE model discretized with equalgrid and time steps such that the error created by the coupling can be highlighted. The paperpresents numerical results both for 1D and 2D hybrid domains where part of the LBM domainis replaced by the macroscopic PDE. In both cases the error associated with the coupling can bemade smaller than the modeling error, related to the PDE approximation of the LBM.This paper reports on our initial efforts where we have focused on a simple model problem withseveral limiting assumptions. In the model we have assumed an equilibrium distribution functionthat depends only on the local density, while in general it also depends on the local momentum andtemperature. This limitation can be easily alleviated by considering a Chapman-Enskog expansionwith a more general equilibrium function.A further assumption is that we used the same time and space grid for the PDE and theLBM. This choice was made to highlight the error made by the coupling mechanism, the ease ofimplementation and to eliminate the error due to the different discretizations. However, there isno reason to prohibit different grid and time spacings. Extra care is then needed to interpolatebetween time and grid spacings. In practice, the grid of the PDE can be further coarsened,depending on local discretization errors. Ideally, the hybrid model is embedded in an adaptivemesh refinement simulation, where at the finest level a LBM is used.We have also kept the boundary between the LBM and the PDE domain fixed during thesimulation at an arbitrary position. In the future, this boundary should be moved adaptivelyusing an accuracy requirement based on the properties of the lifting operator.For the model problem with periodic boundary conditions studied in this paper, the Chapman-Enskog expansion exists everywhere and we could in principle put the boundary between the PDEand the Boltzmann model at any location, provided that we lift accurately. For general Boltzmannmodels, with complicated collision integral operators, such a Chapman-Enskog expansion mightnot exist everywhere in the domain. Then a hybrid model can be constructed where a PDE canreplace the Boltzmann model only in the regions where the Chapman-Enskog expansion is knownto exist.This situation appears in the modelling of laser ablation where a laser heats a surface thatconsequently melts and evaporates. The escaping plasma plume can be described by a Boltzmannequation. Close to the melting surface a complicated non-equilibrium situation appears whereescaping particles evaporate but particles that impinge on the melted surface condensate. Thereis no Chapman-Enskog expansion that can describe this situation close to the surface. Only awayfrom the surface the plasma reaches an equilibrium situation. A hybrid model will then use a fullBoltzmann model near melt while a reduced PDE model can be used away from the surface.27 cknowledgments This work is supported by research project
Hybrid macroscopic and microscopic modelling of laserevaporation and expansion , G.017008N, funded by ‘Fonds Wetenschappelijk Onderzoek’ togetherwith an ‘ID-beurs’ of the University of Antwerp.
References [1]
C.K. Aidun, J.R. Clausen , Lattice-Boltzmann Method for Complex Flows, Annu. Rev. Fluid Mech., 42(2010), pp. 439-472.[2]
M. Banda, A. Klar, L. Pareschi, M. Sea¨ıd , Lattice-Boltzmann type relaxation systems and high orderrelaxation schemes for the incompressible Navier-Stokes equations, Math. Comp., 77 (2008), pp. 943-965.[3]
R. Benzi, S. Succi, M. Vergassola , The lattice Boltzmann equation: theory and applications, Phys. Rep.,222 (1992), pp. 145-197.[4]
J.F. Bourgat, P. Le Tallec, B. Perthame, Y. Qiu , Coupling Boltzmann and Euler equations withoutoverlapping, Contemp. Math., 157 (1994), pp. 377-398.[5]
C. Cercignani , The Boltzmann equation and its applications, Springer, Berlin, 1988.[6]
S. Chapman, T.G. Cowling , The mathematical theory of non-uniform gases, Cambridge University Press,Cambridge, 1953.[7]
B. Chopard, M. Droz , Cellular Automata Modeling of Physical Systems, Cambridge University Press, Cam-bridge, 1998.[8]
P. Degond, G. Dimarco, L. Mieussens , A moving interface method for dynamic kinetic-fluid coupling,J. Comput. Phys., 227 (2007), pp. 1176-1208.[9]
P. Degond, G. Dimarco, L. Mieussens , A multiscale kinetic-fluid solver with dynamic localization of kineticeffects, J. Comput. Phys., 229 (2010), pp. 4907-4933.[10]
G. Dimarco, L. Pareschi , Hybrid multiscale methods II. Kinetic equations, Multiscale Model. Simul., 6(2008), pp. 1169-1197.[11]
A. Dupuis, E.M. Kotsalis, P. Koumoutsakos , Coupling lattice Boltzmann and molecular dynamics modelsfor dense fluids, Phys. Rev. E (3), 75 (2007), pp. 046704.[12]
E.G. Flekkøy, G. Wagner, J. Feder , Hybrid model for combined particle and continuum dynamics, EPL(Europhysics Letters), 52 (2000), pp. 271-276.[13]
A.L. Garcia, J.B. Bell, W.Y. Crutchfield, B.J. Alder , Adaptive mesh and algorithm refinement usingdirect simulation Monte Carlo, J. Comput. Phys., 154 (1999), pp. 134-155.[14]
C.W. Gear, T.J. Kaper, I.G. Kevrekidis, A. Zagaris , Projecting to a slow manifold: singularly perturbedsystems and legacy codes, SIAM J. Appl. Dyn. Syst., 4 (2005), pp. 711-732.[15]
M. Junk, A. Klar, L. Luo , Asymptotic analysis of the lattice Boltzmann equation, J. Comput. Phys., 210(2005), pp. 676-704.[16]
B.D. Kandhai , Large Scale lattice-Boltzmann simulations: computational methods and applications, PhDthesis, Universiteit van Amsterdam, 1999.[17]
I.G. Kevrekidis, C.W. Gear, J.M. Hyman, P.G. Kevrekidis, O. Runborg, C. Theodoropoulos ,Equation-free, coarse-grained multiscale computation: enabling microscopic simulators to perform system-level analysis, Commun. Math. Sci., 1 (2003), pp. 715-762.[18]
J. Latt, B. Chopard, P. Albuquerque , Spatial coupling of a lattice Boltzmann fluid model with a FiniteDifference Navier-Stokes solver, arXiv:physics/0511243v1 [physics.comp-ph], (2008), pp. 1-10.[19]
P. Le Tallec, F. Mallinger , Coupling Boltzmann and Navier-Stokes equations by half fluxes, J. Com-put. Phys., 136 (1997), pp. 51-67.[20]
M. Mehl, T. Neckel, P. Neumann , Navier-Stokes and Lattice-Boltzmann on octree-like grids in the Peanoframework, Internat. J. Numer. Methods Fluids, 65 (2011), pp. 6786. R. Mei, L. Luo, P. Lallemand, D. d’Humi`eres , Consistent initial conditions for lattice Boltzmann simula-tions, Comput. & Fluids, 35 (2006), pp. 855-862.[22]
M.L. Parks, R.B. Lehoucq , Atomistic-to-Continuum Coupling, SIAM NEWS, 24/9/2006.[23]
M. Rheinl¨ander , On the stability structure for lattice Boltzmann schemes, Comput. Math. Appl., 59 (2010),pp. 2150-2167.[24]
J.D. Sterling, S. Chen , Stability analysis of lattice Boltzmann methods, J. Comput. Phys., 123 (1996),pp. 196-206.[25]
S. Succi , The lattice Boltzmann equation for fluid dynamics and beyond, Oxford University Press, Oxford,2001.[26]
S. Succi, O. Filippova, G. Smith, E. Kaxiras , Applying the lattice Boltzmann equation to multiscale fluidproblems, Computing in Science & Engineering, 3 (2001), pp. 26-37.[27]
R. Van der Sman , Introduction to the Lattice Boltzmann method, University of Wageningen, 2004.[28]
P. Van Leemput , Multiscale and equation-free computing for lattice Boltzmann models, PhD thesis, K.U. Leu-ven, 2007.[29]
P. Van Leemput, C. Vandekerckhove, W. Vanroose, D. Roose , Accuracy of hybrid lattice Boltz-mann/finite difference schemes for reaction-diffusion systems, Multiscale Model. Simul., 6 (2007), pp. 838-857.[30]
P. Van Leemput, W. Vanroose, D. Roose , Mesoscale analysis of the equation-free Constrained Runs ini-tialization scheme, Multiscale Model. Simul., 6 (2007), pp. 1234-1255.[31]
C. Vandekerckhove , Macroscopic simulation of multiscale systems within the equation-free framework, PhDthesis, K.U. Leuven, 2008.[32]
C. Vandekerckhove, I. Kevrekidis, D. Roose , An efficient Newton-Krylov implementation of the Con-strained Runs scheme for initializing on a slow manifold, J. Sci. Comput., 39 (2009), pp. 167-188.[33]
Y. Vanderhoydonc, W. Vanroose , Lifting in hybrid lattice Boltzmann and PDE models, Comput. Vis. Sci.,14 (2011), pp. 67-78.[34]
D.A. Wolf - Gladrow , Lattice-gas cellular automata and lattice Boltzmann models, Springer, Berlin, 2000., Lattice-gas cellular automata and lattice Boltzmann models, Springer, Berlin, 2000.