Kévin Guillon

Romane Hélie

Philippe Helluy

STABILITY ANALYSIS OF THE VECTORIAL LATTICE-BOLTZMANN METHOD

Keywords: vectorial lattice-Boltzmann, entropy stability, equivalent equation analysis

. It is a variant of the LBM with extended stability features: it allows handling compressible flows with shock waves, while the LBM is limited to low-Mach number regime. The stability analysis is based on the Legendre transform theory. We also propose a new tool: the equivalent system analysis, which we conjecture to contains both the stability and the consistency of the VLBM.

Introduction

The Vectorial Lattice-Boltzmann Method (VLBM) is a variation of the Lattice-Boltzmann Method (LBM) proposed by several authors [START_REF] Bouchut | Construction of BGK models with a family of kinetic entropies for a given system of conservation laws[END_REF][START_REF] Aregba | Discrete kinetic schemes for multidimensional systems of conservation laws[END_REF][START_REF] Paul | Lattice kinetic schemes for magnetohydrodynamics[END_REF][START_REF] Dubois | Simulation of strong nonlinear waves with vectorial lattice Boltzmann schemes[END_REF][START_REF] Coulette | High-order implicit palindromic discontinuous Galerkin method for kinetic-relaxation approximation[END_REF][START_REF] Baty | A robust and efficient solver based on kinetic schemes for magnetohydrodynamics (MHD) equations[END_REF] for solving systems of conservation laws. Compared to the traditional LBM, which uses a scalar distribution function, the VLBM utilizes a vectorial distribution function and offers several advantages. For instance, it allows for a fully rigorous entropy stability analysis [START_REF] Bouchut | Construction of BGK models with a family of kinetic entropies for a given system of conservation laws[END_REF][START_REF] Dubois | Simulation of strong nonlinear waves with vectorial lattice Boltzmann schemes[END_REF], which ensures that the resulting schemes can compute shock waves [START_REF] Baty | A robust and efficient solver based on kinetic schemes for magnetohydrodynamics (MHD) equations[END_REF] while the standard LBM is only stable on low-Mach number flows.

Both the LBM and the VLBM are based on a kinetic representation of the system of conservation laws, which consists of a set of transport equations coupled through a stiff relaxation term with a small parameter ε. Solving the transport and relaxation separately through a splitting method results in simple and efficient schemes.

However, this methodology raises several natural questions, three of which are listed below:

(1) Is the kinetic representation stable as ε → 0?

(2) Is the kinetic representation consistent with the system of conservation laws as ε → 0?

(3) Is the splitting scheme consistent with the kinetic representation and/or the system of conservation laws as ε → 0? The first question can be answered by an entropy analysis, which has been explored in prior work [START_REF] Bouchut | Construction of BGK models with a family of kinetic entropies for a given system of conservation laws[END_REF][START_REF] Dubois | Simulation of strong nonlinear waves with vectorial lattice Boltzmann schemes[END_REF]. However, the fully rigorous theory has not yet been achieved, especially for general systems of conservation laws in high-dimensional spaces. This is related to challenging questions, such as the mathematical theory of Navier-Stokes equations [START_REF] Saint-Raymond | Hydrodynamic limits of the Boltzmann equation[END_REF], strange behavior of the general solutions of conservation laws [START_REF] De | The H-principle and the equations of fluid dynamics[END_REF], and uncertainty quantification [START_REF] Ulrik S Fjordholm | Construction of approximate entropy measure-valued solutions for hyperbolic systems of conservation laws[END_REF].

The second question regarding whether or not the kinetic model converges towards the approximated conservation laws can be fully answered in some specific cases. However, the general case is also difficult. There are some semi-heuristic arguments, based on the Chapman-Enskog expansion [START_REF] Chapman | The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases[END_REF], that provide explanations for the convergence of the stiff relaxation approximation. Among many works in this direction, [START_REF] Paul | Bulk and shear viscosities in lattice Boltzmann equations[END_REF] presents one approach.

The third question is related to the fact that the stiff relaxation is never resolved in practice. Instead, a time splitting scheme with over-relaxation is used. In most situations, this splitting algorithm does not provide a good approximation of the relaxation model, but rather it approximates the initial conservation laws directly, which is the ultimate objective. Theoretically, this approach is justified with an equivalent equation analysis, performed directly on the time splitting scheme without the relaxation intermediary. The analysis relies on a combination of Taylor and Chapman-Enskog expansions [START_REF] Paul | An interpretation and derivation of the lattice Boltzmann method using Strang splitting[END_REF] and generally provides relevant information on the consistency of the LBM or VLBM method.

However, the stability of the equivalent equation is generally not sufficient to ensure the stability of the VLBM [START_REF] Bouchut | Stability of relaxation models for conservation laws[END_REF]. Thus, separate studies are necessary for both the stability and the consistency of the method in practice.

In this work, we aim to provide a more comprehensive justification for the stability and consistency of the VLBM method. We begin by reviewing the entropy stability of the VLBM, which can be proven using the Legendre transform theory. We demonstrate that this approach can be naturally extended to prove the stability of the VLBM when an over-relaxed splitting scheme is used for time integration.

To analyze consistency, we propose a new algorithmic approach for constructing a system of equivalent Partial Differential Equations (PDEs) that includes not only the conservative variables but also the equilibrium deviation variables. In several examples, we show that the analysis of the equivalent system and the entropy analysis yield the same stability condition.

We then utilize an additional Chapman-Enskog expansion to remove the equilibrium deviation variables, under a smallness hypothesis. This enables us to recover the traditional equivalent equation provided by other authors [START_REF] Dubois | Equivalent partial differential equations of a lattice Boltzmann scheme[END_REF][START_REF] Graille | Approximation of mono-dimensional hyperbolic systems: A lattice Boltzmann scheme as a relaxation method[END_REF][START_REF] Courtès | Vectorial kinetic relaxation model with central velocity. application to implicit relaxations schemes[END_REF].

Finally, we present numerical experiments to verify that the equivalent system with deviation variables provides more accurate information about the stability of the VLBM than the standard equivalent equation.

Hyperbolic conservation law and duality

2.1. Hyperbolic systems of conservation laws. In this work, we consider the numerical resolution of a system of conservation laws (2.1)

∂ t W + d i=1 ∂ i Q i (W ) = 0,
where the unknown is a vector of m conserved quantities W (X, t) ∈ R m , depending on a space variable X ∈ R d and a time variable t ∈ R. Let N = (N 1 , . . . , N d ) ∈ R d be a space vector. The flux of the system is defined by

Q(W, N ) = d i=1 Q i (W )N i .
We assume that the system of conservation laws admits a Lax entropy W → s(W ). Therefore, there is an entropy flux g i (W )N i such that

∂ t s(W ) + d i=1 ∂ i g i (W ) = 0,
whenever W is a smooth solution of (2.1). This imposes that (2.2)

D W s(W )D W Q i (W ) = D W g i (W ),
where we have denoted by D W g(W ) the Jacobian of g(W ). Let us recall that the Jacobian is the transpose of the gradient

(2.3) D W s(W ) = ∇ W s(W ) ⊺ .
Thus D W s(W ) is a row vector, while ∇ W s(W ) is a column vector.

In addition, we assume a convexity hypothesis of s on a closed convex cone K ⊂ R m :

s is strictly convex on K.

Finally, following a standard convex analysis convention, we extend s by an infinite value outside K:

(2.4)

s(W ) = +∞, W ∈ ∁K.
Remark 1. For working in a fully practical framework, it is indeed necessary to consider cases where K ̸ = R d . However, this leads to mathematical questions that we have not yet investigated. Therefore, in the following, the reader can suppose that K = R d .

Entropy and symmetrization.

According to the Mock theorem [START_REF] Peter | Hyperbolic systems of conservation laws and the mathematical theory of shock waves[END_REF][START_REF] Michael | Systems of conservation laws of mixed type[END_REF], the system (2.1) is symmetrizable and thus hyperbolic: for all W ∈ K and all N ∈ R d the Jacobian of the flux

A(W, N ) = D W Q(W, N ) is diagonalizable with real eigenvalues.
Let us recall how to prove this result, because it will be useful. For this we introduce the conjugate of the entropy [START_REF] Hiriart-Urruty | Fundamentals of convex analysis[END_REF][START_REF] Hiriart-Urruty | Convex analysis and minimization algorithms I: Fundamentals[END_REF] defined by

(2.5) s * (W * ) = max W (W * • W -s(W )),
where we denote by • the usual dot product:

W * • W = W * ⊺ W = W ⊺ W * .
Thanks to (2.4) we can also write

(2.6) s * (W * ) = max W ∈K (W * • W -s(W )),
The components of W * are called the dual variables, or entropy variables [START_REF] Harten | On the symmetric form of systems of conservation laws with entropy[END_REF][START_REF] Bourdel | On the approximation of k-diagonalizable hyperbolic systems by finite elements-applications to the Euler equations and to gaseous mixtures[END_REF][START_REF] Croisille | Contribution à l'étude théorique et à l'approximation par éléments finis du système hyperbolique de la dynamique des gaz multidimensionnelle et multiespèces[END_REF]]. The function s * is called the conjugate or the dual entropy. The definition of the conjugate (2.5) applies to functions that are not necessarily smooth or convex. In the regular case, when s is smooth and strictly convex on R m , for instance, s * is defined implicitly by the following two relations (2.7)

W * = ∇s(W (W * )), (2.8) s * (W * ) = W * • W (W * ) -s(W (W * )).
The formula (2.7) determines uniquely the map between the conservative variables W and the dual variables W * . In the convex case, it can also be shown that s * * = s. Thus, we have the reverse relations:

W = ∇s * (W * (W )), s(W ) = W * (W ) • W -s * (W * (W )).
We can also define the dual entropy flux by the relation (2.9)

