Next Article in Journal
Reanalysis of Soil Moisture Used for Rainfall Thresholds for Rainfall-Induced Landslides: The Italian Case Study
Previous Article in Journal
Simulations of the Soil Evaporation and Crop Transpiration Beneath a Maize Crop Canopy in a Humid Area
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Modeling Water-Induced Base Particle Migration in Loaded Granular Filters Using Discrete Element Method

1
School of Civil and Hydraulics Engineering, University of Ningxia, Yinchuan 750021, China
2
Institue of Solid Mechanics, School of Physics and Elecronic-Electrical Engineering, Ningxia University, Yinchuan 750021, China
3
Department of Civil Engineering, University of Engineering and Technology Lahore, Punjab 54890, Pakistan
*
Author to whom correspondence should be addressed.
Water 2021, 13(14), 1976; https://doi.org/10.3390/w13141976
Submission received: 21 June 2021 / Revised: 12 July 2021 / Accepted: 15 July 2021 / Published: 19 July 2021
(This article belongs to the Section Hydraulics and Hydrodynamics)

Abstract

:
Results are reported from a series of filtration tests simulated using coupled computational fluid dynamics and the discrete element method (CCFD-DEM) to investigate the factors controlling the mechanism of base particle erosion and their subsequent capture in loaded granular filters. Apart from geometrical factors such as particle and void sizes, the filter effectiveness was found to be controlled by the magnitudes of the hydraulic gradients and the effective stresses. The results of numerical simulations revealed that the base soils exhibit significant stress reduction that reduces further due to seepage, and the base particles migrate into the filter, bearing very low effective stresses (i.e., localized piping in base soil). Based on the limit equilibrium of hydraulic and mechanical constraints, a linear hydromechanical model has been proposed that governs the migration and capture of base particles in the filter (i.e., filter effectiveness avoiding piping) for cases simulated in this study. Nevertheless, the proposed model agrees closely with the simulation results of this study and those adopted from other published works, thereby showing a reasonable possibility of using the proposed model as a measure of retention capacity of loaded protective filters.

1. Introduction

In practice, granular filters are installed to protect base soils such as dam cores and subgrades from erosion. Recent advances in geotechnical practices have generated an interest in applying granular filters in geo-environmental and transportation engineering [1,2,3], where filtration occurs in the presence of static and cyclic mechanical loading, e.g., railways and highway filters. Unlike large upstream heads in hydropower dams, the hydraulic excitation in railway sub-structures mainly stems from the development of pore pressure at the interface between natural and engineered fills. The filters within these structures function under complex stress states, which significantly influence the erosion-capture mechanism of fine particles such as base soil migration, suffusion, and internal erosion [4,5]. Thus far, various researchers have studied these processes and proposed numerous concepts and definitions of governing threshold hydraulic gradients, e.g., the critical hydraulic gradient for piping [6,7], segregation piping [8], complete erosion of base soil through filters [9,10], and suffusion or washout failure in granular soils [11,12,13]. The erosion of base soils may occasionally favor the development of self-filtering layers that helps prevent any further erosion and increases the longevity of filters. However, the strong upward seepage could incur significant structural instabilities in the forms of piping, suffusion, and heave development. Alternatively, the reduction in permeability due to base soil eroding into the overlying filter layers would result in progressive clogging, especially under cyclic loading conditions [14].
In this study, the results of filtration tests simulated using a coupled CFD-DEM model are reported, based on which we propose a semi-empirical model governing the migration of base particles into loaded filters to quantify filter effectiveness in retaining base particles while avoiding piping. Following its development, the proposed model is compared with the existing theories, and the results are presented. The results of numerical simulations from the current study plus those adopted from published DEM studies on the topic are used to validate the proposed model. Notably, this study focuses on the analysis of loaded granular filters subject to one-dimensional upward seepage and the erosion of base particles from the horizontal interface between base and filter layers, e.g., drainage layers in transportation embankments (sub-ballast and subbase) and inverted filters at the downstream of embankment dams. The scales of the laboratory and numerical simulations reported here are not comparable with the practical dimensions of the field problems, and this is always a limitation in the modeling of most geotechnical processes.

2. Numerical Model for Base Soil Erosion

In this study, the base and filter particles are modeled as spheres using the discrete element method (DEM) wherein their motions are governed by the force–displacement law and Newton’s second law of motion in tandem [15]. The Hertz–Mindlin non-linear contact model is adopted to simulate particle–particle and wall–particle interactions, while the fluid flow is governed by Navier–Stokes’ (N-S) equations incorporating the effect of particles [16]. For brevity, the basic coupled CFD-DEM modeling principals are outlined in Appendix A.

Model and Simulation Scheme

