Numerical method for slope stabilization design with piles

In this paper, a detailed numerical approach is introduced for the slope stabilization design with piles considering the pile capability as well as the slope stability. Two simulation methods for stabilizing piles, i.e., the solid element simulation and the composite element simulation, are introduced and discussed about their applicability. The solid element simulation is found to be potential due to good accuracy and computation efficiency. Interface is employed for pile-soil interaction. Together with strength reduction method, the numerical techniques are applied in the reinforcement for a three-dimensional slope with one row of piles. By varying pile length, pile spacing and pile location, the optimal arrangement of pile reinforcement for the slope is discussed considering both slope stability and pile mechanical behavior. Therefore, an economical as well as effective pile stabilization design is able to be achieved through the suggested practical numerical method in this paper.


Introduction
Stabilizing piles are widely used in both natural landslide treatment and excavated slope reinforcement due to their permanent stabilization and easily construction. The design of pile-reinforced slopes involves two major issues, i.e., the effective improvement of the slope stability and the appropriate arrangement of the piles considering their bear capacity. The three-dimensional (3D) pile-soil interaction has long been recognized as very complicated [1] . Numerous design methods are proposed for slope stabilization design with piles for decades. However, due to the limitation of computation capability in the past, decoupled methods [2][3][4] are mainly proposed that the piles are designed to provide the desired resistance for the slope stability enhancement. In such cases, either the lateral soil pressure or the lateral soil displacement is utilized for the evaluation of the resistance force provided by the proper pile arrangement. The procedures are generally complicated and may require assumptions, for instance, only the soil around the piles is assumed to be in plastic equilibrium state in Ito and Matsui [5] to estimate the lateral soil pressure and a uniformly distributed soil displacement is adopted in Kourkoulis et al [4] to evaluate the lateral load on the piles.
Owing to the rapid development of user-friendly software and high-performance processors, it is able to perform 3D numerical simulations for pile-reinforced slopes at a satisfactory efficiency. Therefore, the coupled method that can simultaneously evaluate the slope stability and pile mechanical behavior using numerical techniques is quite potential and desirable for the further design of stabilizing piles. In this paper, two pile simulation techniques are introduced and discussed for efficiently and accurately analyzing pile mechanical behavior. Interface is adopted for pile-soil interaction. Together with strength reduction method for slope stability evaluation, the numerical techniques are successfully applied in an illustrative example of the stabilizing pile design. Figure 1 displayed the geometry of a pile-reinforced slope, where Lx is the horizontal distance from the pile install location to the slope toe, Lp is the pile length and S is the pile center spacing. The slope is of 12 m in height with a gradient of 1 V:1.5 H. 3D numerical simulation techniques are introduced in the following to estimate pile mechanical behavior as well as the slope stability.

Two simulation methods for stabilizing piles
Concrete piles with a square section are used for slope stabilization in this paper as shown in Figure 1, where W is the square section width. The piles are treated as the linear elastic material. To consider the 3D geometry effects of piles as well as obtaining the internal force, two simulation techniques are presented in the following.

Solid element simulation.
It is straightforward to simulate the 3D geometry effects of piles with solid elements. While the deflection of the solid-element pile model can be captured through the displacement of the grid points, the challenge lies in obtaining the pile internal force. In practice, the bending moments in a pile can be calculated by the integral of the normal stress multiplying the distance to the neutral axis over the cross-sectional area, because the bending moment is the resultant of the forces due to the normal stress acting on the section [6] . The neutral axis is usually regarded to be perpendicular to the bending direction and through the centroid of the cross section. For a cross section in Figure 2, where the bending direction is along X-axis, the bending moment M can be calculated as follows where σzz is the normal stress acting on the cross section, dA is an element of the cross-sectional area and δ is the distance from the centroid of the element to the neutral axis of the cross section. In this paper, the numerical simulation is carried out in FLAC 3D , where σzz can be conveniently derived by the FISH function zone.field.get with the extrapolation method being "polynomial" and be employed for the bending moment calculation through user-defined FISH codes.

Composite element simulation.
To avoid the programming for bending moment calculation in the solid element simulation approach, a composite simulation approach using the couple of solid and structural elements can be adopted for an easy access of pile internal force. Solid elements are still employed for 3D effects, however, the pile stiffness is introduced by beam elements at the center of the pile, i.e., the pile stiffness is assigned to the central beam elements while nearly zero stiffness is assigned to the solid elements. Meanwhile, each solid element grid point around the central beam 3 elements is connected to the central beam element node at the same height through another stiff beam element as shown in Figure 3. The aim is to allow the external forces applied on the solid element grid points to be fully delivered to the central beam elements. To make sure the beam elements and the solid elements deform compatibly, the rigid connection should be set between the solid element grid points and the beam element nodes at the same position by specifying link properties in FLAC 3D . In this way, the composite element pile model is established, where the 3D effects are considered due to the presence of the solid elements and the beam elements acting like the skeleton to bear load so that the pile mechanical behavior is captured directly via the central beam elements.

