Energy-Based Approach: Analysis of a Laterally Loaded Pile in Multi-Layered Non-Linear Elastic Soil Strata

: Several studies have been reported in published literature on analytical solutions for a laterally loaded pile installed in a homogeneous single soil layer. However, piles are rarely installed in an ideal homogeneous single soil layer. The present study describes a new continuum-based analysis or energy-based approach for predicting the pile displacement responses subjected to static lateral loads and moments considering the soil non-linearity. This analytical analysis treats the pile as an elastic Euler–Bernoulli beam and the soil as a three-dimensional (3D) continuum in which the non-linear elastic properties are described by a modulus degradation relationship. The principle of virtual work was applied to the energy equation of a pile–soil system in order to obtain the governing differential equation for the pile and soil displacements. An iterative procedure was adopted to solve the equations numerically using a ﬁnite difference method (FDM). The pile displacement response was obtained using the software MATLAB R2021a, and the results from the energy-based method were compared with those obtained from the ﬁeld test data as well as the ﬁnite element analysis (FEA) based on the software ANSYS Workbench 2021R1. The present study investigated the effect of explicit incorporation of soil properties and layering through a parametric study in order to understand the importance of predicting appropriate pile displacement responses in a linear elastic soil system. The responses indicated that the effect of soil layers and their thicknesses, pile properties and the variation in soil moduli have a direct impact on the displacements of piles subjected to lateral loading. Hence, a proper emphasis has to be given to account for the soil non-linearity. Considering the effect of soil non-linearity, it is observed that the results obtained from the energy-based method agreed well with the ﬁeld measured values and those obtained from the FEA. The results indicated a difference of approximately less than 7% between the proposed method and the FEA. The approach presented in this study can be further extended to piles embedded in multi-layered soil strata subjected to the combined action of axial loads, lateral loads and moments. Furthermore, the same approach can be extended to study the response of the soil to group piles.


Introduction
The growing importance to analyze the structures, such as high-rise buildings, bridges, offshore platforms, etc., resting on pile foundations and acted upon by various horizonal forces (wind, wave, currents and seismic events) has led to various analysis methods over time.An extensive literature review has been conducted by the authors Moussa and Christou [1] who summarized and grouped the various analysis methods of laterally loaded single pile under static loading into four categories: (a) ultimate limit state method (ULS); (b) subgrade reaction approach; (c) finite element method (FEM); and (d) continuum method.Several researchers developed different types of the ultimate limit state (ULS) methods, such as the Blum method [2], Brinch Hansen's method [3] and Brom's method [4][5][6], which aim at obtaining the maximum lateral loads that are supported by the piles.However, the major limitation for the ULS methods is their unrealistic approach to the pile-soil interaction problem, as the mathematical formulations of such methods consider only the pile lateral displacements (no soil deformations).The different types of methods in the subgrade reaction approach developed and studied by several researchers over time include Winkler's approach [7][8][9][10][11][12][13][14][15][16][17][18], the P-y curve method [19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38] and the strain wedge (SW) models [39].These methods have been used to calculate the allowable lateral pile deflection assuming that the soil is represented by a series of linear or non-linear springs.Although the soil idealization in Winkler's hypothesis based on linear independent springs overestimates the soil continuity, shear coupling between the springs and the effect of strength characteristics of the piles on the subgrade reaction, the method is still used because of its simplicity and acceptable results, but it is not recommended to be used when the soil profile is highly non-linear.The p-y curves, on the other hand, are back calculated from empirical test results; therefore, they are highly dependent on the empirical test environment and the variation in soil and pile properties.Hence, it is critical to select the most appropriate p-y curve in the analysis of laterally loaded piles to obtain accurate and realistic results.Furthermore, a review for the applicability of the selected p-y curve is crucial, since they require both the p and y values from the field.The ability to consider the 3D nature of the problem, the soil continuity as well as the soil-pile interaction makes the use of the SW method very appealing to the practicing engineers.Conversely, the two main limitations of this method are: (a) few empirical data were used to develop the stress-strain relationship, and (b) the determination of strain wedge depth and the value of the subgrade reaction modulus below the strain wedge is not a simple task.The finite element method (FEM) is a numerical technique that is based on the concept of the continuum approach developed by various researchers [40][41][42][43][44][45][46][47][48][49][50][51][52][53] to represent the complexity of the loading and the soil boundary conditions.The advancement in the computational power of high-speed computers facilitated the use of FEM in the analysis of laterally loaded piles.However, the limitation lies in its time-consuming nature and also points to the need to validate the finite element model prior to using the results for the design.Finally, the continuum approach is used to obtain the allowable lateral displacement by idealization of the soil as a linear/non-linear elastic infinite stratum.The boundary element method (BEM) and simplified continuum models are the two different mediums used to solve such models.The BEM was implemented in several studies by Refs.[54][55][56][57][58][59][60][61][62][63].The simplified continuum models were further categorized as (i) energy and variational calculus method [64][65][66][67][68][69][70][71] and (ii) modified Reissner models [72][73][74][75][76][77][78].In general, the continuum methods involve complex mathematical formulations to model the continuity of the soil, its non-linearity and the boundary conditions, and require the need to determine a suitable soil modulus.In addition, the complexity of a 3D continuum often requires the use of numerical methods, such as the boundary integral/element method, the FEM and the finite difference method (FDM) [68].
Most studies on laterally loaded piles do not account for layering within the soil medium [79,80].In the p-y method of analysis, soil heterogeneity with depth is approximately considered by assuming different p-y curves with depth [13].Variation of soil properties with depth has been approximately considered in the continuum approach as well by assuming (typically) a linear variation of soil modulus with depth [58,62,63,66].Most solutions that have been used to predict the deformation of a pile under different loads treat the soil as a linear elastic material in order to simplify the calculations.However, in reality, the behavior of soil is highly non-linear, and soil stiffness depends on the stress and strain levels.Therefore, it is highly important to consider the effect of soil non-linearity in multi-layered soil strata.
The present study is based on an energy-based approach (simplified continuum model) to study the displacement profile of a laterally loaded pile in multi-layered soil strata considering soil non-linear elastic behavior.The pile head is subjected to a static force and/or moment with the pile being free at the head and fixed at the tip.Soil moduli (λ and G) are not treated as constants but assumed to vary in radial, circumference and depth directions according to the strain and stress levels.A simple power law from the published literature is adopted, which accounts for the soil non-linear behavior by the degradation of the soil modulus (G) with strain.The pile displacement response is obtained using the software MATLAB R2021a.The present study is based on the FDM, since it yields results faster when compared to the FEM.The results from the energy-based method are compared with those from the field test data as well as the FEA based on the software ANSYS Workbench 2021R1.The results from the energy-based method are in good agreement with both the field test data and the FEA.The present approach is a continuation of the work performed on axially loaded piles by Arvan and Arockiasamy [81].

