Nonlinear fourth order Taylor expansion of lattice Boltzmann schemes

1. Introduction

The lattice Boltzmann schemes in their modern form have been developed with the contributions of d’Humières, Lallemand, Succi [13,37,43,50] and many others. An underlying lattice Boltzmann equation is discretised on a cartesian grid with a finite set of velocities chosen in such a way that during one time step, an exact transport is done between two vertices of the mesh. A lattice Boltzmann scheme is composed with two steps: a nonlinear local relaxation step, followed by a linear advection scheme coupling a given vertex with a given family of neighbours. The relaxation step follows in general the “BGK” approximation introduced by Bhatnagar, Gross and Krook [4]. This numerical method allows the simulation of an important number of physical phenomena as isothermal flows, compressible flows with heat transfer, non-ideal fluids, multiphase and multi-component flows, microscale gas flows, soft-matter flows,…up to quantum mechanics. For the actual status of the method and the various applications, we refer e.g. to the books of Guo and Shu [33], Krüger et al. [41] and to the prospective article of Succi [56]. The lattice Boltzmann method is founded on a mesoscopic Boltzmann model but is not able in general to solve the Boltzmann equation or associated kinetic models. It is admitted that essentially macroscopic models are approximated with the lattice Boltzmann schemes.

The link between mesoscopic and macroscopic models is never straightforward and an analysis is necessary. The usual method of formal analysis is founded on the Chapman–Enskog [8] expansion in the spirit of the book of Huang [38]. The classical approach for the Boltzmann equation and other mesoscopic models posed in a continuum space-time. The Chapman–Enskog method has been adapted to lattice Boltzmann schemes with discrete space and time equations by Chen and Doolen [9] and by Qian and Zhou [51]. It allows the derivation the equivalent partial differential equations founded on a multiscale time analysis. This classical Chapman Enskog expansion is popular in the lattive Boltzmann community, as reported in the contributions of Lallemand and Luo [43], Philippi and Hegele [48], Shan et al. [53] and Otomo et al. [46] among others. A main drawback of this approach is the fact that multiscale expansions are used without a clear mathematical signification of the various variables and associated functions.

Independently of this framework, we have proposed in [15,16] the Taylor expansion method for isothermal fluid problems to obtain formally equivalent partial differential equations at second order accuracy. This work was inspired by the method of equivalent equations for finite difference schemes, in a spirit proposed by Lerat and Peyret [45] and Warming and Hyett [57]. In this approach, the infinitesimal variable is simply the time step (proportional to the space step with the acoustic scaling) and there is only one infinitesimal scale for the analysis. This approach has been experimentaly validated with Pierre Lallemand in various contributions [22,23,25,42]. It was extended to third order accuracy for thermal and fluid problems in [17], automated at an arbitrary order in the linear case in [3], extended with an external drift in [26] and to relative velocities in [19,20]. Our point of view is a priori not to consider the third and fourth order terms as approximations of Burnett of super-Burnett equations (see e.g. [1,55]) but simply as errors of the lattice Boltzmann scheme when solving partial differential equations of second order type. Nevertheless, progress in the approximation of fourth order partial differential equations have been obtained in [47]. Here, we have no particular fluid hypothesis. It can be applied to scalar or to vectorial [18,30] lattice Bolzmann schemes. This Taylor expansion method has been also revisited with the Maxwell iterations of Yong et al. [58].

In this contribution, we first recall in Section 2 classic results related to the Boltzmann equation with discrete velocities. The link with multiple relaxation times lattice Boltzmann schemes as initially proposed by d’Humières [13] is presented in Section 3. Then we present the Taylor expansion method in a very general way and the result is explicited with compact formulas. Our main result (Section 4) concerns the derivation of asymptotic partial differential equations in the nonlinear case, up to fourth order accuracy. The proof for first and second orders is presented in Section 5, with a validation for the isothermal fluid case in the specific case of the D2Q9 scheme. In Section 6, the expansion up to third and fourth orders is detailed, with a confrontation with our previous work [17]. The proofs of these results need specific algebraic developments, presented in this contribution. The linear case with constant coefficients is described in Section 7.

