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

Abstract

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 finite 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 field test data as well as the finite 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 field 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.

1. 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].

2. 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₀ accompanied with a bending moment M₀. 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 H₁, H₂, H₃… 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. 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.

At each point in the soil domain, G and λ were calculated (Figure 3).

2.1. 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⁻⁵. 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].

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:

$$q=A{\left({\epsilon }_q\right)}^n$$

(1)

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₀ is constant; the second region comprising small strains starts from ε₀ 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]:

$$G=G_0{\left(\frac{{\epsilon }_q}{{\epsilon }_{q0}}\right)}^n$$

(2)

$$G=a {\epsilon }_q^n$$

(3)

where $a=\frac{G_0}{{\epsilon }_{q0}^n}$ 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). ${\epsilon }_q=\sqrt{\frac{2}{9}\left[{\left({\epsilon }_{rr}-{\epsilon }_{\theta \theta }\right)}^2+{\left({\epsilon }_{\theta \theta }-{\epsilon }_{zz}\right)}^2+{\left({\epsilon }_{zz}-{\epsilon }_{rr}\right)}^2\right]+\frac{4}{3}\left({\epsilon }_{r\theta }^2+{\epsilon }_{\theta z}^2+{\epsilon }_{zr}^2\right)}$ represents the deviatoric strain; and ε_q0 is the maximum deviatoric strain with linear elastic behavior, which is equal to 10⁻⁵. Soil stiffness G is estimated by calculating the strain at each location followed by the power law.

2.2. 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].

$$u_r=u\left(z\right) {\varphi }_r\left(r\right) cos\theta$$

(4a)

$$u_{\theta }=-u\left(z\right) {\varphi }_{\theta }\left(r\right) sin\theta$$

(4b)

$$u_z=0$$

(4c)

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_{\theta }$ 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₀) that contains the applied force vector F₀ (Figure 2). The above functions are set as ϕ_r(r) = 1 and ϕ_θ(r) = 1 at r = r ₀ (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).

2.3. Governing Differential Equation

In the derivation of the governing differential equations that can capture the non-linear 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

$$\Pi =\frac{1}{2} E_p{I}_p {{\displaystyle \int }}_0^L{\left(\frac{d^{2}u}{d{z}^2}\right)}^{2}dz+{{\displaystyle \int }}_0^{\infty }{{\displaystyle \int }}_0^{2\pi }{{\displaystyle \int }}_{r_p}^{\infty }\frac{1}{2}{\sigma }_{ij }{\epsilon }_{ij} r dr d\theta dz+{{\displaystyle \int }}_L^{\infty }{{\displaystyle \int }}_0^{2\pi }{{\displaystyle \int }}_0^{r_p}\frac{1}{2}{\sigma }_{ij }{\epsilon }_{ij} r dr d\theta dz-{F_{o}u|}_{z=0}+M_o{\frac{du}{dz}|}_{z=0}$$

(5)

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

$$\left[\begin{array}{c}{\sigma }_{rr}\\ {\sigma }_{\theta \theta }\\ {\sigma }_{zz}\\ {\tau }_{r\theta }\\ {\tau }_{rz}\\ {\tau }_{\theta z}\end{array}\right]=\left[\begin{array}{cccccc}\lambda +2G& \lambda & \lambda & 0& 0& 0\\ \lambda & \lambda +2G& \lambda & 0& 0& 0\\ \lambda & \lambda & \lambda +2G& 0& 0& 0\\ 0& 0& 0& G& 0& 0\\ 0& 0& 0& 0& G& 0\\ 0& 0& 0& 0& 0& G\end{array}\right] \left[\begin{array}{c}{\epsilon }_{rr}\\ {\epsilon }_{\theta \theta }\\ {\epsilon }_{zz}\\ {\gamma }_{r\theta }\\ {\gamma }_{rz}\\ {\gamma }_{\theta z}\end{array}\right]$$