P i,⋆ (W * ) = W * • Q i (W (W * )) -g i (W (W * )).
We do not use exactly the same symbol for the dual entropy ( * ) and the dual flux (⋆), because the definitions are slightly different. An important fact is that the knowledge of s * (W * ) and P i,⋆ (W * ) is sufficient to reconstruct the system of conservation laws (2.1). Indeed,

∇ W * P i,⋆ (W * ) = Q i (W (W * )) + W * • D W Q i (W (W * ))D W * W (W * ) -D W g i D W * W (W * ), = Q i (W (W * )) + D W s(W )D W Q i (W (W * ))D W * W (W * ) -D W g i D W * W (W * ),
(because of (2.3) and (2.2))

= Q i (W (W * )) + D W g i (W (W * ))D W * W (W * ) -D W g i (W (W * ))D W * W (W * ),
(because of (2.7)), and thus (2.10)

∇ W * P i,⋆ (W * ) = Q i (W (W * )).
In short: the gradient of the dual entropy gives the conservative variables and the gradient of the dual flux gives the flux.

Theorem 2. The change of variables W * → W (W * ) symmetrizes the system of conservation laws (2.1).

Proof.

Indeed ∂ t W + i ∂ i Q i (W ) = 0, can be rewritten ∂ t ∇ W * s * (W * ) + i ∂ i ∇ W * P i,⋆ (W * ) = 0,
thanks to the above remark (2.7) and (2.10), or

(2.11) D 2 W * W * s * (W * )∂ t W * + i D 2 W * W * P i,⋆ (W * )∂ i W * = 0,
where the Hessian matrices D We have just proven the Mock theorem [START_REF] Michael | Systems of conservation laws of mixed type[END_REF]. In the following sections, we shall often try to guess directly a symmetrization of the system of conservation laws, rather than the full dual entropy theory to obtain the stability conditions of the relaxation scheme.

Remark 3. In the following, the practical stability will derive from a convexity condition imposed on s * . This may sound strange, because the conjugate of any function is always convex, according to the standard definition (2.5). We need to be a little bit more subtle. Actually, it is possible to define two different duality transformations. Let us denote by Fenchel transform the s * given by the first formula (2.5). The Fenchel transform derives from an optimization problem. Let us denote by Legendre transformation the s * given by formula (2.7) and (2.8). The Legendre transform is thus defined by purely algebraic relations. It is unambiguously defined, as soon as W → W * = ∇ W s(W ) is an invertible map. In this algebraic definition, s * is not necessarily convex, but the convexity of s * is equivalent to the convexity of s. The Legendre and Fenchel transform are different in the general case. They coincide in the convex case.

Remark 4. In the most general case, the Legendre transform is a multivalued function. It has to be defined from differential geometry tools. The interested reader can refer to [START_REF] Ekeland | Legendre duality in nonconvex optimization and calculus of variations[END_REF] for an introduction to this topic.

Vectorial kinetic model

3.1. Direct construction. In this section, we recall the entropy theory of the kinetic representation. This theory has a long history, see for instance [START_REF] Lax | Shock waves and entropy[END_REF][START_REF] Michael | Systems of conservation laws of mixed type[END_REF][START_REF] Harten | On the symmetric form of systems of conservation laws with entropy[END_REF][START_REF] Suresh | A second-order accurate kinetic-theory-based method for inviscid compressible flows[END_REF][START_REF] Bourdel | On the approximation of k-diagonalizable hyperbolic systems by finite elements-applications to the Euler equations and to gaseous mixtures[END_REF][START_REF] Croisille | Contribution à l'étude théorique et à l'approximation par éléments finis du système hyperbolique de la dynamique des gaz multidimensionnelle et multiespèces[END_REF][START_REF] Perthame | Boltzmann type schemes for gas dynamics and the entropy property[END_REF]. In our context it has been analyzed by Bouchut in [START_REF] Bouchut | Construction of BGK models with a family of kinetic entropies for a given system of conservation laws[END_REF][START_REF] Bouchut | Entropy satisfying flux vector splittings and kinetic BGK models[END_REF]. The ideas have been rephrased with a systematic use of the Legendre transform in the work of Dubois [START_REF] Dubois | Simulation of strong nonlinear waves with vectorial lattice Boltzmann schemes[END_REF]. Let us now recall the theory.

We consider a formal vectorial kinetic representation of system (2.1)

(3.1) ∂ t F k + V k • ∇F k = 1 ε (F eq k -F k ), k = 1 . . . n v .
The approximate conservative vector is the sum of the kinetic vectors

(3.2) W = nv k=1 F k ,
and the kinetic equilibrium vectors (or Maxwellians) are functions of the conservative data

(3.3) F eq k = F eq k (W
). The kinetic velocities V k are n v constant and given vectors of R d . In practice, it is interesting to introduce a positive parameter λ, whose purpose is to change the size of the kinetic velocities. So we shall often take

(3.4) V k = λ Ṽk ,
where the directions Ṽk of the kinetic velocities are fixed. In this way, the kinetic velocities can be dilated.

For practical reasons we also introduce the kinetic vectors F and F eq , which are column vectors made of all the stacked kinetic data:

F = (F ⊺ 1 , . . . , F ⊺ nv ) ⊺ , F eq = (F eq 1 ) ⊺ , . . . , F eq nv ⊺ ⊺ .
With the help of the following diagonal matrices (1 m is the identity matrix of size m × m)

V i =    V i 1 1 m . . . V i nv 1 m    ,
the kinetic system can also be written in the compact form

∂ t F + d i=1 ∂ i V i F = 1 ε (F eq -F ).
When ε → 0, we expect that

(3.5) F k ≃ F eq k . If we assume that (3.6) W = nv k=1 F eq k (W ), Q i (W ) = nv k=1 V i k F eq k (W ),
then, summing (3.1) on k and using (3.5) we formally obtain

∂ t W + d i=1 ∂ i Q i (W ) ≃ 0,
and we have obtained an approximation of the initial system of conservation laws.

Remark 5. For being more rigorous, we should have made more explicit in (3.1) and (3.2) the dependency of F k and W = k F k with respect to ε. For instance, by using the notations

F (ε) k
and W (ε) . But we prefer to lighten the notations, and in the following the dependency with ε will be implicit. This means that from now on, W = k F k is not an exact solution of (2.1) but an approximate one.

For the moment, instead of relating the equilibrium (3.3) to the conservation laws (2.1), we assume that the equilibrium is obtained from an entropy optimization principle. For this we introduce a microscopic entropy

σ(F ) = nv k=1 s k (F k ),
where the kinetic entropies s k are strictly convex functions of the F k on K. The macroscopic entropy is obtained from the resolution of the following constrained optimization problem

(3.7) s(W ) = min W = k F k σ(F ).
In optimization theory, this operation is known as an inf-convolution operation [START_REF] Hiriart-Urruty | Convex analysis and minimization algorithms I: Fundamentals[END_REF]. The macroscopic entropy is the inf-convolution of the kinetic entropies. In many works, the inf-convolution operator is denoted with a □. We thus have:

s = s 1 □s 2 . . . □s nv .
We denote by F eq k (W ) the (supposed to be unique) values of F k that achieve the minimum

s(W ) = nv k=1 s k (F eq k (W )).
If we introduce the Lagrangian

L(F, Λ) = nv k=1 s k (F k ) + Λ • W - nv k=1 F k .
The minimizer F eq (W ) and the Lagrange multiplier Λ(W ) are characterized by

(3.8) ∇ F k s k (F eq k ) = Λ, nv k=1 F eq k = W.
These relations are simply obtained by deriving the Lagrangian with respect to F k or Λ.

An essential property of the inf-convolution is that the Fenchel transform changes it into a sum. We thus have

(3.9) s * (W * ) = nv k=1 s * k (W * ).
Taking the Legendre transform of (3.8) we see that the couple (F eq k (W ), Λ(W )) is a solution to (3.10)

F eq k = ∇ Λ s * k (Λ(W )), nv k=1 F eq k = W.
Summing over k we also have the Lagrange multiplier in an easier way nv k=1

F eq k = nv k=1 ∇ Λ s * k (Λ(W )), W = ∇ Λ nv k=1 s * k (Λ(W )), = ∇ Λ s * (Λ(W )
), and thus Λ(W ) = W * (W ) = ∇s(W ). The Lagrange multiplier of the constrained optimization problem (3.8) is simply the gradient of the macroscopic entropy.

Let us now assume the additional property

(3.11) P i,⋆ (W * ) = nv k=1 V i k s * k (W * ).
Then we have

∇ W * P i,⋆ (W * ) = nv k=1 V i k ∇ W * s * k (W * ), ∇ W * P i,⋆ (W * ) = nv k=1 V i k F eq k (W ), (from (3
.10)) and thus, from (2.10)

(3.12) Q i (W ) = nv k=1 V i k F eq k (W ).
In the flux form, we can also write

Q(W, N ) = nv k=1 (V k • N )F eq k (W ).
We recover the relations (3.6) that impose the consistency of the kinetic model (3.1) with the system of conservation laws (2.1). But now the consistency derives from an entropy optimization principle. We can sum up the direct construction of a kinetic model (3.1) that is entropy consistent with (2.1):

(1) Compute the Legendre transform s * of the entropy s by (2.7), (2.8).

(2) Compute the dual fluxes P i,⋆ by (2.9).

(3) Find n v strictly convex functions s * k satisfying the consistency relations written in the dual variables:

(3.13) nv k=1 s * k = s * , nv k=1 V i k s * k = P i,⋆ . ( 4
)
The equilibrium is given by F eq k = ∇s * k . Remark 6. If n v ⩾ d + 1, there is generally at least one solution to the algebraic system (3.13). The difficulty is to ensure that the s * k are convex. For this property to hold, the scaling (3.4) is useful. It gives an additional degree of freedom for ensuring convexity.

Remark 7. The main advantage of constructing the kinetic entropy by the dual analysis is theoretical. In the end the algebraic consistency relations given in (3.6) are recovered. The dual analysis ensures that the equilibrium kinetic vector F eq k corresponds to a minimum of the kinetic entropy.

