Numerical Simulation of Ground Subsidence Factors Resulting from Unpressurized Pipeline Rupture Below the Water Table

Abstract

The rupture of an unpressurized pipeline below the water table can lead to the leakage of groundwater along with soil particles into the pipeline. This not only causes blockages in the pipeline but, more critically, can result in ground subsidence. Understanding the factors influencing this phenomenon is a subject of great interest. To delve into this matter, this study employs the DEM-CFD methodology to synergistically encompass particle dynamics and interactions within the flow domain. It introduces an innovative framework for simulating water and soil erosion subsequent to the rupture of subaqueous unpressurized pipelines. This pioneering approach introduces a novel modeling and simulation paradigm catering to the analysis of intricate phenomena of this nature. Upon validating the flow field, our investigation specifically focused on three key factors: particle friction coefficient, groundwater level, and particle size distribution. We conducted a thorough examination of the process and mechanism of water and soil loss at the pipeline leakage point and the subsequent development of stratum subsidence. Our results indicate that particles with a friction coefficient of 0.6 had a reduced maximum displacement by 8.9%, compared to particles with a friction coefficient of 0.3. Similarly, a groundwater depth of 2 m resulted in a 29.6% decrease in maximum displacement compared to a 4 m depth, with a corresponding 160.9% increase in maximum force chain strength. Discontinuous particle gradation, in contrast to continuous gradation, yielded a notable 40.3% reduction in maximum displacement and a substantial 495.1% increase in maximum force chain strength. This underscores the noteworthy influence of particle friction coefficient, groundwater table elevation, and soil particle diameter on the stability of the overlying soil strata in the vicinity of a compromised unpressurized conduit.

1. Introduction

Urban water transportation heavily relies on dense underground networks of pipelines. However, these pipelines often suffer from leakage due to various factors, leading to soil erosion and ground subsidence [1,2,3,4]. According to a statistical analysis by Hu et al. [5], 55% of urban subsidence incidents in China are caused by pipeline damage, and the likelihood of leakage significantly increases when the pipeline’s service life exceeds 25 years [6]. This escalating issue poses a serious geotechnical hazard, endangering residents, buildings, and infrastructure. It has become increasingly prevalent across various regions, resulting in significant socio-economic losses and becoming a destabilizing factor.

Currently, most attention is directed towards the rupture of pressurized pipelines. Some studies suggest that ground subsidence originates from the erosive action of high-pressure jets on the soil, creating liquefaction phenomena in the surrounding soil layers due to high-pressure fluid scouring [7,8,9]. However, researchers have also highlighted the complexity of the interaction between pipelines and soil, owing to changes in the geometry of unconfined flow regimes, piping, and hydraulic fracturing [10].

Unpressurized pipelines for urban rainwater and wastewater drainage are equally critical infrastructure for modern cities. When these pipelines experience leakage, they cause damage to the surrounding soil. Davies et al. [11] summarized the damage process caused by seepage leakage in sewer pipelines, where the soil at the leakage point becomes fluidized and eroded, forming an erosion zone. Extensive studies have also been conducted on the infiltration rate of water through the pipeline leakage point [12,13,14,15]. The erosion process may last for an extended period, forming a low-strength zone at the leakage point, with groundwater gradually propagating towards the surface [16]. Serpente [17] identified four variables influencing the rate of soil loss after damage to unpressurized pipelines: (1) size of the leakage point; (2) groundwater level; (3) amplitude and frequency of surcharge; and (4) soil properties.

In coastal and river areas, or after continuous rainfall, cities often encounter high groundwater levels. Samani et al. [18,19] conducted an in-depth investigation into the intricacies of groundwater models. with unpressurized pipelines typically located below the water table. In such cases, pipeline leakage leads to the gradual formation of a cavity as soil particles are washed into the pipeline with groundwater [20,21,22,23]. The formation of cavities not only causes ground subsidence but also weakens the support for the pipeline, resulting in structural damage [24,25,26,27].

Previous studies have predominantly relied on macroscopic physical model experiments. Currently, the DEM-CFD approach has been extensively applied across various fields for fluid–solid interaction problems. Kloss et al. [28] validated open-source DEM and DEM-CFD coupling models and algorithms. Li et al. [29,30] established mechanics-based models incorporating fluid volume and discrete element coupling, along with a fluid–solid coupling model elucidating the MFSV vibration mechanism, investigating multiphase mixed transport mechanisms. Nigmetova et al. [31] explored the Two-Fluid Model (TFM) approach. Li et al. [32] utilized the DEM-CFD method to investigate the impact of landslides blocking river channels and shock waves. Yan et al. [33] analyzed the kinematic behavior of underwater concrete during pouring, while Krzaczek [34] investigated rock fracturing under the concurrent influences of thermal and hydraulic actions. Concurrently, the advancement in granular mechanics has laid the foundation for a finer-scale exploration of soil failure resulting from unpressurized pipeline rupture [35]. Guo [36] compared three prevalent DEM-CFD coupling models in soil erosion simulation, and Ibrahim [37] conducted detailed numerical investigations into the influences of crack size, defect angle, and sand filling height above the pipeline on erosion.

This study builds upon previous research and, for the first time, employs the DEM-CFD approach to analyze the factors influencing stratum subsidence caused by the rupture of unpressurized pipelines located below the water table. The investigation focuses on particle size distribution, particle friction coefficient, and groundwater level, and it discusses and summarizes the displacement changes in soil particles, the morphological evolution of the stratum, and the variations in contact force chains among particles under different factors. This approach provides a more intuitive understanding of the micro-mechanisms of soil erosion and subsidence. The study aims to offer design guidelines for future underground pipeline installations in cities with high groundwater levels. The numerical computation tool utilized in this research is the Particle Flow Code version 5.0 (PFC 5.0) software developed by Itasca Inc. (Minneapolis, MN, USA).

2. The Computational Principle of the Particle Flow Model

2.1. Principle of Force in Computational Models

2.1.1. The Law of Force–Displacement

In PFC, prevalent contacts comprise sphere–sphere interactions and sphere–wall interactions. As shown in Figure 1.

Let $F_i^n$ denote the normal contact force between two entities, and $U^n$ represent the overlap displacement. In this case, both entities follow the force–displacement law [38]:

(1)$F_i^n=K^n{U}^n{n}^i$

In the given equation, $K^n$ stands for the normal contact stiffness, and $n^i$ represents the normal direction of contact.

The initial shear force exerted on the particles is zero and is dependent on their motion state. The shear force is computed in an incremental manner:

(2)$\Delta F_i^n=-k^s\Delta U^s$

In the equation, $\Delta F_i^n$ stands for the incremental shear force, $k^s$ represents the tangential contact stiffness coefficient, and $\Delta U^s$ corresponds to the normal stiffness displacement.

2.1.2. Newton’s Second Law of Motion

Based on the force–displacement law, the resultant force and moment on each particle can be calculated. By applying Newton’s second law of motion, the acceleration and angular acceleration of the particles can be determined.

