Damien David
email: damien.david@insa-lyon.fr

Frederic Kuznik

Kevyn Johannes

Lucie Merlier

Numerical analysis of truncation error, consistency, and axis boundary condition for axis-symmetric flow simulations via the radius weighted lattice Boltzmann model

Keywords: Lattice Boltzmann models, Axis-symmetric flows, Truncation Error Analysis, Boundary conditions

This paper is the continuity of Guo's work [Phys.Rev., 2009, 046708] about a lattice Boltzmann model for axis-symmetric flows. This is a radius weighted LBM model: all the moments are proportional to the radial coordinate r. A Taylor series analysis is performed on the discrete Boltzmann Model in order to assess the consistency. The truncation error terms do not indicate any error increase along the radial direction, but they contain a spurious term u r /r which considerably reduces the model accuracy when radial velocity is non-null: simulation results highlight a reductions of the order of accuracy (from second to first). A new axis boundary condition was created to improve the model accuracy and to cope with the spurious error term. It was assessed though several test cases.

The lattice Boltzmann method (LB method) is a well understood method for the simulation of fluid flows. Although it is limited to second order accuracy, it can compete with traditional fluid simulation methods in terms of computational performances, and simplicity of the implementation.

Historical origins of the method are found in the kinetic theory: as a continuity of Lattice Gas Cellular Automata (Frisch et al.

[1], McNamara and Zanetti

[2]), or as a special discretization of the Boltzmann equation (He and Luo [START_REF] He | Theory of the lattice boltzmann method: From the boltzmann equation to the lattice boltzmann equation[END_REF]).

The use of lattice Boltzmann models is not restricted to the simulation of the microscopic behaviour of fluids. The method is also a powerful numerical tool for the simulation of partial differential equations. Thus, LB models exist for various differential equations: Full Navier Stokes equations (Qian et al. [START_REF] Qian | Lattice bgk models for navier-stokes equation[END_REF]),

Poisson equation (Chai and Shi [START_REF] Chai | A novel lattice boltzmann model for the poisson equation[END_REF]), Fourier equation (Mishra [START_REF] Subhash | Application of the lattice boltzmann method for solving the energy equation of a 2-d transient conduction-radiation problem[END_REF]), etc... For each LBM model, it is essential that the consistency is expressed at least in one journal paper. Moreover, a truncation error analysis might be helpful for the determination of the optimal model parameters, and for the interpretation of the model behaviour under certain circumstances.

The consistency is assessed when the numerical equations tends toward the appropriate partial differential equations as the discretization parameters ∆X and ∆t decrease. The consistency of LB models is traditionally established through the Chapman Enskog expansion. This method is an heritage from the kinetic theory: it involves a scale separation (between acoustic and diffusive transport) on the distribution functions and the differential operators. Although the method is commonly used for recovering the zero order terms of the simulated equations, it might become a very long and painful experience to use it for recovering the truncation error terms unless restricted scaling hypothesis are adopted (Zhao [7]).

An alternative method was proposed by Holdych [START_REF] David | Truncation error analysis of lattice boltzmann methods[END_REF]: the so-called "Taylor series analysis". The method is an extension of Junk's work [START_REF] Junk | A finite difference interpretation of the lattice boltzmann method[END_REF] for τ = 1 to any value of the relaxation parameter. It does not involve any scale separation that might be confusing for LBM neophytes. It offers a clear an simple technique to get to the zero order equations, and does not require much more effort to reach the whole expression of truncation error terms. Despite its inherent simplicity, the method is still not popular among the LBM community (Ginzburg [START_REF] Ginzburg | Truncation errors, exact and heuristic stability analysis of two-relaxation-times lattice boltzmann schemes for anisotropic advectiondiffusion equation[END_REF]).

The present work was motivated by the need of a precise and reliable axissymmetric flow model, to simulate low-Rayleigh flows inside some cylindrical capsules that contain phase change materials. The model which was selected is the model proposed by Guo [START_REF] Guo | Theory of the lattice boltzmann equation: Lattice boltzmann model for axisymmetric flows[END_REF]. Most of the axis-symmetric LB models were obtained by modifying existing Cartesian LB models in order to recover the additional terms in cylindrical Navier-Stokes equations. Halliday et al. [START_REF] Halliday | Lattice boltzmann equation hydrodynamics[END_REF],

Reis and Philipps [START_REF] Reis | Modified lattice boltzmann model for axisymmetric flows[END_REF][START_REF] Reis | Erratum: Modified lattice boltzmann model for axisymmetric flows[END_REF], and Chen et al. [START_REF] Chen | Lattice boltzmann model for incompressible axisymmetric flows[END_REF][START_REF] Chen | Simulation of buoyancydriven flows in a vertical cylinder using a simple lattice boltzmann model[END_REF][START_REF] Chen | Numerical simulation of fluid flow and heat transfer inside a rotating disk-cylinder configuration by a lattice boltzmann model[END_REF][START_REF] Chen | Simple lattice boltzmann subgrid-scale model for convectional flows with high rayleigh numbers within an enclosed circular annular cavity[END_REF] used virtual source terms. Zhou [START_REF] Guo | Axisymmetric lattice boltzmann method[END_REF][START_REF] Guo | Axisymmetric lattice boltzmann method revised[END_REF] used a radius-dependent relaxation time and a virtual source term.

The way Guo built his model is radically different. The author started from the continuous Boltzmann equation in cylindrical coordinates from Bergers [START_REF] Bergers | Kinetic model solution for axisymmetric flow by the method of discrete ordinates[END_REF].

He discretized the equations and obtained a new kind of LB model, where the radius coordinate appears in all the moments definitions. From now, we refer to this model as the Radius-Weighted Lattice Boltzmann model: RW-LB model.

Despite subsequent extensions of Guo's model to thermal flows (Zheng et al. [START_REF] Zheng | Lattice boltzmann equation for axisymmetric thermal flows[END_REF]), and multi-relaxation time collision operator (Wang et al. [START_REF] Wang | Multi-relaxation-time lattice boltzmann model for axisymmetric flows[END_REF]), the consistency of the discrete model has never been clearly assessed. Guo [START_REF] Guo | Theory of the lattice boltzmann equation: Lattice boltzmann model for axisymmetric flows[END_REF] applied the Chapman Enskog expansion to the continuous Boltzmann equations from which he built the discrete RW-LB model. Moreover, the presence of the radial coordinate in the moments definitions and the force term might modify the truncation error terms: its effect on the model precision must be clearly stated. This is the first purpose of the present paper, which is treated through a Taylor series analysis.

The second purpose of the present paper concerns the boundary condition near the symmetry axis.

Macroscopic equations

Cylindrical formulations of the Navier-Stokes equations are found in every fluid dynamics textbook. The following formulation facilitates the comparison with the macroscopic equations obtained through the Taylor series analysis.

∂ β (ru β ) = 0 (1) ∂t(ruα) + ∂ αβ (ruαu β ) = - 1 ρ ∂α(rp) + p ρ δαr . . . • • • + ν ∂ ββ (ruα) -∂r(uα) - uα r δαr + rF0α ρ (2)
In the present paper, we adopt the Einstein formulation for all the equations.

In equations 1 and 2, the summation index is β and α is a free index.

RW-LB model

The RW-LB model is a D2Q9 BGK model, where the moments and the equilibrium distribution functions are multiples of the radial coordinate r. The simulation domain is discretized into square grids. The grid size is ∆X. The time step is ∆t. Each grid point is connected with itself and the eight neighbouring nodes through nine particle velocities:

-→ c 0 = (0; 0)∆X/∆t -→ c 1-4 = (0; ±1)∆X/∆t, (±1, 0)∆X/∆t -→ c 5-9 = (±1; ±1)∆X/∆t (3)
The distribution functions f 0-8 follow the governing equation:

f k ( -→ x + -→ c k ∆t, t + ∆t) -f k ( -→ x , t) = - 1 τ (f k -f eq k ) + ∆t 1 - 1 2τ R k (4)
The discrete relaxation time τ is related to the viscosity of the fluid through

ν = c 2 s ∆t(τ -0, 5), where c 2 s = ∆X 2 / √ 3∆t 2 is the sound speed.
f eq k is the discretized Maxwellian equilibrium distribution function multiplied by r. It is given by Equation 5, where ω k is a weighting which depends on the associated particle velocity: ω 0 = 4/9, ω 1-4 = 1/9, ω 5-9 = 1/36.

