Mouloud Mansouri

Moulay Said

El Youssoufi

François Nicot

Numerical simulation of the quicksand phenomenon by a 3D coupled Discrete Element -Lattice Boltzmann hydromechanical model

Keywords: granular soils, discrete element method, mollecular dynamics, lattice Boltzmann, fluid-solid interaction, quicksand

This paper deals with the numerical simulation of the quicksand phenomenon using a coupled Discrete Elements -Lattice Boltzmann hydromechanical model. After the presentation of the developed numerical model, simulations of ascending fluid flow through granular deposits are performed. The simulations show that the quicksand actually triggers for a hydraulic gradient very close to the critical hydraulic gradient calculated from the global analysis of classical soil mechanics, that is, when the resultant of the applied external pressure balances submerged weight of the deposit. Moreover, they point out that the quicksand phenomenon does not occur only for hydraulic gradients above the critical hydraulic gradient, but also in some cases with slightly lower gradients. In such cases, a more permeable zone is first gradually built at the bottom of the deposit through a grain rearrangement, which increases the hydraulic gradient in the upper zones and triggers the phenomenon.

INTRODUCTION

When a sand deposit is submitted to an upward hydraulic gradient, the grains undergo hydrodynamic forces essentially oriented in the direction of the applied gradient. These forces grow with increasing hydraulic gradient and may lead to the destabilization of the deposit by carrying the grains upward. This is the so-called quicksand phenomenon. Classical soil mechanics have been interested in this phenomenon in order to estimate the critical hydraulic gradient which is the minimum gradient that triggers quicksand. This critical gradient is thus estimated by expressing the balance between the resultant of external pressures and the submerged weight of the deposit.

The numerical simulation could bring more understanding on the phenomenon as it allows the visualization of its development for different hydraulic gradients. The most appropriate methods in this context are those based on the discrete element method (DEM) due to the discrete nature of the sand material. This method allows the modeling of the mechanical behavior of granular materials through the consideration of external forces acting on the grains and mutual granular interactions. For dynamic phenomena in saturated sands as the case of quicksand, the main external forces are obviously hydrodynamic forces due to fluid flows. Therefore, the DEM should be coupled with other appropriate methods for the evaluation of such forces. Within this context, several numerical approaches have been developed for the simulation of phenomena relative to saturated granular materials. The proposed methods may be classified in two essential categories; the first category includes methods in which the cells resulting from the discretization of the fluid domain are much larger than the grain size, the information within each cell is treated as local average and used to estimate the hydrodynamic forces [START_REF] Tsuji | Discrete particle simulation of two-dimensional fluidized bed[END_REF][START_REF] Nakase | Simulation study on liquefaction using DEM[END_REF][START_REF] Zeghal | Liquefaction of saturated loose and cemented granular soils[END_REF]. Although these methods are developed in order to reduce the computational cost, they may fail on description details of intergranular flows, which can have significant effects on the overall behavior. The methods of the second category are based on the detailed discretization of intergranular spaces and are usually named DNS-DEM approaches. These methods use commonly the Navier-Stokes equations discretized through finite element or finite volume techniques [START_REF] Sharma | A fast computation technique for the direct numerical simulation of rigid particulate flows[END_REF][START_REF] Wachs | A DEM-DLM/FD method for direct numerical simulation of particulate ows : Sedimentation of polygonal isometric particles in a Newtonian uid with collisions[END_REF][START_REF] Catalano | Pore-scale modeling of fluid-particles interaction and emerging poromechanical effects[END_REF]. Simulations based on these methods seem to be very good, however, 3D modeling using such discretizations remains computation demanding and is therefore very expensive.

In this work we use the DEM coupled with the Lattice Boltzmann method (LBM) for the fluid's modeling. This method is based on a detailed discretization of intergranular spaces and seems to be less expensive. Moreover, it has shown considerable capacity in the modeling of flows through porous media such as soils [START_REF] Mansouri | Numerical model for the computation of permeability of a cemented granular material[END_REF]. Although the use of the LBM is convenient for flows with straight boundaries, curved and moving boundaries that are usually involved in granular materials modeling, require special treatments. Therefore, different models are proposed to simulate particle-fluid interactions. The earliest models [START_REF] Ladd | Numerical simulation of particular suspensions via a discretized Boltzmann equation, Part 2, Numerical results[END_REF][START_REF] Aidun | Lattice-boltzmann simulation of solid particles suspended in fluid[END_REF], are based on the 'Standard Bounce Back' (SBB) rule to describe the non-slip boundary condition at the fluid-solid interface. The momentum exchange technique is then used to compute the interaction fluid-solid force. By using the standard bounce back rule the solid boundary is enforced to be located at the midpoints of the lattice links, which are cut by the solid particle boundary. Therefore, for a moving solid the boundary does not move continuously and smoothly in space; instead it will jump from one midpoint to another causing fluctuations on the computed interaction forces. Despite this shortcoming, this method has been successfully applied to the analysis of internal erosion in granular media [START_REF] Lominé | Modeling of fluid-solid interaction in granular media with coupled lattice Boltzmann/discrete element methods: Application to piping erosion[END_REF][START_REF] Sibille | Internal erosion in granular media: direct numerical simulations and energy interpretation[END_REF].

Another approach to consider fluid-solid interaction was proposed by Noble and Torczynski [START_REF] Noble | A lattice Boltzmann method for partially saturated computational cells[END_REF] and based on a modified collision operator for the fluid nodes that are partially covered by the solid. This method involves the fluid filling ratio of the voxel surrounding the fluid boundary node which is computed through the partial voxel's volume covered by the solid. This method has been widely used to model moving solids in a fluid [START_REF] Cook | A direct simulation method for particle-fluid systems[END_REF][START_REF] Feng | The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problems[END_REF][START_REF] Strack | Three-dimensional immersed boundary conditions for moving solids in the Lattice Boltzmann method[END_REF]. It was also used in the coupled DEM-LBM model proposed by Feng et al. [START_REF] Feng | Coupled Lattice Boltzmann and Discrete Element modelling of fluid-particle interaction problems[END_REF].