The simulated systems comprised of a 0.03 m thick layer of gap-graded base soil NB (Figure 1), which is a widely used core material in embankment dams that exhibits greater potential for seepage-induced piping and two uniform filters NF1 and NF2 of thickness of 0.06 m each. Two separate base-filter systems consisting of spherical particles are modeled in a cuboid (0.12 m × 0.12 m × 0.09 m) within five fixed walls (Figure 2). These thicknesses were selected based on a twofold rationale, namely, (1) keeping the size ratio between the smallest cell dimension and the largest particle size well above 5 to avoid potential boundary effects, and (2) in line with the existing studies to accommodate enough particles. For instance, the base layer thickness is 30 mm, while the largest erodible base particle size is 2.8 mm, which yields a size ratio of 10.9 for this study. Note that the conventional definition of relative density (Rd) does not apply to the DEM assemblies [17]; therefore, the base-filter systems were modeled at a target porosity of 50% to indicate a medium dense and stable assembly for which Rd ≈ 50% could be assumed. Prior geometrical assessments using retention criteria of Terzaghi [18] ( ( D 15 / d 85 ) 4 ), NRCS [19] ( D 15 4 × d 85 R ), and Indraratna et al. [20] ( ( D c 35 / d 85 S A ) < 1 ) have been made. Here, D 15 and D c 35 represent the particle size at 15% finer by mass and controlling the constriction size at 35% finer by surface area for the filter soil, respectively. d 85 , d 85 R , and d 85 S A represent base soil particle sizes at 85% finer by mass, regraded base particle size at 85% finer by mass, and base particle size at 85% finer by surface area, respectively. It revealed that both the filters NF1 and NF2 were ineffective in retaining the base soil-NB (potential for ineffectiveness; NF2 > NF1), but they were individually internally stable [21,22]. The base particles were assigned different colors according to their group sizes to visualize their positions during the simulations. Similarly, the z-positions of the eroding particles were monitored to determine the particle infiltration depths and hence the retention capacity of the filters.
Table 1 presents the properties and parameters of all model components including balls, walls, and fluid. The simulation domain consisted of 8 × 8 × 8 fluids grid along the x, y, and z-directions, respectively. Keeping the impermeable sidewalls as the slip boundaries, the hydraulic pressure is applied at the bottom (upstream/inflow boundary), while the top of the model (permeable downstream/outflow boundary) is set at zero pressure. A total of 24 filtration cases were simulated for the base-filter systems NF1-NB and NF2-NB under σ v t = 0, 10, 20, 30, 40 and 50 kPa up to a calculation time of 0.3 s. The porosity variations during simulations have been monitored via five measurement spheres: one installed at 0.015 m, three installed at 0.045 m, and only one at 0.06 m from the bottom (Figure 2b). Nevertheless, the calculation time influences the process of erosion; the current time was deemed sufficient when compared with 0.04 s [23] and 0.3 s [24] from published studies. The target σ v t is obtained by a servo mechanism that controls the boundary wall velocities ( u ˙ w ) using the following algorithm [21]:
u ˙ ( w ) = G Δ σ
u ˙ ( w ) = G ( σ c u r r e n t σ t a r g e t )
where G is the gain parameter that evolves the following stress increment each cycle:
Δ σ = ( k   n   w × N   c × u ˙ ( w ) × Δ t ) / A
where N   c ,     k   n   w ,     Δ t ,     and     A define the number of wall–particle contacts, their average stiffness, time step duration, and the wall area, respectively. A relaxation factor α sets the stability requirements given by Equations (4) and (5), which calculates G and updates Equation (2) during each cycle:
( k   n   w × N   c × u ˙ ( w ) × Δ t / A ) < α | Δ σ |
G = α A / ( k   n   w × N   c × Δ t ) .

3. Results and Discussion

