Franck Lominé
email: franck.lomine@univ-nantes.fr

Luc Scholtès

Luc Sibille

Philippe Poullain

Modeling of fluid-solid interaction in granular media with coupled lattice Boltzmann/discrete element methods: application to piping erosion

Keywords: particle-fluid system, lattice Boltzmann method, discrete element method, fluid flow, soil piping erosion, hole erosion test

In this article, we present a numerical method to deal with fluid-solid interactions and simulate particle-fluid systems as encountered in soils. This method is based on a coupling between two methods, now widely used in mechanics of granular media and fluid dynamics respectively: the discrete element (DE) method and the lattice Boltzmann (LB) method. The DE method is employed to model interactions between particles, whereas the LB method is used to describe an interstitial Newtonian fluid flow. The coupling presented here is a full one in the sense that particle motions act on fluid flow and reciprocally. This article presents in details each of the two methods and the principle of the coupling scheme. Determination of hydrodynamic forces and torques is also detailed, and the treatment of boundaries is explained. The coupled method is finally illustrated on a simple example of piping erosion, which puts in evidence that the combined LB-DE scheme constitutes a promising tool to study coupled problems in geomechanics.

INTRODUCTION

The discrete element (DE) method [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] is now widely used for the investigation of the constitutive mechanical behavior of granular media such as soils, dealing, for instance, with questions about nonassociative flow rule [START_REF] Bardet | Numerical simulations of the incremental responses of idealized granular materials[END_REF][START_REF] Calvetti | A numerical investigation of the incremental behavior of granular soils[END_REF], shear band development [START_REF] Iwashita | Micro-deformation mechanism of shear banding process based on modified distinct element method[END_REF][START_REF] Tordesillas | Force cycles and force chains[END_REF], and failure [START_REF] Sibille | From bifurcation to failure in a granular material, a dem analysis[END_REF]. Although the computation cost constitutes a limit for the DE method, it has also been applied for the study of engineering problems (e.g., soil-inclusion interactions [START_REF] Chareyre | Dynamic spar elements and discrete element methods in two dimensions for the modeling of soil-inclusion problems[END_REF] or rock impacts on embankments [START_REF] Plassiard | Rockfall impact parameters on embankments: a discrete element method analysis[END_REF]).

A challenging point nowadays is to extend the application of the DE method to multiphysic issues. Concerning the framework of the fluid-solid interactions, the main efforts aim, on one hand, to take into account interparticle capillary forces in partially saturated granular media [START_REF] Scholtès | On the capillary stress tensor in wet granular materials[END_REF][START_REF] Richefeu | A model of capillary cohesion for numerical simulations of 3D polydisperse granular media[END_REF], and on the other hand, to describe the fluid dynamics in the pore space of soils and to establish a full mechanical coupling between the solid and the fluid phases.

Concerning this latter improvement, several attempts have been made where different scales were considered for the description of the fluid dynamics. On one hand, the fluid flow can be described at the scale of a cell including several solid particles [START_REF] Zeghal | A continuum-discrete hydromechanical analysis of granular deposit liquefaction[END_REF]; hence, the computational cost is reduced because the fluid dynamics is not resolved in each pore. However, because the method gives only access to an average flow over a subset of particles, the computation of drag forces acting on each solid particle requires strong hypotheses or semiempirical relations [START_REF] Zeghal | A continuum-discrete hydromechanical analysis of granular deposit liquefaction[END_REF]. On the other hand, other approaches aim to describe the fluid flow at a scale smaller than the pore or particle size [START_REF] Araújo | Critical role of gravity in filters[END_REF]. The objective was therefore to limit as much as possible the hypotheses made on the fluid flow and on the hydrodynamic forces acting on grains. In this latter case, the computation cost can rapidly become prohibitive.

The solid-fluid coupled method developed and presented in this article is designed to analyze and to improve the understanding of soil erosion. More particularly, the main objective is to investigate, at the scale of soil grains or aggregates, the relevant parameters involved in the phenomenon of solid particle detachment from the granular skeleton, under a water seepage. In this context, hydrodynamic forces acting on solid particles have to be accurately computed with limited empirical assumptions. One has therefore to properly describe the fluid flow through the interparticle space.

One of the main difficulties in modeling the fluid flow within the pore space of a granular assembly is to deal with the complex moving boundaries, constituted by the deformable solid granular skeleton. Such moving boundaries can be quite easily taken into account with the lattice Boltzmann (LB) method [START_REF] Succi | The Lattice Boltzmann Equation for Fluid Dynamics and Beyond[END_REF], constituting an important advantage of this method over other computational fluid dynamics (CFD) methods. Moreover, as for the DE method, the numerical implementation of the LB method is quite straightforward. Finally, because the LB method is versatile and can be adapted to model multiphase, free surface, thermal or reactive flows, the treatment of a wide range of future developments and applications is possible. Previous research works have presented such a coupling between the LB method and particle methods, mainly for the case of dilute suspension [START_REF] Komnik | Transport phenomena and structuring in shear flow of suspensions near solid walls[END_REF][START_REF] Harting | Lattice-boltzmann simulations of transport phenomena and structuring in suspensions[END_REF][START_REF] Feng | Coupled lattice Boltzmann method and discrete element modelling of particle transport in turbulent fluid flows: Computational issues[END_REF] and generally without detailing all the aspects and difficulties encountered in the implementation of a coupled LB-DE method scheme.

This article aims at presenting and describing the necessary details that should be taken into account to achieve an efficient coupling between the DE and the LB methods. For the sake of simplicity, the developments are here limited to a two-dimensional case. However, all the aspects of the method can be directly generalized to the three-dimensional case. In the first section of this article, the principle of the LB method is presented, with a special emphasis put on the definition of the boundary conditions and on their adaptation to deal with moving solid particles. The second section is devoted to the DE method, to point out the features of this method conditioning the solidfluid coupling, with, afterward, an emphasis on the numerical issues likely to be encountered when achieving the full coupling between the LB and DE methods. To highlight the potentiality of the proposed coupled method, we present in a final section an application concerning piping erosion in a granular assembly similar to the laboratory hole erosion test (HET) [START_REF] Wan | Investigation of internal erosion and piping of soils in embankment dams by the slot erosion test and the hole erosion test[END_REF].

In the different sections, the term particle can refer to several objects. In the DE method, a particle is a solid object interacting with other ones through contact forces. In the LB method, a particle is a microscopic entity carrying state information about the fluid system. These latter particles are submitted to statistical physical interactions. The differences between the two types of particles should always be kept in mind for a good understanding of further explanations.