Following an alternative approach, other researchers have proposed the Interpolated Bounce Back (IBB) conditions [START_REF] Filippova | Grid refinement for lattice-BGK models[END_REF][START_REF] Mei | An accurate curved boundary treatment in the lattice Boltzmann method[END_REF][START_REF] Yu | Viscous Flow Computations with the Lattice Boltzmann Equation Method[END_REF][START_REF] Bouzidi | Momentum transfer of a Boltzmann-lattice fluid with boundaries[END_REF][START_REF] Kao | An investigation into curved and moving boundary treatments in the lattice Boltzmann method[END_REF]. Such conditions use the parts covered by the solid of the lattice links cut by the boundary. Because they are based on a computation of lengths (parts of links covered by the solid) rather than volumes as in the method of Noble and Torczynski [START_REF] Noble | A lattice Boltzmann method for partially saturated computational cells[END_REF], these conditions are easier and more effective in practice. Using the IBB condition proposed by Yu et al. [START_REF] Yu | Viscous Flow Computations with the Lattice Boltzmann Equation Method[END_REF], a coupled DEM-LBM 2D model have been developed to analyze the surface erosion in a particle bed [START_REF] Abdelhamid | Pore-scale modeling of surface erosion in a particle bed[END_REF].

In the present work, we have develop a 3D coupled LBM-DEM model based on the IBB rule proposed by Bouzidi et al. [START_REF] Bouzidi | Momentum transfer of a Boltzmann-lattice fluid with boundaries[END_REF]. The model was first validated and then used to simulate the quicksand phenomenon.

In the remainder of this paper, the methods DEM and LBM are first described, then the procedure of their coupling is presented together with some technical aspects for the implementation. After the validation part, simulations of the quicksand phenomenon are presented. The paper closes by discussing the different obtained results and drawing some concluding remarks.

A DISCRETE ELEMENT METHOD

The Discrete Element approach models the grains of a granular material by independents elements, each element interacts with its neighbors at the contact points. The overall deformation of the medium is mainly due to the relative movements of the grains assumed as rigid bodies. Therefore, the behavior of the medium can be described through the integration of the dynamic equations applied to each element. Such equations are written by taking all external forces into account, such as the contact forces, the gravitational forces and the hydrodynamic forces. Because these forces can change abruptly with time, the integration should be performed in an incremental manner using small enough time steps.

The dynamic of a grain i is governed by Newton's second law, consequently, the equations of translational and rotational motions are ´mi E R

x i D P j E F contact ij C E F hyd i C m i E g I i E R ' i D P j E M contact ij C E M hyd i (1)
where m i and I i are the mass and the mass moment of inertia of the grain, E R x i and E R ' i are accelerations of translation and rotation, E F contact ij is the interaction force applied by a grain j on the grain i, E F hyd i is the hydrodynamic force applied on the grain i and E g is the acceleration of gravity. E M contact ij and E M hyd i are the torques caused by the contact forces and the hydrodynamic force on the grain i, respectively.

The molecular dynamics method originally proposed by Cundall and Strack [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] is implemented in this work. This method allows a slight overlap of the grains which is used to calculate the contact forces through a contact low. The grains are assumed spherical in 3D and circular in 2D modeling. At each time step, the contacts are first detected by examining the normal distance between grains D n (Figure 1). Two grains are in contact if they overlap (D n Ä 0).

A contact force applied by a grain j on a grain i can be decomposed into normal and tangential components. The normal force is calculated here using the viscoelastic linear model (Figure 2a): where V n is the normal velocity of the grain j relative to the grain i, k n ,and n are the elastic and the viscous damping constants, respectively, and n E is the unit vector pointing from i to j :

E F n D . k n D n n V n /:E n (2)
E n D . E x j E x i /=kE x j E x i k (3)
in which E x i and E x j are the positions vectors of the grains. The tangential force is computed using the simple frictional Coulomb model (Figure 2b):

E F s D d F n E V s k E V s k (4)
where E

V s is the tangential velocity of the grain j relative to the grain i and d is the dynamic coefficient of friction.

It should be noted that the choice of the elastic and the viscous damping constants (k n and n ) is not arbitrary, in fact k n must be high enough to avoid substantial overlap that affects the overall behavior and n that controls the damping in the material, is selected such that the restitution coefficient n defined as the ratio of the normal velocities at the start and the end of the contact is between 0 and 1. Hence, it is usual in practice to select n and calculate n from the relation [START_REF] Pöschel | Computational Granular Dynamics -Models and Algorithms[END_REF]:

n D 2log. n / p k n m eff q 2 C .ln n / 2 (5)
where m eff D m i m j =.m i C m j /, m i and m j are the masses of the grains in contact.

Once the contact forces and the other external forces acting on a grain are obtained at a discrete time t , accelerations are calculated through the dynamic equations (Eqs. ( 1)), then integrated to obtain the kinematic variables of the grain at time .t C t/. Considering that accelerations are not constant during a time increment t, special algorithms are used to predict accurately the kinematic variables, such as the Velocity-Verlet algorithm [START_REF] Allen | Computer simulation of liquids[END_REF].

Finally, it is important to note that in order to correctly describe the evolution of the contact force, the time step t should be sufficiently small compared with the contact duration t c . It is usually taken as t t c =10. For the viscoelastic model the contact duration can be approximated by t c D p m eff =k n , so the maximum time step is taken as:

t max 0:1 p m=k n ( 6
)
where m is the smallest effective mass in the system.

LATTICE BOLTZMANN METHOD FOR FLUID FLOWS MODELING

Standard formulation