2. Boltzmann equation with discrete velocities

In the space R d of dimension d, we consider a finite set of q discrete velocities v j ∈ V with components v j α for 1 ⩽ α ⩽ d . The unknowns of the Boltzmann equation are the particle densities f j . They are functions of space x, time t and discrete velocities v j : f j = f j ( x , t ) , x ∈ R d , t ⩾ 0 , 0 ⩽ j < q . The vector f ( x , t ) ∈ R q is constructed with the numbers f j ( x , t ) for 0 ⩽ j < q . We introduce a collision vector R q ∋ f ⟼ Q ( f ) ∈ R q of components Q j ( f ) for 0 ⩽ j < q with given regular functions Q j . We introduce also a small parameter ε > 0 that can be interpreted as a Knudsen number in the context of gas dynamics (see e.g. [8]).

The Boltzmann model with discrete velocities takes the form (1) ∂ f j ∂ t + ∑ 1 ⩽ α ⩽ d v j α ∂ f j ∂ x α = 1 ε Q j ( f ) , 0 ⩽ j < q . It has been proposed by Carleman [7], Gross [32], Broadwell [6] and intensively developed by Gatignol [28] and her co-workers.

∙  Moments

Our framework concerns multi relaxation times: we introduce a constant invertible matrix M called “d’Humières matrix” [13] in this contribution. This matrix defines the vector of moments m ∈ R q by a simple product: (2) m k ≡ ∑ 0 ⩽ j < q M k j f j . We introduce the number of conservation laws N ( 1 ⩽ N < q ) such that the N first moments of the collision kernel are equal to zero: (3) ∑ 0 ⩽ j < q M k j Q j ( f ) = 0 , ∀ f ∈ R q , 0 ⩽ k < N . Then it is natural to divide the vector of moments into two families: (4) m ≡ W Y . The conserved moments or macroscopic variables W constitute a linear space of dimension N and if the Boltzmann model with discrete velocities (1) is satisfied, we have N conservation laws due to the cancellation (3): ∂ W k ∂ t + ∑ 1 ⩽ α ⩽ d ∑ 0 ⩽ j < q M k j v j α ∂ f j ∂ x α = 0 , 0 ⩽ k < N . Observe that the nonconserved moments or microscopic variables Y generate a linear space of dimension q − N .

∙  Equilibrium states

We suppose that the equilibrim states f eq defined by the conditions Q ( f eq ) = 0 are characterized with the help of a regular nonlinear vector field Φ : R N ⟶ R q − N such that (5) f eq = M − 1 W Φ ( W ) . In other words, the vector field W ⟼ Y eq ≡ Φ ( W ) defines the set of equilibrium states.

We suppose moreover that the jacobian matrix d Q ( f eq ) at equilibrium is diagonalizable with real eigenvalues and real eigenvectors. More precisely, taking into account the hypothesis (3), we suppose that there exists a diagonal matrix Z = diag ( 1 τ 1 , … , 1 τ q − N ) of order q − N with strictly positive coefficients τ ℓ such that (6) M d Q ( f eq ) M − 1 = − 0 0 0 Z . In other words, the jacobian operator d Q ( f eq ) admits the matrix M − 1 as a matrix of right eigenvectors and the associated eigenvalues are all non positive.

∙  Momentum-velocity operator

For a linear space of dimension d, we introduce the momentum-velocity operator matrix Λ defined by the relation (7) Λ = M diag ( ∑ 1 ⩽ α ⩽ d v α ∂ α ) M − 1 . It is a q × q operator matrix composed by first-order space differential operators. It is obtained by conjugation of the first order advection operator v . ∇ by the d’Humières matrix M. The operator matrix Λ is nothing else than the advection operator seen in the basis of moments. We introduce a block decomposition of the momentum-velocity operator matrix associated to the decomposition (4) of the moments. We define a N × N operator matrix A, a N × ( q − N ) operator matrix B, a ( q − N ) × N operator matrix C and a ( q − N ) × ( q − N ) operator matrix D according to (8) Λ ≡ A B C D . Remember that in the following, the matrices A, B, C and D are matrices composed with first order space operators.