f eq k = ρrω k 1 + -→ c k . -→ u c 2 s + ( -→ c k . -→ u ) 2 2.c 4 s - -→ u 2 2.c 2 s ( 5
)
The distribution functions f k are related to the macroscopic variables through the first and second order moments:

k f k = rρ (6) k c kα .f k = r ρu α -∆t F α 2 (7)
For the source term R k , two formulations are available. The first formulation (Equation 8) is the He-type formulation (He et al. [START_REF] He | Discrete boltzmann equation model for nonideal gases[END_REF]) proposed by Guo [START_REF] Guo | Theory of the lattice boltzmann equation: Lattice boltzmann model for axisymmetric flows[END_REF] for the original RW-LB model. The second formulation (Equation 9) is the Guo-type formulation from Guo et al. [START_REF] Guo | Discrete lattice effects on the forcing term in the lattice boltzmann method[END_REF]. Since both formulations have the same zero-th, first, and second order moments, they lead to the same macroscopic equations and second order truncation error terms . Both formulations are equivalent for this study. The second formulation is chosen to get lighter expressions.

R He k = ( -→ c k --→ u ) . -→ F ρrc 2 s f eq k He-type (8)
R k = ω k -→ c k --→ u c 2 s + -→ c k . -→ u c 4 s -→ c k . r -→ F Guo-type (9)
The force F α is the sum of the "real" force applied to the fluid F 0α , and a corrective pseudo-force required to recover the cylindrical Navier Stokes equations

F ′ α : F α = F 0α + F ′ α = F 0α + ρc 2 s r -ρν 2u r r 2 δ αr (10)
Please note that Equation 10 is an exact reformulation of the force term given by Guo. The formulation given by Wang et al. [START_REF] Wang | Multi-relaxation-time lattice boltzmann model for axisymmetric flows[END_REF] and Zheng et al. [START_REF] Zheng | Lattice boltzmann equation for axisymmetric thermal flows[END_REF] is wrong: they confused the continuous model relaxation time τ * = ν/c 2 s with the discrete model relaxation time τ = τ * /∆t + 0.5

Consistency and truncation error

Scaling strategy

In LB models, the quantities ν, τ , ∆X and ∆t are connected through the relation:

ν = ∆X 2 3∆t 2τ -1 2 (11)
Because of this relationship, it is impossible to modify ∆X (or ∆t) without modifying other parameters of the model. Thus, whenever a study involves mesh refinement of a LBM model (error analysis, consistency, local mesh refinement...), it is important to clarify the condition of the mesh refinement: which parameter remains constant? How does ∆t vary with ∆X?

The choice of the constant parameter concerns two reasonable options: either the relaxation time τ is constant, or the Mach number M a is constant. Both parameters participate to the stability and the accuracy of the model. It is now common knowledge that acceptable values of τ lies in a narrow range between 0.5 and 2, and that the Mach number should be as low as possible to avoid compressibility effects.

The literature proposes two common scaling options. Let ∆X = 1/N X be the reference small parameter. Acoustic scaling involves ∆t ∝ ∆X, and diffusive scaling involves ∆t ∝ ∆X 2 (Asinari et al. [START_REF] Asinari | Link-wise artificial compressibility method[END_REF], Junk [START_REF] Junk | A finite difference interpretation of the lattice boltzmann method[END_REF]).

The different mesh refinement strategies were evaluated by computing the following quantities: reference velocity U 0 , model viscosity ν, Mach number M a, relaxation parameter τ and number of time steps N t . Those calculation were based on a fixed Reynolds number Re and non dimension simulated time κ:

M a = U 0 ∆X/( √ 3∆t) Re = U 0 N X ∆X ν κ = U 0 N t ∆t N X ∆X (12)
When the Mach number is fixed, we get τ -0, 5 ∝ N X . We might reach unreasonably high values of τ by refining the mesh. On the contrary, by fixing τ , we get M a ∝ 1/N X : we reduce the Mach number value by refining the mesh.

That's why we choose to fix the value of τ during the mesh refinement.

Since the aim of the paper is not to compare the relative order of convergence obtained through either acoustic or diffusive scaling, there was no such a strict criterion for choosing between the scaling strategies. None of them leads to a faster convergence, we always get

N t = κN 2 X Re(τ -0.5) (13)
We chose diffusive scaling based on recommendations from Holdych et al. [START_REF] David | Truncation error analysis of lattice boltzmann methods[END_REF]. With this scaling strategy, the velocity does not depend on N X , and the viscosity can be adjusted through the ratio ∆X 2 /∆t (Equation 11).

Taylor series analysis

A rigorous application of the Taylor series analysis to the RW-LB model is detailed in appendix Appendix A. Here, we present the major results of the calculations.

The Taylor series analysis leads to the modified equations of the numerical scheme. For the continuity equation, we get:

∂ α (ru α ) = (2τ -1) 2 12ν ∆X 2 - r ρν (∂ t p + u α ∂ α p) + ν∂ r u r r + O(∆X 4 ) (14)
and for the momentum equation, we get:

0 = -∂ t (ru α ) -∂ β (ru α u β ) - ∂ α (rp) ρ + p ρ δ αr . . . • • • + ν ∂ ββ (ru α ) -∂ r (u α ) -δ αr u r r + r F 0α ρ + ∆X 2 M E + O(∆X 4 ) (15)
where the momentum truncation error M E is expressed in Equation A.37 of appendix Appendix A.

The zero order elements of Equations If the radial velocity u r is null, the error terms depending on u r /r are null and there is no singularity near the axis. However, with non-zero radial velocities, the term u r /r reaches extreme values near r = 0, the error terms can not be neglected along the symmetry axis. For that reason, we conclude that the consistency of the model is assessed with no restriction if u r is null; it is assessed in all the model except at the axis location when u r is not null. This restriction is not limiting: axis-symmetric simulations are performed on a half-section of the cylindrical domain, and the symmetry axis is treated with proper boundary conditions. The accuracy of the model near the axis depends on the quality of the boundary condition.

Finally, let's consider the evolution of the total mass contained in the simulation domain. The mass is locally preserved during collisions, because the zero order momentum of the collision operator is null. During the streaming process, the quantity ρr is propagated: the total mass being proportional to i,j (ρr) ij , it is exactly preserved by the streaming operator. Thus, despite the presence of truncation error terms in the continuity equation, the total mass of the system is exactly preserved by the model, provided that the boundary conditions are mass conservative.

Axis boundary condition

The axis boundary condition is performed on the axis boundary nodes.

Those boundary nodes are located a half step ahead from the symmetry axis. To analyse the dynamic that is modelled by the axis boundary condition, we define the ghost nodes as the nodes located symmetrically to the boundary nodes with respect to the axis (Figure 1). Thus, the radius is r = ∆X/2 on the boundary nodes, it is r = -∆X/2 on the ghost nodes.

The boundary condition should act as if the macroscopic flow was exactly symmetric on both sides of the axis. Let ρ, (u r , u z ), and (a r , a z ) be the density, velocity, and acceleration on a boundary node. Let ρ, (ũ r , ũz ), and (ã r , ãz ) be the same quantities on the corresponding ghost node. The boundary condition should act as if those quantities were equal to ρ, (-u r , u z ), and (-a r , a z ) (Figure 2).

The symmetry of the acceleration is already ensured by the definition of the axis force in Equation 10. Indeed the real force F 0α is symmetric by definition, and the virtual force F ′ r is an odd function of the radius r and radial velocity u r . If the symmetry is ensured for the density and the velocity, we get ãr = -a r , and ãz = a z .

The boundary condition is applied during the streaming process. Before the streaming process, distribution functions on the boundary node are post-

Symmetry axis

Simulation domain Ghost nodes j j+1 j-1

               f * k = ∆X 2 ρ c kr f * k = ∆X 2 ρ u r + ∆t 2 a r c kz f * k = ∆X 2 ρ u z + ∆t 2 a z (16)

Mirror boundary conditions