In the LBM, one solves the kinetic equation for the particle distribution function f .E x; E ; t/, which depends on the spatial position E x, the velocity of particules E , at the time t . The macroscopic quantities of interest such as mass density and momentum density E u are weighted averages of the distribution function:

D Z f d E .a/ E u D Z E f d E .b/ (7)
A popular kinetic model adopted for the method is the so-called BGK (Bhatnagar, Gross, and Krook) model. In this model the collisions term in the Boltzmann equation is simplified using the simple relaxation time approximation:

@f @t C E : E rf D 1 f f .0/ Á ( 8
)
where f .0/ is the equilibrium distribution function (Maxwell-Boltzmann equilibrium function) and is the relaxation time. To solve for f numerically, Eq. ( 8) is discretized twice: a first discretization with respect to the time involving a time increment t and a second discretization based on the velocity space by choosing a finite set of velocity vectors that particles can have. The continuous particle distribution function f.Ex; E ; t/ becomes therefore, a set of discrete distributions f i .Ex; t/ associated with the chosen velocity vectors e E i .

As an illustration, we present hereafter the D2Q9 model (2 Dimensions, 9 Velocity vectors) which is widely used for 2D simulations. The detailed model's derivation can be found in [START_REF] He | Theory of the lattice Boltzmann method: From the Boltzmann equation to Lattice Boltzmann equation[END_REF]. As sketched in Figure 3a The discretizations lead to the LBGK (Lattice BGK) equation, that describes the incremental evolution of the discrete particle distributions f i :

f i . E x C E e i t; t C t/ f i . E x; t/ D 1 f i . E x; t/ f .eq/ i . E x; t/ Á ( 9
)
where f .eq/ i is the discrete equilibrium distribution given as:

f .eq/ i D w i Ä 1 C 3 c 2 E e i E :u C 9 2c 4 E e i E :u 2 3 2c 2 E u E :u (10)
with w i being the weighting factors; w 0 D 4=9, w 2;4;6;8 D 1=9, w 1;3;5;7 D 1=36, is the dimensionless relaxation time such as:

1 D t (11)
Thus, the mass and the momentum densities (Eqs. ( 7)) can be rewritten as:

D Q 1 X iD0 f i .a/ E u D Q 1 X iD0 E e i f i .b/ ( 12
)
where Q is the number of the chosen discrete velocity vectors. Equation ( 9) describes the incremental evolution of the discrete particle distributions f i at the nodes of a regular lattice having a space step x D ct (Figure 3b). Hence, if the discretization's step of the space is selected, the characteristic speed of the model is defined by: The discrete density distributions f i are calculated at each time step according to Eq. ( 9) as follows

c D x t (13)
collision stepW f out i . E x; t/ D f i . E x; t/ 1 f i . E x; t/ f .eq/ i . E x; t/ Á streaming stepW f i . E x C E e i t; t C t/ D f out i . E x; t/
where f out i represents the post-collision density distribution. After the streaming step, the discrete distributions f i are obtained at each node, consequently, the mass and momentum densities can be calculated at each node too. The fluid pressure p can be computed from the mass density through the equation of state p D c 2 s , where c s is the sound speed of the model related to the lattice model's speed c as c s D c= p 3.

Discretization parameters for incompressible fluid flow simulations

It can be shown through the Chapman-Enskog analysis [START_REF] Luo | Theory of the lattice Boltzmann method: Lattice Boltzmann models for nonideal gases[END_REF] that LB models recover the incompressible Navier-Stokes equations when the density fluctuation of the fluid is assumed to be negligible and the equivalent kinematic viscosity is given by:

D 1 3 cx  1 2 à (14)
Therefore, in order to correctly simulate an incompressible fluid flow, one must ensure that the density fluctuation is sufficiently small. This can be achieved using a model in which the sound speed c s is larger enough than the maximum velocity of the simulated flow u max , that is, when the 'computational' Mach number defined as M a D u max c s is sufficiently small. In practice, M a should be maintained, smaller than 0.1. The discretization parameters are x, t, and , if the viscosity of the fluid is given, only two of these parameters can be chosen independently because they are related through Eq. [START_REF] Feng | The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problems[END_REF]. In practice, it is often convenient to choose and x as two independent parameters and t is derived from Eq. ( 14). This is due to the fact that is largely responsible for the numerical stability of LB simulations and x is often dictated by the space description precision. The BGK-LBM is theoritically convergent for 0:5 < < 1, in fact, this range of corresponds to positive viscosities (Eq. ( 14)). In practice is typically chosen in the range 0:5 ; 3.

Solid moving boundary treatment

Boundary conditions are introduced in LBM in terms of distribution functions f i that are constructed from the imposed physical boundary conditions such as pressure and velocity. We present in the following the treatment of a solid moving boundary, because it is the principal one encountered in granular materials hydromechanics. Details about other boundary conditions in the LBM can be found in the literature [START_REF] Yu | Viscous flow computations with the method of lattice Boltzmann equation[END_REF]. A simple way to represent solid obstacles such as solid particles in a LB dicretization, is to consider them as composed of grouped sets of pixels (voxels in 3D) whose centers are lattice nodes (Figure 4). We denote 'fluid nodes' the nodes that are in the fluid domain and 'solid nodes' the ones covered by solid obstacles.

A solid boundary node is a solid node having at least one link with fluid nodes. The fluid-solid interaction is taken into account through the boundary nodes. The intenal solid nodes are excluded from the computation process. It is assumed that there is no slip of the fluid along a solid surface, that is, the fluid particles in contact with a solid surface have the velocity of this surface at the contact point. This condition could be imposed in LBM using the bounce back scheme. This means that any incoming fluid particle from a fluid node to a solid node will be reflected back to the node it comes from. The momentum of the reflected particle may be different from the momentum of the streamed one, it depends on the position and the velocity of the solid boundary. Therefore, at the solid boundary nodes, instead of the collision step, the distributions that will be used in the streaming step are constructed using the post-collision distributions at neighboring fluid nodes. In this work, we use the IBB scheme proposed by Bouzidi et al. [START_REF] Bouzidi | Momentum transfer of a Boltzmann-lattice fluid with boundaries[END_REF]. For a linear interpolation, the post-collision distributions to be assigned to the solid boundary nodes before the streaming step are 8 < :