LATTICE BOLTZMANN (LB) METHOD

The LB method is based on Boltzmann's equation [START_REF] Boltzmann | Weitere Studien über das Wärmegleichgewicht unter Gas-molekülen[END_REF], derived from the gas kinetic theory. In this work, Boltzmann introduced a binary collision operator to describe the time and spatial evolution of a distribution function of particles. Later, Enskog and Chapman [START_REF] Enskog | Kinetische theorie der vorgänge in mässig verdünnten gasen[END_REF][START_REF] Chapman | The Mathematical Theory of Non-uniform Gases: An Account of the Kinetic Theory of Viscosity[END_REF] showed direct links between the Boltzmann equation [START_REF] Boltzmann | Weitere Studien über das Wärmegleichgewicht unter Gas-molekülen[END_REF] and the Navier-Stokes equations, using a more complex collision operator to deal with multiple collisions. In 1954, a simplified collision operator, called BGK, was introduced by Bhatnagar-Gross-Krook [START_REF] Bhatnagar | A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems[END_REF]. The BGK operator is based on a time relaxation of the distribution function of particles toward an equilibrium function to describe particle collisions.

Besides, in 1973, the first lattice gas cellular automata (LGCA) model was introduced by Hardy and coworkers [START_REF] Hardy | Time evolution of a two-dimensional model system. I. Invariant states and time correlation functions[END_REF][START_REF] Hardy | Molecular dynamics of a classical lattice gas: Transport properties and time correlation functions[END_REF]. The basic idea of the LGCA models is to simulate the macroscopic behavior of a fluid flow with a very simplified model of the microscopic interactions between particles, where time, space, and particle velocities are discretized. The principle of the LGCA model consists of two sequential steps: collision and propagation. The collision step makes the direction of particle velocities change, whereas the propagation step lets particles move to the nearest node in the direction of its velocity. Links between LGCA and Boltzmann equation were established [START_REF]Theory and Applications of Cellular Automata[END_REF][START_REF] Frisch | Lattice Gas Hydrodynamics in Two and Three Dimensions[END_REF] and lead to the introduction of the lattice Boltzmann equation (LBE) [START_REF] Mcnamara | Use of the Boltzmann Equation to Simulate Lattice-Gas Automata[END_REF]. The traditional LB model was then born by adopting the BGK collision operator [START_REF] Chen | Recovery of the Navier-Stokes equations using a lattice-gas Boltzmann method[END_REF][START_REF] Qian | Lattice BGK Models for Navier-Stokes Equation[END_REF], which can lead to Navier-Stokes equations using the Chapman-Enskog expansion [START_REF] Succi | The Lattice Boltzmann Equation for Fluid Dynamics and Beyond[END_REF][START_REF] Chin | Lattice Boltzmann simulation of the flow of binary immiscible fluids with different viscosities using the Shan-Chen microscopic interaction model[END_REF]. Later, He and Luo [START_REF] He | A priori derivation of the lattice Boltzmann equation[END_REF] directly derived the LBE from the Boltzmann equation. This a priori derivation placed the LBE on a rigorous theoretical foundation with no use of the lattice gas automata. We describe in this section the application of the BGK approximation of the LBE and the implementation of boundary conditions (fluid pressure or velocity conditions, no-slip condition on moving particles) for its coupling with the DE method.

Principle of the LB method

In the LB method, the whole domain (solid and fluid) is discretized using a regular lattice.

Here, we employ the widely used two-dimensional D2Q9 model [START_REF] Qian | Lattice BGK Models for Navier-Stokes Equation[END_REF] based on the Face Centered HyperCube model of D'Humières et al. [START_REF] Humières | Lattice Gas Models for 3D Hydrodynamics[END_REF]. The fluid particles located at each node of a square lattice are allowed to move only along eight directions, toward their eight nearest neighbors. As illustrated in Figure 1, these displacements can occur with nine different discrete velocities e i with i = 0, 1, 2 . . . 8, where the case i = 0 corresponds to still particles. According to i, the velocity vectors in the D2Q9 model are given by [START_REF] Ladd | Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation[END_REF][START_REF] Behrend | Solid-fluid boundaries in particle suspension simulations via the lattice Boltzmann method[END_REF][START_REF] Mei | Force evaluation in the lattice Boltzmann method involving curved geometry[END_REF][START_REF] Succi | The Lattice Boltzmann Equation for Fluid Dynamics and Beyond[END_REF]

e i =              (0, 0) if i = 0 C cos π(i-1)

2

, sin π(i-1)

2 for i = 1, . . . , 4 √ 2C cos π(2i-9)

4

, sin π(2i-9)

4 for i = 5, . . . , 8 (1)
where C = h/dt is the lattice velocity, with dt as the time step and h as the lattice spacing. For each direction of every node of the lattice, a particle distribution function or density function f i is associated. f i ( x, t) defines the proportion of particles, moving with velocity e i , along the ith direction of the node located at position x at time t.

Collision step.

An LB algorithm is an iterative process over time. Each time step can be divided into two essential steps, namely, collision and propagation. During the collision phase, each particle distribution function f i evolves toward an equilibrium distribution function f eq i . For a Newtonian fluid, the equilibrium functions are given by [START_REF] Qian | Lattice BGK Models for Navier-Stokes Equation[END_REF]

f eq i = w i ρ 1 + 3 C 2 e i . v + 9 2C 4 ( e i . v) 2 - 3 2C 2 v. v for i = 0, . . . , 8 (2)
where w 0 = 4/9, w 1,2,3,4 = 1/9 and w 5,6,7,8 = 1/36 [START_REF] Qian | Lattice BGK Models for Navier-Stokes Equation[END_REF]. v and ρ are respectively the macroscopic fluid velocity and density of the considered node.

In the rest of this article, overlined characters refer to dimensionless variable. Using ρ = ρ/ρ 0 , ¯ v = v/C and ¯ e i = e i /C, we can write

f eq i = f eq i ρ 0 = w i ρ 1 + 3 ¯ e i . ¯ v + 9 2 ( ¯ e i . ¯ v) 2 - 3 2 ¯ v. ¯ v for i = 0, . . . , 8 (3)
where ρ 0 is the fluid density.

As said previously, this relaxation is performed here with the single-relaxation BGK model [START_REF] Bhatnagar | A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems[END_REF]:

f i ( x, t + ) = f i ( x, t) - 1 τ [f i ( x, t) -f eq i ( x, t)] (4)
The right side of Equation ( 4) is usually called the post-collision term f i ( x, t + ). t + is an intermediate time just following the collision process (t < t + < t + dt). The parameter τ is a dimensionless relaxation time. τ is related to the kinematic viscosity of the fluid ν, the lattice velocity C, and the lattice spacing h through the following relation:

