# On the External Failure Surface in PSEs: Numerical and Theoretical Methods

Authors: Jiaqing Liu; Zeyu Xu; Shuai Zhou
Journal: Advances in Civil Engineering
 (2022)

---

## Abstract

The evolution of soil arching, involving internal and external failure surfaces, is of significance to the load transfer mechanism of pile-supported embankments (PSEs), and external failure surfaces are generally observed in cases of greater embankment height. In this study, the evolution of external failure surfaces was investigated by both numerical and theoretical methods. To begin with, numerical simulations of trapdoor tests were carried out by a two-dimensional discrete element method. The influences of two key parameters (i.e., embankment height and net spacing between piles) on the development of external failure surfaces were emphasized. The measured coefficient of lateral earth pressure around the external failure surface was close to the predicted coefficient of active earth pressure. Then, a theoretical solution considering inclined external failure surfaces, more realistic compared to those adopted in Terzaghi’s method, was proposed. Compared with Terzaghi’s method, the proposed solution exhibited better consistency with laboratory observations, especially when external failure surfaces were significantly observed.


---

## Body

## 1. Introduction

Pile-supported embankments have been widely used to treat the highly compressible ground in the construction of expressways [1–3] and high-speed railways [4, 5] owing to their superiority in reducing postconstruction settlement and shortening the construction period. Soil arching is the most fundamental load transfer mechanism of the pile-supported embankment. To analyze the soil arching effect, numerous analysis models have been proposed, and generally, they can be divided into three groups. The first group is based on the Marston theory [6] and the Terzaghi theory [7], considering the equilibrium of the soil mass. The second group considers the limit equilibrium theory of soil arch [8]. The third group [9] assumes a fixed arch in the embankment, where the embankment loads above the arch are directly transmitted to the pile and the embankment loads under the arch are applied to the soft soil. However, the results of the same project using different models vary greatly, indicating that the analytical models heavily rely on the assumed patterns of arching evolution.The mechanism of arching evolution has been commonly investigated by trapdoor tests [10–14], originally designed by Terzaghi [7]. By this means, soil arching is generally recognized to be evolved with various key parameters, such as embankment height (h), pile diameter (b), and pile spacing (s) [15–19]. In some cases, the failure surface may extend beyond the trapdoor to form an external failure surface. Since vertical failure surfaces were assumed for soil-pile stress distribution in Terzaghi’s [7] were oversimplified, a series of trapdoor tests were performed by Ladanyi and Hoyaux [20] to check Terzaghi’s method. A great discrepancy between theoretical results and laboratory observations can be observed when the external failure surfaces are formed at a great magnitude of trapdoor displacement. For this reason, different failure surfaces in embankments, especially the external failure surface at great trapdoor displacement, can be further investigated.Several experimental investigations have been carried out on the external failure surface [21–24]. In particular, Costa et al. [21] reported that the inclination of external failure surfaces to the horizontal would not be less than the repose angle of the soil, which represents the angle of shearing resistance of soil at its loosest state. Theoretical models [25–27] with sliding surfaces of arbitrary inclinations are also proposed to examine the effect of the sliding surface on vertical stress. Although numerous numerical simulations have been presented on soil arching in PSEs [18, 19, 28, 29], the discrete element method (DEM) has scarcely been adopted to study the microscopic characteristics of the external failure surfaces.In the present study, trapdoor tests were simulated by the particle flow software PFC2D, coded by the discrete element method (DEM) [30]. The DEM can provide better insight into the microbehavior of the embankment as a result of the discrete nature of the granular embankment fill materials. On this basis, more attention will be paid to the influences of embankment height (h) and pile net spacing (s-b) on the evolution of external failure surfaces and the stress state near the external failure surfaces. In parallel, a simple theoretical model considering the external failure surfaces would be proposed, and then it would be compared with Terzaghi’s method.


## 2. The Numerical Model

### 2.1. Description of Reference Case

A numerical model by PFC 2D has been established in this study based on the laboratory tests by Jenck et al. [31] and the configuration of the laboratory model was shown in Figure 1. The embankment fill is represented by steel rods of 3, 4, and 5 mm in diameter (1 : 1:1 by volume) and the porosity of the assembly of steel bars is equal to 0.18. The unit weight and internal friction angle of the assembly are 62 kN/m3 and 24°, respectively. The rigid pile and the soft soil are simulated by rigid elements and foam, respectively.Figure 1
Configuration diagram of the reference case (dimensions in mm).


### 2.2. Numerical Simulation

#### 2.2.1. Basic Kinetics Used in PFC

The sphere (particle) element and the wall element are the basic components used in the PFC model, and the interactions between them are represented by the contact force. The contact force applied to the particles will produce corresponding particle movement as per the motion equation; however, the one applied to the wall would not activate the motion of the wall.Without the constraint of deformation coordination applied between particles, the equilibrium equation governs the behavior of the PFC model. If the resultant force or moment of the particles is not zero, the unbalanced force or moment will lead to the movement of particles according to Newton’s law. The movement of particles is thus dependent on the contact forces of adjacent particles, and the relationship between contact force and resultant displacement is called the physical equation. The calculation of PFC is based on Newton's second law and force-displacement law, where Newton’s second law applied to each particle is used to update the position of the particle and the wall. The force-displacement law is used to update the contact force among particles at each contact. Each particle will undergo a step-by-step iteration with updated position and contact force as per the kinetics, and the iterations will be terminated when the force equilibrium is reached.(1)
Force-displacement law (physical equation)The force-displacement law is applied at the contact to update the contact force. The contact force, however, only works when it is less than the activity distance, and the contact gap between particles should be identified first. The contact forceFi can be decomposed into the normal component Fin perpendicular to the contact plane and the tangential component force Fis in the contact plane.(2)
Newton’s second law (equation of motion)If the resultant force or moment acting on the particles is not zero, translational motion or rotation will be produced accordingly. The translational motion of the particle is described by the linear velocity of the particle center, while the rotation of the particle is described by the angular velocity. The above motions are governed by the following equations:(1)Fi=mx¨−gi,(2)Mi=Iω˙i,where Fi is the unbalanced resultant force acting on the particles; m is the particle mass; x¨ is the acceleration of the particle; gi is the gravitational acceleration; Mi is the unbalanced resultant moment on the particles; I is the moment of inertia of the sphere; and ω˙i is the angular acceleration.Subsequently, the accelerations of particles can be obtained by solving equations (1) and (2) at the time step, Δt, and are then used to update the velocity and angular velocity of the particles of the next time step.


#### 2.2.2. The Contact Constitutive Model