∙  Boltmann-BGK system with discrete velocities

We approximate now the previous discrete model with a BGK [4] type hypothesis: the state f is close to equilibrium and we approach the collision kernel Q by its first order expansion around the equilibrium state: Q ( f ) ≃ Q ( f eq ) + d Q ( f eq ) . ( f − f eq ) . Due to the condition Q ( f eq ) = 0 , we obtain with this approximation that we call “BGK” in this contribution, even it is not exactly the hypothesis done in the original article [4], a new system of partial differential equations: (9) ∂ f j ∂ t + ∑ 1 ⩽ α ⩽ d v j α ∂ f j ∂ x α = 1 ε d Q j ( f eq ) . ( f − f eq ) , 0 ⩽ j < q . Observe that f eq is a function of the particle distribution f through the relations (2), (4) and (5). Taking into account the hypothesis (6), the definition (7) and the notation (8), we can write the Boltmann-BGK system with discrete velocities (9) under the form of a system of two coupled partial differential equations: (10) ∂ W ∂ t + A W + B Y = 0 , ∂ Y ∂ t + C W + D Y = − 1 ε Z ( Y − Φ ( W ) ) .

In applications to fluid dynamics, the first order term Γ 1 ( W ) is associated to the Euler equations of gas dynamics and the second order term Γ 2 ( W ) to the viscous terms of the Navier–Stokes equations. The operators Γ 3 ( W ) and Γ 4 ( W ) represent the Burnett and super-Burnett terms respectively. With this Proposition 1, we essentially reformulate the classic Chapman Enskog expansion [8] in the discrete velocity case, as done e.g. in the book of Gatignol [28].

Proof.

Proof of Proposition 1 First we inject the representation (12) inside the first equation of (10). It comes ∂ W ∂ t + A W + B ( Φ ( W ) + ε Z − 1 Ψ 1 ( W ) + ε 2 Z − 1 Ψ 2 ( W ) + ε 3 Z − 1 Ψ 2 ( W ) ) = O ( ε 4 ) then after rearranging the terms, ∂ W ∂ t + A W + B Φ ( W ) + ε B Z − 1 Ψ 1 ( W ) + ε 2 B Z − 1 Ψ 2 ( W ) + ε 3 B Z − 1 Ψ 3 ( W ) = O ( ε 4 ) . The confrontation of this partial differential equation with the Ansatz (11) shows that Γ 1 ( W ) = A W + B Φ ( W ) Γ 2 ( W ) = B Z − 1 Ψ 1 ( W ) Γ 3 ( W ) = B Z − 1 Ψ 2 ( W ) Γ 4 ( W ) = B Z − 1 Ψ 3 ( W ) as proposed in (13).