Problem Definition
A pile subjected to lateral loading can be modeled as a Euler-Bernoulli beam, especially for long and slender piles in which the pile shear deformation can be neglected for large slenderness ratios of L/D > 10 [82].A laterally loaded pile in an isotropic non-linear elastic multi-layered soil medium is shown in Figure 1.This study considers a pile of length L with circular cross section of radius r p .The pile is embedded in n horizontal soil layers, and the pile head is subjected to a horizontal force F 0 accompanied with a bending moment M 0 .The horizontal soil layers extend to infinity in the radial direction, and the bottom n th layer extends to infinity in the vertical direction.
Geotechnics 2022, 2, FOR PEER REVIEW 3 force and/or moment with the pile being free at the head and fixed at the tip.Soil moduli (λ and G) are not treated as constants but assumed to vary in radial, circumference and depth directions according to the strain and stress levels.A simple power law from the published literature is adopted, which accounts for the soil non-linear behavior by the degradation of the soil modulus (G) with strain.The pile displacement response is obtained using the software MATLAB R2021a.The present study is based on the FDM, since it yields results faster when compared to the FEM.The results from the energy-based method are compared with those from the field test data as well as the FEA based on the software ANSYS Workbench 2021R1.The results from the energy-based method are in good agreement with both the field test data and the FEA.The present approach is a continuation of the work performed on axially loaded piles by Arvan and Arockiasamy [81].

Problem Definition
A pile subjected to lateral loading can be modeled as a Euler-Bernoulli beam, especially for long and slender piles in which the pile shear deformation can be neglected for large slenderness ratios of L/D > 10 [82].A laterally loaded pile in an isotropic non-linear elastic multi-layered soil medium is shown in Figure 1.This study considers a pile of length L with circular cross section of radius rp.The pile is embedded in n horizontal soil layers, and the pile head is subjected to a horizontal force F0 accompanied with a bending moment M0.The horizontal soil layers extend to infinity in the radial direction, and the bottom n th layer extends to infinity in the vertical direction.The terms H1, H2, H3… Hn−1 denote the vertical height from the ground surface to the bottom of any layer i.Therefore, the thickness of any layer i is Hi-Hi−1 with H0 = 0. Due to the axisymmetric problem behavior, a system of cylindrical coordinates (r-θ-z) is chosen, with the origin coinciding with the center of the pile cross section at the pile head, and the z axis coinciding with the pile axis.The pile head is considered to be free, and the tip of the pile is clamped.Another important assumption to be noted is that there is no slippage The terms H 1 , H 2 , H 3 . . .H n−1 denote the vertical height from the ground surface to the bottom of any layer i.Therefore, the thickness of any layer i is H i -H i−1 with H 0 = 0. Due to the axisymmetric problem behavior, a system of cylindrical coordinates (r-θ-z) is chosen, with the origin coinciding with the center of the pile cross section at the pile head, and the z axis coinciding with the pile axis.The pile head is considered to be free, and the tip of the pile is clamped.Another important assumption to be noted is that there is no slippage or separation between the pile and the surrounding soil and between soil layers.The stresses and the displacement within a soil continuum are shown below in Figure 2a,b.
Geotechnics 2022, 2, FOR PEER REVIEW 4 the axisymmetric problem behavior, a system of cylindrical coordinates (r-θ-z) is chosen, with the origin coinciding with the center of the pile cross section at the pile head, and the z axis coinciding with the pile axis.The pile head is considered to be free, and the tip of the pile is clamped.Another important assumption to be noted is that there is no slippage or separation between the pile and the surrounding soil and between soil layers.The stresses and the displacement within a soil continuum are shown below in Figure 2a,b.At each point in the soil domain, G and λ were calculated (Figure 3).

Soil Non-Linearity
Soil behaves as a linear elastic material at an extremely low range of strain, and the shear modulus starts degrading at a strain as low as 10 −5 .Since the soil particles constantly change their position during the application of a load, the resistance offered by the soil mass against deformation also changes; this results in the change in the value of soil modulus with an increase in strain.A non-linear stress-strain curve can be incorporated in an elastic analysis by properly estimating the secant modulus for a given level of strain At each point in the soil domain, G and λ were calculated (Figure 3).At each point in the soil domain, G and λ were calculated (Figure 3).