Guo [START_REF] Guo | Theory of the lattice boltzmann equation: Lattice boltzmann model for axisymmetric flows[END_REF] proposes a "mirror" type boundary condition (mirror BC) at the axis location. The mirror boundary condition allocates the outgoing distribution functions f * 3 , f * 6 and f * 7 to the surrounding unknown distribution functions (in the nodes j -1, j and j + 1) as follow: The mirror boundary condition acts as if, before the streaming, the ghost node distribution functions were a rotation around the vertical axis of tbhe boundary node distribution functions. This behaviour is represented in Figure 3. Now we can compute the macroscopic quantities ρ and (ũ r , ũz ) at the ghost node location before the streaming process. Recalling that the radius is r = -∆X/2, we get:

           f + 5 j+1 = (f * 6 ) j f + 1 j = (f * 3 ) j f + 8 j-1 = (f * 7 ) j (17)
                 ρ = f * k -∆X/2 = -ρ ũr + ∆t 2 ãr = -c kr f * k (-∆X/2)ρ = -u r + ∆t 2 a r ũz + ∆t 2 ãz = c kz f * k (-∆X/2)ρ = u z + ∆t 2 a z ( 18
)
The symmetry is broken for the density. Although the ghost node velocities are correct, the mirror BC simulates opposite densities on both sides of the symmetry axis. It is not suited for the RW-LB model

The radius weighted mirror boundary condition

The radius weighted mirror boundary condition (RW Mirror BC) is a modified version of the mirror boundary condition which ensures the symmetry of the whole set of macroscopic quantities. To do so, we multiply by -1 the value of the distribution functions that cross the boundary;

             f + 5 j+1 = (-f * 6 ) j f + 1 j = (-f * 3 ) j f + 8 j-1 = (-f * 7 ) j (19)
The RW-mirror BC acts as if, before the streaming, the ghost node distribution functions were a rotation around the vertical axis of the boundary node distribution functions, and if all the values were multiplied by -1. This behaviour is represented in Figure 3.

The macroscopic density and velocities now correspond to the targeted values:

                 ρ = (-f * k ) -∆X/2 = ρ ũr + ∆t 2 ãr = -c kr (-f * k ) (-∆X/2)ρ = -u r + ∆t 2 a r ũz + ∆t 2 ãz = c kz (-f * k ) (-∆X/2)ρ = u z + ∆t 2 a z (20)
x-1 However, the formulas from Equation 19are not sufficient to get a good quality boundary conditions. Indeed, during the streaming process the domain

looses f * 3 + f * 6 + f * 7 from each boundary node and gets in return -f * 3 -f * 6 -f * 7 . There is a net mass loss ρ ′ r = 2(f * 3 +f * 6 +f * 7 )
. This quantity must be re-injected into the simulation domain to ensure the mass conservation.

The mass correction is performed in such a way that the post-streaming velocities and accelerations are not affected. Just after the streaming, the mass at the boundary node is ρr = f + k . This mass must be modified to (ρ + ρ ′ )r. We use:

f + k = f + k × 1 + ρ ′ r ρr (21)
The RW-Mirror boundary condition involves two steps: a streaming process (Equation 19) and a mass correction (Equation 21). It ensures the symmetry of the macroscopic quantities, and the mass conservation in the simulation domain.

Additional observations

Axis symmetric simulations are normally performed on a simulation domain that covers a half section [0, R] of the cylindrical domain, with proper boundary conditions near the symmetry axis. To observe the way the distribution functions are naturally arranged by the RW model around the axis location, we performed simulations on a whole section [-R, R] of the cylindrical domain. Moreover, from the full section simulations, we could observe the effect of the error singularity as we approach the axis at r = 0. The solution of the Poiseuille Hagen flow has a zero radial velocity: the singular term ∂ r (u r /r) in equation 14 is null. The simulation results were in good agreement with the analytical results, and no mass deviation was observed. However, the Side driven cavity flow implies radial velocity components: the singular term ∂ r (u r /r) in equation 14 is non null. We observed from the simulations a strong deviation of the total mass.

Symmetry axis

Wall

Periodic boundary condition

Thus, the mass correction of the RW-mirror BC compensates the natural mass deviation caused by the singular error term in the continuity equation.

Model implementation

Dellar [START_REF] Paul | An interpretation and derivation of the lattice boltzmann method using strang splitting[END_REF] interprets the 2 nd order LB scheme as a Strang Splitting procedure. Let S be the streaming operator, and C the collision operator. Dellar recommends the following sequence of operations for the computation of N t time-steps:

f k (n.∆t) = C 1/2 S (Nt-1)×C.S. CS . . . CS C 1/2 (f k (0)) (22)
The initial density, velocity, and force are respectively ρ ini , -→ u ini and -→ F 0ini .

At t = 0, the distribution functions are at the equilibrium state:

f k (0) = f eq k (ρ ini , -→ u ini ).
After the first half-collision, we get:

C 1/2 (f k (0)) = f eq k (ρ ini , -→ u ini ) + ∆t 2 R k ( -→ u ini , -→ F 0ini + -→ F ′ ini ) (23)
where -→ F ′ ini is calculated, using Equation 10, with the initial density and radial velocity.

Then, the streaming and collision operators are applied successively N t -1 times. The streaming operator S splits-up into two steps:

• Normal streaming inside the computational domain

• Correction at the boundary nodes

We used half-way bounce-back boundary conditions along the solid walls.

The collision operator C splits-up into three steps:

• Calculation of the macroscopic quantities from pre-collision distribution functions.

• Calculation of the equilibrium distribution functions (Equation 5), the virtual correcting force -→ F ′ (Equation 10), and the forcing terms R k (Equation 9)

• Collision (RHS of of Equation 4).

The determination of the macroscopic quantities ρ and -→ u is performed from pre-collision distribution functions f k . The density is directly deduced from the zero-th order moment:

ρ = k f k r (24)
For the velocity, it is slightly more complicated because the force term

-→ F = -→ F 0 +
-→ F ′ depends on the radial velocity u r . The first order moment is inverted to get the velocities in the r and z directions:

u r = 1 rρ + ∆t ρν r k c kr f k + ∆t 2 rF 0r + ρc 2 s u z = 1 rρ k c kz f k + ∆t 2 rF 0z (25)
At the end of the simulation, the available distribution functions are post collision distribution functions f * k . Their first order moment is different from Equation 7:

k c kα f * k = r ρu α + ∆t F α 2 (26)
Thus, the macroscopic velocities are obtained through:

u r = 1 rρ -∆t ρν r k c kr f * k - ∆t 2 rF 0r + ρc 2 s u z = 1 rρ k c kz f * k - ∆t 2 rF 0z (27)

Test Cases

The RW-LBM model is confronted to four test cases. The primary goal of those simulations is to compare the performances of the mirror boundary condition and the RW-mirror boundary conditions. The influence of the parameter τ on the model accuracy is investigated with the fourth test case.

The first test case is the Poiseuille Hagen flow. Then, we test the Wormersley flow because of its transitional behaviour. Those two test cases have analytical solutions with zero radial velocities. However, we emphasized the importance of the radial velocity on the error near the symmetry axis. So we investigate two more test cases: the Rayleigh Benard flow, and a made-up test case: the side driven cylinder.

Stopping criterion

The convergence error ε c is the error between two axial velocity fields separated by 500 time steps:

ε c = 1 N Z .N R i j |u z,ij (t) -u z,ij (t -500∆t)| 2 |u z,ij (t)| 2 (28)
Each simulation of a steady state system ended up with a convergence error below 10 -10 . For each of those simulations, the number of time steps we used is indicated through a non-dimension simulation time value κ. The definition of κ depends on the simulated system.

The initial density is: ρ 0 = 1. The total mass deviation ∆M is expressed in 29. It did not exceed 10 -15 during all the simulations. The mass conservation is ensured by the model.

∆M = i j ρ ij r i ∆X 2 R 2 H -1 (29)

Poiseuille Hagen Flow

The Poiseuille Hagen flow is a flow inside a vertical cylinder, driven by a homogeneous vertical force

-→ F 0 = -F 0z .
-→ e z . The simulated system is drawn in Figure 6. The symmetry axis is on the left. Half-way bounce-back boundary condition is applied on the right. Periodic boundary conditions are applied between the top and bottom limits of the domain. The analytical solution is ũr (r) = 0, and ũz (r) given by Equation 30, where U 0 is the velocity at the axis location (Figure 7).

Symmetry axis