ν = 1 3 τ - 1 2 Ch (5)
Practically, for given values of fluid viscosity and parameter τ , dt is defined according to the chosen lattice spacing by

dt = 1 3ν τ - 1 2 h 2 (6)
Macroscopic fluid properties like density ρ and velocity v can be retrieved at each node by

ρ = 8 i=0 f i (7) v = 1 ρ 8 i=0 f i e i (8)
To determine the fluid pressure p, the fluid is assumed to be slightly compressible. We use the following state equation p = c 2 s ρ, where c s is the sound celerity and is defined in the D2Q9 model as c s = C √ 3 . Macroscopic quantities converge to the solution of the incompressible Navier-Stokes equation with order M 2 , where M is the computational Mach number defined by

M = v max C , (9)
where v max is the highest velocity in the flow. To ensure good result accuracy, it is therefore required that

M 1. (10)
Practically, it is convenient to use a value of M 0.1, which leads at worst to a value of 1% for M 2 , which is representative of the compressibility error.

Propagation step.

The collision step, presented previously, consists in computing the right-hand side of Equation (4). After this step, the propagation phase aims at propagating the post-collision distribution functions f i ( x, t + ) over the lattice grid. This step can be represented by the following equation:

f i ( x + e i dt, t + dt) = f i ( x, t + ) (11)
In Equation [START_REF] Zeghal | A continuum-discrete hydromechanical analysis of granular deposit liquefaction[END_REF], for a node located at x, x + e i dt is the position of the nearest neighbor node along the ith direction.

During one time step of the LB algorithm, the two successive collision and propagation phases permit the particle distribution functions to tend toward equilibrium functions, which depend on both the density and the velocity at the node, before they are propagated along each direction. Then, new values of the equilibrium functions are calculated with respect to the new velocity and density at each node, and a new time step begins. Initially, every particle distribution functions f i are set equal to their equilibrium function with the fluid at rest for all nodes (i.e., including distribution functions involved at fluid and solid boundary nodes as they will be defined in Section 2.2.1).

2.1.3. Setting parameters of the LB method.

As seen previously, the LB method comes with the definition of three parameters: h, τ and dt. These parameters are linked together as well as to the fluid viscosity through Equation [START_REF] Sibille | From bifurcation to failure in a granular material, a dem analysis[END_REF]. Thus, they have to be chosen in balance between the desired fluid viscosity, the numerical solution accuracy, and the computational cost. More precisely, the lattice spacing h has a direct impact on the central processing unit (CPU) time used for each step of the LB algorithm. For this reason, h has not to be too small. Nevertheless, the spatial discretization error is directly related to h. The time step dt affects the total number of time step required to reach the end of a simulation. However, in the same time, dt has to be relatively small to minimize time discretization error. The last, and not the least parameter, is the relaxation time τ . τ is the time scale of molecular collisions between fluid particles and is linked to ν, h, and dt through Equation ( 5). Therefore, for a given fluid viscosity, these three parameters have to be set in balance with the spatial discretization and the minimization of the compressibility error and with a reasonable computation time. Practically, we independently choose τ and h and derived dt through Equation ( 6). The choice of adequate values for h and τ constitutes the major issue, but unfortunately, there is no a priori way to estimate them. A testcorrection procedure can be used until the compressibility criterion (Equation ( 10)) is satisfied, and an accurate solution is obtained. The only restriction on τ value comes from Equation [START_REF] Tordesillas | Force cycles and force chains[END_REF]. Indeed, Equation [START_REF] Tordesillas | Force cycles and force chains[END_REF] implies that τ > 0.5, but no upper limit, derived from LB theory, exists. Then, τ can be set to any desired value, even greater than unity [START_REF] Inamuro | A non-slip boundary condition for lattice Boltzmann simulations[END_REF][START_REF] Chen | On boundary conditions in lattice Boltzmann methods[END_REF][START_REF] Gallivan | An evaluation of the bounce-back boundary condition for lattice Boltzmann simulations[END_REF]. Smaller values of τ imply smaller time steps and therefore increase computational cost. As discussed in Ref. [START_REF] Gallivan | An evaluation of the bounce-back boundary condition for lattice Boltzmann simulations[END_REF], results are acceptable for τ < 1.5, especially with the boundary conditions used here (see the following section).

Boundary conditions on obstacles and system borders

The two-step scheme collision-propagation presented previously constitutes the basis of the LB method. Nevertheless, specific treatments need to be achieved for both system boundaries and moving solid obstacles.

Boundary conditions on solid obstacles.

Solid particles, or any other sort of obstacles, can be handled with the LB method by positioning them on the lattice. The nodes belonging to an obstacle which are the closest ones to its external surface are called solid boundary nodes (SB nodes). All the fluid nodes located outside of the obstacle, with at least one direction pointing toward an SB node, are called fluid boundary nodes (FB nodes). All other nodes act as fluid nodes and will be named like this. A link between an SB and an FB node is named a boundary link. For example, Figure 2 illustrates the distinction between SB and FB nodes for a solid circular obstacle. The simplest scheme to calculate momentum transfer between an SB and an FB node is the bounce-back boundary condition on the boundary links. With this method, the effective solid boundary is assumed to be located at the middle of the boundary links as presented in Figure 2. More accurate boundary representation can also be implemented [START_REF] Bouzidi | Momentum transfer of a Boltzmann-lattice fluid with boundaries[END_REF][START_REF] Strack | Three-dimensional immersed boundary conditions for moving solids in the lattice-Boltzmann method[END_REF], but such improvements are not used in the present article. The bounce-back scheme has been validated by comparing the velocity profiles obtained respectively with our LB method and with the commercial CFD software FLUENT for a fluid flowing in a channel around a fixed cylinder (Figure 3).

The LB simulations were run for different resolutions d/h of the lattice, where d is the cylinder diameter. Figure 3 illustrates that a good agreement can be obtained for d/h > 10.

In this work, we choose to use the mid-plane bounce-back scheme in a first approximation for the sake of simplicity and to limit the computation cost.

Fluid particles cannot cross the solid border line that limits the outer solid domain. Therefore, particle distribution functions coming from boundary nodes (SB or FB) are reflected along a boundary link. For a stationary solid surface, the distribution functions at t + dt on a boundary link σ coming from direction i (see Figure 4) are given by