3.2. Reverse construction. Now that we have recalled the entropy theory of the kinetic representation, we can proceed in the reverse way. We choose the equilibrium F eq in such a way that the consistency relation (3.6) is satisfied. From the above theory, we expect that F eq k (W (W * )) is a gradient, when it is expressed in the entropy variables W * . We can thus find dual kinetic entropies

s * k (W * ) such that F eq k (W (W * )) = ∇ W * s * k (W * ).
By Legendre transform, we can (in principle) compute the kinetic entropies s k (F k ) and this gives us the microscopic entropy

σ(F ) = nv k=1 s k (F k ).
The main point in the reverse construction is to ensure that the strict convexity is preserved. In practice, we will see that the microscopic entropy is convex under a condition that λ is large enough. This will be the subcharacteristic condition.

Splitting and over-relaxation

4.1. Splitting scheme. In this section we recall how to solve practically in time the kinetic system (3.1). Indeed, in (3.1) all the transport equations are coupled in a non-linear way. We introduce a splitting method for separating the transport equations. In this numerical method, the approximation of the kinetic data is computed at fixed times t n = n∆t, where ∆t is the time step. The approximation is not continuous at time t n , i.e. we distinguish between the values of

F k (X, t - n ) and of F k (X, t + n ).
Suppose that we know the kinetic data F k (X, t + n-1 ) at the end of time step n -1. Computing the next time step consists first in solving the homogenous transport equations

∂ t F k + V k • ∇F k = 0, for a duration of ∆t, which defines F k (X, t - n ).
Indeed, using the characteristic method (and avoiding for the moment difficulties arising from the boundaries), we get the explicit formula

(4.1) F k (X, t - n ) = F k (X -∆tV k , t + n-1
). In this way, we obtain the new conservative data by

W (X, t n ) = nv k=1 F k (X, t - n ).
We then apply the following over-relaxation formula

(4.2) F k (X, t + n ) = ωF eq k (W (X, t n )) + (1 -ω)F k (X, t - n ).
Several choices can be made for the relaxation parameter ω. The choice ω = 1 corresponds to a projection of the kinetic data on the equilibrium at the end of each time step. The choice ω = 2 (over-relaxation) is interesting, because, in this case, it can be shown that the time integration is second-order accurate (see for instance [START_REF] Coulette | High-order implicit palindromic discontinuous Galerkin method for kinetic-relaxation approximation[END_REF] and included references).

4.2.

Operator notations. The splitting scheme is very simple to implement in a computer program and no additional information is needed. However, we introduce here some additional mathematical notations that will be useful for deriving the equivalent equation analysis.

The interesting output of the kinetic scheme is obviously the conservative variables vector of size m, given by

W = nv k=1 F k .
However, the kinetic representation introduces n v kinetic vectors F k instead of one vector W . For the analysis, it is necessary to supplement W with additional quantities. The equivalent equation could be written in an arbitrary set of variables Y ∈ R mnv . We have done a particular choice, which helps the understanding, to our opinion.

The first components of Y are made of the conservative variables W . The next components will be made of the flux errors

nv k=1 V i k (F k -F eq k ) = nv k=1 V i k F k -Q i , 1 ⩽ i ⩽ d. If n v = d + 1, we have enough variables. If n v > d + 1 we still have to supplement Y with n c = (n v -d -1) independent linear combinations of the kinetic data nv k=1 β ℓ k (F k -F eq k ) , 1 ⩽ ℓ ⩽ n v -d -1 = n c .
It is convenient to choose the coefficients β ℓ k of the linear combinations in order to cancel the contributions of F eq (4.3)

nv k=1 β ℓ k F eq k = 0.
Finally, Y is of the form

Y =             W k V 1 k F k . . . k V d k F k k β 1 k F k . . . k β nv-d-1 k F k             -             0 k V 1 k F eq k . . . k V d k F eq k 0 . . . 0             .
While the m first components of Y are an approximate solution to (2.1)

Y 1•••m = W,
we expect the next components to be small

Y m+1•••mnv ≃ 0.
In the following, we denote the sub-vector Y m+1•••mnv as the flux error (even if the last components of Y , as we have seen above, are not necessarily all related to fluxes).

We can write the Y representation in a matrix form

(4.4) Y (F ) = M F -CM F eq (W (F )), W (F ) = nv k=1 F k .
We call the matrix M the matrix of moments. It is supposed to be invertible. The matrix C is the diagonal matrix

C =   0 m 1 md 0 mnc   ,
where we denote by 1 r the diagonal identity matrix of size r × r and by 0 r the null matrix of size r × r. The nonlinear mapping F → Y can be inverted explicitly

(4.5) F (Y ) = M -1 Y + M -1 CM F eq (W (Y )), W (Y ) = (Y 1 , . . . , Y m ) ⊺ .
The notations are tedious but simple. For making them clearer, let us consider three examples that we will use below.

4.2.1. D1Q2 case. In this case, the space dimension d = 1 and we take two kinetic velocities

V 1 = -λ, V 2 = λ.
Because n v = d + 1 and from the consistency relation (3.12), the equilibrium is necessarily given by

F eq 1 (W ) = 1 2 W - 1 2λ Q(W ), F eq 2 (W ) = 1 2 W + 1 2λ Q(W ). Then M = 1 m 1 m -λ1 m λ1 m .
4.2.2. D2Q3 case. In this case, the space dimension d = 2 and we take n v = 3 kinetic velocities

V 1 = λ 1 0 , V 2 = λ 2 -1 √ 3 , V 3 = λ 2 -1 - √ 3 .
Because of n v = d + 1 the equilibrium is uniquely defined. It is given by

F eq = M -1   W Q 1 Q 2   , with M =   1 m 1 m 1 m λ1 m -1 2 λ1 m -λ 2 1 m 0 m √ 3 2 λ1 m - √ 3 2 λ1 m   .
We find

F eq k = W 3 + 2 3λ 2 V k • Q 1 Q 2 .
4.2.3. D2Q4 case. In this case, the space dimension d = 2 and we take four kinetic velocities

V 1 = λ 1 0 , V 2 = λ -1 0 , V 3 = λ 0 1 , V 4 = λ 0 -1 .
Because n v > d + 1, we have one degree of freedom for the computation of the equilibrium. We can take, by analogy with the D1Q2 model:

F eq 1,2 (W ) = 1 4 W ± 1 2λ Q 1 (W ), F eq 3,4 (W ) = 1 4 W ± 1 2λ Q 2 (W ).
Then the matrix of moments is given by

M =     1 m 1 m 1 m 1 m λ1 m -λ1 m 0 m 0 m 0 m 0 m λ1 m -λ1 m λ 2 1 m λ 2 1 m -λ 2 1 m -λ 2 1 m     .
For the D2Q4 model, other choices are possible for the last row of the matrix of moments. We have made the choice proposed in [START_REF] Février | Extension et analyse des schémas de Boltzmann sur réseau: les schémas à vitesse relative[END_REF], but we could have replaced λ 2 by one, for instance. Actually, we could have supplemented by any row that ensures that M is invertible, but then the condition (4.3) is generally not satisfied.

4.2.4.

Transport operator in the Y variables. With these notations, we can rewrite the transport operator in the Y variables. Let us define the τ k (∆t) shift operators, applied on an arbitrary function X → f (X) by

(τ k (∆t)f )(X) = f (X -∆tV k ).
The vectorial shift operator is then given by

τ (∆t) =    τ 1 (∆t)1 m 0 . . . 0 τ nv (∆t)1 m    .
The initial data Y (•, t + n-1 ) at the end of the time step n-1 are supposed to be known. The transport map T (∆t) is the procedure that takes Y (•, t + n-1 ) and produces Y (•, t - n ) before the relaxation step:

Y (•, t - n ) = T (∆t)Y (•, t + n-1
). The procedure is as follows:

(1) Express the initial data Y in the kinetic variables thanks to (4.5)

F (•, t + n-1 ) = F (Y (•, t + n-1 )). (2) Apply the transport operator F (•, t - n ) = τ (∆t)F (•, t + n-1 ). (3) Go back to the Y variables thanks to (4.4) Y (•, t - n ) = Y (F (•, t - n )).
The transport is nonlinear functional operator, made of shift operators and nonlinear algebraic operations. We do not give all the details. 4.2.5. Relaxation operator in the Y variables. The relaxation operator in the kinetic variables is given by (4.2). In vectorial form it reads

(4.6) F (•, t + n ) = ωF eq (W (•, t - n )) + (1 -ω)F (•, t - n ).
Using the fact that W = k F k does not change during the relaxation, we simply obtain

(4.7) Y (•, t + n ) = R ω Y (•, t - n )
, where the relaxation operator R ω is the diagonal operator

R ω = 1 m 0 0 (1 -ω)1 m(nv-1)
.

In other words, the relaxation operator leaves the conservative variables unchanged:

(R ω Y ) i = Y i = W i , 1 ⩽ i ⩽ m,
and the other components of Y are multiplied by 1 -ω:

(R ω Y ) i = (1 -ω)Y i , i > m.
When ω ≃ 2, 1 -ω ≃ -1 and thus, these components are numerically fluctuating around 0 at each time iteration. There is no damping of the oscillations when ω = 2. The frequency of the fluctuations is 1/∆t and tends to infinity when ∆t tends to zero.

4.2.6. Formal numerical scheme. Now that we have introduced the operator notation, we can give a more formal definition of the split kinetic approximation. For solving (2.1) with the initial data

W (X, 0) = W 0 (X),
we start with

Y + 0 = W 0 (X) 0 .
Assume that we have reached Y + n-1 at time step n -1. Then the next value is given by

Y + n = S(∆t)Y + n-1
, where the split kinetic operator is

(4.8) S(∆t) = R ω • T (∆t).
Then we extract the first m components of Y + n and we expect to obtain an approximation of W at time n:

