# New Mathematical Modelling of Stabilizing Pile with Prestressed Tieback Anchors

Authors: Cheng Huang; Wei-zhong Ren; Ling-wei Kong
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

This paper presents a novel mathematical modelling for analyzing stabilizing piles with prestressed tieback anchors. The new differential equations governing the mechanical response of the stabilizing pile are formulated and the boundary conditions considering the tie-back anchors are mathematically specified. Then, the system of differential equations is numerically solved by the high-accuracy Runge-Kutta finite difference method. A simple computer program has been written on the platform of MATLAB to run the procedure of the proposed algorithm. This approach is entirely different from the traditional finite element method used to design the anchored piles. The FEM is employed to verify the feasibility of the developed method. The comparative case study indicates that the proposed method has more higher modeling and computing efficiency than the FEM and can be an alternative method for designing the anchored pile used for slope stabilization.


---

## Body

## 1. Introduction

For many years, the stabilizing piles built in China to stabilize existing or potential landslides were mostly cantilever. Later, prestressed ground anchors were added to reduce the number and size of the piles. It is shown that anchor prestressing plays an important role in limiting maximum bending moment on the piles and controlling deflection of the pile head. So anchored stabilizing piles offer considerable economic advantages over the conventional cantilevered piles [1–4]. Combined pile tieback anchor systems will find more extensive use as slope stabilization systems in the future.Existing methods for the analysis of stabilizing piles with prestressed tie-back anchor can be generally classified into the following two categories [5–9]: (a) coupled methods (also called continuum analysis) that simultaneously solve pile response and slope stability [10–13] and (b) uncoupled methods which deal with pile and slope separately. In the uncoupled method, the pile response is analyzed using the beam on elastic foundation approach or FEM and pile-soil interaction is represented by equivalent Winkler or p-y springs [14–20].The finite element method is certainly the most comprehensive coupled method to study pile-slope stability. However, its use generally requires high numerical costs and accurate measurements of material properties. This makes the use of coupled approach rather unattractive for practical applications.To date, in practical engineering applications, the uncoupled method is the most widely used approach to design the 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 anchored pile subjected to lateral loading is analysed by FEM modeling it as a beam resting on linear or nonlinear soil/rock spring supports. And we know the conventional subgrade reaction solution available in literature can not take into account the tie-back anchors [15, 18, 19].In this paper, a new uncoupled method to compute the response of anchored piles subjected to lateral earth pressure loading based on new boundary-value problem approach is introduced. This approach is entirely different from the traditional finite element method using the beam elements and spring elements to model the anchored piles. In the following, first, the new governing differential equations including six variables (three internal forces and three displacements) are formulated and the boundary condition considering the tie-back anchor is given. 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, displacements, and the anchor’s internal force. A program for anchored 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 aim of this work is to present an alternative approach based on new governing differential equations to analyze the response of anchored piles used for slope stabilization or earth retaining. The efficiency and accuracy of this developed method are demonstrated through comparative case study.


## 2. Differential Equations Governing the Response of Stabilizing Pile

As known, in order to solve the complicated engineering problem of the response of stabilizing pile with prestressed tie-back anchors by using an accurate mathematical method, it is often needed 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.


### 2.1. Conventional Governing Differential Equation for Pile Deflection

Before deducing the new differential equations governing the arching mechanism of stabilizing piles, we need to review the conventional governing differential equations for pile deflections which can be found in many of the available literature [1, 15]. It is reported here only for the sake of completeness.We assume a planar cartesian coordinate systemw-z, as shown in Figure 1, with its origin at the center of the pile head such that the z-axis coincides with the pile axis and the w-z plane coincides with the plane of the paper. A lateral force Q and a moment M are assumed to be applied at the pile head.Figure 1
The coordinate system of laterally loaded pile model.The pile-soil interaction response can be idealised as a vertical beam placed in a Winkler spring medium. The pile material is assumed to follow linear elastic behaviour.We know that vertical piles resist lateral loads or moments by deflecting until the necessary reaction in the surrounding soil is mobilized.Neglecting the friction force along the pile-soil interface, the deflection response of piles subjected to lateral load and constrained by surrounding Winkler springs is governed by the following differential equation:(1)EId4w(z)dz4+k(z)w(z)b=0,
where z is the position coordinate, w(z) is the deflection function of the pile (it has a unit of length), E is the pile’s Young’s modulus, I is the second moment of inertia of the pile’s cross-section (EI is called the flexural rigidity of the pile), k is the spring constant, also called the modulus of horizontal subgrade reaction (it has a unit of force/length3), which can be assumed to vary linearly with depth for soil or to remain unchanged for rock. Much of the available engineering experience and assumptions can be used to determine its value, and b is the section width.Based on (1), the deflection function w(z) can be obtained by power series solution, finite difference method, and FEM in terms of bar and beam. The detailed solution scheme can be found in many of the available literature. Once the deflection function w(z) is established, the bending moment, shear force, and soil reaction force at various depths along the pile can be deduced by differentiating the deflection function as follows:
(2)φ=dwdz,M=-EId2wdz2,Q=-EId3wdz3.


### 2.2. Derivation of the New Governing Equations

Under the scheme of uncoupled analysis of the pile (as shown in Figure2), 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.Figure 2
The uncoupled analysis model for the stabilizing pile with prestressed tieback anchor.


#### 2.2.1. The Loading Condition

The external forces acting on the stabilizing piles include landslide thrust force, active earth pressure, and concentrated pulling force exerted by tension anchors 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. The subgrade reaction force is acting on the pile segment below the slip surface. These external forces considered are in equilibrium.


#### 2.2.2. The Pile-Soil Interaction Model

Due to its simplicity and efficiency of use, the Winkler soil model is also adopted in the current analysis to describe the pile-soil interaction behavior. Winkler method assumed that the substratum is composed of independent horizontal springs. Under the Winkler hypothesis, the soil reaction pressures (p) acting on the pile can be modeled by discrete independent linear or nonlinear springs in form of the following equation:
(3)p=k·w,
where p is the horizontal soil reaction pressure (it has a unit of force/length2).


#### 2.2.3. The New Equilibrium Differential Equations

Considering an isolated free portion of pile, as shown in Figure3, having an infinitesimal length of ds and acted upon by external distributed normal load qn and tangential load qτ. The free segment can be imagined to be cut out of the pile and the internal forces (M, N, and Q) in the original pile may become external forces on the isolated free portion.Figure 3
The free body diagram of infinitesimal isolated segment of pile.We define sign conventions so that the six variables as shown in Figure3 are positive. The sign convention adopted for forces is that positive sign   indicates tensile axial force N, the positive shearing force Q should be directed so that they will tend to rotate the element counterclockwise and the positive bending moment M will tend to make the element be concave leftward. The sign convention adopted for displacements is that the positive normal displacement v points outward normal, the positive tangential displacement u points right when facing outward normal, the positive φ is counterclockwise.Thus, considering the equilibrium of the above infinitesimal pile segment under the action of the applied loads (shown in Figure3), we arrive at two force equilibrium equations in the direction of u and v and one moment equilibrium equation:
(4)(N+dNdsds)-N-qτds-ksuDds=0,(Q+dQdsds)-Q-qnds-knvbds=0,(M+dMdsds)-M+(Q+dQdsds)ds2+Qds2=0.Simplifying the above equations and neglecting the higher order terms and the terms with the square of the differential, the equilibrium equations now take the following form:(5)dNds=ksDu+qτ,dQds=knbv+qn,dMds=-Q,
where  S   is the position coordinate and D is perimeter of the cross-section. b is the width of the cross-section, ks is the Winkler modulus of vertical subgrade reaction, and kn is the Winkler modulus of horizontal subgrade reaction.For the sake of convenience of formula deducing, let(6)X=[NQM];Z=[uvφ];P=[qτqn0].The set of equations of equilibrium (5) can be rewritten as the following matrix form:
(7)dXds=B·X+L·Z+P,
where B=[0000000-10]; L=[ksD000knb0000].


