Numerical insight into the micromechanics of jet erosion of a cohesive granular material

. Here we investigate the physical mechanisms behind the surface erosion of a cohesive granular soil induced by an impinging jet by means of numerical simulations coupling fluid and grains at the microscale. The 2D numerical model combines the Discrete Element and Lattice Boltzmann methods (DEM-LBM) and accounts for the granular cohesion with a contact model featuring a paraboloidal yield surface. Here we review first the hydrodynamical conditions imposed by the fluid jet on a solid granular packing, turning then the attention to the impact of cohesion on the erosion kinetics. Finally, the use of an additional subcritical debonding damage model based on the work of Silvani and co-workers provides a novel insight into the internal solicitation of the cohesive granular sample by the impinging jet.


Introduction
The physical phenomenon of surface erosion afflicts often the earthen hydraulic constructions such as earth-dams and levees [1]. In this context, the Jet Erosion Test (JET) has proved useful for the assessment of the sensitivity of soils to the occurrence of surface erosion. The interpretation of results can be made in terms of the critical threshold of hydrodynamic shear stress τc and the erosion modulus kd that quantifies the erosion kinetics [2]. A simple erosion law is thereby commonly adopted as ‫ܧ‬ = ݇ ௗ (߬ − ߬ ) where E is the erosion rate and τ is the hydraulic shear stress. This implies a description of the erosion evolution with one single variable, the shear stress, which is averaged over time and space aiming to represent the hydrodynamic conditions at the fluid-solid interface. However, this is a rough simplification of the complex conditions at the surface under an impinging jet, where actually the shear stress should be zero right at the impingement point.
This paper aims to provide a micromechanical insight into the mechanisms taking place during the jet erosion of a cohesive granular material. Firstly, we introduce a numerical model considering on the one hand a model for cohesive intergranular bonds and then a suitable extension for transient subcritical debonding processes (i.e. damage). Then, the results from a parametric study varying the bond strength are put forward. Finally, the erosion mechanisms within the granular sample at the onset of erosion are briefly discussed in the light of preliminary results with the extended damage model.

Numerical methods
The analysis of fluid-solid interactions at the micro-scale is performed here numerically, combining the Lattice Boltzmann Method (LBM) for the fluid analysis with the Discrete Element Method (DEM) for the solid particles (see further applications for instance in [3,4]).

Solid mechanics
A Molecular Dynamics method has been used here to describe the granular soil as a two-dimensional assembly of round particles whose trajectories are governed by Newton's equations of motion. The interaction between particles for the case of purely frictional contacts is formulated in terms of an interaction force F with normal and shear components and an interaction moment M applied at the common contact point. The normal force Fn depends on the local interpenetration δn through a viscoelastic relationship featuring a normal stiffness kn and damping coefficient ηn. On the other hand, a viscousregularized form of Coulomb's law is used here to compute the shear force arising at a frictional contact in dependence of the sliding velocity ߜ௦ by means of a static friction coefficient μ and the viscous coefficient of shear regularization ks. Finally, the interaction moment is defined by the shear force with the particle's radius as lever arm and a rolling friction component that depends on the relative velocity of rotation through a rolling friction coefficient μω and a coefficient of regularization kω.

Intergranular cohesion
For the analysis of cohesive granular samples we consider that all particles initially at contact are bound by a solid bridge with an elastic rheology characterized by the normal and shear bond stiffnesses kn,b and ks,b. An elastoplastic model with paraboloidal yield surface gu in the space of contact forces (a three-dimensional space in terms of Fn, Fs and M) provides here the limits of cohesion allowing for tensile normal forces as long as they remain in the interior of the yield surface ( Figure 1). Whenever the contact forces reach or trespass the yield surface, the cohesive bond is broken and the contact becomes purely frictional. For convenience, the single thresholds for normal, shear and moment interactions Cn, Cs and Mb have all been set to depend here only on a single parameter C = Cn = 2Cs = Mb/(2Dmean), thereby fixing the shape of the paraboloid for different degrees of cohesion. This way, the parameter C represents the strength, or degree of cementation, of the solid bond and permits the characterization of the relative bond strength in a polydisperse assembly by defining a dimensionless number ‫ܤ‬ = ‫ߩ∆(/ܥ‬ ݃ ܸ) as the ratio of the bond cohesion C to the particle's own buoyant weight.

Fig. 1.
Section of the yield surface of cohesive bonds in the plane of normal and shear interaction forces, after [5]. The third dimension of interaction moments is not depicted here.
In general, it can be noted that B ≥ 0 and that solid bonds with B < 1 tend to be unstable and short-lived since any slight rearrangement of the assembly under its own weight has the potential to cause bond rupture. Here it was further observed that the granular assemblies remained completely bonded under gravity whenever B ≥ 3 for all particles in the sample.