f -σi ( x F B , t + dt) = f σi ( x F B , t + ) f σi ( x SB , t + dt) = f -σi ( x SB , t + ) (12)
In Equation ( 12), -i simply denotes the opposite direction of i.

When the solid boundary is moving, the surface velocity needs to be taken into account in the bounce-back rule to ensure the no-slip boundary condition on the solid surface. This is achieved by using the modified bounce-back rule [START_REF] Ladd | Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation[END_REF]:

f -σi ( x F B , t + dt) = f σi ( x F B , t + ) -2α i V b . e i f σi ( x SB , t + dt) = f -σi ( x SB , t + ) + 2α i V b . e i (13) with x SB = x F B + e σi dt,
where α i is defined by α i = 3w i ρ/C 2 and V b is the boundary velocity at the middle of the boundary link σ:

V b = V c + w ∧ r c (14) with
• V c and w as the translational and angular velocities respectively at the center of mass of the solid obstacle, • r c as the vector joining the center of mass x c of the solid obstacle to the middle of the boundary link:

Figure 2. A solid circular obstacle mapped on the LB lattice with the different node species.

r c = x F B + e σi dt 2 -x c (15)
2.2.2. Boundary conditions on system boundaries.

In the LB method, pressure or velocity boundary conditions cannot be directly imposed. As pressure and velocities are derived from the distribution functions f i , distribution functions have to be defined to match the desired boundary condition. We consider here a rectangular channel as presented in Figure 5. In this work, we use boundary conditions similar to those defined by Zou and He [START_REF] Zou | On pressure and velocity boundary conditions for the lattice Boltzmann BGK model[END_REF]. These conditions are based on the idea of the bounce-back scheme of the nonequilibrium part of the particle distribution functions. The channel presented in Figure 5 requires defining eight sets of distribution functions denoted respectively from a to g. We consider first a node at the inlet of the system (case a in Figure 5). After the propagation step, f 2 , f 3 , f 4 , f 6 , f 7 are known and f 1 , f 5 , f 8 need to be determined to impose the boundary condition. We assume that the bounce-back rule is still correct for the nonequilibrium part of the particle distribution function normal to the boundary. Thus, f 1f eq 1 = f 3f eq 3 . Using the expression of equilibrium distribution functions given by Equation (3), it comes with dimensionless variables for the sake of simplicity: By projecting Equation ( 8) on the two space directions, we can write

e -σi e σi x F B x SB = x F B + e σi dt
f1 + f5 + f8 = ρv x + f3 + f6 + f7 f5 -f8 = ρv y + -f2 + f4 -f6 + f7 (17)
Then, combining Equations ( 17) and ( 16) leads to

       f1 = f3 + 2 3 ρv x f5 = f7 -1/2( f2 -f4 ) + 1/6ρv x + 1/2ρv y f8 = f6 + 1/2( f2 -f4 ) + 1/6ρv x -1/2ρv y (18)
The same procedure can be used for the cases b, c and d of Figure 5. The relations used to determine the unknown particle distribution functions in these cases are given in appendix. Due to a different number of unknown distribution functions on the corner nodes (cases e, f, g and h), a specific procedure is needed. For example, let us take the bottom left corner node (case e). At this node, f 3 , f 4 and f 7 are known after the propagation step, but f 1 , f 2 , f 5 , f 6 and f 8 need to be determined. Using the bounce-back rule for the nonequilibrium part, we find that

     f1 = f3 + ( f eq 1 -f eq 3 ) = f3 + 2 3 ρv x f2 = f4 + ( f eq 2 -f eq 4 ) = f4 + 2 3 ρv y (19)
In addition, by developing Equations ( 7) and ( 8), we have

           f5 -f6 + f8 = ρv x -( f1 -f3 -f7 ) f5 + f6 -f8 = ρv y -( f2 -f7 -f4 ) f6 + f8 = ρ -( f0 + f1 + f2 + f3 + f4 + f5 + f7 ) (20)
Then, by combining Equation [START_REF] Enskog | Kinetische theorie der vorgänge in mässig verdünnten gasen[END_REF] with Equation ( 20), we obtain the expressions of f5 , f6 and f8 :

           f5 = f7 + 1/6ρ( vx + vy ) f6 = 1/2[ρ(1 -vx -2/3v y ) -f0 -2( f3 + f4 + f7 )] f8 = 1/2[ρ(1 -2/3v x -vy ) -f0 -2( f3 + f4 + f7 )] (21)
Here, again, the same procedure can be applied to cases f, g and h of Figure 5. All the relationships between distribution functions are listed in appendix for the D2Q9 model. From a practical point of view, to apply the boundary condition defined by Equation ( 18) for the case a at the inlet of the channel, ρ, vx and vy have to be fixed. If a pressure condition p in is imposed, the inlet density ρ in is deduced from the relation p in = c 2 s ρ in . For convenience, we define v y = 0. Using Equation ( 7) it comes

f1 + f5 + f8 = ρin -( f0 + f2 + f3 + f4 + f6 + f7 ) (22)
by combining with Equation ( 17) with ρ = ρ in , we obtain the velocity at the inlet:

vx = 1 - f0 + f2 + f4 + 2( f3 + f6 + f7 ) ρin (23)
The same procedure can be applied to impose velocity boundary conditions. In such a case, considering the inlet (case a), vx and vy are known at the inlet. Then, the pressure (density) needs to be determined first prior to the determination of f1 , f5 and f8 ( f2 , f3 , f4 , f6 and f7 are known after streaming). Using Equation ( 23), ρin can be determined with

ρin = 1 1 -vx f0 + f2 + f4 + 2 f3 + f6 + f7 (24)
Then, the determination of unknown distribution functions is similar to the one used for pressure boundary conditions from Equation [START_REF] Boltzmann | Weitere Studien über das Wärmegleichgewicht unter Gas-molekülen[END_REF]. A similar reasoning can be done for every channel boundary, and we will present in Section 5 an example of use of such boundary conditions applied to a channel flow. The treatment of boundary nodes, including those on the corner, will be detailed.

DISCRETE ELEMENT (DE) METHOD

To model solid particle interactions, we use the YADE open source software [START_REF] Kozicki | A new open-source software developed for numerical simulations using discrete modeling methods[END_REF][START_REF] Yade | Open source discrete element method[END_REF] based on the DE method, as introduced by Cundall and Strack [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]. Particles move according to Newton's law, and two particles in contact are allowed to slightly overlap among each other. Contact interaction forces are computed for each couple of particles in contact according to a cohesive frictional interaction law.

In the direction normal to the tangent contact plane, the interaction law is elastic and characterized by the normal stiffness k n :