Figure 3 shows the initial and final effective stress magnitudes in the base and filter at simulation times t = 0, 0.2, and 0.3 s. The base soil exhibited significant stress reduction due to seepage that could eventually neutralize to approximately zero (i.e., piping) for the cases with σ v t ≤ 30 kPa. This observation endorsed the previous experimental observations that a relatively higher magnitude of seepage stress would be required to induce any disturbance to the base soil (e.g., piping) when a greater magnitude of σ v t is applied on the upper filter layer. Notably, the stress distributions in base and filter soils for both samples NF1-NB and NF2-NB almost remained unchanged when the simulation time increased from 0.2 to 0.3 s. This clearly depicts that a simulation time of 0.2 s would be sufficient, beyond which the stress conditions become steady.
Figure 4 shows the contact force fc distributions for specimen NF1-NB when subjected to σ v t = 10–50 kPa. The initial contact force, fc, and hence the contact stress σ c distribution within the filter layer is almost uniform for the given σ v t , but it is reduced significantly inside the base soil layer, although there too it remained uniform. This substantial σ c reduction during transfer from the filter to the base soil clearly indicated that the magnitudes of stress reduction factor ( β = σ b a s e / σ f i l t e r ) in Appendix B for the current base-filter systems were less than unity. Similarly, the magnitudes of σ c inside both the base and the filter layers markedly increased with the corresponding increase in σ v t . The final fc distributions (at t = 0.2 s) in the filter remained almost the same compared to the base, wherein the seepage stresses markedly reduced the base particle contact stresses.
In this study, the hydraulic gradients i a = 20 and 35 were applied on both specimens NF1-NB and NF2-NB. Notably, these high hydraulic gradients were chosen to ensure that the base particles are fully mobilized. For no external loading, the magnitude of minimum hydraulic gradient for the complete erosion of base soil through an ineffective granular filter could be computed using model of Indraratna and Radampola [10] and the current magnitudes have been kept slightly larger. The erosion ratio R e was computed to quantify the extent of base soil erosion and hence the filter effectiveness [23]:
R e = 100 × j = 1 j = n j r   j   3 / k = 1 k = n k r   k   3
where n j ,     n k ,   and   r j ,   r k define the number of eroded and original base particles and their respective radii. The filter was divided into six equal layers of 0.01 m each, and the quantities of eroded base fines and their respective R e values were calculated. Based on the published DEM studies, a filter exhibiting R e < 2.25% could be considered effective [24]. Figure 5a presents the correlation between R e and σ v t , NF1-NB exhibited relatively smaller R e compared to NF2-NB at σ v t = 0 and their R e magnitudes decreased significantly with the increase in the magnitudes of σ v t . For instance, at i a = 35, the specimen NF1-NB exhibited reductions in R e from 6.3% (ineffective) at σ v t = 0 to 2.0% (effective) at σ v t = 40 kPa, which reduced further to almost 0.95% (effective) at σ v t = 50 kPa. Similarly, at i a = 20, NF2-NB showed reductions in R e from 7.5% (ineffective) at σ v t = 0 to 1.4% (effective) at σ v t = 50 kPa. Notably, when sample NF2-NB was subjected to i a = 35, a subtle increase in the magnitude of R e was observed. However, subjecting NF1-NB to a relatively smaller i a = 20 at the same stress levels could substantially reduce the R e values to almost half of that at i a = 35. It clearly depicted that the erosion and hence the effectiveness of a base-filter system would be significantly influenced by the magnitude of applied hydraulic gradient.
Figure 5b,c show the spatial distributions of eroded base fines captured at various depths inside the filters at t = 0.2 s and 0.3 s, respectively. In both cases, the computed R e values were higher near the interface, because the hydraulic gradients were higher at the interface than elsewhere, and R e gradually decreased with the filter depth. Notably for the given i a values, a significant amount of base fines could reach the downstream filter boundary at relatively smaller magnitudes of applied effective stresses, i.e., σ v t ≤ 20 kPa. These fines could have eroded further had either the magnitude of i a increased or that of σ v t reduced. Interestingly, the magnitudes of eroded fines ( R e values) for effective base-filter systems were insensitive to the simulation time because the steady-state conditions could be reached until t = 0.2 s. In contrast, the continuous erosion from ineffective base-filter systems at t > 0.2 s resulted in markedly increased R e values, as shown by red lines in Figure 5b,c.
The expressions of hydraulic gradients for a unit displacement of a base particle inside the filter layer and for complete erosion through the entire filter layer read [10]:
i c 0 = [ d b 2 γ ( 2 + k f c tan ) 3 γ w [ 3 8 D c m 2 + d b δ z ] ] i c r = [ γ d b h f ( 2 + k f c tan ) 3 γ w [ 3 8 D c m 2 + d b δ z ] ] i c p = [ σ v t γ w Δ y i ] × [ d b 2 γ ( 2 + k f c tan ) 3 γ w [ 3 8 D c m 2 + d b δ z ] ] .
where γ and γ w represent unit weights of soil and water, respectively, is the soil’s angle of internal friction in degrees, k f c is the particle contact factor, d b is the base particle size, δ z is the unit displacement of base particle, D c m is the mean constriction size, and h f is the filter depth.
Figure 6a,b plot the critical hydraulic gradient values for unit displacement and complete erosion of base particles versus those predicted by the above Equation (7), respectively. Not surprisingly, the magnitudes of i c 0 computed through Equation (7) agree closely with those recorded during the numerical simulations of this study, and it may be because both models consider the particles to be spheres. Nevertheless, the values of the critical hydraulic gradients would be expectedly higher for actual granular soil particles (i.e., non-spheres). It clearly establishes that the values of i c 0 obtained from Equation (7) would be conservative when applied to the actual field conditions. Therefore, this study adopts both Equation (7) to capture the unit displacement and the complete erosion of base particles through the filter layer (i.e., unit displacement approaching h f ), respectively. Furthermore, the additional hydraulic gradient ( i c p ) arising due to the normal loading needs to be overcome for once by hydrodynamic forces to dislodge an erodible fine and induce unit displacement (i.e., initiation of erosion) and can be given by (Appendix B):
Thus, the expression for the complete erosion of base particles (i.e., both initiation and development) from the filter layer can be given by:
i c t = [ γ d b 2 ( 2 + k f c tan ) 3 γ w [ 3 8 D c m 2 + d b 2 ] ] × ( σ v t γ w h f ) + [ γ d b h f ( 2 + k f c tan ) 3 γ w [ 3 8 D c m 2 + d b δ z ] ] .
Figure 7 presents the base-filter particle distributions for select specimens at the end of simulations i.e., at t = 0.2 s for filter NF2 subjected to 0, 10 kPa, 30 kPa, and 50 kPa. As shown, the magnitude of erosion has been significantly reduced with the increase in the magnitude of normal stress. Similarly, the filter’s ability of capturing the eroded fines has been markedly enhanced with the increase in stress magnitude. This clearly depicts that the filter effectiveness could be increased under higher stress magnitudes. Nevertheless, it is fully consistent with our understanding of increasing the stability and effectiveness of inverted filters by constructing berms and providing additional surcharge at downstream of hydraulic structures.

4. Validation of Proposed Model