Secondly, we write the second equation of (10) under the form Y = Φ − ε Z − 1 ( ∂ t Y + C W + D Y ) . We observe that − ∂ t Y = d Y . ( − ∂ t W ) = d ( Φ ( W ) + ε Z − 1 Ψ 1 ( W ) + ε 2 Z − 1 Ψ 2 ( W ) ) . ( Γ 1 + ε Γ 2 ( W ) + ε 2 Γ 3 ( W ) ) + O ( ε 3 ) = d Φ ( W ) . Γ 1 + ε ( Z − 1 d Ψ 1 ( W ) . Γ 1 + d Φ ( W ) . Γ 2 ) + ε 2 ( Z − 1 d Ψ 2 ( W ) . Γ 1 + Z − 1 d Ψ 1 ( W ) . Γ 2 + d Φ ( W ) . Γ 3 ) + O ( ε 3 ) . We report this expression inside the value of Y: Y = Φ + ε Z − 1 ( − ∂ t Y − C W − D Y ) = Φ + ε Z − 1 [ d Φ ( W ) . Γ 1 + ε ( Z − 1 d Ψ 1 ( W ) . Γ 1 + d Φ ( W ) . Γ 2 ) + ε 2 ( Z − 1 d Ψ 2 ( W ) . Γ 1 + Z − 1 d Ψ 1 ( W ) . Γ 2 + d Φ ( W ) . Γ 3 ) − C W − D ( Φ ( W ) + ε Z − 1 Ψ 1 ( W ) + ε 2 Z − 1 Ψ 2 ( W ) ) ] + O ( ε 4 ) = Φ + ε Z − 1 ( d Φ ( W ) . Γ 1 − C W − D Φ ) + ε 2 Z − 1 ( Z − 1 d Ψ 1 ( W ) . Γ 1 + d Φ ( W ) . Γ 2 − D Z − 1 Ψ 1 ) + ε 3 Z − 1 ( Z − 1 d Ψ 2 ( W ) . Γ 1 + Z − 1 d Ψ 1 ( W ) . Γ 2 + d Φ ( W ) . Γ 3 − D Z − 1 Ψ 2 ) + O ( ε 4 ) . Comparatively to the relation (12), we obtain Ψ 1 ( W ) = d Φ ( W ) . Γ 1 − ( C W + D Φ ) Ψ 2 ( W ) = Z − 1 d Ψ 1 ( W ) . Γ 1 + d Φ ( W ) . Γ 2 ( W ) − D Z − 1 Ψ 1 ( W ) Ψ 3 ( W ) = Z − 1 d Ψ 1 ( W ) . Γ 2 + d Φ ( W ) . Γ 3 − D Z − 1 Ψ 2 + Z − 1 d Ψ 2 ( W ) . Γ 1 and the last relations of (13) are explicited. The proof is completed. □

3. Multiple relaxation times lattice Boltzmann schemes

We introduce now a regular cartesian lattice L composed by vertices x separated by distances that are simple expressions of the space step Δ x . A discrete time t is supposed to be an integer multiple of a time step Δ t > 0 . The unknowns of a lattice Boltzmann scheme with q discrete velocities are the particle densities f j . They are now functions of discrete space x, discrete time t and discrete velocities v j for 0 ⩽ j < q : f j = f j ( x , t ) , x ∈ L , t = n Δ t , n ∈ N , v j ∈ V .

The discrete velocity set V does not depend on space or time.

Multiple relaxation times lattice Boltzmann scheme considered in this contribution consists in a simple Lie–Trotter splitting applied to the Boltzmann-BGK model with discrete velocities (9). First a relaxation step f ( t ) ⟶ R ( Δ t ) f ( t ) ≡ f ∗ for the approximation of the system of ordinary differential equations ∂ f ∂ t = 1 ε d Q ( f eq ) . ( f − f eq ) and secondly a transport step f ∗ ⟶ T ( Δ t ) f ∗ ≡ f ( t + Δ t ) for the free advection ∂ f j ∂ t + ∑ 1 ⩽ α ⩽ d v j α ∂ f j ∂ x α = 0 , 0 ⩽ j < q .

∙  Relaxation step

Due to the formulation (10), the relaxation scheme is simply described in the moment representation and we have ∂ W ∂ t = 0 , ∂ Y ∂ t = − 1 ε Z ( Y − Φ ( W ) ) . After relaxation, the vector of moments m is transformed into a new vector m ∗ : (14) m ∗ ≡ W ∗ Y ∗ . The conserved moments are not modified during this relaxation step: (15) W ∗ = W . As remarked in [16], we use an explicit first order Euler scheme for the microscopic moments: Y ∗ − Y Δ t = − 1 ε Z ( Y − Φ ( W ) ) , id est (16) Y ∗ = Y + S ( Φ ( W ) − Y ) with (17) S = Δ t ε Z . Usually, as pointed by Lallemand and Luo [43], the relaxation matrix S is a diagonal matrix, in coherence with the hypothesis concerning the diagonalization of the operator d Q ( f eq ) presented in the previous section: S = diag ( s 1 , … , s q − N ) . Moreover, the matrix S is a fixed invertible square matrix of dimension ( q − N ) × ( q − N ) . In other words, the time step Δ t is if the order of the Knudsen number ε in our asymptotic study. Finally the particle distribution f ∗ after relaxation is given according to (18) f ∗ = M − 1 m ∗ .