(6)

where G and λ are the elastic constants of the soil. The strain–displacement relationship is given by

$$\left[\begin{array}{c}{\epsilon }_{rr}\\ {\epsilon }_{\theta \theta }\\ {\epsilon }_{zz}\\ {\gamma }_{r\theta }\\ {\gamma }_{rz}\\ {\gamma }_{\theta z}\end{array}\right]=\left[\begin{array}{c}-\frac{\partial {\mathrm{u}}_{\mathrm{r}}}{\partial \mathrm{r}}\\ -\frac{u_r}{r}-\frac{1}{r}\frac{\partial {\mathrm{u}}_{\theta }}{\partial \theta }\\ -\frac{\partial {\mathrm{u}}_{\mathrm{z}}}{\partial \mathrm{z}}\\ -\frac{1}{r}\frac{\partial {\mathrm{u}}_{\mathrm{r}}}{\partial \theta }-\frac{\partial {\mathrm{u}}_{\theta }}{\partial \mathrm{r}}+\frac{u_{\theta }}{r}\\ -\frac{\partial {\mathrm{u}}_{\mathrm{z}}}{\partial \mathrm{r}}-\frac{\partial {\mathrm{u}}_{\mathrm{r}}}{\partial \mathrm{z}}\\ -\frac{1}{r}\frac{\partial {\mathrm{u}}_{\mathrm{z}}}{\partial \theta }-\frac{\partial {\mathrm{u}}_{\theta }}{\partial \mathrm{z}}\end{array}\right]=\left[\begin{array}{c}-u(z)\frac{d{\varphi }_r(r)}{dr}\cos \theta \\ -u(z)\frac{{\varphi }_r(r)-{\varphi }_{\theta }(r)}{r}\cos \theta \\ u(z)\left\{\frac{{\varphi }_r(r)-{\varphi }_{\theta }(r)}{r}+\frac{d{\varphi }_{\theta }(r)}{dr}\right\}\sin \theta \\ -\frac{du(z)}{dz}{\varphi }_r(r)\cos \theta \\ \frac{du(z)}{dz}{\varphi }_{\theta }(r)\sin \theta \end{array}\right]$$

(7)

Combining Equations (6) and (7) gives the strain energy density within any layer:

$$\begin{array}{ll}\frac{1}{2}{\sigma }_{ij}{\epsilon }_{ij}=& \frac{1}{2}[(\lambda +2G)u^2{\left(\frac{d{\varphi }_r}{dr}\right)}^2{\cos }^2\theta +2\lambda u^2\frac{d{\varphi }_r}{dr}\frac{{\varphi }_r-{\varphi }_{\theta }}{r}{\cos }^2\theta \\ & +(\lambda +2G)u^2\frac{{\left({\varphi }_r-{\varphi }_{\theta }\right)}^2}{r^2}{\cos }^2\theta +G{u}^2\frac{{\left({\varphi }_r-{\varphi }_{\theta }\right)}^2}{r^2}+G{u}^2{\left(\frac{d{\varphi }_{\theta }}{dr}\right)}^2{\sin }^2\theta \\ & +2G{u}^2\frac{{\varphi }_r-{\varphi }_{\theta }}{r}\frac{d{\varphi }_{\theta }}{dr}{\sin }^2\theta +G{\left(\frac{du}{dz}\right)}^2{\varphi }_r^2{\cos }^2\theta +G{\left(\frac{du}{dz}\right)}^2{\varphi }_{\theta }^2{\sin }^2\theta ]\end{array}$$

(8)

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