In PFC simulations, the macroscopic mechanical properties of materials are mainly governed by the contact between microparticles, or the particle contact constitutive model. The linear contact model is used in this study because it is applicable for most analysis scenarios and also by using this model the required microparameters can be easily obtained. The mechanical model of particle contact is shown in Figure2.Figure 2
A schematic diagram of the linear contact model.(3)kn=knAknBknA+knBks=ksAksBksA+ksB.Normal and tangential stiffness are, respectively, defined in PFC, and the stiffness corresponds with the particle stiffnesses acting in series. The normal stiffness,kn, and tangential (shear) stiffness, ks, can thus be estimated as follows: where knA and knB are the normal stiffness of particles A and B, respectively, whereas ksA and ksB are the tangential stiffness of particles A and B, respectively.


#### 2.2.3. Damping

In this model, the local damping is always active with the default value of 0.7.


#### 2.2.4. Servo Control

Because the stress cannot be directly endowed to activate the wall in PFC, confining pressure or loading is usually applied by displacing the wall through the servo control.


#### 2.2.5. Calibration of Embankment Fills

The linear-based model was adopted in the DEM simulation to represent the contact behaviors among embankment fills. Three major parameters were required for the linear contact model: the normal stiffness of particles,kn; the shear stiffness of particles, ks; and the friction coefficient, μ_s. Their magnitudes were first roughly estimated based on the recommendations by Stahl and Konietzky [32] and Han et al. [28], and then they were precisely calibrated by using a numerical biaxial test before modeling trapdoor tests. A calibration chamber with a size of 200 mm × 200 mm was built. According to the observation [26], samples were simulated by disk particles with three different diameters (3, 4, and 5 mm) at a volume ratio of 1 : 1:1 and a porosity of 0.18. The servo-controlled confining stress was applied to the embankment fill in the horizontal direction, and the vertical stress was controlled by displacing the horizontal wall at a constant velocity. As shown in Figure 3(a), the curves of deviatoric stress versus the axial strain are obtained at three different confining stresses (i.e., 20, 30, and 40 kPa). It can be observed that numerical simulations and experimental results were in good agreement. Shear stress versus normal stress is shown in Figure 3(b). The internal friction angle of the embankment fill φ can be calculated as follows:(4)sinϕ=σ1−σ3σ1+σ3,where σ1 and σ3 are maximum and minimum principal stresses, respectively. The value of the internal friction angle is 24.4°, which is close to the experimental parameters (i.e., 24°), and the other micromechanical parameters of the particles are listed in Table 1.Figure 3
Calibration of embankment fills: (a) curves of the deviatoric stress versus the axial strain and (b) the internal friction angle of the embankment fills.
(a)(b)Table 1
Micromechanical properties.
ParameterValueEmbankment fillsNormal stiffness of particles,kn2.25 × 107 N/mShear stiffness of particles,ks1.50 × 107 N/mFriction coefficient,μ_s0.7The density of the particle,ρ7600 kg/m3


#### 2.2.6. Numerical Filling Technique

Numerical simulation techniques to generate embankment fills have been improved gradually in recent years [28, 29, 33]. The embankments primitively generated by the “radius expansion” and “falling rain” methods were unrealistic. The distribution of the fills was nonuniform and the vertical stresses did not vary with the density gradient along the depth. To ensure the natural stress state, Bhandari and Han [33], and Han et al. [28] numerically compacted the embankment fills through the servo mechanism of a moving wall. Lai et al. [29] proposed a multilayer compaction method that can better simulate practical embankment filling processes and can guarantee the natural stress state. A new problem, however, would inherently arise if the embankment fills were compacted to the target porosity. Then, the density of the region near the moving wall would be denser than that of the region away from the wall. In this study, overcompaction and release were used to improve the multilayer compaction method. In this study, the initial porosity ni and target nr were 0.59 and 0.18, respectively, and the thickness of each compaction layer of the embankment was 100 mm, which was consistent with the compaction process in experimental tests [31]. The flowchart of the improved multilayer compaction method is shown in Figure 4, and the detailed steps are as follows.(1)
Particle generation. Particles are generated in an area (area =s×200 mm) according to the calculated porosity ni and the cycle has sufficient time for uniform distribution of particles.(2)
Overcompression. The particles are compressed by moving down the top platen to a height of less than 100 mm, and the particles in the area are thus compressed to be overlapped within a transitional porosityng smaller than nr.(3)
Normal compression. The platen is moved up to the designed position (h = 100 mm). The particles are then bounced off by the contact force between the overlapped particles.(4)
The second layer is generated by following the steps from 1 to 3.(5)
The platen between the layers is removed and it is cycled to the equilibrium.Figure 4
Flowchart of an improved multilayer compaction method.In this way, the embankment can be generated in lifts without gravity (g = 0) and friction (μs = 0).The embankment generated by this technique would be uniform, compacted, and can achieve the target porosity. After imposing gravity on the generated embankment, Figure4 illustrates the accuracy of the embankment stress state generated by the improved multilayer compaction method. By fitting the numerical prediction in Figure 5, the slope of the numerical prediction curve, indicating the unit weight used in the simulation is 64.33, which is close to the magnitude of the unit weight in the experiment (γ = 62 kN/m3).Figure 5
Verification of the embankment filling method.


#### 2.2.7. The Simulation Model for Reference Cases

A schematic diagram of the numerical model is shown in Figure6, where the setups are consistent with the experiments. The pile diameter (b) of 100 mm and three different pile spacing (s) of 320, 450, and 650 mm were considered in reference cases, and the embankment height (h) increased from 100 mm to 600 mm for each pile spacing. A movable wall located at the bottom of the embankment simulates the trapdoor, and its movement was denoted as trapdoor displacement (δ). In the numerical simulation, the moving rate should be slow enough to guarantee the quasistatic equilibrium of the sample during the test [34, 35], and a constant vertical velocity of 1 mm/s was applied to the trapdoor. In this model, the default time step for each calculation cycle was set as 2.1 × 10−5 s and the corresponding moving rate was equivalent to 2.1 × 10−5 mm/step, so that it takes more than 46,000 steps to move 1 mm. At the same time, a series of measurement circles were set to evaluate the performance of the embankment.Figure 6
A schematic diagram of the numerical model (dimensions in mm).


### 2.3. Validation of the Numerical Model

The efficacy of loads distribution,E, is used to characterize the soil arching effect, and it can be calculated as follows [8]:(5)E=FpW,where Fp is the load carried by piles and W is the total weight of the embankment. The results of the DEM simulations were compared with the experimental results for validation. As shown in Figure 7, the efficacy E increases with the embankment height and decreases with the increase of pile-net spacing. It can be found that the results of the DEM simulations match well with the experimental results.Figure 7
Variation of efficacy with embankment height for various pile-net spacings.


### 2.4. Numerical Results

#### 2.4.1. Evolution of External Failure Surfaces

