François Dubois
email: francois.dubois@math.u-psud.fr

Pierre Lallemand
email: pierre.lallemand1@free.fr

On triangular lattice Boltzmann schemes for scalar problems

Keywords: Laplacian operator, heat equation, d'Humières scheme, D2T4, D2T7 AMS classification: 65-05, 65Q99, 82C20

We propose to extend the d'Humières version of the lattice Boltzmann scheme to triangular meshes. We use Bravais lattices or more general lattices with the property that the degree of each internal vertex is supposed to be constant. On such meshes, it is possible to define the lattice Boltzmann scheme as a discrete particle method, without need of finite volume formulation or Delaunay-Voronoi hypothesis for the lattice. We test this idea for the heat equation and perform an asymptotic analysis with the Taylor expansion method for two schemes named D2T4 and D2T7. The results show a convergence up to second order accuracy and set new questions concerning a possible super-convergence.

1) Introduction

The importance of extending the lattice Boltzmann scheme from square type regular meshes to unstructured triangulations has been recognized during the last years of 20th century [START_REF] Chen | Volumetric formulation of the lattice Boltzmann method for fluid dynamics: Basic concept[END_REF][START_REF] Karlin | Lattice Boltzmann Method for Irregular Grids[END_REF][START_REF] Peng | A Finite Volume Scheme for the Lattice Boltzmann Method on Unstructured Meshes[END_REF]. In particular the "volumetric formulation" of Chen [START_REF] Chen | Volumetric formulation of the lattice Boltzmann method for fluid dynamics: Basic concept[END_REF] makes a link with finite volumes, using control volumes around each vertex (the "Inria cells" [START_REF] Vijayasundaram | Transonic flow simulations using an upstream centered scheme of Godunov in finite elements[END_REF]) of a finite element type triangulation. This method is still under active development with the work of Succi, Ubertini and co-workers [START_REF] Pontrelli | The unstructured lattice Boltzmann method for non-Newtonian flows[END_REF][START_REF] Ubertini | Lattice Boltzmann method on unstructured grids: Further developments[END_REF][START_REF] Ubertini | Lattice Boltzmann schemes without coordinates[END_REF]. In a dual way, van der Sman [START_REF] Van Der Sman | Diffusion Lattice Boltzmann scheme on an Orthorhombic Lattice[END_REF][START_REF] Van Der Sman | Convection-Diffusion Lattice Boltzmann scheme for Irregular Lattices[END_REF][START_REF] Van Der Sman | Lattice Boltzmann scheme for Diffusion on Triangular Grids[END_REF][START_REF] Van Der Sman | Diffusion on unstructured triangular grids using Lattice Boltzmann[END_REF] uses rectangles and triangles as control volumes with a "cell center" type approach in Roache [START_REF] Roache | Computational Fluid Dynamics[END_REF] denomination. He has developed an approximation of diffusion equation with Delaunay-Voronoi meshes for a BGK variant of the lattice Boltzmann scheme.

In a previous contribution [START_REF] Dubois | On lattice Boltzmann scheme, finite volumes and boundary conditions[END_REF], we have observed that for usual lattice Boltzmann schemes (as for example the well known D2Q9), several (two for D2Q9) families of finite volumes are naturally associated with the scheme. As a consequence, we consider now the lattice Boltzmann scheme essentially as a "particle" method on a given (a priori fixed) mesh with discrete velocities. Recall that the "Particle In Cell" method has been first proposed in 1964 by Harlow et al. [START_REF] Harlow | The particle-in-cell computing method for fluid dynamics[END_REF] and has been analyzed in the eighties by Beale and Majda [START_REF] Beale | Vortex methods II: Higher order accuracy in two and three dimensions[END_REF], Raviart, Cottet and Mas Gallic [START_REF] Cottet | A particle method to solve the Navier-Stokes system[END_REF][START_REF] Gallic | A Particle Method for First order Symmetric Systems[END_REF][START_REF] Raviart | An analysis of particle method[END_REF] among others. We remark that this particle method does not suppose a priori the existence of a given lattice. The surrounding cells are recomputed at each time step in order to make the particle interact. Dynamic triangulation is an alternative to the previous methodology. It has been developed recently by Cianci, Klales, Love and co-workers [START_REF] Klales | Lattice gas simulations of dynamical geometry in two dimensions[END_REF][START_REF] Love | From the Boltzmann equation to fluid mechanics on a manifold[END_REF] in the context of lattice gas automata.

In this contribution, instead of adopting the volumetric formulation or a Delaunay-Voronoi hypothesis, we develop the framework of lattice Boltzmann schemes as a variant of the particle method. We propose an extension of the approach of d'Humières [START_REF] Humières | Generalized Lattice-Boltzmann Equations[END_REF] to triangular meshes and we restrict this first tentative to scalar problems like the heat equation without advection.

The outline of the contribution is the following. We first recall the classic D2T7 lattice Boltzmann scheme in the next section. At this occasion, we put in evidence a property of symmetry of Bravais lattices. It is possible to adapt the Taylor expansion analysis [START_REF] Dubois | Equivalent partial differential equations of a lattice Boltzmann scheme[END_REF] to this triangular lattice, with a diffusive scaling. This development is presented in Section 3 and applied to the D2T7 scheme. Several simulations with the D2T7 lattice Boltzmann scheme for the heat equation are presented in Section 4. In Section 5, we set the question of defining a discrete particle method on a finite element type triangular lattice. We propose a partial answer when each vertex of the lattice has a constant number of neighbours. This framework is applied in Section 6 to define a D2T4 lattice Boltzmann scheme for the heat equation. We repeat in Section 7 with this new scheme "D2T4" the simulations presented in Section 4. This work validates the potential of applications of our proposal. The conclusion (section 8) serves also as a discussion concerning encountered difficulties.



On triangular lattice Boltzmann schemes for scalar problems

2) D2T7 lattice Boltzmann scheme

We consider a Bravais lattice L connecting nodes labelled by the letter x and parametrized by a typical space scale ∆x. The neighbour vertex number j of the node x ∈ L is denoted by x j and we set [START_REF] Arnoldi | The principle of minimized iteration in the solution of the matrix eigenvalue problem[END_REF] x j = x + ξ j ∆x .

For each x ∈ L and each direction ξ j linking two vertices, the "opposite node" with number σ(j) defined according to

(2)

x σ(j) ≡ x -ξ j ∆x , ξ j + ξ σ(j) ≡ 0
is also a vertex of the lattice L (i.e. x σ(j) ∈ L). In the following, we emphasize this property satisfied by Bravais lattices and qualify it as a symmetric property. Most "DdQq" schemes (with a notation introduced by Qian et al. [START_REF] Qian | Lattice BGK Models for Navier-Stokes Equation[END_REF]) presented in the literature use a Bravais lattice. This symmetry property is also mandatory e.g. to define "two relaxation times" lattice Boltzmann schemes as proposed by Ginzburg et al. [START_REF] Ginzburg | Two-relaxation-time lattice Boltzmann scheme: about parametrization, velocity, pressure and mixed boundary conditions[END_REF].

The lattice Boltzmann scheme with multiple relaxation times is defined in a classical manner. Consider a vertex x that belongs to the lattice L. Then the j o direction of propagation is defined with a vector ξ j and ξ j ∈ V, set of directions that define the vicinity of the vertex x. The j o density of particles at vertex x and time t is denoted by f j (x, t). After a local step of relaxation, the j o density of particles is named f * j (x, t). Because a Bravais lattice L is symmetric, the neighbouring vertex x σ(j) defined in (2) in the direction opposite to the j o direction of propagation belongs to the lattice L. The lattice Boltzmann scheme can be completely defined:

(3) f j (x, t + ∆t) = f * j (x -ξ j ∆x, t) . Moreover the basic iteration (3) of a lattice Boltzmann scheme supposes explicitly that the lattice is symmetric, as illustrated in Figure 1 (left). Typical stencil of a lattice Boltzmann scheme for a Bravais lattice (left) ; both opposite directions ξ j and -ξ j connect two vertices of the mesh. Local numbering of the six neighbours (right) of the D2T7 lattice Boltzmann scheme on triangles.

The D2T7 lattice Boltzmann uses equilateral triangles as suggested by Frisch, Hasslacher and Pomeau in 1986 [START_REF] Frisch | Lattice gas automata for the Navier Stokes equation[END_REF] in the context of lattice gas automata. We precise the parameters that we have to consider. A vertex x has a total of six neighbours (seven including itself) ξ j (j = 0, . . . , 6) as in Figure 1 (right). Following d'Humières approach [START_REF] Humières | Generalized Lattice-Boltzmann Equations[END_REF] we introduce moments m k as linear functions of the particle distribution f :

(4) m k = j M kj f j .


François Dubois and Pierre Lallemand

We restrict our study to the simple case of only one conservation (thermal problem). Following [START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability[END_REF], we introduce a family P of polynomials p k for k = 0, • • • , 6:

(5)

P = 1, X, Y, X 2 + Y 2 , 4 √ 3 X Y, 2 (X 2 -Y 2 ), 3 Y -4 Y 3 .
The coefficients of the matrix M introduced at relation ( 4) are simply given by a nodal value in the velocity space:

(6)

M kj = p k (ξ j ) , 0 ≤ j, k ≤ 6 .
We remark that M αj = ξ α j for α = 1, 2. We have only one conserved moment ρ ≡ m 0 = m eq 0 = m * 0 = j f j and the other moments at equilibrium follow the relations m eq 1 = m eq 2 = 0, m eq 3 = a 3 ρ, m eq 4 = m eq 5 = m eq 6 = 0. The relaxation of moments out of equilibrium is also very simple:

(7) m * k = m k + s k (m eq k -m k ) , k = 1, • • • , 6 .
with s 1 = s 2 and s 4 = s 5 to enforce isotropy.

3) Taylor expansion with diffusive scaling

We can analyse the D2T7 lattice Boltzmann scheme with the Taylor expansion method [START_REF] Dubois | Equivalent partial differential equations of a lattice Boltzmann scheme[END_REF]. We consider one time step of iteration (3) and we replace the particle distribution in the right hand side by the moments after relaxation:

f j (x, t + ∆t) = ℓ M -1 jℓ m * ℓ (x -ξ j ∆x, t).
In consequence, we have the formal expansion in the moment space