Validation and discussion of pile simulation methods
As shown in Figure 4, two types of beams with different boundary conditions, i.e., the cantilever beam and the clamped beam, are presented for the validation of the two suggested simulation methods. The theoretical solutions of bending moment M and deflection f for the two types of beams according to Euler-Bernoulli beam theory can be written as follows where q is a uniform load, L is length of the beams, h is the distance of from the bottom of the beams, ξ = (L-h)/ L is a parameter and EI is the bending stiffness of the beams. Equations (2.1) and (2.2) are for the cantilever beam, while equations (3.1) and (3.2) are for the clamped beam. In the following, the two types of beams also with a square section are analyzed using the two simulation methods as well as Euler-Bernoulli beam theory under the conditions that q = 100 kN/m, L = 10 m, E = 30 GPa and W = 0.5 m. Under the above conditions, the deflection and bending moment in the cantilever beam using the two suggested simulation methods are of good agreement with those based on Euler-Bernoulli beam theory according to Figure 5(a) and Figure 5(b). And the same is for the clamped beam in Figure 5(c) and Figure 5(d). To determine applicability of the two simulation methods, the influence of beam size and mesh size are further discussed in the following.

Influence of beam size.
Euler-Bernoulli beam theory works well for thin beams with length to thickness ratio of the order 20 or more, to which both the two simulation methods are applicable according to the analysis above. In the following, thick beams with W = 2 m are studied keeping the other conditions the same, and the results are exhibited in Figure 6. For the cantilever beam, the bending moments calculated via both two simulation methods are coincident with that based on Euler-Bernoulli beam theory as can be observed from Figure 6 (a) and Figure 6 (b), while the deflection based on solid element simulation is slightly larger than that based on Euler-Bernoulli beam theory and the deflection using composite element simulation agrees with that based on Euler-Bernoulli beam theory. And the results for the clamped beam is similar, however, the deflection based on solid element simulation is significantly larger than those based on Euler-Bernoulli beam theory and composite element simulation according to Figure 6 (c) and Figure 6 (d).
From the above analysis, it proves the beam elements in FLAC 3D are based on Euler-Bernoulli beam theory. Since the deformation obtained from solid element simulation is supposed to be close to the accurate solution as long as the mesh is fine enough, the different deflection results indicate the

2.2.2.
Influence of mesh size. The above analyses are conducted using relatively small mesh sizes, i.e., the solid element pile model consists of cubic elements with an edge length of 0.0625 m and the size along the pile length for the composite element pile model is also 0.0625 m with the sizes for other two directions being W/2 since the external force is directly delivered from the solid element grid points through the stiff beam elements to the central beam elements. In the following, the clamped beam model with 2 times and 4 times the original mesh size is analyzed using the two simulation methods. The thin beam conditions above are considered and the results are displayed in Figure 7. Apart from the calculation accuracy, the calculation efficiency should also be considered. The analyses in this paper are performed on a computer with an AMD Ryzen Threadripper 3970X 32-Core processor clocked at 3.69 GHz. Table 1 lists the computation time for each simulation. The computation time decreases with the increase of mesh size for both simulation methods as can be observed from Table 1, however, it is much more time consuming using composite element simulation than using solid element simulation with the equal mesh size. It indicates that the solid element simulation can still be a good choice even fine mesh is required.