Assuming at time $t_0$, the particle experiences a resultant force $F_x$ and a moment $M_x$ in the x-direction, with a corresponding moment of inertia $I_x$, the translational acceleration in the x-direction can be calculated.

(3)$u_x(t_0)=\frac{F_x}{m}$

The rotational acceleration is given by Equation (4):

(4)$w_x(t_0)=\frac{M_x}{I_x}$

At time $t_1=t_0+\frac{\Delta t}{2}$ the translational acceleration of the particle in the x-direction is given by Equation (5) as follows [38]:

(5)$u_x(t_1)=u_x(t_0-\frac{\Delta t}{2})+u_x(t_0)\Delta t$

The rotational acceleration, as described by Equation (6) [38]:

(6)$w_x(t_1)=w_x(t_0-\frac{\Delta t}{2})+w_x(t_0)\Delta t$

$t_0$ represents the initial time, and $\Delta t$ denotes the time step. $t_1=t_0+\Delta t$.

2.2. Contact Model

In this research, the utilized silty sand possesses certain cohesive properties. Based on prior investigations, the linear contact bond model (linearcbond) was chosen. The essential parameters requiring calibration comprise effective modulus, stiffness ratio, tensile strength, shear strength, and friction coefficient.

As shown in Figure 2, the linear contact bonding model provides an ultra-small linear elastic mechanical representation for the contact interface. When the contact surfaces are bonded, they cannot sustain frictional forces, whereas when they are not bonded, they bear frictional forces. Under the condition of bonding, the contact behavior remains consistently linearly elastic. However, upon exceeding the strength limit, the bonding is disrupted, resulting in the separation of the contact surfaces. Even after separation, the contact behavior remains linearly elastic, but it can support sliding when the shear force satisfies the Coulomb criterion.

In the absence of bonding, the linear contact model is effectively analogous to a linear model. Within the linear model, the contact forces can be partitioned into two distinct components: linear and damping. The linear component accounts for linear elastic behavior and frictional response, while the damping component contributes to the viscous behavior. The linear force is generated by springs with constant normal stiffness $k_n$ and shear stiffness $k_s$ The fulfillment of the sliding condition is achieved by applying the Coulomb criterion to the shear force through the friction coefficient.

2.3. The Equation for Fluid Forces

Within the CFD module, the interaction between particles and fluid is incorporated through the inclusion of additional forces. As represented by Equations (7) and (8) [39].

(7)$\frac{\partial \overrightarrow{u}}{\partial t}=\frac{\overrightarrow{f_{mech}}+\overrightarrow{f_{fluid}}}{m}+\overrightarrow{g}$

(8)$\frac{\partial \overrightarrow{w}}{{\partial }_t}=\frac{\overrightarrow{M}}{I}$

In the given equations: $\overrightarrow{u}$ denotes the velocity of the particles; m represents the mass of the particles; $\overrightarrow{f_{fluid}}$ signifies the total force exerted on the particles by the fluid; $\overrightarrow{f_{mech}}$ represents the external forces acting on the particles; g denotes the acceleration due to gravity; w represents the angular velocity of the particles; I denotes the moment of inertia; and M signifies the total torque acting on the particles.

$\overrightarrow{f_{fluid}}$ is composed of two components: the drag force and the fluid pressure gradient force. This is represented by Equation (9) as follows [40]:

(9)$\overrightarrow{f_{fluid}}=\overrightarrow{f_{drag}}+\frac{4}{3}\pi r^3(\overrightarrow{\nabla }p-{\rho }_f\overrightarrow{g})$

The drag force is determined based on the fluid elements containing the particle. Throughout the computation, the fluid–particle interaction is always applied at the center of the particle, thus inducing no bending moment on the particle. The drag force is represented by $\overrightarrow{f_{drag}}$ as follows [40]:

(10)$\overrightarrow{f_{drag}}=\overrightarrow{f_0}{\epsilon }^{-\chi }$

In the equations, $\overrightarrow{f_0}$ denotes the drag force exerted on a single particle, $\epsilon$ represents the porosity of the fluid element in which the particle is situated, and ${\epsilon }^{-\chi }$ is an empirical coefficient considering the local porosity.

The drag force $\overrightarrow{f_0}$ acting on a single particle can be expressed as follows [40]:

(11)$\overrightarrow{f_0}=\frac{1}{2}{C}_d{\rho }_f\pi r^2\left|\overrightarrow{u}-\overrightarrow{v}\right|(\overrightarrow{u}-\overrightarrow{v})$

The empirical coefficient $\chi$ is defined as follows [40]:

(12)$\chi =3.7-0.65\exp \left[-\frac{1}{2}{(1.5-{\log }_{10}{R}_{ep})}^2\right]$

In the given equations, C_d denotes the drag force coefficient, ${\rho }_f$ represents the fluid density, r is the particle radius, $\mu$ stands for the particle velocity, and v denotes the fluid velocity.

The drag coefficient C_d is defined as follows [39]:

(13)$C_d={(0.63+\frac{4.8}{\sqrt{R_{ep}}})}^2$

The Reynolds number R_ep is calculated using the following formula [39]:

(14)$R_{ep}=\frac{2{\rho }_{f}r\left|\overrightarrow{\mu }-\overrightarrow{v}\right|}{{\mu }_f}$

${\mu }_f$ denotes the dynamic viscosity of the fluid.

2.4. The Equation for Pressure Gradient

In porous media, the interaction forces between the fluid and solid can be represented by the pressure gradient equation, as shown in Equation (15) [38]. In this study, Darcy’s law is adopted to describe this interaction.

(15)$\frac{d_p}{d_x}=-\frac{{\rho }_{f}g}{k}{\mu }_{x0}=-\frac{{\mu }_f{\rho }_f}{k}{\mu }_{x0}$

In the equation, ${\mu }_{x0}$ represents the average flow velocity, and k denotes the intrinsic permeability coefficient.

In the context of low Reynolds number flow in porous media, Darcy’s law is formulated as follows [38]:

(16)$\overrightarrow{v}=\frac{K}{\eta \epsilon }\overrightarrow{\nabla }p$

In the given equation, $\overrightarrow{v}$ denotes the flow velocity of the fluid, K represents the permeability matrix, $\eta$ stands for fluid viscosity, $\epsilon$ represents the porosity matrix, and represents the fluid pressure. It is commonly acknowledged that the fluid is considered incompressible, thus:

(17)$\overrightarrow{\nabla }\cdot v=0$

When imposing pressure gradients on the fluid grid, it determines the average flow velocity of the fluid within that grid using the aforementioned equation. The velocity is considered piecewise linear within each element. Therefore, the grid must be adequately refined, considering that the number of grid cells and its complexity will influence both the computation speed and the accuracy of the calculation. Consequently, the selection of an appropriate grid is of paramount importance. Additionally, ensuring a sufficient number of particles within a fluid unit is crucial to prevent particle jumping and maintain smoothness.