m k (x, t + ∆t) = jℓ M kj M -1 jℓ m * ℓ (x -ξ j ∆x, t) = jℓ M kj M -1 jℓ m * ℓ (x, t) -ξ α j ∆x ∂ α m * ℓ + O(∆x 2 ) = m * k -∆x ℓ j M kj ξ α j M -1 jℓ ∂ α m * ℓ + O(∆x 2 ).
We introduced the momentum-velocity tensor introduced in [START_REF] Dubois | Equivalent partial differential equations of a lattice Boltzmann scheme[END_REF]:

Λ ℓ kp ≡ j M kj M pj M -1 jℓ .
Then we have up to third order accuracy

m k (x, t + ∆t) = m * k -∆x Λ ℓ kα ∂ α m * ℓ + 1 2 ∆x 2 Λ p kα Λ ℓ pβ ∂ α ∂ β m * ℓ + O(∆x 3
) and using the relaxation step [START_REF] Ciarlet | General Lagrange and Hermite interpolation in IR n with applications to finite element methods[END_REF],

m * k (x, t + ∆t) = m eq k -∆x 1 -s k s k Λ ℓ kα ∂ α m eq ℓ + O(∆x 2
). We adopt the so-called "diffusive scaling" proposed initially for rarefied flows by Sone [START_REF] Sone | Asymptotic theory of flow of rarefied gas over a smooth boundary I[END_REF] (see an explicit derivation for lattice Boltzmann schemes e.g. in Junk et al. [START_REF] Junk | Asymptotic analysis of the lattice Boltzmann equation[END_REF]) [START_REF] Cottet | A particle method to solve the Navier-Stokes system[END_REF] ∆t

≡ ∆x 2 ζ 
where ζ is a constant for homogeneity of dimensions. We add some advection term by enforcing the relations according to m eq 1 = u ∆x ζ and m eq 2 = v ∆x ζ . After some pages of formal calculus, following the method presented in details in [START_REF] Dubois | Towards higher order lattice Boltzmann schemes[END_REF], we obtain the equivalent partial differential equation : [START_REF] Humières | Generalized Lattice-Boltzmann Equations[END_REF] ∂ρ ∂t

+ u ∂ρ ∂x + v ∂ρ ∂x -µ ∆ρ = Θ ∆x 2 ∆ 2 ρ + ∆x 4 A 6 ρ + O(∆x 6 ) .
Up to second order accuracy, we have an approximation of the heat equation with a diffusivity coefficient µ given according to

µ = 1 2 ζ a 3 σ 1 .
The coefficients σ k for the nonconserved moments are given by the Hénon's relation [START_REF] Hénon | Viscosity of a Lattice Gas[END_REF]

σ k ≡ 1 s k -1 2 .
The coefficient Θ in front of the fourth order term in ( 9) is explicited as follows for u = v = 0 :

(10) Θ = - 1 16 σ 1 a 3 ζ (1 -a 3 ) 1 -4 σ 1 σ 3 -2 σ 1 σ 4 + 4 a 3 σ 2 1 .
In the relation ( 9), A 6 is a sixth order operator. The development of the other moments can also be achieved. In particular, we have

m α = m eq α -a 3 2 s 1 ∆x ∂ α ρ + O(∆x 2 ).
• "Second order", "quartic" and "hexahedric" coefficients We have chosen the following numerical values ζ = 1 , a 3 = 1 4 , s 1 = 0.8 compatible with a diffusivity coefficient µ = 0.09375. In these conditions, the D2T7 lattice Boltzmann scheme is formally equivalent to the heat equation up to order 2 (id est, due to (9) and [START_REF] Dubois | Equivalent partial differential equations of a lattice Boltzmann scheme[END_REF], Θ = 0 and A 6 = 0) when using to fix the ideas the following "second order" coefficients (given here with 15 decimals for a possible implementation): (11) s 3 = 1.428571428571428, s 4 = s 5 = 0.481927710843373, s 6 = 0.476190476190476.

With the following choice of "quartic" relaxation coefficients (12) s 3 = 1.428571428571428, s 4 = s 5 = 0.930232558139534, s 6 = 0.526315789473684, the D2T7 lattice Boltzmann scheme is formally of the order 4 (id est Θ = 0 and A 6 = 0). Last but not least, we can impose Θ ≡ 0 and A 6 ≡ 0 and the D2T7 scheme is of order 6. The "hexahedric" coefficients can be taken as follows: (13) s 3 = 1.086117521785847, s 4 = s 5 = 1.344205296559553, s 6 = 0.647305233773416.

4) Diffusion simulations with the D2T7 scheme

We have done several simulations: a "one point" periodic analysis, a numerical evaluation of the modes for a periodic pipe and a rectangle, the computation of harmonic functions by time asymptotics of the heat equation, the dissipation of a triangular Dirichlet mode and the direct numerical computation of triangular Dirichlet modes.

• One point periodic analysis

The one point analysis can be conducted as follows. We start from the iteration (3) of a lattice Boltzmann scheme. We suppose that the particle field for the neighbouring points of vertex x satisfy the following periodicity condition :

(14) f j x -ξ j ∆x, t) = exp -i k • ξ j ∆x f j x, t)
for some wave vector k = (k cos θ, k sin θ). From ( 14), the evaluation of the right hand side of ( 3) is easy in the context of the d'Humières version of the lattice Boltzmann scheme. The state vector f is then solution of an eigenvalue problem of small dimension q for a general lattice Boltzmann problem with q velocities. In the D2T7 case for thermal problems, we obtain six eigenvalues λ ℓ ≃ 1 -s ℓ for ℓ ≥ 1 and one physical eigenvalue

λ(k) ≃ 1 -µ k 2 .
This eigenvalue has a real meaning for applications to macroscopic physics. A numerical diffusivity µ num ≡ (1 -λ(k))/k 2 can be extracted from the previous relation. In Figure 2, we have plotted the error ǫ ≡ | µ -µ num | as a function of the modulus of the wave vector. With the three versions of the D2T7 scheme detailed in ( 11), ( 12) and ( 13), the errors for the diffusivity have an order of convergence directly predicted by the Taylor expansion analysis. • Periodic pipe and rectangle The analysis for a periodic pipe is conducted by following the same ideas. A D2T7 lattice Boltzmann solver is considered on a simple geometry of nx ≡ 96 by ny ≡ 4 mesh points.

The unknown is now a vector f ∈ IR 7 nx ny . The iteration of the scheme defines a linear operator A and the first eigenvalue of this operator is determined thanks to an Arnoldi algorithm [START_REF] Arnoldi | The principle of minimized iteration in the solution of the matrix eigenvalue problem[END_REF]. The first eigenvalue λ ≡ 1 corresponds to the conservation of mass in the whole domain, including boundary conditions. The second eigenvalue λ 0 corresponds to the smallest wave vector compatible with the computational domain. It is compared with the modulus of the wave vector to evaluate a numerical diffusivity µ num = (1 -λ 0 )/k 2 as  On triangular lattice Boltzmann schemes for scalar problems previously. The different errors ǫ ≡ | µ -µ num | are presented in Figure 3 (left). The first two versions (11)(12) of orders two and four present a coherent numerical convergence. The results are not so clear with the sixth order tuning of the parameters. It seems to be due to the round-off errors for this study involving three orders of magnitude for wave vector.

The analysis is analogous for a rectangle nx ≡ 36 by ny ≡ 52 mesh points. The results are depicted in Figure 3 (right). The lattice Boltzmann scheme has a coherent order of convergence for the "second order" and "fourth order" versions of the scheme. The "sixth order" scheme exhibits now an error numerically evaluated as only fifth order accurate. This fact seems again to due to round-off errors in the Arnoldi process [START_REF] Arnoldi | The principle of minimized iteration in the solution of the matrix eigenvalue problem[END_REF]. • Harmonic polynomials on a triangle We have developed a D2T7 solver for a triangular geometry (see a typical mesh in Figure 4)

The initial condition is a priori identically null. We determine the numerical boundary conditions compatible with a polynomial expression p H (x, y) ≡ x 2 -y 2 on the boundary with an "anti-bounce-back" version of the algorithm of Bouzidi et al. [START_REF] Bouzidi | Momentum transfer of a Boltzmann-lattice fluid with boundaries[END_REF]. The computation converges in time towards the harmonic function introduced above. We present in Figure 4 (right) the numerical result ρ(x, y) ≃ p H (x, y) when we use 61 points on the edge of the triangle (that corresponds to a total of 1891 vertices for the entire mesh). In Figure 5, we have plotted the error field for the three versions (11)(12)(13) of the D2T7 lattice Boltzmann scheme. The results are qualitatively coherent: the more the scheme is theoretically precise, the more the error is reduced. In Figure 6 (left), we observe that the  On triangular lattice Boltzmann schemes for scalar problems L ∞ error is substantially reduced when the parameters induce a better precision. But the order of convergence remains very close to second order even for "fourth order" and "sixth order" versions of the scheme using the set of parameters [START_REF] Dubois | On lattice Boltzmann scheme, finite volumes and boundary conditions[END_REF] or [START_REF] Dubois | Towards higher order lattice Boltzmann schemes[END_REF]. In this case this default can be due to a possible deficit of time steps and to crude boundary conditions. The lattice Boltzmann scheme is explicit and the time iterations (see Figure 6, right) take too much time to reach convergence to the stationary state with a satisfactory reduction of the error.

• Dissipation of a triangular Dirichlet mode

We have also experimented the relaxation of a Dirichlet mode. The first mode is simply a product of three "sinus" functions, as first explicited by Lamé (see McCartin [START_REF] Mccartin | Eigenstructure of the Equilateral Triangle, Part I: The Dirichlet Problem[END_REF]). With our nomenclature, the eigenvalue number ℓ is proportional to 3 (ℓ -1) 2 . The reference value is in consequence equal to 12, 48 and 108 for ℓ equal to 3, 5 and 7 respectively. We present in Figure 7 the results at T = 4/3 and the evolution of the physical field at the center. The asymptotic analysis obtained by successive mesh refinements is presented in Figure 8. We measure the error in time for the center vertex as the mesh size tends to zero and the L ∞ error at the precise time T = 4/3 in the same conditions. The results are correct but not easy to interpret. The "second order" scheme is just a bit better that the order 3/2. The "fourth order" version is of order 3 and the "sixth order" scheme hesitates between the orders 3 and 4. Negative values in blue and positive ones in red. Time and L ∞ space errors for several meshes and several "orders" with the D2T7 lattice Boltzmann scheme. The obtained accuracy is not the one proposed by the Taylor expansion method. The space numerical accuracy is going from 1.7 to 3.0 with a good tuning of the numerical parameters.