$$\begin{array}{ll}\delta \Pi =& {{\displaystyle \int }}_0^L{E}_p{I}_p\frac{d^{2}u}{d{z}^2}\delta \left(\frac{d^{2}u}{d{z}^2}\right)dz\\ & +{{\displaystyle \int }}_0^{\infty }{{\displaystyle \int }}_{r_p}^{\infty }{{\displaystyle \int }}_0^{2\pi }[(\lambda +2G)u\delta u{\left(\frac{d{\varphi }_r}{dr}\right)}^2{\mathit{cos}}^2\theta +(\lambda +2G)u^2{\mathit{cos}}^2\theta \left(\frac{d{\varphi }_r}{dr}\right)\delta \left(\frac{d{\varphi }_r}{dr}\right)\\ & +2\lambda u\delta u{\mathit{cos}}^2\theta \frac{d{\varphi }_r}{dr}\frac{\left({\varphi }_r-{\varphi }_{\theta }\right)}{r}+\lambda u^2{\mathit{cos}}^2\theta \frac{1}{r}\left({\varphi }_r-{\varphi }_{\theta }\right)\delta \left(\frac{d{\varphi }_r}{dr}\right)+\lambda u^2{\mathit{cos}}^2\theta \frac{1}{r}\frac{d{\varphi }_r}{dr}\delta {\varphi }_r\\ & -\lambda u^2{\mathit{cos}}^2\theta \frac{1}{r}\frac{d{\varphi }_r}{dr}\delta {\varphi }_{\theta }+\left(\lambda +2G\right) u \delta u co{s}^2\theta \frac{1}{r^2} {\left({\varphi }_r-{\varphi }_{\theta }\right)}^2+\left(\lambda +2G\right)u^{2}co{s}^2\theta \frac{1}{r^2}\left({\varphi }_r-{\varphi }_{\theta }\right)\delta {\varphi }_r\\ & -\left(\lambda +2G\right)u^{2}co{s}^2\theta \frac{1}{r^2}\left({\varphi }_r-{\varphi }_{\theta }\right)\delta {\varphi }_{\theta }+G u \delta u {\mathit{sin}}^2\theta \frac{1}{r^2} {\left({\varphi }_r-{\varphi }_{\theta }\right)}^2+2G u \delta u {\mathit{sin}}^2\theta \frac{1}{r} \frac{d{\varphi }_{\theta }}{dr} ({\varphi }_r\\ & -{\varphi }_{\theta })+G{u}^2 {\mathit{sin}}^2\theta \frac{1}{r^2} \left({\varphi }_r-{\varphi }_{\theta }\right)\delta {\varphi }_r-G u^2{\mathit{sin}}^2\theta \frac{1}{r^2} \left({\varphi }_r-{\varphi }_{\theta }\right)\delta {\varphi }_{\theta }+G u^2 {\mathit{sin}}^2\theta \frac{1}{r} \frac{d{\varphi }_{\theta }}{dr}\delta {\varphi }_r\\ & - G u^2 {\mathit{sin}}^2\theta \frac{1}{r} \frac{d{\varphi }_{\theta }}{dr}\delta {\varphi }_{\theta }+Gu \delta u {\mathit{sin}}^2\theta {\left(\frac{d{\varphi }_{\theta }}{dr}\right)}^2\\ & +G u^2{\mathit{sin}}^2\theta \frac{d{\varphi }_{\theta }}{dr} \delta \left(\frac{d{\varphi }_{\theta }}{dr}\right)+G u^2{\mathit{sin}}^2\theta \frac{1}{r}\left({\varphi }_r-{\varphi }_{\theta }\right)\delta \left(\frac{d{\varphi }_{\theta }}{dr}\right)+G \frac{du}{dz} \delta \left(\frac{du}{dz}\right){\varphi }_r^2{\mathit{cos}}^2\theta \\ & +G {\left(\frac{du}{dz}\right)}^{2}co{s}^2\theta {\varphi }_r \delta {\varphi }_r+G \frac{du}{dz} \delta \left(\frac{du}{dz}\right){\varphi }_{\theta }^2 {\sin }^2\theta +G {\left(\frac{du}{dz}\right)}^2{\mathit{sin}}^2\theta {\varphi }_{\theta } \delta {\varphi }_{\theta }]rd\theta drdz\\ & +\pi r_p^2{{\displaystyle \int }}_L^{\infty }G\frac{du}{dz} \delta \left(\frac{du}{dz}\right)dz-{F_o\delta u|}_{z=0}+{M_o\delta \left(\frac{du}{dz}\right)|}_{z=0}=0\end{array}$$