f out O i . E x s ; t/ D 2qf out i . E x f ; t/ C .1 2q/f out i . E x f 0 ; t/ C 6w i w E e O i : E u w c 2 ; q 6 1 2 f out O i . E x s ; t/ D 1 2q f out i . E x f ; t/ C 2q 1 2q f out O i . E x f 0 ; t/ C 3 q w i w E e O i : E u w c 2 ; q > 1 2 (15)
where E e O i refers to the opposite direction of the direction E e i (E e O i D E e i ), E u w is the velocity of the boundary, w is the fluid mass density at the boundary and q defines the fraction of the boundary intersected link located in the fluid domain and calculated by referring to Figure 4 as:

q D k E x f E x w k k E x f E x s k (16)

Approximation of the distributions on the nodes regaining the fluid domain due to a solid boundary movement

When a solid boundary moves, there are grid nodes that move from the solid region into the fluid region to become fluid nodes (indicated by in Figure 4). Therefore, one must specify the distribution functions applying to these nodes. There are several techniques to compute values of the unknown distribution functions in this context, including the extrapolation from neighboring nodes, or by approximation as equilibrium distributions, using the velocity of the solid at the considered node just before it leaves the solid region with an averaged mass density. Or alternatively, by systematically updating the distribution functions in the solid nodes by performing collisions as in the fluid nodes. It seems that all these schemes produce similar results [START_REF] Lallemand | Lattice Boltzmann method for moving boundaries[END_REF]. In this work, unknown distributions are approximated as the equilibrium distributions computed using the averaged fluid's density in the whole system and the velocity of the solid particle at the specified node's position just before it leaves the solid region.

Hydrodynamic force on a grain, momentum exchange method

Based on the bounce back scheme, the hydrodynamic forces acting on an obstacle can be obtained through the momentum exchange method. The momentum change of the bounced fluid particle is OEf out

O i . E x s ; t/E e O i f out i . E
x f ; t/E e i x 3 . Then, the average force transmitted over a time step t to the solid particle, along a boundary link `defined by the boundary nodes E

x f and E x s (Figure 4) is:

E F `D x 3 t h f out O i . E x s ; t/E e O i f out i . E x f ; t/E e i i (17)
The total hydrodynamic force exerted on the grain can be calculated by summing up the forces from all the related boundary links:

E F hyd D X `E F `( 18
)
The force E F `induces a torque around the grain's center of mass calculated as:

E M `D E r ` E F `( 19
)
where E r