• Dirichlet modes for a triangle We used the D2T7 lattice Boltzmann scheme (3) to define a linear operator f (t) -→ f (t + ∆t) ≡ A•f (t) where f (t) is the vector of all unknowns for the entire mesh. Then the first eigenvalues of the linear operator A are computed with the Arnoldi algorithm [START_REF] Arnoldi | The principle of minimized iteration in the solution of the matrix eigenvalue problem[END_REF].

Figure 9. Isovalues of the Dirichlet mode number "3" of an equilateral triangle (left) and errors for a D2T7 computation. The exact reference eigenvalue is equal to 12 (in appropriate units). The numerical eigenvalue is equal to 11.99902 with second order parameters (middle) and to 11.99938 with "fourth order" parameters (right). The L ∞ error for the modes is equal to 4 10 -4 at order 2 and 10 -4 at order 4. The figures show the isovalues of the error for both computations with different scales. We observe that the global shape of these errors is similar to isovalues of the reference eigenvector. Some exact reference modes are displayed in the left part of Figures 9 to 11. The numerical approximation is globally of very good quality and we have plotted the errors for different  modes computed on the same lattice in the same figures. We perform the computations for each mode, one with the "second order" accurate version of the D2T7 scheme and the other one with a "fourth order" accurate tuning of numerical parameters. In each case, we compare the theoretical eigenvalue after applying a suitable normalization and the computed eigenvalue by the Arnoldi algorithm. The results are of good quality and the quartic parameters give a better precision for the numerical results. Even if the fourth order convergence is not established, the tuning of parameters improves clearly the numerical quality. Dirichlet mode number "7" of an equilateral triangle (left) and errors for a D2T7 computation. The exact reference eigenvalue is equal to 108. The numerical eigenvalue is equal to 107.90777 at order 2 (middle) and to 107.92705 at order 4 (right). The L ∞ error for the modes with this computation is equal to 1.02 10 -2 at order 2 and 4.2 10 -3 at order 4.

5) Lattice Boltzmann scheme on arbitrary meshes ?

Imagine that we move the vertices in the Bravais lattice presented at Figure 4 (left). We obtain a topologically regular mesh in the sense that the number of edges containing a given internal vertex is constant. An example is proposed at Figure 12. This mesh is a good candidate for future extensions of the lattice Boltzmann scheme. With this kind of classical finite element type mesh, it is possible to use all the engineering tools of automatic meshing in two and three space dimensions as described e.g. in [START_REF] George | Automatic Mesh Generation: Applications to Finite Element Methods[END_REF]. But this goal is still not the purpose of the present contribution. The vertices of the mesh of Figure 12 are now the nodes of a cellular complex and each vertex has a constant number of neighbours. In other terms, the degree of each vertex is constant. We denote by x j ≡ x + ξ j (x) ∆x the vertex belonging to the lattice L with a local neighbouring number j relative to the vertex x. Remark that the vertex x is also a neighbour of the vertex x j with a local number ℓ ≡ n j (x). We have the obvious relation ξ j (x) + ξ ℓ (x j ) ≡ 0 and in other terms the identity [START_REF] George | Automatic Mesh Generation: Applications to Finite Element Methods[END_REF] ξ j (x) + ξ n j (x) (x j ) ≡ 0 .

As previously, we denote by f j (x) the density of particles going from vertex x towards vertex x j . Moreover, the outgoing particles from vertex x j are also ingoing particles "into" vertex x with an index denoted by ℓ. We precise the previous notation. If f j (x) is the density of particles from the vertex x towards the vertex x j , f * j (x) denotes the same quantity after relaxation. In a dual vision, we denote by f * ℓ (x j ) the density of particles going from the vertex x j in the direction of the vertex x after relaxation. We have also to consider the density f j (x) of particles going from the vertex x j towards the vertex x. The lattice Boltzmann scheme is a particle method. The flight of particles between the vertex x j and the vertex x takes exactly one time step : f j (x, t + ∆t) = f * ℓ (x j , t). If we replace the notation ℓ for the index of vertex x relative to its neighbour x j by the notation n j (x) introduced previously at relation [START_REF] George | Automatic Mesh Generation: Applications to Finite Element Methods[END_REF], the lattice Boltzmann scheme takes the form [START_REF] Ginzburg | Two-relaxation-time lattice Boltzmann scheme: about parametrization, velocity, pressure and mixed boundary conditions[END_REF] f j (x, t + ∆t) = f * n j (x) (x j , t) .

In the case of general meshes, the relation ( 16) replaces the initial formula (3), correct only for Bravais lattices, as illustrated in Figure 13.

x + ξ ∆ x = x x j j j n (x) * f (x ) j f (x) j Figure 13
. Iteration of a lattice Boltzmann scheme: the ingoing particles f j (x) into vertex x are coming from the neighbouring vertex x j after a relaxation step.

We precise now how to compare the ingoing particles f j (x), the outgoing particles f * j (x) emitted from the vertex x and the associated moments. We introduce first a matrix M (x) in order to compute the moments m k (x) from the ingoing particles. As previously (see the relation ( 5)), we suppose given a family P of polynomials p k . In an analogous way suggested by the relation ( 6), we just reverse the direction of velocities and we have M (x) kj = p k -ξ j , with p k ∈ P. If the polynomials 1, X and Y are the first polynomials of the family P, we have as in the previous studies

M (x) 0j = 1 , M (x) αj = -ξ α j (x) , 1 ≤ α ≤ d.
The moments are evaluated for the incoming particles with the natural relation [START_REF] Griffiths | On the scope of the method of modified equations[END_REF] m k (x)

≡ j M (x) kj f j (x) , 0 ≤ k ≤ q -1 , x ∈ L .
The relaxation step is essentially unchanged. The moments can be seen as the eigenvectors of the jacobian of the relaxation matrix (see e.g. [START_REF] Dubois | Introduction au Schéma de Boltzmann sur Réseau[END_REF]) and this operator is diagonal with this representation:

(18) m * k (x) = m k (x) + s k (x) m eq k (x) -m k (x)
, where the index k in relation ( 18) is running on all nonconserved moments. The outgoing particles after relaxation are supposed to be a linear functional of the moments:

(19) f * j (x) ≡ k P (x) jk m * k (x) , 0 ≤ j ≤ q -1 , x ∈ L .
The question is now to determine the matrix P . We have the following property.

Proposition 1.

Transition matrix from moments to particle distribution. If the Taylor expansion approach is valid at the order zero and if each internal node of lattice L is of constant degree (the number of neighbours of each vertex does not depend of the vertex x ∈ L), the matrix P (x) of relation ( 19) is given by the relation [START_REF] Junk | Asymptotic analysis of the lattice Boltzmann equation[END_REF] P (x

) i ℓ = M (x i ) -1 n i (x) ℓ . Proof of Proposition 1.
The proof can be conducted as follows. We start from the time iteration (3) of the lattice Boltzmann scheme. Then after multiplication by the matrix M (x), with the help of ( 17), ( 16) and ( 19), we have

m k (x, t + ∆t) = j M (x) kj f j (x, t + ∆t) = j M (x) kj f * n j (x) (x j , t)
 François Dubois and Pierre Lallemand

= j M(x) kj ℓ P (x j ) n j (x) ℓ m * ℓ (x j , t) = ℓ j M(x) kj P (x j ) n j (x) ℓ m * ℓ (x+ξ j ∆x, t).
We expand this relation at order one. Due to relaxation, we just have a small perturbation between m and m * :

m k (x) + O(∆t) = m * k (x) + O(∆x). In consequence, j M (x) kj P (x j ) n j (x) ℓ ≡ δ k ℓ
and in other terms,

P (x j ) n j (x) ℓ = M (x) -1 j ℓ .
We change the names of the vertices. We replace the letter x j by the letter x. Then we replace the index n j (x) by some neighbor i of vertex x and the index j is now equal to n i (x). With this change of notation, we obtain P (x) i ℓ = M (x i ) -1 n i (x) ℓ which is exactly the relation [START_REF] Junk | Asymptotic analysis of the lattice Boltzmann equation[END_REF].

We can now make explicit the d'Humières lattice Boltzmann scheme on an arbitrary mesh where the degree of each vertex is constant. When all the outgoing densities of particles f * j (x, t) are known for all the vertices of the lattice at some discrete time t, the ingoing densities f j (x, t + ∆t) at the new time step are simply evaluated by a free flight [START_REF] Ginzburg | Two-relaxation-time lattice Boltzmann scheme: about parametrization, velocity, pressure and mixed boundary conditions[END_REF] during one time step. Then the moments m k are a local linear transform of the particle densities thanks to [START_REF] Griffiths | On the scope of the method of modified equations[END_REF]. The first moments compose a set W (x) of conserved variables and the equilibrium moments m eq are a given (in general nonlinear) function G(W ) of this field: m eq k (x) = G k (W (x)), x ∈ L. The relaxation of moments follow the relation [START_REF] Harlow | The particle-in-cell computing method for fluid dynamics[END_REF]. Note that in general the coefficients s k (x) now depend a priori explicitly on the vertex x. Last but not least, the outgoing particles at the new time step from the vertex x follow the local linear transform [START_REF] Hénon | Viscosity of a Lattice Gas[END_REF].

6) D2T4 scheme for equilateral triangles