F n = k n δ n ( 25
)
where F n is the normal contact force and δ n is the overlap between the couple of particles in contact. Due to the introduction of a normal cohesion C n at contact, tensile normal forces are allowed for F n C n (compressive forces F n are counted positively, then C n takes usually negative values). If F n < C n , failure occurs and the contact is lost.

In the direction included in the tangent contact plane, the interaction law is elastic-plastic and follows the Coulomb friction law. In the elastic regime, k s represents the shear stiffness:

∆F s = -k s ∆u s (26)
where ∆F s is the shear force increment and ∆u s is the increment of the tangential relative displacement of the particles at contact. Tangential sliding occurs if

|F s | > F n tan(Φ) + C s ( 27
)
with Φ is the contact friction angle and C s is the shear cohesion. When sliding occurs, normal C n and shear C s cohesions are set to zero. Consequently, the contact becomes then purely frictional. It means that sliding occurs for F s > F n tan(Φ) and contact is lost for F n 0, as presented in Figure 6. At each time step of the DE method, updated position and velocity of each solid particle are deduced from the explicit time integration of Newton's second law, taking into account the interparticle contact force F c and the hydrodynamic force F h acting on the considered particle (see Section 4 for the definition of F h ):

Fs C s Fn C n Φ Φ ×Ð Ò Ò Ó × ÓÒ ÖÓ Ò C n = C s = 0 ÓÒØ Ø ÐÓ×Ø
F c + F h = m ẍc ( 28
)
where m is the mass of the solid particle and x c is its position. Rotation of particles is computed in the same way by integrating the equation:

T c + T h = J ẇ (29)
where T c is the torque resulting from the interparticle forces, T h is the hydrodynamic torque, J is the moment of inertia (for circular or spherical solid particles), and w is the angular velocity. The stability of the explicit scheme for the time integration holds if the time step of the DE method dt DE is sufficiently small to describe the wave propagation in the mass-spring system constituted by the solid granular assembly. Hence, the critical time step dt cr DE is related to the period of oscillation of this assembly. For a single one-dimensional mass-spring oscillator of mass m and stiffness k, the period of oscillation and thus dt cr DE are equal to 2π m/k. In YADE, dt cr DE is computed from a generalization of this definition for all the degrees of freedom of the particles, in rotation and in translation [START_REF] Šmilauer | Yade Documentation[END_REF][START_REF] Šmilauer | Yade dem formulation[END_REF]. A ratio of 10 between dt DE and dt cr DE is usually enough to avoid any instability of the numerical method.

COUPLING BETWEEN LB AND DE METHODS

The coupling between the LB and DE methods requires the exchange of information from one method to the other. Hydrodynamic forces and torques are deduced from the LB computation and used in the DE cycle to compute new positions and velocities of solid particles, which in turn define the new boundary conditions for the fluid dynamic problem. This has to be done by taking care of the fact that each method is characterized by a different time discretization. Issues concerning the computation of hydrodynamic forces, numerical time step adjustment, and modification of the fluid domain due to new particle positions (resulting in the affectation of new lattice nodes in the fluid phase) are discussed hereafter.

Forces acting on solid obstacles

The fluid force acting on an obstacle can be evaluated from the momentum exchange between fluid and solid particle surface [START_REF] Ladd | Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation[END_REF][START_REF] Mei | Force evaluation in the lattice Boltzmann method involving curved geometry[END_REF][START_REF] Bouzidi | Momentum transfer of a Boltzmann-lattice fluid with boundaries[END_REF] along two opposite propagation directions. If we consider a boundary link σ with direction i pointing from an FB node to an SB node (see Figure 4), using Equation ( 13), the momentum exchange Π is

Π( x, t + 1 2 dt) = e σi f σi ( x, t + ) -e -σi f -σi ( x, t + dt) = e σi [f σi ( x, t + ) + f -σi ( x, t + dt)] = 2 e σi f σi ( x, t + ) -α i V b . e i (30)
Then, the corresponding force and torque are respectively

F σ ( x, t + 1 2 dt) = 2 Ω dt f σi ( x, t + ) -α i V b . e i e σi (31)
T σ ( x, t + 1 2 dt) = r c × F σ ( x, t + 1 2 dt) ( 32
)
where Ω is the cell lattice volume in three dimensions or the cell lattice surface in two dimensions. The total hydrodynamic force and torque exerted on the solid particle are obtained by summing Equations ( 31) and ( 32) over all the boundary links of the particle:

F h (t + 1 2 dt) = σ F σ ( x, t + 1 2 dt) (33)
T h (t + 1 2 dt) = σ T σ ( x, t + 1 2 dt) (34)

Numerical coupling

The idea behind the coupling of the DE and LB methods is a subcycling scheme.

dt DE = dt n , ( 35
)
where n is the nearest integer greater or equal to dt/dt DE . It has been chosen here to reduce the DE time step to be sure to satisfy automatically the critical time step. Thus, in one computational step of the LB method, n substeps of time integration with the DE method are performed. During this subcycling, the hydrodynamic forces and torques remain unchanged. As the time step of the LB method is larger than the DE time step, it means that the LB method can describe physical phenomena characterized by a lower frequency than those that could be described with the DE method. Similarly, the highest frequency of the solid-fluid interaction described by such a coupling should depend essentially on the time step of the LB method. Although further investigations are required to understand clearly how large n can be without altering LB-DE coupling, we choose n values smaller than 10. This condition gives another important constraint on the determination of an appropriate value of the relaxation time τ (as discussed in Section 2.1.3). To be time consistent, hydrodynamic forces and torques applied on obstacles are obtained by an average over two consecutive time steps of the LB method:

F h (t) = 1 2 F h (t + 1 2 dt) + F h (t - 1 2 dt) (36)
T h (t) = 1 2 T h (t + 1 2 dt) + T h (t - 1 2 dt) (37)
The DE method software YADE is a modular program that lets the user add some features with plug-ins. We have developed a plug-in to resolve fluid flow with the LB method and to perform a full coupling with the DE method implemented in YADE. Due to the fact that our plug-in execution is driven by that DE loop, our LB plug-in is activated every n DE time steps to fulfill the subcycling condition. In addition, from a given time t, the execution of the LB time step and the set of n DE time steps are performed sequentially. First, the LB plug-in is run for one time step to compute the fluid flow at t + dt. Consequently, hydrodynamic forces and torques at time t can be deduced from Equations ( 36) and [START_REF] Gallivan | An evaluation of the bounce-back boundary condition for lattice Boltzmann simulations[END_REF]. Then, the DE method is run for n time steps to compute the positions of solid particles at t + dt. Hydrodynamic forces and torques, previously computed, are applied to solid particles in addition to the intergranular contact forces at every DE time step, during these n steps.

