# Boundary Value Problem for Analysis of Portal Double-Row Stabilizing Piles

Authors: Cheng Huang
Journal: Journal of Applied Mathematics
 (2013)

---

## 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.


---

## Body

## 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 Figure1), 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 1
The uncoupled analysis model for the double-row portal stabilizing piles.


### 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 (p) acting on the pile can be modeled by discrete independent linear or nonlinear springs in the form of the following equation:
(1)p=k·w,
where k 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. p is the horizontal soil reaction pressure (it has a unit of force/length2). w 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 Figure2, having a 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,Q) in the original pile may become external forces on the isolated free portion.Figure 2
The free-body diagram of infinitesimal isolated segment of pile.Figure2 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 Figure2 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 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, and the positive φ is counterclockwise. The lateral pressure qn is considered positive when applied from left to right. The lateral pressure qτ 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 Figure2), we arrive at two force equilibrium equations in the directions of u and v and one moment equilibrium equation:
(2)(N+dNdsds)-N-qτds-ksHuDds=0,(Q+dQdsds)-Q-qnds-knHvbds=0,(M+dMdsds)-M+(Q+dQdsds)ds2+Qds2=0.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:(3)dNds=ksHDu+qτ,dQds=knHbv+qn,dMds=-Q,
where H={1,contact,0,noncontact, the “one” indicates that Winkler reaction force exists and the “zero” indicates that Winkler reaction force does not exist. S is the position coordinate; 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(4)X=[NQM];Z=[uvφ];P=[qτqn0].The set of equations of equilibrium (3) can be rewritten in the following matrix form:
(5)dXds=B·X+L·Z+P,
where B=[0000000-10]; L=[ksHD000knHb0000].


### 2.4. Geometric and Constitutive Equations

When the deformation (du-, dv-, dφ-) of the differential element (shown in Figure 2) induced by the internal forces (N,Q,M) is considered given, the corresponding strains can be expressed as
(6)dU-ds=(du-ds,dv-ds,dφ-ds).
According to the related theory of elastic beam, the internal forces (N,Q,M) can be related to strains as in the following linear constitutive equation:
(7)dU-ds=[du-dsdv-dsdφ-ds]=[1EA000αGA0001EI]·[NQM],
where α is the 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; EI is the flexural rigidity of the pile’s cross-section.Considering the deformation,dU- 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 displacement onto this direction which takes the form B-Zds, where B- is the a undetermined third-order square matrix. Then, the deformation dU- can be expressed as
(8)dU-=dZ+B-Zds.
Applying the principle of virtual work to the isolated differential element of pile (shown in Figure 2), B- can be determined. We suppose that each point of the body is given an infinitesimal virtual displacement δZ satisfying displacement boundary conditions where 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 (5) and (8) 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 which means that B-=BT. At last, we develop the following geometric and constitutive equations:
(9)duds=NEA,dvds=αQGA+φ,dφds=MEI.


### 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:
(10)dNds=ksHDu+qτ,dQds=knHbv+qn,dMds=-Q,duds=NEA,dvds=αQGA+φ,dφds=MEJ.(11)LetK=[000HksD000000Hknb00-100001EA000000αGA0001001EI000],X={NQMuvφ}T,p={qτqn0000}T.Then, we can use matrix notation to present the above equations (10) in the form:
(12)dXds=KX+p=F(X,s),
where K is the coefficients square matrix of six order and X and p 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).


## 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 (p) acting on the pile can be modeled by discrete independent linear or nonlinear springs in the form of the following equation:
(1)p=k·w,
where k 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. p is the horizontal soil reaction pressure (it has a unit of force/length2). w 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 Figure2, having a 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,Q) in the original pile may become external forces on the isolated free portion.Figure 2
The free-body diagram of infinitesimal isolated segment of pile.Figure2 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 Figure2 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 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, and the positive φ is counterclockwise. The lateral pressure qn is considered positive when applied from left to right. The lateral pressure qτ 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 Figure2), we arrive at two force equilibrium equations in the directions of u and v and one moment equilibrium equation:
(2)(N+dNdsds)-N-qτds-ksHuDds=0,(Q+dQdsds)-Q-qnds-knHvbds=0,(M+dMdsds)-M+(Q+dQdsds)ds2+Qds2=0.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:(3)dNds=ksHDu+qτ,dQds=knHbv+qn,dMds=-Q,
where H={1,contact,0,noncontact, the “one” indicates that Winkler reaction force exists and the “zero” indicates that Winkler reaction force does not exist. S is the position coordinate; 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(4)X=[NQM];Z=[uvφ];P=[qτqn0].The set of equations of equilibrium (3) can be rewritten in the following matrix form:
(5)dXds=B·X+L·Z+P,
where B=[0000000-10]; L=[ksHD000knHb0000].


## 2.4. Geometric and Constitutive Equations

When the deformation (du-, dv-, dφ-) of the differential element (shown in Figure 2) induced by the internal forces (N,Q,M) is considered given, the corresponding strains can be expressed as
(6)dU-ds=(du-ds,dv-ds,dφ-ds).
According to the related theory of elastic beam, the internal forces (N,Q,M) can be related to strains as in the following linear constitutive equation:
(7)dU-ds=[du-dsdv-dsdφ-ds]=[1EA000αGA0001EI]·[NQM],
where α is the 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; EI is the flexural rigidity of the pile’s cross-section.Considering the deformation,dU- 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 displacement onto this direction which takes the form B-Zds, where B- is the a undetermined third-order square matrix. Then, the deformation dU- can be expressed as
(8)dU-=dZ+B-Zds.
Applying the principle of virtual work to the isolated differential element of pile (shown in Figure 2), B- can be determined. We suppose that each point of the body is given an infinitesimal virtual displacement δZ satisfying displacement boundary conditions where 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 (5) and (8) 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 which means that B-=BT. At last, we develop the following geometric and constitutive equations:
(9)duds=NEA,dvds=αQGA+φ,dφds=MEI.


## 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:
(10)dNds=ksHDu+qτ,dQds=knHbv+qn,dMds=-Q,duds=NEA,dvds=αQGA+φ,dφds=MEJ.(11)LetK=[000HksD000000Hknb00-100001EA000000αGA0001001EI000],X={NQMuvφ}T,p={qτqn0000}T.Then, we can use matrix notation to present the above equations (10) in the form:
(12)dXds=KX+p=F(X,s),
where K is the coefficients square matrix of six order and X and p 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 (N,Q,M,u,v,φ). 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:
(13)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 calculation (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 matrices. In this study, the following possible pile end conditions were considered.(1)
Free head (allows both displacement and rotation):N=Q=0, M=0. The corresponding matrix of boundary condition is
(14)[100000010000001000]×{NQMuvφ}=o.(2)
In case of bottom end hinged (allows rotation without displacement):u=v=0, M=0. The corresponding matrix of boundary condition is
(15)[001000000100000010]×{NQMuvφ}=o.(3)
In case of bottom end fixed (allows neither displacement nor rotation):u=v=0, φ=0. The corresponding matrix of boundary condition is
(16)[000100000010000001]×{NQMuvφ}=o.(4)
In case of bottom end partially hinged (allows rotation without vertical displacement):u=0, Q=0, and M=0. The corresponding matrix of boundary condition is
(17)[010000001000000100]×{NQMuvφ}=o.(5)
Elastic vertical support at the bottom end:N=Kv·u·Area, Q=0, and M=0. The corresponding matrix of boundary condition is
(18)[100-Kv·Area00010000001000]×{NQMuvφ}=o,where Kv 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:
(19)dXds=KX+p,CX|S=0=o,DX|S=L=o.
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 (N,Q,M) are so much higher than those of the displacements (u,v,φ) 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 K. We multiply the displacements (u,v,φ) by E and substitute the original displacement variables by the expressions (Eu,Ev,Eφ). So, we redefine two variables as follows:
(20)X~={NQMu-v-φ-}={NQMEuEvEφ},K~=[000HksDE000000HknbE00-100001A000000EαGA0001001I000].Finally, we obtain the following system of ordinary differential equations:(21)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).


### 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 formdX/ds=F(X,s). 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:
(22)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 (22) is a numerical method of fifth order, where δ is the stepsize of difference and F(X,s)=K-X+p. For convenience of computation, this formulation may be rewritten in the following form:
(23)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 represents the beginning point of calculation and n=m represents the end point of calculation.In the above equation,Gn and Hn can be obtained from the following recursion formula:
(24)Kn(j)=K~(sn+δj),Pn(j)=p(sn+δj),α1=α2=α3=12,α4=1,β1=β4=16,β2=β3=13,δj=(∑k=1jγk)δ,γ1=γ3=0,γ2=γ4=12,δj⟹δ1=0,δ2=δ3=12δ,δ4=δ,(25)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,
where I is the identity matrix.


#### 5.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 (23) and imposing the boundary conditions at pile head and base.Considering the recursion formula of (23), Xn can be expressed in terms of X0 as follows:
(26)Xn=D(n)X0+F(n).In the case ofn=0, we have D(0)=I,F(0)=0 and substitute it into the recursion formula of (23). We get
(27)Xn+1=Gn(D(n)X0+F(n))+Hnδ.It can be rewritten as follows:(28)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:
(29)D(0)=I,D(n+1)=GnD(n),F(0)=0,F(n+1)=GnF(n)+Hnδ.Now considering the case of boundary point:Xm=D(m)X0+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:
(30)[CDD(m)]X0=[o-DF(m)].The above set 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 (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 K- and p 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 Dm, Fm, and  Xn, recursively.Figure 3
Scheme of local coordinate system transformation between pile and connection beam.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(31)[N-Q-]=[cosθ-sinθsinθcosθ][NQ].
And the bending moment remains unchanged.According to the principle of vector analysis, we get(32)[u-v-]=[cosθ-sinθsinθcosθ][uv].
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:
(33)[N-nQ-nM-nu-nv-nφ-n]=[cosβ-sinβ0000sinβcosβ0000001000000cosβ-sinβ0000sinβcosβ0000001]×[NnQnMnunvnφn],X-n=Cn·Xn.As we can see, the above connection matrix is the so-called orthogonal matrix whose inverse is its transposed matrix.As shown in Figure3, the node of difference (no. n) is the connection node. First we need to distinguish K- and p in the recursion formula for the two connected segments and then insert the connection matrix for transformation from Xn to X-n when dealing with the Dm, Fm, and Xn, recursively. The detailed procedure is as follows.(1)
For the calculation ofXn,
(34)Xi+1=GiXi+Hiδ,(i=1,2,…,n-1)X-n=CnXn=Cn(Gn-1Xn-1+Hn-1δ), where Cn is the transformation matrix for point n. Thus, Xn+1=GnX-n+Hnδ. Next,
(35)Xi+1=GiXi+Hiδ,(i=n+1,n+2,…m),
where m indicates the end point of calculation.(2)
For the calculation ofDm: because Xn=D(n)X0+F(n), and
(36)Xn+1=GnX-n+HnδXn+1=Gn(CnXn)+HnδXn+1=Gn(Cn(D(n)X0+F(n)))+HnδXn+1=GnCnD(n)X0+GnCnF(n)+HnδXn+1=D(n+1)X0+F(n+1)Xn+1⟹{D(n+1)=GnCnD(n)F(n+1)=GnCnF(n)+Hnδ,
So, the above procedure also applies to calculation of D(m) as follows:
(37)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).(3)
For the calculation ofFm,
(38)F(i+1)=GiF(i)+Hiδ,(i=1,2,…,n-1),F(n+1)=GnCnF(n)+Hnδ,F(i+1)=GiF(i)+Hiδ,(i=n+1,n+2,…m),where m 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 ofGn and Hn using the given Equation (25);(2)
calculatingD(m) and F(m) using the given recursion formula of (29);(3)
calculating the vectorX0 by solving linear algebraic Equations (30);(4)
calculatingXn 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.


## 5.1. Uniformity Preprocessing

The orders of magnitude of the section internal forces (N,Q,M) are so much higher than those of the displacements (u,v,φ) 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 K. We multiply the displacements (u,v,φ) by E and substitute the original displacement variables by the expressions (Eu,Ev,Eφ). So, we redefine two variables as follows:
(20)X~={NQMu-v-φ-}={NQMEuEvEφ},K~=[000HksDE000000HknbE00-100001A000000EαGA0001001I000].Finally, we obtain the following system of ordinary differential equations:(21)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).


## 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 formdX/ds=F(X,s). 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:
(22)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 (22) is a numerical method of fifth order, where δ is the stepsize of difference and F(X,s)=K-X+p. For convenience of computation, this formulation may be rewritten in the following form:
(23)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 represents the beginning point of calculation and n=m represents the end point of calculation.In the above equation,Gn and Hn can be obtained from the following recursion formula:
(24)Kn(j)=K~(sn+δj),Pn(j)=p(sn+δj),α1=α2=α3=12,α4=1,β1=β4=16,β2=β3=13,δj=(∑k=1jγk)δ,γ1=γ3=0,γ2=γ4=12,δj⟹δ1=0,δ2=δ3=12δ,δ4=δ,(25)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,
where I is the identity matrix.


### 5.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 (23) and imposing the boundary conditions at pile head and base.Considering the recursion formula of (23), Xn can be expressed in terms of X0 as follows:
(26)Xn=D(n)X0+F(n).In the case ofn=0, we have D(0)=I,F(0)=0 and substitute it into the recursion formula of (23). We get
(27)Xn+1=Gn(D(n)X0+F(n))+Hnδ.It can be rewritten as follows:(28)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:
(29)D(0)=I,D(n+1)=GnD(n),F(0)=0,F(n+1)=GnF(n)+Hnδ.Now considering the case of boundary point:Xm=D(m)X0+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:
(30)[CDD(m)]X0=[o-DF(m)].The above set 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 (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 K- and p 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 Dm, Fm, and  Xn, recursively.Figure 3
Scheme of local coordinate system transformation between pile and connection beam.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(31)[N-Q-]=[cosθ-sinθsinθcosθ][NQ].
And the bending moment remains unchanged.According to the principle of vector analysis, we get(32)[u-v-]=[cosθ-sinθsinθcosθ][uv].
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:
(33)[N-nQ-nM-nu-nv-nφ-n]=[cosβ-sinβ0000sinβcosβ0000001000000cosβ-sinβ0000sinβcosβ0000001]×[NnQnMnunvnφn],X-n=Cn·Xn.As we can see, the above connection matrix is the so-called orthogonal matrix whose inverse is its transposed matrix.As shown in Figure3, the node of difference (no. n) is the connection node. First we need to distinguish K- and p in the recursion formula for the two connected segments and then insert the connection matrix for transformation from Xn to X-n when dealing with the Dm, Fm, and Xn, recursively. The detailed procedure is as follows.(1)
For the calculation ofXn,
(34)Xi+1=GiXi+Hiδ,(i=1,2,…,n-1)X-n=CnXn=Cn(Gn-1Xn-1+Hn-1δ), where Cn is the transformation matrix for point n. Thus, Xn+1=GnX-n+Hnδ. Next,
(35)Xi+1=GiXi+Hiδ,(i=n+1,n+2,…m),
where m indicates the end point of calculation.(2)
For the calculation ofDm: because Xn=D(n)X0+F(n), and
(36)Xn+1=GnX-n+HnδXn+1=Gn(CnXn)+HnδXn+1=Gn(Cn(D(n)X0+F(n)))+HnδXn+1=GnCnD(n)X0+GnCnF(n)+HnδXn+1=D(n+1)X0+F(n+1)Xn+1⟹{D(n+1)=GnCnD(n)F(n+1)=GnCnF(n)+Hnδ,
So, the above procedure also applies to calculation of D(m) as follows:
(37)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).(3)
For the calculation ofFm,
(38)F(i+1)=GiF(i)+Hiδ,(i=1,2,…,n-1),F(n+1)=GnCnF(n)+Hnδ,F(i+1)=GiF(i)+Hiδ,(i=n+1,n+2,…m),where m 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 ofGn and Hn using the given Equation (25);(2)
calculatingD(m) and F(m) using the given recursion formula of (29);(3)
calculating the vectorX0 by solving linear algebraic Equations (30);(4)
calculatingXn 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.


## 5.2.1. Derivation of the Recursion Formula

The following finite difference formula (22) is one format of the Runge-Kutta methods:
(22)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 (22) is a numerical method of fifth order, where δ is the stepsize of difference and F(X,s)=K-X+p. For convenience of computation, this formulation may be rewritten in the following form:
(23)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 represents the beginning point of calculation and n=m represents the end point of calculation.In the above equation,Gn and Hn can be obtained from the following recursion formula:
(24)Kn(j)=K~(sn+δj),Pn(j)=p(sn+δj),α1=α2=α3=12,α4=1,β1=β4=16,β2=β3=13,δj=(∑k=1jγk)δ,γ1=γ3=0,γ2=γ4=12,δj⟹δ1=0,δ2=δ3=12δ,δ4=δ,(25)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,
where I is the identity matrix.


## 5.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 (23) and imposing the boundary conditions at pile head and base.Considering the recursion formula of (23), Xn can be expressed in terms of X0 as follows:
(26)Xn=D(n)X0+F(n).In the case ofn=0, we have D(0)=I,F(0)=0 and substitute it into the recursion formula of (23). We get
(27)Xn+1=Gn(D(n)X0+F(n))+Hnδ.It can be rewritten as follows:(28)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:
(29)D(0)=I,D(n+1)=GnD(n),F(0)=0,F(n+1)=GnF(n)+Hnδ.Now considering the case of boundary point:Xm=D(m)X0+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:
(30)[CDD(m)]X0=[o-DF(m)].The above set 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 (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 K- and p 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 Dm, Fm, and  Xn, recursively.Figure 3
Scheme of local coordinate system transformation between pile and connection beam.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(31)[N-Q-]=[cosθ-sinθsinθcosθ][NQ].
And the bending moment remains unchanged.According to the principle of vector analysis, we get(32)[u-v-]=[cosθ-sinθsinθcosθ][uv].
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:
(33)[N-nQ-nM-nu-nv-nφ-n]=[cosβ-sinβ0000sinβcosβ0000001000000cosβ-sinβ0000sinβcosβ0000001]×[NnQnMnunvnφn],X-n=Cn·Xn.As we can see, the above connection matrix is the so-called orthogonal matrix whose inverse is its transposed matrix.As shown in Figure3, the node of difference (no. n) is the connection node. First we need to distinguish K- and p in the recursion formula for the two connected segments and then insert the connection matrix for transformation from Xn to X-n when dealing with the Dm, Fm, and Xn, recursively. The detailed procedure is as follows.(1)
For the calculation ofXn,
(34)Xi+1=GiXi+Hiδ,(i=1,2,…,n-1)X-n=CnXn=Cn(Gn-1Xn-1+Hn-1δ), where Cn is the transformation matrix for point n. Thus, Xn+1=GnX-n+Hnδ. Next,
(35)Xi+1=GiXi+Hiδ,(i=n+1,n+2,…m),
where m indicates the end point of calculation.(2)
For the calculation ofDm: because Xn=D(n)X0+F(n), and
(36)Xn+1=GnX-n+HnδXn+1=Gn(CnXn)+HnδXn+1=Gn(Cn(D(n)X0+F(n)))+HnδXn+1=GnCnD(n)X0+GnCnF(n)+HnδXn+1=D(n+1)X0+F(n+1)Xn+1⟹{D(n+1)=GnCnD(n)F(n+1)=GnCnF(n)+Hnδ,
So, the above procedure also applies to calculation of D(m) as follows:
(37)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).(3)
For the calculation ofFm,
(38)F(i+1)=GiF(i)+Hiδ,(i=1,2,…,n-1),F(n+1)=GnCnF(n)+Hnδ,F(i+1)=GiF(i)+Hiδ,(i=n+1,n+2,…m),where m 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 ofGn and Hn using the given Equation (25);(2)
calculatingD(m) and F(m) using the given recursion formula of (29);(3)
calculating the vectorX0 by solving linear algebraic Equations (30);(4)
calculatingXn 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 Figure4) 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 c 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.Figure 4
Scheme of portal double-row piles used to stabilize a landslide (units: m).The front row pile has a lengthL=16.0 m and sectional dimension 1.5×2.0m. The length of the portion embedded into the sliding surface is 6.0 m. The back row pile has a length L=37.0 m and sectional dimension 1.5×2.8m. 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 F=2147×4kN 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 Figure5. The boundary condition at the pile base is considered as free head which allows both lateral displacement and rotation.Figure 5
The conceptual calculation model for portal double-row stabilizing pile.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 Table1 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.Table 1
Parameters of material properties.
Young's modulus(GPa)
Shear modulus(GPa)
Modulus of subgrade reaction(MPa/m)
C30 concrete
30
12
/
Completely decomposed mudstone
/
/
80
Intensely weathered mudstone
/
/
110
Moderately weathered mudstone
/
/
150Then, 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 Figures6 and 7. Complete agreement between them can be observed.Figure 6
Comparison of internal force and deflection of front row pile between the boundary value method and FEM.Figure 7
Comparison of internal force and deflection of back row pile between the boundary value method and FEM.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.



---

*Source: 485632-2013-09-09.xml*