We consider a general two-dimensional mesh L composed by triangles. Note here that a cellular complex is composed by "vertices" in L 0 of dimension zero, by edges in L 1 of dimension one and by triangles of dimension two: x ∈ L 2 . In other words, we adopt a "cell center" framework in the sense proposed in Roache [START_REF] Roache | Computational Fluid Dynamics[END_REF]. We can also locate the degree of freedom x at the center of gravity of the corresponding triangle. Remark that we make here a priori no other regularity hypothesis. Each triangle x has three edges. Each edge inside the border of x is part of the boundary of (at most) two triangles : the triangle x itself and its j th neighbor x j . It is then natural to consider outgoing particles (f j ) 0≤j≤4 going from x towards x j with a local velocity ξ j (x) ∆t chosen in such a way that the centers of both triangles x and x j are joined in exactly one time step of duration ∆t. Of course, the null velocity is not excluded. This remark explains the name "D2T4" of this type of lattice Boltzmann scheme. A typical regular mesh for a D2T4 computation is presented in Figure 14. Typical mesh with equilateral triangles. The four degrees of freedom of D2T4 scheme are located at the center of gravity of each triangle. The links between triangles create the dual hexagonal mesh around the vertices of the triangular mesh.

The degrees of freedom in Figure 14 are the centers of the initial triangular mesh. This "secondary mesh" is no longer a Bravais lattice. We lose the possibility of straight propagation of particles in the lattice and also the symmetry property of Bravais meshes emphasized in Figure 1. But we keep the property that the number of neighbours is constant. And this property is maintained whatever the initial triangulation with cellular complexes. Two types of triangles for the D2T4 scheme with equilaterals. The local numbers are explicited for each edge.

We observe that such a lattice contains only two types of equilateral triangles: the "left" and "right" ' triangles as displayed in Figure 15. We precise now the choices we have done to construct our scheme. The family P of polynomials is simply composed by a restriction of (5) to the first four terms: P = 1, X, Y, X 2 + Y 2 . Because we have two generic triangles, we have two families of neighboring directions

ξ left = -1, 0 , 1 2 , - √ 3 2 , 1 2 , √ 3
2 , ξ right = 1, 0 , -1 2 , √ 3 2 , -1 2 , - √ 3
2 . We observe also that due to the simple numbering of local edges (see the figure 15), we have the simple relations n j (x) ≡ j and ξ left j +ξ right j = 0. In this contribution, we consider only one conserved variable ρ = m 0 ≡ 3 j=0 f j . The moments at equilibrium are simply chosen with m eq 1 = m eq 2 = 0 and m eq 3 = a 3 ρ. Proposition 2.

Transition matrix for the D2T4 lattice Boltzmann scheme For the D2T4 lattice Boltzmann scheme defined previously, we have [START_REF] Junk | Weighted L 2 -Stability of the Lattice Boltzmann Method[END_REF] P left = M left -1 , P right = M right -1 . In this particular case, the relations [START_REF] Junk | Weighted L 2 -Stability of the Lattice Boltzmann Method[END_REF] are exactly analogous to the ones for lattice Boltzmann schemes on Bravais lattices. In some sense, for the D2T4 scheme, the relations (21) remain (too !) simple !

Proof of Proposition 2

Recall that due to [START_REF] Junk | Asymptotic analysis of the lattice Boltzmann equation[END_REF], we have P (x) i ℓ = M (x i ) -1 n i (x) ℓ with n i (x) ≡ i due to our precise choice of numbering (see the figure 15). Then we have the two matrix equalities • Taylor expansion analysis for the D2T4 scheme The analysis can now be conducted without difficulty in the same framework than previously. We adopt the diffusive-scaling [START_REF] Cottet | A particle method to solve the Navier-Stokes system[END_REF]. After some developments with the help of formal calculus (see e.g. [START_REF] Dubois | Towards higher order lattice Boltzmann schemes[END_REF]) we derive the equivalent partial differential equation at the order 6:

P left i ℓ = M right
(22)    ∂ρ ∂t -µ ∆ρ = a 3 ζ 24 (12 σ 2 1 -1) ∆x ∂ 2 x -3 ∂ 2 y ∂ x ρ + Θ 2 ∆x 2 ∆ 2 ρ + Θ 3 ∆x 3 ∂ 2 x -3 ∂ 2 y ∆ ∂ x ρ + ∆x 4 A 6 ρ + O(∆x 6
) . The notation σ k is identical to the one used at Hénon's relation [START_REF] Hénon | Viscosity of a Lattice Gas[END_REF]. The diffusion coefficient µ satisfies the relation µ = ζ a 3 σ 1 .

• "First order", "second order", "third order" and "quartic" coefficients We have chosen ζ = 1. For first order simulations, we have taken the following numerical values [START_REF] Klales | Lattice gas simulations of dynamical geometry in two dimensions[END_REF] a 3 = 0.216506350946109, s 1 = 1.2, s 3 = 0.750796078775233 compatible with a diffusion coefficient µ = 0.0721687836487032 = 1 4 √

12 . With the choice σ 1 = 1 √ 12 the scheme is at least second order accurate (see the right hand side of ( 22)) and we take parameters to fit the previous choice of the diffusion coefficient:

id est σ 3 = √ 3
2 , we have Θ 2 = Θ 3 = 0. With these conditions, the D2T4 scheme is theoretically fourth order accurate.



On triangular lattice Boltzmann schemes for scalar problems

7) Diffusion simulations with the D2T4 scheme

We have done essentially the same simulations as performed with the D2T7 lattice Boltzmann scheme (see Section 4).

• One point periodic analysis The results are presented in Figure 16. The theoretical orders with the four choices of parameters proposed previously are exactly the one proposed by the Taylor expansion analysis. A defect of isotropy for the numerical diffusivity is clearly visible for parameters that lead to a first order and third order schemes with this D2T4 simulator. 23), ( 24), ( 25), ( 26) and ( 27) lead to schemes of several orders. The measured orders with a linear regression are displayed in the right column.

• Periodic pipe and rectangle

We have tested the fourth order version (24) (27) of the D2T4 lattice Boltzmann scheme on two simple periodic geometries presented in Section 4. The numerical results (Figure 17) show that the scheme is convergent, but simply at second order accuracy. The predicted coefficients for the order 4 define a second order scheme in this particular case.

• Harmonic polynomials on a triangle

The numerical computation of Laplace equation with non-homogeneous boundary conditions has been also performed by integrating the heat equation and taking the limit for time large enough. Our simulation (Figure 18) shows that the D2T4 scheme is convergent with second order accuracy. Nevertheless the second order and fourth order versions of the scheme give essentially the same results. Figure 18. Two-dimensional computation of the harmonic function p H (x, y) = x 2 -y 2 on a triangle with the D2T4 lattice Boltzmann scheme (left). Convergence of the L ∞ error for several meshes. The D2T4 lattice Boltzmann scheme remains of order 2 even if quartic parameters are used in the simulation. Dissipation of the first Dirichlet mode for an equilateral triangle (right). Time and L ∞ space errors for several meshes and several "orders". The numerical accuracy is equal to 2 for all the parameters. With "quartic" parameters the absolute level of the error is substantially reduced.

• Dissipation of a triangular Dirichlet mode

The dissipation of the first mode described at the figure 7 for the D2T7 scheme has been constructed without difficulty. Now the two main versions of the scheme (second and fourth orders) converge with second order accuracy as shown in Figure 18. We observe that even if no extra order of convergence has been obtained, the results with quartic parameters give a better precision.

• Dirichlet modes for a triangle

The simulations done with the D2T7 lattice Boltzmann scheme have been compared with a D2T4 simulator. The results (Figure 19) explicit this comparison. A first result is that the level of error for D2T4 is comparable with D2T7 results at order two. If we look precisely to the error fields, distinguished contribution is due to the boundary conditions.

8) Conclusion

We have proposed an extension of the lattice Boltzmann method for triangular meshes. Our first step concerns a single conservation law and we made numerical simulations for the heat equation. For an extension of the discrete particle method, we have considered non Bravais lattices such that the degree of each vertex is constant. Our formulation does not need any finite volume or Delaunay-Voronoi triangulation hypothesis as in the previous contributions. We have used the Taylor expansion analysis with a diffusive scaling to explicit some parameters of the d'Humières scheme. With this method, it is possible to get formally a better accuracy. Our simulations show that this extra accuracy can be obtained with very fundamental one point periodic hypothesis. In more realistic cases, this extra-accuracy is in general not observed.

We think that triangular meshes explicit the limit of validity of the Taylor expansion analysis. In fact when we write the lattice Boltzmann scheme with the relation (3) or ( 16) and when we perform the Taylor expansion, we suppose that there exists a very regular function f (x, t) of space and time that support the definition of the scheme. In particular, this function is supposed to be independent of the lattice ! This last Ansatz is in defect for triangular meshes on nonsymmetric lattices as D2T4. Note that this kind of remark recover other critics [START_REF] Chang | A critical analysis of the modified equation technique of Warming and Hyett[END_REF][START_REF] Griffiths | On the scope of the method of modified equations[END_REF] relative to this kind of symptotic analysis [START_REF] Lerat | Noncentered Schemes and Shock Propagation Problems[END_REF][START_REF] Warming | The modified equation approach to the stability and accuracy analysis of finite difference methods[END_REF].

Two directions of research are natural in the continuation of the present contribution. First we can try to develop a true mathematical analysis of the lattice Boltzmann scheme, following e.g. previous work of Junk and Yong [START_REF] Junk | Weighted L 2 -Stability of the Lattice Boltzmann Method[END_REF] with appropriate mathematical tools, as done typically by Ciarlet and Raviart for finite elements [START_REF] Ciarlet | General Lagrange and Hermite interpolation in IR n with applications to finite element methods[END_REF] or Gallouët and coworkers for finite volumes [START_REF] Champier | Convergence of an Upstream Finite Volume Scheme for a Nonlinear Hyperbolic Equation on a Triangular Mesh[END_REF]. Second we can extend triangular lattice Boltzmann schemes to systems with other conservation laws for acoustics and fluid flow applications, revisiting the breakthrough of Frisch, Hasslacher and Pomeau [START_REF] Frisch | Lattice gas automata for the Navier Stokes equation[END_REF]. Preliminary results have been obtained for D2T10, which are not described in this article due to unnecessary complications.

  Figure 1.Typical stencil of a lattice Boltzmann scheme for a Bravais lattice (left) ; both opposite directions ξ j and -ξ j connect two vertices of the mesh. Local numbering of the six neighbours (right) of the D2T7 lattice Boltzmann scheme on triangles.