Y + n (X) = W n (X) . . . , W (X, t n ) ≃ W n (X).
As stated above, the relaxation operator generates time fluctuations with frequency 1/∆t. For the analysis, it is convenient to suppress these non-relevant time fluctuations. This is done by considering only an even number of time steps of (4.8) for instance. Another point is that the split operator (4.8) is not symmetric in time. In order to remove the fast time oscillations and the lack of symmetry, we can consider the modified split kinetic operator (4.9)

S(∆t) = T ∆t 4 • R ω • T ∆t 2 • R ω • T ∆t 4 .
In the case ω = 2, this modified operator is indeed time symmetric in the sense that (4.10)

S(0) = I d , S(∆t) -1 = S(-∆t).
This is this property that ensures that the split scheme is second order accurate when ω = 2. We refer, for instance to [START_REF] Coulette | High-order implicit palindromic discontinuous Galerkin method for kinetic-relaxation approximation[END_REF] and to the consistency analysis given below. The symmetric operator and the non-symmetric operators sampled at the even time steps are similar. They differ in the initial and final steps by a term of order ∆t. This difference has in practice very few effect on the order of the scheme on the conservative variables (see Section 6.1.1 in [START_REF] Coulette | High-order implicit palindromic discontinuous Galerkin method for kinetic-relaxation approximation[END_REF]).

Applications

In this section, we apply the entropy theory to practical examples.

5.1. Application to the transport equation . We first consider the construction of the D1Q2 model for the one-dimensional transport equation

∂ t W + c∂ x W = 0,
where the velocity c is supposed to be constant and X = x is the one-dimensional space variable.

In this case

Q(W ) = cW, V 1 = -λ, V 2 =
λ, and we have no free choice for choosing the equilibrium kinetic data, which are given by

F eq 1 (W ) = W 2 - cW 2λ , F eq 2 (W ) = W 2 + cW 2λ .
For this simple linear conservation law we can take the entropy associated to the L 2 norm

s(W ) = W 2 2 .
The dual entropy is simply

s * (W * ) = W * 2 2 ,
and the entropy variable is

W * = ∇ W s(W ) = W. Thus F eq 1 (W (W * )) = W * 2 - cW * 2λ , F eq 2 (W (W * )) = W * 2 + cW * 2λ .
From (3.10) we deduce the dual kinetic entropies

s * 1 (W * ) = 1 4 (1 -c/λ)W * 2 , s * 2 (W * ) = 1 4 (1 + c/λ)W * 2 .
They are strictly convex under the subcharacteristic condition (5.1) λ > |c| .

We can then compute the kinetic entropies

s 1 (F 1 ) = λ λ -c F 2 1 , s 2 (F 2 ) = λ λ + c F 2 2 .
The microscopic entropy is then

σ(F 1 , F 2 ) = λ λ -c F 2 1 + λ λ + c F 2 2 .
As expected, it is a diagonal quadratic form in the (F 1 , F 2 ) variables.

Let us express the microscopic entropy with respect to the Y = (W, y) ⊺ variables. We have

W = F 1 + F 2 ,
and

y = -λF 1 + λF 2 -Q(W ), = -λF 1 + λF 2 -c(F 1 + F 2 ).
After simple computations, we find that the microscopic entropy is also

(5.2) σ(W, y) = σ(F 1 , F 2 ) = W 2 2 + y 2 2(λ 2 -c 2 )
.

It is a convex function of W and y under condition (5.1). As expected, it is minimal when the flux error y vanishes. In addition, in the relaxation step, the entropy is exactly conserved when ω = 2 because

(5.3) σ(W, (1 -ω)y) = σ(W, -y) = σ(W, y).
We are now in a position to prove the entropy stability of the over-relaxation scheme when

1 ⩽ ω ⩽ 2.
Theorem 8. With periodic boundary conditions, or in an infinite domain, the over-relaxation scheme is entropy stable under the sub-characteristic condition (5.1) when 1 ⩽ ω ⩽ 2.

Proof. It is sufficient to prove the decrease of the entropy

S(t) = x σ(F 1 (x, t), F 2 (x, t)) = x s 1 (F 1 ) + s 2 (F 2 ),
for a single time step. In the transport step, one solves

∂ t F k + V k ∂ x F k = 0,
and thus the microscopic entropies

s k (t) = x s k (F k ), are separately conserved s k (t + ∆t -) = s k (t +
). In the relaxation step, W is not changed and y(x, t + ∆t + ) = (1 -ω)y(x, t + ∆t -), because |1 -ω| ⩽ 1 we see from the expression (5.2) of the entropy in the (W, y) variable that the microscopic entropy decreases pointwise, at each x. Therefore, S(t + ∆t + ) ⩽ S(t + ∆t -). □ 5.2. Application to the shallow water equations. In order to show that the approach still works for non-linear systems of conservation laws, we try now to apply the above method to the shallow water model where the unknowns are the water height h(x, t) and the velocity u(x, t). It reads

∂ t W + ∂ x Q(W ) = 0, with W = h hu , Q(W ) = hu hu 2 + gh 2 /2 .
We define the primitive variables v = h u .

For smooth solutions, we also have

∂ t v + B(v)∂ x v = 0, with B(v) = u h g u .
Assume that the Lax entropy s(W ) = H(v) is expressed in the primitive variables, and that the entropy flux G(W ) = R(v). Then we must have

D v H(v)B(v) = D v R(v).
Denoting the partial derivatives with indices we obtain

H h H u u h g u = R h R u .
We search H under the form

H(h, u) = h u 2 2 + e(h).
Because

H h = u 2 2 + e ′ (h), H u = hu,
this gives u 3 2 + ue ′ + ghu = R h , 3hu 2 2 + he ′ = R u .
We take

R = h u 3 2 + ue + gu h 2 2 .
Then

R u = 3hu 2 2 + e + g h 2 2 = 3hu 2 2 + he ′ .
e(h) is then solution of the differential equation e -he ′ + gh 2 /2 = 0.

A solution is

e(h) = gh 2 2 .
In the end, we find

s(W ) = h u 2 2 + gh 2 2 , G(W ) = h u 3 2 + ugh 2 .
This allows us to compute the entropy variables (5.4)

W * 1 = gh - u 2 2
, W * 2 = u, and the reverse change of variables

h = 2W * 1 + W * 2 2 2g , u = W * 2 .
The equilibrium kinetic vectors are given by

F eq 1 = W 2 - Q(W ) 2λ , F eq 2 = W 2 + Q(W ) 2λ .
After some calculations, we can express this equilibrium in the entropy variables

F eq 1 = (W * 2 2 +2W * 1 )(λ-W * 2 ) 4gλ - (W * 2 2 +2W * 1 )(-4W * 2 λ+5W * 2 2 +2W * 1 ) 16gλ ⊺ , F eq 2 = (W * 2 2 +2W * 1 )(λ+W * 2 ) 4gλ (W * 2 2 +2W * 1 )(4W * 2 λ+5W * 2 2 +2W * 1 ) 16gλ ⊺ .
From the above theory, we know that

F eq k = ∇ W * s * k ,
for some dual kinetic entropies s * k . This is indeed the case and after more calculations we find

s * 1 = (λ -W * 2 ) W * 2 2 + 2W * 1 2 16gλ , s * 2 = W * 2 2 + 2W * 1 2 (λ + W * 2 ) 16gλ
.

It is then possible to compute the Hessians of s * k and express them in the (h, u) variables with (5.4). We find

D W * W * s * 1 =   λ-u 2gλ -gh+λu-u 2 2gλ -gh+λu-u 2 2gλ (gh+u 2 )λ-3hug-u 3 2gλ   , D W * W * s * 2 =   λ+u 2gλ gh+λu+u 2 2gλ gh+λu+u 2 2gλ (gh+u 2 )λ+3hug+u 3 2gλ   .
The two matrices are positive definite iff the first diagonal terms and the determinants are positive. This is equivalent to λ > |u| , (λ -u) 2 -gh > 0, (λ + u) 2 -gh > 0, which is again equivalent to λ > |u| + gh. This is the expected sub-characteristic condition. It is difficult to go farther because the Legendre transforms s 1 and s 2 of s * 1 and s * 2 cannot be computed explicitly. However, we can reproduce the stability property of the linear case. The microscopic entropy is given by

σ(F 1 , F 2 ) = s 1 (F 1 ) + s 2 (F 2 ).

Using the relations

W = F 1 + F 2 , y = -λF 1 + λF 2 -Q(F 1 + F 2 ), we deduce F 1 = W 2 - Q(W ) 2λ - y 2λ , F 2 = W 2 + Q(W ) 2λ + y 2λ ,
and the microscopic entropy can be expressed in function of W and y

σ(W, y) = s 1 W 2 - Q(W ) 2λ - y 2λ + s 2 W 2 + Q(W ) 2λ + y 2λ .
For a fixed W the minimum is achieved for y = 0, therefore the macroscopic entropy is

s(W ) = σ(W, 0),
and ∇ y σ(W, 0) = 0. Then, with a Taylor expansion near to y = 0, we get σ(W, y) = σ(W, -y) + O( |y|3 ).

The relation (5.3) thus still holds but with a third-order term in y. This means that the relaxation scheme with ω = 2 is entropy preserving up to third order in y. In principle, it is also possible to construct a scheme that preserves exactly the entropy in the non-linear case. It is sufficient to choose the relaxation parameter ω = ω(W, y) in such way that (5.5) σ(W, (1 -ω(W, y))y) = σ(W, y).

In practice, this would not be very interesting, one would get ω(W, y) ≃ 2

and ω(W, y) would have to be computed numerically by first computing s 1 and s 2 numerically and then by solving (5.5) also numerically. What is interesting, however, is that the reasoning ensures the existence of a relaxation parameter ω(W, y) ≃ 2, such that the whole scheme is entropy preserving. And if the scheme is run with a smaller relaxation parameter, it is ensured to be entropy stable. Similar ideas for controlling entropy dissipation in the over-relaxation step have been developed by Karlin in the standard LBM for Navier-Stokes. See [START_REF] Ilya V Karlin | Perfect entropy functions of the lattice Boltzmann method[END_REF] or [START_REF] Seyed | Entropic lattice Boltzmann methods: A review[END_REF] for a recent review article. See also [START_REF] Brownlee | Stability and stabilization of the lattice Boltzmann method[END_REF]. 5.3. Application to isothermal Euler equations. Finally, we also apply the theory to the isothermal Euler model, which reads

