Numerical Analysis of Laterally Loaded Long Piles in Cohesionless Soil

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


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, earthretaining 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 beamon-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 fullscale 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 = n j=0 c j P j (η). By applying the Galerkin method to minimize the residual yields n j=0 c j A P j (η) − f , P r (η) = 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 multilayered soil with the modulus of subgrade reaction varies with any functions with the depth.

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 E p and I p , 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   Assuming all forces acting horizontally and applying equilibrium conditions leads to According to that, Summation of moments about axis at end of element yields The relation between bending moment and the shear force is given by Differentiating Eq. (4) with respect to z and substituting it into Eq. (2) provide The basic moment-curvature relationship of elementary pile can be written as The differential equation of a laterally loaded pile can take the form where y(z) is the horizontal pile deflection at depth z and E p I p is the flexural stiffness of pile.
Assuming a linear variation of soil stiffness with soil depth, k can be written as follows [29] k = n h z, where here n h is the coefficient of subgrade reaction that represents the rate of increase of subgrade reaction modulus with depth z.
Eqs. (7) and (8) yield The solution of Eq. (9) requires an iterative procedure to achieve convergence of the relative stiffness factor T as follows [30] Defining a dimensionless variable x = z/T and introducing it into Eq. (9) yields By substituting Eq. (10) into Eq. (11), the governing equation for the soil-pile system is rewritten the form

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.
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. 2180 CMC, 2022, vol.71, no.2

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

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

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, P n (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].
with the orthogonality on the interval [−1, 1], i.e., Lemma [32]: If the arbitrary constants q, k, m and j are positive integers then and where Theorem: If the arbitrary constants k, l, q and m are positive integers, then where C q (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.

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 So, Eq. (12) could be rewritten as subject to the conditions of free-head pile Moreover, the conditions of fixed-head pile are: 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 Reducing Eq. (25) is performed by orthogonalizing the residual with respect to the basic functions as follows where the inner product · , · is defined as By substituting Eq. (28) into Eq. (29), a discrete system of n + 1 unknowns is given in the form c n−1 c n ] T and a j, r = P j (η), P r (η) + 2 L 3 (η + 1)P j (η), P r (η) . All values of a j, r could be calculated using the given theorem. The vector b based on the boundary conditions as follows.

Free-Head Piles
In this case, the vector b is written as The linear system in Eq. (31) By using the inverse linear transformation η = 2 L x−1 and x = z T , the solutions of the governing differential equation for a laterally loaded pile on elastic subgrade, Eq. (12) (35)

Fixed-Head Piles
In this case, the vector b is written as b = 0 0 0 · · ·  (38)
If the pile top condition is a free-head one, then the system of equations given in Eq. (39) is subjected to the initial conditions where the constants s 1 and s 2 are unknowns to be evaluated during the process of solution.
On the other hand, if the pile top condition is a fixed-head one, then the system of equations given in Eq. (39) is subjected to the initial conditions where the constants s 3 and s 4 are unknowns to be determined as a part of the numerical solution.
The system of equations Eq. (39) can be rewritten in the vector form as where ] and x max is used as an adaptation for x → ∞.
To obtain local error estimates for adaptive step-size control effectively, consider two Runge-Kutta formulas of different orders p andp = p + 1. A Runge-Kutta process generates a sequence Y n as an approximation of Y(x n ) for 0 = x 0 < x 1 < · · · < x n = x max . In the interval from x n to x n+1 = x n + h n , there are two approximations of Y(x n+1 ) called Y n+1 and Y n+1 for p andp = p + 1, respectively. Their forms are where s is the number of stages and, Considering p = 4 and s = 7, the coefficients a i,j , b i ,b i and c i = i−1 j=0 a i,j , i = 1, 2, · · · , s, can be evaluated by Bogacki and Shampine [24] to produce the efficient pair of formulas Y n+1 (4 th order formula) and Y n+1 (5 th order formula).
is generated based on the unknowns s 1 and s 2 for free-head pile and the unknowns s 3 and s 4 for fixed-head pile cases. The values of all these unknowns are determined upon solving the equations y(x max ) = y 1 (x max ) = y 2 (x max ) = y 3 (x max ) = 0 to satisfy the given boundary conditions given in Eqs. (14) and (15). The value of x max is adapted by gradual increase till the maximum difference between any values of s 1 , s 2 , s 3 or s 4 in succession is less than the required accuracy τ = 10 −6 .
The error between the two numerical solutions Y n+1 and Y n+1 is calculated as In case of e n+1 ≤ τ , one can use Y n+1 or Y n+1 as the final approximate value of Y(x n+1 ). On the other hand, if e n+1 > τ, the error e n+1 is used to adapt the step size h n toĥ n as followŝ The adapted step size is used to estimate the new values of Y n+1 and Y n+1 until achieving e n+1 ≤ τ .