Figure 2 .

 2 Figure 2. D2T7 lattice Boltzmann scheme for the heat equation. One point periodic analysis. Error ǫ ≡ | µ -µ num | between exact and numerical diffusivities. Note that the exact dispersion equation can be obtained, and when solved by successive approximations in powers of k, it leads to the same results, obtained from the successive equivalent equations.

Figure 3 .

 3 Figure 3. D2T7 lattice Boltzmann scheme for the heat equation. Periodic modes for a pipe with nx = 96 and ny = 4 mesh points (left) an periodic modes for a rectangle of nx = 36 by ny = 52 points (right). Error ǫ ≡ | µ -µ num | between exact and numerical diffusivities. The hexahedric predicted coefficients define a fourth or fifth order scheme.

Figure 4 .

 4 Figure 4. Typical two-dimensional mesh for a D2T7 computation on a triangle (left). Two-dimensional computation of the harmonic function p H (x, y) = x 2 -y 2 on a triangle with the D2T7 lattice Boltzmann scheme (right). The iso-contours are composed by discrete hyperbolas.

Figure 5 .

 5 Figure 5. Two-dimensional computation of the harmonic function p H (x, y) = x 2 -y 2 on a triangle with the D2T7 lattice Boltzmann scheme. Iso-contours of the errors for three sets of parameters presented at relations[START_REF] Dubois | Introduction au Schéma de Boltzmann sur Réseau[END_REF],[START_REF] Dubois | On lattice Boltzmann scheme, finite volumes and boundary conditions[END_REF] and[START_REF] Dubois | Towards higher order lattice Boltzmann schemes[END_REF]. Negative values are in blue and positive ones in red. The maximal errors are equal to 8.14 10 -4 (left), 2.36 10 -4 (middle) and 4.47 10 -5 (right) when using a D2T7 scheme with formal order of 2 (left), 4 (middle) and 6 (right).

Figure 6 .

 6 Figure 6. Two-dimensional computation of the harmonic function p H (x, y) = x 2 -y 2 on a triangle with the D2T7 lattice Boltzmann scheme. No extra order is observed for the L ∞ error when refining the mesh (left). The evolution in time is very slow (right), even initializing the computation with the exact solution !

Figure 7 .FrançoisFigure 8 .

 78 Figure 7. Isovalues of the first Dirichlet mode for an equilateral triangle (left). Dissipation of this mode by time evolution: D2T7 solution at time T = 4/3 for 76 points on the edge (2926 vertices, middle). Exponential decay at the center of the mesh (55 vertices, right).

Figure 10 .

 10 Figure10. Isovalues of the Dirichlet mode number "5" of an equilateral triangle (left) and errors for a D2T7 computation. The exact reference eigenvalue is equal to 48. The numerical eigenvalue is equal to 47.98339 at order 2 (middle) and to 47.98842 at order 4 (right). The L ∞ error for the modes is equal to 3. 10 -2 at order 2 and 1.1 10 -3 at order 4.

Figure 11 .

 11 Figure 11.Dirichlet mode number "7" of an equilateral triangle (left) and errors for a D2T7 computation. The exact reference eigenvalue is equal to 108. The numerical eigenvalue is equal to 107.90777 at order 2 (middle) and to 107.92705 at order 4 (right). The L ∞ error for the modes with this computation is equal to 1.02 10 -2 at order 2 and 4.2 10 -3 at order 4.

Figure 12 .

 12 Figure 12. Triangular lattice obtained from a little random displacement of the vertices of an equilateral triangular mesh.

Figure 14 .

 14 Figure 14.Typical mesh with equilateral triangles. The four degrees of freedom of D2T4 scheme are located at the center of gravity of each triangle. The links between triangles create the dual hexagonal mesh around the vertices of the triangular mesh.

  Figure 15.Two types of triangles for the D2T4 scheme with equilaterals. The local numbers are explicited for each edge.