∂ t W + ∂ x Q(W ) = 0, with W = ρ ρu , Q(W ) = ρu ρu 2 + c 2 ρ .
In primitive variables v = (ρ, u) ⊺ , the system reads

∂ t v + B(v)∂ x v = 0, with B(v) = u ρ c 2 ρ u .
First, let us find a Lax entropy H(v) and an entropy flux R(v). We must have

(H ρ , H u ) u ρ c 2 ρ u = (R ρ , R u ).
We search for H(v) in the form

H = ρ u 2 2 + e(ρ). Then R ρ = u 3 2 + e ′ (ρ)u + c 2 u, R u = ρ 3u 2 2 + e ′ (ρ)ρ.
Integrating the first equation with respect to ρ we get

R = ρ u 3 2 + e(ρ)u + ρc 2 u + C(ρ),
(we can take C(ρ) = 0). And deriving with respect to u we get

ρ 3u 2 2 + e ′ (ρ)ρ = ρ 3u 2 2 + e(ρ) + ρc 2 .
Therefore e(ρ) is a solution of the differential equation

e ′ (ρ)ρ = e(ρ) + ρc 2 .
We can take e(ρ) = c 2 ρ(ln ρ -1).

Finally

s(W ) = ρ u 2 2 + c 2 ρ(ln ρ -1), G(W ) = u(s(W ) + c 2 ρ).
The entropy variables are

W * 1 = - u 2 2 + c 2 ln ρ, W * 2 = u.
The reverse change of variables is

ρ = W 1 = exp 2W * 1 + W * 2 2 2c 2 , ρu = W * 2 exp 2W * 1 + W * 2 2 2c 2 .
The equilibrium distribution is

F eq 1 = W 2 - Q(W ) 2λ = 1 2 W 1 -W 2 /λ W 2 -(W 2 2 /W 1 + c 2 W 1 )/λ . F eq 2 = W 2 + Q(W ) 2λ = 1 2 W 1 + W 2 /λ W 2 + (W 2 2 /W 1 + c 2 W 1 )/λ .
In the dual variables we get

∇s * 1 = exp 2W * 1 +W * 2 2 2c 2 2λ λ -W * 2 -W * 2 2 + λW * 2 -c 2 , ∇s * 2 = exp 2W * 1 +W * 2 2 2c 2 2λ λ + W * 2 W * 2 2 + λW * 2 + c 2 , and finally s * 1 = exp 2W * 1 +W * 2 2 2c 2 2λ (λ -W * 2 ), s * 2 = exp 2W * 1 +W * 2 2 2c 2 2λ (λ + W * 2 ).
With similar calculations as for the shallow water system we find the following sub-characteristic condition

λ > c + |u| .
As for the shallow water system, it is difficult to compute explicitly s 1 and s 2 .

Equivalent equation analysis

The entropy analysis of the over-relaxation scheme ensures the stability of the scheme as soon as the sub-characteristic condition is satisfied. However, for the moment it is not obvious that the scheme provides an approximation of the system (2.1). Indeed, the consistency of the kinetic approximation with the system of conservation laws is ensured, as soon as F (X, t) ≃ F eq (W (X, t)), and the over-relaxation formula (4.6) enforces F = F eq at the end of the time step only when ω = 1.

The objective of the equivalent system analysis is to provide a consistency theory in the case 1 < ω ⩽ 2. This consistency analysis is based on two ingredients: a Taylor expansion followed by a Chapman-Enskog analysis.

6.1. Taylor expansion: equivalent system in Y . The first ingredient is a Taylor expansion of S(∆t) -S -1 (∆t) 2∆t

when ∆t tends to zero. The Taylor expansion provides an equivalent system of Partial Differential Equations (PDE) expressed on Y . This system involves the whole vector Y , which contains both W and the flux error. We denote it by the equivalent system in the following. It takes the general form (6.1)

∂ t Y + r(ω) ∆t 0 Y m+1••• + 1⩽i⩽d A i (Y, ω)∂ i Y -∆t 1⩽i,j⩽d B i,j (Y, ω)∂ i,j Y = O(∆t 2 ).
In the case ω = 2, the formula is simpler and we shall find that (6.2) r(2) = 0, B i,j (Y, 2) = 0, and that the matrices A i are of the form

(6.3) A i (Y, 2) = D W Q i (W ) 0 m 0 m(nv-1) × .
The Taylor expansion is thus sufficient to get the second-order time consistency of the split scheme with the initial system of conservation law when ω = 2. This consistency holds even when the flux error is large. In addition, up to third order terms, the evolution of W is uncoupled from the evolution of the flux error. This surprising result relies essentially on the symmetry property (4.10), which ensures that the stiff terms in ∆t -1 vanish in the Taylor expansion.

6.2. Asymptotic analysis: equivalent equation in W . In the case ω ̸ = 2, the equivalent system contains stiff terms. And the evolution of W is no more uncoupled from the evolution of the flux error. In order to remove this coupling, the second step is to perform an additional asymptotic analysis (similar to the Hilbert or Chapman-Enskog expansion), with the assumption that the flux error is of order O(∆t). From (6.1) we can deduce an algebraic relation between the flux error

Y m+1••• and the gradient of W Y m+1••• = - ∆t r(ω) T 1⩽i⩽d A i (Y, ω)∂ i W 0 + O(∆t 2 ), T = 0 m 0 m(nv-1) 0 m(nv-1) 1 m(nv-1)
.

Reinjecting this approximation in the first row of (6.1) provides a simpler system involving only W . We denote it by the equivalent equation in the following. It takes the form (6.4)

∂ t W + 1⩽i⩽d ∂ i Q i (W ) -∆t 1⩽i,j⩽d ∂ i D i,j (W, ω)∂ j W = O(∆t 2 ).
6.3. Stability conditions. Once we have obtained the equivalent system (6.1) and the equivalent equation ( 6.4) we can study their stability.

6.3.1. Equivalent system. From the theory developed in Section (3), we have stability under the subcharacteristic condition. This condition ensures dissipation of the kinetic entropy σ

(Y ) = σ(W, Y m+1•••mnv ) = k s k (F k ).
We also know that σ(W, 0) = s(W ), corresponds to the minimum of σ with respect to the flux error Y m+1•••mnv . It is therefore expected that the change of variables Y → ∇ Y σ(Y ) symmetrizes the equivalent system (6.1), which is thus hyperbolic. We have also seen that the computation of the kinetic entropies s k is not necessarily easy. That is why in the following we try to find directly a symmetrization for the first order part of the equivalent system, rather than computing the kinetic entropy. We introduce the following definition.

Definition 9. We shall say that the equivalent system (6.1) is hyperbolic iff we can find a matrix P 0 (Y ), 0 ⩽ i ⩽ d, such that P 0 (Y ) is symmetric positive definite and

P 0 (Y )A i (Y, ω) is symmetric for 1 ⩽ i ⩽ d.
Remark 10. Our definition of hyperbolicity is stronger than the usual one, which only states that

d i=1 N i A i (Y, ω)
is diagonalizable with real eigenvalues for all directions N . We expect that

P 0 (Y ) = ∇ 2 σ(Y ),
at least in the case of the linear transport equation discussed in Section 5.1. But we don't know if the stability condition of Definition 9 is equivalent to the convexity of the dual kinetic entropies. Maybe that it is true. For a discussion around these questions, we refer to [START_REF] Bouchut | Stability of relaxation models for conservation laws[END_REF].

Equivalent equation.

For the equivalent equation, we already know that its first order part is hyperbolic. The stability thus depends on the sign of the second-order terms. This gives another stability criterion. For obtaining this criterion, we multiply (6.4) on the left by D W s(W ) and we integrate by part the second-order term. The entropy is dissipated if

ij,k,ℓ ∂ 2 k,ℓ s∂ i W k D i,j ∂ j W ℓ ⩾ 0.
We thus introduce the quadratic form acting on a second-order tensor α k i :

α → q(α) = i,j,k,ℓ ∂ 2 k,ℓ s(W )D i,j (W, ω)α k i α l j .
Definition 11. The equivalent equation (6.4) is dissipative iff the quadratic form q(α) is positive.

Remark 12. This definition amounts to checking that the Hessian matrix ∇ 2 α q(α) is positive. This (symmetric) matrix is of size md × md. In the scalar case m = 1, which we study below, the condition is simpler. It simply states that the quadratic form

x → i,j D i,j x i x j is positive.

Applications to the transport equation

Now we apply the equivalent system analysis and the equivalent equation analysis to the simple scalar transport equation (thus m = 1)

∂ t w + d i=1 ∂ i (v i w) = 0,
where the velocity vector (v 1 , . . . , v d ) is supposed to be constant. We consider the cases d = 1 or d = 2 and the D1Q2, D2Q3 and D2Q4 models.

For each model, we compute the equivalent system. The expansion of

S(∆t) -S -1 (∆t) 2∆t Y (•, t) = ∂ t Y (•, t) + O(∆t 2 ),
is performed with the Computer Algebra System (CAS) Maple. This is done by entering the explicit definition of the symmetric operator given in Section 4.2, step by step. Without a CAS, the calculations would be extremely tedious... 7.1. D1Q2. For the D1Q2 model, we use the notations

Y = w y , v 1 = v.
The equivalent system reads