Table 2 presents the results of numerical data of this study plus those adopted from the published DEM studies for validation. In Figure 8, the magnitudes of maximum applied hydraulic gradients during DEM simulations, i   a ,   m a x , and those predicted by the proposed model, i   c t (Equation (8)), were plotted. The line of equality demarcates a visible boundary between effective and ineffective base-filter systems, whereby the ones plotting on its right were deemed ineffective and vice versa. The assessment results of filter effectiveness shown in Figure 8 have been quantified in Table 2 in the form of effectiveness of a filter in fully protecting the base soil. For brevity, the current model assesses a filter as ineffective should the applied hydraulic gradient be greater than its anticipated i   c t value (Equation (8)). Then, this assessment is compared with the actual DEM response observed in this study or that reported in the adopted works. As Table 2 shows, out of 40 simulation results from three independent DEM studies including the current work, only five incorrect assessments (i.e., two conservative assessments including N3 and N10 and three unsafe including N11, N25, and N39) are obtained by the currently proposed theoretical model. This acceptably smaller discrepancy in numerical results and theoretical assessments of this study could be because the proposed model is an explicit solution using a continuum approach, while CCFD-DEM simulations of this study capture micro-mechanical responses at the fluid–particle level.

Theoretical Envelope

A graphical illustration of the proposed theoretical hydro-mechanical envelopes (Equation (8)) governing the migration of base particles in loaded filters is presented in Figure 9. It followed a straight-line relationship with the dimensionless mechanical constraint Ɲ m   ( = σ v t / γ w Δ y i ), having a slope equal to the unit i c 0 (Equation (7)). This definition of slope is more realistic because the proposed envelope for a given base-filter system tends to be unique by the virtue of its physical properties. For an internally stable filter under σ v t = 0 (i.e., with no additional surcharge on the filter layer except the self-weight), the i c 0 can be estimated by the ordinate intercept (OA) of the theoretical envelope in Figure 9. Almost all the inter-particle contacts are mechanically active in an internally stable soil [25]; thus, any value of σ v t > 0 can help a filter resist the erosion of particles by virtue of the increased mechanical and frictional constraints. In turn, this requires larger i c r values to trigger the migration of base particles (path BCD), thus imitating Terzaghi’s [6] recommendation of using inverted loaded filters to avoid quicksand conditions at the toe of embankment dams.
Nevertheless, a different response would be expected with internally unstable soils whereby the inter-particle contacts are not mechanically significant, and hence, the load transfer is not efficient. This would not increase the inter-particle friction sufficiently, and hence the theoretical envelope tends not to obey the proposed law. The seepage forces accompany erodible finer fractions from the filter itself, which progressively increases its porosity (thus, the constriction sizes), and then, the base particles require less i c r for erosion through the filter. Hence, a prior assessment of internal stability through an acceptable geometrical criterion is imperative to select internally stable and effective filters.
In Figure 10, the applied i a and the Ɲ m values for the numerical simulation results of this study are plotted with the corresponding theoretical envelopes for the base-filter systems NF1-NB and NF2-NB. Apart from a single incorrect assessment NF1-NB-30 (i.e., test N20 in Table 2), the rest of the base-filter systems are correctly plotted in their respective effective and ineffective zones by the proposed model. Nonetheless, based on the geometrical and mechanical factors, the proposed theoretical envelopes governing base particle erosion for a given base-filter system should be unique.

5. Conclusions

The following conclusions could be drawn from this study.
The results of CCFD-DEM numerical modeling of this study plus those adopted from the published literature showed good agreements with both the experimental observations and the model predictions. Unlike the available particle size distribution (PSD)-based geometrical criteria for assessing the effectiveness of filters, a performance-based hydro-geo-mechanical approach that considers the level of hydro-mechanical excitation (ia and σ v t ) and the geometrical and physical characteristics of filter media (e.g., PSD, Rd, and ∅’ etc.) was presented.
No significant erosion occurs from an effective base-filter system at extended simulation times (t = 0.2 s to t = 0.3 s), whereby steady-state conditions prevail due to the formation of self-filtering layers inside the filter. In contrast, the erosion of base particles continues through an ineffective filter should the simulation time increase. Furthermore, the changes in the magnitudes of applied hydraulic gradients markedly affect the magnitude of erosion.
An enhanced alternative method for quantifying filter effectiveness through the factor of safety against base particle erosion instead of traditional “effective or ineffective” bifurcation may be more efficient and cost effective in practice; for instance, designing filters for properly monitored low-risk structures such as drainage layers subject to low magnitude hydraulic excitation in slopes, embankments, levees, railway, and highway sub-structures.
A novel theoretical model has been proposed to quantify filter effectiveness under static loading condition that could also be presented graphically. Nevertheless, the current simulation results and those adopted from the published literature could sufficiently verify the proposed model, thus enhancing the user confidence in adopting it for various preliminary analyses and assessment purposes.
Needless to mention, the findings from this study only govern internally stable and non-cohesive base-filter systems, but the authors do envisage extending the scope to internally unstable granular filters and cohesive base soils.

Author Contributions

Conceptualization, G.Z. and J.I.; methodology, G.Z.; software, J.I.; validation, J.I., W.M. and H.W.; formal analysis, G.Z.; investigation, J.I.; resources, H.W.; data curation, W.M.; writing—original draft preparation, J.I.; writing—review and editing, G.Z. and H.W.; visualization, W.M.; supervision, J.I.; project administration, G.Z.; funding acquisition, G.Z., W.M. and H.W. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Key R&D Program of Ningxia Hui Autonomous Region, China (No. 2018BFH03010), National Natural Science Foundation of China (No. 51768059) and Helan Mountain Research Scholar Program of Ningxia University (2020–2025).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

Supports received by the second/corresponding author from the Helan Mountain Research Scholar Program, Faculty Development Program and International Postgraduate Tuition Award from Ningxia University (China), University of Engineering and Technology Lahore (Pakistan) and University of Wollongong, respectively, are gratefully appreciated. Private discussions and training of corresponding author on DEM with distinguished professor Buddhima Indraratna, previously at University of Wollongong Australia, are heartily acknowledged.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Basic Coupled CFD-DEM Modeling Principles