(9)

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.

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

$$\begin{array}{ll}{{\displaystyle \int }}_0^L{E}_p{I}_p\frac{d^{2}u}{d{z}^2}\delta \left(\frac{d^{2}u}{d{z}^2}\right)dz& +\pi r_p^2{{\displaystyle \int }}_L^{\infty }G\frac{du}{dz} \delta \left(\frac{du}{dz}\right)dz\\ & +{{\displaystyle \int }}_0^{\infty }{{\displaystyle \int }}_{r_p}^{\infty }{{\displaystyle \int }}_0^{2\pi }[\left(\lambda +2G\right) u \delta u{\left(\frac{d{\varphi }_r}{dr} \right)}^2{\mathit{cos}}^2\theta +2\lambda u \delta u co{s}^2\theta \frac{d{\varphi }_r}{dr}\frac{({\varphi }_r-{\varphi }_{\theta })}{r}\\ & +Gu \delta u {\mathit{sin}}^2\theta {\left(\frac{d{\varphi }_{\theta }}{dr}\right)}^2+G \frac{du}{dz} \delta \left(\frac{du}{dz}\right){\varphi }_r^2{\mathit{cos}}^2\theta +G \frac{du}{dz} \delta \left(\frac{du}{dz}\right){\varphi }_{\theta }^2{\mathit{sin}}^2\theta \\ & +G u \delta u {\mathit{sin}}^2\theta \frac{1}{r^2} {\left({\varphi }_r-{\varphi }_{\theta }\right)}^2+\left(\lambda +2G\right) u \delta u co{s}^2\theta \frac{1}{r^2} {\left({\varphi }_r-{\varphi }_{\theta }\right)}^2]rd\theta drdz-F_o\delta u{|}_{z=0}\\ & +{M_o\delta \left(\frac{du}{dz}\right)|}_{z=0}=0 \end{array}$$

(10)

The governing equation becomes

$$E_p{I}_p\frac{d^{4}u}{d{z}^4}-C\frac{d^{2}u}{d{z}^2}+ku=0$$

(11)

where

$$C={{\displaystyle \int }}_{r_p}^{\infty }{{\displaystyle \int }}_0^{2\pi }{G}_{si}({\varphi }_{r }^2 {\cos }^2\theta +{\varphi }_{\theta }^2 {\sin }^2\theta )rd\theta$$

(12)

where i denotes the number of the layer from the surface to the length of the pile

$$C_s={{\displaystyle \int }}_{r_p}^{\infty }{{\displaystyle \int }}_0^{2\pi }{G}_{sj}({\varphi }_{r }^2 {\cos }^2\theta +{\varphi }_{\theta }^2 {\sin }^2\theta )rd\theta$$

(13)

and where j represents the layer number from pile length to infinity

$$k={{\displaystyle \int }}_{r_p}^{\infty }{{\displaystyle \int }}_0^{2\pi }(\left(\lambda +G\right){\left(\frac{d{\varphi }_r}{dr} \right)}^2{\cos }^2\theta +2\lambda \frac{1}{r}\frac{d{\varphi }_r}{dr}({\varphi }_r-{\varphi }_{\theta })co{s}^2\theta +\left(\lambda +G\right)\frac{{\left({\varphi }_r-{\varphi }_{\theta }\right)}^2}{r^2} co{s}^2\theta +G \frac{{\left({\varphi }_r-{\varphi }_{\theta }\right)}^2}{r^2} si{n}^2\theta +2G\frac{({\varphi }_r-{\varphi }_{\theta })}{r} \frac{d{\varphi }_{\theta }}{dr} si{n}^2\theta +G {\left(\frac{d{\varphi }_{\theta }}{dr}\right)}^{2}si{n}^2\theta ) rd\theta dr$$