Wall

Periodic boundary condition

ũz (r) = -U 0 1 - r 2 R 2 = - F 0z R 2 4ρν 1 - r 2 R 2 (30)
The dynamics of a Poiseuille Hagen flow is fully characterized by the Reynolds number Re = 2RU 0 /ν. We used the following settings: Re = 10, τ = 0.75, κ = (N t ∆t)U 0 /R = 100, R = 1 and ∆t/∆X 2 = 1. Density variations and horizontal velocities were close to the machine precision for all simulations.

The bottom axis of Figure 7 show the repartition of the axial velocity error along the radial coordinate, for the Mirror and the RW-Mirror boundary conditions. The error values are identical near the solid wall. They are due to the half way bounce-back scheme which, when used with a BGK model, causes a slight deviation of the wall location. This phenomenon is explained by Noble et al. [START_REF] Noble | A consistent hydrodynamic boundary condition for the lattice boltzmann method[END_REF] and Hecht and Harting [START_REF] Hecht | Implementation of on-site velocity boundary conditions for d3q19 lattice boltzmann simulations[END_REF]. The error repartition which is obtained using the Mirror-BC is much worse.

The error deviates strongly when approaching the axis.

To determine the degree of accuracy of the models, we compute the error ε on the axis axial velocity:

ε = U 0 -u k z0 U 0 (31)
The axial velocity u z0 = u z (r = 0) is extrapolated from the two neighbouring velocity values u z,0 and u z,1 , with a second order extrapolation scheme assuming ∂u z /∂r(0) = 0:

u z0 = 9 8 u z,0 - 1 8 u z,1 (32)
The error ε is represented as a function of the radial resolution N R in Figure 8. For both boundary conditions, we get ε ∝ ∆X 2 , the model is second order accurate. The axis velocity predictions are more precise with the RW-Mirror BC for all mesh sizes.

Womersley Flow

The Womersley flow is the flow inside a tube, driven by an oscillatory axial force:

-→ F 0 = -F 0z cos(ω.t) -→ e z (33)
The reference velocity is U 0 = F 0z R 2 /(4ρν). The dynamic of the systems depends on the Reynolds number Re = 2RU 0 /ν and the Womersley number α = R ω/ν. The analytical solution has a zero radial velocity, and an axial velocity given by Equation 34, where J 0 is the zero order complex Bessel function.

ũz (r, t) = abs

  4U 0 α 2 i   1 - J 0 iω/ν.r J 0 iω/ν.R   exp(iωt)   (34)
The parameters we used for the simulations are: Re = 1200, α = 20, τ = 0.6, R = 1 and ∆t/∆X 2 = 1. The simulations were run over 400 periods to reach the periodic regime.

Simulated axial velocity profiles are drawn for four moments of the oscillatory process in Figures 9 (mirror boundary conditions) and 10 (radius weighted mirror boundary conditions). In each figure, the simulation results are confronted to the theoretical profiles. In Figure 9, the mirror boundary causes a deviation of the velocity profile near r = 0 (circled area). This deviation is null with the Radius weighted boundary condition.

(N R = 32) ε = 1 N t T N tT k ũz (0, t k ) -u k z0 |ũ z (0, t k )| (35)
It is represented as a function of the radial resolution N R in Figure 11. It is difficult to extract any order accuracy from those curves. The precision is generally improved during the mesh refinement. Axis velocity predictions are on average 500 times more precise with the RW-mirror BC.

Rayleigh Benard convection

The choice of the Rayleigh Benard convection flow is motivated by the presence of non-zero radial velocity, and the existence of reference results in the literature.

The fluid, initially at rest, is enclosed inside a vertical cylinder, with a hot temperature T h at the bottom, a cold temperature T c at the top, and adiabatic conditions on the surrounding wall (Figure 12). The Boussinesq force

-→ F 0 = ρgβ(T -T m )
-→ e z drags a toroidal flow inside the cylinder. The spatial derivatives ∂r , ∂z , ∂2 r and ∂2 z are centred discrete derivatives. They are second order accurate in ∆X.

Symmetry axis

T k+1 = T k + ∆t -u r ∂r T k -u z ∂z T k + α ∂2 r T k + ∂2 z T k + ∂r T k r (36)
The reference velocity is U 0 = gβ(T h -T c )H. The following criteria were fulfilled in order to ensure the stability of the energy conservation model (see Tannehill et al. [START_REF] Tannehill | Computational Fluid Mechanics and Heat Transfer[END_REF]):

• Diffusive criterion: α∆t ∆X 2 ≤ 1 4 • Two dimensions Courant Friedrichs Lewi criterion: U0∆t ∆X ≤ 1 2
• Mixed criterion:

U 2 0 ∆t α ≪ 1
The dynamics of the system is fully characterized by three dimensionless numbers: Rayleigh Ra = gβ(T h -Tc)H 3 αν , Prandtl P r = ν α , and shape factor Γ = R/H. The parameters we used were: Ra = 5000, P r = 0.7, Γ = 1,

τ = 0.55, R = 1, ∆t/∆X 2 = 1 and κ = αN t ∆t/R 2 = 14.
The reference result given by Lemembre and Petit [START_REF] Lemembre | Laminar natural convection in a laterally heated and upper cooled vertical cylindrical enclosure[END_REF] is the maximum dimensionless axis velocity: u z0,max /U 0 = 0.351. The small number of digits of this reference result is due to the weakly oscillatory nature of the toroidal flow. The simulation results are compared to the reference value in Figure 13.

The models with both boundary conditions converge toward the right velocity value. However, the precision is far better with the RW-mirror BC than with the mirror BC. We fix z = 0 at the bottom of the system. A weighting is applied on the side wall velocity to avoid discontinuity issues near the corners.

U ′ 0 (z) = U 0 × 1 - 2z H -1 18 2 (37)
The dynamic of the system is entirely defined by two dimensionless numbers:

Reynolds number Re = U 0 H/ν, and a shape factor: Γ = H/R. We performed the simulations for Re = 50 and Γ = 2, with the following settings: τ = 0.6, R = 1, ∆t/∆X 2 = 1, and κ = N t ∆tU 0 /H = 30.

We got the reference velocity values from much finer simulations performed on a finite volume model based on the axis-symmetric Navier-Stokes equations.

For those simulations, we used a steady state model, with a pressure-velocity coupling scheme which offers a second order accuracy for the pressure, and a third order accuracy for the momentum, through the MUSCLE scheme. The reference results (ũ r , ũz ) were obtained with a regular square mesh of 1024×2048 cells. A fourth order spline interpolation technique was used to move reference and LBM velocity values to the same locations for comparisons.

The lines z = 0.5, z = 1 and z = 1.5 are identified in Figure 15. The radial and axial velocity profiles along those lines are represented in Figure 16, where the reference results are drawn in black. Absolute error repartitions are drawn in the bottom axis.

The error values are much greater for the side driven cavity flow (around 10-15% with N R = 64) than for the Poiseuille Hagen flow (less than 0.1%, see fig. 7). The present test case is stiffer, due to non-null radial velocity, and velocity gradients in both radial and axial directions.

The distinction between Mirror and RW-mirror BC is perceptible near the axis (circled areas). With the RW-Mirror BC, the radial velocity error tends toward zero and the axial velocity error stabilizes. With the Mirror BC, there is a distortion of the error profiles near the axis. The corresponding error are reported on the curves below.

The global simulation error ε is computed from the axial velocity profiles at the axis location r = 0 (Equation 38), where u z0,j are extrapolated values (Equation 32).

ε = 1 N Z j |u z0,j -ũz0,j | U 0 (38)
The error ε is drawn in figure 17 for several mesh sizes. The convergence is now first order. This is an issue, because the truncation error analysis does not highlight any first order error term. It seems that the spurious terms u r /r become so great near the axis that they influence the convergence rate of the model. The effect of the relaxation parameter τ is investigated in Figure 18. We set R = 1, U 0 = 1, N R = 32 and ν = 2/Re, with Re = 5, Re = 10 or Re = 100. Since ν is fixed and τ is varying, the ratio ∆t/∆X 2 is varying from one simulation to another.