#### 2.2.4. The Geometric and Constitutive Equations

When the deformation (du-, dv-, and dφ-) of the differential element (shown in Figure 3) induced by the internal forces (N, Q, and M) is considered given, the corresponding strains can be expressed as
(8)dU-ds=(du-ds,dv-ds,dφ-ds).According to the related theory of elastic beam, the internal forces (N, Q, and M) can be related to strains as the following linear constitutive equation:
(9)dU-ds=[du-dsdv-dsdφ-ds]=[1EA000αGA0001EI]·[NQM],
where α is the shearing constant related to the shape of pile cross section (α=6/5 for rectangular cross section; α=10/9 for circular cross section), A is the cross-sectional area, and EI is the flexural rigidity of the pile’s cross section.The deformationdU- can be decomposed into two parts. One part is the dZ induced by the displacement on its direction and another part is the projection of other displacements on to this direction which takes the form B-Zds. Where B- is the undetermined third-order square matrix. Then, the deformation dU- can be expressed as
(10)dU-=dZ+B-Zds.Applying the principle of virtual work to the isolated differential element of pile (shown in Figure3), B- can be determined. We suppose that each point of the body is given an infinitesimal virtual displacement δZ satisfying displacement boundary conditions which were prescribed. The virtual deformation associated with the infinitesimal virtual displacement is δU-. The virtual work of the external surface forces is -∫P-T(δZ)ds, where P-=P+LZ. The virtual work of the internal forces is ∫XTd(δU-). By equating the external work to the internal work, we have -∫P-T(δZ)ds=∫XTd(δU-). Substituting (7) and (10) into the above equation and simplifying yields ∫XT(BT-B-)δZds=0. Since this equation is satisfied for arbitrary δZ, the terms in the brackets in the integral must vanish at every point, that means B-=BT. At last, we developed the following geometric and constitutive equations:
(11)duds=NEA,dvds=αQGA+φ,dφds=MEI.


#### 2.2.5. Matrix Form of the Governing Equations

For the sake of convenience of problem solving, combining the three equilibrium differential equations (5) and the three geometric and constitutive equations (11) led to a system of six equations as follows:
(12)dNds=ksDu+qτ,dQds=knbv+qn,dMds=-Q,duds=NEA,dvds=αQGA+φ,dφds=MEJ.
Letting
(13)K=[000ksD000000knb00-100001EA000000αGA0001001EI000];X={NQMuvφ}T;p={qτqn0000}T.The matrix notation is used to present (12) in thefollowing form:
(14)dXds=KX+p=F(X,s),
where K is the coefficients square matrix of six order and X and p represent two column matrix.This system of six independent differential equations can be solved for six unknown functions (three independent forces functions and three independent displacements functions).


## 2.1. Conventional Governing Differential Equation for Pile Deflection

Before deducing the new differential equations governing the arching mechanism of stabilizing piles, we need to review the conventional governing differential equations for pile deflections which can be found in many of the available literature [1, 15]. It is reported here only for the sake of completeness.We assume a planar cartesian coordinate systemw-z, as shown in Figure 1, with its origin at the center of the pile head such that the z-axis coincides with the pile axis and the w-z plane coincides with the plane of the paper. A lateral force Q and a moment M are assumed to be applied at the pile head.Figure 1
The coordinate system of laterally loaded pile model.The pile-soil interaction response can be idealised as a vertical beam placed in a Winkler spring medium. The pile material is assumed to follow linear elastic behaviour.We know that vertical piles resist lateral loads or moments by deflecting until the necessary reaction in the surrounding soil is mobilized.Neglecting the friction force along the pile-soil interface, the deflection response of piles subjected to lateral load and constrained by surrounding Winkler springs is governed by the following differential equation:(1)EId4w(z)dz4+k(z)w(z)b=0,
where z is the position coordinate, w(z) is the deflection function of the pile (it has a unit of length), E is the pile’s Young’s modulus, I is the second moment of inertia of the pile’s cross-section (EI is called the flexural rigidity of the pile), k is the spring constant, also called the modulus of horizontal subgrade reaction (it has a unit of force/length3), which can be assumed to vary linearly with depth for soil or to remain unchanged for rock. Much of the available engineering experience and assumptions can be used to determine its value, and b is the section width.Based on (1), the deflection function w(z) can be obtained by power series solution, finite difference method, and FEM in terms of bar and beam. The detailed solution scheme can be found in many of the available literature. Once the deflection function w(z) is established, the bending moment, shear force, and soil reaction force at various depths along the pile can be deduced by differentiating the deflection function as follows:
(2)φ=dwdz,M=-EId2wdz2,Q=-EId3wdz3.


## 2.2. Derivation of the New Governing Equations

Under the scheme of uncoupled analysis of the pile (as shown in Figure2), 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.Figure 2
The uncoupled analysis model for the stabilizing pile with prestressed tieback anchor.


### 2.2.1. The Loading Condition

The external forces acting on the stabilizing piles include landslide thrust force, active earth pressure, and concentrated pulling force exerted by tension anchors 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. The subgrade reaction force is acting on the pile segment below the slip surface. These external forces considered are in equilibrium.


### 2.2.2. The Pile-Soil Interaction Model

Due to its simplicity and efficiency of use, the Winkler soil model is also adopted in the current analysis to describe the pile-soil interaction behavior. Winkler method assumed that the substratum is composed of independent horizontal springs. Under the Winkler hypothesis, the soil reaction pressures (p) acting on the pile can be modeled by discrete independent linear or nonlinear springs in form of the following equation:
(3)p=k·w,
where p is the horizontal soil reaction pressure (it has a unit of force/length2).


### 2.2.3. The New Equilibrium Differential Equations