The drag forces due to fluid flow induce particle movements by overcoming the inertial and contact (normal and tangential) forces, while the moving particles re-interact with the fluid flow further reducing the magnitude of drag forces; thus, the fluid–particle coupled interaction is modeled. The motion of spheres is governed by Newton’s second law of motion and force–displacement law, which is given by the following momentum equations [26]:
m   p v ˙   p = c f   c + f   d + f   g
I   p ω ˙   p = c r   c × f   c
where m   p ,     v ˙   p ,     ω ˙   p ,     f   c ,     f   d ,     I   p ,     r   c   and     f   g ( = m   p g ) define the mass of a soil particle, its interstitial and angular velocity vectors, contact force, drag force exerted by the fluid on soil particles, moment of inertia, vector connecting the centers of particles in contact, and gravitational force, respectively. The drag force evolves from buoyancy and fluid particle interaction terms and reads:
f   d = ( f   i ¯ 1 n p f ) V   p
where f   i ¯ ,     n ,     , p f   and   V   p are the average fluid–particle interaction vector (i.e., f   i ¯ ( x ,   t ) ), local porosity of soil, gradient operator ( = { / x 1 ,     / x 2 ,     / x 3 } ), average fluid pressure, and average particle volume, respectively. To determine the force f   c between two or more contacting particles, the no-slip Hertz–Mindlin contact model with a linear spring-dashpot system is adopted that consists of tangential and normal contact components. If K   s ,   K   n   and   C   s , C   n define the spring and dashpot constants, d s r ,     d n r     and   v s r ,     v n r define the relative displacement and velocity vectors for solid phase (where, subscripts n and s indicate normal and tangential components, respectively); the normal and shear forces at the contacts read [27]:
f   n = K   n d n r + C   n v n r
f   s = K   s d s r + C   s v s r
The generalized N-S equations described by the two-phase mass and momentum conservation equations could capture the coupled response of fluid (i.e., pressure and velocities) and solid phases. The governing equations could be modified to include the effect of a particulate solid mixed into the fluid phase for non-Darcy flow regimes, i.e., R n ≫ 1 [16]. For instance, the average effects of the coupling force f   b and the porosity n over a finite number of particles can be expressed by the following forms of N-S equations:
ρ f [ ( n v ) t + v . ( n v ) ] + n p f = μ 2 ( n v ) + f   b + n ρ f f   g
( n v ) t + ( n v ) = 0
where ρ f = fluid density, μ = dynamic viscosity of fluid, f   b = body force per unit volume, and nv = Darcy velocity. In this study, a volume-averaged, two-way coupling was implemented to model the fluid–particle interaction, which requires updating the effect of movement of particles on the fluid conditions (velocities and pressures) and vice versa [27]. As shown previously, the governing equations for the fluid phase were formulated with the porosity term n to account for the presence of particles (Equations (A6) and (A7)). The forces on particles applied by the fluid (e.g., f   d ) are assigned locally to each particle, depending on the conditions within the fluid elements that the particles occupy. Based on how porosity is determined, a particle may intersect more than one fluid cell at a time, which must be accounted for in calculations for accurate results. In this study, forces were distributed based on the fractional overlap between the particles and the fluid cells. As a result, a body force f   b , proportional to the particle volume evolves on each particle. For spherical soil particles with radius r i assumed in this study, the drag force acting on each particle reads:
f   d = 4 3 π r i 3 f   b ( 1 n ) .
The total fluid force on an individual sphere reads:
f   f = f   d + 4 3 π r i 3 ρ f g
To consider the effects of externally applied mechanical loads (i.e., forces or stresses) and the motion of solid particles in fluid in tandem, an additional forcing term f   m is added to Equation (A1) to account for the interaction between fluid and solid phases:
m   p v ˙   p = c f   c + f   d + f   g + f   m .
In differential form, Equation (A10) is modified to read:
v   p t = f   f + f   m m   p + g
where f   f   ( = c f   c + f   d ) and g define the combined effect of fluid forces on the particles and acceleration due to gravity, respectively. Notably, the resultant fluid force f   f is applied at the centroids of particles. The coupling is implemented in the form of data exchanges during the calculation cycle that involve solving Equations (A1)–(A5) for the updated particle positions, contact forces, and porosity of fluid elements in PFC3D and then Equations (A6)–(A11) are solved for the fluid velocity and pressure gradient in each fluid element [27]. The fluid force f   f moves the solid particles, resulting in porosity deteriorations for different fluid cells, which subsequently affect the fluid flow rate of the grids. The magnitudes of effective stresses in the base ( σ   b a s e ) and filter ( σ   f i l t e r ) could be quantified using the following relationships [22]:
σ   f i l t e r = 1 V   f 1 c r   f [ f   n + f   s ] σ   b a s e = 1 V   b 1 c r   b [ f   n + f   s ]
where V   f ,   V   b = volume fractions and r   f ,   r   b = particle radii of filter and base, respectively, and f   n ,   f   s = normal and tangential vector components of contact forces.

Appendix B. Equivalent Hydraulic Gradient for the Mechanical Constraint