∂ t w y - ω(ω -2)(ω 2 -2ω + 2) 2∆t(ω -1) 2 0 y + v γ 1 (λ 2 -v 2 )γ 1 -v(ω 4 -4ω 3 +6ω 2 -4ω+2) 2(ω-1) 2 ∂ x w y + -(ω 2 -6ω + 6)(λ 2 -v 2 ) 3v(ω 2 -2ω + 2) 3v(λ 2 -v 2 )(ω 2 -2ω + 2) -5v 2 ω 2 -3λ 2 ω 2 + 6v 2 ω + 10λ 2 ω -6v 2 -10λ 2 × ∆tω(ω -2) 32(ω -1) 2 ∂ xx w y = O(∆t 2 ), (7.1)
with

γ 1 = (ω-2) 2 (ω 2 -2ω+2) 8(ω-1) 2
. We can check that when ω = 2, we indeed obtain an equivalent system with the simplification (6.2) and (6.3). The stiff term vanishes and the evolution of w is uncoupled from that of y at order 2. This means that the consistency is achieved even when the flux error y is large.

For obtaining the equivalent equation we assume that we have y = O(∆t). Let us write y = ∆tỹ. We obtain

ω(ω -2)(ω 2 -2ω + 2) 2(ω -1) 2 ỹ = (λ 2 -v 2 )(ω -2) 2 (ω 2 -2ω + 2) 8(ω -1) 2 ∂ x w + O(∆t).
By simplifying, we have

(7.2) y = (λ 2 -v 2 )(ω -2) 4ω ∆t∂ x w + O(∆t 2 ).
By reinjecting this expression of y in the first equation of equivalent system 7.1, we obtain the equivalent equation on w

∂ t w + v∂ x w + (ω -2) 2 (ω 2 -2ω + 2) 8(ω -1) 2 (λ 2 -v 2 )(ω -2) 4ω ∆t∂ xx w - ∆tω(ω -2) 32(ω -1) 2 (ω 2 -6ω + 6)(λ 2 -v 2 )∂ xx w = O(∆t 2 ),
which can be simplified in

(7.3) ∂ t w + v∂ x w = 1 2 1 ω - 1 2 (λ 2 -v 2 )∆t∂ xx w + O(∆t 2 ).
We can notice that we recover the equivalent equation given in [START_REF] Dubois | Equivalent partial differential equations of a lattice Boltzmann scheme[END_REF][START_REF] Graille | Approximation of mono-dimensional hyperbolic systems: A lattice Boltzmann scheme as a relaxation method[END_REF][START_REF] Courtès | Vectorial kinetic relaxation model with central velocity. application to implicit relaxations schemes[END_REF].

Theorem 13. When 1 ⩽ ω < 2, the sub-characteristic diffusive stability condition of the D1Q2 model is |v| < λ.
Proof. The equivalent equation on w of the D1Q2 model (7.3) is stable if the diffusion term is positive. As ω ∈ [1, 2), the term 1 ω -1 2 is positive. The positivity of the diffusion term is then equivalent to

λ 2 -v 2 > 0,
which gives us the stability condition |v| < λ. □ Remark 14. When ω = 2, the diffusion term of the equivalent equation of the D1Q2 model disappears, which gives us ∂ t w + v∂ x w = O(∆t 2 ). We obtain that the solution given by the D1Q2 model is an approximation of order 2 of the solution of the initial equation.

Theorem 15. The matrix

P = 1 0 0 1 λ 2 -v 2
, symmetrizes the equivalent system of the D1Q2 model (7.1), if the diffusive sub-characteristic stability condition is satisfied. Consequently, the equivalent system (7.1) is hyperbolic if

|v| < λ.
Proof. We search a matrix P = p 1 p 2 p 2 p 3 such as P A is symmetric and P is symmetric positive definite. We have