Considering an isolated free portion of pile, as shown in Figure3, having an infinitesimal length of ds and acted upon by external distributed normal load qn and tangential load qτ. The free segment can be imagined to be cut out of the pile and the internal forces (M, N, and Q) in the original pile may become external forces on the isolated free portion.Figure 3
The free body diagram of infinitesimal isolated segment of pile.We define sign conventions so that the six variables as shown in Figure3 are positive. The sign convention adopted for forces is that positive sign   indicates tensile axial force N, the positive shearing force Q should be directed so that they will tend to rotate the element counterclockwise and the positive bending moment M will tend to make the element be concave leftward. The sign convention adopted for displacements is that the positive normal displacement v points outward normal, the positive tangential displacement u points right when facing outward normal, the positive φ is counterclockwise.Thus, considering the equilibrium of the above infinitesimal pile segment under the action of the applied loads (shown in Figure3), we arrive at two force equilibrium equations in the direction of u and v and one moment equilibrium equation:
(4)(N+dNdsds)-N-qτds-ksuDds=0,(Q+dQdsds)-Q-qnds-knvbds=0,(M+dMdsds)-M+(Q+dQdsds)ds2+Qds2=0.Simplifying the above equations and neglecting the higher order terms and the terms with the square of the differential, the equilibrium equations now take the following form:(5)dNds=ksDu+qτ,dQds=knbv+qn,dMds=-Q,
where  S   is the position coordinate and D is perimeter of the cross-section. b is the width of the cross-section, ks is the Winkler modulus of vertical subgrade reaction, and kn is the Winkler modulus of horizontal subgrade reaction.For the sake of convenience of formula deducing, let(6)X=[NQM];Z=[uvφ];P=[qτqn0].The set of equations of equilibrium (5) can be rewritten as the following matrix form:
(7)dXds=B·X+L·Z+P,
where B=[0000000-10]; L=[ksD000knb0000].


### 2.2.4. The Geometric and Constitutive Equations

When the deformation (du-, dv-, and dφ-) of the differential element (shown in Figure 3) induced by the internal forces (N, Q, and M) is considered given, the corresponding strains can be expressed as
(8)dU-ds=(du-ds,dv-ds,dφ-ds).According to the related theory of elastic beam, the internal forces (N, Q, and M) can be related to strains as the following linear constitutive equation:
(9)dU-ds=[du-dsdv-dsdφ-ds]=[1EA000αGA0001EI]·[NQM],
where α is the shearing constant related to the shape of pile cross section (α=6/5 for rectangular cross section; α=10/9 for circular cross section), A is the cross-sectional area, and EI is the flexural rigidity of the pile’s cross section.The deformationdU- can be decomposed into two parts. One part is the dZ induced by the displacement on its direction and another part is the projection of other displacements on to this direction which takes the form B-Zds. Where B- is the undetermined third-order square matrix. Then, the deformation dU- can be expressed as
(10)dU-=dZ+B-Zds.Applying the principle of virtual work to the isolated differential element of pile (shown in Figure3), B- can be determined. We suppose that each point of the body is given an infinitesimal virtual displacement δZ satisfying displacement boundary conditions which were prescribed. The virtual deformation associated with the infinitesimal virtual displacement is δU-. The virtual work of the external surface forces is -∫P-T(δZ)ds, where P-=P+LZ. The virtual work of the internal forces is ∫XTd(δU-). By equating the external work to the internal work, we have -∫P-T(δZ)ds=∫XTd(δU-). Substituting (7) and (10) into the above equation and simplifying yields ∫XT(BT-B-)δZds=0. Since this equation is satisfied for arbitrary δZ, the terms in the brackets in the integral must vanish at every point, that means B-=BT. At last, we developed the following geometric and constitutive equations:
(11)duds=NEA,dvds=αQGA+φ,dφds=MEI.


### 2.2.5. Matrix Form of the Governing Equations

For the sake of convenience of problem solving, combining the three equilibrium differential equations (5) and the three geometric and constitutive equations (11) led to a system of six equations as follows:
(12)dNds=ksDu+qτ,dQds=knbv+qn,dMds=-Q,duds=NEA,dvds=αQGA+φ,dφds=MEJ.
Letting
(13)K=[000ksD000000knb00-100001EA000000αGA0001001EI000];X={NQMuvφ}T;p={qτqn0000}T.The matrix notation is used to present (12) in thefollowing form:
(14)dXds=KX+p=F(X,s),
where K is the coefficients square matrix of six order and X and p represent two column matrix.This system of six independent differential equations can be solved for six unknown functions (three independent forces functions and three independent displacements functions).


## 2.2.1. The Loading Condition

The external forces acting on the stabilizing piles include landslide thrust force, active earth pressure, and concentrated pulling force exerted by tension anchors 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. The subgrade reaction force is acting on the pile segment below the slip surface. These external forces considered are in equilibrium.


## 2.2.2. The Pile-Soil Interaction Model

Due to its simplicity and efficiency of use, the Winkler soil model is also adopted in the current analysis to describe the pile-soil interaction behavior. Winkler method assumed that the substratum is composed of independent horizontal springs. Under the Winkler hypothesis, the soil reaction pressures (p) acting on the pile can be modeled by discrete independent linear or nonlinear springs in form of the following equation:
(3)p=k·w,
where p is the horizontal soil reaction pressure (it has a unit of force/length2).


## 2.2.3. The New Equilibrium Differential Equations

Considering an isolated free portion of pile, as shown in Figure3, having an infinitesimal length of ds and acted upon by external distributed normal load qn and tangential load qτ. The free segment can be imagined to be cut out of the pile and the internal forces (M, N, and Q) in the original pile may become external forces on the isolated free portion.Figure 3
The free body diagram of infinitesimal isolated segment of pile.We define sign conventions so that the six variables as shown in Figure3 are positive. The sign convention adopted for forces is that positive sign   indicates tensile axial force N, the positive shearing force Q should be directed so that they will tend to rotate the element counterclockwise and the positive bending moment M will tend to make the element be concave leftward. The sign convention adopted for displacements is that the positive normal displacement v points outward normal, the positive tangential displacement u points right when facing outward normal, the positive φ is counterclockwise.Thus, considering the equilibrium of the above infinitesimal pile segment under the action of the applied loads (shown in Figure3), we arrive at two force equilibrium equations in the direction of u and v and one moment equilibrium equation:
(4)(N+dNdsds)-N-qτds-ksuDds=0,(Q+dQdsds)-Q-qnds-knvbds=0,(M+dMdsds)-M+(Q+dQdsds)ds2+Qds2=0.Simplifying the above equations and neglecting the higher order terms and the terms with the square of the differential, the equilibrium equations now take the following form:(5)dNds=ksDu+qτ,dQds=knbv+qn,dMds=-Q,
where  S   is the position coordinate and D is perimeter of the cross-section. b is the width of the cross-section, ks is the Winkler modulus of vertical subgrade reaction, and kn is the Winkler modulus of horizontal subgrade reaction.For the sake of convenience of formula deducing, let(6)X=[NQM];Z=[uvφ];P=[qτqn0].The set of equations of equilibrium (5) can be rewritten as the following matrix form:
(7)dXds=B·X+L·Z+P,
where B=[0000000-10]; L=[ksD000knb0000].


## 2.2.4. The Geometric and Constitutive Equations