For a granular medium of arbitrary depth Δ y i (from D f   to   h f ) under an effective stress σ t and subject to upward flow due to i a (Figure A1), the bottom effective stress σ b reads:
σ b = σ t + γ . Δ y i γ w . Δ y i
where;
Δ y i = [ Δ y ] D f h f .
Figure A1. Vertical effective stresses on soil volume subject to the upward seepage flow.
Figure A1. Vertical effective stresses on soil volume subject to the upward seepage flow.
Water 13 01976 g0a1
At the onset of erosion into a filter, the stress σ b on the eroding base particles should be zero:
i c p = ( i c t , 0 + σ v t γ w . Δ y i ) .
Equation (A15) governs the uniform soils where all particles participate equally in transferring σ v t . However, with a non-uniform soil, the share of coarse particles due to their larger surface areas is more than with fine particles. The contribution of fines declines further as the seepage forces increase, and it becomes zero when they begin to erode (quick condition). Thus, a stress reduction factor β (= i c 0 / i c r , 0 )) is introduced with i c 0 denoting the actual hydraulic gradient governing the particle erosion at σ b = 0 , and i c r , 0 ( = γ / γ w ) is the critical hydraulic gradient for quick condition in base soil [6]. The factor β has physical meanings, whereby the particles will not erode individually until β→1 ( i c 0 i c r , 0 ). However, for the ineffective base-filter systems, the base particles may start to erode even when β < 1 and the filter particles are still intact, i.e., at σ b a s e 0 , it is likely that σ f i l t e r > 0 .
i c p = β ( i c r , 0 + σ v t γ w . Δ y i ) = i c 0 i c r , 0 ( i c r , 0 + σ v t γ w . Δ y i )
At the inception of erosion in internally stable base soils, the term i c r , 0 1 (fines in their loosest state). Similarly, the term i c r , 0 in Equation (A16) represents the contribution of buoyant weight that could be relevant for full-scale problems with non-negligible filter thickness and is negligible in laboratory tests. Therefore, ignoring the first term in Equation (A17) when compared with the magnitude of the second term for a 150 mm thick sub-ballast layer or a 500 mm thick filter layers under more than 30 kPa [28] and 100 kPa surcharges, respectively:
i c p = i c 0 σ v t / γ w . Δ y i

References

  1. Fourie, A.B.; Copeland, A.M.; Barrett, A.J. Optimization of the as-placed properties of hydraulic backfill. J. S. Afr. Inst. Min. Met. 1994, 94, 199–210. [Google Scholar]
  2. Rice, J.D.; Duncan, J.M. Deformation and cracking of seepage barriers in dams due to changes in the pore pressure regime. J. Geotech. Geoenviron. Eng. 2010, 136, 16–25. [Google Scholar] [CrossRef]
  3. Aziz, M. Using grain size to predict engineering properties of natural sands in Pakistan. Geomech. Eng. 2020, 22, 165–171. [Google Scholar]
  4. Chang, D.; Zhang, L. Critical Hydraulic Gradients of Internal Erosion under Complex Stress States. J. Geotech. Geoenviron. Eng. 2013, 139, 1454–1467. [Google Scholar] [CrossRef]
  5. Xiao, M.; Shwiyhat, N. Experimental investigation of the effects of suffusion on physical and geomechanic characteristics of sandy soils. Geotech. Test. J. 2012, 35, 890–900. [Google Scholar] [CrossRef] [Green Version]
  6. Terzaghi, K. Failure of dam foundations by piping and means for preventing it. Wasserkr. Z. Gesamte Wasserwirtsch. 1922, 17, 445–449. (In German) [Google Scholar]
  7. Dirkx, W.; Beek, R.V.; Bierkens, M. Piping the Influence of Grain Size Distribution on the Hydraulic Gradient for Initiating Backward Erosion. Water 2020, 12, 2644. [Google Scholar] [CrossRef]
  8. Skempton, A.W.; Brogan, J.M. Experiments on piping in sandy gravels. Geotechnique 1994, 44, 449–460. [Google Scholar] [CrossRef]
  9. Irfan, M.; Uchimura, T. Modified triaxial apparatus for determination of elastic wave velocities during infiltration tests on unsaturated soils. KSCE J. Civ. Eng. 2016, 20, 197–207. [Google Scholar] [CrossRef]
  10. Indraratna, B.; Radampola, S. Analysis of critical hydraulic gradient for particle movement in filtration. J. Geotech. Geoenviron. Eng. 2002, 128, 347–350. [Google Scholar] [CrossRef]
  11. Li, M.; Fannin, R.J. Comparison of two criteria for internal stability of granular soil. Can. Geotech. J. 2008, 45, 1303–1309. [Google Scholar] [CrossRef]
  12. Indraratna, B.; Israr, J.; Rujikiatkamjorn, C. Geometrical method for evaluating the internal instability of granular filters based on constriction size distribution. J. Geotech. Geoenviron. Eng. 2015, 141. [Google Scholar] [CrossRef] [Green Version]
  13. Feng, Q.; Ho, H.; Man, T.; Wen, J.; Jie, Y.; Fu, X. Suffusion Internal Stability Evaluation of Soils. Water 2019, 11, 1439. [Google Scholar] [CrossRef] [Green Version]
  14. Trani, L.D.O.; Indraratna, B. Assessment of Subballast Filtration under Cyclic Loading. J. Geotech. Geoenviron. Eng. 2010, 136, 1519–1528. [Google Scholar] [CrossRef]
  15. Cundall, P.; Strack, O. A discrete numerical model for granular assemblies. Geotechnique 1979, 29, 47–65. [Google Scholar] [CrossRef]
  16. Tsuji, Y.; Kawaguchi, T.; Tanaka, T. Discrete particle simulation of two-dimensional fluidized bed. Powder Technol. 1993, 77, 79–87. [Google Scholar] [CrossRef]
  17. Shire, T.; O’Sullivan, C. Experiments on piping in sandy gravels. Acta Geotech. 2013, 8, 81–90. [Google Scholar] [CrossRef] [Green Version]
  18. Terzaghi, K. Soil mechanics—A new chapter in engineering science. J. Inst. Civ. Eng. 1939, 12, 106–142. [Google Scholar] [CrossRef]
  19. Natural Resources Conservation Services (NRCS). Gradation Design of Sand and Gravel Filters. In National Engineering Handbook; USDA: Washington, DC, USA, 1994; Chapter 26, Part 633. [Google Scholar]
  20. Indraratna, B.; Raut, A.K.; Khabbaz, H. Constriction-based retention criterion for granular filter design. J. Geotech. Geoenviron. Eng. 2007, 133, 266–276. [Google Scholar] [CrossRef] [Green Version]
  21. Kenney, T.C.; Lau, D. Internal stability of granular filters. Can. Geotech. J. 1985, 22, 215–225. [Google Scholar] [CrossRef]
  22. Raut, A.K.; Indraratna, B. Further advancement in filtration criteria through constriction-based techniques. J. Geotech. Geoenviron. Eng. 2008, 134, 883–887. [Google Scholar] [CrossRef]
  23. Zou, Y.; Chen, Q.; Chen, X.; Cui, P. Discrete numerical modelling of particle transport in granular filters. Comput. Geotech. 2013, 32, 340–357. [Google Scholar]
  24. Qing-Fu, H.; Mei-Li, Z.; Jin-Chang, S.; Yu-Long, L.; Bao-Yu, S. Investigation of fluid flow induced particle migration in granular filters using a DEM-CFD method. J. Hydrodyn. 2014, 26, 406–415. [Google Scholar]
  25. Langroudi, M.F.; Soroush, A.; Tabatabaie, S.P.; Shafipour, R. Stress transmission in internally unstable gap-graded soils using discrete element modelling. Powder Technol. 2013, 247, 161–171. [Google Scholar] [CrossRef]
  26. Shamy, U.E.; Zeghal, M. Coupled continuum-discrete model for saturated granular soils. J. Geotech. Geoenviron. Eng. 2004, 131, 413–426. [Google Scholar] [CrossRef]
  27. Itasca. Particle Flow Code, PFC3D, Release 4.0.; Itasca Consulting Group Inc.: Minneapolis, MN, USA, 2004. [Google Scholar]
  28. Christie, D. Bulli field trial: Vertical and lateral pressure measurement. In Proceedings of the Rail CRC Seminar; University of Wollongong: Wollongong, NSW, Australia.