Due to the movement of solid particles on the LB lattice, lattice nodes initially covered by solid particles may become fluid nodes during simulation. To avoid the creation of new distribution functions when new fluid nodes are uncovered [START_REF] Iglberger | Simulation of moving particles in 3D with the Lattice Boltzmann method[END_REF], all nodes are initialized at the simulation startup, even those included in solid particles. Thus, a solid particle consists of a nutshell of mass and inertia, equal to the mass and inertia of the whole solid particle, filled with fluid of same density as the bulk fluid. The presence of inner fluid guarantees mass conservation locally [START_REF] Aidun | Lattice Boltzmann simulation of solid particles suspended in fluid[END_REF][START_REF] Ladd | Lattice-Boltzmann Simulations of Particle-Fluid Suspensions[END_REF] and counterbalances error in force and torque calculations. Therefore, no errors are induced by this method [START_REF] Clausen | Galilean invariance in the lattice-Boltzmann method and its effect on the calculation of rheological properties in suspensions[END_REF].

APPLICATION TO PIPING EROSION

Piping erosion and HET

The HET [START_REF] Wan | Investigation of internal erosion and piping of soils in embankment dams by the slot erosion test and the hole erosion test[END_REF][START_REF] Pham | Erosion et dispersion des sols argileux par un fluide[END_REF][START_REF] Regazzoni | Confrontation et analyse d'érodimètres et caractérisation de la sensibilité à l'érosion d'interface[END_REF] is used in civil engineering to characterize erodability of soils with respect to piping erosion. A cylindrical compacted soil sample, drilled along its axis, is placed in a channel. This channel is saturated with water, and a pressure gradient is applied between the inlet and the outlet (see Figure 7). Due to hydraulic loading, soil particles can be detached and transported through the soil pipe.

Classically, the interpretation of results from an HET experiment [START_REF] Wan | Investigation of internal erosion and piping of soils in embankment dams by the slot erosion test and the hole erosion test[END_REF] is based on the erosion law suggested by Shields [START_REF] Shields | Anwendung der Aechichkeits-mechanik und der turbuleng forschung aufdir Geschiebewegung[END_REF]:

˙ = k d (τ h -τ c ) if τ h > τ c , ( 38
)
where ˙ is the rate of eroded mass per unit pipe area, k d is the erosion coefficient, τ h is the fluid shear stress, and τ c is a critical shear stress. To deduce the erosion parameters (k d and τ c ) from such an experiment, the erosion rate ˙ and the shear stress τ h have to be estimated. However, generally, the evolution of the eroded mass is not followed during the experiment time † . Consequently, the change in times of the soil pipe diameter is deduced from the initial and final diameters, then the erosion rate is estimated from the rate of the pipe diameter growth. The fluid shear stress is computed from the fluid force summation on the boundaries of the pipe of radius R at time t. In the case of a two-dimensional pipe (case studied hereafter) of height d h and length l h , this latter expression reads

τ h = d h ∆P 2l h , ( 39
)
where ∆P is the pressure gradient between the inlet and the outlet of the apparatus. In this section, a two-dimensional simulation of piping erosion inspired by HET is presented using the developed LB-DE method previously presented. The objective is to verify if the linear relation between the erosion rate and the fluid shear stress can be described by the coupled model.

Piping erosion modeling with the coupled LB-DE model

To create the numerical model, a packing of cohesive frictional disks is generated and compacted (by enlarging disk diameters) in a rectangular box of length l, as illustrated in Figure 7. This packing is made of monosized grains of diameter d with a size dispersion of 30% and a mean diameter d mean . The cohesive frictional contact law, used in this work, is the one presented previously in Figure 6. The granular packing is drilled by removing grains, which are not in the hatched regions of Figure 7, and the left and right walls are removed. This system can be seen as a simplified model of HET initial configuration with a hole of diameter d h and without taking gravity into account. Finally, the system is spatially discretized with a square lattice to set up the LB algorithm and to let the voids fill up with water. In all the simulations, we use the parameters listed in Table I. Therefore, in these simulations, the mean number of node along a particle radius is d mean /h ≈ 17. Because of the twodimensional description adopted here, there is no water flow between two particles in contact and consequently through the bulk of the granular assembly in the hatched regions (Figure 7). This point is certainly the main limitation with respect to a three-dimensional model. However, this hypothesis seems reasonable in such a case of piping erosion, where water flow through the bulk of the soil is generally negligible with respect to the water flow in the pipe.

Geometry DE Method LB Method

l = 0.03 m and L = l/3 C n /d = C s /d = 0.253 N m -1 N x = 1003, N y = 335 d mean = 5.06 × 10 -4 m Φ = 20 • ν = 1 × 10 -6 m 2 s -1 d h = 0.2 × l/3 k n /d = 150 × 10 6 N m -2 ρ 0 = 1 000 kg m -3 l h = 0.8 × l k s /k n = 0.4 dt = 1.807 × 10 -4 s, n = 3, τ = 1.1
Table I. Table of parameters used in our simulations. Nx and Ny are respectively the number of nodes along the x and y directions.

To initiate fluid flow throughout the pipe, a pressure gradient is imposed between the inlet and the outlet of the system:

∆P = p in -p out = c 2 s ρ 0 (ρ in -ρout ) (40)
with ρout = 1. Results presented in this article will be expressed in terms of hydraulic gradient i = ∆P/(ρ 0 gl h ) = 10∆P/(ρ 0 g8l), where g = 9.81 m.s -2 is the gravitational acceleration. To apply the pressure gradient, we use pressure boundary conditions as described in Section 2.2.2. On the inlet and the outlet, we suppose that v y = 0. Then v x is set according to Equation [START_REF] Hardy | Molecular dynamics of a classical lattice gas: Transport properties and time correlation functions[END_REF]. For the bottom and upper fixed walls, velocity boundary conditions (v x = v y = 0) are imposed. ρ value is fixed according to Equation [START_REF]Theory and Applications of Cellular Automata[END_REF]. Here, corner nodes are nodes belonging to fixed walls; therefore, they are initialized in a similar way to wall nodes (v x = 0 and v y = 0). Nevertheless, on corner nodes, the locally available information is insufficient for the evaluation of the density. In such a case, it is common to take the density of corner nodes equal to the one of the nearest node [START_REF] Zou | On pressure and velocity boundary conditions for the lattice Boltzmann BGK model[END_REF]. For each corner node, we set ρ equal to the density of the nearest node, which belongs to the same wall.