When the deformation (du-, dv-, and dφ-) of the differential element (shown in Figure 3) induced by the internal forces (N, Q, and M) is considered given, the corresponding strains can be expressed as
(8)dU-ds=(du-ds,dv-ds,dφ-ds).According to the related theory of elastic beam, the internal forces (N, Q, and M) can be related to strains as the following linear constitutive equation:
(9)dU-ds=[du-dsdv-dsdφ-ds]=[1EA000αGA0001EI]·[NQM],
where α is the shearing constant related to the shape of pile cross section (α=6/5 for rectangular cross section; α=10/9 for circular cross section), A is the cross-sectional area, and EI is the flexural rigidity of the pile’s cross section.The deformationdU- can be decomposed into two parts. One part is the dZ induced by the displacement on its direction and another part is the projection of other displacements on to this direction which takes the form B-Zds. Where B- is the undetermined third-order square matrix. Then, the deformation dU- can be expressed as
(10)dU-=dZ+B-Zds.Applying the principle of virtual work to the isolated differential element of pile (shown in Figure3), B- can be determined. We suppose that each point of the body is given an infinitesimal virtual displacement δZ satisfying displacement boundary conditions which were prescribed. The virtual deformation associated with the infinitesimal virtual displacement is δU-. The virtual work of the external surface forces is -∫P-T(δZ)ds, where P-=P+LZ. The virtual work of the internal forces is ∫XTd(δU-). By equating the external work to the internal work, we have -∫P-T(δZ)ds=∫XTd(δU-). Substituting (7) and (10) into the above equation and simplifying yields ∫XT(BT-B-)δZds=0. Since this equation is satisfied for arbitrary δZ, the terms in the brackets in the integral must vanish at every point, that means B-=BT. At last, we developed the following geometric and constitutive equations:
(11)duds=NEA,dvds=αQGA+φ,dφds=MEI.


## 2.2.5. Matrix Form of the Governing Equations

For the sake of convenience of problem solving, combining the three equilibrium differential equations (5) and the three geometric and constitutive equations (11) led to a system of six equations as follows:
(12)dNds=ksDu+qτ,dQds=knbv+qn,dMds=-Q,duds=NEA,dvds=αQGA+φ,dφds=MEJ.
Letting
(13)K=[000ksD000000knb00-100001EA000000αGA0001001EI000];X={NQMuvφ}T;p={qτqn0000}T.The matrix notation is used to present (12) in thefollowing form:
(14)dXds=KX+p=F(X,s),
where K is the coefficients square matrix of six order and X and p represent two column matrix.This system of six independent differential equations can be solved for six unknown functions (three independent forces functions and three independent displacements functions).


## 3. Boundary Conditions

As mentioned above, there is a total of six unknown functions to be determined (N, Q, M, u, v, and φ). Therefore, six boundary conditions are needed for the problem solving.The boundary conditions for (14) are determined according to the way in which the pile’s endpoints 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:
(15)CX-|S=0=o,DX-|S=L=o,
where S=0 indicates the beginning point of calculation (the pile’s base point), S=L indicates the end point of caculation (the pile’s head point), C is the matrix of boundary condition on the beginning point, and D is the matrix of boundary condition on the end point. They are 3 × 6 matrix. In this study, the following possible pile head or base conditions were considered.


### 3.1. Top End Condition

