Discrete modelling of front propagation in backward piping erosion

A preliminary discrete numerical model of a REV at the front region of an erosion pipe in a cohesive granular soil is briefly presented. The results reported herein refer to a simulation carried out by coupling the Discrete Element Method (DEM) with the Lattice Boltzmann Method (LBM) for the representation of the granular and fluid phases, respectively. The numerical specimen, consisiting of bonded grains, is tested under fully-saturated conditions and increasing pressure difference between the inlet (confined) and the outlet (unconfined) flow regions. The key role of compression arches of force chains that transversely cross the sample and carry most part of the hydrodynamic actions is pointed out. These arches partition the REV into an upstream region that remains almost intact and a downstream region that gradually degrades and is subsequently eroded in the form of a cluster. Eventually, the collapse of the compression arches causes the upstream region to be also eroded, abruptly, as a whole. A complete presentation of the numerical model and of the results of the simulation can be found in [12].


Introduction
Piping erosion within embankment dams and dykes or in their foundations is a frequent hazard that may lead these structures to failure [1,2].The process is generally described through four ordered phases of evolution: initiation, continuation, progression and breaching [3].Throughout these different stages, two different erosion mechanisms can be singled out: the enlargement of the pipe is driven by tangential erosion caused by the turbulent flow, while backward erosion induced by Darcy flow prevails at the upstream-propagating pipe tip.Several experimental, analytical and numerical studies on the tangential erosion mechanism have led to substantial advances in the understanding and modelling of the pipe enlargement kinetics but few models for localised backward erosion processes in porous media have been developed [4].This paper provides a brief account of a first, discrete numerical investigation of the backward erosion process at the pipe tip.The Reader is referred to [12] for a more extended description of this study.Modelling of both the cohesive granular soil and the pore fluid phase in the front region of the pipe is obtained herein by coupling the Discrete Element Method (DEM) with the Lattice Bolzmann Method (LBM), respectively.Recently, a similar DEM-LBM coupling scheme has been used to study the tangential erosion process at the walls of the erosion conduit [5,6].

DEM-LBM coupling
An in-house 2D DEM code was used to investigate the micromechanical (grain-scale) processes underwent by the granular skeleton in the pipe tip region.The code was developed following a standard molecular dynamics approach [7].Grains are represented by circular discs and their interactions are modelled by enabling normal (F n ) and tangential (F t ) forces at the contact points.They can result from a bilateral (bond between grains) or a unilateral contact (no bond).The failure envelope for the bonds is characterised by: where the strength parameters A > 0 (bond adhesion), C > 0 (bond cohesion) and μ > 0 (friction coefficient) are chosen such that C ≥ μA.As the bond breaks, the normal interaction becomes unilateral and the tangential force limit is computed on a purely frictional basis.Gravity is neglected.The same model is employed for the interaction between a grain and a rigid confining wall.
As for the fluid phase, LBM is a kinetic-theory-based numerical approach to fluid dynamics problems.The method is hinged upon a discrete form of the Boltzmann equation and therefore on the determination of the distribution function.The latter is defined as the probability density for the presence of a fluid particle with velocity c at position x (as random variables) and time t (as a parameter).In the in-house code developed for this study, the Boltzmann equation is discretised according to the ninevelocity square lattice model D2Q9 [13].The Multiple Re-laxation Times (MRT) approximation [10] is used for the collision operator of the Boltzmann equation, to overcome some known drawbacks of the standard Bhatnagar-Gross-Krook (BGK) approximation [8] and improve numerical stability [14].A detailed description of the in-house LBM code and its validation can be found in [12].
No-slip conditions at stationary as well as moving solid boundaries, straight or curved, are imposed through an interpolated "bounce-back" scheme [9].The actions exchanged between the fluid and the grains are computed by the LBM algorithm (in terms of momentum transfer from the fluid nodes) and summed, in the DEM algorithm, to the resultant forces and force moments acting on each grain due to contact interactions.
Fluid transport is not possible through a 2D granular packing, due to the pore network not being connected.To overcome this limitation and enable, in the 2D model, typical permeabilities and drag forces of 3D assemblies, the hydraulic radius of a grain in the LBM model is set (reduced) to 0.8 times the reference radius of the same grain in the DEM model [11].

Backward erosion test 3.1 Preparation of the granular specimen
The arrangement in Figure 1a represents the numerical specimen employed to model a granular soil REV located on the upstream side of the pipe face (see Figure 2).It consists of 800 circular grains with a core mass density of 2.65 • 10 3 kg m −2 and radii randomly dispersed in the range from 0.75 mm to 0.95 mm.The specimen was obtained by a "dry" preparation procedure.A preliminary "lubricated" (i.e., with null contact friction) isotropic compaction was performed, starting from an initial randomly-dispersed configuration, up to a confining pressure of 30 kN m −1 .The normal and tangential contact stiffnesses were both set to 5.4 • 10 7 N m −1 .
The normal and tangential contact damping coefficients were set at about 80% of the critical value for an harmonic oscillator characterised by the same mass as the average grain.The contact friction coefficient μ (which is actually the same parameter appearing in the bond failure criterion) was chosen equal to 0.5.The contact adhesion was set as A = 1 N, which is consistent with an extremely weak macroscopic tensile traction resistance (of order of magnitude 10 2 N m −1 ) and the contact cohesion was set as a multiple of the contact adhesion (C = 4 N).New contacts created after the initialisation of the bonds were considered by default as unilateral frictional-visco-elastic.The time step of the DEM algorithm was set to 2.5 • 10 −7 s.