P A = p 1 p 2 p 2 p 3 v γ 1 (λ 2 -v 2 )γ 1 -vγ 2 , = vp 1 + (λ 2 -v 2 )γ 1 p 2 γ 1 p 1 -vγ 2 p 2 vp 2 + (λ 2 -v 2 1 p 3 γ 1 p 2 -vγ 2 p 3 .
with

γ 1 = (ω-2) 2 (ω 2 -2ω+2) 8(ω-1) 2
, and

γ 2 = (ω 4 -4ω 3 +6ω 2 -4ω+2) 2(ω-1) 2
. As we want P A to be symmetric, we need to satisfy the condition

γ 1 p 1 -vγ 2 p 2 = vp 2 + (λ 2 -v 2 )γ 1 p 3 ,
which is equivalent to

p 3 = 1 (λ 2 -v 2 ) p 1 -v 1 + γ 2 (λ 2 -v 2 )γ 1 p 2 .
Let us choose p 2 = 0 and p 1 = 1. We obtain

P = 1 0 0 1 λ 2 -v 2 .
As its eigenvalues are 1 and 1 λ 2 -v 2 , P is definite positive if |v| < λ.

□

Remark 16. We obtain the same condition on v and λ as for the diffusive stability condition given in Proposition 13. In this case, the diffusive analysis and the hyperbolicity analysis give the same stability condition.

D2Q3.

For the D2Q3, we use the notations

Y =   w y 1 y 2   , v 1 = a, v 2 = b.
The equivalent system of the D2Q3 model is

∂ t   w y 1 y 2   - ω(ω -2)(ω 2 -2ω + 2) 2∆t(ω -1) 2   0 y 1 y 2   +   a -2γ 1 0 γ 1 (2a + λ)(a -λ) γ 2 (2a -λ) 0 γ 1 b(2a + λ) 2bγ 2 γ 2 λ   ∂ 1   w y 1 y 2   +   b 0 -2γ 1 γ 1 b(2a + λ) 0 γ 2 (2a + λ) γ 1 (aλ + 2b 2 -λ 2 ) γ 2 λ 2bγ 2   ∂ 2   w y 1 y 2   = O(∆t), (7.4)
with

γ 1 = -1 16 (ω 2 -2ω+2)(ω-2) 2 (ω-1) 2 and γ 2 = -1 4 ω 4 -4ω 3 +6ω 2 -4ω+2 (ω-1) 2
. For getting the equivalent equation, we assume that (y 1 , y 2 ) = O(∆t). As above, we express the flux error in function of the gradient of w up to order 2. This gives:

y 1 = ∆t 1 ω - 1 2 (2a + λ) 2 ((a -λ) ∂ 1 w + b∂ 2 w) + O(∆t 2 ),
and

y 2 = ∆t 2 1 ω - 1 2 (2ab + bλ) ∂ 1 w + λa + 2b 2 -λ 2 ∂ 2 w + O(∆t 2 ).
We reinject these expressions of y 1 and y 2 in the first equation of the equivalent system (7.4), we obtain (7.5)

∂ t w + ∂ 1 (aw) + ∂ 2 (bw) = ∆t 2 1 ω - 1 2 ∇ • (D 3 ∇w) + O(∆t 2 ),
with the diffusion matrix

D 3 = λ 2 (λ + a) -a 2 -λ 2 b -ab -λ 2 b -ab λ 2 (λ -a) -b 2 . Theorem 17.
The sub-characteristic stability condition of the D2Q3 model is

(7.6) λ 2 -a 2 -b 2 -(a 2 + b 2 ) 2 + λ(-2a 3 + 6ab 2 ) + λ 2 (a 2 + b 2 ) > 0.
Remark 18. This condition has a geometric interpretation. It also states that the velocity vector (a, b) has to be inside the triangle formed by the kinetic velocities, which are here (λ, 0), (-λ/2, λ √ 3/2) and (-λ/2, -λ √ 3/2). It is not easy to guess it from the inequality (7.6), but it becomes obvious if we plot numerically the stability region (see Figure 7.1).

Proof. Indeed, with a linear flux, we have

∂ t w + ∂ 1 (aw) + ∂ 2 (bw) = ∆t 2 1 ω - 1 2 ∇ • (D 3 ∇w) + O(∆t 2 ),
with the diffusion matrix

D 3 = λ 2 (λ + a) -a 2 -λ 2 b -ab -λ 2 b -ab λ 2 (λ -a) -b 2 .
The eigenvalues of this diffusion matrix are

d 1,2 = 1 2 λ 2 -a 2 -b 2 ± (a 2 + b 2 ) 2 + λ(-2a 3 + 6ab 2 ) + λ 2 (a 2 + b 2 ) .
Finally, the model D2Q3 is stable if D 3 is positive definite, namely if d 1 > 0 and d 2 > 0. □ Theorem 19. The matrix

P =   λ 2 (a 2 -2aλ -3b 2 + λ 2 )(2a + λ) 0 0 0 -(aλ + 2b 2 -λ 2 ) b(2a + λ) 0 b(2a + λ) -(a -λ)(2a + λ)   .
symmetrizes the equivalent system of the D2Q3 model (7.4), if the diffusive sub-characteristic stability condition (17) is verified. Consequently, the equivalent system (7.4) is hyperbolic if

λ 2 -a 2 -b 2 -(a 2 + b 2 ) 2 + λ(-2a 3 + 6ab 2 ) + λ 2 (a 2 + b 2 ) > 0.
Proof.

We are searching for a matrix P =   p 1 p 2 p 3 p 2 p 4 p 5 p 3 p 5 p 6   such as P A 1 and P A 2 are symmetric and P is symmetric positive definite. When we compute the matrices P A 1 and P A 2 , the symmetry imposes 6 equations on the unknown p 1 , p 2 , p 3 , p 4 , p 5 , p 6 . This gives us the matrix

P =    λ 2 (a 2 -2aλ-3b 2 +λ 2 )p5 b 0 0 0 -(aλ+2b 2 -λ 2 )p5 b(2a+λ) p 5 0 p 5 -(a-λ)p5 b    ,
where p 5 must be chosen. We choose p 5 = b(2a + λ). We obtain

P =   λ 2 (a 2 -2aλ -3b 2 + λ 2 )(2a + λ) 0 0 0 -(aλ + 2b 2 -λ 2 ) b(2a + λ) 0 b(2a + λ) -(a -λ)(2a + λ)   .
The eigenvalues of P are

e 1 = λ 2 (a 2 -2aλ -3b 2 + λ 2 )(2a + λ),
Y =     w y 1 y 2 z 3     , v 1 = a, v 2 = b.
We can also compute the equivalent system on (w, y 1 , y 2 , z 3 ) of the D2Q4 model. We obtain

∂ t     w y 1 y 2 z 3     - ω(ω -2)(ω 2 -2ω + 2) 2∆t(ω -1) 2     0 y 1 y 2 z 3     +     a 2γ 1 0 0 γ 1 (λ 2 -2a 2 ) -2aγ 2 0 γ 2 -2abγ 1 -2bγ 2 0 0 2λ 2 aγ 1 2λ 2 γ 2 0 0     ∂ 1     w y 1 y 2 z 3     +     b 0 2γ 1 0 -2abγ 1 0 -2aγ 2 0 γ 1 (λ 2 -2b 2 ) 0 -2bγ 2 -γ 2 -2λ 2 bγ 1 0 -2λ 2 γ 2 0     ∂ 2     w y 1 y 2 z 3     = O(∆t), (7.7)
with γ 1 = (ω-2) 2 (ω 2 -2ω+2)

16(ω-1) 2

and γ 2 = ω 4 -4ω 3 +6ω 2 -4ω+2 4(ω-1) 2

. With the same method as above, we derive the equivalent equation

∂ t w + ∂ 1 (aw) + ∂ 2 (bw) = ∆t 2 1 ω - 1 2 ∇ • (D 4 ∇w) + O(∆t 2 ),
with the diffusion matrix

D 4 = λ 2 2 -a 2 -ab -ab λ 2 2 -b 2 . Theorem 21. The D2Q4 model is stable if a 2 + b 2 ⩽ λ 2 2 . Proof. We have (7.8) ∂ t w + ∂ 1 (aw) + ∂ 2 (bw) = ∆t 2 1 ω - 1 2 ∇ • (D 4 ∇w) + O(∆t), with D 4 = λ 2 2 -a 2 -ab -ab λ 2 2 -b 2 .
The model is stable if the diffusion matrix D 4 is positive. Its eigenvalues are:

e 1 = 1 2 λ 2 -a 2 -b 2 -(λ 2 -a 2 -b 2 ) 2 -λ 4 + 2λ 2 b 2 + 2λ 2 a 2 and e 2 = 1 2 λ 2 -a 2 -b 2 + (λ 2 -a 2 -b 2 ) 2 -λ 4 + 2λ 2 b 2 + 2λ 2 a 2 .
As e 1 ⩽ e 2 , the eigenvalues are both positive if e 1 ⩾ 0, which means if

a 2 + b 2 ⩽ λ 2 2 .
□ Theorem 22. The matrix

P =     λ 2 (4a 2 -λ 2 )(4b 2 -λ 2 ) 0 0 0 0 -2λ 2 (4b 2 -λ 2 ) 0 2a(4b 2 -λ 2 ) 0 0 -2λ 2 (4a 2 -λ 2 ) -2b(4a 2 -λ 2 ) 0 2a(4b 2 -λ 2 ) -2b(4a 2 -λ 2 ) -2a 2 -2b 2 + λ 2     ,
symmetrizes the equivalent system of the D2Q4 model (7.7), if

(7.9) 4 max(a 2 , b 2 ) < λ 2 .
Consequently, under this condition, the equivalent system (7.7) is hyperbolic.

Proof. We are searching for a matrix such as P A 1 and P A 2 are symmetric and P is symmetric positive definite. We can compute P A 1 and P A 2 . As we want these matrices to be symmetric, we obtain conditions on the coefficients p i . We deduce that

P =    
P =      1 2a λ 2 (2a -λ)(2a + λ)p 7 0 0 0 0 -p 7 λ 2 /a 0 p 7 0 0 -p7λ 2 (2a-λ)(2a+λ) (a(2b-λ)(2b+λ)) -bp7(2a-λ)(2a+λ) (a(2b-λ)(2b+λ)) 0 p 7 -bp7(2a-λ)(2a+λ) (a(2b-λ)(2b+λ)) -1(2a 2 +2b 2 -λ 2 )p7 2(a(2b-λ)(2b+λ))     
.

By choosing p 7 = 2a(2b -λ)(2b + λ), we obtain

(7.10) P =     λ 2 (4a 2 -λ 2 )(4b 2 -λ 2 ) 0 0 0 0 -2λ 2 (4b 2 -λ 2 ) 0 2a(4b 2 -λ 2 ) 0 0 -2λ 2 (4a 2 -λ 2 ) -2b(4a 2 -λ 2 ) 0 2a(4b 2 -λ 2 ) -2b(4a 2 -λ 2 ) -2a 2 -2b 2 + λ 2     .
As P is symmetric, according to the Sylvester's criterion, P is positive definite if and only if all the leading principal minors are positive, that is to say if the following conditions are satisfied

       |P 1 | = λ 2 (4b 2 -λ 2 )(4a 2 -λ 2 ) > 0, |P 2 | = -2λ 4 (4a 2 -λ 2 )(4b 2 -λ 2 ) 2 > 0, |P 3 | = 4λ 6 (4a 2 -λ 2 ) 2 (4b 2 -λ 2 ) 2 > 0, |P 4 | = 4λ 4 (4a 2 -λ 2 ) 3 (4b 2 -λ 2 ) 3 > 0.
This is equivalent to

4a 2 < λ 2 , 4b 2 < λ 2 ,
which can be rewritten 2 max (|a|, |b|) < λ. □ Remark 23. The hyperbolicity condition obtained is more restrictive than the diffusive stability condition obtained in Proposition 21. We can see in Figure 7.2 the values of a/λ and b/λ for which the diffusive stability condition is verified, the circle colored in yellow, are included in the blue square, for which the hyperbolicity condition is checked. This is coherent with the review of stability conditions given by Bouchut in [START_REF] Bouchut | Stability of relaxation models for conservation laws[END_REF].

Some numerical results

8.1. D1Q2: consistency error. Now that we have obtained the equivalent equations, we wish to quantify numerically how they are close to the kinetic equations. We shall compute analytic solutions of the equivalent equations and compare them with the solutions of the kinetic equation. We shall also compute the error between the two solutions.

For the numerical experiments, we use particular solutions of the form (8.

+ ivkw = - ∆t 2 1 ω - 1 2 k 2 λ 2 -v 2 w.
It gives us the value of γ with respect to k and v

γ = - ∆t 2 1 ω - 1 2 k 2 λ 2 -v 2 -ivk.
A particular solution of the equivalent equation (7.3) is then

w = w 0 e -( ∆t 2 ( 1 ω -1 2 )k 2 (λ 2 -v 2 )+ivk)t e ikx .
In order to deal with real solutions, we compute the real part of this particular solution, that we denote w eqeq and which is still a solution of (7.3)

w eqeq = ℜ(w) = w 0 e -∆t 2 ( 1 ω -1 2 )k 2 (λ 2 -v 2 )t cos (k(x -vt)) .
To compute the equivalent equation, we assume that we have the relation between y and ∂ x w given by (7.2)

(8.2) y = (λ 2 -v 2 )(ω -2) 4ω ∆t∂ x w.
We denote y eqeq the real part of y

y eqeq = ℜ(y), = - (λ 2 -v 2 )(ω -2) 4ω ∆tkw 0 e -∆t 2 ( 1 ω -1 2 )k 2 (λ 2 -v 2 )t sin (k(x -vt)) .
8.1.2. Particular solution of the equivalent system. Now, we inject the expression of the particular solution (8.1) in the equivalent system (7.1). We obtain

γI 2 + 1 ∆t R(ω) + ikA(Y, ω) + ∆tk 2 B(Y, ω) w y = 0,
with R(ω) = 0 0 0 r(ω) .

The previous system admits two solutions γ 1 (k) and γ 2 (k) depending on k, which are the eigenvalues of -1 ∆t R(ω) -ikA(Y, ω) -∆tk 2 B(Y, ω). We have

γ 1 (k) = 1 ∆t 16ω 4 -64ω 3 + 96ω 2 -64ω 32(ω -1) 2 + O(∆t 0 ), and
γ 2 (k) = -ikD(W, ω) - D(W, ω) 2 (λ 2 ω 2 -2) + 2λ 2 4ω k 2 ∆t + O(∆t 2 ).
One of the solutions, γ 1 (k), behaves as O( 1 ∆t ) when ∆t → 0, and the real part of the other solution γ 2 (k) behaves as O(∆t) when ∆t → 0. If we compute the particular solution (8.1) with the eigenvalue γ 1 in O( 1 ∆t ), we observe that y decreases rapidly toward 0. If we consider instead, the solution given by the second eigenvalue γ 2 , y stays small and has slower variations. We choose to keep this eigenvalue γ 2 for a relevant comparison with the expected behavior.

A particular solution of the equivalent system (7.1) is then

w y (x, t) = w 0 y 0 e γ2t e ikx .
To test, we compute the real part of w, that we denote w syseq w syseq = ℜ(w) = w 0 e ℜ(γ2)t cos (ℑ(γ 2 )t + kx) , and the real part of y, denoted by y syseq y syseq = ℜ(y) = e ℜ(γ2)t (ℜ(y 0 ) cos (ℑ(γ 2 )t + kx) -ℑ(y 0 ) sin (ℑ(γ 2 )t + kx)) .

8.1.3. Numerical comparison of w. We take k = 2. We choose w 0 = 1, and we take y 0 such as (w 0 , y 0 ) belong to the kernel of the matrix -1 ∆t R(ω) -ikA(Y, ω) -∆tk 2 B(Y, ω). We denote w LB the solution given by the D1Q2 model with the initialization w LB y LB (x, 0) = w 0 y 0 cos(kx).

We compute the relative L 2 error between the solution of the equivalent equation w eqeq and the solution given by the D1Q2 model w LB at the final time

N x i=0 w i,N t LB -w i,N t eqeq 2 N x i=0 w i,N t LB 2 ,
and the relative L 2 error between the solution of the equivalent system w syseq and w LB

N x i=0 w i,N t LB -w i,N t syseq 2 N x i=0 w i,N t LB 2 .
We compute the solution for different amounts of time steps N t = 16, 32, 64, 128, 256, 512, 1024 and 2048, which gives us different time steps ∆t = T N t , with T = π. We obtain the relative errors of Figure 1, for different relaxation parameters ω. The equivalent equation and the equivalent system both converge at the order 2 toward the solution given by the D1Q2 model. When ω ∈ [1.8, 2], the equivalent equation and the equivalent system give similar accuracy. When ω ∈ [1.5, 1.8], the equivalent system is a better approximation of the solution given by the D1Q2 model, while when ω ⩽ 1.4, the equivalent equation is more accurate. 8.1.4. Numerical comparison of y. Now, we want to compute the error on the flux error y.

We can compute the relative L 2 errors between y LB and the flux error y eqeq that we assume to have in order to compute the equivalent equation and between y LB and the solution of the equivalent system y syseq N x i=0 y i,N t LB -y i,N t We obtain the Figure 2. We can observe that the flux error y eqeq given by the equivalent equation converges at the order 1 toward the y LB given by the D1Q2 model, while the y syseq given by the equivalent system converges at the order 2.

Remark 24. When ω = 2, the error between the flux error y given by the equivalent equation and the one given by the D1Q2 model is constant. This is due to the fact that y eqeq is given by (8.2), which is equal to 0 when ω = 2. Indeed, when ω = 2, w and y are independent, so we do not have to assume the smallness hypothesis y = O(∆t) to deduce the equivalent equation from the equivalent system. 8.2. D2Q4: numerical stability. As we can see in Figure 7.2, for some choice of velocity (v 1 , v 2 ) = (a, b) and norm of the kinetic velocity λ, the diffusive stability condition can be satisfied, but not the hyperbolicity condition. We want to test numerically what happened when we are in this case. Table 1. Relative L 2 error on w with respect to the time step ∆t, for different relaxation parameters ω.

We consider a square geometry [0, 1] × [0, 1] with periodic boundary conditions. We consider N x = 200 space steps in both directions. We denote by ∆x = 1/N x the grid step. We initialized w with a Gaussian function centered in the middle of the square w(x, y, 0) = e -80((x-0.5) 2 +(y-0.5) 2 ) .

Let us choose (a, b) = (1, 0). The stability condition is satisfied if

λ > 2(a 2 + b 2 ) = √ 2.
The hyperbolicity condition is satisfied if

λ > 2 max(|a|, |b|) = 2.
We are going to compare the solution obtained with λ = 1.6, that is when the diffusive stability condition is satisfied, but not the hyperbolicity condition, and λ = 2.2, namely when both the diffusive stability and the hyperbolicity conditions are satisfied. We draw the solutions w(x, y, T ) at time T = 1, for different values of the relaxation parameter: ω = 1.2, ω = 1.6 and ω = 2.

As we are solving the transport step of time step ∆t 4 with a Lattice-Boltzmann method, we need to have the relation between the time and space step ∆t = 4∆x λ = 4 λN x .

Consequently, the number of time step is

N t = T ∆t = λN x 4 ,
and which depends on the λ chosen: we do N t = 80 time steps when λ = 1.6 and N t = 110 steps when λ = 2.2. Table 2. Relative L 2 error on y with respect to the time step ∆t, for different relaxation parameters ω.

We obtained the solution w of the Table 3. When λ = 2.2, that is to say when both the diffusive stability and the hyperbolicity conditions are verified, we obtained a Gaussian centered in the middle of the square, as expected. However, the closer the relaxation parameter ω is to 1, the more the Gaussian function dampens due to the relaxation step. When ω = 2 or ω = 1.6, the solutions obtained with λ = 1.6, namely when the diffusive stability condition is satisfied but not the hyperbolicity condition, are not stable. Oscillations appear and grow over time. When ω = 1.2 and λ = 1.6, we obtained a solution close to the expected Gaussian function, but a little distorted. Moreover, this solution is stable, we do not observe any oscillations.

Conclusion

In this work we have provided a general methodology for studying the stability and the consistency of the Vectorial Lattice-Boltzmann Method (VLBM). We have first shown that the dual entropy analysis of [START_REF] Bouchut | Construction of BGK models with a family of kinetic entropies for a given system of conservation laws[END_REF][START_REF] Dubois | Simulation of strong nonlinear waves with vectorial lattice Boltzmann schemes[END_REF] can be applied for a direct and rigorous proof of the stability of the over-relaxed time splitting algorithm. It is not necessary to pass through the stiff relaxation intermediary.

Secondly, we have proposed an automatic way to construct an equivalent system of PDE, consistent with the VLBM. This equivalent system contains stiff terms in ∆t. The classical equivalent equation can be derived from the equivalent system by a Chapman-Enskog analysis when ∆t is small and the kinetic data close to equilibrium. It seems, but it is a conjecture, that the hyperbolicity condition of the equivalent system is exactly the entropy stability condition.

In future works, we plan to investigate this conjecture, perform additional numerical experiments on truly non-linear systems, beyond the simple transport equation. Finally, an important question is to extend the stability analysis for handling boundary conditions correctly (see [START_REF] Hélie | Schéma de relaxation pour la simulation de plasmas dans les tokamaks[END_REF] for preliminary results). 3. Solutions w at time T = 1 for different values of λ and ω.

e 2 =

 2 λ 2 -a 2 -b 2 + (a 2 + b 2 ) 2 + λ(-2a 3 + 6ab 2 ) + λ 2 (a 2 + b 2 ), and e 3 = λ 2 -a 2 -b 2 -(a 2 + b 2 ) 2 + λ(-2a 3 + 6ab 2 ) + λ 2 (a 2 + b 2 ). By noticing that e 2 > e 3 and e 2 e 3 = 2e 1 , we deduce that P is definite positive if e 3 > 0. □ Remark 20. The hyperbolicity condition on a, b and λ is the same as the diffusive stability condition given in the Proposition 17. Here again, the diffusive analysis and the hyperbolicity analysis are equivalent.

Figure 7 . 1 .

 71 Figure 7.1. Graphic representation of the diffusive stability and hyperbolicity condition of the D2Q3 model for λ = 1. The stable region is the blue triangle, whose vertices are the kinetic velocities.

Figure 7 . 2 .

 72 Figure 7.2. Graphic representation of the diffusive stability and hyperbolicity condition of the D2Q4 model.

  p 1 p 2 p 3 p 4 p 2 p 5 p 6 p 7 p 3 p 6 p 8 p 9 p 4 p 7 p 9 p 10

	
	   ,

  Particular solution of the equivalent equation. If we inject this particular solution (8.1) in the equivalent equation (7.3) on w, we obtain γw

	1)	w y	(x, t) =	w 0 y 0	e

γt e ikx , with k ∈ N and γ ∈ C.

8.1.1.


Leila Galai

Emilie Bérard

Solenn Reguer

Eddy Foy

Cristian Mocuta

Sebastian Schöder

Enrique Vega

Caroline Toffolon-Masclet

Thomas Guilbert

Ivan Guillot

New experimental approach for the metallurgical study of medieval ferrous artefact: the case of arms and armour manufacturing

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

New experimental approach for the metallurgical study of medieval ferrous artefact: the case of arms and armour manufacturing

Content

The metallurgical investigations of medieval armours are of great interest the history of techniques [1]. Studying the nature and microstructure of the alloys provides crucial insights into the manufacturing processes used by the armourers. Indeed, the nature and heat treatment of the metal play a key role on its mechanical properties, and thus on the effectiveness of an armour. This approach has already brought new perspectives in the study of medieval armaments. However, metallographic studies require to invasively collect samples, which is rarely allowed on museum artefacts. Moreover, it only provides localized information that could be not representative of the whole armour which microstructure can be heterogenous.

For this reason, non-invasive approaches have been developed. X-ray diffraction (XRD) using Synchrotron Radiation (SR) is a well-suited tool to identify and quantify mineral phases, commonly found in ancient ferrous alloys, like cementite, ferrite and martensite [2]. The high flux, the good collimation and monochromaticity of the photon beam favors the detection of minor phases with low volume fraction (such as cementite) and can provide information on residual strain and texture. Combined with the fast diffractograms acquisition, it allows to multiply the analysis over extended regions of the samples to assess the heterogeneity of the metal on large museum artefact.

The aim of the present study is to set up a specific experimental approach for the analysis of museum armour. It has involved, in a first step, the fabrication and characterization of Fe-C model alloys at targeted carbon concentration ([C] < 0.8 wt%) and typical of microstructure found in ancient ferrous alloys such as perlite, ferrite, martensite and bainite. Thermal treatment including different cooling rates and monitored by quenching dilatometry were applied to the alloys to obtain the desired microstructure. Then, of characterization techniques (Optical Microscopy, Backscattering electron (BSE) imaging using the scanning electron microscope (SEM) and microhardness measurements) was carried out to identify and attest the obtained microstructure. SR-XRD measurements were carried out at the DiffAbs and PUMA beamline (SOLEIL Synchrotron) on the model alloys at different energies to assess the in-depth heterogeneity of the metal, depending on the microstructure. The measurements made on the model samples were then used to select the best operating condition to perform SR-XRD analysis on complete piece of armour. In parlicular, these measurements concern an Armour made by Valentin Siebenbürger's workshop. This armour belongs to the "Musée des Beaux Arts de Rennes"(Rennes, France, Inv. n°1096/4850 et 4851 DEP) and is currently on display at the "Musée de l'Armée"(Paris, France).

The data collected on the model alloys and late medieval European armour preserved at the Musée de l'Armée (Paris, France) will contribute to the elaboration of a reference table on the ferrous alloys used in the field of medieval armament but also for other ancient ferrous artefacts. This will bring new insights on the technical gestures and manufacturing processes and will pave the way for large-scale interdisciplinary in this field.