(
1
) In case the anchor is locked against the top end of the pile, pile head is tied back by prestressed anchor (as shown in Figure 4). The corresponding matrix of this type of boundary condition can be deduced as followsFigure 4
The scheme of pile with top end tied back by a prestressed anchor.As shown in Figure4, the local coordinate system (u, v) at the pile head rotates θ degree to become the (u′,v′) coordinate system along the anchor line. According to the principle of coordinate transformation, the following relationship exists:
(16)[u′v′]=[cosθ-sinθsinθcosθ][uv]⟹u′=cosθ·u-sinθ·v,[N′Q′]=[cosθ-sinθsinθcosθ][NQ]⟹{N′=cosθ·N-sinθ·Q=-u′LEaAa+TPQ′=sinθ·N+cosθ·Q=0,
where θ is the rotating angle of anchor as shown in Figure 4, L is the length of unbonded segment of the anchor, EaAa is the product of the cross-sectional area and the elastic young’s modulus of anchor material, and  TP  is the prepull axial force exerted in the anchor.Then, another conditionM=0 is added. So we have in total three conditions. Substituting the expression of u′ into the above equation and simplifying yields
(17)[cosθ-sinθ0cosθEaAaL-sinθEaAaL0sinθcosθ0000001000]×{NQMuvφ}={TP00}.In case the anchors are not joined to the pile at the top end of the pile, the pile head is free and the anchor pulling point needs to be processed using other special mathematical techniques.(
2
) Free head (allows both displacement and rotation), N=Q=0, M=0. The corresponding matrix of boundary condition is
(18)[100000010000001000]×{NQMuvφ}=o.


### 3.2. Bottom End Condition

(
1
) In case of the bottom end hinged (allows rotation without displacement), u=v=0, M=0. The corresponding matrix of boundary condition is
(19)[001000000100000010]×{NQMuvφ}=o.(
2
) In case of the bottom end fixed (allows neither displacement nor rotation), u=v=0, φ=0. The corresponding matrix of boundary condition is
(20)[000100000010000001]×{NQMuvφ}=o.(
3
) In case of the bottom end partially hinged (allows rotation without vertical displacement), u=0, Q=0, M=0. The corresponding matrix of boundary condition is
(21)[010000001000000100]×{NQMuvφ}=o.(
4
) Elastic vertical support on the bottom end, N=Kv·u·area, Q=0, M=0. The corresponding matrix of boundary condition is
(22)[100-Kv·Area00010000001000]×{NQMuvφ}=o,
where Kv is the modulus of vertical compressibility. Area is the cross-sectional area of the pile.Next, we impose the boundary conditions (15) at the pile top and base point upon the derived new governing differential equation (14) to define a boundary-value problem of the following equations:
(23)dXds=KX+p,CX|S=0=o,DX|S=L=o.To this end, the response of anchored stabilizing pile is mathematically idealised as the boundary-value problem of (23). Thus, many numerical methods to solve the ordinary differential equations can be adopted to solve (23).It should be pointed out that the existence and uniqueness of a solution for (23) should be mathematically proved. But this matter is outside the scope of the authors’ major. We can only imagine that the solution exists and is unique according to the physical character of the problem. The solution can be validated through comparative studies.


## 3.1. Top End Condition

(
1
) In case the anchor is locked against the top end of the pile, pile head is tied back by prestressed anchor (as shown in Figure 4). The corresponding matrix of this type of boundary condition can be deduced as followsFigure 4
The scheme of pile with top end tied back by a prestressed anchor.As shown in Figure4, the local coordinate system (u, v) at the pile head rotates θ degree to become the (u′,v′) coordinate system along the anchor line. According to the principle of coordinate transformation, the following relationship exists:
(16)[u′v′]=[cosθ-sinθsinθcosθ][uv]⟹u′=cosθ·u-sinθ·v,[N′Q′]=[cosθ-sinθsinθcosθ][NQ]⟹{N′=cosθ·N-sinθ·Q=-u′LEaAa+TPQ′=sinθ·N+cosθ·Q=0,
where θ is the rotating angle of anchor as shown in Figure 4, L is the length of unbonded segment of the anchor, EaAa is the product of the cross-sectional area and the elastic young’s modulus of anchor material, and  TP  is the prepull axial force exerted in the anchor.Then, another conditionM=0 is added. So we have in total three conditions. Substituting the expression of u′ into the above equation and simplifying yields
(17)[cosθ-sinθ0cosθEaAaL-sinθEaAaL0sinθcosθ0000001000]×{NQMuvφ}={TP00}.In case the anchors are not joined to the pile at the top end of the pile, the pile head is free and the anchor pulling point needs to be processed using other special mathematical techniques.(
2
) Free head (allows both displacement and rotation), N=Q=0, M=0. The corresponding matrix of boundary condition is
(18)[100000010000001000]×{NQMuvφ}=o.


## 3.2. Bottom End Condition

(
1
) In case of the bottom end hinged (allows rotation without displacement), u=v=0, M=0. The corresponding matrix of boundary condition is
(19)[001000000100000010]×{NQMuvφ}=o.(
2
) In case of the bottom end fixed (allows neither displacement nor rotation), u=v=0, φ=0. The corresponding matrix of boundary condition is
(20)[000100000010000001]×{NQMuvφ}=o.(
3
) In case of the bottom end partially hinged (allows rotation without vertical displacement), u=0, Q=0, M=0. The corresponding matrix of boundary condition is
(21)[010000001000000100]×{NQMuvφ}=o.(
4
) Elastic vertical support on the bottom end, N=Kv·u·area, Q=0, M=0. The corresponding matrix of boundary condition is
(22)[100-Kv·Area00010000001000]×{NQMuvφ}=o,
where Kv is the modulus of vertical compressibility. Area is the cross-sectional area of the pile.Next, we impose the boundary conditions (15) at the pile top and base point upon the derived new governing differential equation (14) to define a boundary-value problem of the following equations:
(23)dXds=KX+p,CX|S=0=o,DX|S=L=o.To this end, the response of anchored stabilizing pile is mathematically idealised as the boundary-value problem of (23). Thus, many numerical methods to solve the ordinary differential equations can be adopted to solve (23).It should be pointed out that the existence and uniqueness of a solution for (23) should be mathematically proved. But this matter is outside the scope of the authors’ major. We can only imagine that the solution exists and is unique according to the physical character of the problem. The solution can be validated through comparative studies.


## 4. Numerical Solution Scheme

### 4.1. Uniformity Preprocessing for the Order of Magnitude of the Element in the Coefficient MatrixK

The magnitude order of the section internal forces (N, Q, and M) are so much higher than that of the displacements (u, v, and φ) that the numerical solving of the equations may meet singularity difficulty. So for reasons of numerical stability, it is needed to perform uniformity preprocessing for the magnitude order of the element in the coefficient matrix K. We multiply the displacements (u, v, and φ) by E and substitute the original displacement variables by the expressions (Eu, Ev, Eφ). Then, we can redefine two variables as follows:
(24)X~={NQMu-v-φ-}={NQMEuEvEφ};K~=[000HksDE000000HknbE00-100001A000000EαGA0001001I000].Finally, we arrive at the following system of ordinary differential equations:(25)dX~ds=K~X~+p,C·X~|S=0=o,D·X~|S=L=o.This system of six independent differential equations subjected to boundary conditions can be numerically solved for six unknown functions (three forces and three displacements).


### 4.2. The Finite Difference Method

According to the theory of numerical analysis, the Runge-Kutta algorithm is a convenient, powerful, and high-accuracy procedure for solving the ordinary differential equation of the formdX/ds=F(X,s). So, it is adopted to solve (25).


#### 4.2.1. Derivation of the Recursion Formula

The following finite difference formula (26) is one format of the Runge-Kutta methods:
(26)Xn+1=Xn+16(K1+2K2+2K3+K4),K1=δ·F(Xn,sn),K2=δ·F(Xn+12K1,sn+12δ),K3=δ·F(Xn+12K2,sn+12δ),K4=δ·F(Xn+K3,sn+δ).The Runge-Kutta algorithm of this type (26) is a numerical method of fifth order. Where δ is the step size of difference and F(X,s)=K-X+p. For convenience of computation, this formulation may be rewritten in the following form:
(27)Xn+1=GnXn+Hnδ.As we can see, the value of function at point (n+1) can be determined from the value of function at point (n). where n=0 represent the beginning point of caculation and n=m represent the end point of caculation.The following recursion formula can be used to obtain the value ofGn and Hn in (27):
(28)Kn(j)=K~(sn+δj),Pn(j)=p(sn+δj),δ1=0,δ2=δ3=12δ,δ4=δ,(29)Gn=I+∑j=14βjG(j),Hn=∑j=14βjH(j),G(j)=(δK(j))(I+αjG(j-1)),G(0)=0,H(j)=(δK(j))αjH(j-1)+P(j),H(0)=0,α1=α2=α3=12,α4=1,β1=β4=16,β2=β3=13,
where I is the identity matrix.


#### 4.2.2. Determination of the Initial VectorX0

The initial value is the start point of the recursion formula. Now we discuss in the following how to obtain the initial vectorX0 by using the recursion formula of (27) and imposing the boundary conditions at piles head and base.Considering the recursion formula of (27), the Xn can be expressed in terms of X0 as follows:
(30)Xn=D(n)X0+F(n).In the case ofn=0, we have D(0)=I, F(0)=0 and substituting it into the recursion formula of (27). We get
(31)Xn+1=Gn(D(n)X0+F(n))+Hnδ.It can be rewritten as follows:(32)Xn+1=(GnD(n))X0+(GnF(n)+Hnδ),Xn=D(n)X0+F(n)⟹Xn+1=D(n+1)X0+F(n+1),Xn+1=(GnD(n))X0+(GnF(n)+Hnδ),Xn+1=D(n+1)X0+F(n+1).Comparing the above two equations, the recursion formula forD(n) and F(n) is obtained as follows:
(33)D(0)=I,D(n+1)=GnD(n),F(0)=o,F(n+1)=GnF(n)+Hnδ.Now considering the case of boundary point,X
m
=
D
(
m
)
X
0
+
F
(
m
), we substitute the boundary conditions at end point CX0=o, DXm=o into the above equation. This leads to the equation to solve for X0(34)[CDD(m)]X0=[o-DF(m)].The above system of linear algebraic equations can be solved forX0 by using the method of Gaussian elimination with pivot selection. Once X0 is known, Xn can be obtained in sequence using the recursion formula of (27).


#### 4.2.3. Procedure to Deal with the Intersection Point of Pile and Anchors

In some cases where the pile is tied back with multiple level anchors, the intersection point of pile and anchor is not located at the pile top (as shown in Figure5). So the effect of anchor’s pull cannot be considered as boundary condition. The pile’s internal force discontinuity induced by anchor’s oblique pull loads needs to be processed in the flow process of the finite difference method. The procedure is described as follows.Figure 5
Scheme of local coordinate system transformation.We establish two local coordinate systems on the intersection point (as shown in Figure5) to help calculate the axial force of the anchors.According to Figure5, the coordinate system transformation relationship can be expressed as
(35)[u′v′]=[cosθ-sinθsinθcosθ][uv]⟹u′=cosθ·u-sinθ·v.
So the axial tension force of anchors can be expressed as
(36)N′=-u′LEA+TP=(cosθ·u-sinθ·v)(-EaAa)L+TP,
where θ is the rotating angle of anchor as shown in Figure 5 and TP is the prepull force exerted in the anchor.In the node numbering system,n is assumed to be the intersection point of pile and anchors (as shown in Figure 6). The procedure to go across node n is described as follows.Figure 6
Scheme of local coordinate system transformation across the points of discontinuity.Because noden is subjected to oblique pull loads of anchor (as shown in Figure 6), the discontinuity of shear force, axial force across node n exists which can be expressed as follows:
(37)N-n=Nn+N′·cos(π-θ),Q-n=Qn+N′·sin(π-θ).Because the anchor cannot resist the bending moment, the transfer of moment at noden is not affected by the anchor. The deformation transfer across node n is also continuous. Substituting the expression of N′ into the above equation and simplifying yields
(38)[N-nQ-nM-nu-nv-nφ-n]=Cn×[NnQnMnunvnφn]+F→P,
where
(39)Cn=[100cosθ-EALcos(π-θ)-sinθ-EALcos(π-θ)0010cosθ-EALsin(π-θ)-sinθ-EALsin(π-θ)0001000000100000010000001];F→P=[TPcos(π-θ)TPsin(π-θ)0000].For convenience of formular expressing, (38) may be abbreviated as
(40)X-n=Cn·Xn+F→P,
where Cn is termed as connecting matrix. The recursion process relating to node n is described as follows.When using the recursion formular to calculateD(m), F(m), and Xi on node n, the connecting matrix Cn must be used as in the following steps:(
1
) the step of recursion caculation of Xi(41)Xi+1=GiXi+Hiδ(i=1,2,…,n-1),X-n=CnXn+F⃑P=Cn(Gn-1Xn-1+Hn-1δ),Xn+1=GnX-n+Hnδ.
The following process is the same as other nodes’:
(42)Xi+1=GiXi+Hiδ(i=n+1,n+2,…m).(
2
) The step of recursion caculation of D(m) and F(m).BecauseXn=D(n)X0+F(n), then we get
(43)Xn+1=GnX-n+Hnδ=Gn(CnXn+F→P)+Hnδ=Gn[Cn(D(n)X0+F(n))+F→P]+Hnδ=GnCnD(n)X0+GnCnF(n)+GnF⃑P+Hnδ=D(n+1)X0+F(n+1)⟹{D(n+1)=GnCnD(n)F(n+1)=GnCnF(n)+GnF⃑P+Hnδ.The recursion steps to calculateD(m) are as follows:
(44)D(i+1)=GiD(i)(i=1,2,…,n-1),D(n+1)=GnCnD(n),D(i+1)=GiD(i)(i=n+1,n+2,…m).The recursion steps to calculateF(m) is as follows:
(45)F(i+1)=GiF(i)+Hiδ(i=1,2,…,n-1),F(n+1)=GnCnF(n)+GnF⃑P+Hnδ,F(i+1)=GiF(i)+Hiδ(i=n+1,n+2,…m).


#### 4.2.4. The Solution Flow Process

In short, the proposed solution procedure involves the following four main steps:(1)
calculating the value ofGn and Hn using the given (29),(2)
calculatingD(m) and F(m) using the given recursion formula of (33),(3)
calculating the vectorX0 by solving linear algebraic equation (34),(4)
calculatingXn using the given recursion formula of (27).Because the equations and solution formula are all given in form of matrix. 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.


## 4.1. Uniformity Preprocessing for the Order of Magnitude of the Element in the Coefficient MatrixK

The magnitude order of the section internal forces (N, Q, and M) are so much higher than that of the displacements (u, v, and φ) that the numerical solving of the equations may meet singularity difficulty. So for reasons of numerical stability, it is needed to perform uniformity preprocessing for the magnitude order of the element in the coefficient matrix K. We multiply the displacements (u, v, and φ) by E and substitute the original displacement variables by the expressions (Eu, Ev, Eφ). Then, we can redefine two variables as follows:
(24)X~={NQMu-v-φ-}={NQMEuEvEφ};K~=[000HksDE000000HknbE00-100001A000000EαGA0001001I000].Finally, we arrive at the following system of ordinary differential equations:(25)dX~ds=K~X~+p,C·X~|S=0=o,D·X~|S=L=o.This system of six independent differential equations subjected to boundary conditions can be numerically solved for six unknown functions (three forces and three displacements).


## 4.2. The Finite Difference Method

According to the theory of numerical analysis, the Runge-Kutta algorithm is a convenient, powerful, and high-accuracy procedure for solving the ordinary differential equation of the formdX/ds=F(X,s). So, it is adopted to solve (25).


### 4.2.1. Derivation of the Recursion Formula

The following finite difference formula (26) is one format of the Runge-Kutta methods:
(26)Xn+1=Xn+16(K1+2K2+2K3+K4),K1=δ·F(Xn,sn),K2=δ·F(Xn+12K1,sn+12δ),K3=δ·F(Xn+12K2,sn+12δ),K4=δ·F(Xn+K3,sn+δ).The Runge-Kutta algorithm of this type (26) is a numerical method of fifth order. Where δ is the step size of difference and F(X,s)=K-X+p. For convenience of computation, this formulation may be rewritten in the following form:
(27)Xn+1=GnXn+Hnδ.As we can see, the value of function at point (n+1) can be determined from the value of function at point (n). where n=0 represent the beginning point of caculation and n=m represent the end point of caculation.The following recursion formula can be used to obtain the value ofGn and Hn in (27):
(28)Kn(j)=K~(sn+δj),Pn(j)=p(sn+δj),δ1=0,δ2=δ3=12δ,δ4=δ,(29)Gn=I+∑j=14βjG(j),Hn=∑j=14βjH(j),G(j)=(δK(j))(I+αjG(j-1)),G(0)=0,H(j)=(δK(j))αjH(j-1)+P(j),H(0)=0,α1=α2=α3=12,α4=1,β1=β4=16,β2=β3=13,
where I is the identity matrix.


### 4.2.2. Determination of the Initial VectorX0

The initial value is the start point of the recursion formula. Now we discuss in the following how to obtain the initial vectorX0 by using the recursion formula of (27) and imposing the boundary conditions at piles head and base.Considering the recursion formula of (27), the Xn can be expressed in terms of X0 as follows:
(30)Xn=D(n)X0+F(n).In the case ofn=0, we have D(0)=I, F(0)=0 and substituting it into the recursion formula of (27). We get
(31)Xn+1=Gn(D(n)X0+F(n))+Hnδ.It can be rewritten as follows:(32)Xn+1=(GnD(n))X0+(GnF(n)+Hnδ),Xn=D(n)X0+F(n)⟹Xn+1=D(n+1)X0+F(n+1),Xn+1=(GnD(n))X0+(GnF(n)+Hnδ),Xn+1=D(n+1)X0+F(n+1).Comparing the above two equations, the recursion formula forD(n) and F(n) is obtained as follows:
(33)D(0)=I,D(n+1)=GnD(n),F(0)=o,F(n+1)=GnF(n)+Hnδ.Now considering the case of boundary point,X
m
=
D
(
m
)
X
0
+
F
(
m
), we substitute the boundary conditions at end point CX0=o, DXm=o into the above equation. This leads to the equation to solve for X0(34)[CDD(m)]X0=[o-DF(m)].The above system of linear algebraic equations can be solved forX0 by using the method of Gaussian elimination with pivot selection. Once X0 is known, Xn can be obtained in sequence using the recursion formula of (27).


### 4.2.3. Procedure to Deal with the Intersection Point of Pile and Anchors

In some cases where the pile is tied back with multiple level anchors, the intersection point of pile and anchor is not located at the pile top (as shown in Figure5). So the effect of anchor’s pull cannot be considered as boundary condition. The pile’s internal force discontinuity induced by anchor’s oblique pull loads needs to be processed in the flow process of the finite difference method. The procedure is described as follows.Figure 5
Scheme of local coordinate system transformation.We establish two local coordinate systems on the intersection point (as shown in Figure5) to help calculate the axial force of the anchors.According to Figure5, the coordinate system transformation relationship can be expressed as
(35)[u′v′]=[cosθ-sinθsinθcosθ][uv]⟹u′=cosθ·u-sinθ·v.
So the axial tension force of anchors can be expressed as
(36)N′=-u′LEA+TP=(cosθ·u-sinθ·v)(-EaAa)L+TP,
where θ is the rotating angle of anchor as shown in Figure 5 and TP is the prepull force exerted in the anchor.In the node numbering system,n is assumed to be the intersection point of pile and anchors (as shown in Figure 6). The procedure to go across node n is described as follows.Figure 6
Scheme of local coordinate system transformation across the points of discontinuity.Because noden is subjected to oblique pull loads of anchor (as shown in Figure 6), the discontinuity of shear force, axial force across node n exists which can be expressed as follows:
(37)N-n=Nn+N′·cos(π-θ),Q-n=Qn+N′·sin(π-θ).Because the anchor cannot resist the bending moment, the transfer of moment at noden is not affected by the anchor. The deformation transfer across node n is also continuous. Substituting the expression of N′ into the above equation and simplifying yields
(38)[N-nQ-nM-nu-nv-nφ-n]=Cn×[NnQnMnunvnφn]+F→P,
where
(39)Cn=[100cosθ-EALcos(π-θ)-sinθ-EALcos(π-θ)0010cosθ-EALsin(π-θ)-sinθ-EALsin(π-θ)0001000000100000010000001];F→P=[TPcos(π-θ)TPsin(π-θ)0000].For convenience of formular expressing, (38) may be abbreviated as
(40)X-n=Cn·Xn+F→P,
where Cn is termed as connecting matrix. The recursion process relating to node n is described as follows.When using the recursion formular to calculateD(m), F(m), and Xi on node n, the connecting matrix Cn must be used as in the following steps:(
1
) the step of recursion caculation of Xi(41)Xi+1=GiXi+Hiδ(i=1,2,…,n-1),X-n=CnXn+F⃑P=Cn(Gn-1Xn-1+Hn-1δ),Xn+1=GnX-n+Hnδ.
The following process is the same as other nodes’:
(42)Xi+1=GiXi+Hiδ(i=n+1,n+2,…m).(
2
) The step of recursion caculation of D(m) and F(m).BecauseXn=D(n)X0+F(n), then we get
(43)Xn+1=GnX-n+Hnδ=Gn(CnXn+F→P)+Hnδ=Gn[Cn(D(n)X0+F(n))+F→P]+Hnδ=GnCnD(n)X0+GnCnF(n)+GnF⃑P+Hnδ=D(n+1)X0+F(n+1)⟹{D(n+1)=GnCnD(n)F(n+1)=GnCnF(n)+GnF⃑P+Hnδ.The recursion steps to calculateD(m) are as follows:
(44)D(i+1)=GiD(i)(i=1,2,…,n-1),D(n+1)=GnCnD(n),D(i+1)=GiD(i)(i=n+1,n+2,…m).The recursion steps to calculateF(m) is as follows:
(45)F(i+1)=GiF(i)+Hiδ(i=1,2,…,n-1),F(n+1)=GnCnF(n)+GnF⃑P+Hnδ,F(i+1)=GiF(i)+Hiδ(i=n+1,n+2,…m).


### 4.2.4. The Solution Flow Process

In short, the proposed solution procedure involves the following four main steps:(1)
calculating the value ofGn and Hn using the given (29),(2)
calculatingD(m) and F(m) using the given recursion formula of (33),(3)
calculating the vectorX0 by solving linear algebraic equation (34),(4)
calculatingXn using the given recursion formula of (27).Because the equations and solution formula are all given in form of matrix. 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.


## 4.2.1. Derivation of the Recursion Formula

The following finite difference formula (26) is one format of the Runge-Kutta methods:
(26)Xn+1=Xn+16(K1+2K2+2K3+K4),K1=δ·F(Xn,sn),K2=δ·F(Xn+12K1,sn+12δ),K3=δ·F(Xn+12K2,sn+12δ),K4=δ·F(Xn+K3,sn+δ).The Runge-Kutta algorithm of this type (26) is a numerical method of fifth order. Where δ is the step size of difference and F(X,s)=K-X+p. For convenience of computation, this formulation may be rewritten in the following form:
(27)Xn+1=GnXn+Hnδ.As we can see, the value of function at point (n+1) can be determined from the value of function at point (n). where n=0 represent the beginning point of caculation and n=m represent the end point of caculation.The following recursion formula can be used to obtain the value ofGn and Hn in (27):
(28)Kn(j)=K~(sn+δj),Pn(j)=p(sn+δj),δ1=0,δ2=δ3=12δ,δ4=δ,(29)Gn=I+∑j=14βjG(j),Hn=∑j=14βjH(j),G(j)=(δK(j))(I+αjG(j-1)),G(0)=0,H(j)=(δK(j))αjH(j-1)+P(j),H(0)=0,α1=α2=α3=12,α4=1,β1=β4=16,β2=β3=13,
where I is the identity matrix.


## 4.2.2. Determination of the Initial VectorX0

The initial value is the start point of the recursion formula. Now we discuss in the following how to obtain the initial vectorX0 by using the recursion formula of (27) and imposing the boundary conditions at piles head and base.Considering the recursion formula of (27), the Xn can be expressed in terms of X0 as follows:
(30)Xn=D(n)X0+F(n).In the case ofn=0, we have D(0)=I, F(0)=0 and substituting it into the recursion formula of (27). We get
(31)Xn+1=Gn(D(n)X0+F(n))+Hnδ.It can be rewritten as follows:(32)Xn+1=(GnD(n))X0+(GnF(n)+Hnδ),Xn=D(n)X0+F(n)⟹Xn+1=D(n+1)X0+F(n+1),Xn+1=(GnD(n))X0+(GnF(n)+Hnδ),Xn+1=D(n+1)X0+F(n+1).Comparing the above two equations, the recursion formula forD(n) and F(n) is obtained as follows:
(33)D(0)=I,D(n+1)=GnD(n),F(0)=o,F(n+1)=GnF(n)+Hnδ.Now considering the case of boundary point,X
m
=
D
(
m
)
X
0
+
F
(
m
), we substitute the boundary conditions at end point CX0=o, DXm=o into the above equation. This leads to the equation to solve for X0(34)[CDD(m)]X0=[o-DF(m)].The above system of linear algebraic equations can be solved forX0 by using the method of Gaussian elimination with pivot selection. Once X0 is known, Xn can be obtained in sequence using the recursion formula of (27).


## 4.2.3. Procedure to Deal with the Intersection Point of Pile and Anchors

In some cases where the pile is tied back with multiple level anchors, the intersection point of pile and anchor is not located at the pile top (as shown in Figure5). So the effect of anchor’s pull cannot be considered as boundary condition. The pile’s internal force discontinuity induced by anchor’s oblique pull loads needs to be processed in the flow process of the finite difference method. The procedure is described as follows.Figure 5
Scheme of local coordinate system transformation.We establish two local coordinate systems on the intersection point (as shown in Figure5) to help calculate the axial force of the anchors.According to Figure5, the coordinate system transformation relationship can be expressed as
(35)[u′v′]=[cosθ-sinθsinθcosθ][uv]⟹u′=cosθ·u-sinθ·v.
So the axial tension force of anchors can be expressed as
(36)N′=-u′LEA+TP=(cosθ·u-sinθ·v)(-EaAa)L+TP,
where θ is the rotating angle of anchor as shown in Figure 5 and TP is the prepull force exerted in the anchor.In the node numbering system,n is assumed to be the intersection point of pile and anchors (as shown in Figure 6). The procedure to go across node n is described as follows.Figure 6
Scheme of local coordinate system transformation across the points of discontinuity.Because noden is subjected to oblique pull loads of anchor (as shown in Figure 6), the discontinuity of shear force, axial force across node n exists which can be expressed as follows:
(37)N-n=Nn+N′·cos(π-θ),Q-n=Qn+N′·sin(π-θ).Because the anchor cannot resist the bending moment, the transfer of moment at noden is not affected by the anchor. The deformation transfer across node n is also continuous. Substituting the expression of N′ into the above equation and simplifying yields
(38)[N-nQ-nM-nu-nv-nφ-n]=Cn×[NnQnMnunvnφn]+F→P,
where
(39)Cn=[100cosθ-EALcos(π-θ)-sinθ-EALcos(π-θ)0010cosθ-EALsin(π-θ)-sinθ-EALsin(π-θ)0001000000100000010000001];F→P=[TPcos(π-θ)TPsin(π-θ)0000].For convenience of formular expressing, (38) may be abbreviated as
(40)X-n=Cn·Xn+F→P,
where Cn is termed as connecting matrix. The recursion process relating to node n is described as follows.When using the recursion formular to calculateD(m), F(m), and Xi on node n, the connecting matrix Cn must be used as in the following steps:(
1
) the step of recursion caculation of Xi(41)Xi+1=GiXi+Hiδ(i=1,2,…,n-1),X-n=CnXn+F⃑P=Cn(Gn-1Xn-1+Hn-1δ),Xn+1=GnX-n+Hnδ.
The following process is the same as other nodes’:
(42)Xi+1=GiXi+Hiδ(i=n+1,n+2,…m).(
2
) The step of recursion caculation of D(m) and F(m).BecauseXn=D(n)X0+F(n), then we get
(43)Xn+1=GnX-n+Hnδ=Gn(CnXn+F→P)+Hnδ=Gn[Cn(D(n)X0+F(n))+F→P]+Hnδ=GnCnD(n)X0+GnCnF(n)+GnF⃑P+Hnδ=D(n+1)X0+F(n+1)⟹{D(n+1)=GnCnD(n)F(n+1)=GnCnF(n)+GnF⃑P+Hnδ.The recursion steps to calculateD(m) are as follows:
(44)D(i+1)=GiD(i)(i=1,2,…,n-1),D(n+1)=GnCnD(n),D(i+1)=GiD(i)(i=n+1,n+2,…m).The recursion steps to calculateF(m) is as follows:
(45)F(i+1)=GiF(i)+Hiδ(i=1,2,…,n-1),F(n+1)=GnCnF(n)+GnF⃑P+Hnδ,F(i+1)=GiF(i)+Hiδ(i=n+1,n+2,…m).


## 4.2.4. The Solution Flow Process

In short, the proposed solution procedure involves the following four main steps:(1)
calculating the value ofGn and Hn using the given (29),(2)
calculatingD(m) and F(m) using the given recursion formula of (33),(3)
calculating the vectorX0 by solving linear algebraic equation (34),(4)
calculatingXn using the given recursion formula of (27).Because the equations and solution formula are all given in form of matrix. 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.


## 5. Verification

The practical examples of a tied back retaining piles used in roadway cuts are considered herein to verify the developed numerical calculation techniques. Manual digging discrete reinforced concrete piles with prestressed tie-back anchors were used to stabilize a cut slope beside a highway roadbed at mountainous area, the city of ChongQing (shown in Figure7). 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 of slope stability is estimated to be 800 kN/m using the limit equilibrium method. This force is transferred by the pile to the stable layer below the slip surface. The piles were designed to be installed along the side of the road at a spacing of 5 m to increase the factor of safety of the whole slope to the required value of 1.3.Figure 7
A pile retaining wall with prestressed tie-back anchors (design case I).The stabilizing pile has a length ofL=30.0 m, sectional dimension 2.0 × 3.0 m. The height of the portion embedded into the sliding surface is 20.0 m. The pile is constructed using C30 concrete (assuming concrete does not crack). A lateral force F=800×5 kN is assumed to act upon the segment above the sliding surface. The pressure distribution is considered to have a rectangular shape, as proposed by Chinese design code.Two design schemes were presented for selection. For design case I (as shown in Figure7), the stabilizing pile is tied back at the top with one row of anchor. The anchor is prestressed to have an axial force of 550 kN. The unbonded length of the anchor is 28 m. Anchor inclination is 15 degrees from the horizontal. The anchors are prestressed at the time of installation. For design case II (as shown in Figure 8), in order to reduce the long-term relaxation of anchor prestress, the pile is tied back with two rows of anchors and the anchor’s axial preforce is designed to be 300 kN which is smaller than 550 kN. The unbonded length of the anchors is 28 m for the upper one and 25 m for the lower one. The upper anchor inclination is 15 degrees from the horizontal and the lower one’s inclination is 20 degrees from the horizontal.Figure 8
A pile retaining wall with prestressed tie-back anchors (design case II).The conceptual calculation model used to simulate the lateral response of the pile is shown in Figures9 and 10 respectively. The boundary condition at the pile base is considered as free head which allows both lateral displacement and rotation.Figure 9
The conceptual calculation model for stabilizing pile with one row of tie-back anchor (design case I).Figure 10
The conceptual calculation model for stabilizing pile with two rows of tie-back anchors (design case II).The modulus of horizontal subgrade reaction for the deposit and bedrock below the slip surface is assumed to follow the following distribution types, respectively.(1)
For the bedrock,k is assumed to be a constant of 100.0 × 106 Pa/m.(2)
For the deposit,k is assumed to vary linearly with depth as the following formular, k=A+mz. Where A is 30.0 × 106 Pa/m and m is 7.0 × 106 Pa/m2.For simplicity, the friction force along the pile-soil interface can be neglected.Then, the developed method is applied to analyze the lateral response of the anchored pile. In order to compare the results with those from the FEM, we built the corresponding finite element model which considers the pile as a series of connected elastic beam elements with its node supported by horizontal spring element representing the surrounding soil. Comparisons of shear force, bending moment, and deflection of the pile between boundary value method (BVM) and FEM are presented in Figures11 and 12. Complete agreement between them can be observed.Figure 11
Comparison of internal forces and deflection of pile between BVM and FEM (design case I).Figure 12
Comparison of internal forces and deflection of stabilizing pile between BVM and FEM (design case II).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 stabilizing pile with prestress tie-back anchors.


## 6. Summary and Conclusions

In this work, a new numerical uncoupled method for calculating the response of stabilizing piles with prestressed tie-back anchors is proposed. The detailed derivation of the governing differential equations and the proposed numerical solution scheme are presented. 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 anchored pile used for slope stabilization.



---

*Source: 601508-2013-09-24.xml*