For high Reynolds numbers (Re = 50, Re = 100), the trend is clear: the closer is the relaxation parameter to 0.5, the more precise is the model. Here, the value ν = 1/Re being very small, the dominant error terms in Equations 14 and A.37 are multiples of (2τ -1). The error is minimized when τ approaches 0.5. For Re = 5, the minimal error occurs near τ = 6.5. The dominant error terms in 14 and A.37 are no more the multiples of (2τ -1), because the terms in 6τ 2 -6τ + 1 (root τ ≈ 0.789) and 12τ 2 -12τ + 1 (root τ ≈ 0.908) are getting bigger. The optimal value of τ is then shifted upward. This observation is in concordance with the analysis of the classical BGK model made by Holdych [START_REF] David | Truncation error analysis of lattice boltzmann methods[END_REF].

Conclusions

The present paper offers a deep and critical analysis of the RW-LB model, based the expression of its generalized equation. The generalized equation assesses the consistency of the model. All the error terms are multiples of ∆X 2

or more: the model is second order accurate. However, we identified spurious error terms u r /r which may considerably reduce the model precision when u r is not null.

The fact that the spurious terms might reach extreme values near the axis motivated an investigation about the axis boundary condition. This investigation resulted in the definition of a new boundary condition which mimics the natural behaviour of the model when u r = 0, and compensates the mass loss

f k (1) -f k = - 1 τ (f k -f eq k ) + ∆t 1 - 1 2τ R k (A.2)
By simply shifting the time continuum by -∆t, we get:

f k = 1 - 1 τ f k (-1) + 1 τ f eq k (-1) + ∆t 1 - 1 2τ R k (-1) (A.3)
This expression is applied recursively over an infinite number of time steps.

By this process, all the unknown distribution functions f k (-n) are replaced by equilibrium distribution functions f eq k (-n) and source terms R k (-n):

f k = 1 τ ∞ n=1 1 - 1 τ n-1 f eq k (-n)+ ∆t 1 - 1 2τ ∞ n=1 1 - 1 τ n-1 R k (-n) (A.4)
f eq k and R k are directly related to the macroscopic quantities ρ, p, -→ u , and -→ F . Equation A.4 expresses f k ( -→ x , t) as a function of the macroscopic quantities at the previous time steps. We use Taylor series expansion to express them at the present time:

f eq k (-n) = f eq k + ∞ m=1 1 m! -n∆t∂ t -n∆X -→ w k . -→ ∇ m f eq k R k (-n) = R k + ∞ m=1 1 m! -n∆t∂ t -n∆X -→ w k . -→ ∇ m R k (A.5)
Those expressions are introduced in A.4:

f k = 1 τ ∞ n=1 1 - 1 τ n-1 f eq k + ∆t 1 - 1 2τ ∞ n=1 ∞ m=1 1 - 1 τ n-1 R k + 1 τ ∞ n=1 ∞ m=1 1 - 1 τ n-1 (-n) m m! ∆t∂ t + ∆X -→ w k . -→ ∇ m f eq k + ∆t 1 - 1 2τ ∞ n=1 1 - 1 τ n-1 (-n) m m! ∆t∂ t + ∆X -→ w k . -→ ∇ m R k (A.6)
We introduce the function p[τ, m] which is defined as:

p[τ, m] = 1 τ 2 ∞ n=1 1 - 1 τ n-1 (-n) m (A.7)
The function p[τ, m] has a formal expression for the first values of m:

p[τ, 0] = 1/τ p[τ, 1] = -1 p[τ, 2] = 2τ -1 p[τ, 3] = -6τ 2 + 6τ -1 p[τ, 4] = 24τ 3 -36τ 2 + 14τ -1 (A.8)
Equation A.6 is then reduced to the following expression:

f k =f eq k + ∆t τ - 1 2 R k + τ ∞ m=1 p[τ, m] m! ∆t∂ t + ∆X -→ w k . -→ ∇ m f eq k + ∆tτ (2τ -1) 2 ∞ m=1 p[τ, m] m! ∆t∂ t + ∆X -→ w k . -→ ∇ m R k (A.9)
Finally, we expand the global differential operators:

∆t∂ t + ∆X -→ w k . -→ ∇ m = m q=0 m! q!(m -q)! ∆t m-q ∆X q [∂ t ] m-q -→ w k . -→ ∇ q (A.10)
The previous equation becomes:

f k = f eq k + ∆t (2τ -1) 2 R k + ∆tτ ∞ m=1 m q=0 p[τ, m] q!(m -q)! ∆t m-q-1 ∆X q [∂ t ] m-q -→ w k . -→ ∇ q f eq k + ∆tτ ∞ m=1 m q=0 p[τ, m] q!(m -q)! 2τ -1 2 ∆t m-q ∆X q [∂ t ] m-q -→ w k . -→ ∇ q R k (A.11)
In Equation A.11, the distribution functions f k ( -→ x , t) are entirely defined from the macroscopic quantities at the time t and the coordinate -→ x . We kept one term τ ∆t at the beginning of each sum to simplify subsequent calculations.

Appendix A.2. Explicit expressions of the equilibrium distribution functions and force terms

To select the dominant terms in Equation A.11, we must expand the expression of f eq k and R k as functions of ∆t and ∆X. For that purpose, we consider c 2 s = ∆X 2 /3∆t 2 and -→ c k = (∆X/∆t) -→ w k .

The first element of f eq k is directly expressed as a function of the pressure p through the density / pressure relationship for incompressible flows:

ρ = p c 2 s = 3∆t 2 ∆X 2 p (A.12)
With this artefact, the pressure p appears naturally in the momentum equation. We expand f eq k and R k :

f eq k = ω k r p c 2 s + ρ -→ c k . -→ u c 2 s + ρ ( -→ c k . -→ u ) 2 2c 4 s -ρ -→ u 2 2c 2 s = ω k r 3∆t 2 ∆X 2 p + 3∆t ∆X ρ -→ w k . -→ u + 9∆t 2 2∆X 2 ρ( -→ w k . -→ u ) 2 - 3∆t 2 2∆X 2 ρ -→ u 2 = ∆t 2 ∆X 2 3ω k rp + 9ω k 2 ρr( -→ w k . -→ u ) 2 - 3ω k 2 ρr -→ u 2 + ∆t ∆X [3ω k ρr -→ w k . -→ u ] = ∆t 2 ∆X 2 [A f k ] + ∆t ∆X [B f k ] (A.13) R k = ω k r -→ c k . -→ F c 2 s + ( -→ c k . -→ u )( -→ c k . -→ F ) c 4 s - -→ u . -→ F c 2 s = ω k r 3∆t ∆X -→ w k . -→ F + 9∆t 2 ∆X 2 ( -→ w k . -→ u )( -→ w k . -→ F ) - 3∆t 2 ∆X 2 -→ u . -→ F = ∆t 2 ∆X 2 9ω k r( -→ w k . -→ u )( -→ w k . -→ F ) -3ω k r -→ u . -→ F + ∆t ∆X 3ω k r -→ w k . -→ F = ∆t 2 ∆X 2 [A Rk ] + ∆t ∆X [B Rk ] (A.14)
In both expanded expressions, we get some that are multiples of ∆t 2 /∆X 2 , and on term which is a multiple of ∆t/∆X. The ∆t 2 /∆X 2 elements contain an even number of -→ w k . The ∆t/∆X element contain and odd number of -→ w k . This distinction is important given that the odd-order momentums in -→ w k are null. For the even order moments: We get from A.11:

k ω k = 1 k ω k w kα w kβ =    1/3 if 2r or 2z 0 if rz k ω k w kα w kβ w kγ w kτ =          1/3 if
f k = f eq k + ∆t (2τ -1) 2 R k + ∆tτ ∞ m=1 m q=0 p[τ, m] q!(m -q)! ∆t m-q+1 ∆X q-2 [∂ t ] m-q -→ w k . -→ ∇ q A f k +∆tτ ∞ m=1 m q=0 p[τ, m] q!(m -q)! 2τ -1 2 ∆t m-q+2 ∆X q-2 [∂ t ] m-q -→ w k . -→ ∇ q A Rk + ∆tτ ∞ m=1 m q=0 p[τ, m] q!(m -q)! ∆t m-q ∆X q-1 [∂ t ] m-q -→ w k . -→ ∇ q B f k +∆tτ ∞ m=1 m q=0 p[τ, m] q!(m -q)! 2τ -1 2 ∆t m-q+1 ∆X q-1 [∂ t ] m-q -→ w k . -→ ∇ q B Rk (A.16) Appendix A.3. Scaling hypothesis (A.17)
The scaling strategy is justified in section 4.1. The space step ∆X is the small parameter that is reduced when refining the mesh. The time step is adjusted through diffusive scaling: ∆t ∝ ∆X 2 . The ratio (2τ -1)/6ν is fixed.