The series of snapshots presented in Figure 8 are taken from one simulation run for i = 12.7 × 10 -4 . Figure 8 shows that particles are detached under the hydraulic loading. Then, detached particles are transported by the flow toward the channel outlet. We denote R max the largest radius of the particles. If a particle mass center comes closer than a distance R max to the channel output, the particle is simply removed from the simulation domain, and it is then considered as being an eroded one. When removing one particle, no specific treatment is needed to deal with the resulting fluid nodes because they were already treated as fluid nodes when they belonged to the solid particle. The status of the solid SB and fluid FB boundary nodes, on each side of the solid particle perimeter, is simply switched to classical fluid nodes.

During a simulation, we record the evolution of the eroded particle mass M e (t) and define the ratio M e (t)/M 0 , where M 0 is the total mass of the initial granular assembly. Before studying dependency of eroded mass on pressure gradient applied, let us consider some basic results to validate our simulations. Figure 9(a) shows the total hydrodynamic force acting on the solid phase. It illustrates that the total force along x direction is conserved and equal to L∆P . In a more rigorous way, the fluid compressibility should be taken into account to retrieve equality between the total hydrodynamic force and L∆P . Figure 9(b) presents the total mass variation during the simulation, determined with a summation of density over all nodes. It can be observed that, during a simulation, the maximum mass variation is about 1.7%. As explained previously, our model implies that fluid is considered as slightly compressible. Nevertheless, Figure 9(b) puts in evidence that, thanks to a careful choice of parameters (see Section 2.1.3), the compressibility of the fluid is negligible, and therefore, that mass conservation is achieved. with time for different values of the hydraulic gradient i. The ratio M e (t)/M 0 increases with time to reach unity, which corresponds to the case where all particles are eroded. Moreover, Figure 10 illustrates that erosion rate is higher when hydraulic gradient increases. For the lowest hydraulic gradient, detachment of particles is not triggered, and the mass of eroded particles is nil all along the simulation. Using fluid velocity profile along the direction perpendicular to the cylinder axis, the shear stress τ h at the vicinity of the hole boundary can be computed as

τ h = νρ 0 dV x dy (41)
Here, the velocity gradient dVx dy at the vicinity of the hole boundary is deduced from a velocity cross profile averaged along the length of the hole. It is then assumed that maximum velocity gradient, identified on both sides of this averaged velocity profile, results from the presence of the hole boundaries. Therefore, this latter maximum value of dVx dy has been chosen to compute the shear stress according to Equation [START_REF] Kozicki | A new open-source software developed for numerical simulations using discrete modeling methods[END_REF]. τ h can be calculated based on hole diameter measurement as defined in Section 5.1. Contrary to previous experiments performed by Wan and Fell [17], our simulations let us access to the hole diameter at any time during the piping erosion process. Here, the mean hole diameter is determined from an averaged cross-section density profile over the whole packing. Then, τ h is calculated using Equation [START_REF] Strack | Three-dimensional immersed boundary conditions for moving solids in the lattice-Boltzmann method[END_REF]. To overcome wall influence when estimating τ h and τ h , we do not consider shear stresses when ratio between hole diameter d h and channel height is greater than 75%.

The way to compute the eroded mass M e , as the cumulative sum of the detached particles reaching the outlet of the channel, introduces a delay in the detection of the particles effectively eroded, which corresponds to the time of transfer from the location of detachment of a particle to the outlet of the channel. Nevertheless, the characteristic transfer time of detached particles is less than 3% of the total simulated time to erode the whole granular packing and corresponds to a maximum hole diameter increase of 5%. Therefore, the computation of the eroded mass as it is done here at a time t can be assumed as concomitant with the hole diameter and the hydraulic shear stress determined at this same time t. Figures 11 and 12 present respectively dependency of ˙ on shear stresses τ h and τ h for simulations of piping erosion with different hydraulic gradients. In both Figures 11 and12, a least squares linear adjustment is represented. It illustrates that the linear erosion law introduced by Shields [START_REF] Shields | Anwendung der Aechichkeits-mechanik und der turbuleng forschung aufdir Geschiebewegung[END_REF] and used by Wan and Fell [START_REF] Wan | Investigation of internal erosion and piping of soils in embankment dams by the slot erosion test and the hole erosion test[END_REF] can be fairly retrieved with this coupled LB-DE model in a first approximation. With a more accurate observation, in particular of Figure 11, the ˙τ h relation slightly deviates from a linear tendency. However, this result is not in contradiction with experimental results obtained by Pham [START_REF] Pham | Erosion et dispersion des sols argileux par un fluide[END_REF] for some soil fabrics, where a slightly nonlinear ˙τ h relation was exhibited. The erosion coefficient k d corresponds to the slopes of linear adjustments. We can observe that two different values of k d were obtained according to the method used in determining the shear stress. The erosion coefficient obtained with τ h (k d = 4.7 s.m -1 ) is about twice small as the one obtained with τ h (k d = 10.4 s.m -1 ). The pressure gradient used in Equation [START_REF] Strack | Three-dimensional immersed boundary conditions for moving solids in the lattice-Boltzmann method[END_REF] gradient between the inlet and the outlet of the channel. It has been shown experimentally [START_REF] Marot | An Energy Based Method for Providing Soil Surface Erodibility Rankings[END_REF] that water head losses resulting from the narrowing (at the entrance of the channel) and the enlargement (at the exit) of the water flow width cannot be neglected. Hence, Marot et al. [START_REF] Marot | An Energy Based Method for Providing Soil Surface Erodibility Rankings[END_REF] estimated, for their HET experiments, that only about 30% of the pressure gradient applied between the inlet and the outlet of the channel was effectively applied between the entrance and the exit of the soil pipe. Such water head losses could explain the lower value found for τ h than for τ h . Nevertheless, the ratio between the initial pipe diameter d h and the channel diameter L (d h /L = 0.2) is much larger in the numerical model than for the real classical HET (where the ratio is typically about 0.06) and not any meaningful pressure drop has been observed at the entrance and exit of the soil pipe. Beyond that point, the main geometrical difference between this simplified numerical model and classical HET experiments on soils is the ratio between the typical size of soil particles and the pipe width.