∙  Necessary conditions for stability

We deduce from the relation (16) Y ∗ − Φ ( W ) = ( I − S ) ( Y − Φ ( W ) ) and a natural stability condition takes the form ‖ I − S ‖ ⩽ 1 , id est 0 ⩽ s j ⩽ 2 , 0 ⩽ j < q − N . The case of over-relaxation 1 < s j < 2 is widely used when implementing lattice Boltzmann schemes to take into account very law viscosities in fluid mechanics. We will emphasize this point in the next section.

∙  Advection step

In the second step of advection, the method of characteristics is applied in the very particular case when it is exact. In one time step each particle distribution f j is exactly transported from the note x − v j Δ t of the lattice L to the vertex x ∈ L . In other words, the advection step occurs with a Courant–Friedrichs–Lewy number c f l ≡ 1 for all discrete velocities v j .

∙  Lattice Boltzmann scheme as a Lie–Trotter splitting scheme

The iteration of the lattice Boltzmann scheme is finally given by the relation (19) f ( t + Δ t ) = T ( Δ t ) R ( Δ t ) f ( t ) and it is usefull to explicit it for each component: (20) f j ( x , t + Δ t ) = f j ∗ ( x − v j Δ t , t ) , v j ∈ V , x ∈ L . Observe that this framework is very general. The so-called “BGK version” of the lattice Boltzmann scheme [5,37,50] corresponds to the choice M = I and Z = 1 τ I . Usual multi-relaxation lattice Boltzmann schemes [13,21,22,43] suppose in general some orthogonality property for the d’Humières matrix. This paradigm permits also the introduction of two particle distributions [2,21,34,40] or eventually more with the vectorial schemes developed by Graille [30] or a variant proposed in [18]. In those last cases, the numbering j ⟼ v j of the velocities is simply non injective.

It is naturally possible to think to a second order accurate splitting scheme. A natural idea is the Strang splitting [54]: f ( t + Δ t ) = R ( Δ t 2 ) T ( Δ t ) R ( Δ t 2 ) f ( t ) . In this case, as put in evidence by Dellar [11], one can write f ˜ ( t ) ≡ R ( Δ t 2 ) − 1 f ( t ) and the previous Strang scheme can be written as f ˜ ( t + Δ t ) = T ( Δ t ) R ( Δ t 2 ) R ( Δ t 2 ) f ˜ ( t ) . This scheme, up to the shift of one half relaxation step, is identical to the Lie–Trotter splitting (19), except that the product R ( Δ t 2 ) R ( Δ t 2 ) plays the role of R ( Δ t ) . A first possibility is to divide the time step Δ t by a factor 2 in the relations (16), (17): Y ‾ = Y + 1 2 S ( Φ − Y ) , Y ∗ = Y ‾ + 1 2 S ( Φ − Y ‾ ) . Then Y ∗ = ( 1 − S + 1 4 S 2 ) Y + ( S − 1 4 S 2 ) Φ ( W ) and over-relaxation is no more possible because s j − 1 4 s j 2 ⩽ 1 for all s j ∈ R . So such a numerical scheme is not used by the engineers at our knowledge.

A second possibility is to enforce the relation R ( Δ t 2 ) = R ( Δ t ) , id est Y ‾ = I − S Y + ( I − I − S ) Φ ( W ) . But this relation has a sense only if s j ⩽ 1 and over-relaxation is again excluded.

An other idea is to exchange the roles of advection and relaxation in the Strang scheme: f ( t + Δ t ) = T ( Δ t 2 ) R ( Δ t ) T ( Δ t 2 ) f ( t ) or to use more elaborated splittings as proposed by Drui et al. [14]. In this case, the transport step is nomore at a Courant number equal to unity and an interpolation procedure, costly in terms of numerical viscosity or complex to implement as the semi-Lagrangian method [27], is not usually chosen by experts developping lattice Boltzmann schemes.

In conclusion of this sub-section, the naive Lie–Trotter splitting scheme (19) gives very good numerical results. It is important to put in evidence precise qualities and defects of this scheme through the following formal asymptotic analysis.