Damage model. Subcritical debonding
The cohesion model presented so far is suitable to reproduce the cohesive response to instantaneous solicitations such as dynamic impacts and transient peaks of applied forces for instance, but it lacks itself of transience beyond the dual state-disjunctive (from intact to broken bond state). A characteristic time can however be introduced for instance by means of the subcritical debonding concept (see e.g. [6]) which permits the possibility of a progressive degradation of the cohesive strength for solicitations contained within the yield surface (i.e. subcritical solicitations). To this end, an additional surface, the damage surface, is defined in the interior of the failure surface discriminating the bond solicitations that induce damage (space of interaction forces bounded by the yield and damage surfaces, depicted in pink colour in Figure 2) from those that cause no bond degradation (space contained within the damage surface, shown in green colour). All other solicitations, outside the yield surface, cause an immediate rupture of the cohesive bond, as in the original model. Fig. 2. Application of the subcritical debonding concept after [6] to the cohesion model of Delenne and coworkers [5]. The third dimension (moments) has been omitted in the figure.
The transience can then be introduced into the model by any suitable definition of a damage variable d and its time derivative, for instance as follows: where 〈•〉 denotes the MacCauley brackets ‫〉ݔ〈(‬ = ‫ݔ‬ if x ≥ 0; ‫〉ݔ〈‬ = 0 if x < 0), η is a characteristic time, C0 stands for the initial damage threshold under pure tensile forces, g0(Fn,Fs,M,d) is the damage criterion and dc is the ultimate value of damage, which depends both on the material and loading parameters [6]. It can be noted that, since the yield and damage surfaces gu and g0 include now a negative dependency on d, this formulation implies that both failure and damage surfaces are actually displaced towards the origin as the damage variable grows, thus increasing the susceptibility of the bond to further damage or rupture.

Fluid dynamics
The the components sα of the diagonal relaxation matrix S (inverse of the different relaxation times) and the fluid material parameters of density ρf and kinematic viscosity ν. Further details can be found in [3]. The method for momentum exchange by Bouzidi and coworkers [8] provides here the coupling between the fluid and solid phases and permits the computation of the hydrodynamic forces on the discrete particles. Since the size of time steps required for the fluid computation is larger than for the discrete solid phase, a sub-cycling time integration technique has been employed, performing two DEM subcycles for each LBM step [9]. In order to retrieve a non-zero permeability through the two-dimensional assembly of solid grains, a reduced "hydraulic" radius is introduced [10]. Figure 3 illustrates the hydrodynamic conditions derived by the consideration of a fluid jet with prescribed velocity at the injection nozzle (nozzle diameter b) which impinges orthogonally on the surface of a granular packing located at a distance H from the nozzle, as described more thoroughly in [11].

Jet hydrodynamics at the soil surface
The geometrical, material and rheological parameters employed for the simulation have been chosen here either for convenience or based on usual values from the literature (see e.g. [12]) and are summarized in Table 1. The flow may be described here as an inertial laminar one (jet's Reynolds number Rej ~ 50 to 200) in the transition from a laminar to a turbulent regime. The flow is confined between the solid surface and the upper, closed boundary, while the lateral boundaries have been left open, and so the appearance of two convective cells can be observed. The profiles of fluid velocity, pressure and shear stress right over the solid surface ( Figure 4) show the characteristic "M"-shape profile of fluid velocity with the stagnation point right under the jet's axis and its complementary maximum of fluid pressure. However, the maxima of shear stress, as well as the highest pressure gradients, are located actually right on the spots of the most prominent grains (the most exposed ones) of the irregular granular surface.  This seems to imply that the topology of erosion, at least at its onset, can be dictated by the irregularities of the surface in the relevant impingement area, approximately in the range (-H, H) of radial distance from the jet's axis.

Soil erosion. Simple cohesion
When the soil grains are allowed to move under the action of the fluid jet with a velocity above the erosion critical threshold, a scouring crater will generally form. To quantify the erosion kinetics, it is necessary to specify an erosion criterion filtering the eroded from the non-eroded particles. Here, we classify a particle as eroded if it attains at any moment a kinetic energy above 2x10 -5 Joules [11]. This was useful to discriminate the grains at the debonding front separating the cohesive assembly from the detached and re-settled particles. Figure 5 shows the time evolution of the relative eroded mass (i.e. the proportion of eroded particles compared to the total mass of the granular assembly) for different strengths of the cohesive bonds and a jet velocity of 1.5 m/s. The purely frictional sample (B = 0) shows a sharp increase of eroded mass when the fluid jet reaches the soil surface and, after 5 seconds, half of the assembly has already been eroded. This proportion is significantly reduced by the introduction of cohesion, while the complete absence of erosion is achieved when the bond strength is B = 25, indicating the critical cementation degree for the beginning of erosion.

Bond damage caused by an impinging jet
The extended subcritical debonding model can provide an insight into the micromechanics of erosion even before the erosion onset. Figure 6 shows the network of cohesive bonds among the soil particles (Dmean = 1mm) in the impingement area after 5 seconds of simulated jet flow with C = 1N, C0 = 0.01N and η = 0.1s. The figure appears to show a pattern of damage with a clear preferential directionality along a tangential polar coordinate centred on the impingement point. A closer inspection of the damaged bonds (Fig. 7) shows as well a second preferential damage direction deviating 30° from the tangential direction. These results, along with the fact that the most damaged bonds are located below the exposed surface, appear to suggest that the fluid pressure gradients are at least a major agent (if not the driving one) behind the jet erosion. Fig. 7. Angular distribution of damaged cohesive bonds. Angular deflections from polar tangent centred on the jet's impingement point. The French regional administration Provence-Alpes-Côte d'Azur has kindly provided the funding for these investigations in the frame of the ESCAPE research project.