In order to account for the influence of particle motion on fluid flow, the permeability coefficient is computed based on the porosity within the PFC3D model. The permeability estimation is carried out using the Kozeny–Carman equation [38]:

(18)$K(\epsilon )=\{\begin{array}{cc}\frac{1}{180}\frac{{\epsilon }^3}{{(1-\epsilon )}^2}(2r_e)& \epsilon \le 0.7\\ K(0.7)& \epsilon >0.7\end{array}$

In the given equation, r_e denotes the particle radius, and $\epsilon$ represents the default minimum porosity in the system, assumed to be 0.3. To compute the permeability coefficient, the upper limit of porosity is set at 0.7. Consequently, when the porosity exceeds 0.7, the permeability coefficient is considered to be the value corresponding to a porosity of 0.7. This approach ensures the consistency of the permeability coefficient even when the porosity exceeds the specified upper limit.

2.5. Coupled Fluid–Solid Interaction Process

In the context of fluid–particle interaction, the bidirectional coupling process involves a sequence of data exchanges between the fluid solver and PFC3D. Upon activating the CFD module, the fluid–particle interaction forces $\overrightarrow{f_{fluid}}$ are applied to the PFC particles during the solving process. During the iterative computations within PFC3D, the fluid–particle interaction forces $\overrightarrow{f_{fluid}}$ and the porosity ${\epsilon }^i$ of each fluid cell are continuously computed and rigorously updated at specified time intervals according to the PFC commands.

Regarding the fluid–particle interaction, the procedures are depicted in Figure 3 as follows:

3. Establishment of Numerical Model

3.1. Calibration of Material Mesoscopic Parameters

In this study, the linear contact model requires calibration of five parameters: effective modulus, stiffness ratio, tensile strength, shear strength, and friction coefficient.

The calibration of material parameters was conducted by contrasting laboratory direct shear test outcomes with numerical model calculations for direct shear tests. Three sets of soil samples were taken, and direct shear tests were carried out under normal stresses of 100 kPa, 200 kPa, and 300 kPa, respectively. The soil employed for the tests was silty clay with a water content of 15.67% and a density of 1.57 g/cm³. Subsequently, a direct shear test model was established in the PFC software. Various parameters were continually adjusted based on the test results to compute shear stress-shear displacement curves under different normal stresses. These calculated curves were subsequently compared to the experimental results, as depicted in Figure 4.

From the graph, it can be observed that for normal stresses of 100 kPa and 200 kPa, there is a good agreement between the computed results and the experimental results. At 300 kPa, while the alignment is slightly less compared to the previous two groups, the overall data deviation is acceptable. Moreover, the trends are similar, which implies that the parameter values in the numerical model are reasonable. The final calibrated results for the microscale material parameters in numerical calculations are presented in Table 1.

3.2. Model Dimensions and Particle Morphology

In a cuboid region of dimensions 15 m ∗ 6 m ∗ 2 m (X ∗ Y ∗ Z), particles are generated. The generation of particles based on their true sizes results in a particle count that significantly exceeds the computational capacity of the computer. To address this limitation, the soil particles are proportionally enlarged, ensuring that the simulation captures the essential characteristics of the soil particles within an acceptable error range. Following the proportional enlargement, the simulated particles exhibit an average diameter of 0.12 m, assuming a regular spherical shape, and their interparticle contacts are modeled using a linear approach. To facilitate clear visualization of displacement and deformation in different soil layers during the failure process, the soil model is discretized into 12 layers along the Y-direction, each represented with distinct colors.

3.3. Assumptions in the Model

During the course of this study, several assumptions have been made to enable more efficient numerical simulations. However, these assumptions may introduce certain limitations:

(1) The assumption of rigid spherical particles is made. This assumption serves to simplify the geometric configuration of simulations and reduce computational complexity. However, it significantly diverges from the irregular shapes found in actual soil particles. Irregularly shaped soil particles can enhance shear strength through interlocking, and their motion and rotation induce internal structural adjustments, thus affecting the macroscopic shear strength and deformation characteristics of the soil. This can lead to an elevation of the dynamic shear modulus and dynamic damping ratio.

(2) The hypothesis of elastic collisions is employed. It assumes that collisions between particles are purely elastic, neglecting energy dissipation post-collision. Such an assumption is reasonable in many instances, particularly for low-speed collisions or minor deformations. Nevertheless, real particle collisions might entail energy loss and plastic deformation, especially when high-speed impacts or intricate interparticle interactions are present.

(3) The supposition of an absence of internal structure is invoked. It presumes that particles lack internal structure, rendering mechanical and thermal analyses more straightforward. This postulation reduces computational costs in simulations. Yet, in scenarios where the internal structure of particles significantly affects mechanical and thermal properties, accurate emulation of actual phenomena might be compromised.

3.4. Boundary Conditions

To accurately reflect pipeline leakage, this study sets the inlet as a velocity inlet oriented perpendicular to the upper boundary, and the outlet as a pressure outlet oriented perpendicular to the lower boundary. The upper boundary is treated as a free boundary, while the lower boundary is configured as a composite boundary featuring a rectangular aperture measuring 2 m in length and 0.9 m in width, flanked by solid walls on its left and right sides. Additionally, the front, back, left, and right boundaries of the model are all defined as solid walls to constrain the fluid and solid regions.

The upper boundary is utilized to simulate the surface environment, while the lower boundary is designed to investigate the impact of underground pipeline damage on ground subsidence. Since this work primarily focuses on the effects of pipeline breakage without considering pipeline deformation, the lower boundary is modeled with solid walls to represent a typical pipeline boundary. The rectangular aperture in the lower boundary is employed to simulate a pressure-less environment within the pipeline.

This configuration will have certain impacts on the results, primarily manifested in the following aspects:

(1) Constrained Flow Field Distribution: The rectangular aperture condition at the lower boundary simulates a pressure-less environment within the pipeline, yet it constrains the distribution of the flow field, preventing a complete representation of the fluid’s flow characteristics near the damaged pipeline, especially in proximity to the boundaries.

(2) Pipeline Deformation Not Considered: This omission leads to the simulation results of soil erosion and sediment transport not including the effects caused by pipeline deformation, which is a significant factor in real-world scenarios.

(3) Boundary Reflection Effects: Setting the front, back, left, and right boundaries as solid walls serves to confine the fluid and solid regions, but it also introduces boundary reflection effects, impacting the evolution and distribution of the flow field.

Despite these influences, the research outcomes still contribute significantly to the understanding of soil deformation settlement and the development of soil particle migration processes.

3.5. The Establishment of the Flow Field

The flow field is established by importing it through the Python function “itasca.cfdarray.create_mesh”. The fluid mesh consists of a regular arrangement of 1 m ∗ 1 m ∗ 1 m (X ∗ Y ∗ Z) hexahedral cells, which are aligned with the walls of the constructed soil model. A total of 120 fluid hexahedral cells are generated to simulate the flow field corresponding to the groundwater levels at 2 m and 4 m underground, as well as the pipe failure point at a depth of 6 m, resulting from groundwater seepage.

Using the Python programming language, a computational fluid field is generated by implementing the DarcyFlowSolution class. This solution involves importing the numpy, pylab, fipy, and itasca libraries. The class defines a method for establishing the computational fluid grid and solving the liquid pressure equation. By passing parameters such as pressure values, locations where pressure is applied, and liquid viscosity coefficients, the pressure distribution of the fluid is computed. Additionally, the interaction between the fluid flow and particles is calculated.

An opening is positioned at the bottom boundary, where a pressure of 0 Pa is maintained. The fluid pressure field is formulated without accounting for the hydrostatic component, necessitating the explicit inclusion of buoyancy effects. Following this configuration, once the flow field achieves stability, the maximum fluid velocity reaches 1.7 m/s. This velocity significantly exceeds the laminar flow velocities typically predicted by Darcy’s law. Consequently, for the computational analysis, the actual fluid velocity should be utilized as the reference standard. To address this disparity, in this research, the hydrostatic pressure at the top boundary for both groundwater levels will be reduced, ultimately regulating the fluid flow velocity within the range of 0.001 to 0.0037 m/s.

The CFD grids employed in this study are illustrated in Figure 5. Sensitivity tests were conducted within the flow field region using two different resolutions: 120 grids and 960 grids, respectively. Changes in porosity above the opening at the bottom of the model were recorded. The results, as depicted in Figure 6, show that the numerical values from both grid resolutions are close, and their trends are similar, indicating the reliability of the data. However, the high-resolution grid exhibits oscillations and numerical instability in the initial stages of computation. Moreover, the high-resolution grid requires significantly more computation time compared to the low-resolution grid. Consequently, considering the optimization of computational efficiency, the division into 120 grids is deemed reasonable.

The mesoscopic parameters of the flow field are presented in Table 2. Figure 7 illustrates the flow field model.

3.6. The Arrangement of the Measurement System

Considering the dimensions of the current study, measurement circles with a diameter of 2 m were selected and arranged in a tightly packed regular pattern within the specimen. A total of 27 measurement circles were generated. The “measure” command was employed to monitor the variations in several parameters of the soil sample during the permeation process. Figure 8 illustrates the arrangement of the measurement circles.

3.7. Flow Field Model Verification for Authenticity

In order to validate the feasibility of the interaction between the flow field and particles as described above, a simulated experiment of particle free settling within the fluid is conducted to verify the accuracy of the bidirectional fluid–particle coupling method.

In the simulated experiment of underwater particle free settling, three particles with different radii but consistent density and centroid positions are placed within the flow field to allow for their unconstrained settling. The vertical velocities of these particles are monitored during the experiment. The data utilized in the experiment are provided in Table 3.

Given the premise of a particle undergoing unimpeded sedimentation within a fluid medium, the initial velocity of the particle relative to the fluid is zero. As the particle begins to move, according to Stokes’ law, the viscous drag force acting on the particle by the fluid is given by [38]:

(19)$F=6\pi \eta \mu r$

In Equation (19), $\eta$ represents the viscosity coefficient of the fluid, $\mu$ denotes the settling velocity of particles, and $r$ stands for the particle radius.

The gravitational force acting on the particle is given by:

(20)$G=\frac{4}{3}\pi r^3{\rho }_{p}g$

The buoyant force acting on the particle is given by:

(21)$F_f=\frac{4}{3}\pi r^3{\rho }_{f}g$

Initially, the particle undergoes accelerated descent, and the viscous drag force increases. When the particle reaches a certain value of acceleration, the sum of the viscous drag force and the buoyant force equals the gravitational force. At this point, the particle experiences a net external force of zero, and its velocity becomes constant. The equilibrium equation at this stage is as follows:

(22)$G=F_f+F$

(23)$\frac{4}{3}\pi r^3{\rho }_{p}g=\frac{4}{3}\pi r^3{\rho }_{f}g+6\pi \eta \mu r$

In a state of equilibrium, the particle’s descending velocity can be calculated as follows:

(24)$\mu =\frac{2r^2({\rho }_p-{\rho }_f)g}{9\eta }$

By substituting the data from Table 3 into the Equation (24), the free settling velocities for the three types of particles are determined as follows: −1.45 mm/s, −3.27 mm/s, and −5.81 mm/s.

The settling velocities of the three types of particles obtained through the PFC-CFD method are shown in Figure 9. Due to variations in particle sizes, they attain different free settling velocities. Based on the measured data, the final settling velocities for the three particles are −1.51 mm/s, −3.35 mm/s, and −5.86 mm/s, respectively. The simulated results closely agree with the calculated values based on Stokes’ law.

Consequently, it is inferred that the PFC-CFD computational approach employed in this research exhibits a notable degree of precision.

3.8. The Design of the Numerical Simulation Scheme

The design scheme for this study consists of 5 groups, and the specific details are presented in Table 4:

With Scheme 1 and Scheme 2 dedicated to investigating the influence of friction coefficient, while Scheme 3 and Scheme 4 focus on examining the impact of groundwater level, and finally, Scheme 4 and Scheme 5 delve into the analysis of particle gradation.

4. Simulation Results and Discussions

The computer processor model used in this simulation is Intel(R) Xeon(R) CPU E5-2630 v4 2.20 GHz (dual processor), with 10 cores and 20 threads. The machine is equipped with RAM (112 GB available), and the duration for each computational model ranges from 2 to 3 h.

4.1. The influence of Particle Friction Coefficient

The particle friction coefficient refers to the resistance between two particles when they come into contact due to the interaction between their surfaces. It determines the resistance magnitude when particles slide or attempt to slide relative to each other. As shown in Figure 10. At the initial computation stage, the soil with friction coefficients of 0.3 and 0.6 showed minimal changes, and each soil layer had an orderly arrangement with clear boundaries. Importantly, no significant particle loss occurred in the bottom layer of the soil. However, as the computation progressed, Figure 11 illustrates both soil types experienced particle loss from the damaged pipeline section. Specifically, the soil with a friction coefficient of 0.3 lost particles from six layers above the damage site, while the soil with a friction coefficient of 0.6 only lost particles from four layers above. This particle loss caused a funnel-shaped deformation in the lower soil layers, with smaller friction coefficients associated with narrower funnel openings. Additionally, it is worth noting that neither soil type experienced surface settlement, regardless of their different friction coefficients.

The particle displacement contour map reveals that different particle friction coefficients result in varying ranges of particle displacement within the soil, even under the same computation duration. In Figure 12, during the initial stage, both friction coefficients of 0.3 and 0.6 exhibit similar influences on soil displacement, with the numerical values falling within a comparable range. However, as the simulation progresses to 400,000 steps (Figure 13), the soil with a friction coefficient of 0.3 experiences a greater range of particle displacement compared to the soil with a friction coefficient of 0.6. Despite this difference, the numerical values of displacement for both cases still remain approximately within the same range, and there is no observed displacement of soil particles on either side of the model.

Figure 14 and Figure 15 display force chain patterns in two types of soil with varying friction coefficients, both at the beginning and end of the computation. While their forms show no significant differences, the disparity lies in the magnitude of force chain values, with higher friction coefficients leading to stronger force chains among particles. This effect becomes more pronounced as the computation steps increase. Frictional forces emerge when external forces induce relative motion among particles, causing the transmission of internal force chains. This process transfers frictional forces from the surface to the interior of the material, impeding relative sliding between surfaces. The magnitude of frictional forces is influenced by the product of the normal stress on particle surfaces and the friction coefficient. Larger friction coefficients enable more efficient force chain transmission, resulting in the generation of larger frictional forces during relative sliding. Consequently, this leads to slower relative movement between objects.

Based on the analysis of the computational results presented earlier, the specific effects of particle friction coefficient on the process of soil particle slip and discharge can be summarized in the following three points:

(1) Sliding Resistance: Higher friction coefficients result in increased frictional resistance between soil particles, hindering relative motion.

(2) Stability of Water and Soil Particles: A higher particle friction coefficient counteracts water’s wetting effect on soil particles, improving soil stability.

(3) Critical Slope: The particle friction coefficient influences the critical slope at which soil begins to slide. Larger friction coefficients correspond to higher critical slopes, providing greater stability even at steeper angles.

4.2. The Influence of Water Table Level

Considering that particle friction coefficients above 0.3 yielded soil displacement that was not significantly different, a uniform particle friction coefficient of 0.1 was adopted for subsequent investigations.

This choice allowed for a comparative analysis between Figure 16 and Figure 17, highlighting the significant influence of groundwater levels on soil loss. By the 400,000th simulation step, the high groundwater scenario exhibited pronounced soil subsidence, with substantial displacements observed in the middle of the top layer, forming a distinctive funnel-shaped configuration. The surrounding soil on both sides of the leakage opening symmetrically slid towards it, while the soil from the ninth layer above the leakage opening was also affected. In contrast, the low groundwater scenario exclusively affected the seventh layer.

In regions with high water table levels, a greater proportion of soil particles are located within the saturated zone. Buoyancy effects weaken the interlocking behavior among particles, and the higher hydraulic gradient induced by the elevated water level increases their susceptibility to displacement. Consequently, during the initial simulation stages, particles in the high groundwater scenario experience significantly larger displacements compared to those in the low groundwater scenario. Upon comparing the computed values of soil particles in Figure 18 and Figure 19 at the end of the simulation, it becomes evident that a considerable amount of soil particles has been lost. The contact resistance between particles further diminishes, and under the influence of high hydraulic gradients, a stable leakage channel forms within the soil mass. As a result, the displacement of particles in both scenarios is significantly higher at the end of the simulation compared to the initial stages.

Upon comparing Figure 20 and Figure 21, it is evident that at the end of the simulation, the contact force chains between soil particles in both groundwater level scenarios have become sparser and weaker than in the initial stages. In the high groundwater level scenario, a substantial number of particles have been lost due to soil subsidence, resulting in a considerable loosening of the force chains in the top region. At this point, most soil particles are submerged in the groundwater, rendering them saturated. Consequently, under the influence of buoyancy and hydrodynamic forces, the interparticle friction gradually weakens, leading to a significantly weaker force chain strength compared to the low groundwater level condition.

The computed results clearly indicate that groundwater levels significantly affect soil particle displacement and stratum subsidence. When a pipeline rupture occurs, it disrupts the previously stable groundwater flow field, creating a noticeable hydraulic gradient between the unpressurized pipeline’s interior and the surrounding soil. This hydraulic gradient becomes more pronounced with higher groundwater levels, leading to a higher flow velocity of groundwater into the pipeline. Consequently, the water flow exerts a stronger erosive action on the soil particles. Furthermore, in areas with high groundwater levels, a larger portion of the subsurface becomes saturated, weakening the interlocking behavior between particles. When combined with the increased hydraulic gradient, the soil becomes more susceptible to significant damage and erosion.

4.3. The Influence of Particle Size Distribution

The section includes numerical experiments conducted on continuous gradation soil, which are identical to the previous numerical experiments involving high groundwater levels. Therefore, the focus is solely on presenting the computational results for discontinuous gradation soil and conducting a comparative analysis with the findings obtained from the continuous gradation soil calculations.

At 5000 simulation steps, the computational model with 12 layers of particles exhibits not obvious displacement, as shown in Figure 16a and Figure 22. The boundaries between adjacent layers remain distinct, and only a few particles near the damaged section of the lower pipe have escaped, indicating the system’s stability. However, as the simulation progresses to 400,000 steps, the continuous gradation soil experiences substantial downward displacement in the central region, leading to evident surface subsidence. Particle loss has occurred in all layers above the leakage point, with the 9th layer being displaced to the leakage point. Despite varying displacement between layers, clear boundaries persist, and there is no apparent mixing of particles among different layers. In contrast, the gradation soil shows a similar pattern of significant displacement in the central region. However, due to the presence of larger-sized particles, the soil arching effect becomes more pronounced near the leakage point. As a result, only the bottom four layers of particles have leaked out, disrupting the boundaries between layers. Smaller particles migrate more extensively downward along the interstices between the skeleton particles, particularly noticeable near the leakage point. As shown in Figure 23.

Figure 18a and Figure 24 show that both soil under different gradations initially remain stable during the simulations. However, near the leakage point, a small area experiences particle displacement. Continuous gradation particles exhibit a larger variation in displacement area compared to discontinuous gradation particles. Specifically, continuous gradation particles show more even displacement throughout the area, while discontinuous gradation particles primarily experience displacement around the larger particles near the leakage point, resulting in a smaller overall displacement area. As the simulations progress to 400,000 steps, continuous gradation particles display a larger overall displacement area, and the displacement values on both sides of the leakage point axis become more symmetrical. This observation highlights a significant correlation between particle displacement and the distance from the leakage point, with displacement diminishing as the distance increases, resulting in mostly stationary particles in the edge region. In contrast, Figure 25 demonstrates the displacement of discontinuous gradation particles lacks the axis symmetry observed in continuous gradation particles. The occurrence and magnitude of particle displacement become more unpredictable. As particles in the vicinity of the leakage point and central region experience overall displacement, fine particles at the edge region also undergo displacement, providing numerical evidence that the finer particles are transferring through the channels formed between the skeleton particles.

When particle contacts are loose, particles are prone to sliding under external forces. Conversely, densely packed particles form stable force chains through close interlocking and mutual compression. At the beginning of the simulations, continuous gradation particles have more force chains due to their smaller size, resulting in a denser arrangement. However, As shown in Figure 26. These force chains are relatively short and weaker compared to those formed by discontinuous gradation particles. As the simulations progress to 400,000 steps, the density of force chains in the continuous gradation soil significantly decreases compared to the initial stage. This loosening becomes more pronounced, indicating that particle loss initiates from the central region at the lower part, causing the upper soil layers to lose strength. For discontinuous gradation particles, as shown in Figure 27. Force chain density decreases, and values reduce above the leakage point due to particle loss. Nevertheless, the overall skeleton force chain network remains intact, and though force chain strength is lower than at the beginning of the simulations, it is still noticeably higher compared to continuous gradation particles. Force chains align with the direction of applied pressure, and stronger force chains lead to less particle sliding, as they collectively deform to resist external forces.

The size and gradation of soil particles significantly influence their movement towards a leakage pipe under high groundwater levels. Several factors contribute to this phenomenon:

(1) Erosion differences: Finer soil particles, due to their smaller radius, are more susceptible to erosion in saturated conditions, leading to more pronounced movement in water. Consequently, finer soil particles are washed into the pipe more frequently.

(2) Percolation pathways: Soils with different particle gradations form distinct percolation pathways. Coarser particles create larger pores, facilitating easier water flow. The skeleton effect between larger particles provides greater resistance to deformation, although smaller particles may be lost through these larger pores. On the other hand, soils with finer particles create numerous small percolation channels. As smaller particles erode more easily, these channels gradually widen, accelerating the erosion process.

(3) Chain reaction in continuous gradation: Soils with continuous gradation have particle sizes distributed continuously. Smaller particles fill the voids between larger particles, leading to interlocking and embedding effects. However, when saturated and subjected to hydraulic forces, erosion of smaller particles weakens the interlocking effect between medium-sized and larger particles, making it more susceptible to erosion and failure.

4.4. Discussion of Different Soil Types

The study at hand encompasses the simulation of silty clay conditions, a soil type characterized by diminutive clay particles manifesting elevated cohesion and internal friction angles. This attributes to heightened resistance to particle displacement within the context of hydraulic forces. The cohesive nature of clay particles precipitates the establishment of a robust interparticle adhesion framework when subjected to aqueous flow, thereby attenuating the pace of particle detachment. The intricacies of the clay’s pore structure and its permeability orchestrate the dynamics of water infiltration and interstitial fluid flow, ultimately wielding a discernible imprint upon the mode of soil erosion.

In contrast, sandy soil, distinguished by its coarser particle constitution, is rendered susceptible to pronounced displacement and dislodgement in response to the hydraulic scouring phenomenon. A reflection of its relatively diminished cohesive and internal friction attributes, sand particles demonstrate a proclivity to forfeit resistance in the presence of water flux, thereby rendering their susceptibility to entrainment and consequent amplification of the rate of erosional attrition. Beyond this, the augmented permeability characterizing sandy soil engenders active groundwater percolation and soil particle mobilization, thereby heralding ground subsidence and soil mass migration mechanisms that diverge distinctly from the tenets observed in clayey contexts.

5. Conclusions

Numerical calculations explored soil loss and development processes following the rupture of unpressurized pipelines below the water table. Factors considered were particle friction coefficient, groundwater level, and particle size distribution, with a focus on their impact on stratum subsidence. The conclusions are summarized below:

(1) A higher particle friction coefficient results in slower soil loss around unpressurized pipelines below the groundwater level, ensuring greater stability. In practical engineering, the particle friction coefficient equals the tangent of the soil’s internal friction angle. Adjusting the internal friction angle during underground pipeline installation, through methods such as soil replacement or compaction changes, is recommended in areas with high groundwater levels.

(2) The influence of groundwater level on particle loss and stratum subsidence is significant. Higher groundwater levels create increased pressure differentials, making soil particle loss into unpressurized pipelines more likely. To address this, selecting suitable dewatering methods during pipeline installation becomes crucial. Additionally, implementing effective waterproofing measures for the pipelines during later use is essential to prevent potential issues.

(3) The particle size distribution also significantly affects particle transport. According to Stokes’ law, particle movement speed in a fluid depends on their radius. Hence, fine and continuous particle size distributions promote soil particle loss. During pipeline installation and backfilling, considering materials with larger particles for backfilling can help reduce soil loss.

In this study, only a small portion of factors contributing to soil erosion resulting from unpressurized pipeline rupture below the groundwater level were examined, primarily focusing on particle displacement and stratum subsidence. In future research, the analysis could be extended to include other factors such as pipeline materials and broader environmental conditions. Moreover, a more detailed exploration of particle interactions, including aspects such as particle aggregation and settling dynamics, could provide a more comprehensive and systematic understanding of this complex issue. Finally, the lack of actual engineering data for comparison with our research results is a limitation of this study, which restricts the direct applicability of the research findings in real-world engineering projects.

Author Contributions

J.-Y.C. organized the textual work of the entire paper, participated in the establishment of the computational model, and summarized and analyzed the computational results. F.-Y.L. initially proposed this idea, modified the computational model, and provided theoretical support during the analysis process. Y.-W.G. contributed to the establishment of the computational model and collected the computational results. C.P. provided technical guidance during the model construction. M.M. offered more reasonable suggestions for the entire paper. All authors made contributions to the writing of this paper. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Schematic diagram of contact.

Figure 2. Contact model.

Figure 3. Diagram of fluid–particle interaction process.

Figure 4. Comparison between computed results and experimental results.

Figure 5. The CFD grids.

Figure 6. Variations in porosity at different resolutions.

Figure 7. The flow field model.

Figure 8. The arrangement of the measurement circles.

Figure 9. The numerical variation curves of particle settling velocities.

Figure 10. The subsurface displacements at 5000 simulation steps.

Figure 11. The subsurface displacements at 400,000 simulation steps.

Figure 12. The subsurface particle displacements at 5000 simulation steps.

Figure 13. The subsurface particle displacements at 400,000 simulation steps.

Figure 14. The internal force chains within the subsurface at 5000 simulation steps.

Figure 15. The internal force chains within the subsurface at 400,000 simulation steps.

Figure 16. The subsurface displacements at 5000 simulation steps.

Figure 17. The subsurface displacements at 400,000 simulation steps.

Figure 18. The subsurface particle displacements at 5000 simulation steps.

Figure 19. The subsurface particle displacements at 400,000 simulation steps.

Figure 20. The internal force chains within the subsurface at 5000 simulation steps.

Figure 21. The internal force chains within the subsurface at 400,000 simulation steps.

Figure 22. The subsurface displacements at 5000 simulation steps.

Figure 23. The subsurface displacements at 400,000 simulation steps.

Figure 24. The subsurface particle displacements at 5000 simulation steps.

Figure 25. The subsurface particle displacements at 400,000 simulation steps.

Figure 26. The internal force chains within the subsurface at 5000 simulation steps.

Figure 27. The internal force chains within the subsurface at 400,000 simulation steps.

Table 1

The mesoscopic calibration results of the model.

Soil Types	Effective Modulus	Stiffness Ratio	Coefficient of Friction	Tensile Strength	Compressive Strength	Particle Number	Particle Volume Fraction
silty clay	3.5 × 106	1.2	0.3	3.0 × 103	2.2 × 103	19,860	96.5%

Table 2

The mesoscopic parameters of the flow field.

Parameter Name	Symbol	Unit	Numerical Value
Flow field density	${\rho }_f$	${\mathrm{kg}/\mathrm{m}}^3$	1000
dynamic viscosity coefficient	${\mu }_f$	1	1.0
topmost hydraulic pressure	Pu	kPa	20, 40
hydraulic pressure at the bottom aperture	Pb	kPa	0

Table 3

The mesoscopic parameters utilized in simulation experiments of particle settling in water.

Parameter Name	Symbol	Unit	Numerical Value
particle radius	r	m	0.001\0.0015\0.002
particle number	id	/	1\2\3
Particle density	${\rho }_p$	g/cm3	2000
flow field density	${\rho }_f$	g/cm3	1000
viscosity coefficient	μ	1	1.5
gravitational acceleration	g	m/s2	9.81

Table 4

Scheme Design.

Scheme Number	Friction Coefficient	Water Table Height	Particle Size Distribution
1	0.3	2 m	Continuous gradation(0.08–0.16 m)
2	0.6	2 m	Continuous gradation(0.08–0.16 m)
3	0.1	2 m	Continuous gradation(0.08–0.16 m)
4	0.1	4 m	Continuous gradation(0.08–0.16 m)
5	0.1	4 m	Discontinuous gradation(0.04 m and 0.30 m)

References

1. Balkaya, M.; Moore, I.D.; Sağlamer, A. Study of nonuniform bedding support because of erosion under cast iron water distribution pipes. J. Geotech. Geoenviron. Eng.; 2012; 138, pp. 1247-1256. [DOI: https://dx.doi.org/10.1061/(ASCE)GT.1943-5606.0000689]

2. Toshifumi, M.; Naoko, K.; Jun, O. Image analysis of soil failure on defective underground pipe due to cyclic water supply and drainage using X-ray CT. Front. Struct. Civ. Eng.; 2012; 6, pp. 85-100. [DOI: https://dx.doi.org/10.1007/s11709-012-0159-5]

3. Van Zyl, J.E.; Alsaydalani, M.O.; Clayton, C.R.; Bird, T.; Dennis, A. Soil fluidisation outside leaks in water distribution pipes—Preliminary observations. J. Water Manag.; 2013; 166, pp. 546-555. [DOI: https://dx.doi.org/10.1680/wama.11.00119]

4. Wen, S.J.; Cheng, W.C.; Li, D.F.; Hu, W.L. Immobilizing lead using loess and nanoscale zerovalent iron (nZVI)-amended loess: Insights from macroscopic and microscopic tests. Environ. Technol. Innov.; 2023; 31, 103228. [DOI: https://dx.doi.org/10.1016/j.eti.2023.103228]

5. Hu, Y.H.; Bai, Y.C.; Xu, H.J. Analysis of reasons for urban road collapse and prevention and control countermeasures in recent decade of China. Highway; 2016; 61, pp. 130-135.

6. Kuwano, R.; Horii, T.; Kohashi, H.; Yamauchi, K. Defects of sewer pipes causing cave-in’s in the road. Proceedings of the 5th International Symposium on New Technologies for Urban Safety of Mega Cities in Asia; Phuket, Thailand, 16–17 November 2006.

7. He, Y.; Zhu, D.Z.; Zhang, T.; Shao, Y.; Yu, T. Experimental observations on the initiation of sand-bed erosion by an upward water jet. J. Hydraul. Eng.; 2017; 143, 06017007. [DOI: https://dx.doi.org/10.1061/(ASCE)HY.1943-7900.0001302]

8. Tang, Y.; Chan, D.H.; Zhu, D.Z. Numerical investigation of sand-bed erosion by an upward water jet. J. Eng. Mech.; 2017; 143, 04017104. [DOI: https://dx.doi.org/10.1061/(ASCE)EM.1943-7889.0001319]

9. Wen, S.J.; Cheng, W.C.; Li, D.F.; Hu, W.L. Evaluating gas breakthrough pressure and gas permeability in a landfill cover layer for mitigation of hazardous gas emissions. J. Environ. Manag.; 2023; 336, 117617. [DOI: https://dx.doi.org/10.1016/j.jenvman.2023.117617]

10. Greyvenstein, B.; van Zyl, J.E. An experimental investigation into the pressure leakage relationship of some failed water pipes. Aqua-J. Water Supply; 2007; 56, pp. 117-124. [DOI: https://dx.doi.org/10.2166/aqua.2007.065]

11. Davies, J.; Clarke, B.; Whiter, J.; Cunningham, R. Factors influencing the structural deterioration and collapse of rigid sewer pipes. Urban Water; 2001; 3, pp. 73-89. [DOI: https://dx.doi.org/10.1016/S1462-0758(01)00017-6]

12. Guo, S.; Zhang, T.; Zhang, Y.; Zhu, D.Z. An approximate solution for two-dimensional groundwater infiltration in sewer systems. Water Sci. Technol.; 2013; 67, pp. 347-352. [DOI: https://dx.doi.org/10.2166/wst.2012.568] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23168634]

13. Lambert, A. What do we know about pressure-leakage relation-ships in distribution systems? In Proceedings of the IWA Conference on Systems Approach to Leakage Control and Water Distribution System Management, Citeseer, Brno, Czech Republic, 16–18 May 2001.

14. Ssozi, E.; Reddy, B.; van Zyl, J. Numerical investigation of the influence of viscoelastic deformation on the pressure-leakage behavior of plastic pipes. J. Hydraul. Eng.; 2015; 142, 04015057. [DOI: https://dx.doi.org/10.1061/(ASCE)HY.1943-7900.0001095]

15. Yang, Y.; Zhu, D.Z.; Zhang, T.; Liu, W.; Guo, S. Improved model for contaminant intrusion induced by negative pressure events in water distribution systems. J. Hydraul. Eng.; 2016; 142, 06016012. [DOI: https://dx.doi.org/10.1061/(ASCE)HY.1943-7900.0001176]

16. Guo, S.; Shao, Y.; Zhang, T.; Zhu, D.; Zhang, Y. Physical Modeling on Sand Erosion around Defective Sewer Pipes under the Influence of Groundwater. J. Hydraul. Eng.; 2013; 139, pp. 1247-1257. [DOI: https://dx.doi.org/10.1061/(ASCE)HY.1943-7900.0000785]

17. Serpente, P.E. Understanding the modes of failure for sewer. Urban Drainage Rehabilitation Programs and Techique Selected Papers on Urban Drainage Rehabilitation from 1988–1993; Macaitis, W.A. ASCE: New York, NY, USA, 1994.

18. Samani, S.; Ye, M.; Zhang, F.; Pei, Y.Z.; Tang, G.P.; Elshall, A.; Moghaddam, A.A. Impacts of prior parameter distributions on Bayesian evaluation of groundwater model complexity. Water Sci. Eng.; 2018; 11, pp. 89-100. [DOI: https://dx.doi.org/10.1016/j.wse.2018.06.001]

19. Samani, S.; Moghaddam, A.A.; Ye, M. Investigating the effect of complexity on groundwater flow modeling uncertainty. Stoch. Env. Res. Risk. A; 2018; 32, pp. 643-659. [DOI: https://dx.doi.org/10.1007/s00477-017-1436-6]

20. Bertrand-Krajewski, J.; Cardoso, M.; Ellis, B.; Frehmann, T.; Giulianelli, M.; Gujer, W.; Krebs, P.; Pliska, J.; Pryl, K. Towards a better knowledge and management of infiltration and exfiltration in sewer systems: The APUSS project. Water Pract. Technol.; 2006; 1, wpt2006021. [DOI: https://dx.doi.org/10.2166/wpt.2006.021]

21. Cardoso, A.; Prigiobbe, V.; Giulianelli, M.; Baer, E.; Bénédittis, J.D.; Coelho, S.T. Assessing the impact of infiltration and exfiltration in sewer systems using performance indicators: Case studies of the APUSS project. Water Pract. Technol.; 2006; 1, wpt2006019. [DOI: https://dx.doi.org/10.2166/wpt.2006019]

22. Karpf, C.; Hoeft, S.; Scheffer, C.; Fuchs, L.; Krebs, P. Groundwater infiltration, surface water inflow and sewerage exfiltration considering hydrodynamic conditions in sewer systems. Water Sci. Technol.; 2011; 63, pp. 1841-1848. [DOI: https://dx.doi.org/10.2166/wst.2011.388]

23. Meguid, M.A.; Dang, H.K. The effect of erosion voids on existing tunnel linings. Tunn. Undergr. Space Technol.; 2009; 24, pp. 278-286. [DOI: https://dx.doi.org/10.1016/j.tust.2008.09.002]

24. Fenner, R.A. The influence of sewer bedding arrangements on infiltration rates and soil migration. Munic. Eng.; 1991; 8, pp. 105-117.

25. Rogers, C.J. Sewer Deterioration Studies: The Background to the Structural Assessment Procedure in the Sewerage Rehabilitation Mannual; 2nd ed. Water Research Centre: London, UK, 1986.

26. Abraham, D.M.; Wirahadikusumah, R. Development of prediction models for sewer deterioration. Proceedings of the 8th Conference on the Durability of Building Materials and Components; Vancouver, BC, Canada, 30 May–3 June 1999; pp. 1257-1267.

27. Moore, I. Assessment of damage to rigid sewer pipes and erosion voids in the soil, and implications for design of liners. Proceedings of the 2008 No-Dig Conference & Exhibition, North American Society for Trenchless Technology; Dallas, TX, USA, 3–6 June 2008.

28. Kloss, C.; Goniva, C.; Hager, A.; Amberger, S.; Pirker, S. Models, algorithms and validation for opensource DEM and CFD–DEM. Progress in Computational Fluid Dynamics. Prog. Comput. Fluid Dyn.; 2012; 12, pp. 140-152. [DOI: https://dx.doi.org/10.1504/PCFD.2012.047457]

29. Li, L.; Gu, Z.; Xu, W.; Tan, Y.; Fan, X.; Tan, D. Mixing mass transfer mechanism and dynamic control of gas-liquid-solid multiphase flow based on VOF-DEM coupling. Energy; 2023; 272, 127015. [DOI: https://dx.doi.org/10.1016/j.energy.2023.127015]

30. Li, L.; Xu, W.; Tan, Y.; Yang, Y.; Yang, J.; Tan, D. Fluid-induced vibration evolution mechanism of multiphase free sink vortex and the multi-source vibration sensing method. Mech. Syst. Signal. Pract.; 2023; 189, 110058. [DOI: https://dx.doi.org/10.1016/j.ymssp.2022.110058]

31. Nigmetova, A.; Masi, E.; Simonin, O.; Dufresne, Y.; Moureau, V. Three-dimensional DEM-CFD simulation of a lab-scale fluidized bed to support the development of two-fluid model approach. Int. J. Multiphas. Flow; 2022; 156, 104189. [DOI: https://dx.doi.org/10.1016/j.ijmultiphaseflow.2022.104189]

32. Li, D.; Nian, T.; Tiong, R.L.K.; Shen, Y.; Shao, Z. River blockage and impulse wave evolution of the Baige landslide in October 2018: Insights from coupled DEM-CFD analyses. Eng. Geol.; 2023; 321, 107169. [DOI: https://dx.doi.org/10.1016/j.enggeo.2023.107169]

33. Yan, W.; Cui, W.; Qi, L. Simulation of underwater concrete movement in flowing water using DEM-CFD coupling method. Constr. Build. Mater.; 2022; 319, 126134. [DOI: https://dx.doi.org/10.1016/j.conbuildmat.2021.126134]

34. Krzaczek, M.; Tejchman, J. Hydraulic fracturing process in rocks–small-scale simulations with a novel fully coupled DEM/CFD-based thermo-hydro-mechanical approach. Eng. Fract. Mech.; 2023; 289, 109424. [DOI: https://dx.doi.org/10.1016/j.engfracmech.2023.109424]

35. Rao, K.K.; Nott, P.R. An Introduction to Granular Flow; Cambridge University Press: New York, NY, USA, 2008.

36. Guo, Y.; Yu, X.B. Comparison of the implementation of three common types of coupled CFD-DEM model for simulating soil surface erosion. Int. J. Multiph. Flow; 2017; 91, pp. 89-100. [DOI: https://dx.doi.org/10.1016/j.ijmultiphaseflow.2017.01.006]

37. Ibrahim, A.; Meguid, M.A. CFD-DEM simulation of sand erosion into defective gravity pipes under constant groundwater table. Tunn. Undergr. Space Technol.; 2023; 131, 104823. [DOI: https://dx.doi.org/10.1016/j.tust.2022.104823]

38. Shi, C.; Zhang, Q.; Wang, S.N. Numerical Simulation Technology and Application with Particle Flow Code(PFC5.0); 1st ed. China Architecture & Building Press: Beijing, China, 2018.

39. Li, X.J.; Lu, Y.; Wu, Y.J. Application of CFD-DEM coupling method in seepage-induced ground subsidence. J. Shanghai Univ. (Nat. Sci. Ed.); 2020; 26, pp. 842-852.

40. Wang, Y.L.; Lu, Y. CFD-DEM simulation of disaster mechanism of pavement collapse caused by pipeline leakage in urban areas. J. Civil. Environ. Eng.; 2021; 43, pp. 60-67.