(14)

The fourth-order differential Equation (11) can be solved using a central finite difference scheme, represented as

$$E_p{I}_p\left(\frac{u_{i-2}-4u_{i-1}+6u_i-4u_{i+1}+u_{i+2}}{\mathsf{\Delta }{z}^4} \right)-C\left(\frac{-u_{i-1}+u_{i+1}}{\mathsf{\Delta }z}\right)+k{u}_i=0$$

(15)

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₀ ≤ r ≤∞ and equating the summation to zero yields

$$\begin{array}{ll}{{\displaystyle \int }}_0^{\infty }{{\displaystyle \int }}_{r_p}^{\infty }{{\displaystyle \int }}_0^{2\pi }& [\left(\lambda +2G\right)u^2{\cos }^2\theta \left(\frac{d{\varphi }_r}{dr}\right)\delta \left(\frac{d{\varphi }_r}{dr}\right)+G{u}^2{\sin }^2\theta \frac{1}{r}\left(\frac{d{\varphi }_{\theta }}{dr}\right)\delta {\varphi }_r+\lambda u^2{\cos }^2\theta \frac{({\varphi }_r-{\varphi }_{\theta })}{r}\delta \left(\frac{d{\varphi }_r}{dr}\right)\\ & +\lambda u^2{\cos }^2\theta \frac{1}{r}\left(\frac{d{\varphi }_r}{dr}\right)\delta {\varphi }_r+G{u}^2{\sin }^2\theta \frac{1}{r^2}({\varphi }_r-{\varphi }_{\theta })\delta {\varphi }_r+\left(\lambda +2G\right)u^2{\cos }^2\theta \frac{1}{r^2} ({\varphi }_r-{\varphi }_{\theta }) \delta {\varphi }_r\\ & +G{\cos }^2\theta {\left(\frac{du}{dz}\right)}^2{\varphi }_r\delta {\varphi }_r] rdrd\theta dz=0\end{array}$$

(16)

$$\begin{array}{ll}{{\displaystyle \int }}_0^{\infty }{{\displaystyle \int }}_{r_p}^{\infty }{{\displaystyle \int }}_0^{2\pi }& [-G{u}^2{\sin }^2\theta \left(\frac{d{\varphi }_{\theta }}{dr}\right)\delta {\varphi }_{\theta }-\lambda u^2{\cos }^2\theta \frac{1}{r}\left(\frac{d{\varphi }_r}{dr}\right)\delta {\varphi }_{\theta }-G{u}^2{\sin }^2\theta \frac{1}{r^2}({\varphi }_r-{\varphi }_{\theta })\delta {\varphi }_{\theta }\\ & -\left(\lambda +2G\right)u^2{\cos }^2\theta \frac{1}{r^2} ({\varphi }_r-{\varphi }_{\theta }) \delta {\varphi }_{\theta }\\ & +G{u}^2{\sin }^2\theta \left(\frac{d{\varphi }_{\theta }}{dr}\right)\delta \left(\frac{d{\varphi }_{\theta }}{dr}\right)+G{u}^2{\sin }^2\theta \frac{({\varphi }_r-{\varphi }_{\theta })}{r}\delta \left(\frac{d{\varphi }_{\theta }}{dr}\right)+G{\sin }^2\theta {\left(\frac{du}{dz}\right)}^2{\varphi }_{\theta }\delta {\varphi }_{\theta }] rdrd\theta dz\\ & =0\end{array}$$

(17)

Simplifying Equations (16) and (17) gives