Thus, we get from Equation 11:

∆t = 2τ -1 6ν ∆X 2 (A.18)
This expression is introduced in Equation A.16

f k = f eq k + ∆t (2τ -1) 2 R k + ∆tτ ∞ m=1 m q=0 p[τ, m] q!(m -q)! 2τ -1 6ν m-q+1 ∆X 2m-q [∂ t ] m-q -→ w k . -→ ∇ q A f k +∆tτ ∞ m=1 m q=0 p[τ, m] q!(m -q)! 3ν 2τ -1 6ν m-q+3 ∆X 2m-q+2 [∂ t ] m-q -→ w k . -→ ∇ q A Rk + ∆tτ ∞ m=1 m q=0 p[τ, m] q!(m -q)! 2τ -1 6ν m-q ∆X 2m-q-1 [∂ t ] m-q -→ w k -→ ∇ q B f k +∆tτ ∞ m=1 m q=0 p[τ, m] q!(m -q)! 3ν 2τ -1 6ν m-q+2 ∆X 2m-q+1 [∂ t ] m-q -→ w k . -→ ∇ q B Rk (A.19)
Now, there is only one small parameter (∆X) in the native equation. The order of magnitude of the different terms will be determined from the power of this parameter.

Appendix A.4. Zero order moment: continuity equation

The continuity equation is obtained by computing the zero order momentum of Equation A. [START_REF] Guo | Axisymmetric lattice boltzmann method[END_REF]. We know that k f k = k f eq k = ρr and k R k = 0: all the first line terms are eliminated. The coefficient τ ∆t at the begining of the remaining terms is removed:

0 = ∞ m=1 m q=0 p[τ, m] q!(m -q)! 2τ -1 6ν m-q+1 ∆X 2m-q [∂ t ] m-q k -→ w k . -→ ∇ q A f k ∞ m=1 m q=0 p[τ, m] q!(m -q)! 3ν 2τ -1 6ν m-q+3 ∆X 2m-q+2 [∂ t ] m-q k -→ w k . -→ ∇ q A Rk ∞ m=1 m q=0 p[τ, m] q!(m -q)! 2τ -1 6ν m-q ∆X 2m-q-1 [∂ t ] m-q k -→ w k . -→ ∇ q B f k ∞ m=1 m q=0 p[τ, m] q!(m -q)! 3ν 2τ -1 6ν m-q+2 ∆X 2m-q+1 [∂ t ] m-q k -→ w k . -→ ∇ q B Rk (A.20) ∆X 0 ∆X 2 line 1 (A f k ) - (1, 0), (2, 2) line 2 (A Rk ) - - line 3 (B f k ) (1, 1) (2, 1), (3, 3) line 4 (B Rk ) - (1, 1)
Table A.1: Selected combination of (m, q) for the continuity equation

Here, despite the infinite number of terms in the sum, only a few of them have an order of ∆X lower or equal to 2. Moreover, for the moments to be non zero, the coefficient q must be even in the two first lines, and odd in the two last lines.

The combinations (m, q) that lead to a ∆X order lower or equal to 2, are given in Table A.1 for each line of the former equation. After computing the related moments, we get: Table A.2: Selected combination of (m, q) for the momentum equation

