Coupled DEM-CFD Analysis of the Initiation of Internal Instability in a Gap-Graded Granular Embankment Filter

Internal instability is a form of internal erosion that can occur in embankment dams or flood embankments where the finer fraction of the material is washed out under the action of seepage flow; if undetected this process can progress to cause embankment collapse. Gap-graded materials are particularly susceptible. Skempton and Brogan [1] proposed that a key contributor to instability is the reduced stress transmitted by the finer fraction and that the magnitude of this reduced stress could be inferred from the hydraulic gradients observed at the initiation of particle migration in experiments. Here Skempton and Brogan’s hypothesis is assessed at the particle scale using a discrete element method (DEM) model coupled with computational fluid dynamics (CFD). This contribution discusses validation of the coupled DEM-CFD software prior to describing the simulation of a permeameter experiment. The simulation generated particlescale data at the initiation of instability by considering a gap-graded sample subject to at a hydraulic gradient of 1.0 (upward flow). The results provide insight into the instability mechanism, most notably showing that while the particles that move under seepage flow do indeed transmit relatively small effective stress, a finite proportion of the particles that move transfer relatively large stresses.


Introduction
Internal instability is a form of internal erosion which can occur in water-retaining earth structures such as embankment dams and flood levees.Internally unstable soils typically have fairly distinct coarse and fine fractions (i.e.gap or broadly-graded soils), and it is this fine fraction which can be preferentially eroded by seepage flow [2].While initially this phenomenon might not be associated with any noticeable deformation of the embankment, the migration of particles can cause changes in local permeability and the overall flow field and can progress to a scenario involving significant deformations, possibly causing embankment breach.In recent years there have been a number of experimental research studies considering this mechanism (e.g.[1,[3][4][5]) which have hypothesized the fundamental particle-scale interactions that trigger this instability.A particular hazard of internal instability is the erosion of the fine particles which can initiate at much lower hydraulic gradients than would be predicted by the continuum analyses used in engineering practice which assume a homogenous soil.Skempton and Brogan [1] attributed this phenomenon to the lower proportion of the effective stress transmitted by the finer particles in comparison with the coarser particles.They proposed that this could be quantified using the stress reduction factor, Į: where ıƍ fine is the effective stress transferred by the fines and ıƍ is effective overburden stress.Skempton and Brogan [1] inferred Į by comparing the hydraulic gradient at which erosion was observed in a permeameter test with no surface loading (i c ) to the the hydraulic gradient required for zero effective stress from Terzaghi's classical theory (i c,heave ): In order for the fines to carry reduced effective stress they must be insufficient to fill the voids between the coarser particles.Depending on sample porosity the fines should just fill the voids between the coarse particle at a fines content of between 24% and 29% by volume [1,6].Shire et al. [6] carried out a micro-scale study of the stress distribution within gap-graded soils using the discrete element method (DEM).By directly measuring the stress within coarse and fine fractions of the soil they demonstrated that Į is a function of fines content, packing density and the size ratio between coarse and fine particles.However, as no water was applied they were unable to directly link Į to the hydraulic gradient at which erosion initiates, i.e. they could not confirm that Equation 2 is valid.
In this paper the initiation of internal instability was modelled using a coupled discrete element methodcomputational fluid dynamics (DEM-CFD) approach, where DEM was used to model the soil particles and CFD the seepage of water.The simulation approach and validation are first introduced, then results for a representative sample are analysed in detail.

Simulation approach 2.1 Methodology
The simulations were carried out using PFC3D with the CCFD add on as detailed in [7].Using the DEM data the stress tensor for each particle, p ij V was calculated using the method described in [8].The average stress tensor for the entire volume, ij V is then: where V = total volume (particles + voids); N p = number of particles; p ij V = average stress in particle p; and p V = volume of particle p.The stress tensor for the finer fraction where n =sample porosity; N p,fine = number of fine particles.As an isotropic stress state is applied here, the stress reduction factor, Į, is given in terms of mean normal stress pƍ = 1/3(ı ii ƍ): The coarse-grid approach to DEM-CFD coupling proposed by [9,10] was adopted here.The CCFD software discretises the fluid domain into a regular fluid grid, which is overlain on the DEM domain.The averaged Navier-Stokes equations are solved to obtain average values of fluid pressure and velocity for each cell within the grid and the empirical Di Felice [11] expression is used to calculate the fluid-particle interaction force, which is added as an additional force term in the DEM simulation.Full details of the coupling scheme are available in [7].

Validation using a single column of particles
The coupled software system was validated using the case proposed by Suzuki et al. [12] in which 100 monodisperse spherical particles were placed at diameter intervals along the centre line of a vertical square column, which was filled with water.Free slip, prescribed velocity and prescribed pressure conditions were imposed at the lateral, base and top boundaries respectively.Suzuki et al. [12] derived an expression for the distribution of pressure in the system using the Ergun expression for drag.
Following their approach a similar expression was derived for the Di Felice drag expression as detailed in [13].Representative data comparing the values obtained using this analytical expression with DEM simulation results are presented in Figure 1.For the case illustrated, the model parameters are identical to those used by Suzuki et al. [12], with the exception of the fluid cell length, which was twice the length that used by [12].After the particles were generated, a gravitational body force and buoyancy were applied and the system was cycled until the particles came to rest.Then, an upward fluid velocity of v f = 0.005 m/s was applied at the bottom boundary and a pressure of P = 0 Pa was applied at the top boundary.The distributions of water pressure in excess of hydrostatic pressure are compared in Figure 1.An exact match between the analytical and experimental data was not expected as the analytical expression is based on application of the heat equation to a continuum material.Consequently, the agreement is considered to be reasonable as the DEM-CFD generated data were within 5% of the analytical values.The excess pore water pressure values obtained from the analytical expression derived with the Ergun drag expression are also included to show the sensitivity of the predicted system response to the drag expression used.
Figure 1.Distribution of excess fluid pressure from the validation case

Permeameter test simulation
For the main simulations a Hertz-Mindlin contact model was used.The shear modulus was 27.0 GPa, Poisson's ratio was 0.3, particle density was 2,670 kg/m 3 and coefficient of interparticle friction was 0.3.The particle timestep was 2.0×10 -8 s and the fluid timestep was 1.0×10 -3 s.A sample of 20,000 spherical particles was generated by creating a non-contacting cloud of particles to a specified particle size distribution (PSD) within a rigid-walled cubical cell.The particle size distribution is shown in Figure 2. The particles were isotropically compressed to a mean normal stress of pƍ = 50 kPa.Gravity was then applied to the sample.The PSD is classified as borderline internally unstable according to the empirical guidelines proposed by Kézdi [14], with a Dƍ Once the sample was in equilibrium under gravity the coupled DEM-CFD simulations were initiated, details of which can be found in [13].Free slip and prescribed pressure conditions were used on the lateral and top/bottom boundaries respectively.The sample was allowed to come into equilibrium under gravity and buoyancy.A hydraulic gradient of i = 1.0 was then applied across the top and bottom boundaries to induce an upward flow of water.

Soil fabric
Referring to Table 1 the soil fabric was quantified by considering porosity, average coordination number and Į.Data from Shire et al. [6] are also included relate to a second DEM sample isotropically compressed using the open-source DEM code Granular LAMMPS [15] with periodic boundaries (no gravity) and having with the same PSD and input variables.Despite the differences in sample preparation and boundary conditions the porosity and stress-reduction factor values agree well.However the coordination number for the current study is much higher than that from [6], in which gravity was neglected so that non-stress transmitting particles remained suspended within the voids.
Figure 3 shows the cumulative distribution of normalised particle mean normal stresses തതതത ᇱ .It can be seen that the fine particles (~25% by volume) carry minimal stress for both datasets.These data along with the D values indicate that there are insufficient fine particles to fill the voids between the coarse particles, leaving the fine particles understressed within these voids.The coarser particles (75% by volume) have a heterogeneous stress distribution, which is also independent of sample preparation method.

Fluid-particle interaction
Upon application of a hydraulic gradient of i = 1.0, the permeability was determined to be 2.3×10 -4 m/s, similar to the value which would be obtained from the empirical Hazen formula [16].Figure 4 shows the change in coordination number with time following the application of the hydraulic gradient.There was a large initial drop as understressed fine particles began to move due to the upward flow of water.The coordination number then increased as these mobile fine particles became trapped at the void boundaries between coarse particles.However, after around 0.45 s the coordination number became stable at a value lower than the initial value, indicating that some particles were not retained within the voids and are therefore potentially erodible.Figure 5 presents the relationship between the net change in connectivity and displacement.The histogram represents how many particles had a given connectivity change.Generally, particles which had a net reduction in connectivity moved further than those which gained contacts, demonstrating that fine particles which lose contacts are less likely to be retained by the coarser particles.The relationship between particle stress and displacement is shown in Figure 6.It can be clearly seen that particles with lower initial stress undergo much more displacement than those with higher stress.However, some particles with higher initial stress undergo significant displacement.These highly stressed but erodible particles could be significant as their erosion could lead to buckling of force chains and volumetric strain of the soil.

Conclusions
This paper presents the results of a coupled DEM-CFD analysis of internal instability.Following software validation, a permeameter test on a borderline internally unstable gap-graded soil was simulated.The conclusions are: (a) The coupled DEM-CFD simulations were successfully validated using a case proposed by [12].Permeability values were also close to those suggested by the empirical Hazen rule

Figure 2 .
Figure 2. Particle size distribution of sample

Figure 3 .
Figure 3. Cumulative distribution of particle stresses

Figure 4 .Figure 5 .
Figure 4. Change in coordination number following application of fluid flow
[16].(b) In coupled DEM-CFD simulations the values of parameters such as fluid pressure are quite sensitive to the choice of empirical drag expression, e.g.[11,17].(c) The stress transfer characteristics of the DEM sample were insensitive to preparation method.(d) On application of fluid flow many fine particles were mobilised, resulting in a large drop in coordination number.(e) More mobile fine particles tended to have a larger net drop in connectivity, suggesting they were less likely to be retained by the coarse particles.(f) Most particles which underwent significant displacement had low initial stress.However, some highly stressed particles were displaced, and this could cause volumetric strain in a real soil.