$${{\displaystyle \int }}_{r_p}^{\infty }(m_{s1}\left(\frac{d{\varphi }_r}{dr}\right)\delta \left(\frac{d{\varphi }_r}{dr}\right)+m_{s3}\frac{1}{r}\left(\frac{d{\varphi }_r}{dr}\right)\delta {\varphi }_r+m_{s3}\frac{({\varphi }_r-{\varphi }_{\theta })}{r}\delta \left(\frac{d{\varphi }_r}{dr}\right)+m_{s1}\frac{1}{r^2} ({\varphi }_r-{\varphi }_{\theta }) \delta {\varphi }_r+m_{s2}\frac{1}{r^2} ({\varphi }_r-{\varphi }_{\theta }) \delta {\varphi }_r+m_{s2}\frac{1}{r}\left(\frac{d{\varphi }_{\theta }}{dr}\right)\delta {\varphi }_r+n_{s1}{\varphi }_r\delta {\varphi }_r) dr=0$$

(18)

$${{\displaystyle \int }}_p^{\infty }(-m_{s2}\left(\frac{d{\varphi }_{\theta }}{dr}\right)\delta \left(\frac{d{\varphi }_{\theta }}{dr}\right)-m_{s3}\frac{1}{r}\left(\frac{d{\varphi }_r}{dr}\right)\delta {\varphi }_{\theta }-m_{s2}\frac{1}{r^2} ({\varphi }_r-{\varphi }_{\theta }) \delta {\varphi }_{\theta }-m_{s1}\frac{1}{r^2} ({\varphi }_r-{\varphi }_{\theta }) \delta {\varphi }_{\theta }+m_{s2}\left(\frac{d{\varphi }_{\theta }}{dr}\right)\delta \left(\frac{d{\varphi }_{\theta }}{dr}\right)+m_{s2}\frac{({\varphi }_r-{\varphi }_{\theta })}{r}\delta \left(\frac{d{\varphi }_{\theta }}{dr}\right)n_{s2}{\varphi }_{\theta }\delta {\varphi }_{\theta }) dr=0$$

(19)

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

$$m_{s1}\left(r\right)={{\displaystyle \int }}_0^{\infty }{{\displaystyle \int }}_0^{2\pi }\left({\lambda }_s+2G_s\right)u^2{\cos }^2\theta rd\theta dz={{\displaystyle \sum }}_{i=n}^n{{\displaystyle \int }}_{H_{i-1}}^{H_i}{{\displaystyle \int }}_0^{2\pi }\left({\lambda }_{si}+2G_{si}\right)u_i^2{\cos }^2\theta rd\theta dz$$

(20a)

$$m_{s2}\left(r\right)={{\displaystyle \int }}_0^{\infty }{{\displaystyle \int }}_0^{2\pi }{G}_s{u}^2{\sin }^2\theta rd\theta dz={{\displaystyle \sum }}_{i=n}^n{{\displaystyle \int }}_{H_{i-1}}^{H_i}{{\displaystyle \int }}_0^{2\pi }{G}_{si}{u}_i^2{\sin }^2\theta rd\theta dz$$

(20b)

$$m_{s3}\left(r\right)={{\displaystyle \int }}_0^{\infty }{{\displaystyle \int }}_0^{2\pi }{\lambda }_s{u}^2{\cos }^2\theta rd\theta dz={{\displaystyle \sum }}_{i=n}^n{{\displaystyle \int }}_{H_{i-1}}^{H_i}{{\displaystyle \int }}_0^{2\pi }{\lambda }_{si}{u}_i^2{\cos }^2\theta rd\theta dz$$

(20c)

$$n_{s1}\left(r\right)={{\displaystyle \int }}_0^{\infty }{{\displaystyle \int }}_0^{2\pi }{G}_s{\left(\frac{du}{dz}\right)}^2{\cos }^2\theta rd\theta dz={{\displaystyle \sum }}_{i=n}^n{{\displaystyle \int }}_{H_{i-1}}^{H_i}{{\displaystyle \int }}_0^{2\pi }{G}_{si}{\left(\frac{d{u}_i}{dz}\right)}^2{\cos }^2\theta rd\theta dz$$