Soil Non-Linearity
Soil behaves as a linear elastic material at an extremely low range of strain, and the shear modulus starts degrading at a strain as low as 10 −5 .Since the soil particles constantly change their position during the application of a load, the resistance offered by the soil mass against deformation also changes; this results in the change in the value of soil modulus with an increase in strain.A non-linear stress-strain curve can be incorporated in an elastic analysis by properly estimating the secant modulus for a given level of strain (or stress).Soil non-linearity is taken into the account by estimating the decay of soil stiffness with strain [71,82,83].

Soil Non-Linearity
Soil behaves as a linear elastic material at an extremely low range of strain, and the shear modulus starts degrading at a strain as low as 10 −5 .Since the soil particles constantly change their position during the application of a load, the resistance offered by the soil mass against deformation also changes; this results in the change in the value of soil modulus with an increase in strain.A non-linear stress-strain curve can be incorporated in an elastic analysis by properly estimating the secant modulus for a given level of strain (or stress).Soil non-linearity is taken into the account by estimating the decay of soil stiffness with strain [71,82,83].

574
The variation of the shear stress with strain can be described using two parameters, A and n, which were obtained experimentally using a pressuremeter test, as shown in the equation below: where q represents the equivalent shear stress, ε q is the deviator shear strain.Atkinson [84] shows that the decay curves of soil stiffness with strain can be divided into three regions, as shown in Figure 4.The first region in Figure 4 represents the very small strain where the stiffness G 0 is constant; the second region comprising small strains starts from ε 0 until ε = 0.1%, and the third region exceeds ε = 0.1%, indicative of large strains.
Geotechnics 2022, 2, FOR PEER REVIEW 5 The variation of the shear stress with strain can be described using two parameters, A and n, which were obtained experimentally using a pressuremeter test, as shown in the equation below: where q represents the equivalent shear stress, εq is the deviator shear strain.Atkinson [84] shows that the decay curves of soil stiffness with strain can be divided into three regions, as shown in Figure 4.The first region in Figure 4 represents the very small strain where the stiffness G0 is constant; the second region comprising small strains starts from ε0 until ε = 0.1%, and the third region exceeds ε = 0.1%, indicative of large strains.In the second region, the stiffness decays rapidly, and in the third region, with large strain levels, the stiffness is the smallest, which concludes that soil stiffness is high at small strain and decreases with large strain [84].Figure 5 shows the degradation of soil stiffness with increasing strains for different clay types.In the second region, the stiffness decays rapidly, and in the third region, with large strain levels, the stiffness is the smallest, which concludes that soil stiffness is high at small strain and decreases with large strain [84].Figure 5 shows the degradation of soil stiffness with increasing strains for different clay types.
Geotechnics 2022, 2, FOR PEER REVIEW 5 The variation of the shear stress with strain can be described using two parameters, A and n, which were obtained experimentally using a pressuremeter test, as shown in the equation below: where q represents the equivalent shear stress, εq is the deviator shear strain.Atkinson [84] shows that the decay curves of soil stiffness with strain can be divided into three regions, as shown in Figure 4.The first region in Figure 4 represents the very small strain where the stiffness G0 is constant; the second region comprising small strains starts from ε0 until ε = 0.1%, and the third region exceeds ε = 0.1%, indicative of large strains.In the second region, the stiffness decays rapidly, and in the third region, with large strain levels, the stiffness is the smallest, which concludes that soil stiffness is high at small strain and decreases with large strain [84].Figure 5 shows the degradation of soil stiffness with increasing strains for different clay types.The present study is based on a non-linear elastic model developed by Osman et al. [86], which assumes the decay of soil stiffness with strain using a power law to describe the stress-strain behavior of soil [87,88]: where a = G 0 ε n q0 is a constant determined empirically; n describes soil non-linearity, which is equal to (−0.5), according to the experimental data analyzed by Osman et al. [86] (Figure 6).
zr represents the deviatoric strain; and ε q0 is the maximum deviatoric strain with linear elastic behavior, which is equal to 10 −5 .Soil stiffness G is estimated by calculating the strain at each location followed by the power law.
Geotechnics 2022, 2, FOR PEER REVIEW 6 The present study is based on a non-linear elastic model developed by Osman et al. [86], which assumes the decay of soil stiffness with strain using a power law to describe the stress-strain behavior of soil [87,88]: where  = is a constant determined empirically; n describes soil non-linearity, which is equal to (−0.5), according to the experimental data analyzed by Osman et al. [86] (Figure 6). = [( −  ) + ( −  ) + ( −  ) ] + ( +  +  ) represents the deviatoric strain; and 0 is the maximum deviatoric strain with linear elastic behavior, which is equal to 10 −5 .Soil stiffness G is estimated by calculating the strain at each location followed by the power law.

Basic Assumptions
For a laterally loaded pile, the displacement at any point within the continuum (Figure 2) can be expressed as a product of three separable functions, each accounting for one of the dimensions [64,89].
where u(z) is a displacement function (with a dimension of length) varying with depth z, representing the deflection of the pile axis;   ,   and   are the soil displacements in

Basic Assumptions
For a laterally loaded pile, the displacement at any point within the continuum (Figure 2) can be expressed as a product of three separable functions, each accounting for one of the dimensions [64,89].
where u(z) is a displacement function (with a dimension of length) varying with depth z, representing the deflection of the pile axis; u r , u θ and u z are the soil displacements in the direction r, θ and z; φ r (r) and φ θ (r) are dimensionless soil displacement functions varying with the radial coordinate r, and θ is measured from a vertical reference section (r = r 0 ) that contains the applied force vector F 0 (Figure 2).The above functions are set as φ r (r) = 1 and φ θ (r) = 1 at r = r 0 (i.e., at the pile-soil interface) and φ r (r) = 0 and φ θ (r) = 0 at r = ∞.Such an assumption ensures the decrease in the displacement (owing to pile deflection) within the soil mass with increasing radial distance from the pile.Thus, φ r and φ θ vary between 1 at the pile-soil interface to 0 at an infinite radial distance from the pile [82,83,90,91].It is also reasonable to assume that the vertical displacement of the pile owing to the lateral load and moment is negligible, and this justifies Equation (4c).