Figure 1. Current base and filter soil gradations simulated in this study.
Figure 1. Current base and filter soil gradations simulated in this study.
Water 13 01976 g001
Figure 2. Illustration of (a) simulation cell, (b) numerical model specimen, and (c) the Hertz–Mindlin contact model.
Figure 2. Illustration of (a) simulation cell, (b) numerical model specimen, and (c) the Hertz–Mindlin contact model.
Water 13 01976 g002
Figure 3. Effective stress distributions in specimens; (a) NF1-NB at t = 0 s, (b) NF1-NB at t = 0.2 s, (c) NF1-NB at t = 0.3 s, (d) NF2-NB at t = 0 s, (e) NF2-NB at t = 0.2 s, and (f) NF2-NB at t = 0.3 s.
Figure 3. Effective stress distributions in specimens; (a) NF1-NB at t = 0 s, (b) NF1-NB at t = 0.2 s, (c) NF1-NB at t = 0.3 s, (d) NF2-NB at t = 0 s, (e) NF2-NB at t = 0.2 s, and (f) NF2-NB at t = 0.3 s.
Water 13 01976 g003
Figure 4. Initial contact force distribution at σ v t = (a) 10 kPa, (b) 20 kPa, (c) 30 kPa, (d) 40 kPa, and (e) 50 kPa for specimen NF1-NB (Note: time t = 0 s and t = 0.2 s represent conditions before and after the coupled simulations, respectively).
Figure 4. Initial contact force distribution at σ v t = (a) 10 kPa, (b) 20 kPa, (c) 30 kPa, (d) 40 kPa, and (e) 50 kPa for specimen NF1-NB (Note: time t = 0 s and t = 0.2 s represent conditions before and after the coupled simulations, respectively).
Water 13 01976 g004
Figure 5. (a) Observed relationship between erosion ratio Re and effective stress magnitude in base soil, (b) spatial distributions of particles eroded and captured at various depths within the filter medium at t = 0.2 s, and (c) spatial distributions of particles eroded and captured at various depths within the filter medium at t = 0.3 s. (Note: Red lines in Figure 5b,c represent ineffective filters).
Figure 5. (a) Observed relationship between erosion ratio Re and effective stress magnitude in base soil, (b) spatial distributions of particles eroded and captured at various depths within the filter medium at t = 0.2 s, and (c) spatial distributions of particles eroded and captured at various depths within the filter medium at t = 0.3 s. (Note: Red lines in Figure 5b,c represent ineffective filters).
Water 13 01976 g005
Figure 6. (a) Plot of hydraulic gradient for unit displacement of base particle i c 0 observed [10] and that estimated by Equation (7) versus mean constriction size D c m values, and (b) comparison between experimentally observed and theoretically estimated i c r data for the complete erosion of base soil through filter.
Figure 6. (a) Plot of hydraulic gradient for unit displacement of base particle i c 0 observed [10] and that estimated by Equation (7) versus mean constriction size D c m values, and (b) comparison between experimentally observed and theoretically estimated i c r data for the complete erosion of base soil through filter.
Water 13 01976 g006
Figure 7. Illustration of base particle erosion and capture inside the filter NF2 under σ v t =; (a) 0 kPa, (b) 10 kPa, (c) 30 kPa, and (d) 50 kPa.
Figure 7. Illustration of base particle erosion and capture inside the filter NF2 under σ v t =; (a) 0 kPa, (b) 10 kPa, (c) 30 kPa, and (d) 50 kPa.
Water 13 01976 g007
Figure 8. The comparison between observed and estimated i c r values by the proposed model for numerical modeling data in Table 2.
Figure 8. The comparison between observed and estimated i c r values by the proposed model for numerical modeling data in Table 2.
Water 13 01976 g008
Figure 9. Proposed theoretical envelope to quantify the effectiveness of loaded granular filters.
Figure 9. Proposed theoretical envelope to quantify the effectiveness of loaded granular filters.
Water 13 01976 g009
Figure 10. Validation through coupled CFD-DEM simulation results (test series numbers in Table 2 for data).
Figure 10. Validation through coupled CFD-DEM simulation results (test series numbers in Table 2 for data).
Water 13 01976 g010
Table 1. Summary of parameters for the numerical model.
Table 1. Summary of parameters for the numerical model.
ParametersBallsWallsFluidOthers
Number:
NF1-NB16,7356--
NF2-NB16,2956--
Density ρ (kg/m3)2700-1000-
Friction coefficient0.50.5--
Porosity ( n )0.50.5--
Stiffness (N/m):
Normal ( k n ) 1 × 10 6 1 × 10 6 --
Shear ( k s ) 1 × 10 6 1 × 10 6 --
Time step (s):
DEM--- 1 × 10 6
CFD--- 1 × 10 3
Viscous coefficient μ (N-s/m2)--0.001-
Table 2. Summary of numerical data used for validation.
Table 2. Summary of numerical data used for validation.
Test
Series
Number
Sample
ID
n   f
( % )
i   c t i a , m a x Filter EffectivenessReferences
ObservedCurrent
Model
N1B-F15099.787.5YesYesQing-fu
et al. [23]
(Upward
Flow)
N2B-F25031187.5YesNo
N3B-F3506.187.5NoNo
N4B-F4503.887.5NoNo
N5Base-F15099.76YesYesZou et al.
[24]
(Horizontal
Flow)
N6Base-F15099.730YesYes
N7Base-F15099.783YesYes
N8Base-F250316YesYes
N9Base-F2503130YesYes
N10Base-F2503183YesNo
N11Base-F3506.16NoYes
N12Base-F3506.130NoNo
N13Base-F3506.183NoNo
N14Base-F4503.86NoNo
N15Base-F4503.830NoNo
N16Base-F4503.883NoNo
N17NF1-NB-0509.135NoNoCurrent Study
(Upward
Flow)
N18NF1-NB-105019.235NoNo
N19NF1-NB-205029.335NoNo
N20NF1-NB-305039.435YesYes
N21NF1-NB-405049.535YesYes
N22NF1-NB-505059.635YesYes
N23NF2-NB-0506.120NoNo
N24NF2-NB-105013.520NoNo
N25NF2-NB-205020.920NoYes
N26NF2-NB-305028.320YesYes
N27NF2-NB-405035.720YesYes
N28NF2-NB-505043.120YesYes
N29NF1-NB-0509.120NoNo
N30NF1-NB-105019.220NoNo
N31NF1-NB-205029.320YesYes
N32NF1-NB-305039.420YesYes
N33NF1-NB-405049.520YesYes
N34NF1-NB-505059.620YesYes
N35NF2-NB-0506.135NoNo
N36NF2-NB-105013.535NoNo
N37NF2-NB-205020.935NoNo
N38NF2-NB-305028.335NoNo
N39NF2-NB-405035.735NoYes
N40NF2-NB-505043.135YesYes
Note: Here n f , i a , m a x , and i c t define filter porosity, maximum applied hydraulic gradient, and hydraulic gradient governing particle migration estimated from Equation (8), respectively.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Zhang, G.; Israr, J.; Ma, W.; Wang, H. Modeling Water-Induced Base Particle Migration in Loaded Granular Filters Using Discrete Element Method. Water 2021, 13, 1976. https://doi.org/10.3390/w13141976

AMA Style

Zhang G, Israr J, Ma W, Wang H. Modeling Water-Induced Base Particle Migration in Loaded Granular Filters Using Discrete Element Method. Water. 2021; 13(14):1976. https://doi.org/10.3390/w13141976

Chicago/Turabian Style

Zhang, Gang, Jahanzaib Israr, Wenguo Ma, and Hongyu Wang. 2021. "Modeling Water-Induced Base Particle Migration in Loaded Granular Filters Using Discrete Element Method" Water 13, no. 14: 1976. https://doi.org/10.3390/w13141976

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop