Abstract
This paper presents a new numerical approach for computing the internal force and displacement of portal double-row piles used to stabilize potential landslide. First, the new differential equations governing the mechanical behaviour of the stabilizing pile are formulated and the boundary conditions are mathematically specified. Then, the problem is numerically solved by the high-accuracy Runge-Kutta finite difference method. A program package has been developed in MATLAB depending on the proposed algorithm. Illustrative examples are presented to demonstrate the validity of the developed program. In short, the proposed approach is a practical new idea for analyzing the portal double-row stabilizing pile as a useful supplement to traditional methods such as FEM.
1. Introduction
Double-row stabilizing piles have been widely used in slope reinforcement engineering and treatment of landslide geological disasters, which have some advantages such as larger rigidity, less displacement in the top of the piles, and large resisting force. Existing methods for the analysis of double-row stabilizing piles can be generally classified into the following two categories [1–6]: (1) coupled method (continuum analysis) that simultaneously solves pile response and slope stability [7]; (3) uncoupled method which deals with pile and slope separately. In the uncoupled method, pile-soil interaction is commonly represented by equivalent Winkler or p-y springs [8–13].
As for coupled method, the finite element method is certainly the most comprehensive approach to study pile-slope stability. However, its use generally requires high numerical costs and accurate measurements of material properties. This makes the use of this method rather unattractive for practical applications [6].
To date, in practical engineering applications, the uncoupled method is the most widely used approach to design the double-row reinforcing piles to increase slope stability due to its simplicity of use. First, the lateral force acting on the pile segment above the slip surface due to soil movement is evaluated usually by the limit equilibrium method. Second, the response of the double-row pile subjected to lateral loading is analysed by FEM modeling it as a beam resting on linear or nonlinear soil/rock spring supports. The FEM modeling is reasonably accurate but complicated and time consuming.
In this paper, a new uncoupled method to compute the response of portal double-row piles subjected to lateral earth pressure loading based on new boundary value problem approach is introduced. First, the new governing differential equations including six variables (three internal forces and three displacements) are formulated and the boundary condition is specified. Second, the high-accuracy Runge-Kutta differential method is used to solve the corresponding system of differential equations to obtain the pile’s internal forces and displacements. A program for pile response analysis and graphics edit is developed. At last, the program was verified against the FEM analysis results in terms of pile deflection, bending moment, and shear force along the length of the pile.
The objective of this study is to provide an alternative method for the design of portal double-row pile used for slope stabilization or earth retaining.
2. Derivation of Governing Differential Equations
In order to solve the complicated engineering problem of the response of double-row stabilizing pile under laterally loading by using accurately mathematical model, we often need to define its boundary value problem which involves the governing differential equations and corresponding boundary conditions. Then, closed-form or numerical solutions for the engineering problem can be obtained by many appropriate mathematical methods.
Under the scheme of uncoupled analysis of the pile (as shown in Figure 1), the new governing differential equations for stabilizing piles embedded in slope will be developed in the following according to the general principles of solid and structural mechanics including static force equilibrium, deformation compatibility, and constitutive relationship.
2.1. The Loading Condition
Below the slip surface, the sliding bed is assumed to be stable and cannot move. Before the active force induced by the sliding mass act upon the pile segment above the slip surface, the earth pressure acting at the front and back sides of the pile segment below the slip surface is in equilibrium. So it can be neglected. Only the active force induced by the movement of the sliding mass is considered in the calculation model, which are loaded on the pile segment above the slip surface. Their distribution along the pile shaft can be assumed as uniform, triangular, trapezoidal, and rectangular profiles. Then, the reaction force acting at the pile segment below the slip surface is calculated. These active forces and reaction forces considered are in equilibrium [14].
2.2. The Soil-Pile Interaction Model
Due to its simplicity and reasonable accuracy, the Winkler foundation model is adopted in current analysis to describe the pile-soil interaction behavior. The Winkler method assumed that the substratum is composed of independent horizontal springs. Under the Winkler hypothesis, the soil reaction pressures () acting on the pile can be modeled by discrete independent linear or nonlinear springs in the form of the following equation: where is the spring constant, also called the modulus of horizontal subgrade reaction (it has a unit of force/length3). The main difference between the different Winkler foundation models available is in the selection of the foundation stiffness coefficients. is the horizontal soil reaction pressure (it has a unit of force/length2). is the horizontal displacement (it has a unit of length).
2.3. The New Equilibrium Differential Equations
Let us consider an isolated free portion of pile, as shown in Figure 2, having a infinitesimal length of and acted upon by external distributed normal load and tangential load . The free segment can be imagined to be cut out of the pile, and the internal forces () in the original pile may become external forces on the isolated free portion.
Figure 2 also shows the coordinate system used in this paper for introducing the governing equations.
We adopt sign conventions so that the six variables as shown in Figure 2 are positive. The sign convention adopted for forces is that positive sign indicates tensile axial force , the positive shearing force should be directed so that they will tend to rotate the element counterclockwise, and the positive bending moment will tend to make the element concave leftward. The sign convention adopted for displacements is that the positive normal displacement points outward normal, the positive tangential displacement points right when facing outward normal, and the positive is counterclockwise. The lateral pressure is considered positive when applied from left to right. The lateral pressure is considered positive when applied from up to down.
Thus, considering the equilibrium of the above infinitesimal pile segment, as it bends under the action of the applied loads (shown in Figure 2), we arrive at two force equilibrium equations in the directions of and and one moment equilibrium equation:
Simplifying the above equations and neglecting the higher-order term and the term with the square of the differential, the equilibrium equations now take the following forms: where the “one” indicates that Winkler reaction force exists and the “zero” indicates that Winkler reaction force does not exist. is the position coordinate; is perimeter of the cross-section. is the width of the cross-section, is the Winkler modulus of vertical subgrade reaction, and is the Winkler modulus of horizontal subgrade reaction.
For the sake of convenience of formula deducing, let
The set of equations of equilibrium (3) can be rewritten in the following matrix form: where ; .
2.4. Geometric and Constitutive Equations
When the deformation (, , ) of the differential element (shown in Figure 2) induced by the internal forces () is considered given, the corresponding strains can be expressed as According to the related theory of elastic beam, the internal forces () can be related to strains as in the following linear constitutive equation: where is the constant related to the shape of pile cross-section ( for rectangular cross-section; for circular cross-section); is the cross-sectional area; is the flexural rigidity of the pile’s cross-section.
Considering the deformation, can be decomposed into two parts. One part is the induced by the displacement on its direction and another part is the projection of other displacement onto this direction which takes the form , where is the a undetermined third-order square matrix. Then, the deformation can be expressed as Applying the principle of virtual work to the isolated differential element of pile (shown in Figure 2), can be determined. We suppose that each point of the body is given an infinitesimal virtual displacement satisfying displacement boundary conditions where prescribed. The virtual deformation associated with the infinitesimal virtual displacement is . The virtual work of the external surface forces is , where . The virtual work of the internal forces is . By equating the external work to the internal work, we have . Substituting (5) and (8) into the above equation and simplifying yields . Since this equation is satisfied for arbitrary , the terms in the brackets in the integral must vanish at every point which means that . At last, we develop the following geometric and constitutive equations:
2.5. Matrix Form of the Governing Equations
For the sake of convenience of problem solving, combining the three equilibrium differential equations (3) and the three geometric and constitutive equations (9) leads to a system of six equations:
Then, we can use matrix notation to present the above equations (10) in the form: where is the coefficients square matrix of six order and and represent two column matrices.
This system of six independent differential equations can be solved for six unknown functions (three independent forces and three independent displacements).
3. Boundary Conditions
As mentioned above, there are total six unknowns to be determined (). Therefore, six boundary conditions are needed for the problem solving.
The boundary conditions for (12) are determined according to the way in which the pile’s head and base are supported or restrained. There are three conditions at the base point and three conditions at the head point. We use matrix notation to present these boundary conditions in the following form: where indicates the beginning point of calculation (the pile’s base point), indicates the end point of calculation (the pile’s head point), is the matrix of boundary condition on the beginning point, and is the matrix of boundary condition on the end point. They are matrices. In this study, the following possible pile end conditions were considered.(1)Free head (allows both displacement and rotation): , . The corresponding matrix of boundary condition is (2)In case of bottom end hinged (allows rotation without displacement): , . The corresponding matrix of boundary condition is (3)In case of bottom end fixed (allows neither displacement nor rotation): , . The corresponding matrix of boundary condition is (4)In case of bottom end partially hinged (allows rotation without vertical displacement): , , and . The corresponding matrix of boundary condition is (5)Elastic vertical support at the bottom end: , , and . The corresponding matrix of boundary condition is where is the modulus of vertical compressibility. Area is the cross-sectional area of the pile.
4. Definition of Boundary Value Problem
Next, we impose the boundary conditions (13) at the pile head and base upon the derived new governing differential equations (12) to define a boundary value problem of the following equations: To this end, the response of double-row portal stabilizing pile is mathematically idealised as the boundary value problem of (19).
Thus, many numerical methods to solve the ordinary differential equations can be adopted to solve the boundary value problem of (19).
It should be noted that the existence and uniqueness of solution for the boundary value problem of (19) should be mathematically proved. This matter is however outside the scope of the writer’s major. According to the physical character of the problem, we can imagine that the solution exists and is unique. The solution can be validated through comparative studies.
5. Method of Solution
5.1. Uniformity Preprocessing
The orders of magnitude of the section internal forces () are so much higher than those of the displacements () that numerical solving of the equations may meet singularity difficulty. So for reasons of numerical stability, it is necessary to perform uniformity preprocessing for the order of magnitude of the element in the coefficient matrix . We multiply the displacements () by and substitute the original displacement variables by the expressions (). So, we redefine two variables as follows:
Finally, we obtain the following system of ordinary differential equations:
This system of six independent differential equations subjected to boundary conditions can be numerically solved for six unknown functions (three forces and three displacements).
5.2. The Runge-Kutta Finite Difference Algorithm
The Runge-Kutta algorithm is commonly used for the solution of the ordinary differential equation of the form . So, it is chosen to solve (21).
5.2.1. Derivation of the Recursion Formula
The following finite difference formula (22) is one format of the Runge-Kutta methods:
The Runge-Kutta algorithm of this type (22) is a numerical method of fifth order, where is the stepsize of difference and . For convenience of computation, this formulation may be rewritten in the following form:
As we can see, the value of function at point () can be determined from the value of function at point (), where represents the beginning point of calculation and represents the end point of calculation.
In the above equation, and can be obtained from the following recursion formula: where is the identity matrix.
5.2.2. Determination of the Initial Vector
The initial value is the start point of the recursion formula. Now, we discuss in the following how to obtain the initial vector by using the recursion formula of (23) and imposing the boundary conditions at pile head and base.
Considering the recursion formula of (23), can be expressed in terms of as follows:
In the case of , we have and substitute it into the recursion formula of (23). We get
It can be rewritten as follows:
Comparing the above two equations, the recursion formula for and is obtained as follows:
Now considering the case of boundary point: , we substitute the boundary conditions at end point , into the above equation. This leads to the equation to solve for :
The above set of linear algebraic equations can be solved for by using the method of Gaussian elimination with pivot selection. Once is known, can be obtained in sequence using the recursion formula of (23).
5.2.3. Coordinate Transformation between Pile and Beam
When solving the problem of double-row portal piles using the recursion formula of (23), due to the direction change of the axes of the pile and the beam at the connection point (shown in Figure 3), first we need to distinguish and in the recursion formula for the two connected segments. And then in order to satisfy the equilibrium of internal forces and maintain the continuity of displacements at the connection point, we need to introduce the so-called connection matrix to the corresponding formula when dealing with , , and , recursively.
The centroidal axes’ rotation from pile to beam at connection point means mathematically that the local coordinate system rotates clockwise by degree at the connection point (shown in Figure 3).
According to the principle of equilibrium of internal forces, we get And the bending moment remains unchanged.
According to the principle of vector analysis, we get And the rotation displacement remains unchanged.
Combining the above equations, the corresponding transformation matrix for the local coordinate system rotating clockwise by degree can be expressed as follows:
As we can see, the above connection matrix is the so-called orthogonal matrix whose inverse is its transposed matrix.
As shown in Figure 3, the node of difference (no. ) is the connection node. First we need to distinguish and in the recursion formula for the two connected segments and then insert the connection matrix for transformation from to when dealing with the , , and , recursively. The detailed procedure is as follows.(1)For the calculation of , , where is the transformation matrix for point . Thus, . Next, where indicates the end point of calculation.(2)For the calculation of : because , and So, the above procedure also applies to calculation of as follows: (3)For the calculation of , where denotes the end point of calculation.
5.2.4. The Solution Flow Process
In short, the proposed solution procedure involves four main steps:(1)calculating the value of and using the given Equation (25);(2)calculating and using the given recursion formula of (29);(3)calculating the vector by solving linear algebraic Equations (30);(4)calculating using the given recursion formula of (23).
Because the equations and solution formula are all given in form of matrices, a simple computer program has been written on the platform of MATLAB to run this procedure. At last, we can get the shear, bending-moment, and deflection diagram along the pile.
6. Verification
The practical examples of portal double-row piles used to stabilize an potential landslide (shown in Figure 4) are considered herein to verify the developed numerical calculation techniques. Soil strength parameters used in the stability analysis are from laboratory shear testing on the undisturbed soil samples. The resisting (shear) force required to achieve the desired safety factor and transferred by the pile is estimated to be 2147 kN/m. Before the pile is installed, the slope is approaching limit state and the safety factor can be assumed to be one. The values of and can be determined based on experiment data and satisfying this limit state condition. Then, the required resisting force 2147 kN/m can be obtained using back analysis. When this additional force is applied to the specified place of the landslide, the safety factor calculated using the limit equilibrium method can achieve the desired value 1.3. The manual digging discrete reinforced concrete piles were designed to be installed at a spacing of 4 m to increase the factor of safety of the whole slope to the required value of 1.3.
The front row pile has a length m and sectional dimension . The length of the portion embedded into the sliding surface is 6.0 m. The back row pile has a length m and sectional dimension . The length of the portion embedded into the sliding surface is 18.0 m. The sectional dimension of the connection beam is 0.8 0.8 m. The pile is constructed using C30 concrete (assuming that concrete does not crack during working). A lateral force is assumed to act upon the pile segment above the sliding surface. The pressure distribution is considered to have a rectangular shape, as proposed by the Chinese design code (Code for design on retaining structures of railway subgrade no. TB10025-2006).
The conceptual calculation model used to simulate the lateral response of the pile is shown in Figure 5. The boundary condition at the pile base is considered as free head which allows both lateral displacement and rotation.
For simplicity of program editing, the friction force along the pile-soil interface can be neglected. The modulus of horizontal subgrade reaction for the mudstone below the slip surface is listed in Table 1 for the three different mudstone layers. The coefficient of subgrade reaction was determined according to the suggestions for the mudstone in the related Chinese design code and in situ tests. The related Chinese design code provides detailed experiment procedure to determine the modulus of horizontal subgrade reaction.
Then, the developed method is applied to analyze the lateral response of the double-row pile which is also computed by the FEM program we developed. We employ an elastic beam column element to model the pile and horizontal spring element to represent the reactions of the surrounding soil in the FEM model. Comparisons of shear force, bending moment, and deflection of the pile between boundary value method (BVM) and FEM are presented in Figures 6 and 7. Complete agreement between them can be observed.
Through the above comparative studies, it has been found that the program we developed works very well and can replace the existing numerical methods that have been used to design the portal double-row stabilizing pile.
7. Summary and Conclusions
In this paper, a new numerical uncoupled method for calculating the response of portal double-row stabilizing piles is proposed. The theoretical background and a detailed derivation of the proposed numerical solution scheme are described. The feasibility of the method developed was verified using the comparative case study. The proposed method has more higher modeling and computing efficiency than the FEM and can be an alternative method for analyzing the behavior of portal double-row piles used for slope stabilization.
Acknowledgments
This study was jointly financially supported by grants from the National Natural Science Foundation of China (Grant no. 51008298) and Key Research Program of Chinese Academy of Sciences (no. KZZD-EW-05). The support is gratefully acknowledged. The author also thanks the anonymous referees whose comments helped to improve the work and the presentation of this paper.