Validation of the Proposed Methods
The performance and capability of the proposed methods to predict the behavior of laterally loaded piles in cohesion1ess soil have been demonstrated by comparing the obtained numerical results and the observed results from field experiments in full-scale lateral load tests reported by Cox. et al. [33]. In these tests, the flexible free-head steel tube pipe pile of 0.61 m in diameter, 21 m in length, 9.525 mm in wall thickness and 163000 kN·m 2 bending rigidity, was embedded in a deposit of submerged sand. The soil profile at the site is composed of a uniformly fine-graded sand with internal angle of friction, φ = 39 • and submerged density of γ =10.4 kN/m 3 . To investigate the soil properties below the ground surface, the Standard Penetration Test, SPT, was performed by [33] and reported as 18 blows per 30 cm. Accordingly, the recommended value of the horizontal subgrade reaction constant n h for submerged sand is estimated to be 15 MN/m 3 , [34]. Moreover, an analytical power series solution is also used for comparison, including the solution proposed by Fayun et al. [11], which reported a simplified analytical solution for laterally loaded long piles based on  The design of pile foundation under lateral loads is extensively bound to study both the pile head deflection and the maximum bending moment. Legendre-Galerkin method and Runge-Kutta formulas of order four and five were employed to solve the flexural equation of long piles embedded in homogeneous cohesionless soil with a modulus of subgrade reaction increases linearly with depth. In the numerical simulation, the pile is subjected to horizontal force H 0 = 100 kN at distance e = 0.305 m above-soil level, as shown in Fig. 1. The considered parameters in terms of pile length L, pile bending rigidity E P I P and coefficient of subgrade reaction n h , are set to be 21 m, 163000 kN/m 2 and 15000 kN/m 3 respectively. In order to ensure that our results produce reliable pile deflection and bending moment not only at the pile head but also for the entire pile length, and the analysis is equally applicable for different pile head conditions, the obtained results are compared with the results of simplified solution proposed by Fayun et al. [11]. The responses of free-head pile foundation in terms of lateral deflection and bending moment distributions are shown in Fig. 5. In addition, the corresponding responses in the case of fixed-head pile foundation are also presented in Fig. 6. The plotted curves reveal complete overlaps between Legendre-Galerkin solution and Runge-Kutta solution. Furthermore, the comparison charts demonstrate a very good agreement between the numerical results estimated via Legendre-Galerkin and Runge-Kutta and the corresponding procedure introduced by Fayun et al. [11]. These observations suggest that the lateral deflection and bending moment profiles can be represented accurately by the proposed method. However, insignificant differences between the location of peak deflection and bending moment from proposed methods and solution proposed by Fayun et al. can also be observed regardless the type of head condition. In the present study, Legendre-Galerkin and Runge-Kutta formulas of order four and five methods have been introduced to obtain simplified numerical approaches for understanding the behaviour of single piles against lateral loads. For the purpose of analysis and design of laterally loaded piles crossing sandy soil, simple expressions for the pile lateral deflection and bending moment can be evaluated. The procedure is programmed with the most computational software program Mathematica, which is considered as the world's leading computational software. The numerically computed pile responses are compared with the results from the full-scale lateral load tests. The proposed approaches are well validated. Moreover, these proposed approaches provide evidence that high precision can be achieved with a small amount of computational work. It has been found from the study that the Legendre-Galerkin solution almost coincides with the Runge-Kutta solution for both free-head and fixed-head piles. The suggested numerical expressions obtained in this study can be reasonably applied to analyze and design laterally loaded long piles conveniently. In addition, these techniques can also modify to design/analyze laterally loaded long piles in soil with the modulus of subgrade reaction in any functions of the depth. The proposed approaches capture the long pile behaviour. Furthermore, the proposed solutions are aimed at providing an effective and convenient method for engineers to predict the responses of the entire pile length under the applied lateral load.

Funding Statement:
The authors received no specific funding for this study.

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study.