Soil-pile interaction simulation
Interface elements are employed for soil-pile interaction simulation. Figure 8 gives the diagram for the interface constitutive model that is defined by a slider, a tensile strength, a normal stiffness kn and a shear stiffness ks. The interface normal behavior is defined by kn and the tensile strength that is usually set to take no effective tension; the shear behavior of interface elements is limited by ks and the slider that limits the maximum shear force fsmax for an interface node through the Coulomb shear-strength criterion as follows [ where cs and φs are the cohesion and the friction angle for the interface elements, respectively; a is the representative area associated with the interface node and fn is the normal force acting on the node. Once the criterion is satisfied, sliding begins. kn and ks are recommended to be 10 times the equivalent stiffness of the stiffest neighboring zone [7] that are given as follows where K and G are the bulk modulus and the shear modulus, respectively; ΔZmin is the smallest dimension of an adjoining zone in the normal direction.

Strength reduction method
Strength reduction method (SRM) is a powerful tool for slope stability analysis based on 3D numerical model, by which the factor of safety (FOS) of the slope can be estimated and the critical slip surface is automatically identified as well [8] . where cmin and φmin are the values for c and φ at limit state, respectively. For a linear elastic-perfectly plastic constitutive model, e.g., Mohr-Coulomb yield condition, the bracketing and bisection method can be adopted to implement SRM for the FOS because a sudden transition from elastic to plastic behavior of the material can be captured [9] . In the following, we will illustrate the application of the above simulation techniques for the design of slope stabilization with piles.

Illustrative example
For the pile-reinforced slope presented in Figure 1, one row of piles with W = 1.5 m and E = 30 GPa are used. The soil is modeled with the Mohr-Coulomb yield condition considering full tension cut-off and the non-associated rule. The cohesion, friction angle and Young's modulus of the soil are c = 10 kPa, φ = 20 o and Es = 40 MPa, respectively. The semi-symmetric finite difference model is used for computation efficiency as presented in Figure 9. In the following, the design of pile reinforcement is illustrated by considering the pile bending moment distribution as well as the FOS.

Application of the two pile simulation techniques
Consider the pile reinforcement conditions with Lp = 15 m, Lx = 9 m and S/W = 4, the applicability of the two pile simulation techniques is further determined in following. The critical slip surfaces determined using SRM can be clearly defined by the maximum shear strain in the slope [7,10] . Figure 10 exhibits the critical slip surfaces based on the maximum shear strain as well as the corresponding FOS using the two simulation techniques. Figure 11 presents the pile bending moment distribution using the two simulation techniques.  According to Figure 10, the critical slip surfaces based on the solid element simulation and the solid element simulation are quite close as well as the corresponding FOS. Figure 11 presents the close results of pile bending moment distribution using the two pile simulation methods. Nevertheless, the computation efficiency differs much that it takes 43.29 min using the solid element simulation while it costs 455.57 min using the composite simulation. In such a case, the solid element simulation technique is adopted for the further parametric analysis.

Pile reinforcement parametric analysis
The design of pile reinforcement is supposed to be economical as well as effective. In the following, the pile length, pile spacing and pile location are discussed for the optimal pile configuration. Figure 12 depicts the critical slip surfaces developing at the pile center with Lp = 12 m and Lp = 18 m when Lx = 9 m and S/W = 4 as well as the corresponding FOS. When Lp = 12 m, the pile reinforcement is relatively insufficient because a deep critical slip surface nearly goes through the pile bottom and the corresponding FOS is quite smaller than that when Lp = 15 m although the pile bending moments stay low as shown in Figure 13. When Lp = 18 m, the critical slip surface is nearly the same with that when Lp = 15 m that is shallow and divided into two parts due to the pile reinforcement; the corresponding FOS increases slowly when Lp increases from 15 m to 18 m but the pile bending moments grow quite a bit when it goes deeper than 6 m along the pile length as presented in Figure 13. This indicates Lp larger than 15 may be unnecessary that may cause additional structural requirements and construction work.  Figure 12. Critical slip surface of the slope with using different pile length When pile spacing is relatively small, i.e., S/W = 2, the critical slip surface mainly develops locally above the pile. When S/W increases to 4, the upper and lower parts of the critical slip surface almost connected between the piles. And the critical failure completely goes through the toe to the top between the piles when S/W = 6. The corresponding FOS decreases slowly when S/W increases from 2 to 4 and drops significantly when S/W > 4. Figure 15 displays the pile bending moment distribution with different pile spacing. The bending moments grow obviously when S/W increases from 2 to 4 and increase slowly when S/W increases from 4 to 6. Overall, S/W ≤ 4 is recommended for the slope for an effective reinforcement but a proper pile spacing may be further determined for the balance of the structural requirements for enough bearing capacity and the involved construction work.   Figure 16 presents the critical slip surfaces developing at the pile center with Lx = 4.5 m and Lx = 13.5 m when Lp = 15 and S/W = 4 as well as the corresponding FOS. When Lx = 4.5 m, the critical slip surface develops above the pile; when Lx = 13.5 m, the critical slip surface develops in front of the pile. In both cases, the FOS are quite small indicating inadequate pile reinforcement. Figure 17 shows that relatively small bending moments develops in both cases that can hardly make full use of the pile bear capacity. Therefore, the row of piles are preferred to be installed at around the middle of the slope.

Conclusions
This paper presents a detailed approach for the slope stabilization design with piles using numerical simulations. Two simulation techniques for stabilizing piles are introduced and discussed. The solid element simulation method is recommended for good accuracy and computation efficiency. Interface is employed for pile-soil interaction. Together with strength reduction method, the numerical techniques are applied in the reinforcement for a three-dimensional slope with one row of piles. The illustrative example proves the suggested approach can provide a useful and practical way for the design of pile-reinforced slopes considering not only the slope stability but also the pile capability.