0 = -∂ α (ρru α ) + (2τ -1) 2 12ν ∆X 2 - r ν ∂ t p + 2∂ tα (ρru α ) + ∂ αβ (rpδ αβ + ρru α u β ) -∂ α (rF α ) + -6τ 2 + 6τ -1 6 ∆X 2 [∂ ββα (ρru α )] + O(∆X
For the other terms, we use the following development:

k c kα [-] = ∆X ∆t k w kα [-] = 2τ -1 6ν -1 ∆X -1 k w kα [-] (A.23)
The only surviving term in the first line of the first order momentum equation is ∆tτ (rF α ). The coefficient τ ∆t at the beginning of each line is eliminated.

We get:

0 = rF α + ∞ m=1 m q=0 p[τ, m] q!(m -q)! 2τ -1 6ν m-q ∆X 2m-q-1 [∂ t ] m-q k -→ w k . -→ ∇ q w kα A f k ∞ m=1 m q=0 p[τ, m] q!(m -q)! 3ν 2τ -1 6ν m-q+2 ∆X 2m-q+1 [∂ t ] m-q k -→ w k . -→ ∇ q w kα A Rk ∞ m=1 m q=0 p[τ, m] q!(m -q)! 2τ -1 6ν m-q-1 ∆X 2m-q-2 [∂ t ] m-q k -→ w k . -→ ∇ q w kα B f k ∞ m=1 m q=0 p[τ, m] q!(m -q)! 3ν 2τ -1 6ν m-q+1 ∆X 2m-q [∂ t ] m-q k -→ w k . -→ ∇ q w kα B Rk (A.24)
The coefficient q must be odd to get non zero moments in the two first lines.

It must be even for lines 3 and 4. The combinations (m, p) for an order of ∆X equal or lower to 2 are given in Thus, when extracting the density from the zero order derivative ∂ α (ρru α ), we must take into account the corresponding second order element in the error terms. However, when we extract the density from derivatives that are already within the second order terms, the equation is not modified (the additional terms would be fourth order). We get the modified continuity equation: We recover on the zero order terms the Navier Stokes momentum after sub-stituting the force term F α by its expression in A.27.

∂ α (ru α ) = (2τ -

  14 and 15 are identical to the cylindrical Navier Stokes Equations 1 and 2. The consistency of the model is assessed wherever the truncation error do not reach infinite values. One concern that might arise with the RW-LB model is the repartition of the error along the radial coordinate, due to the presence of r in the definition of the moments. The expression of the error terms can give us some answers about this concern. There are two kinds of error terms in Equations 14 and A.37. The first ones are functions of rp or r -→ u . The other ones are functions of u r /r. The first kind of error terms are functions of quantities that are proportional to r. However, Equations 14 and 15 govern the evolution of rρ and rρ -→ u , which are also proportional to r. The r-proportionality in the error and the main equation terms compensate each other. The first kind of error terms should not involve any increasing trend along the radius coordinate.

Figure 1 :

 1 Figure 1: Boundary node and ghost node near the symmetry axis

Figure 2 :

 2 Figure 2: Target behaviour for the axis boundary condition

Figure 3 :

 3 Figure 3: Ghost node distribution function for the mirror boundary condition

Figure 4 :

 4 Figure 4: Ghost node distribution function for the RW-mirror boundary condition

Figure 5 :

 5 Figure 5: Simulation domain for the preliminary simulations of the Poiseuille flow

Figure 6 :

 6 Figure 6: Simulation domain of the Poiseuille Hagen flow

Figure 7 :

 7 Figure 7: Poiseuille Hagen flow velocity profile for N R = 64 and Re = 10. Bottom axis: distribution of the axial velocity error along the radial coordinate

Figure 8 :

 8 Figure 8: Grid dependency on the axis velocity for the Poiseuille flow

Figure 9 : 22 -Figure 10 :

 92210 Figure 9: Typical solution for the wormersley flow with α = 20. Mirror boundary conditions (N R = 32)

Figure 11 :

 11 Figure 11: Error on the axis velocity for the Wormersley flow

Figure 12 :

 12 Figure 12: Simulated system for the Rayleigh Benard convection

Figure 13 :

 13 Figure 13: Maximum axis velocity for Rayleigh Benard convection

Figure 14 :Figure 15 :

 1415 Figure 14: Velocity streamlines for Rayleigh-Benard convection:14(a) initial cold temperature, 14(b) initial hot temperature

Figure 16 :

 16 Figure 16: Axial and radial velocity profiles along the lines y = 0.5, y = 1 and y = 1.5, for N R = 64 and Re = 50. Reference profiles are drawn in black, simulated profiles with the RW-mirror and mirror boundary conditions are respectively in solid lines and dashed lines.

Figure 17 :

 17 Figure 17: Axis axial velocity error for the side driven cylinder flow

Figure 18 :

 18 Figure 18: Effect of the relaxation parameter τ on the model precision.

4 )(A. 21 )c kα f k = ρru α - ∆t 2 rF α k c kα f eq k = ρru α k c kα R k = rF α (

 421α Appendix A.[START_REF] Chai | A novel lattice boltzmann model for the poisson equation[END_REF]. First order moment: momentum equationThe momentum equation is obtained from the first order moments of Equation A.19. For the terms in the first line, we use:

  Table A.2. The momentum equation at this ββα(ρru α ) + 2∂ ααβ (ρru β ) -∂ rα (ρu α ) + δ αr ∂ α -∂ β (ρru β )It is clear that the divergence terms ∂ α (ρru α ) have an order of magnitude ∆X 2 . They are removed from the second order elements:

	0 =(eq.A.25) +	(2τ -1) 2 12ν	∆X 2 ∂ α (eq.A.31)
	0 = -∂ α (ρru α ) + (2τ -1) 2 12ν	∆X 2 -	r ν	∂ t (p) + ∂ tα (ρru α )
	+	(2τ -1) 2 12ν	∆X 2 ν ∂ r	+	ρu r r
	+	-6τ 2 + 6τ -1 6	∆X 2 [∂ ββα (ρru α )] + O(∆X 4 )
								(A.32)
	∂ α (ρru α ) =	(2τ -1) 2 12ν	∆X 2 -	r ν	∂ t (p) + ν∂ r	ρu r r	(A.33)
	The density variations are second order terms:
	∂ρ =	∂p c 2 s	=	(2τ -1) 2 12ν	∆X 2 ∂p ν	(A.34)
	stage turns into:

  For the momentum equation, we begin the simplifications by substituting the divergence terms using Equation A.33:0 = -∂ t (ρru α ) -∂ β (rpδ αβ + ρru α u β ) + ν ∂ ββ (ρru α ) -∂ r (ρu α ) + δ αr [∂ tt (ρru α ) -∂ t (rF α ) + 2∂ tβ (rpδ αβ + ρru α u β ) -∂ β (ru α F 0β + ru β F 0α )] ∂ tββ (ρru α ) + 2∂ tαβ (ρru β ) + ∂ αββ (rp -[∂ tt (ru α ) + 2∂ tβ (ru α u β )] (rF α ) -∂ β (ru α F 0β + ru β F 0α )] ∂ tββ (ru α ) + 2∂ tαβ (ru β ) + ∂ αββ (

									ρu r r	+ rF α
	+ ∆X 2 + (2τ -1) 2 12ν (2τ -1) 2 12ν ∆X 2 [δ αr ∂ t (p) -2∂ tα (rp)]
	+	(2τ -1) 2 12ν	∆X 2 ν ∂ ββ (rF α ) + 2∂ αβ (rF β ) + 2∂ r	ρu α u r r	+ 2∂ β	ρu β u r r	δ αr
	+	(2τ -1) 2 12ν	∆X 2 ν 2 2∂ αr	ρu r r	-	δ αr r	∂ r	ρu r r
	+	-6τ 2 + 6τ -1 6	∆X 2 ρru γγ 2	) +	⊕ α
	+	12τ 2 -12τ + 1 4	ν∆X 2	⊗ u r r	+	rF α ρ
	+ [-ru + (2τ -1) 2 12ν ∆X 2 ρν (2τ -1) 2 12ν ∆X 2 + (2τ -1) 2 12ν ∆X 2 ρ [δ αr ∂ + (2τ -1) 2 12ν ∆X 2 ρ [-∂ t + (2τ -1) 2 12ν ∆X 2 ν ∂ ββ (rF α ) + 2∂ αβ (rF β ) ρ	+ 2∂ r	u α u r r	+ 2∂ β	u β u r r	δ αr
	+ +	1) 2 12ν ∆X 2 ν 2 2∂ αr ∆X 2 -u r r -6τ 2 + 6τ -1 (2τ -1) 2 12ν 6 ∆X 2 rp r ρν (∂ t p + u α ∂ α p) + ν∂ r -δ αr r u r ∂ r r u r r ρ	-	ru γγ 2	(A.35) ) + ⊕ α /ρ
	+	12τ 2 -12τ + 1 4	ν∆X 2	⊗ α /ρ + O(∆X 4 )
									(A.37)

α + O(∆X 4 ) (A.36)

The density is taken out of the derivatives:

0 = -∂ t (ru α ) -∂ β (ru α u β ) -∂ r (rp) ρ + ν ∂ ββ (ru α ) -∂ r (u α ) + δ αr α ∂ t (p) -ru α u β ∂ β (p)] t (p) + 2∂ β (ru α )∂ β (p) + ru α ∂ ββ (p) -u α ∂ r (p)]

Acknowledgments

We would like to thank the French agency ANR for its financial support through the project Stock-E MICMCP.

Appendix A. Taylor series analysis

when u r = 0. The new boundary condition clearly provides more accurate simulation results.

The spurious error terms have no effect on the simulations results when u r = 0. We recovered the second order accuracy that we expected from the generalized equation. However, our last simulation results showed that strong values of u r reduced the order of accuracy of the model. We still don't know how to quantify this reduction, this is a work to be continued.

Anyway, despite those limitations, the RW-LB model still has the great advantage over any other model that the collision process is purely local. It remains an alternative to any Navier Stokes discretization when an exceptional precision is not necessarily needed.

Appendix A.1. Native Expression

The first step of the Taylor series analysis is to build the native expression of the model. The native expression is when the distribution functions f k ( -→ x , t) are expressed as a function of macroscopic quantities ρ, p, -→ u , and -→ F at the time t and coordinate -→ x .

To distinguish the time step ∆t and space step ∆X, the particle velocities -→ c k are formulated as the product of an amplitude ∆X/∆t and a direction vector -→ w k . For the sake of clarity, we adopt the following formulation:

We start with the LBM Equation 4:

Where:

. Extraction of the additional momentum terms

Here come the tricky part of the development. In Equation A.25, only the first line seems to be zero order, the following lines being multiples of ∆X 2 .

However, due to the radius coordinate appearing within space derivatives, and due to the definition of the force term, one additional zero order term arises.

The definition of the force term is:

(line 2):

The cross terms u α p and u β p are of order -2. Indeed:

They are extracted from the second order elements of Equation A.25, providing additional zero order elements in the momentum equation:

The momentum equation then turns into:

First, we use the zero order terms of Equation A.31 to simplify the continuity equation:

In section Appendix A.6, the additional zero order element was extracted from the force term F ′ α . We shall point out that, even if the definition of F ′ α did not contain the pressure p, the same additional zero order element would appear from the double time derivative ∂ tt (ρru α ) in Equation A. [START_REF] Guo | Discrete lattice effects on the forcing term in the lattice boltzmann method[END_REF].

The procedure for extracting this term is not so simple. First, we express the zero order part of Equation A. [START_REF] Guo | Discrete lattice effects on the forcing term in the lattice boltzmann method[END_REF], by taking into account the additional zero order term through the unknown variable K α :

Then, Equation A.38 is used to substitute the inner time derivative of ∂ t (∂ t (ρru α )). One of the resulting terms is ∂ tβ (ρru α u β ). We develop this term in the following manner:


J Jumel

F Lepoutre

F Enguehard

D Rochais

M Cataldi

Microscopic thermal characterization of C/C composites

Introduction

Both carbon fibers and matrix of C/C composites (figure 1) are made of graphite nanoelements, whose texture rules the local thermal properties, and consequently the macroscopic ones of the composite [START_REF] Legendre | Le matériau carbone[END_REF]. Thus, precise characterizations of fibers and matrix are required to improve these composites and to optimize their processes. It has been demonstrated that photoreflectance microscopy is able to provide the thermal properties measurement of complex materials such as graphite at a microscopic scale [START_REF] Pottier | [END_REF]3]. Nevertheless, additional efforts are necessary to remove the measurement artifacts, encountered when characterizing C/C composites local properties [4].

Microscopic thermal characterization by photoreflectance microscopy