-1 i ℓ and P right i ℓ =) . But ξ left j + ξ right j = 0

 ℓ0 M left -1 i ℓ . We remark also that M left kj = p k (ξ left j ) and M left kj = p k (-ξ left j ). Analogously M right kj = p k (ξ right j ) and M right kj = p k (-ξ right j , then P left = M right -1 = (M left ) -1 and for the other family of triangles P right = M left -1 = (M right ) -1 . The relation (21) is established.

( 24 )

 24 a 3 = 0.25, s 1 = 1.267949192431122 With the particular value (25) s 3 = 0.422649730810374 , we have Θ 2 = 0 and the D2T4 scheme is formally second order accurate. With (26) s 3 = 0.758775495823486 , we have Θ 2 = 0, Θ 3 = 0 and the D2T4 scheme is formally third order accurate. With the choice

Figure 16 .

 16 Figure[START_REF] Ginzburg | Two-relaxation-time lattice Boltzmann scheme: about parametrization, velocity, pressure and mixed boundary conditions[END_REF]. One point periodic analysis. Error ǫ ≡ | µ -µ num | between numerical and theoretical diffusivities. Four sets of parameters defined at relations (23), (24), (25), (26) and (27) lead to schemes of several orders. The measured orders with a linear regression are displayed in the right column.

Figure 17 .

 17 Figure 17.D2T4 lattice Boltzmann scheme for the heat equation. Periodic modes for a pipe with nx = 96 and ny = 4 mesh points (left) and a rectangle of nx = 36 by ny = 52 points (right). Error ǫ ≡ | µ -µ num | between exact and numerical diffusivities. The predicted coefficients for the order 4 define a second order scheme in this particular case.

Figure 19 .

 19 Figure 19. Isovalues of the errors for D2T4 computation of Dirichlet modes of an equilateral triangle. The exact reference eigenvalues are equal respectively to 12, 48 and 108 and the computed ones to 11.97493, 47.59816 and 105.95870. This figure can be compared to second order accurate D2T7 results at Figures 9, 10 and 11 respectively.

Acknowledgments

The authors thank the "LaBS project" (Lattice Boltzmann Solver, www.labs-project.org), funded by the French "FUI8 research program", for supporting this contribution. Last but not least, the authors thank the referees for very constructive remarks. Some of them have been incorporated into the present edition of the article.


François Dubois
email: francois.dubois@math.u-psud.fr

Pierre Lallemand
email: pierre.lal@free.fr

Mahdi Tekitek
email: mohamed-mahdi.tekitek@math.u-psud.fr

On a superconvergent lattice Boltzmann boundary scheme

Keywords: Lattice Boltzmann scheme, boundary conditions, Taylor expansion method AMS classification: 65-05, 65Q99, 82C20

 Ginzburg and Adler analyzed the bounce-back boundary conditions for the lattice Boltzmann scheme and showed that it could be made exact to second order for the Poiseuille flow if some expressions depending upon the parameters of the method were satisfied, thus defining so-called "magic parameters". Using the Taylor expansion method that one of us developed, we analyze a series of simple situations (1D and 2D) for diffusion and for linear fluid problems using bounce-back and "anti bounce-back" numerical boundary conditions. The result is that "magic parameters" depend upon the detailed choice of the moments and of their equilibrium values. They may also depend upon the way the flow is driven.

Introduction

The theoretical analysis of the lattice Boltzmann scheme [START_REF] Frisch | Lattice gas hydrodynamics in two and three dimensions[END_REF][START_REF] Higuera | Boltzmann approach to lattice gas simulations[END_REF][START_REF] Higuera | Lattice gas dynamics with enhanced collisions[END_REF][START_REF] Qian | Lattice BGK models for Navier-Stokes equation[END_REF][START_REF] Humières | Generalized Lattice-Boltzmann Equations[END_REF][START_REF] Karlin | Maximum Entropy Principle for Lattice Kinetic Equations[END_REF][START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability[END_REF] is an active subject of research. Recall that the method was first analyzed by d'Humières [START_REF] Humières | Generalized Lattice-Boltzmann Equations[END_REF] with a Chapman Enskog expansion coming from statistical physics; we also refer to Asinari and Ohwada [START_REF] Asinari | Connection between kinetic methods for fluid-dynamic equations and macroscopic finite-difference schemes[END_REF] for a method of analysis based on the Grad moment system. A fruitful idea followed by Junk et al [START_REF] Junk | Asymptotic analysis of the lattice Boltzmann equation[END_REF] and [START_REF] Dubois | Une introduction au schéma de Boltzmann sur réseau[END_REF][START_REF] Dubois | Equivalent partial differential equations of a lattice Boltzmann scheme[END_REF] is to use the so-called equivalent equation method derived independently by Lerat-Peyret [START_REF] Lerat | Noncentered Schemes and Shock Propagation Problems[END_REF] and Warming and Hyett [START_REF] Warming | The modified equation approach to the stability and accuracy analysis of finite difference methods[END_REF] (see also [START_REF] Villatoro | On the method of modified equations. V: Asymptotic analysis of and direct-correction and asymptotic successive-correction techniques for the implicit midpoint method[END_REF]). An infinitesimal parameter is introduced and the finite differences operators are expanded into a family of equivalent partial differential equations. The main goal of this study is to use the Taylor expansion method [START_REF] Dubois | Une introduction au schéma de Boltzmann sur réseau[END_REF][START_REF] Dubois | Equivalent partial differential equations of a lattice Boltzmann scheme[END_REF] in order to increase the accuracy of boundary conditions for simple problems with analytical solutions. We first consider a one-dimensional (1D) diffusion problem and study the influence of the definition of the moments and of their equilibrium value. We then consider a two-dimensional (2D) Poiseuille flow using several ways to enforce a pressure gradient.

We consider regular lattices parametrized by a space step ∆x. We introduce a time step ∆t and adopt the "acoustic" scaling: the ratio λ ≡ ∆x ∆t is a fixed reference velocity for each study. As a consequence, the parameters ∆x and ∆t are equivalent infinitesimals. Note that as this work is devoted to boundaries, we shall use a particular way to test the accuracy of a numerical scheme as will be discussed later.

Diffusion problem in one space dimension

We consider the classical Lattice Boltzmann model D1Q3 with three discrete velocities and one conservation law to model diffusion problems. We choose the velocities v i (0 ≤ i ≤ 2) such that v 0 = 0, v 1 = λ, v 2 = -λ. At each mesh point, there are three functions {f j } that can be interpreted as populations of fictitious particles. These populations evolve according to the lattice Boltzmann scheme which we write as in [START_REF] Dubois | Une introduction au schéma de Boltzmann sur réseau[END_REF]:

(1)

f j (x, t + ∆t) = f * j (x -v j ∆t, t) , 0 ≤ j ≤ 2,
where the superscript * denotes post-collision quantities and x a vertex of the lattice. Therefore during each time increment ∆t there are two fundamental steps: advection and collision. The advection step describes the motion of a particle which has undergone collision at node x -v j ∆t and goes to the jth neighbouring node. Following d'Humières [START_REF] Humières | Generalized Lattice-Boltzmann Equations[END_REF], the collision step is defined in the space of moments. For D1Q3 three moments {m ℓ } are obtained by a linear transformation of vectors f j :

(2)

m 0 = f 0 + f 1 + f 2 ≡ ρ (density), m 1 = λ(f 1 -f 2 ), m 2 = λ 2 2 (f 1 + f 2 ).
In consequence, we introduce a matrix of moments M to represent moments like (2); it takes the form

(3) M =    1 1 1 0 λ -λ 0 λ 2 2 λ 2 2
and the relations (2) can be simply written as m = M f . To simulate diffusion problems, we conserve only the density moment ρ in the collision step and obtain one macroscopic scalar equation. The other quantities (non-conserved moments) are assumed to relax towards equilibrium values (m eq 1 , m eq 2 ) following:

(4) m * ℓ = (1 -s ℓ ) m ℓ + s ℓ m eq ℓ , 1 ≤ ℓ ≤ 2,
where s ℓ (0 < s ℓ < 2, for ℓ = 1, 2) are relaxation rates, not necessarily equal to a single value as in the BGK case [START_REF] Qian | Lattice BGK models for Navier-Stokes equation[END_REF]. The equilibrium values m eq ℓ of the non conserved moments in equation (4) determine the macroscopic behavior of the scheme. Indeed with the following choice of equilibrium values (neglecting non-linear contributions):

(5)

m eq 1 = 0 , m eq 2 = ζ λ 2 2 ρ
and using the Taylor expansion method we find (see e.g. [START_REF] Dubois | Third order equivalent equation of lattice Boltzmann scheme[END_REF]) that the equivalent partial differential equation of the numerical scheme up to order three in ∆x is a diffusion equation:

(6) ∂ρ ∂t -κ ∂ 2 ρ ∂x 2 = O(σ 1 ∆x 3 ) .
The value of the diffusivity κ is given according to

(7) κ = ∆t λ 2 σ 1 ζ
where σ ℓ ≡ 1 s ℓ -1 2 , ℓ = 1, 2. Remark that the thermal diffusivity κ is imposed by the Physics. Moreover the scale velocity λ is fixed and the coefficient ζ is also imposed. When we refine the mesh, the coefficient σ 1 must be chosen in order to enforce relation [START_REF] Frisch | Lattice gas hydrodynamics in two and three dimensions[END_REF]. In other terms, the product σ 1 ∆t must be maintained constant. Then the right hand side of relation (6) exhibits a second order truncation error of the lattice Boltzmann scheme for a given thermal diffusivity κ. Associated with stability properties (see Junk and Yong [START_REF] Junk | Weighted L 2 Stability of the Lattice Boltzmann Method[END_REF]), convergence properties of lattice Boltzmann scheme can be established, as in [START_REF] Junk | Convergence of Lattice Boltzmann Methods for Navier-Stokes Flows in Periodic and Bounded Domains[END_REF].

Localization of a one-dimensional boundary

Let us introduce a constant c and consider the following one-dimensional Poisson problem:

(8) -K d 2 ρ dx 2 = c in ]0, 1[, ρ(0) = ρ(1) = 0 .
We take an "anti bounce-back" numerical boundary condition at x = 0: to model the right hand side c of equation [START_REF] Ginzburg | Generic boundary conditions for lattice Boltzmann models and their application to advection and anisotropic dispersion equations[END_REF]. So we can write the lattice Boltzmann scheme as follows: (i) m = M f, (ii) m 0 = m 0 + 1 2 δp, (iii) evaluate the other moments, (iv) relaxation (4) of the other moments, (v) m 0 = m 0 + 1 2 δp, (vi) f = M -1 m, (vii) advection step (1) and boundary conditions. The exact solution of problem (8) is elementary:

(9) f 1 (x b , t + ∆t) = -f 2 (x e , t + ∆t) = -f * 2 (x b , t) ,
u(x) = c x (1-x) 2 K .
We analyze the behavior of the discrete model to show whether it can be tuned so that the location of the "numerical boundary" can be fixed at mid-point as expected from "anti bounce-back". Thus we shall use as criterion for accuracy the difference between the imposed boundary and the "numerically determined" boundary.

From a theoretical point of view, we suppose that the discrete fields f j (x, t) vary slowly in space and time in order to be able to use Taylor expansions. We analyse the lattice Boltzmann scheme in terms of equivalent partial differential equations and formal developments. It is well known (see e.g. Griffiths and Sanz-Serna [START_REF] Griffiths | On the scope of the method of modified equations[END_REF] or Chang [START_REF] Chang | A critical analysis of the modified equation technique of Warming and Hyett[END_REF]) that this method of analysis fails a priori to predict boundary effects properly if this hypothesis is not satisfied. We keep in mind this restriction in our numerical experiments. Nevertheless, this elementary tool can produce nontrivial results, as we will see hereafter.

We say in the following that a boundary scheme (such as [START_REF] Ginzbourg | Boundary flow condition analysis for three-dimensional lattice Boltzmann model[END_REF] to fix the ideas) is of order p at location ∆q relative to homogeneous Dirichlet boundary condition (that are present in (8)) if the numerical boundary condition implies [START_REF] Ginzburg | Second order boundary method for Lattice Boltzmann model[END_REF] ρ

(x b -∆q) = O(∆x p )
for the continuous conserved field issued from the particle field f j (x, t) according to (2).

We have the following result:

Proposition 1.
Superconvergent relation between parameters For the D1Q3 lattice Boltzmann scheme (1) (2) (4) (5), the "anti bounce-back" numerical boundary condition ( 9) is of order 3 at location ∆q = ∆x 2 relative to the homogeneous Dirichlet boundary condition of problem [START_REF] Ginzburg | Generic boundary conditions for lattice Boltzmann models and their application to advection and anisotropic dispersion equations[END_REF] if and only if the following condition [START_REF] Ginzburg | Multireflection boundary conditions for lattice Boltzmann models[END_REF] σ

1 σ 2 = 1 8 is satisfied.
Relation [START_REF] Ginzburg | Multireflection boundary conditions for lattice Boltzmann models[END_REF] defines superconvergent parameters σ 1 and σ 2 . Recall that they have been called "magic" by I. Ginzbourg and P.M. Alder [START_REF] Ginzbourg | Boundary flow condition analysis for three-dimensional lattice Boltzmann model[END_REF].

Proof of Proposition 1.

We have introduced in [START_REF] Dubois | Une introduction au schéma de Boltzmann sur réseau[END_REF][START_REF] Dubois | Equivalent partial differential equations of a lattice Boltzmann scheme[END_REF] the "tensor of momentum velocities" and the so-called "defects of conservation" which are defined respectively by

Λ ℓ kα ≡ j M kj M αj (M -1 ) jℓ (12) θ k ≡ ∂ t m eq k + Λ ℓ kα ∂ α m eq ℓ , k ≥ 1.
For the D1Q3 lattice Boltzmann scheme applied to diffusion problem like ( 8), ( 13)

θ 1 = ζλ 2 ∂ρ ∂x , θ 2 = ζ λ 2 2 ∂ρ ∂t .
Then we obtain the following development of non-equilibrium moments at third order (as described in [START_REF] Dubois | Third order equivalent equation of lattice Boltzmann scheme[END_REF]):

(14) m * k = m eq k + ∆t 1 2 -σ k θ k -∆t(σ k ∂ t θ k + σ ℓ Λ ℓ kα ∂ α θ ℓ ) + O(∆x 3 ), k ≥ 1.
Thus for k = 1:

m * 1 = ∆t 1 2 -σ 1 θ 1 -∆t(σ 1 ∂ t θ 1 + σ 1 Λ 1 11 ∂ x θ 1 + σ 2 Λ 2 11 ∂ x θ 2 ) + O(∆x 3 ) .
Using (13) and ∂ t θ 1 = O(∆t 2 ), Λ 1 11 = 0 the above equation becomes:

(15) m * 1 = ∆tλ 2 1 2 -σ 1 ζ ∂ρ ∂x + O(∆x 3 ).
For k = 2, we use expression (5) of m eq 2 , together with θ 2 = 0 and Λ 1 21 = λ 2 2 to obtain from equation ( 14):

(16) m * 2 = λ 2 ζ 2 ρ -∆t 2 λ 4 ζ 2 σ 1 1 2 -σ 2 ∂ 2 ρ ∂x 2 + O(∆x 3 )
Using the inverse moment matrix M -1 , the post-collision f are given by:

f * 1 = 1 2λ 2 [2m * 2 + λm * 1 ] , f * 2 = 1 2λ 2 [2m * 2 -λm * 1 ] . ( 17
)
At the boundary, due to ( 9) and (1), we consider the following quantity:

(18) f * 1 (x e ) + f * 2 (x b ) = 1 2λ 2 [2(m * 2 (x e ) + m * 2 (x b )) + λ(m * 1 (x e ) -m * 1 (x b ))] .
Using relations [START_REF] Humières | Generalized Lattice-Boltzmann Equations[END_REF] and (15) we obtain respectively:

(19)      m * 2 (x e ) + m * 2 (x b ) = λ 2 ζ 2 [ρ(x e ) + ρ(x b )] -∆t 2 λ 4 ζ 2 σ 1 1 2 -σ 2 ∂ 2 ρ ∂x 2 (x e ) + ∂ 2 ρ ∂x 2 (x b ) + O(∆x 3 )  (20) m * 1 (x e ) -m * 1 (x b ) = ∆t λ 2 ζ ( 1 2 -σ 1 ) ∂ρ ∂x (x e ) - ∂ρ ∂x (x b ) + O(∆x 3 ).
With the help of classical Taylor expansion we have, with the notation

x i ≡ 1 2 (x b + x e ): ρ(x e ) + ρ(x b ) = 2ρ(x i ) + ∆x 2 4 ∂ 2 ρ ∂x 2 (x i ) + O(∆x 3 ), (21) ∂ρ ∂x (x e ) - ∂ρ ∂x (x b ) = -∆x ∂ 2 ρ ∂x 2 (x i ) + O(∆x 3 ), (22)
∂ 2 ρ ∂x 2 (x e ) + ∂ 2 ρ ∂x 2 (x b ) = 2 ∂ 2 ρ ∂x 2 (x i ) + O(∆x 2 ). ( 23
)
Considering equation [START_REF] Junk | Convergence of Lattice Boltzmann Methods for Navier-Stokes Flows in Periodic and Bounded Domains[END_REF], together with [START_REF] Junk | Weighted L 2 Stability of the Lattice Boltzmann Method[END_REF], [START_REF] Karlin | Maximum Entropy Principle for Lattice Kinetic Equations[END_REF] and taking into account relations ( 21), ( 22) and ( 23) we obtain:

(24) f * 1 (x e ) + f * 2 (x b ) = ζ ρ(x i ) + ζ ∆x 2 σ 1 σ 2 - 1 8
∂ 2 ρ ∂x 2 (x i ) + O(∆x 3 ).
Due to the simple fundamental expression (1) of a lattice Boltzmann scheme, the left hand side of ( 24) is identically null when the numerical boundary condition (9) occurs. Due to the relation

x i = 1 2 (x b + x e )
, the condition ( 10) is satisfied with ∆q = ∆x 2 and p = 3 if and only if σ 1 σ 2 = 1 8 .

• Let us now consider the effect of using a different moment matrix for the D1Q3 case:

(25) M =   1 1 1 0 λ -λ -2λ 2 λ 2 λ 2   .
obtained from (3) by a Gram-Schmidt orthogonalization algorithm as usual with the lattice Boltzmann scheme (see e.g. [START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability[END_REF]). The moments at equilibrium are now given by m eq 1 = 0 and m eq 2 = λ 2 ζ ρ. We remark that the matrix of moments [START_REF] Warming | The modified equation approach to the stability and accuracy analysis of finite difference methods[END_REF] leads to an equivalent macroscopic conservation law of type (6) with a diffusivity κ which is now given by κ = ∆t

λ 2 σ 1 2+ ζ 3 . Proposition 2.
Third order at the boundary The D1Q3 lattice Boltzmann scheme (1) (25) (4) (5) associated to the "anti bounce-back" numerical boundary condition (9) is of order 3 at location ∆q = ∆x 2 for the homogeneous Dirichlet boundary condition of problem [START_REF] Ginzburg | Generic boundary conditions for lattice Boltzmann models and their application to advection and anisotropic dispersion equations[END_REF] if and only if σ 1 σ 2 = 3 8 . Remark that the superconvergent parameters satisfying the relation [START_REF] Ginzburg | Multireflection boundary conditions for lattice Boltzmann models[END_REF] emerging from Proposition 1 with the choice of transformation matrix M given by (3) are different from those obtained in the case with matrix [START_REF] Warming | The modified equation approach to the stability and accuracy analysis of finite difference methods[END_REF].



On a superconvergent lattice Boltzmann boundary scheme

Proof of Proposition 2.

For this model we have the following Taylor development of non-conserved moments up to order 2 on ∆x:

m * 1 = ∆t λ 2 1 2 -σ 1 2 + ζ 3 ∂ρ ∂x + O(∆x 3 ), m * 2 = λ 2 ζ ρ -∆t 2 λ 4 σ 1 1 2 -σ 2 2 + ζ 3
∂ 2 ρ ∂x 2 + O(∆x 3 ).
With the help of the matrix moments [START_REF] Warming | The modified equation approach to the stability and accuracy analysis of finite difference methods[END_REF] we have:

f * 1 = 1 3 ρ + 1 6λ 2 m * 2 + 3λm * 1 , f * 2 = 1 3 ρ + 1 6λ 2 m * 2 -3λm * 1 .
As f 2 (x b ) is internal to the domain we add to ρ(x b ) a body source δρ = -2+ ζ 3 ∂ 2 ρ ∂x 2 . Now by using the same method as in the proof of Proposition 1, we obtain:

f * 1 (x e ) + f * 2 (x b ) = 2 + ζ 3 ρ(x i ) + ∆x 2 8σ 1 σ 2 -3 2 + ζ 72 ∂ 2 ρ ∂x 2 (x i ) + O(∆x 3 ).
The conclusion is a direct consequence of the above calculus. To illustrate the preceeding discussion, we perform a numerical simulation of the two lattice Boltzmann models and analyze (after a suitable number of iterations to reach steady state) the "Poiseuille" parabolic profile. We measure the numerical error in terms of a precise location of the boundary for Dirichlet type boundary condition. We follow  François Dubois, Pierre Lallemand and Mahdi Tekitek a method proposed by Ginzburg and d'Humières [START_REF] Ginzburg | Second order boundary method for Lattice Boltzmann model[END_REF]: from the numerical discrete field u LB (j ∆x) we determine by least squares a parabola that fit at best the data. Then we calculate where this approximation of the numerical solution u LB is equal to zero. We interpret this location as the "experimental numerical location" of the solid wall. We find experimentally that the extrapolated location of the Dirichlet boundary condition is located between x b and x e and this exact solid wall location is parametrized under the form x b -∆q, with 0 ≤ ∆q ≤ ∆x. The results obtained for several values of σ 1 and σ 2 are shown in Figure 2 to depend only upon the product σ 1 σ 2 and go through 0 respectively for 1 8 or 3 8 , in complete coherence with the Taylor expansion method developed in Propositions 1 and 2.

The two-dimensional Poiseuille flow

We consider here the classical D2Q9 model (see e.g. [START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability[END_REF]). We study a Poiseuille flow (in linear regime), first with an imposed uniform body force and periodic boundary condition at the inlet and oulet of the channel. Then we consider the same flow with an imposed difference of pressure between inlet and outlet. The evolution of the lattice Boltzmann scheme is given by equation [START_REF] Asinari | Connection between kinetic methods for fluid-dynamic equations and macroscopic finite-difference schemes[END_REF]. The corresponding moments have an explicit physical significance: m 0 ≡ ρ is the density, m 1 ≡ j x and m 2 ≡ j y are x and y components of momentum, m 3 is the energy, m 4 is related to square energy, m 5 , m 6 are x and y components of heat flux and m 7 , m 8 are diagonal stress and off-diagonal stress. A Gram-Schmidt orthogonalization method is also used and the matrix of moments is exactly that used in [START_REF] Dubois | Une introduction au schéma de Boltzmann sur réseau[END_REF][START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability[END_REF]. The collision is described in the moments space as:

(26) m * ℓ = (1 -s ℓ ) m ℓ + s ℓ m eq ℓ , 3 ≤ ℓ ≤ 8,
where the equilibrium values m eq k are given by:

(27) m eq 3 = α ρ , m eq 4 = β ρ , m eq 5 = -j x λ , m eq 6 = -j y λ , m eq 7 = 0 , m eq 8 = 0 .

 • The Poiseuille flow

We introduce a two-dimensional domain Ω =]0, L[×]0, H[ (see Figure 3). Let u(t, x, y) ≡ (u, v) be the velocity of fluid and p the pressure solution of the "Poiseuille" Stokes system:

(28) -ν ∆u + 1 ρ ∇p = 0 , div u = 0 in Ω
with the following boundary conditions:

(29) p(0, y) = -p(L, y) = δp for 0 ≤ y ≤ H , u(x, 0) = u(x, H) = 0 for 0 ≤ x ≤ L .
The solution of the above problem is classical: u(x, y) ≡ (u(y) = K y (H -y), v = 0) , p(x, y) = 2 ρ ν K x+ P 0 , where P 0 is a given constant. We note that the problem (28) ( 29) is equivalent to a flow resulting from the action of a constant external force F between two arallel walls with periodic boundary conditions on the inlet and the outlet of the channel (i.e. in the Ox direction). So the problem (28) (29) becomes:

(30) -ν ∆u = F , u(x, 0) = u(x, H) = 0,
where F = (F x , 0) is the external force. The solution is given by u

(x, y) = u(y) = Fx 2 ν y (H -y) , v = 0 .
• A first lattice Boltzmann scheme We use the D2Q9 lattice Boltzmann scheme to model the Poiseuille flow described by equation (30). We use the equilibrium (27) with s 7 = s 8 = 1 2 + 3ν λ 2 ∆t -1 to have the exact viscosity ν present in equation (30). The implementation of the lattice Boltzmann scheme is conducted as follows for an arbitrary mesh vertex x of the lattice.

At initial time step t = t 0 we set the vectors f (x, t 0 ) = 0. For any given time t, we first determine the moments m k using the relation m ≡ M f. Then we change velocity j x before the collision step by adding a half of the external force F x :  x = j x + ∆t 2 F x . Thus the macroscopic moments (density and velocity) are evaluated. Then we perform the collision step in moments space according to relation (26) and we add half of the external force F x to the conserved velocity  x : j x (t + ∆t) =  x + ∆t 2 F x . Using the matrix M -1 we compute the particle distributions f * α (x, t). We perform advection through a relation analogous to (1) and we obtain the vector f α (x + v α ∆t, t + ∆t) for 0 ≤ α ≤ 8, if x + v α ∆t is a lattice node. For a boundary node as x b of Figure 3, we use (with the usual numbering of the degrees of freedom for D2Q9 scheme [START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability[END_REF]) the following bounce-back boundary condition (31)

   f 2 (x b , t + ∆t) = f 4 (x e , t + ∆t) = f * 4 (x b , t) f 5 (x b , t + ∆t) = f 7 (x c , t + ∆t) = f * 7 (x b , t) f 6 (x b , t + ∆t) = f 8 (x d , t + ∆t) = f * 8 (x b , t) .
Periodic boundary conditions are considered in the longitudinal direction for abscissae equal to 0 and L. We repeat those steps until convergence to a steady state.

 Proposition 3.

Order three for bounce-back For the D2Q9 lattice Boltzmann scheme (1) (26) (27 ) the bounce-back numerical boundary condition (31) is of order 3 at location ∆q = ∆x 2 for the Dirichlet boundary condition u = 0 if and only if σ 5 σ 8 = 3 8 . Proof of Proposition 3. We calculate the defects of conservation θ k defined by [START_REF] Ginzburg | Two-relaxation-time lattice Boltzmann scheme: About parametrization, velocity, pressure and mixed boundary conditions[END_REF] for k > 3:

θ 3 = ∂ t (αρ) = -α divj + O(∆x 2 ) , θ 4 = ∂ t (βρ) -divj = -(β + 1) divj + O(∆x 2 ) , θ 5 = - ∂ t j x λ + λ(α + β) 3 ∂ x ρ = λ 6 (4 + 3α + 2β) ∂ x ρ + O(∆x) , θ 6 = - ∂ t j y λ + λ(α + β) 3 ∂ y ρ = λ 6 (4 + 3α + 2β) ∂ y ρ + O(∆x) , θ 7 = 2 3 (∂ x j x -∂ y j y ) θ 8 = 1 3 (∂ y j x + ∂ x j y ) .
Nonequilibrium moments at second order are given by the expansion ( 14) (justified in [START_REF] Dubois | Third order equivalent equation of lattice Boltzmann scheme[END_REF][START_REF] Dubois | Using the Lattice Boltzmann Scheme for Anisotropic Diffusion Problems, Finite Volumes for complex Applications V[END_REF]). Then we have:

m * 3 = αρ + ∆t 1 2 -σ 3 -α divj + ∆t λ 2 6 (σ 3 α(4 + α) + σ 5 (4 + 3α + 2β))∆ρ + O(∆x 3 ), m * 4 = βρ+∆t ( 1 2 -σ 4 ) -(β+1)divj+∆t λ 2 6
(σ 4 (β+1)(4+α)+σ 5 (4+3α+2β))∆ρ +O(∆x 3 ),

m * 5 = - j x λ + ∆t 1 2 -σ 5 λ (4 + 3α + 2β) 6 ∂ x ρ + ∆t λ 3 (4 + 3α + 2β) 2 σ 5 ∂ x divj +ασ 3 ∂ x divj + (β + 1)σ 4 ∂ x divj + 2σ 8 ∂ x (∂ x j x -∂ y j y ) -σ 8 ∂ y (∂ y j x + ∂ x j y ) + O(∆x 3 ), m * 6 = - j y λ + ∆t 1 2 -σ 5 λ (4 + 3α + 2β) 6 ∂ y ρ + ∆t λ 3 (4 + 3α + 2β) 2 σ 5 ∂ y divj +ασ 3 ∂ y divj + (β + 1)σ 4 ∂ y divj -2σ 8 ∂ y (∂ x j x -∂ y j y ) -σ 8 ∂ x (∂ y j x + ∂ x j y ) + O(∆x 3 ), m * 7 = ∆t 1 2 -σ 8 2 3 (∂ x j x -∂ y j y ) + ∆t λ 2 9 σ 8 (4 + α) + σ 5 2 (4 + 3α + 2β) (∂ 2 x ρ -∂ 2 y ρ) + O(∆x 3 ), m * 8 = ∆t 1 2 -σ 8 1 3 (∂ y j x + ∂ x j y ) + ∆t λ 2 9 (σ 8 (4 + α) -σ 5 (4 + 3α + 2β))∂ xy ρ + O(∆x 3 ) .
We have j y = 0, j x = j x0 + y∂ y j x + y 2 2 ∂ 2 y j x and ρ = constant. We evaluate the non conserved moments m * k (x b ) and add m

1 (x b ) = j x (x b ) = j x (x i ) + ∆x 2 3 σ 8 ∂ 2 y j x (x i ) + O(∆x 3 ). We compute moments m * k (x e ), m * k (x c
) and m * k (x d ) at the "external nodes" depicted in Figure 3. Using the matrix M -1 we evaluate f

* k (x b ), f * k (x c ), f * k (x d ) and f * k (x e ). Finally we obtain (32) f * 5 (x c ) -f * 7 (x b ) = 1 6 j x (x i ) + ∆x 2 144 (8σ 5 σ 8 -3) ∂ 2 j x ∂y 2 (x i ) + O(∆x 3 )
verse profile for longitudinal momentum. the algebra then follows what is presented for proposition 3.

• We then perform simulations of the two situations discussed above. For this we consider a domain of size N x = 100, N y = 21 and analyze the flow in the steady state.

For several values of σ 5 and of σ 8 , we determine a parabola by best fit with the velocity profile in the middle section of the channel. We verified that the domain was long enough in order to reduce to a negligible level the errors due to mismatch in the end boundary conditions for links that intersect both a solid boundary (imposed flux) and the input boundary (imposed pressure), that would require a more sophisticated treatment.

As in relation [START_REF] Ginzburg | Second order boundary method for Lattice Boltzmann model[END_REF], we define ∆q as the experimental point where the parabola goes through zero. The results (Figure 4) depend only upon the product σ 5 σ 8 and are coherent with the theoretical results established in propositions 3 and 4. For α = -2 and β = 1, the superconvergent accuracy is obtained "experimentally exactly" at the boundary for

Conclusion

The "magic" parameters introduced by Ginzburg and Adler [START_REF] Ginzbourg | Boundary flow condition analysis for three-dimensional lattice Boltzmann model[END_REF] which allow to increase the accuracy of lattice Boltzmann simulations in the presence of solid boundaries have been considered for a few simple situations. We have shown that they depend upon the choice of moments and of their equilibrium values. In addition they depend upon the way the flow is driven. The analysis requires the determination of the non equilibrium moments up to second order accuracy and this explicitation is obtained in the framework of the Taylor expansion method. Note that the same results could be obtained with the Chapman-Enskog procedure. The work described here can easily be extended to more complicated lattice Boltzmann schemes for boundaries parallel to one of the velocities of the model. In all cases that we considered, the results can be expressed in terms of products of the type σ i σ j , where σ i corresponds to the relevant transport coefficient (diffusivity or shear viscosity) and σ j to other moments of opposite symmetry (i.e. odd order moments of f , "energy flux" and higher order terms of the same symmetry for models with a large enough number of velocities), and thus the "magic" conditions are the same as those presented in the comprehensive paper of Ginzburg, Verhaeghe and d'Humières [START_REF] Ginzburg | Two-relaxation-time lattice Boltzmann scheme: About parametrization, velocity, pressure and mixed boundary conditions[END_REF]. They are also valid for special BGK situations that we have in addition to the "magic" conditions, σ i = σ j .

FrançoisFigure 1 :

 1 Figure 1: A boundary surface cutting the link between a fluid node x b and a fictitious outside node x e ≡ x b -∆x.

Figure 2 :

 2 Figure 2: "Experimental numerical location" of the solid wall ∆q versus σ 1 σ 2 . The first D1Q3 model with symbol: superconvergent parameters such that σ 1 σ 2 = 1 8 . The second D1Q3 model with symbol: superconvergent parameters satisfying σ 1 σ 2 = 3 8 .

Figure 3 :

 3 Figure 3: Domain Ω = ]0, L[ × ]0, H[ (left) and notations for the numerical treatment of a boundary vertex x b at the bottom of domain Ω (right).

  σ 5 σ 8 = 3 16 . When α = -β = -52 , the same observation occurs for σ 5 σ 8 = 3 8 .

Figure 4 : 2 for σ 5 σ 8

 428 Figure 4: Product σ 5 σ 8 versus solid wall location ∆q with N y = 21 for the second lattice Boltzmann boundary scheme (33). Parameters α = -2, β = 1 with symbol, parameters α = -2.5, β = 2.5 with symbol. The boundary is experimentally located at ∆q = ∆x 2 for σ 5 σ 8 = 3 16 in the first case and σ 5 σ 8 = 3 8 in the second, as suggested in Proposition 4.

Acknowledgments

The referees conveyed to the authors very interesting remarks that have been incorporated into the present edition of the article.

and similar relations for f * 2 (x e ) -f * 4 (x b ) and f * 6 (x d ) -f * 8 (x b ). The conclusion is clear: when the left hand side of (32) is identically null due to the boundary condition (31), the momentum j x (x i ) on the surface located at ∆q = 1 2 ∆x is null "up to third order accuracy" as defined in [START_REF] Ginzburg | Second order boundary method for Lattice Boltzmann model[END_REF] if and only if the relation 8σ 5 σ 8 -3 = 0 occurs.

• We remark that if we apply the body force following the algorithm (i) m = M f, (ii) collision, (iii) f = M -1 m, (iv) apply the body force following the precise relations for transformation of particle distribution f -→ f :

λ and m k = m k for the other moments, the solid wall for the Poiseuille problem is "numerically located" at ∆q = ∆x 2 up to third order accuracy if the relation σ 5 σ 8 = 3 16 is satisfied between the relaxation parameters, as proposed by Ginzburg and d'Humières [START_REF] Ginzburg | Generic boundary conditions for lattice Boltzmann models and their application to advection and anisotropic dispersion equations[END_REF][START_REF] Ginzburg | Multireflection boundary conditions for lattice Boltzmann models[END_REF][START_REF] Ginzburg | Two-relaxation-time lattice Boltzmann scheme: About parametrization, velocity, pressure and mixed boundary conditions[END_REF].

• A second lattice Boltzmann scheme We can also model the Poiseuille flow described by (28) (29) with the explicit introduction of a pressure gradient δp. So the scheme (26) (27) has the same steps as the preceeding scheme with F x ≡ 0 and the wall boundary conditions are still given by (31). We consider the boundary condition for nodes x ≡ (k ∆x, ℓ ∆x) at the entrance (k = 1) and at the output (k = N x ) as follows:

(33)

with δρ = δp/c 2 s the density drop corresponding to the pressure step considered in (29), ( c s is the speed of sound) and (α, β) parameters for equilibrium introduced at the relation (27). Note that these expressions may be called "anti bounce-back" with an imposed scalar quantity (similar to what is used when the lattice Boltzmann scheme is set to simulate diffusion problems). In this case we perform the same proof as for proposition 3, we take F x = 0 and the exact solution is given by a linear longitudinal profile for density and a parabolic trans-