`D E x f C E x s 2 E x G .
Likewise, the total hydrodynamic torque is obtained by summing over all the boundary links:

E M hyd D X `E M `(20)

DEM AND LBM COUPLING

Coupling procedure

The coupling of the two methods involves the use of the DEM by introducing the hydrodynamic forces as external forces in the dynamic equations of the grains (Eqs. ( 1)). These forces are obviously dependent of grains motion, therefore, they should be re-evaluated continuously. The space LB discretization is dictated in this case by the diameters of the granular material particles. That is, solid particles should be discretized with a sufficient resolution. For moderate Reynolds numbers, Yu et al. [START_REF] Yu | Viscous flow computations with the method of lattice Boltzmann equation[END_REF] have shown that in order to obtain an accurate evaluation of the hydrodynamic force on a grain, the discretization's step (x) must be lower than about one-tenth (1=10) of the grain diameter in the 2D modeling (circular grain) and about one-seventh (1=7) of the grain diameter in the 3D case (spherical grain).

As mentioned in the previous section, the time step in LBM (noted in the following t LB ) depends on the other discretization parameters and calculated from Eq. ( 14). t LB is often larger than the maximum value DEM time step noted here t DE max and calculated from Eq. ( 6). Therefore, one should perform a number n d of DEM computation steps then perform one LB computation step. This can be carried out by selecting the DEM time step t DE 6 t DE max , such that n d t DE D t LB , with the integer number n d may be computed as:

n d D I nt  t LB t DE max à C 1 (21)
then the DEM time step is set:

t DE D t LB n d (22)
It should noted finally that the model presented in this manuscript does not account for any gravity effect on the fluid. The effect of buoyancy could be adequately computed by re-scaling the specific weight of the grains without changing their inertial properties. This can be achieved by multiplying the gravitationnal acceleration applied to the grains by the coefficient .1 w = s /, where w and s are the fluid and solid densities, respectively.

Numerical inspection

The applications presented here are focused on the validation of the hydrodynamic forces computation procedure using 2D or 3D models. For 2D cases the D2Q9 Lattice Boltzmann model is used, whereas in 3D we employ D3Q19 model whose the discrete velocity vectors are (Figure 5):

where c is the characteristic speed of discretization. The discrete equilibrium distributions f .eq/ i have the same expressions as in the D2Q9 model (Eq. ( 10)) with the weight factors are: 6(a)) and the expressions of hydrodynamic force and torque applied on the sphere are given as follows [START_REF] Clague | The hydrodynamic force and torque on a bounded sphere in Poiseuille flow[END_REF]:

w i D 8
F hyd D 6 rU h " 1 1 9 r h 2 1 0:6526 r h C 0:3160 r h 3 0:2420 r h 4 # (25)
M hyd D 8 3 r 2 U h r h Á Ä 1 C 0:0758 r h Á C 0:0490 r h Á 2 (26)
where is the fluid's dynamic viscosity, r is the sphere's radius, U h is the fluid velocity at position h (indicated in Figure 6) for a Poiseuille flow (without the sphere), it is expressed:

U h D 1 2 hH h 2  dp dx à ( 27
)
where dp dx is the pressure gradient in the x direction. In this application, the flow is driven by introducing an acceleration a x in the x direction. The equivalent pressure gradient is then dp dx D f a x , f being the fluid density. The acceleration is introduced according to the procedure given [START_REF] Yu | Viscous flow computations with the method of lattice Boltzmann equation[END_REF] and selected to get the desired particulate Reynolds number defined as:

Re D U h d (28)
where d D 2r is the sphere's diameter and is the fluid's kinematic viscosity.

As the problem is infinite, it is convenient to implement periodic boundary conditions in x and y directions. The simulated problem is then a periodic cell whose dimensions are the periods in the two directions D x and D y (Figure 6(b)).

The effect of the discretization size on the computed forces is also examined in these simulations. The discretization size is given by the sphere's radius in LB units N r D r= x . Three series of simulations are then performed for N r D 5:50, N r D 3:52 and N r D 3:08. The geometrical parameters used in the simulation are set to r=H D 0:11, h=H D 0:25. Periods are the same in both directions and set D x D D y D D. The period D is increased gradually in order to show its effects on the computed force and torque. The dimensionless relaxation time is set in all simulations to D 0:7 and the Reynolds number is chosen equal to 0.0006, which is in the range of Stokes flow.

Figures (7 and 8) show the evolution of the computed normalized drag force and normalized torque when enlarging the cell size with respect to the sphere's radius. The normalized drag force and torque are defined as:

N F hyd D F hyd 6 rU h (29) N M hyd D M hyd 8 3 r 2 U h r h (30)
The dashed-line curves through the data are provided to guide the eye. It is clear that the computed drag force and torque increase with increasing distance between periodic spheres and tend to Wakiya's solution for large periods. The reduction of the discretization size has almost no effect on the computed force, however, for the computed torque, this reduction causes a slight overestimation. Neverthless, the resulting errors remain low for the tested sizes. This is consistent with the conclusion drawn by Mei et al. [START_REF] Mei | Lattice Boltzmann method for 3-D flows with curved boundary[END_REF] namely the computed force on a sphere is accurate when the radius is greater than 3.5 lattice spaces.

Hydrodynamic force on a moving sphere.

In the previous section, we have checked the validity of the model for a fixed sphere in a Poiseuille flow. This application, however, does not involve the movement of the solid across the lattice. Thus, in this section, we compute the hydrodynamic force on a moving spherical grain in a channel bounded by two parallel plates. The geometrical configuration of the problem is similar to that of the previous section, with the specifications r=H D 0:11, h=H D 0:5, D y =r D 6 and D x =r D 45. At time t D 0, the fluid in the channel is at rest and the sphere is impulsively started to move with a constant speed u s D U 0 in x-direction. The sphere's initial position is .x 0 ; y 0 ; ´0/ D .D x =10; D y =2; H=2/.

This problem should be equivalent to the following one: the channel with the fluid is moving at a uniform speed U x , at time t D 0, the sphere is placed at its position and maintained at rest, the speed of the lower and upper plates being maintained u W D U 0 .

Therefore, in the first problem the sphere is moving relatively to the lattice grid, whereas in the second one the fluid is flowing but the sphere is at a fixed position relatively to the lattice grid. The following dimensionless parameters are used for the simulations: the sphere's radius in LB units is N r D 5:5, the dimensionless imposed velocity is N U x D U x =c D 0:04 and the dimensionless relaxation time is D 0:63. The flow's characteristic Reynolds number is then R e D 10. Figure 9 shows the drag force acting on the sphere over time. The plot is presented in a dimensionless form where the dimensionless drag force is computed as N F hyd D F hyd =. 2 x c 2 / and the time is given in LB units, that is, the number of time steps performed. Two curves are shown, the first one is a fluctuating line (thin line) corresponding to the case of an impulsively started sphere movement at a constant velocity N u G D 0:04 with the channel being at rest ( N u W D 0), the second one is a smouth line (thick line), corresponding to the inverse case where the sphere is at rest ( N u G D 0) and the channel starts moving with the velocity N u W D 0:04. These curves show that the drag force variation is the same for the two cases. On the whole, the fluctuations observed on the computed drag force in the case of the moving sphere are small compared with the mean value of the force at any time. These fluctuations can be attributed to the approximation used in the computation of the distributions at the nodes joining the fluid domain due to the grain movement (as discussed in section 3.4).

Sedimentation of two particles 'drafting, kissing, and tumbling'.

The two validations performed earlier do not involve the effective DEM-LBM coupling. In this application, we simulate the sedimentation of two solid particles in a column of fluid. It is known that two particles dropped close to each other in a Newtonian fluid will undergo drafting, kissing, and tumbling as demonstrated experimentally by Fortes et al. [START_REF] Fortes | Nonlinear mechanics of fluidization of beds of spherical particles[END_REF]. This phenomenon has also been captured by numerical simulations using different methods in two and three dimensions [START_REF] Feng | The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problems[END_REF][START_REF] Strack | Three-dimensional immersed boundary conditions for moving solids in the Lattice Boltzmann method[END_REF][START_REF] Qi | Lattice-boltzmann simulations of particles in non-zero-Reynolds-number flows[END_REF][START_REF] Patankar | A new formulation of the distributed Lagrange multiplier/fictitious domain method for particulate flows[END_REF][START_REF] Patankar | A formulation for fast computations of rigid particulate flows[END_REF]. Because of its clarity in visualization, a 2D simulation is performed here. We replicate a problem solved by Patankar et al. [START_REF] Patankar | A new formulation of the distributed Lagrange multiplier/fictitious domain method for particulate flows[END_REF][START_REF] Patankar | A formulation for fast computations of rigid particulate flows[END_REF], using the finite method based on the domain method. This problem was resolved later by Feng and Michaelides [START_REF] Feng | The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problems[END_REF] in order to validate their immersed boundary-lattice Boltzmann method.

The channel is 2 cm wide (x-direction) and 8 cm high (y-direction). The fluid has the same properties as water with viscosity 0:01 g=cm:s and density 1 g=cm 3 . The particles density is 1:01 g=cm 3 , and the radii of the particles are 0:1 cm. Initially, the first particle is shifted with 0:001 cm from the channel's center and at a height of 7:2 cm, the second particle is at the channel center at a height of 6:8 cm. The two particles start settling in the y-direction, due to the gravity force. The relaxation time is set to D 0:75 and the computational domain is discretized to 200 800 lattice units, the radius of each particle is then 10 lattice units.

Figure 10 shows the positions of the particles at different times. Fixed straigth lines on the particles are represented to show their rotations during the sedimentation. It can be observed that the drafting, kissing, and tumbling are quite well captured by the LBM-DEM simulation. Qualitatively, the sedimentation process is almost the same as that observed by Feng and Michaelides [START_REF] Feng | The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problems[END_REF].

The longitudinal velocities during the sedimentation process are plotted in Figure 11 together with those obtained by Feng and Michaelides [START_REF] Feng | The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problems[END_REF] and Patnakar [START_REF] Patankar | A formulation for fast computations of rigid particulate flows[END_REF]. This graph shows a good agreement in the velocity evolution obtained by the three methods. It can be noted that the different phases of the phenomenon occur at nearly identical times for the three methods. From a quantitative point of view, the graph shows small differences between the results of the three methods during the two first phases, that is, the drafting and the kissing phases, but from the start of the tumbling phase (at time t 2:4 s on the graph) the difference becomes more significant. Feng and Michaelides [START_REF] Feng | The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problems[END_REF] in a comparison of their results with those of Patnakar [START_REF] Patankar | A formulation for fast computations of rigid particulate flows[END_REF], attributed the notable difference that arises from the start of the tumbling phase to the physical instability of the tumbling's initiation.

QUICKSAND PHENOMENON SIMULATIONS

Before analyzing the simulations, it should be noted that the presented model was previously used in a simplified form to estimate the critical hydraulic gradient corresponding to the quicksand onset [START_REF] Mansouri | A 3D DEM-LBM approach for the assessment of the quick condition for sands[END_REF][START_REF] Mansouri | Hydromécanique Des Milieux Granulaires Saturés Par Une Approche Micromécanique Couplée éléments Discrets / Lattice Boltzmann -Application à La Géotechnique[END_REF]. This gradient was estimated through the analysis of the evolution of the vertical stress in a deposit of grains subjected to an increasing hydraulic gradient. The obtained results were very close to the ones calculated through the global analysis of classical soil mechanics. In that work, the hydrodynamic forces acting on the grains were assumed to be proportional to the applied gradient. As a consequence, they are computed for a single value of the hydraulic gradient, then obtained by proportionality for any applied gradient. For this purpose, no relative movement between the grains is allowed until the lifting of the deposit occurs. Therefore, it was not possible to visualize the deposit after the boiling onset.

To better understand the micromechanical mechanisms that occur in the quicksand phenomenon, we present the simulations of grain deposits subjected to ascending water flows under hydraulic gradients near the critical hydraulic gradient.

3D quicksand simulations

As the phenomenon is 3D, the first simulations are performed through a 3D modeling. We construct polydisperse granular deposits with periodic boundaries in both horizontal directions using the DEM. The deposits are composed of spherical grains having diameters that follow the cumulative beta distribution [START_REF] Voivret | Space-filling properties of polydisperse granular media[END_REF], with maximum and minimum diameters of 0.4 and 0:8 mm, respectively. Once a grain deposit is built, it is placed into a water column of the same section and of larger height, then subjected to the hydraulic gradient by imposing a pressure difference between the horizontal lower and upper boundaries. For this purpose, we employ the pressure boundary conditions developed by Zou & He [START_REF] Zou | On pressure and velocity boundary conditions for the Lattice Boltzmann BGK model[END_REF] and based on the assumption the bounce back of non-equilibrium distribution. Note that the details of these conditions for the 3D case are not given in this reference, their implementation procedure for the D3Q19 model is therefore presented in Appendix. On both horizontal directions, periodic conditions are implemented on the water flow.

The unit weight and kinematic viscosity of water are taken as w D 10 kN=m 3 and D 10 6 m 2 , respectively, while the unit weight of solid grains s is set to 26 kN=m 3 .

Starting from its classical expression, the hydraulic gradient could be written as:

i D 0 w pS W 0 ( 31
)
where 0 is the submerged unit weight of the soil, p is the imposed pressure drop between the lower and upper horizontal surfaces of the deposit, S and W 0 are the horizontal section and the submerged weight of the deposit, respectively. The product pS represents the resultant of external pressures applied to the deposit. With increasing p, the quicksand phenomenon triggers when the resultant of external pressures balances the submerged weight of the deposit, that is, pS D W 0 , therefore, the critical hydraulic gradient is i c D 0 w and the applied hydraulic gradient could be written in terms of the critical hydraulic gradient as

i D i c pS W 0 (32)
As S and W 0 are fixed (and characterize the deposit), it results from Eq. ( 32) that in order to impose a fraction of the critical hydraulic gradient, the only parameter to select is the pressure drop p. It should be noted that W 0 is calculated as the sum of the submerged grain weights.

Hereafter an example of a 3D quicksand simulation is described. The deposit consists of 100 grains, periodically disposed in both horizontal directions with a period of 2:4 mm. The discretization steps in the horizontal directions are 40x 40x, so that the diameter of the smallest grain is discretized in 6.67 space steps (x).

In order to analyze the deposit's behavior when varying the applied hydraulic gradient, the deposit is subjected to different gradients near the critical hydraulic gradient. It is observed that: For i i c no grain movement is observed, the deposit remains in its initial state. For i slightly less than i c , a relaxation occurs in some cases starting from the bottom of the deposit. Then the grains undergo a carrying upwards. Figure 12 shows this behavior for an applied gradient i 0:92i c . For i > i c , the deposit undergoes an upward movement with its preserved form as soon as the gradient is applied, but it loosens from the bottom during the movement. Moreover, it is noted that the movement is faster when the gradient is increased.

It is worth noting that for this deposit, the hydraulic gradient i 0:92i c is the minimum gradient that triggers the quicksand. These observations indicate that quicksand can occur in cases where i < i c . Similar results are obtained for other samples with larger numbers of grains, but the minimum gradient below the critical gradient which causes the quicksand is variable, depending on the initial compactness of the deposit and on the applied pressure rate. From a theoretical point of view, the quicksand is expected for i > i c , but not for i < i c . Therefore, further investigations are required. For this purpose, it is convenient to visualize the movement of deposit grains, nevertheless the visualization is difficult for 3D simulations. This is why 2D simulations of quicksand were carried out, making it easier to visualize the grain movements.

2D quicksand simulations

First and foremost, it is worth noting that 2D discrete element modeling of a granular deposit gives rise to a porous medium with non interconnected pores. This completely changes the hydraulic properties of a saturated granular material. Boutt et al. [START_REF] Boutt | Coupled fluid-solid model for problems in geomechanics: Application to sand production[END_REF], in order to simulate hydromechanics of saturated granular materials with a 2D model, have used a reduced grain diameter in LB modeling while keeping the actual diameter in DEM modeling. The model simulates therefore a connected pore network. In addition, they proposed an approximate expression of the permeability including the average radius and a radius reduction multiplier, in order to extrapolate the 2D simulation results to 3D problems with comparable permeability and average grain size.

An inverse approach is proposed in this manuscript. The diameters of the grains are increased to compute the contact forces while conserving the actual diameters in the calculation of the grain weights and drag forces. This method is mainly adopted in order to calculate the drag force with the actual diameter. By doing so, the contact forces between the grains develop before the contacts actually occur and a connected pore network is created. In this study, the distance added to the actual radius in the discrete elements modeling is set equal to one spatial step in the Lattice-Boltzmann discretization .x/, consequently, the space between each two grains contains at least one node of the Lattice-Boltzmann grid. The simulated deposit is composed of 300 grains, periodic in the horizontal direction with a period length of 14:4 mm. The resulting average height is about 8:8 mm and the total submerged weight is 1:265 N . The D2Q9 Lattice Boltzmann model is used, the discretization in the horizontal direction is 400x, so that the smallest grain diameter is discretized into about 11x.

The deposit is first subjected to a vertical flow under a pressure difference equal to 100N=m 2 applied slowly. The resultant force of the applied pressure is then of 1:44 N (pS ) and the corresponding hydraulic gradient is i 0:90i c . The applied pressure was maintained for a long time but nothing happened, the deposit has kept its initial configuration. It was verified that when the steady state flow is established, the resultant of the computed hydrodynamic forces is very close to the resultant of applied pressure (1:44 N ), it is only about 0:018% larger.

After that, the applied gradient was gradually increased, but nothing happened until i 0; 97i c . Figures 13,14, and 15 show the deposit evolution over time for three applied hydraulic gradients i D 1; 1i c , i D i c and i 0; 97i c .

The time increment used to plot these figures (ıt D 0:065s) was selected to compare conveniently the three cases. Figures 13 and14 corresponding to the gradients i D 1:1i c and i D i c , respectively, show that quicksand occurs in a similar manner for gradients i > i c , in fact the deposit is fully raised initially, but for a low gradient (i D i c ) it loosens quickly while for a high gradient it may be raised like a shutter. This can be attributed to the high pressure applied at the bottom of the deposit which prevents loosening. Furthermore, it can be observed that the lifting is more quicker for the high gradient, because the grains reach the end of the column of water in a shorter time; 3ıt for i D 1:1i c against 4:7ıt for i D i c .

For the case i 0:97i c , Figure 15 shows that quicksand develops step by step; first, there is a grain rearrangement in some areas at the bottom of the deposit that allows the creation of large channels where there is no pressure drop (surrounded areas in Figure 15). Accordingly, the hydraulic gradient in the area overlying the channel becomes larger. This increase in gradient initiates an uprising of this zone which leads to the loosening of the lower zone. This loosening grows progressively until the quicksand onset. It should be noted that in this case (i 0:97i c ) the quicksand process is much slower than in the case with i > i c , and the initial stage of grain rearrangement is the most time consuming stage. These results confirm the previous results obtained in the 3D simulations. In addition, they bring some comprehensive elements on the quicksand triggering for cases where i<i c .Infact, for loose deposits, the water flows induce a grain rearrangement at the bottom of the deposit which leads to an increased hydraulic gradient at the top causing the onset of quicksand.

CONCLUDING REMARKS

We have presented a 3D coupled Lattice-Boltzmann-Discrete element hydromechanical model to investigate the micro-mechanical dynamic response of saturated granular materials. The developed model was first checked and validated from different heuristic applications, including the computation of the hydrodynamic forces on fixed and moving grains and the simulation of the 'Drafing-Kissing-Tumbling' phenomenon. The obtained results are very close to those reported in the literature including numerical and experimental data. Then, the model was used to simulate the quicksand phenomenon. The performed simulations show that the quicksand actually triggers for a hydraulic gradient very close to the critical hydraulic gradient calculated from the standard soil mechanics analysis, that is, when the resultant of the applied external pressure balances submerged weight of the deposit. Moreover, they reveal that this phenomenon does not occur only for hydraulic gradients above the critical hydraulic gradient, but also at lower gradients when a more permeable zone can develop at the bottom of the deposit due to a grain rearrangement.

Finally, as in the presented quicksand simulations the developed model could be applied in simulation of more complex phenomena related to saturated granular media, such as sand liquefaction and internal erosion of soils. This could shed light on the micro-mechanical mechanisms involved. After the streaming step, the unknown distributions are f 7 , f 9 , f 12 , f 15 and f 17 . Assuming that their non-equilibrium parts are equal to that in the opposite directions (i.e., f i f .eq/ i D f O i f .eq/ O i

), therefore they can be expressed as:

Figure 1 .

 1 Figure 1. Normal distance between two grains.

Figure 2 .

 2 Figure 2. Modeling of the contact interaction forces; (a) the normal force model, (b) the tangential force model.

  , the discrete velocity vectors are chosen as follows : e E 0 D c.0; 0/; e E 1 D c. 1; 1/; e E 2 D c. 1; 0/; e E 3 D c. 1; 1/; e E 4 D c.0; 1/; e E 5 D c.1; 1/; e E 6 D c.1; 0/; e E 7 D c.1; 1/; e E 8 D c.0; 1/, where c is a characteristic speed of the model.

Figure 3 .

 3 Figure 3. (a) D2Q9 model, (b) Flow domain discretization using the D2Q9 model.

Figure 4 .

 4 Figure 4. Representation of a moving grain on the lattice grid.

Figure 5 .

 5 Figure 5. D3Q19 Lattice Boltzmann model, each arrow corresponds to a velocity between the central node and one of the 18 neighboring nodes.

E e 1 D

 1 c. 1; 0; 0/; E e 2 D c.0; 1; 0/; E e 3 D c.0; 0 1; / E e 4 D c. 1; 1; 0/; E e 5 D c. 1; 1; 0/; E e 6 D c. 1; 0; 1/ E e 7 D c. 1; 0; 1/; E e 8 D c.0; 1; 1/; E e 9 D c.0; 1; 1/ E e 10 D c.1; 0; 0/; E e 11 D c.0; 1; 0/; E e 12 D c.0; 0; 1/ E e 13 D c.1; 1; 0/; E e 14 D c.1; 1; 0/; E e 15 D c.1; 0; 1/ E e 16 D c.1; 0; 1/; E e 17 D c.0; 1; 1/; E e 18 D c.0; 1; 1/

Figure 6 .

 6 Figure 6. (a) Fixed spherical grain in a 3D Poiseuille flow, (b) Periodic cells used in the modeling.

Figure 7 .

 7 Figure 7. Hydrodynamic force vs. cell's size.

Figure 8 .

 8 Figure 8. Hydrodynamic torque vs. cell's size.

Figure 9 .

 9 Figure 9. Drag force vs. time.

Figure 10 .

 10 Figure 10. Numerical simulation of two circular particles sedimentation at different time stages.

Figure 11 .

 11 Figure 11. Evolution of the particles longitudinal velocities vs. time.

Figure 12 .

 12 Figure 12. Snapshots in a chronological order, during the 3D quicksand simulation for a deposit submitted to a gradient i 0:92i c . The picture (1) corresponds to the initial state of the deposit not subjected to the flow.

Figure 13 .

 13 Figure 13. Snapshots in a chronological order, during the 2D simulation of the quicksand of a grain deposit for i D 1:1i c , ıt D 0:065 s.

Figure 14 .

 14 Figure 14. Snapshots in in a chronological order, during the 2D simulation of the quicksand of a grain deposit for i D i c , ıt D 0:065 s.

Figure 15 .

 15 Figure 15. Snapshots in a chronological order, during the 2D simulation of the quicksand of a grain deposit for i 0:968i c , ıt D 0:065 s.

  APPENDIX A: PRESSURE AND VELOCITY BOUNDARY CONDITIONS FOR THE D3Q19 LATTICE BOLTZMANN MODEL We consider that the pressure (or the velocity) is imposed on the boundary of a discretized domain by the D3Q19 model (Figure (A.1)).

Figure A. 1 .

 1 Figure A.1. Boundary of a domain discretized using the D3Q19 Lattice Boltzmann model.

  In this application, we compute the drag force and the torque applied on a stationary spherical grain in a flow between two parallel plates at low Reynolds Number. In addition, analytical approximate solutions established by Wakiya (1957) are used. These solutions give the drag force and the torque applied on a fixed sphere in a creeping Poiseuille flow, that is, a flow between two infinite parallel plates at very low Reynolds number.The configuration of the problem is shown in figure (Figure

< : 1=3 pour i D 0; 1=18 for short velocity vectors like : c.1; 0; 0/; 1=36 for long velocity vectors like : c.1; 1; 0/: (24) 4.2.1. Force and torque on a sphere in a Poiseuille flow.

ACKNOWLEDGEMENTS

The authors would like to express their sincere thanks to the French Research Network MeGe (Multiscale and multi-physics couplings in geo-environmental mechanics GDR CNRS 3176/2340, 2008-2015) for having supported this work.

<

:

According to Eq. (A.1), the momentum densities in the three directions (Eq. (12b)) can be given as:

where e ix , e iy , and e i´a re the three components of the discrete velocity E e i . It can be observed that these expressions (A.1) of unknown distributions lead to inexact momentum densities in the x and y directions, therefore, it is necessary to introduce corrections.

Let err x and err y be the errors on the momentum densities in the x and y directions resulting from equations (Eq. (A.2)) and (Eq. (A.3)):

These errors can be corrected by introducing additional terms in the distribution expressions (Eqs. (A.1)) as follows: add and subtract err x 2c to f 15 et f 7 respectively add and subtract err y 2c to f 17 et f 9 respectively The corrected expressions of the unknown distributions are finally: 8 <

:

It shall be noted that the fluid's inplane velocities in a pressure boundary are commonly assumed zero, that is, for the case of figure (