As the diameter of solid particles is here important with respect to the pipe width, first, the assumption of flat pipe surface does not hold, and second, the influence on water flow of particles detached and being carried away cannot be neglected. As pipe surface is quite tortuous, the action of fluid on pipe walls cannot be resumed to a single shear stress in the direction of the main fluid flow (i.e., the horizontal direction in Figure 7 or 8). In addition, due to the relative important size of solid particles transported in the pipe, the action of fluid on these particles (involving viscous shear stress or even important pressure gradient when the pipe is clogged) can be important. Consequently, the computation of the shear stress τ h from Equation (41) can be seen as a rather good estimation of the shear stress on the pipe walls, in a plane parallel to the main flow direction. However, this shear stress is only a partial view of the action of fluid on solid particles constituting the pipe surface (shear stress in direction perpendicular to the main flow direction can occur). On the other hand, the estimation of the shear stress τ h through Equation ( 39) is probably overestimated because it assumes that all the action of fluid on solid takes place on pipe walls along the main flow direction and discards transported particles. To resume, the estimation of a typical value of hydraulic shear stress applied on pipe walls is not straightforward (at least in our case), and computation of τ h and τ h probably give only boundary values.

CONCLUSION

In this article, we present a full two-way coupling scheme to deal with fluid-solid interaction modeling. Principles of the DE method and of the LB method are presented, and a complete description of the coupling implementation is provided. The coupled LB-DE method is very powerful and can be used in a wide range of engineering applications, especially when simulations of fluid-solid interactions at the grain scale are required. We have presented and used the coupling scheme only in two dimensions. Nevertheless, it can be easily extended to three dimensions for both the LB and DE parts.

To illustrate the method, we have performed simulations of piping erosion inspired by HET experiments, where both detachment and transport of particles can be modeled. The erosion law deduced from these simulations seems to be in agreement with the law verified and used experimentally. We found that the rate of erosion per unit area is proportional to the shear stress. Nevertheless, complementary simulations are required to accurately characterize dependency of erosion coefficient and critical shear stress on grain assembly properties. Such simulations can also be useful to identify relevant parameters responsible of particle detachment from the granular assembly.

In our knowledge, this simple example is one of the first attempts to apply the LB-DE method to such dense cases. Indeed, coupling method using LB method to simulate particle-fluid systems was used mainly to study suspensions or dilute systems contrary to systems encountered in geomechanics. This article and its presented application example highlight the abilities and efficiency of the method, which proves to be a promising tool in soil mechanics and more generally in civil engineering.

-case g) -known distributions: f 0 , f 2 , f 3 , f 6 -unknown distributions:

Figure 1 .

 1 Figure 1. Discrete velocity directions in the D2Q9 model.

Figure 3 .

 3 Figure 3. Velocity profiles of water flow around a circular obstacle. The evolution of velocity is presented along two perpendicular axes passing through the disk center (see insert). Three LB simulation results are presented for different values of d/h, where d is the diameter of the solid disk. Velocity profiles obtained with the commercial software FLUENT are also represented.

Figure 4 .Figure 5 .

 45 Figure 4. A boundary link σ. The solid boundary seen by the LB method and the real solid boundary are represented respectively with solid lines and dashed curve.

Figure 6 .

 6 Figure 6. Cohesive contact law used in the DE method.

10 L = l/ 3 lh = 8l/ 10 Figure 7 .

 103107 Figure 7. Schematic drawing of the numerical model. Hatched regions correspond to the locations of the cohesive material after drilling of the granular assembly.

Figure 8 .

 8 Figure 8. Series of snapshot from a simulation of piping erosion with the coupled LB-DE method for i = 12.7 × 10 -4 . Color gradient is proportional to fluid velocity.

Figure 9 .

 9 Figure 9. (a) Total hydrodynamic force acting on solid obstacles. The two components (along x and y) are represented. Dashed line represents the product of the pressure gradient ∆P = 0.03 with the channel height L = 0.01. (b) Total mass inside the channel during a simulation with i = 1.3 × 10 -4 .

4 Figure 10 .

 410 Figure 10. Evolution of the ratio Me/M 0 with time for different values of the hydraulic gradient i.

4 Figure 12 .

 412 Figure 12. Evolution of erosion rate ˙ with shear stress τ h for four simulations with different hydraulic gradients.

  f1 = f3 + 2/3 ρv x f4 = f2 -2/3 ρv y f5 = 1/2 [ρ(1 -2/3 vx + vy ) -f0 -2( f2 + f3 + f6 )] f7 = 1/2 [ρ(1vx + 2/3 vy ) -f0 -2( f2 + f3 + f6 )] f8 = f6 + 1/6 ρ(v xvy ) -case h) -known distributions: f 0 , f 1 , f 4 , f 8 -unknown distributions: f2 = f4 + 2/3 ρv y f3 = f1 -2/3 ρv x f5 = 1/2 [ρ(1 + vx -2/3 vy ) -f0 -2( f1 + f4 + f8 )] f6 = f8 -1/6 ρ(v xvy ) f7 = 1/2 [ρ(1 + 2/3 vxvy ) -f0 -2( f1 + f4 + f8 )]

  Indeed, generally, the DE time step dt DE (fixed with respect to a critical time step, see Section 3) is smaller than the LB time step dt. Consequently, to ensure the correspondence between the two physical times simulated by each method respectively, an integer number of DE loops should be included in a single LB cycle. Because it is easier to adjust dt DE than dt (depending on the viscosity, the relaxation time and the space discretization, through Equation (6)), the DE time step dt DE is reduced to a new value dt DE such as

  to compute τ h , is the pressure Figure11. Evolution of erosion rate ˙ with shear stress τ h for four simulations with different hydraulic gradients.

	(kg s -1 m -2 )	0.4
	. erosion rate per unit area ε	0.1 0.2 0.3							i= 8.3e-4 i=12.5e-4 i=16.6e-4 i=20.8e-4
		0 0	0.01	0.02	0.03	0.04	0.05	0.06	0.07
				shear stress τ	h (N m -2 )

† In recent research works, the experimental setup has been improved to achieve a direct measurement of the erosion rate[START_REF] Pham | Erosion et dispersion des sols argileux par un fluide[END_REF].

ACKNOWLEDGEMENTS

The scientific context and the financial support provided by the French Region Pays de la Loire through the EMERMOD project are gratefully acknowledged.

Contract/grant sponsor: Region Pays de la Loire -EMERMOD project

APPENDIX: BOUNDARY CONDITION RELATIONSHIPS BETWEEN DISTRIBUTION FUNCTIONS

This appendix presents relations used to determine the unknown distribution functions after the propagations. As explained in Section 2.2.2 this set of relations is used to impose a specific boundary condition with the D2Q9 model.

-case a) -known distributions: