Computers, Materials & Continua DOI:10.32604/cmc.2022.021899
Article

Numerical Analysis of Laterally Loaded Long Piles in Cohesionless Soil

Ayman Abd-Elhamed1,2,*, Mohamed Fathy3 and Khaled M. Abdelgaber1

1Physics and Engineering Mathematics Department, Faculty of Engineering-Mattaria, Helwan University, Cairo, 11718, Egypt
2Faculty of Engineering, King Salman International University, South Sinai, Ras Sedr, Egypt
3Basic and Applied Science Department, College of Engineering and Technology, Arab Academy for Science, Technology and Maritime Transport, Cairo, 11799, Egypt
*Corresponding Author: Ayman Abd-Elhamed. Email: moh_fathy79_6@aast.edu
Received: 19 July 2021; Accepted: 30 August 2021

Abstract: The capability of piles to withstand horizontal loads is a major design issue. The current research work aims to investigate numerically the responses of laterally loaded piles at working load employing the concept of a beam-on-Winkler-foundation model. The governing differential equation for a laterally loaded pile on elastic subgrade is derived. Based on Legendre-Galerkin method and Runge-Kutta formulas of order four and five, the flexural equation of long piles embedded in homogeneous sandy soils with modulus of subgrade reaction linearly variable with depth is solved for both free- and fixed-headed piles. Mathematica, as one of the world's leading computational software, was employed for the implementation of solutions. The proposed numerical techniques provide the responses for the entire pile length under the applied lateral load. The utilized numerical approaches are validated against experimental and analytical results of previously published works showing a more accurate estimation of the response of laterally loaded piles. Therefore, the proposed approaches can maintain both mathematical simplicity and comparable accuracy with the experimental results.

Keywords: Numerical solution; laterally loaded pile; cohesionless soil; Legendre-Galerkin; Runge-Kutta

1  Introduction

Pile foundations are frequently used, especially in weaker soils, to support various structures subjected to lateral loads such as high-rise buildings, communication towers, wind turbines, earth-retaining structures, bridges, tanks and offshore structures. Lateral loads owing to wind, wave, dredging, traffic and seismic events are considered significant on these structures since they are eventually transmitted to the piles [1,2]. As a result, the piles have been analyzed by considering a concentrated force and/or moment acting on the top of the pile. Over many decades, several methods have been proposed for designing and analysis of piles subjected to lateral loads including the subgrade reaction approach [3], the p-y approach [4,5], the finite element approach [6], the finite difference approach [7], and the analytical method [8]. Of these approaches, elastic solutions based on beam-on-Winkler-foundation model, albeit approximate, are probably the most widely used in engineering practice due to their simplicity, as well as they provide satisfactory results. Winkler's model is a particularly attractive approach used to reliably capture the soil-pile interaction. In this model the pile is simulated as a flexural beam connected to a series of narrowly spaced independent and continuous Winkler springs and dashpots distributed along the pile shaft. To investigate the mechanical behavior of laterally loaded piles in clay, several numerical investigations using the LPILE software were performed by Moayed et al. [9]. Moreover, Chang [10] derived an analytic solution to get the responses of laterally loaded long piles in cohesive soil considering constant subgrade reaction modulus. Furthermore, Different techniques have been proposed to analyze the behaviour of piles frequently subjected to lateral loads in sandy soil with different boundary conditions at their ends, including power series solution [11], finite element method [12], and finite difference method [3,13]. Further full-scale tests were performed to explore the behavior of laterally loaded piles either in sand [14] or in clay [15]. Numerical simulations considering both theoretical predictions and experimental validations are the common powerful tools of analysis in the field of geotechnical engineering particularly complex applications.

Recently, many researchers used Legendre polynomials in different methods to construct various mathematical models. These methods can solve Lane-Emden type of differential equation [16], differential equation with second and fourth order [17], the equation of Cahn-Hilliard [18], integral equation of Fredholm type [19], Helmholtz equation [20], Volterra integral equations in the second kind [21], integral-differential of Fredholm type in linear form [22] and Abels integral equation [23]. Based on Legendre-Galerkin method, the pile flexural equation can be written as Ay(η)=f, where A is a differential operator. The solution of our problem can be approximated in Legendre series as y=∑j=0ncjPj(η). By applying the Galerkin method to minimize the residual yields ∑j=0ncj⟨APj(η)−f,Pr(η)⟩=0. The differential equation is converted into discrete linear system. This system is solved, and the coefficients are determined for our approximate solution of the differential equation. In addition to the above method, the Bogacki–Shampine method [24] is a powerful numerical solution used for solving ordinary differential equations of the investigated problem. It is a Runge–Kutta method of two successive orders, for example 4th and 5th, with multi stages. An adaptive step size is implemented based on the error estimated between the solutions of the two successive orders. Elbashbeshy et al. [25,26] used this method in flow and heat transfer problems of fluid/nanofluid over different stretching surfaces and found that it is accurate by comparing it with methods used to solve the same problems.

This paper aims to present simple numerical methods to capture the behaviour of single piles under lateral working loads. First, mathematical formulation of the problem and the derivation of the governing differential equations are presented. Thereafter, two techniques are developed for numerical approximations of the derived equations. Then the effect of boundary conditions of the pile head on the behaviors of piles is studied as well. In order to examine the validity of the proposed numerical techniques, the obtained numerical results for deflection and bending moment along the laterally loaded piles are compared with results of previous studies. The proposed numerical techniques in the current study provide a better approach for structural designers to simply solve for the displacement and bending moment responses of laterally loaded piles. Consequently, the techniques can be easily applied in practice as an alternative approach to analyze and design laterally loaded long piles. It is worth noting that, the proposed solutions are applicable only for homogeneous soil with reaction modulus that varies linearly with depth. So, further investigation and analyses are needed for multi-layered soil with the modulus of subgrade reaction varies with any functions with the depth.

2  Mathematical Formulation

The model under examination consists of a pile that is assumed to be perfectly glued to the surrounding soils suggesting that there is no relative movement between the soil and the pile [27]. It is worth noting that the horizontal and vertical loads are rarely coupled. Thereupon, it is common among structural designers to neglect the vertical load effect during handling the laterally loaded pile problem. Accordingly, the pile under consideration is subjected to applied lateral force H0 at distance e above the soil surface, as shown in Fig. 1. Furthermore, the thickness of soil strata, L, coincides with the embedment length of pile. The pile elastic modulus and moment of inertia are Ep and Ip, respectively. The soil is characterized by a subgrade modulus k which linearly increases with the depth as a good approximation for granular soils [28]. For a laterally loaded pile on elastic homogeneous foundation, the governing differential equation employing the subgrade reaction theory can be derived by considering the equilibrium of transverse forces over a differential segment of the pile, as shown in Fig. 2.

Figure 1: Illustration of pile–soil system using Winkler model

Figure 2: Forces acting in differential pile element

Assuming all forces acting horizontally and applying equilibrium conditions leads to

V−(V+dVdzdz)−k(z)y(z)dz=0(1)

According to that,

dVdz+k(z)y(z)=0(2)

Summation of moments about axis at end of element yields

M−(M+dMdzdz)+Vdz=0(3)

The relation between bending moment and the shear force is given by

V=dMdz(4)

Differentiating Eq. (4) with respect to z and substituting it into Eq. (2) provide

d2Mdz2+k(z)y(z)=0(5)

The basic moment-curvature relationship of elementary pile can be written as

M=EpIpd2ydz2(6)

The differential equation of a laterally loaded pile can take the form

EpIpd4ydz4+k(z)y(z)=0,(7)

where y(z) is the horizontal pile deflection at depth z and EpIp is the flexural stiffness of pile.

Assuming a linear variation of soil stiffness with soil depth, k can be written as follows [29]

k=nhz,(8)

where here nh is the coefficient of subgrade reaction that represents the rate of increase of subgrade reaction modulus with depth z.

Eqs. (7) and (8) yield

EpIpd4ydz4+nhzy(z)=0.(9)

The solution of Eq. (9) requires an iterative procedure to achieve convergence of the relative stiffness factor T as follows [30]

T=(EpIpnh)0.2.(10)

Defining a dimensionless variable x=z/T and introducing it into Eq. (9) yields

EpIpnhT5d4ydx4+xy(x)=0.(11)

By substituting Eq. (10) into Eq. (11), the governing equation for the soil-pile system is rewritten the form

d4ydx4+xy(x)=0.(12)

3  Boundary Conditions

As shown in Fig. 3, the solution of Eq. (12) requires boundary conditions to balance the number of equations and the number of unknowns.

Figure 3: Boundary conditions for the pile under consideration

Boundary conditions of infinitely long piles fixed at the bottom can be written as

y(LT)=0,y′(LT)=0,y′′(LT)=0,y′′′(LT)=0.(13)

The boundary conditions at the top of the pile depend on the circumstances of the lateral deflection, slope, bending moment and shear force. These are generalized into the following two categories.

3.1 Free-Head Pile

The pile head is not restrained against rotation and translation. Therefore, the boundary conditions at the pile head can take the form

y′′′(0)=HOT3EpIp,y′′(0)=MOT2EpIp.(14)

3.2 Fixed-Head Pile

The pile head is completely restrained against rotation. The boundary conditions at the pile head are given by

y′′′(0)=HOT3EpIp,y′(0)=0.(15)

4  Legendre-Galerkin Method

4.1 Preliminaries

Orthogonal systems play a vital role in performing mathematical analysis. This can be due to functions belonging to very general classes can be expanded in series of orthogonal functions e.g., Fourier-Bessel series, Fourier series, etc. Orthogonal polynomials of degree n, Pn(x);n=0,1,2,…, represent an important class of orthogonal systems where many of the special functions encountered in the applications, such as Jacobi, Legendre, Laguerre, Hermite and Chebyshev polynomials, are part of that class. The Legendre polynomial is defined as [31].

Pn(x)=2n∑k=0n⁡(nk)(n+k−12n)xk,(16)

with the orthogonality on the interval [−1,1], i.e.,

∫−11⁡Pn(x)Pm(x)dx={22n+1,ifm=n0,ifm≠n.(17)

Lemma [32]: If the arbitrary constants q,k,m and j are positive integers then

Pk(q)(x)=∑i=0k−q⁡Cq(k,i)Pi(x),(k+1)iseven(18)

and

xmPj(x)=∑n=02m⁡Θm,n(j)Pj+m−n(x),P−r(x)=0forr≥1,(19)

where

Cq(k,i)=2q−1(2i+1)Γ(q+k−i2)Γ(q+k−i+12)Γ(q)Γ(2−q+k−i2)Γ(3−q+k−i2),(20)

Θm,n(j)=(−1)n2j+m−nm!(2j+2m−2n+1)Γ(j+1)Γ(j+m−n+1)×∑k=max(0,j−n)j+m−nΓ(j+k+1)2k(n+k−j)!Γ(3j+2m−2n−k+2),×∑ℓ=0j−k(−1)ℓΓ(2j+m−n−k−ℓ+1)Γ(j+m+ℓ−n+1)ℓ!(j−k−ℓ)!Γ(j−ℓ+1)Γ(k+ℓ+1)×2F1(j−k−n,j+m+ℓ;3j+2m−2n−k+2;2),(21)

and

Theorem: If the arbitrary constants k, l, q and m are positive integers, then

∫−11⁡Pk(q)(x)Pl(x)dx={2∑i=0k−qCq(k,i)2i+1,ifl=i,(k+1)iseven0,ifl≠i(22)

and

∫−11⁡xmPk(x)Pl(x)dx={2∑n=02mΘm,n(j)2(j+m−n)+1ifl=j+m−n0,ifl≠j+m−n,(23)

where Cq(k,i) and Θm,n(j) are defined in Eqs. (20) and (21), respectively.

Proof. By using the given lemma and Eq. (17), the theorem is proved.

4.2 Legendre-Galerkin Method

To solve Eq. (12), the solution domain of our problem must be changed from the interval [0,L] where L→∞ to the interval [−1,1] using the linear transformation

x=L2(η+1).(24)

So, Eq. (12) could be rewritten as

y′′′′(η)+(2L)3(η+1)y(η)=0,(25)

subject to the conditions of free-head pile

y′′′(−1)=(2L)3HT3EpIp,y′′(−1)=(2L)2MT2EpIp,y(1)=y′(1)=y′′(1)=y′′′(1)=0,(26)

Moreover, the conditions of fixed-head pile are:

y′′′(−1)=(2L)3HT3EpIp,y′′(−1)=y(1)=y′(1)=y′′(1)=y′′′(1)=0,(27)

The solution of Eq. (25) subject to the conditions given in Eqs. (26) or (27) is written in finite expansion of shifted Legendre function as

yn(η)≈∑j=0n⁡cjPj(η).(28)

Reducing Eq. (25) is performed by orthogonalizing the residual with respect to the basic functions as follows

⟨y′′′′(η),Pr(η)⟩+⟨(2L)3(η+1)y(η),Pr(η)⟩=0,(29)

where the inner product ⟨⋅,⋅⟩ is defined as

⟨ζ(η),χ(η)⟩=∫−11⁡ζ(η)χ(η)dη.(30)

By substituting Eq. (28) into Eq. (29), a discrete system of n+1 unknowns is given in the form

Ac=b,(31)

c=[c0c1c3⋯cn−1cn]T and aj,r=⟨Pj′′′′(η),Pr(η)⟩+⟨(2L)3(η+1)Pj(η),Pr(η)⟩. All values of aj,r could be calculated using the given theorem. The vector b based on the boundary conditions as follows.

4.2.1 Free-Head Piles

In this case, the vector b is written as

b=[000⋯(2L)3HT3EpIp(2L)2MT2EpIp0000]T.(32)

The linear system in Eq. (31) is solved using one of the available numerical techniques. By obtaining the coefficients values of {cj}0n, the approximate solution given in Eq. (28) becomes

y(η)=3.99HEI2/5nh3/5(−1.9×10−6+8.02×10−5η−8.4×10−4η2+2.9×10−3η3+9.6×10−3η4−0.109η5+0.284η6−0.012η7−1.13η8+2.06η9−3.17η10+5.04η11+6.35η12−35.29η13+24.47η14+52.42η15−74.02η16−24.81η17+82.27η18−9.96η19−47.01η20+16.04η21+13.97η22−6.69η23−1.73η24+η25)+2.58MEI2/5nh3/5(9.44×10−7+3.61−5η−8.31−4η2+6.121−3η3−0.010η4−0.079η5+0.376η6−0.373η7−0.591η8+1.02η9−2.26η10+11.16η11−9.75η12−33.95η13+59.65η14+23.97η15−103.32η16+18.82η17+88.61η18−40.16η19−41.67η20+26.40η21+10.35η22−8.13η23−1.06η24+η25).(33)

By using the inverse linear transformation η=(2L)x−1andx=zT, the solutions of the governing differential equation for a laterally loaded pile on elastic subgrade, Eq. (12), are given for deflection and moment as

y(z)=2.429Hnh3/5(EI5−0.67EI24/5znh1/5+0.07EI22/5z3nh3/5−3.05×10−7EI21/5z4nh4/5−8.33×10−3EI4z5nh+1.84×10−3EI19/5z6nh6/5+1.7×10−5EI18/5z7nh7/5−5.7×10−5EI17/5z8nh8/5+1.04×10−5EI16/5z9nh9/5−2.99×10−6EI3z10nh2+1.28×10−6EI14/5z11nh11/5−3.73×10−7EI13/5z12nh12/5+7.28×10−8EI12/5z13nh13/5−1.02×10−8EI11/5z14nh14/5+1.07×10−9EI2z15nh15/5−8.76×10−11EI9/5z16nh16/5+5.62×10−12EI8/5z17nh17/5−2.84×10−13EI7/5z18nh18/5+1.13×10−14EI6/5z19nh19/5−3.52×10−16EIz20nh20/5+8.35×10−18EI4/5z21nh21/5−1.46×10−19EI3/5z22nh22/5+1.79×10−21EI2/5z23nh23/5−1.36×10−23EI1/5z24nh24/5+4.85×10−26z25nh5)+1.619Mnh3/5(EI5−1.08nh1/5EI24/5z+0.31nh2/5EI23/5z2+1.1×10−6nh4/5EI21/5z4−8.34×10−3nhEI4z5+3.01×10−3nh6/5EI19/5z6−3.8×10−4nh7/5EI18/5z7+1.2×10−5nh8/5EI17/5z8−4.8×10−6nh9/5EI16/5z9+2.6×10−6nh2EI3z10−3.97×10−7nh11/5EI14/5z11−2.81×10−8nh12/5EI13/5z12+2.15×10−8nh13/5EI12/5z13−3.97×10−7nh11/5EI14/5z11−2.81×10−8nh12/5EI13/5z12+2.15×10−8nh13/5EI12/5z13−4.47×10−9nh14/5EI11/5z14+5.82×10−10nh3EI2z15−5.44×10−11nh16/5EI9/5z16+3.84×10−12nh17/5EI8/5z17−2.09×10−13nh18/5EI7/5z18+8.86×10−15nh19/5EI6/5z19−2.89×10−16nh4EIz20+7.17×10−18nh21/5EI4/5z21−1.31×10−19nh22/5EI3/5z22+1.65×10−21nh23/5EI2/5z23−1.29×10−23nh24/5EI1/5z24+4.71×10−26nh5z25)(34)

and

M(z)=H(EI27/5z−8.8×10−6nh1/5EI26/5z2−0.41nh2/5EI5z3+0.134nh3/5EI24/5z4+1.8×10−3nh4/5EI23/5z5−7.8×10−3nhEI22/5z6+1.8×10−3nh6/5EI21/5z7−6.5×10−4nh7/5EI4z8+3.4×10−4nh8/5EI19/5z9−1.19×10−4nh9/5EI18/5z10+2.7×10−5nh2EI17/5z11−4.5×10−6nh11/5EI16/5z12+5.48×10−7nh12/5EI3z13−5.11×10−8nh13/5EI14/5z14+3.71×10−9nh14/5EI13/5z15−2.11×10−10nh3EI12/5z16+9.41×10−12nh16/5EI11/5z17−3.24×10−13nh17/5EI2z18+8.52×10−15nh18/5EI9/5z19−1.64×10−16nh19/5EI8/5z20+2.20×10−18nh4EI7/5z21−1.82×10−20nh21/5EI6/5z22+7.07×10−23nh22/5EIz23)+Mnh1/5(EI28/5+2.3×10−5nh2/5EI26/5z2−0.27nh3/5EI5z3+0.15nh4/5EI24/5z4−2.64×10−2nhEI23/5z5+1.1×10−3nh6/5EI22/5z6−5.6×10−4nh7/5EI21/5z7+3.8×10−4nh8/5EI4z8−7×10−5nh9/5EI19/5z9−6.01×10−6nh2EI18/5z10+5.4×10−6nh11/5EI17/5z11−1.3×10−6nh12/5EI16/5z12+1.978×10−7nh13/5EI3z13−2.11×10−8nh14/5EI14/5z14+1.69×10−9nh3EI13/5z15−1.038×10−10nh16/5EI12/5z16+4.91×10−12nh17/5EI11/5z17−1.781×10−13nh18/5EI2z18+4.878×10−15nh19/5EI9/5z19−9.76×10−17nh4EI8/5z20+1.349×10−18nh21/5EI7/5z21−1.151×10−20nh22/5EI6/5z22+4.57×10−23nh23/5EIz23).(35)

4.2.2 Fixed-Head Piles

In this case, the vector b is written as b=[000⋯(2L)3HT3EpIp00000]T. The approximate solution given in Eq. (28) becomes

y(η)=1.599HEI2/5nh3/5(−6.15×10−6+1.45×10−4η−8.656×10−4η2−1.82×10−3η3+3.95×10−2η4−0.153η5+0.147η6+0.527η7−1.94η8+3.62η9−4.53η10−4.10η11+30.45η12−37.298η13−28.146η14+94.987η15−30.193η16−90.098η17+72.802η18+35.205η19−54.997η20+0.553η21+19.385η22−4.528η23−2.727+η25).(36)

By using the inverse linear transformation η=(2L)x−1andx=zT, the solutions of the governing differential equation for a laterally loaded pile on elastic subgrade, Eq. (12), are given for deflection and moment as

y(z)=0.928Hnh3/5(EI5−0.50EI23/5z2nh2/5+0.179EI22/5z3nh3/5−2.72×10−6EI21/5z4nh4/5−8.308×10−3EI4z5nh−6.2×10−5EI19/5z6nh6/5+6.75×10−4EI18/5z7nh7/5−1.7×10−4EI17/5z8nh8/5+3.5×10−5EI16/5z9nh9/5−1.21×10−5EI3z10nh2+4×10−6EI14/5z11nh11/5−9.336×10−7EI13/5z12nh12/5+1.557×10−7EI12/5z13nh13/5−1.944×10−8EI11/5z14nh14/5+1.870×10−9EI2z15nh3−1.414×10−10EI9/5z16nh16/5+8.484×10−12EI8/5z17nh17/5−4.048×10−13EI7/5z18nh18/5+1.529×10−14EI6/5z19nh19/5−4.526×10−16EIz20nh4+1.026×10−17EI4/5z21nh21/5−1.725×10−19EI3/5z22nh22/5+2.025×10−21EI2/5z23nh23/5−1.482×10−23EI1/5z24nh24/5+5.089×10−26z25nh5)(37)