(20d)

$$n_{s1}\left(r\right)={{\displaystyle \int }}_0^{\infty }{{\displaystyle \int }}_0^{2\pi }{G}_s{\left(\frac{du}{dz}\right)}^2{\sin }^2\theta rd\theta dz={{\displaystyle \sum }}_{i=n}^n{{\displaystyle \int }}_{H_{i-1}}^{H_i}{{\displaystyle \int }}_0^{2\pi }{G}_{si}{\left(\frac{d{u}_i}{dz}\right)}^2{\sin }^2\theta rd\theta dz$$

(20e)

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

$$\frac{d^2{\varphi }_r}{d{r}^2}+\frac{1}{m_{s1}}\frac{d{m}_{s1}}{dr}\frac{d{\varphi }_r}{dr}-\left\{\frac{1}{r^2} \frac{m_{s1}+m_{s2}+m_{s3}}{m_{s1}}-\frac{1}{r}\frac{1}{m_{s1}}\frac{d{m}_{s3}}{dr}+\frac{n_{s1}}{m_{s1}}\right\}{\varphi }_r=\frac{m_{s2}+m_{s3}}{m_{s1}}\frac{1}{r}\frac{d{\varphi }_{\theta }}{dr}-\left\{\frac{1}{r^2}\frac{m_{s1}+m_{s2}+m_{s3}}{m_{s1}}-\frac{1}{r}\frac{1}{m_{s1}}\frac{d{m}_{s3}}{dr}\right\}{\varphi }_{\theta }$$

(21)

$$\frac{d^2{\varphi }_{\theta }}{d{r}^2}+\frac{1}{m_{s2}}\frac{d{m}_{s2}}{dr}\frac{d{\varphi }_{\theta }}{dr}-\left\{\frac{1}{r^2} \frac{m_{s1}}{m_{s2}}+\frac{1}{r}\frac{1}{m_{s2}}\frac{d{m}_{s2}}{dr}+\frac{n_{s2}}{m_{s2}}\right\}{\varphi }_{\theta }=\frac{m_{s2}+m_{s3}}{m_{s2}}\frac{1}{r}\frac{d{\varphi }_r}{dr}-\left\{\frac{1}{r^2}\frac{m_{s1}}{m_{s2}}+\frac{1}{r}\frac{1}{m_{s2}}\frac{d{m}_{s2}}{dr}\right\}{\varphi }_r$$

(22)

Considering the soil boundary conditions ϕ_r = 1 at r = r₀; ϕ_r = 0 at r → ∞; ϕ_θ = 1 at r = r₀; and ϕ_θ = 0 at r → ∞, the governing differential equations of soil can be rewritten as

$$\frac{d^2{\varphi }_r}{d{r}^2}+y_1\frac{d{\varphi }_r}{dr}-y_2{\varphi }_r=y_3\frac{d{\varphi }_{\theta }}{dr}-y_4{\varphi }_{\theta }$$

(23)

$$\frac{d^2{\varphi }_{\theta }}{d{r}^2}+y_5\frac{d{\varphi }_{\theta }}{dr}-y_6{\varphi }_{\theta }=-y_7\frac{d{\varphi }_r}{dr}-y_8{\varphi }_r$$

(24)

where

$$y_1=\frac{1}{m_{s1}}\frac{d{m}_{s1}}{dr}$$

(25a)

$$y_2=\frac{1}{r^2}\frac{m_{s1}+m_{s2}+m_{s3}}{m_{s1}}-\frac{1}{r}\frac{1}{m_{s1}}\frac{d{m}_{s3}}{dr}+\frac{n_{s1}}{m_{s1}}$$

(25b)

$$y_3=\frac{1}{r}\frac{m_{s2}+m_{s3}}{m_{s1}}$$