The evolution of soil arching can be described by the development of failure surfaces, which are generally explored by the displacement contour [18, 19, 21, 36], and a typical case is illustrated in Figure 7 (s = 320 mm, h = 400 mm). Figure 8(a) demonstrates the development of soil displacement during the displacement of the trapdoor. When the trapdoor reaches the designated displacement of 17.25 mm, the embankment is divided into three regions according to the displacement contour, which includes the active region, transition region, and stable region. According to Costa et al. [21], the boundary between the transition region and active region (curve OA) represents the internal surface and the boundary between the transition region and stable region (curve OB) represents the external surface. The active region will flow downward with the displacement of the trapdoor, making the transition region move obliquely downward to the active region accordingly.Figure 8
Displacement contour of case1A (s = 320 mm, (h = 400 mm): (a) development of soil displacement and (b) zone division of final displacement contour.
(a)(b)In this study, the failure surfaces are outlined in Figure9. To facilitate the description of the failure surface, simplified straight failure surfaces were adopted to approximately represent the actual failure surfaces, and they share the same starting and ending points. The inclination of the simplified failure surface, θ, is the angle between the line of the failure surface and the horizontal plane, as shown in Figure 9. The external failure surfaces are only observed when the embankment height, h, is larger than the pile-net spacing (s-b). For example, for cases with a net spacing of 220 mm, the external failure surface would occur when the embankment heights are 300, 400, 500, and 600 mm. It is also true for the cases with the net spacing of 350 mm when embankment heights are 400, 500, and 600 mm, and the cases with the net spacing of 550 mm when the embankment height is 600 mm. It should be noted that the inclination of the failure surface (θ) decreases gradually with the embankment height.Figure 9
Evolution of the failure surfaces with various embankment height for different pile-net spacing: (a) (s)-(b) = 220 mm, (b) (s)-(b) = 350 mm, and (c) (s)-(b) = 550 mm.
(a)(b)(c)


#### 2.4.2. Distribution of Vertical Stress

The contact force between particles can be monitored by measurement circles and then can be converted into continuum stress based on the following equation:(6)σ¯ij=1−n∑Np1Sp∑Np1∑Ncsicsipnic,pFjc,where Np1 is the number of particles whose centroids are contained within the measurement circle; Nc is the number of contacts between particles; n is the porosity within the measurement circle; S(p) is the area of a particle; si(p) and si(c) are the locations of a particle centroid and its contact, respectively; ni(c,p) is the unit normal vector directed from a particle centroid to its contact location; Fj(c) is the force acting on the contact; and indices i and j are in the set of (x, z).The vertical stress,σv, can be calculated based on the following equation for i = z and j = z:(7)σv=σzz,where σzz is the average stress in the direction of the z-axis.Figure10 shows the distribution of vertical stress in the embankment before and after displacing the trapdoor. The stress redistribution results in the prominent increase of the loads acting on piles increasing from 25 kPa to 60 kPa, but the reduction of loads on the soil decreasing from 25 kPa to 9.5 kPa with the trapdoor displacement. The areas of increased vertical stress and corresponding equipotential lines of the contour roughly coincide with the stable regions and external failure surfaces in Figure 8, respectively. It implies that more loads of surrounding soil are transferred to the stable regions through the external failure surface, thus increasing the vertical stress on the pile top.Figure 10
Vertical stress contour: (a) before trapdoor displacement and (b) after trapdoor displacement.
(a)(b)The ratio of the horizontal stress to vertical stress, defined as the coefficient of lateral earth pressure (k), is an important parameter and it is calculated before and after displacing the trapdoor as shown in Figure 11. The coefficient of lateral earth pressure around the external failure surface is 0.48, roughly approaching the coefficient of active earth pressure, ka = tan2(45°-φ/2) = 0.42. On this basis, the coefficient of active earth pressure can be used to represent the coefficient of lateral earth pressure around the external failure surface in the proposed method.Figure 11
Coefficient of lateral earth pressure: (a) before trapdoor displacement and (b) after trapdoor displacement.
(a)(b)


## 2.1. Description of Reference Case

A numerical model by PFC 2D has been established in this study based on the laboratory tests by Jenck et al. [31] and the configuration of the laboratory model was shown in Figure 1. The embankment fill is represented by steel rods of 3, 4, and 5 mm in diameter (1 : 1:1 by volume) and the porosity of the assembly of steel bars is equal to 0.18. The unit weight and internal friction angle of the assembly are 62 kN/m3 and 24°, respectively. The rigid pile and the soft soil are simulated by rigid elements and foam, respectively.Figure 1
Configuration diagram of the reference case (dimensions in mm).


## 2.2. Numerical Simulation

### 2.2.1. Basic Kinetics Used in PFC

The sphere (particle) element and the wall element are the basic components used in the PFC model, and the interactions between them are represented by the contact force. The contact force applied to the particles will produce corresponding particle movement as per the motion equation; however, the one applied to the wall would not activate the motion of the wall.Without the constraint of deformation coordination applied between particles, the equilibrium equation governs the behavior of the PFC model. If the resultant force or moment of the particles is not zero, the unbalanced force or moment will lead to the movement of particles according to Newton’s law. The movement of particles is thus dependent on the contact forces of adjacent particles, and the relationship between contact force and resultant displacement is called the physical equation. The calculation of PFC is based on Newton's second law and force-displacement law, where Newton’s second law applied to each particle is used to update the position of the particle and the wall. The force-displacement law is used to update the contact force among particles at each contact. Each particle will undergo a step-by-step iteration with updated position and contact force as per the kinetics, and the iterations will be terminated when the force equilibrium is reached.(1)
Force-displacement law (physical equation)The force-displacement law is applied at the contact to update the contact force. The contact force, however, only works when it is less than the activity distance, and the contact gap between particles should be identified first. The contact forceFi can be decomposed into the normal component Fin perpendicular to the contact plane and the tangential component force Fis in the contact plane.(2)
Newton’s second law (equation of motion)If the resultant force or moment acting on the particles is not zero, translational motion or rotation will be produced accordingly. The translational motion of the particle is described by the linear velocity of the particle center, while the rotation of the particle is described by the angular velocity. The above motions are governed by the following equations:(1)Fi=mx¨−gi,(2)Mi=Iω˙i,where Fi is the unbalanced resultant force acting on the particles; m is the particle mass; x¨ is the acceleration of the particle; gi is the gravitational acceleration; Mi is the unbalanced resultant moment on the particles; I is the moment of inertia of the sphere; and ω˙i is the angular acceleration.Subsequently, the accelerations of particles can be obtained by solving equations (1) and (2) at the time step, Δt, and are then used to update the velocity and angular velocity of the particles of the next time step.


### 2.2.2. The Contact Constitutive Model

In PFC simulations, the macroscopic mechanical properties of materials are mainly governed by the contact between microparticles, or the particle contact constitutive model. The linear contact model is used in this study because it is applicable for most analysis scenarios and also by using this model the required microparameters can be easily obtained. The mechanical model of particle contact is shown in Figure2.Figure 2
A schematic diagram of the linear contact model.(3)kn=knAknBknA+knBks=ksAksBksA+ksB.Normal and tangential stiffness are, respectively, defined in PFC, and the stiffness corresponds with the particle stiffnesses acting in series. The normal stiffness,kn, and tangential (shear) stiffness, ks, can thus be estimated as follows: where knA and knB are the normal stiffness of particles A and B, respectively, whereas ksA and ksB are the tangential stiffness of particles A and B, respectively.


### 2.2.3. Damping

In this model, the local damping is always active with the default value of 0.7.


### 2.2.4. Servo Control

Because the stress cannot be directly endowed to activate the wall in PFC, confining pressure or loading is usually applied by displacing the wall through the servo control.


### 2.2.5. Calibration of Embankment Fills

The linear-based model was adopted in the DEM simulation to represent the contact behaviors among embankment fills. Three major parameters were required for the linear contact model: the normal stiffness of particles,kn; the shear stiffness of particles, ks; and the friction coefficient, μ_s. Their magnitudes were first roughly estimated based on the recommendations by Stahl and Konietzky [32] and Han et al. [28], and then they were precisely calibrated by using a numerical biaxial test before modeling trapdoor tests. A calibration chamber with a size of 200 mm × 200 mm was built. According to the observation [26], samples were simulated by disk particles with three different diameters (3, 4, and 5 mm) at a volume ratio of 1 : 1:1 and a porosity of 0.18. The servo-controlled confining stress was applied to the embankment fill in the horizontal direction, and the vertical stress was controlled by displacing the horizontal wall at a constant velocity. As shown in Figure 3(a), the curves of deviatoric stress versus the axial strain are obtained at three different confining stresses (i.e., 20, 30, and 40 kPa). It can be observed that numerical simulations and experimental results were in good agreement. Shear stress versus normal stress is shown in Figure 3(b). The internal friction angle of the embankment fill φ can be calculated as follows:(4)sinϕ=σ1−σ3σ1+σ3,where σ1 and σ3 are maximum and minimum principal stresses, respectively. The value of the internal friction angle is 24.4°, which is close to the experimental parameters (i.e., 24°), and the other micromechanical parameters of the particles are listed in Table 1.Figure 3
Calibration of embankment fills: (a) curves of the deviatoric stress versus the axial strain and (b) the internal friction angle of the embankment fills.
(a)(b)Table 1
Micromechanical properties.
ParameterValueEmbankment fillsNormal stiffness of particles,kn2.25 × 107 N/mShear stiffness of particles,ks1.50 × 107 N/mFriction coefficient,μ_s0.7The density of the particle,ρ7600 kg/m3


### 2.2.6. Numerical Filling Technique

Numerical simulation techniques to generate embankment fills have been improved gradually in recent years [28, 29, 33]. The embankments primitively generated by the “radius expansion” and “falling rain” methods were unrealistic. The distribution of the fills was nonuniform and the vertical stresses did not vary with the density gradient along the depth. To ensure the natural stress state, Bhandari and Han [33], and Han et al. [28] numerically compacted the embankment fills through the servo mechanism of a moving wall. Lai et al. [29] proposed a multilayer compaction method that can better simulate practical embankment filling processes and can guarantee the natural stress state. A new problem, however, would inherently arise if the embankment fills were compacted to the target porosity. Then, the density of the region near the moving wall would be denser than that of the region away from the wall. In this study, overcompaction and release were used to improve the multilayer compaction method. In this study, the initial porosity ni and target nr were 0.59 and 0.18, respectively, and the thickness of each compaction layer of the embankment was 100 mm, which was consistent with the compaction process in experimental tests [31]. The flowchart of the improved multilayer compaction method is shown in Figure 4, and the detailed steps are as follows.(1)
Particle generation. Particles are generated in an area (area =s×200 mm) according to the calculated porosity ni and the cycle has sufficient time for uniform distribution of particles.(2)
Overcompression. The particles are compressed by moving down the top platen to a height of less than 100 mm, and the particles in the area are thus compressed to be overlapped within a transitional porosityng smaller than nr.(3)
Normal compression. The platen is moved up to the designed position (h = 100 mm). The particles are then bounced off by the contact force between the overlapped particles.(4)
The second layer is generated by following the steps from 1 to 3.(5)
The platen between the layers is removed and it is cycled to the equilibrium.Figure 4
Flowchart of an improved multilayer compaction method.In this way, the embankment can be generated in lifts without gravity (g = 0) and friction (μs = 0).The embankment generated by this technique would be uniform, compacted, and can achieve the target porosity. After imposing gravity on the generated embankment, Figure4 illustrates the accuracy of the embankment stress state generated by the improved multilayer compaction method. By fitting the numerical prediction in Figure 5, the slope of the numerical prediction curve, indicating the unit weight used in the simulation is 64.33, which is close to the magnitude of the unit weight in the experiment (γ = 62 kN/m3).Figure 5
Verification of the embankment filling method.


### 2.2.7. The Simulation Model for Reference Cases

A schematic diagram of the numerical model is shown in Figure6, where the setups are consistent with the experiments. The pile diameter (b) of 100 mm and three different pile spacing (s) of 320, 450, and 650 mm were considered in reference cases, and the embankment height (h) increased from 100 mm to 600 mm for each pile spacing. A movable wall located at the bottom of the embankment simulates the trapdoor, and its movement was denoted as trapdoor displacement (δ). In the numerical simulation, the moving rate should be slow enough to guarantee the quasistatic equilibrium of the sample during the test [34, 35], and a constant vertical velocity of 1 mm/s was applied to the trapdoor. In this model, the default time step for each calculation cycle was set as 2.1 × 10−5 s and the corresponding moving rate was equivalent to 2.1 × 10−5 mm/step, so that it takes more than 46,000 steps to move 1 mm. At the same time, a series of measurement circles were set to evaluate the performance of the embankment.Figure 6
A schematic diagram of the numerical model (dimensions in mm).


## 2.2.1. Basic Kinetics Used in PFC

The sphere (particle) element and the wall element are the basic components used in the PFC model, and the interactions between them are represented by the contact force. The contact force applied to the particles will produce corresponding particle movement as per the motion equation; however, the one applied to the wall would not activate the motion of the wall.Without the constraint of deformation coordination applied between particles, the equilibrium equation governs the behavior of the PFC model. If the resultant force or moment of the particles is not zero, the unbalanced force or moment will lead to the movement of particles according to Newton’s law. The movement of particles is thus dependent on the contact forces of adjacent particles, and the relationship between contact force and resultant displacement is called the physical equation. The calculation of PFC is based on Newton's second law and force-displacement law, where Newton’s second law applied to each particle is used to update the position of the particle and the wall. The force-displacement law is used to update the contact force among particles at each contact. Each particle will undergo a step-by-step iteration with updated position and contact force as per the kinetics, and the iterations will be terminated when the force equilibrium is reached.(1)
Force-displacement law (physical equation)The force-displacement law is applied at the contact to update the contact force. The contact force, however, only works when it is less than the activity distance, and the contact gap between particles should be identified first. The contact forceFi can be decomposed into the normal component Fin perpendicular to the contact plane and the tangential component force Fis in the contact plane.(2)
Newton’s second law (equation of motion)If the resultant force or moment acting on the particles is not zero, translational motion or rotation will be produced accordingly. The translational motion of the particle is described by the linear velocity of the particle center, while the rotation of the particle is described by the angular velocity. The above motions are governed by the following equations:(1)Fi=mx¨−gi,(2)Mi=Iω˙i,where Fi is the unbalanced resultant force acting on the particles; m is the particle mass; x¨ is the acceleration of the particle; gi is the gravitational acceleration; Mi is the unbalanced resultant moment on the particles; I is the moment of inertia of the sphere; and ω˙i is the angular acceleration.Subsequently, the accelerations of particles can be obtained by solving equations (1) and (2) at the time step, Δt, and are then used to update the velocity and angular velocity of the particles of the next time step.


## 2.2.2. The Contact Constitutive Model

In PFC simulations, the macroscopic mechanical properties of materials are mainly governed by the contact between microparticles, or the particle contact constitutive model. The linear contact model is used in this study because it is applicable for most analysis scenarios and also by using this model the required microparameters can be easily obtained. The mechanical model of particle contact is shown in Figure2.Figure 2
A schematic diagram of the linear contact model.(3)kn=knAknBknA+knBks=ksAksBksA+ksB.Normal and tangential stiffness are, respectively, defined in PFC, and the stiffness corresponds with the particle stiffnesses acting in series. The normal stiffness,kn, and tangential (shear) stiffness, ks, can thus be estimated as follows: where knA and knB are the normal stiffness of particles A and B, respectively, whereas ksA and ksB are the tangential stiffness of particles A and B, respectively.


## 2.2.3. Damping

In this model, the local damping is always active with the default value of 0.7.


## 2.2.4. Servo Control

Because the stress cannot be directly endowed to activate the wall in PFC, confining pressure or loading is usually applied by displacing the wall through the servo control.


## 2.2.5. Calibration of Embankment Fills

The linear-based model was adopted in the DEM simulation to represent the contact behaviors among embankment fills. Three major parameters were required for the linear contact model: the normal stiffness of particles,kn; the shear stiffness of particles, ks; and the friction coefficient, μ_s. Their magnitudes were first roughly estimated based on the recommendations by Stahl and Konietzky [32] and Han et al. [28], and then they were precisely calibrated by using a numerical biaxial test before modeling trapdoor tests. A calibration chamber with a size of 200 mm × 200 mm was built. According to the observation [26], samples were simulated by disk particles with three different diameters (3, 4, and 5 mm) at a volume ratio of 1 : 1:1 and a porosity of 0.18. The servo-controlled confining stress was applied to the embankment fill in the horizontal direction, and the vertical stress was controlled by displacing the horizontal wall at a constant velocity. As shown in Figure 3(a), the curves of deviatoric stress versus the axial strain are obtained at three different confining stresses (i.e., 20, 30, and 40 kPa). It can be observed that numerical simulations and experimental results were in good agreement. Shear stress versus normal stress is shown in Figure 3(b). The internal friction angle of the embankment fill φ can be calculated as follows:(4)sinϕ=σ1−σ3σ1+σ3,where σ1 and σ3 are maximum and minimum principal stresses, respectively. The value of the internal friction angle is 24.4°, which is close to the experimental parameters (i.e., 24°), and the other micromechanical parameters of the particles are listed in Table 1.Figure 3
Calibration of embankment fills: (a) curves of the deviatoric stress versus the axial strain and (b) the internal friction angle of the embankment fills.
(a)(b)Table 1
Micromechanical properties.
ParameterValueEmbankment fillsNormal stiffness of particles,kn2.25 × 107 N/mShear stiffness of particles,ks1.50 × 107 N/mFriction coefficient,μ_s0.7The density of the particle,ρ7600 kg/m3


## 2.2.6. Numerical Filling Technique

Numerical simulation techniques to generate embankment fills have been improved gradually in recent years [28, 29, 33]. The embankments primitively generated by the “radius expansion” and “falling rain” methods were unrealistic. The distribution of the fills was nonuniform and the vertical stresses did not vary with the density gradient along the depth. To ensure the natural stress state, Bhandari and Han [33], and Han et al. [28] numerically compacted the embankment fills through the servo mechanism of a moving wall. Lai et al. [29] proposed a multilayer compaction method that can better simulate practical embankment filling processes and can guarantee the natural stress state. A new problem, however, would inherently arise if the embankment fills were compacted to the target porosity. Then, the density of the region near the moving wall would be denser than that of the region away from the wall. In this study, overcompaction and release were used to improve the multilayer compaction method. In this study, the initial porosity ni and target nr were 0.59 and 0.18, respectively, and the thickness of each compaction layer of the embankment was 100 mm, which was consistent with the compaction process in experimental tests [31]. The flowchart of the improved multilayer compaction method is shown in Figure 4, and the detailed steps are as follows.(1)
Particle generation. Particles are generated in an area (area =s×200 mm) according to the calculated porosity ni and the cycle has sufficient time for uniform distribution of particles.(2)
Overcompression. The particles are compressed by moving down the top platen to a height of less than 100 mm, and the particles in the area are thus compressed to be overlapped within a transitional porosityng smaller than nr.(3)
Normal compression. The platen is moved up to the designed position (h = 100 mm). The particles are then bounced off by the contact force between the overlapped particles.(4)
The second layer is generated by following the steps from 1 to 3.(5)
The platen between the layers is removed and it is cycled to the equilibrium.Figure 4
Flowchart of an improved multilayer compaction method.In this way, the embankment can be generated in lifts without gravity (g = 0) and friction (μs = 0).The embankment generated by this technique would be uniform, compacted, and can achieve the target porosity. After imposing gravity on the generated embankment, Figure4 illustrates the accuracy of the embankment stress state generated by the improved multilayer compaction method. By fitting the numerical prediction in Figure 5, the slope of the numerical prediction curve, indicating the unit weight used in the simulation is 64.33, which is close to the magnitude of the unit weight in the experiment (γ = 62 kN/m3).Figure 5
Verification of the embankment filling method.


## 2.2.7. The Simulation Model for Reference Cases

A schematic diagram of the numerical model is shown in Figure6, where the setups are consistent with the experiments. The pile diameter (b) of 100 mm and three different pile spacing (s) of 320, 450, and 650 mm were considered in reference cases, and the embankment height (h) increased from 100 mm to 600 mm for each pile spacing. A movable wall located at the bottom of the embankment simulates the trapdoor, and its movement was denoted as trapdoor displacement (δ). In the numerical simulation, the moving rate should be slow enough to guarantee the quasistatic equilibrium of the sample during the test [34, 35], and a constant vertical velocity of 1 mm/s was applied to the trapdoor. In this model, the default time step for each calculation cycle was set as 2.1 × 10−5 s and the corresponding moving rate was equivalent to 2.1 × 10−5 mm/step, so that it takes more than 46,000 steps to move 1 mm. At the same time, a series of measurement circles were set to evaluate the performance of the embankment.Figure 6
A schematic diagram of the numerical model (dimensions in mm).


## 2.3. Validation of the Numerical Model

The efficacy of loads distribution,E, is used to characterize the soil arching effect, and it can be calculated as follows [8]:(5)E=FpW,where Fp is the load carried by piles and W is the total weight of the embankment. The results of the DEM simulations were compared with the experimental results for validation. As shown in Figure 7, the efficacy E increases with the embankment height and decreases with the increase of pile-net spacing. It can be found that the results of the DEM simulations match well with the experimental results.Figure 7
Variation of efficacy with embankment height for various pile-net spacings.


## 2.4. Numerical Results

### 2.4.1. Evolution of External Failure Surfaces

The evolution of soil arching can be described by the development of failure surfaces, which are generally explored by the displacement contour [18, 19, 21, 36], and a typical case is illustrated in Figure 7 (s = 320 mm, h = 400 mm). Figure 8(a) demonstrates the development of soil displacement during the displacement of the trapdoor. When the trapdoor reaches the designated displacement of 17.25 mm, the embankment is divided into three regions according to the displacement contour, which includes the active region, transition region, and stable region. According to Costa et al. [21], the boundary between the transition region and active region (curve OA) represents the internal surface and the boundary between the transition region and stable region (curve OB) represents the external surface. The active region will flow downward with the displacement of the trapdoor, making the transition region move obliquely downward to the active region accordingly.Figure 8
Displacement contour of case1A (s = 320 mm, (h = 400 mm): (a) development of soil displacement and (b) zone division of final displacement contour.
(a)(b)In this study, the failure surfaces are outlined in Figure9. To facilitate the description of the failure surface, simplified straight failure surfaces were adopted to approximately represent the actual failure surfaces, and they share the same starting and ending points. The inclination of the simplified failure surface, θ, is the angle between the line of the failure surface and the horizontal plane, as shown in Figure 9. The external failure surfaces are only observed when the embankment height, h, is larger than the pile-net spacing (s-b). For example, for cases with a net spacing of 220 mm, the external failure surface would occur when the embankment heights are 300, 400, 500, and 600 mm. It is also true for the cases with the net spacing of 350 mm when embankment heights are 400, 500, and 600 mm, and the cases with the net spacing of 550 mm when the embankment height is 600 mm. It should be noted that the inclination of the failure surface (θ) decreases gradually with the embankment height.Figure 9
Evolution of the failure surfaces with various embankment height for different pile-net spacing: (a) (s)-(b) = 220 mm, (b) (s)-(b) = 350 mm, and (c) (s)-(b) = 550 mm.
(a)(b)(c)


### 2.4.2. Distribution of Vertical Stress

The contact force between particles can be monitored by measurement circles and then can be converted into continuum stress based on the following equation:(6)σ¯ij=1−n∑Np1Sp∑Np1∑Ncsicsipnic,pFjc,where Np1 is the number of particles whose centroids are contained within the measurement circle; Nc is the number of contacts between particles; n is the porosity within the measurement circle; S(p) is the area of a particle; si(p) and si(c) are the locations of a particle centroid and its contact, respectively; ni(c,p) is the unit normal vector directed from a particle centroid to its contact location; Fj(c) is the force acting on the contact; and indices i and j are in the set of (x, z).The vertical stress,σv, can be calculated based on the following equation for i = z and j = z:(7)σv=σzz,where σzz is the average stress in the direction of the z-axis.Figure10 shows the distribution of vertical stress in the embankment before and after displacing the trapdoor. The stress redistribution results in the prominent increase of the loads acting on piles increasing from 25 kPa to 60 kPa, but the reduction of loads on the soil decreasing from 25 kPa to 9.5 kPa with the trapdoor displacement. The areas of increased vertical stress and corresponding equipotential lines of the contour roughly coincide with the stable regions and external failure surfaces in Figure 8, respectively. It implies that more loads of surrounding soil are transferred to the stable regions through the external failure surface, thus increasing the vertical stress on the pile top.Figure 10
Vertical stress contour: (a) before trapdoor displacement and (b) after trapdoor displacement.
(a)(b)The ratio of the horizontal stress to vertical stress, defined as the coefficient of lateral earth pressure (k), is an important parameter and it is calculated before and after displacing the trapdoor as shown in Figure 11. The coefficient of lateral earth pressure around the external failure surface is 0.48, roughly approaching the coefficient of active earth pressure, ka = tan2(45°-φ/2) = 0.42. On this basis, the coefficient of active earth pressure can be used to represent the coefficient of lateral earth pressure around the external failure surface in the proposed method.Figure 11
Coefficient of lateral earth pressure: (a) before trapdoor displacement and (b) after trapdoor displacement.
(a)(b)


## 2.4.1. Evolution of External Failure Surfaces

The evolution of soil arching can be described by the development of failure surfaces, which are generally explored by the displacement contour [18, 19, 21, 36], and a typical case is illustrated in Figure 7 (s = 320 mm, h = 400 mm). Figure 8(a) demonstrates the development of soil displacement during the displacement of the trapdoor. When the trapdoor reaches the designated displacement of 17.25 mm, the embankment is divided into three regions according to the displacement contour, which includes the active region, transition region, and stable region. According to Costa et al. [21], the boundary between the transition region and active region (curve OA) represents the internal surface and the boundary between the transition region and stable region (curve OB) represents the external surface. The active region will flow downward with the displacement of the trapdoor, making the transition region move obliquely downward to the active region accordingly.Figure 8
Displacement contour of case1A (s = 320 mm, (h = 400 mm): (a) development of soil displacement and (b) zone division of final displacement contour.
(a)(b)In this study, the failure surfaces are outlined in Figure9. To facilitate the description of the failure surface, simplified straight failure surfaces were adopted to approximately represent the actual failure surfaces, and they share the same starting and ending points. The inclination of the simplified failure surface, θ, is the angle between the line of the failure surface and the horizontal plane, as shown in Figure 9. The external failure surfaces are only observed when the embankment height, h, is larger than the pile-net spacing (s-b). For example, for cases with a net spacing of 220 mm, the external failure surface would occur when the embankment heights are 300, 400, 500, and 600 mm. It is also true for the cases with the net spacing of 350 mm when embankment heights are 400, 500, and 600 mm, and the cases with the net spacing of 550 mm when the embankment height is 600 mm. It should be noted that the inclination of the failure surface (θ) decreases gradually with the embankment height.Figure 9
Evolution of the failure surfaces with various embankment height for different pile-net spacing: (a) (s)-(b) = 220 mm, (b) (s)-(b) = 350 mm, and (c) (s)-(b) = 550 mm.
(a)(b)(c)


## 2.4.2. Distribution of Vertical Stress

The contact force between particles can be monitored by measurement circles and then can be converted into continuum stress based on the following equation:(6)σ¯ij=1−n∑Np1Sp∑Np1∑Ncsicsipnic,pFjc,where Np1 is the number of particles whose centroids are contained within the measurement circle; Nc is the number of contacts between particles; n is the porosity within the measurement circle; S(p) is the area of a particle; si(p) and si(c) are the locations of a particle centroid and its contact, respectively; ni(c,p) is the unit normal vector directed from a particle centroid to its contact location; Fj(c) is the force acting on the contact; and indices i and j are in the set of (x, z).The vertical stress,σv, can be calculated based on the following equation for i = z and j = z:(7)σv=σzz,where σzz is the average stress in the direction of the z-axis.Figure10 shows the distribution of vertical stress in the embankment before and after displacing the trapdoor. The stress redistribution results in the prominent increase of the loads acting on piles increasing from 25 kPa to 60 kPa, but the reduction of loads on the soil decreasing from 25 kPa to 9.5 kPa with the trapdoor displacement. The areas of increased vertical stress and corresponding equipotential lines of the contour roughly coincide with the stable regions and external failure surfaces in Figure 8, respectively. It implies that more loads of surrounding soil are transferred to the stable regions through the external failure surface, thus increasing the vertical stress on the pile top.Figure 10
Vertical stress contour: (a) before trapdoor displacement and (b) after trapdoor displacement.
(a)(b)The ratio of the horizontal stress to vertical stress, defined as the coefficient of lateral earth pressure (k), is an important parameter and it is calculated before and after displacing the trapdoor as shown in Figure 11. The coefficient of lateral earth pressure around the external failure surface is 0.48, roughly approaching the coefficient of active earth pressure, ka = tan2(45°-φ/2) = 0.42. On this basis, the coefficient of active earth pressure can be used to represent the coefficient of lateral earth pressure around the external failure surface in the proposed method.Figure 11
Coefficient of lateral earth pressure: (a) before trapdoor displacement and (b) after trapdoor displacement.
(a)(b)


## 3. The Theoretical Model

### 3.1. Theoretical Model Establishment

It is essential to establish a theoretical model for calculating the efficacy of pile-supported embankments with external failure surfaces, and the schematic of the theoretical model is shown in Figure12(a). The coordination system is established with the origin located at the center of the pile top plane, and the right direction and upward direction are denoted as positive, respectively. Considering the elastic-plastic behavior of the soil commonly described by the Mohr-Coulomb model [37], two simplified straight failure surfaces AB and A′B′ are assumed to deviate with the displacement direction (i.e., z-axis) at an angle of α (α = π/4 + φ/2), as shown in Figure 12(a). The pile-supported embankment is a three-dimensional (3D) problem; however, the 2D plane-strain model is much easier to solve and can also obtain reasonable accuracy for engineering practices [38, 39]. In this study, the sliding soil over the trapdoor is sliced into segments along the depth, and the coefficient of active earth pressure is used in the proposed method based on the embankment behavior in the numerical simulation.Figure 12
(a) A schematic diagram of the theoretical model and (b) force equilibrium of segmentd(z).
(a)(b)Considering the vertical force equilibrium of the segment with the thickness of dz, shown in Figure 12(b), the following equation can be derived:(8)PzD+2Fnzsinαdzcosα+2Fnztanϕcosαdzcosα=PszD+PszdD+dPszD+dPszdD+γD+D+dD2,where Ps(z) is the vertical stress at the depth of z, Fn(z) is the stress perpendicular to the failure surface at a depth of z, D is the lower width of the soil unit, α is the angle of the failure surface with a z-axis, φ is the internal friction of soil, and γ is the unit weight of soil.For the horizontal force equilibrium of the differential segment dz, the following equation can be obtained:(9)Fnzcosα−Fnztanϕsinα=Pszka,where ka is the coefficient of active earth pressure, equal to tan2 (45°-φ/2).By rearranging equation (10), the following relationship can be obtained:(10)Fnz=Pszkacosα−tanϕsinα=Pszkaa1,where a1 = cosα-tanφsinα.According to geometrical relationships, we can obtain the following formulas:(11)D=s−b+2ztanα,(12)dD=2tanαdz.Substituting equations (11)–(13) into (9), the following equation stands:(13)2Pszkaa1tanα+tanϕ=2Psztanα+s−b+2ztanαPs′z+γs−b+2ztanα.The boundary condition should satisfy(14)Psb2tanα=γh−b2tanα.Ps (z) can then be solved by combining equations (14) and (15).For the equilibrium of vertical stress atz = 0, we obtain the following equation:(15)psz=0×1−m+ppz=0×m=γh,where m is the area replacement ratio (m = b/s), Ps|z = 0 is the embankment load suffered by the soil, and Pp|z = 0 is the embankment load suffered by the pile, which can be solved by using equation (16).Expression of efficacyE can be obtained from equation (17) and is displayed as(16)E=ppz=0×bppz=0×b+psz=0×s−b.


### 3.2. Comparative Analyses

The parameters used in the comparative analyses are consistent with the experimental tests [31]: the pile net spacings (s-b) are 350 mm and 550 mm, the internal friction of soil φ is 24°, and the weight of soil γ is 62 kN/m3. The calculated results of the efficiency E from the proposed method and the comparison between the experimental results and those from Terzaghi’s method [7] are shown in Figure 12. The method for predicting the efficiency E suggested by Terzaghi [7] can be briefly clarified by the following formulas:Whenh ≤ 2s(17)psz=0=γs−b2ktanϕ1−e−2khs−btanϕppz=0=γhs−psz=0s−bb.Whenh > 2s(18)ps|z=0=γs−b2ktanϕ1−e−4ktanϕ+γh−2s−b,ppz=0=γhs−psz=0s−bb,where k, equals 1, as suggested by Terzaghi.As shown in Figure13, the results calculated by Terzaghi’s method show good agreement with the experimental results when the embankment height, h, is less than the pile-net spacing, s-b. As concluded from the numerical simulation, the embankment failure surface is vertical when the embankment height is less than that of the pile-net spacing and the shape of the failure surface is consistent with the assumption proposed in Terzaghi’s method. Nevertheless, the significant differences from the assumption can be observed when the embankment height is larger than the net spacing because the failure surface no more extends vertically along the edge of the trapdoor but it extends beyond the trapdoor. As a result, the results calculated by the proposed method would be more accurate for high embankments with external failure surfaces compared to those by Terzaghi’s method with the simple vertical failure surface. In conclusion, Terzaghi’s method is appropriate when the embankment height is less than the net spacing, whereas the proposed method is appropriate when the embankment height is larger than the net spacing.Figure 13
Comparison of the results of the proposed method, Terzaghi’s method, and the experiment.


## 3.1. Theoretical Model Establishment

It is essential to establish a theoretical model for calculating the efficacy of pile-supported embankments with external failure surfaces, and the schematic of the theoretical model is shown in Figure12(a). The coordination system is established with the origin located at the center of the pile top plane, and the right direction and upward direction are denoted as positive, respectively. Considering the elastic-plastic behavior of the soil commonly described by the Mohr-Coulomb model [37], two simplified straight failure surfaces AB and A′B′ are assumed to deviate with the displacement direction (i.e., z-axis) at an angle of α (α = π/4 + φ/2), as shown in Figure 12(a). The pile-supported embankment is a three-dimensional (3D) problem; however, the 2D plane-strain model is much easier to solve and can also obtain reasonable accuracy for engineering practices [38, 39]. In this study, the sliding soil over the trapdoor is sliced into segments along the depth, and the coefficient of active earth pressure is used in the proposed method based on the embankment behavior in the numerical simulation.Figure 12
(a) A schematic diagram of the theoretical model and (b) force equilibrium of segmentd(z).
(a)(b)Considering the vertical force equilibrium of the segment with the thickness of dz, shown in Figure 12(b), the following equation can be derived:(8)PzD+2Fnzsinαdzcosα+2Fnztanϕcosαdzcosα=PszD+PszdD+dPszD+dPszdD+γD+D+dD2,where Ps(z) is the vertical stress at the depth of z, Fn(z) is the stress perpendicular to the failure surface at a depth of z, D is the lower width of the soil unit, α is the angle of the failure surface with a z-axis, φ is the internal friction of soil, and γ is the unit weight of soil.For the horizontal force equilibrium of the differential segment dz, the following equation can be obtained:(9)Fnzcosα−Fnztanϕsinα=Pszka,where ka is the coefficient of active earth pressure, equal to tan2 (45°-φ/2).By rearranging equation (10), the following relationship can be obtained:(10)Fnz=Pszkacosα−tanϕsinα=Pszkaa1,where a1 = cosα-tanφsinα.According to geometrical relationships, we can obtain the following formulas:(11)D=s−b+2ztanα,(12)dD=2tanαdz.Substituting equations (11)–(13) into (9), the following equation stands:(13)2Pszkaa1tanα+tanϕ=2Psztanα+s−b+2ztanαPs′z+γs−b+2ztanα.The boundary condition should satisfy(14)Psb2tanα=γh−b2tanα.Ps (z) can then be solved by combining equations (14) and (15).For the equilibrium of vertical stress atz = 0, we obtain the following equation:(15)psz=0×1−m+ppz=0×m=γh,where m is the area replacement ratio (m = b/s), Ps|z = 0 is the embankment load suffered by the soil, and Pp|z = 0 is the embankment load suffered by the pile, which can be solved by using equation (16).Expression of efficacyE can be obtained from equation (17) and is displayed as(16)E=ppz=0×bppz=0×b+psz=0×s−b.


## 3.2. Comparative Analyses

The parameters used in the comparative analyses are consistent with the experimental tests [31]: the pile net spacings (s-b) are 350 mm and 550 mm, the internal friction of soil φ is 24°, and the weight of soil γ is 62 kN/m3. The calculated results of the efficiency E from the proposed method and the comparison between the experimental results and those from Terzaghi’s method [7] are shown in Figure 12. The method for predicting the efficiency E suggested by Terzaghi [7] can be briefly clarified by the following formulas:Whenh ≤ 2s(17)psz=0=γs−b2ktanϕ1−e−2khs−btanϕppz=0=γhs−psz=0s−bb.Whenh > 2s(18)ps|z=0=γs−b2ktanϕ1−e−4ktanϕ+γh−2s−b,ppz=0=γhs−psz=0s−bb,where k, equals 1, as suggested by Terzaghi.As shown in Figure13, the results calculated by Terzaghi’s method show good agreement with the experimental results when the embankment height, h, is less than the pile-net spacing, s-b. As concluded from the numerical simulation, the embankment failure surface is vertical when the embankment height is less than that of the pile-net spacing and the shape of the failure surface is consistent with the assumption proposed in Terzaghi’s method. Nevertheless, the significant differences from the assumption can be observed when the embankment height is larger than the net spacing because the failure surface no more extends vertically along the edge of the trapdoor but it extends beyond the trapdoor. As a result, the results calculated by the proposed method would be more accurate for high embankments with external failure surfaces compared to those by Terzaghi’s method with the simple vertical failure surface. In conclusion, Terzaghi’s method is appropriate when the embankment height is less than the net spacing, whereas the proposed method is appropriate when the embankment height is larger than the net spacing.Figure 13
Comparison of the results of the proposed method, Terzaghi’s method, and the experiment.


## 4. Conclusion

The evolution of soil arching in pile-supported embankments plays an important role in the selection of arching models, and it evolves with embankment height, pile diameter, and pile spacing. In addition to internal failure surfaces, there are also external failure surfaces in some cases of the pile-supported embankment. Numerical models were developed in this study to analyze the external failure surfaces. To accurately replicate the initial state of the embankment, an improved multilayer compaction method was adopted so that the porosity of the generated particles in the PFC approaches the realistic embankment fill. The micromechanical parameters of embankment fills were calibrated by numerical biaxial tests at various confining pressures, and the numerical model was then validated by comparing the efficacy with the case observation. It is concluded that the external failure surface will only occur when the embankment height is greater than the pile-net spacing and the embankment load is redistributed towards the pile head through the external failure surface. The coefficient of lateral earth pressure around the external failure surface is 0.48, and it is close to the coefficient of active earth pressure. The calculation results indicate that the proposed solution demonstrates better consistency with experimental results than that of Terzaghi’s method for the case with external failure surfaces.



---

*Source: 7218904-2022-09-13.xml*