Governing Differential Equation
In the derivation of the governing differential equations that can capture the nonlinear soil response, the soil within each layer is assumed to be elastic and isotropic but heterogeneous (with respect to r and θ but not with respect to z), with no sliding or separation between the soil layers or between the pile and the soil.By solving the equations (valid for elastic, heterogeneous soil) for different magnitudes of load (with appropriate values and variations of soil modulus), the analysis can trace the non-linear progression of pile deflection (due to soil non-linearity) with increasing applied load.
The total potential energy of the pile-soil system, including both the internal and external potential energies, is given by where u is the lateral pile deflection, and σ ij , ε ij are the stress and strain tensors in the soil (Figure 2).The first integral represents the internal potential energy of the pile.The second and third integrals represent the internal potential energy of the continuum.The remaining two terms represent the external potential energy.Soil non-linearity is considered by varying the soil elastic parameters (G and λ) at each discretized nodal point in the soil domain.The stress-strain and strain-displacement relationships at any given nodal point in the soil medium are idealized by the following relationships.The stress-strain relationship is expressed as where G and λ are the elastic constants of the soil.The strain-displacement relationship is given by Geotechnics 2022, 2

577
Combining Equations ( 6) and (7) gives the strain energy density within any layer: The equilibrium equations are obtained using the principle of minimum potential energy, according to which δΠ = 0. Substituting Equation (8) into Equation ( 5) and applying δΠ = 0, we obtain The variations δu, δφ r and δφ θ are functions of u(z), φ r (r) and φ θ (r) and are independent.To obtain the pile governing equation, from Equation ( 9), all terms associated with δu, δ(du/dz) will be collected.Furthermore, the terms that are related to δφ r and δφ θ will be collected to obtain the governing equation of the soil.