Photothermal microscopy is based on the measurement and analysis of the periodic temperature increase induced by the absorption of an intensity modulated laser beam (pump beam). By detecting with a secondary laser beam (probe beam), the thermally induced reflection coefficient variations, we measure the temperature increase at the sample surface with a sensitivity better than 10 -3 K.Hz -1/2 . This optical method has a micrometric spatial resolution, and furthermore, is non destructive and contactless.

Thermal waves propagation

When an isotropic, homogeneous medium is heated by a periodic point-like heat source, the periodic temperature increase, also called "thermal wave", at a distance r from the pump location is described by the following equation:

( ) ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + - = r j kr Q x T µ π δ 1 exp 4 (1)
where:

f π α µ= (2)
is called thermal diffusion length, f being the modulation frequency of the harmonic heat source, α=k/ρC the thermal diffusivity of the material, k its thermal conductivity, and ρC its volume specific heat. Thermal waves are damped, the higher the modulation frequency is, the sharper the periodic temperature decrease. At a few hundreds of kHz and for usual thermal diffusivity values (smaller than 1cm 2 /s), the thermal wave is confined in a few cubic micrometers volume. The argument of δT (i.e. its phase) slope along r is equal to µ 1 -, which permits to evaluate the thermal diffusivity of the media, using a simple linear regression. A focused laser beam, is generally used as pump beam. Its energy distribution, which limits the spatial resolution of the measurements [4], must be taken into account.

Photoreflectance measurements principle and experimental set-up

To detect the temperature variations with a microscopic spatial resolution, the most convenient technique is photoreflectance microscopy. It consists in measuring with a photocell the light flux of a continuous laser beam (probe beam) of intensity I 0 , after being reflected by the sample surface in the periodically heated zone. If the reflection coefficient of the material, R, varies linearly with temperature,

] ) ( 1 1 [ ) ( 0 0 t j e r T T R R R r R ω δ ∂ ∂ + = (3)
the periodic reflected probe beam intensity variations are then equal to:

] ) ( 1 .[ . ) ( 0 0 0 t j e r T T R R I R r ω δ φ ∂ ∂ = (4)
These very small periodic intensity variations can be detected by a photocell providing the signal is filtered by a lock-in analyzer which strongly increases the signal to noise ratio. The signal being proportional to the local periodic temperature increase, the phase of the signal only depends on the material thermal diffusivity, while, on the contrary, the amplitude depends on the optical properties as well. The experimental set-up, described on figure (2), contains three main parts. An optical set-up to focus and to position the pump and probe beams, a device to measure the reflected probe beam flux, and several electronic devices to detect the signal and drive the experiment. The pump beam is an Ar + laser whose maximum power is 600 mW. The intensity is modulated with an acousto-optic modulator operating up to 1 MHz. An a.c. waves generator delivers the driving signal. The pump beam is then orientated by a dichroïc mirror and finally focused onto the sample surface with a microscope. The probe beam is a laser diode that goes through a quarter wave plate, the dichroïc mirror and then is focused onto the sample surface with the same microscope. After reflection, it goes through the quarter wave plate again and is then deviated to the photodiode with a beam splitter cube. An optical filter prevents from any return of the pump beam toward the detector. A lock-in amplifier extracts the amplitude and phase of the periodic signal. A personal computer controls the dichroic mirror orientation and consequently the distance r between probe and pump beam positions.

C/C Composites microscopic thermal characterization

In this section, we wish to emphasize the influence of graphite local texture on the microscopic thermal properties of carbon fibres and matrix. At first, heat diffusion in the transverse section of different types of carbon fiber is visualized. In a second step, the three dimensional thermal diffusivity tensor of various morphologies of pyrocarbon matrix are precisely measured.

Heat diffusion in carbon fiber visualization.

As mentioned in section 2, the phase signal evolution around the heat source reveals the way heat diffuses in the media. Then, to point-out the role of fiber morphology in the heat transfer process at a microscopic scale, the heat diffusion in fiber cross section is visualized by focusing the pump beam on the fiber and by measuring the periodic temperature increase around this point. Figure 3. shows typical results in the case of pitch based and PAN based fibers. When the heat source is slightly offset with respect to the pitch based fiber cross section center, the strong thermal anisotropy of heat diffusion immediately appears. The thermal wave essentially propagates along the radial direction, while it is strongly blocked along the circumferential one. These phenomenon are due to the strong radial texture of this fiber as revealed by the fractography presented in figure 4. The same measurement performed in a PAN based fiber confirm the expected transversely isotropic thermal behavior of these fibers due to their crumpled paper like microstructure.

Pyrocarbon matrix characterization.

The pyrocarbon matrix is made of graphitic planes that surround the carbon fiber such as in a leek. Because of the strong thermal anisotropy of the graphite this makes the thermal properties of the matrix orthotropic in cylindrical coordinates. The curvature of the graphitic planes has to be taken into account for precise determination of the matrix circumferential thermal diffusivity. The periodic temperature increase, in the matrix, during the photoreflectance experiment is approximated by: ( )

⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + - = R j R k R T r r 1 exp 4 1 ~µ π δ (8) where ( ) 2 0 2 2 cos 2 ~z k k k k rR R r R z r r + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - - + = θ θ θ (9)
k z , k r and k θ are respectively the thermal conductivities of the media in the longitudinal, radial and circumferential direction. (R,θ 0 ) being the polar coordinates of the heat source with respect to the fibre referential. According to relations (8) and (9), when the graphitic curvature is small, the phase slope along the circumferential direction permits to determine the matrix circumferential thermal diffusivity, as explained in section 2. The results given in table 1 demonstrate, at a microscopic scale, the correlation between graphitization and thermal ansitropy. Furthermore, it gives an experimental support to supply quantitative microscopic thermal properties. Correlated with TEM or polarized light microscopies, these results help material scientists to determine the optimal process parameters and permit physical interpretation of the macroscopic properties evolution of the composite.

Conclusion

Once having recalled the photoreflectance microscopy principle and performance, experimental results showing heat diffusion, at a microscopic scale, in PAN and pitch based carbon fibers were presented, in order to outline the role of fiber morphology during the microscopic heat transfer processes. Using a rigorous thermal parameters extraction scheme, different varieties of pyrocarbon matrix are then analyzed, in order to determine quantitatively their thermal diffusivity tensor. Finally, the results on different samples are analyzed in relation with their known morphologies as determined by polarized light microscopy.

Figure 1 .

 1 Figure 1. Micrograph of a C/C composite showing the locations and the dimensions of carbon fibers and matrix.

Figure 2 .

 2 Figure 2. Description of the photothermal microscope.

Figure 3 .

 3 Figure 3. Visualization of heat diffusion in a carbon fiber cross section. On the left, pitch based fiber, on the right, PAN based fiber, grayscales are in degrees.

Figure 4 .

 4 Figure 4. On the left, pitch based fiber fractography revealing the strong radial morphology of this fiber (typical diameter 10 µm), and, on the right, typical high modulus PAN based carbon fiber "transversely isotropic" microstructure.

Figure 3 .

 3 Figure 3. Experimental configuration for matrix characterization.

Figure 4 .

 4 Figure 4. Photoreflectance measurement in a matrix longitudinal section (f = 1MHz), attenuation gray scale in dB (right), Phase gray scale in degree (left).

Figure 5 .

 5 Figure 5. Photoreflectance measurement in a matrix cross section (f = 1MHz).

Table 1 .

 1 Figures6 and 7show typical aspects of the periodic temperature increase field in matrix Sample F: Sample B, + medium temperature heat treatment. Influence of processes on the matrix thermal properties.

	Sample	A	B	C	D	E	F
	α θ (cm 2 /s) α θ /α r	0,17 6	0,14 7	0,42 20	1,25 22	3,04 48	0,48 17

longitudinal and cross sections. From these experimental data, thermal diffusivities along the longitudinal, radial and circumferential directions are found to be equal respectively to: 8.1x10 -5 m 2 /s, 3.0x10 -6 m 2 /s and 7.8x10 -5 m 2 /s. Figure 7. illustrates the effect of graphitic planes curvature with heat diffusing circularly. Using this method, we studied six different samples of various morphologies determined by polarized light microscopy: Sample A: dark laminar, no heat treatment, Sample B: smooth laminar, no heat treatment, Sample C: rough laminar, no heat treatment, Sample D: Sample C, + medium temperature heat treatment, Sample E: Sample C, + high temperature heat treatment,