(25c)

$$y_4=\frac{1}{r^2}\frac{m_{s1}+m_{s2}+m_{s3}}{m_{s1}}-\frac{1}{r}\frac{1}{m_{s1}}\frac{d{m}_{s3}}{dr}$$

(25d)

$$y_5=\frac{1}{m_{s2}}\frac{d{m}_{s2}}{dr}$$

(25e)

$$y_6=\frac{1}{r^2}\frac{m_{s1}}{m_{s2}}+\frac{1}{r}\frac{1}{m_{s2}}\frac{d{m}_{s2}}{dr}+\frac{n_{s2}}{m_{s2}}$$

(25f)

$$y_7=\frac{1}{r}\frac{m_{s2}+m_{s3}}{m_{s2}}$$

(25g)

$$y_8=\frac{1}{r^2}\frac{m_{s1}}{m_{s2}}+\frac{1}{r}\frac{1}{m_{s2}}\frac{d{m}_{s2}}{dr}$$

(25h)

All the equations mentioned in this study are solved using the software MATLAB 2021a.

3. 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₁, y₂, y₃, y₄, y₅, y₆, y₇ and y₈. 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₁, y₂, y₃, y₄, y₅, y₆, y₇ and y₈, 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₁, y₂, y₃, y₄, y₅, y₆, y₇ and y₈ 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 $\frac{{\gamma }_{old}-{\gamma }_{new}}{{\gamma }_{old}}<0.001$. Figure 7 shows the flow chart of the iterative solution procedure.

4. Results

4.1. 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. Pile and soil properties for understanding the effect of explicit incorporation.

S. No.	Pile Properties	Heights of Soil Layers, Hi (m)	Shear Moduli, G0 (MPa)	Poisson’s Ratio (v)
Short Pile
Case 1	L = 10 m rp = 0.5 m Ep = 25 × 106 kN/m2 F0 = 1000 kN	Homogeneous soil layer	G02/G01 = 1 G01 = 25	0.25
Case 2		Two-layer soil system, H1 = 2 m	G02/G01 = 2 G01 = 25; G02 = 50
Case 3		Two-layer soil system, H1 = 2 m	G02/G01 = 4 G01 = 25; G02 = 100
Case 4		Two-layer soil system, H1 = 2 m	G02/G01 = 0.5 G01 = 50; G02 = 25
Long Pile
Case 1	L = 20 m rp = 0.25 m Ep = 25 × 106 kN/m2 F0 = 1000 kN	Homogeneous soil layer	G01 = G02 = G03 = G04 G01 = 10	0.25
Case 2		Four-layer soil system, H1 = 1 m; H2 = 3 m; H3 = 5 m	G01 = 0.5G02 = 0.25G03 = 0.125G04 G01 = 10; G02 = 20; G03 = 40; G04 = 80
Case 3		Four-layer soil system, H1 = 1 m; H2 = 3 m; H3 = 5 m	G01 = 0.25G02 = 0.25G03 = 0.125G04 G01 = 10; G02 = 40; G03 = 40; G04 = 80

shows the pile and soil details considered for the analyses.

Figure 8 and Figure 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.

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

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.

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. The computational speed of the present energy-based approach and the FEA.

Problems	3D FEA (Seconds)	Energy-Based Method (Seconds)
12 (a)	4343	169
12 (b)	4389	178

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 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₀) 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.

• (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². 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 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].

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². The details of the soil strata are shown in

Table 3. Soil details at the test site investigated by McClelland and Focht [19] (adapted from Basu [90]).

Depth (m)	Extent of Soil Layers (m)	Shear Modulus, G0 (MPa)	Poisson’s Ratio, ν
2	0 to −4	1.6	0.3
6	−4 to −8	4.8	0.3
10	−8 to −12	8	0.3
17.5	−12 to greater depth	14	0.3

.

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.

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

6. Future Work

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