Output Parameters (i) Pile Displacement
The governing equation of piles can be obtained by collecting the terms associated with δu and δ(du/dz) for (0 ≤ z ≤ L) from Equation ( 9) and then equating the summation to zero.
The governing equation becomes where where i denotes the number of the layer from the surface to the length of the pile and where j represents the layer number from pile length to infinity The fourth-order differential Equation ( 11) can be solved using a central finite difference scheme, represented as where i denotes the ith node in z direction, and ∆z is the distance between two nodes.At each point in the soil domain, λ and G were calculated (Figure 3).
(ii) Soil Displacement Similar to the pile governing equation, the soil governing equation can be formulated by considering the variations of φ r (r) and φ θ (r).Collecting the terms associated with φ r and φ θ over the domain r 0 ≤ r ≤∞ and equating the summation to zero yields Simplifying Equations ( 16) and ( 17 where n s1 , n s2 , m s1 , m s2 and m s3 for homogeneous and non-homogeneous soil comprising multi-layer soil summations with i representing the number of layers and n being the last layer number are given below Equations ( 18) and ( 19) can be simplified further using integration by parts of the terms that consist of δ((dφ r )/dr) and δ((dφ θ )/dr), respectively.The governing equation is then obtained by dividing the resulting equations of ( 18) by (−m s1 ) and ( 19) by (−m s2 ).
Considering the soil boundary conditions φ r = 1 at r = r 0 ; φ r = 0 at r → ∞; φ θ = 1 at r = r 0 ; and φ θ = 0 at r → ∞, the governing differential equations of soil can be rewritten as where All the equations mentioned in this study are solved using the software MATLAB 2021a.

Iterative Solution Methodology
The pile displacement Equation ( 11) can be solved when the soil parameters k and C are known.However, these parameters are dependent on the unknown dimensionless soil functions, φ r and φ θ , which can be estimated by calculating y 1 , y 2 , y 3 , y 4 , y 5 , y 6 , y 7 and y 8 .Soil deformation is calculated in radial and circumference directions, then calculated in the depth direction.To obtain the soil displacement, the initial numbers of these values, y 1 , y 2 , y 3 , y 4 , y 5 , y 6 , y 7 and y 8 , must be assumed.Then, they are inserted into Equations ( 23) and ( 24), from which the soil parameters k and C are obtained as a result of the pile displacement.New values of y 1 , y 2 , y 3 , y 4 , y 5 , y 6 , y 7 and y 8 can then be inserted into Equations ( 23) and ( 24) to evaluate φ r and φ θ , which are then inserted into Equation ( 11) to obtain the pile displacement u, so an iteration technique is needed to satisfy the condition γ old −γ new γ old < 0.001.Figure 7 shows the flow chart of the iterative solution procedure.
Geotechnics 2022, 2, FOR PEER REVIEW (23) and (24) to evaluate ϕr and ϕθ, which are then inserted into Equation ( 11) to the pile displacement u, so an iteration technique is needed to satisfy the con < 0.001.Figure 7 shows the flow chart of the iterative solution procedure

Effect of Explicit Incorporation of Soil Characteristics and Layering
The present study investigated the effect of explicit incorporation of soil pro and layering in order to understand the importance of predicting appropriate pi placement responses in a linear elastic soil system.Two different types of piles wi ying soil properties and layers were considered.Table 1 shows the pile and soil considered for the analyses.

Effect of Explicit Incorporation of Soil Characteristics and Layering
The present study investigated the effect of explicit incorporation of soil properties and layering in order to understand the importance of predicting appropriate pile displacement responses in a linear elastic soil system.Two different types of piles with varying soil properties and layers were considered.Table 1 shows the pile and soil details considered for the analyses.Figures 8 and 9 show the pile displacement responses along the depth z for all the cases of the two different pile types discussed above (Table 1).9 show the pile displacement responses along the depth z for all the cases of the two different pile types discussed above (Table 1).The responses indicated that the effect of soil layers and their thicknesses, pile properties and the variation in soil moduli have a direct impact on the displacements of piles subjected to lateral loading.Hence, a proper emphasis has to be given to account for soil non-linearity.The present study is capable enough in predicting accurate pile responses considering proper characterization of soil deposits and explicit accounting for different layers, as it utilizes 3D interactions between the pile and the soil deposits considering the effect of soil non-linearity.The responses indicated that the effect of soil layers and their thicknesses, pile properties and the variation in soil moduli have a direct impact on the displacements of piles subjected to lateral loading.Hence, a proper emphasis has to be given to account for soil non-linearity.The present study is capable enough in predicting accurate pile responses considering proper characterization of soil deposits and explicit accounting for different layers, as it utilizes 3D interactions between the pile and the soil deposits considering the effect of soil non-linearity.

Accuracy of the Model Considering Soil Non-Linearity (a) Energy-based Method versus FEA
The present analytical study is compared with the results of a 3D FEA in order to verify its accuracy and computational efficiency.The FEA analyses were performed using the software ANSYS Workbench 2021R1, which substantially reproduced the realistic behavior of the pile-soil system.ANSYS SOLID65 eight-node elements (three translational degrees of freedom) were used for both the pile and soil.Appropriate separation was set up between both the bodies in order to ensure that there is no slippage/separation.The soil-pile interface was set manually using the contact tool in which the pile was the "target" surface and the soil a "contact" surface.The soil medium was extended to about 30 times the pile diameter in both the horizontal and radial directions.The top surface of the soil medium was at the same level with the pile head and extended to a finite direction.The boundary conditions for the pile were set to being free at the head and fixed at the tip in order to replicate the response obtained by the energy-based assumptions.For the soil, the model included roller supports on the sides and a fixed support at the bottom.Figure 10 shows the initial input properties for the pile and soil and their boundary conditions.
times the pile diameter in both the horizontal and radial directions.The top surface of the soil medium was at the same level with the pile head and extended to a finite direction The boundary conditions for the pile were set to being free at the head and fixed at the tip in order to replicate the response obtained by the energy-based assumptions.For the soil the model included roller supports on the sides and a fixed support at the bottom.Figure 10 shows the initial input properties for the pile and soil and their boundary conditions.A concentrated force and/or moment was applied to the pile head that varied in several fixed increments.The soil non-linear elastic relationship is defined using Equation ( 2).The same curve generated for the variation of the secant shear modulus of soil with strain (Equation ( 2) is given as an input in the material definition section of the FEA A concentrated force and/or moment was applied to the pile head that varied in several fixed increments.The soil non-linear elastic relationship is defined using Equation (2).The same curve generated for the variation of the secant shear modulus of soil with strain (Equation ( 2) is given as an input in the material definition section of the FEA software.A directional deformation function was used to record the load-displacement responses for the different load levels.The discretized finite element model consisted of a linear explicit coarse mesh with 5908 nodes and 3855 elements.A mesh representation of the pile-soil model is shown in Figure 11.
software.A directional deformation function was used to record the load-displacement responses for the different load levels.The discretized finite element model consisted of a linear explicit coarse mesh with 5908 nodes and 3855 elements.A mesh representation of the pile-soil model is shown in Figure 11. Figure 12a,b show the pile head displacements obtained for the applied forces and moments, respectively.The results were compared for the soil conditions modeled as a linear elastic and a non-linear elastic material.The results obtained by the present analytical study (energy-based method) agree reasonably well with the FEA with a difference of approximately less than 7%.For an applied lateral force at the pile head, the non-linear elastic response shows a slight underestimation until an applied force of about 400 kN. Figure 12a,b show the pile head displacements obtained for the applied forces and moments, respectively.The results were compared for the soil conditions modeled as a linear elastic and a non-linear elastic material.The results obtained by the present analytical study (energy-based method) agree reasonably well with the FEA with a difference of approximately less than 7%.For an applied lateral force at the pile head, the non-linear elastic response shows a slight underestimation until an applied force of about 400 kN.
Another advantage of the present study is its computational efficiency over the FEA.Table 2 depicts the computational speed of the present energy-based approach solved using MATLAB Scripts and the FEA using ANSYS 2021R1 for the problems solved with a computer Intel Core i7-9750H CPU @ 2.60 GHz and 16 GB RAM.It is also important to note that the time taken to design a specific model in the FEA software is not defined in the table, which adds an additional advantage for the developed ready-to-use MATLAB Scripts.The pile-soil problem analyzed above is now compared against different values for Poisson's ratio to the 3D FEA to confirm the accuracy of the present study, considering the non-linear elastic constitutive relationship given by Equation (2). Figure 13a,b show the results obtained for pile head lateral displacements for two different values of Poisson's ratio (v = 0.2 and v = 0.49).Such values of Poisson's ratios were chosen, as they are the most adopted values of the previous studies [82,83,[90][91][92].It is evident from the figures that the present study produces responses, which are in good agreement with the FEA.
Figure 12a,b show the pile head displacements obtained for the applied forces and moments, respectively.The results were compared for the soil conditions modeled as a linear elastic and a non-linear elastic material.The results obtained by the present analytical study (energy-based method) agree reasonably well with the FEA with a difference of approximately less than 7%.For an applied lateral force at the pile head, the non-linear elastic response shows a slight underestimation until an applied force of about 400 kN.Another advantage of the present study is its computational efficiency over the FEA.Table 2 depicts the computational speed of the present energy-based approach solved using MATLAB Scripts and the FEA using ANSYS 2021R1 for the problems solved with a computer Intel Core i7-9750H CPU @ 2.60 GHz and 16 GB RAM.It is also important to note that the time taken to design a specific model in the FEA software is not defined in the table, which adds an additional advantage for the developed ready-to-use MATLAB Scripts.The pile-soil problem analyzed above is now compared against different values for Poisson's ratio to the 3D FEA to confirm the accuracy of the present study, considering the non-linear elastic constitutive relationship given by Equation (2). Figure 13a,b show the results obtained for pile head lateral displacements for two different values of Poisson's ratio (v = 0.2 and v = 0.49).Such values of Poisson's ratios were chosen, as they are the most adopted values of the previous studies [82,83,[90][91][92].It is evident from the figures that the present study produces responses, which are in good agreement with the FEA.  Figure 14a-g show the soil and pile responses obtained using the present energybased method and the FEA. Figure 14a,b show the responses of a laterally loaded pile in terms of pile head displacements with respect to the depth and the magnitude of applied forces.The input soil and pile properties are shown in the figures themselves.The results are compared to the 3D FEA and indicate that the present energy-based method is capable of producing responses quite close to the FEA. Figure 14c-g show the soil responses for the analyzed laterally loaded pile.Figure 14c-e represent the normalized radial (ur/rp) and tangential (uθ/rp) displacements along θ = 0° (Figure 14c), θ = 45° (Figure 14d) and θ = 90° (Figure 14e).Figure 14f shows the response of the dimensionless displacement functions ϕr and ϕθ along θ = 45° to further demonstrate the accuracy of the present energy-based method with the FEA. Figure 14g shows the response of the variation of the secant shear modulus (G/G0) along the radial direction for θ = 0°, 45° and 90°. Figure 15a-g show the pile and soil responses for the same problem analyzed but with an applied moment at the pile head.The present study produces both the soil and pile responses accurately when compared with the FEA. Figure 14a-g show the soil and pile responses obtained using the present energy-based method and the FEA. Figure 14a,b show the responses of a laterally loaded pile in terms of pile head displacements with respect to the depth and the magnitude of applied forces.The input soil and pile properties are shown in the figures themselves.The results are compared to the 3D FEA and indicate that the present energy-based method is capable of producing responses quite close to the FEA. Figure 14c-g show the soil responses for the analyzed laterally loaded pile.Figure 14c-e represent the normalized radial (u r /r p ) and tangential (u θ /r p ) displacements along θ = 0 • (Figure 14c), θ = 45 • (Figure 14d) and θ = 90 • (Figure 14e).Figure 14f shows the response of the dimensionless displacement functions φ r and φ θ along θ = 45 • to further demonstrate the accuracy of the present energy-based method with the FEA. Figure 14g shows the response of the variation of the secant shear modulus (G/G 0 ) along the radial direction for θ = 0 • , 45 • and 90 • .Figure 15a-g show the pile and soil responses for the same problem analyzed but with an applied moment at the pile head.The present study produces both the soil and pile responses accurately when compared with the FEA.
(Figure 14e).Figure 14f shows the response of the dimensionless displacement functions ϕr and ϕθ along θ = 45° to further demonstrate the accuracy of the present energy-based method with the FEA. Figure 14g shows the response of the variation of the secant shear modulus (G/G0) along the radial direction for θ = 0°, 45° and 90°. Figure 15a-g show the pile and soil responses for the same problem analyzed but with an applied moment at the pile head.The present study produces both the soil and pile responses accurately when compared with the FEA."A Numerical Study into Lateral Cyclic Nonlinear Soil-Pile Response": Allotey and ElNaggar [92].
A numerical analysis of lateral cyclic non-linear soil-pile response was carried out by Allotey and El Naggar [92] using a beam on non-linear Winkler foundation model (BNWF dynamic model), where the model was compression dominant, used to model the pilesoil interaction.This model was produced to predict the response of pile for four main (b) Energy-based method versus the published numerical analysis and FEA "A Numerical Study into Lateral Cyclic Nonlinear Soil-Pile Response": Allotey and ElNaggar [92].
A numerical analysis of lateral cyclic non-linear soil-pile response was carried out by Allotey and El Naggar [92] using a beam on non-linear Winkler foundation model (BNWF dynamic model), where the model was compression dominant, used to model the pile-soil interaction.This model was produced to predict the response of pile for four main parts: the backbone curve, the standard reload curve, the general unload curve and the direct reload curve.Allotey and El Naggar [92] simulated a reinforced concrete pile of length 12 m and diameter 0.6 m embedded in a single stratum consisting of uniform medium-stiff clay with an undrained shear strength of 50 KPa and unit weight of 19 kN/m 2 .The pile head was laterally loaded with 25 uniform two-way loading cycles with an amplitude of 213 kN and a pile-yield moment of 636 kNm.This study is compared to Case C of the paper, which assumed full connection between the pile and the soil.A FEA was also conducted using the software ANSYS Workbench 2021R1 for the soil-pile problem, as discussed above.A static structural analysis was chosen, and the input data for the pile characteristics, soil characteristics and soil non-linearity were taken from Allotey and El Naggar [92].The soil block's length was chosen as 25 times the pile diameter (25D), which was fixed at the bottom and free on the surroundings.The discretized FEM had a fine mesh of 9194 nodes and 1845 elements.The directional deformation along the x axis was used to calculate the displacements of pile head for the respective force applied.
Figure 16 shows the comparison of load-deflection responses obtained from the published numerical analyses, energy-based solution and the FEA.The pile head was laterally loaded with 25 uniform two-way loading cycles with an amplitude of 213 kN and a pile-yield moment of 636 kNm.This study is compared to Case C of the paper, which assumed full connection between the pile and the soil.A FEA was also conducted using the software ANSYS Workbench 2021R1 for the soil-pile problem, as discussed above.A static structural analysis was chosen, and the input data for the pile characteristics, soil characteristics and soil non-linearity were taken from Allotey and El Naggar [92].The soil block's length was chosen as 25 times the pile diameter (25D), which was fixed at the bottom and free on the surroundings.The discretized FEM had a fine mesh of 9194 nodes and 1845 elements.The directional deformation along the x axis was used to calculate the displacements of pile head for the respective force applied.
Figure 16 shows the comparison of load-deflection responses obtained from the published numerical analyses, energy-based solution and the FEA.The present study makes a comparison of the proposed energy-based method with the pile load test data obtained by McClelland and Focht [19].This study further compares the results to those proposed by several other researchers [63,64,90].The pile was embedded in a normally consolidated clay, and the pile details included a pile length of 23 m and a radius of 0.305 m.The pile was subjected to both lateral force and moment of 300 kN and −265 kNm, respectively.The elastic modulus of the pile was given by Ref. [63] as 68.42 × 106 kN/m 2 .The details of the soil strata are shown in Table 3.The present study makes a comparison of the proposed energy-based method with the pile load test data obtained by McClelland and Focht [19].This study further compares the results to those proposed by several other researchers [63,64,90].The pile was embedded in a normally consolidated clay, and the pile details included a pile length of 23 m and a radius of 0.305 m.The pile was subjected to both lateral force and moment of 300 kN and −265 kNm, respectively.The elastic modulus of the pile was given by Ref. [63] as 68.42 × 106 kN/m 2 .The details of the soil strata are shown in Table 3.The results of the present study show the pile deflection profile closely matches the field data, as shown in Figure 17.The results of the present study show the pile deflection profile closely matches the field data, as shown in Figure 17.The present study considers soil non-linearity, which is more representative of the actual soil behavior, substantiated by the available published field data.The pile deflections considering soil non-linearity appear to be within the bounds of linear elastic soil analyses.

Discussion and Conclusions
The beam-on-non-linear-foundation approach (p-y method) has been used extensively to understand the pile-displacement responses.This approach has the limitation of representing the surrounding soil using non-linear springs requiring realistic user-specified soil characteristics of the p-y curves, which do not represent the three-dimensional pile-soil interaction.Several researchers have studied the pile responses using numerical methods, including BEM, FEM and FDM.Although the FEM using appropriate soil constitutive relationship, elements and domains for the soil and pile gives realistic results, the method requires enormous computation time, and the resources required for such an analysis stand out as a major limitation.This gives rise to the need for a more efficient method, which has both the strength of a rigorous three-dimensional non-linear approach for the pile-soil interaction and potentials for obtaining a faster computational effort.The present study utilizes an analytical solution based on an energy method (discretized The present study considers soil non-linearity, which is more representative of the actual soil behavior, substantiated by the available published field data.The pile deflections considering soil non-linearity appear to be within the bounds of linear elastic soil analyses.

Discussion and Conclusions
The beam-on-non-linear-foundation approach (p-y method) has been used extensively to understand the pile-displacement responses.This approach has the limitation of representing the surrounding soil using non-linear springs requiring realistic user-specified soil characteristics of the p-y curves, which do not represent the three-dimensional pile-soil interaction.Several researchers have studied the pile responses using numerical methods, including BEM, FEM and FDM.Although the FEM using appropriate soil constitutive relationship, elements and domains for the soil and pile gives realistic results, the method requires enormous computation time, and the resources required for such an analysis stand out as a major limitation.This gives rise to the need for a more efficient method, which has both the strength of a rigorous three-dimensional non-linear approach for the pile-soil interaction and potentials for obtaining a faster computational effort.The present study utilizes an analytical solution based on an energy method (discretized continuum approach) to predict pile-soil displacements, where the soil is assumed to be non-linear elastic (soil parameters vary in radial, tangential and vertical directions).The analysis was conducted on a laterally loaded pile embedded in multi-layered soil.The governing equations for pile and soil were obtained by applying the variational principle to the potential energy, which are then solved using the software MATLAB R2021a.Comparisons were made with the published field data and the FEA.It is observed that the energy-based method presented in this study is in good agreement with the field data when compared to the linear elastic solution that does not consider soil non-linearity.The developed analytical model incorporating the Euler-Bernoulli beam theory proved to be an effective way in estimating the load-displacement responses of piles embedded in multi-layered non-linear elastic soil strata.The non-linear elastic constitutive relationship, which described the variation of secant shear modulus with strain through a power law, showed reasonably accurate predictions when compared to the published field test data and the FEA.The results indicated a difference of approximately less than 7% between the proposed method and the FEA.The study illustrates the superior advantage of the fast solutions obtained by the energy-based method over the FEA.

•
The proposed analytical models of the present study are subjected to static loads and can be extended to the effect of dynamic loading.

•
In the present study, the load-displacement responses of single piles were considered.However, the deformation responses of group piles are larger than the displacement of isolated single piles.Hence, the present work needs to be extended to understand the group action of the piles when subjected to several external loads.

Figure 1 .
Figure 1.Laterally loaded pile in layered non-linear elastic medium.

Figure 1 .
Figure 1.Laterally loaded pile in layered non-linear elastic medium.

Figure 2 .
Figure 2. Stresses and displacements within a soil continuum.

Figure 2 .
Figure 2. Stresses and displacements within a soil continuum.

Geotechnics 2022, 2 ,
FOR PEERREVIEW  4    or separation between the pile and the surrounding soil and between soil layers.The stresses and the displacement within a soil continuum are shown below in Figure2a,b.

Figure 2 .
Figure 2. Stresses and displacements within a soil continuum.

Geotechnics 2022, 2 , 80 Figures 8
Figures8 and 9show the pile displacement responses along the depth z for all the cases of the two different pile types discussed above (Table1).

Figure 8 .
Figure 8.Effect of explicit incorporation of soil characteristics and layering (two-layer soil system) for the 10 m short pile.

Figure 9 .
Figure 9.Effect of explicit incorporation of soil characteristics and layering (four-layer soil system)

Figure 8 .
Figure 8.Effect of explicit incorporation of soil characteristics and layering (two-layer soil system) for the 10 m short pile.

Figure 8 .
Figure 8.Effect of explicit incorporation of soil characteristics and layering (two-layer soil system) for the 10 m short pile.

Figure 9 .
Figure 9.Effect of explicit incorporation of soil characteristics and layering (four-layer soil system) for the 20 m long pile.

Figure 9 .
Figure 9.Effect of explicit incorporation of soil characteristics and layering (four-layer soil system) for the 20 m long pile.

Figure 10 .
Figure 10.Input of (a) pile-soil system properties and (b) boundary conditions in FEA software ANSYS 2021R1.

Figure 10 .
Figure 10.Input of (a) pile-soil system properties and (b) boundary conditions in FEA software ANSYS 2021R1.

Figure 11 .
Figure 11.The discretized FEM (a) A mesh representation of the FEM and (b) A wireframe representation of the FEM (ANSYS 2021R1).

Figure 11 .
Figure 11.The discretized FEM (a) A mesh representation of the FEM and (b) A wireframe representation of the FEM (ANSYS 2021R1).

Figure 12 .
Figure 12.Linear elastic and non-linear elastic soil (a) Pile head displacement responses for an applied force, (b) Pile head displacement responses for an applied moment.

Figure 12 .
Figure 12.Linear elastic and non-linear elastic soil (a) Pile head displacement responses for an applied force, (b) Pile head displacement responses for an applied moment.

Figure 13 .
Figure 13.Pile head lateral displacements for two different values of Poisson's ratio v = 0.2 and v = 0.49 with (a) an applied lateral force at the pile head and (b) an applied moment at the pile head.

Figure 13 .
Figure 13.Pile head lateral displacements for two different values of Poisson's ratio v = 0.2 and v = 0.49 with (a) an applied lateral force at the pile head and (b) an applied moment at the pile head.

Figure 15 .
Figure 15.Pile and soil responses obtained from the present analysis and 3D FE analysis using Equation (2) for an applied moment at the pile head (a) pile embedment depth versus the pile head displacement, (b) applied moment versus pile head displacement, (c) normalized radial soil displacement along θ = 0 • at z = 0 m, (d) normalized radial and tangential soil displacement along θ = 45 • at z = 0 m, (e) normalized tangential soil displacement along θ = 90 • at z = 0 m, (f) normalized radial and tangential soil displacement functions φ r and φ θ along θ = 45 • at z = 0 m and (g) graphical representation of the normalized secant shear modulus versus normalized radial distance along θ = 0 • , 45 • and 90 • at z = 0 m at the end of the final load increment.

Geotechnics 2022, 2 ,
FOR PEERREVIEW  24    parts: the backbone curve, the standard reload curve, the general unload curve and the direct reload curve.Allotey and El Naggar[92] simulated a reinforced concrete pile of length 12 m and diameter 0.6 m embedded in a single stratum consisting of uniform medium-stiff clay with an undrained shear strength of 50 KPa and unit weight of 19 kN/m 2 .

Figure 16 .
Figure 16.Pile head deflection with respective lateral loads.The result from the FEA underestimates the pile displacements by about 10% for the load values up to 70 kN.However, there is a good agreement in the results from the backbone curve of Allotey and El Naggar[92], FEA and the proposed energy-based solution.(c)Energy-based method versus the field test data "Soil Modulus for Laterally Loaded Piles": McClelland and Focht[19].

Figure 16 .
Figure 16.Pile head deflection with respective lateral loads.The result from the FEA underestimates the pile displacements by about 10% for the load values up to 70 kN.However, there is a good agreement in the results from the backbone curve of Allotey and El Naggar [92], FEA and the proposed energy-based solution.(c) Energy-based method versus the field test data "Soil Modulus for Laterally Loaded Piles": McClelland and Focht [19].

Table 1 .
Pile and soil properties for understanding the effect of explicit incorporation.

Table 1 .
Pile and soil properties for understanding the effect of explicit incorporation.

Table 2 .
The computational speed of the present energy-based approach and the FEA.

Table 2 .
The computational speed of the present energy-based approach and the FEA.

) Extent of Soil Layers (m) Shear Modulus, G0 (MPa) Poisson's Ratio, ν
The authors acknowledge the support received from the Florida Department of Transportation Grant: DOT-RFP-20-9069-CA, Gangals Foundation Inc. for providing the senior author with International Post Graduate Engineering Scholarship and Florida Atlantic University for providing the senior author with the Presidential Fellowship to successfully finish this study.
Author Contributions: Conceptualization, P.A.A.; methodology, P.A.A.; resources, P.A.A.; data curation, P.A.A.; writing-original draft preparation, P.A.A.; writing-review and editing, P.A.A. and M.A.; visualization, P.A.A.; supervision, M.A.All authors have read and agreed to the published version of the manuscript.Funding: Institutional Review Board Statement: Not applicable.Informed Consent Statement: Not applicable.