Flow-induced erosion of the granular specimen
During the erosion test, the specimen resulting from the preparation procedure was confined by a rigid wall on the left (upstream) wall.The latter was maintained at the same position as at the end of the preparation procedure.The right wall used for confining the specimen during compaction had been carefully removed by the end of the preparation procedure.The right (downstream) side of the specimen is therefore unconfined.The degrees of freedom of the black-coloured grains in Figure 1a were locked in order to form two oblique rough boundaries at the top and the bottom, respectively.The resulting trapezoidal shape of this REV region, where the gray-coloured grains are mobile (erodible grains), is reminiscent of an angular sector at the supposedly curved soil-pipe interface in the front region.Figure 1b represents the hydraulic boundary conditions implemented in the simulation, as well as the "hydraulic configuration" of the specimen at the beginning of the test.An incompressible fluid flow, confined by the impervious walls located at the top and the bottom of the sample, was enforced through the REV by the pressure difference Δp = p out − p in < 0 between the right (outlet) and the left (inlet) boundaries of the fluid domain.The fluid is characterised by a mass density ρ f = 10 3 kg m −2 and a dynamic viscosity of 10 −3 N s m −1 .The evolution of the normalised pressure difference −Δp/(ρ f gL), where g is the standard gravity acceleration, is plotted in Figure 3i for the whole duration of the simulated test.The set of parameters controlling directly or indirectly the numerical implementation of the LBM model was chosen consistently with a number of requirements.In particular, the spacing between neighbouring lattice nodes was set to 5 • 10 −5 m (i.e., 1/12 of the smallest hydraulic radius) and the lattice time step to 10 −5 s, which enables an accurate estimate of the hydrodynamic actions.The complete set of LBM parameter and the relevant criteria can be found in [12].
The conservation of the fluid mass is not an implicit feature of the implemented no-slip condition for moving boundaries.Therefore, the total fluid mass was recorded during the first part of the erosion test, up to time t = 1.55 s, i.e., before the first grain started exiting the fluid domain through the channel outlet.An increase of 0.63% was measured, which is considered consistent with the required accuracy.

Results
The evolution of the erosion process during the test is quantified in Figure 3ii in terms of the eroded mass fraction M er /M 0 , where M er and M 0 are the eroded and erodible masses, respectively.At a given time instant, the eroded mass is defined as the cumulative mass of the grains having crossed the dashed line marking the right (downstream) limit of the configuration of the granular specimen in Figure1.
At the beginning of the test, for low-enough values of the pressure difference, no erosion was observed.Visual information on the process at this stage is provided in Figure 4a, corresponding to the marks with label (a) in Figure 3.Even at very low values, the hydraulic load affects the stress transmission through force chains within the specimen: a concentration of compressive and tensile force chains can be observed roughly at the upstream and downstream sides of the specimen, respectively.It can be inferred, already from Figure 4b and Figure 4c, that the main resistance mechanism developed in response to the action of the hydrodynamic forces consists of a few main arches of compressive force chains.These arches divide the sample into two parts: the compressive upstream region, which is essentially preserved during the process up to Figure 4(e), and the downstream region that is degrading due to bond breakage as the pressure difference increases.Backward erosion was observed for larger values of the driving pressure difference, roughly in the branch from (b) to (d) in Figure 3 (given the considered definition for eroded grains).The eroded mass is detached and transported by the flow in the form of a large, deformable cluster of grains.The increasing branch after point (d) of the eroded mass diagram in Figure 3ii is to be referred to previously mobilised grains progressively exiting the right limit of the measure volume: no significant further detachment of granular material was observed up to point (e).This was due to the resistance provided by the persisting arches of force chains.The erosion kinetics, by taking the form of released clusters of particles, is compatible with what is observed on site.After point (e), the final increasing branch in Figure 3ii corresponds to the final collapse of the persisting arches, resulting in all the remaining material being washed away simultaneously.Finally, the plots in Figure 5(a-e) illustrate the evolution of the norm of the hydrodynamic forces on the grains, during the erosion test, and refer to the same configurations identified in Figure 3. Spatial variability of the hydrodynamic forces is observed at the scale of the specimen, with a tendency for lower values towards the downstream side of the granular region.Some randomness at the scale of a few grain diameters can also be observed, especially in Figures 5d and 5e, i.e., for larger values of the imposed pressure difference.Figures 5b to 5d also illustrate the complexity of the interaction between the fluid and solid phases during the backward erosion process.

Conclusion
A numerical simulation of the backward erosion process for a REV of cohesive granular soil, at the front region of an erosion pipe, has been performed using a coupled DEM-LBM approach.Throughout the test, the most distinct effect that was observed was the marked arching through force chains, as a self-organised response of the contact/bond network to drag forces and couples on the grains.The main resisting arches that transversely crosses the sample partitioned the REV into a compressive upstream region, that remained almost intact during the largest part of the test, and a downstream region that gradually degraded as the pressure difference increased.Eventually, the downstream material was eroded and clusters of particles were released.The collapse of the compressive arches of force chains, by the end of the test, caused all the remaining material to be washed away as a whole.

Figure 1 .
Figure 1.Simulation of backward erosion: geometric configurations (L = 66.8 mm, H = 33 mm), boundary conditions and initial particle arrangements for the coupled models developed with DEM (a) and LBM (b).

Figure 2 .
Figure 2. Sketch of the granular REV at the soil-pipe interface in the front region