Lateral load capacity of flexible piles in the SDMC model with cylindrical cavity expansion theory

Abstract

This paper presents a method for calculating the lateral load capacity of flexible piles based on the state-dependent Mohr-Coulomb (SDMC) model combined with cylindrical cavity expansion theory. The soil surrounding the pile is idealized as a cylindrical cavity, and the stress-strain behavior of the soil during monotonic lateral loading of the pile is described using the cavity expansion model. By incorporating the stress equilibrium equation, displacement compatibility conditions, boundary conditions, and the constitutive model of the soil, a system of governing equations consisting of partial differential equations (PDEs) is derived. These equations are solved to obtain monotonic p-y curves that capture elastic, elastoplastic, and critical processes. Subsequently, the deflection equilibrium differential equation is introduced to transfer the lateral load from the pile top to the tip, and the resulting system of equations is solved numerically using the finite-difference method. The accuracy of the cylindrical cavity expansion theory integrated with the SDMC model is verified through finite element method (FEM) simulations. Furthermore, the validity and capacity of the proposed approach in predicting the monotonic load-displacement response of piles are demonstrated by comparing the results with existing centrifuge test data. The finding confirm that the proposed method effectively captures key phenomena observed in pile load test.

Introduction

In offshore geotechnical engineering, piles are commonly employed as bridge foundations to resist lateral loads induced by wind and waves. Fleming et al.¹ pointed out that current design practices for pile engineering often rely on empirical and semi-empirical correlations, which can lead to inconsistent and scattered predictions. In response, numerous methods have been developed to predict the lateral load capacity of piles. These approaches can be broadly classified into the following categories: (a) The beam on elastic foundation method². This approach models the pile as a beam supported by an elastic foundation, with the load-displacement behavior analyzed using the Winkler foundation model. (b) The p-y curve method³. Widely adopted in practice, this method simulates the nonlinear response of laterally loaded piles by representing the soil as a series of independent nonlinear springs, each with stiffness varying according to depth and deformation. (c) The failure-wedge method⁴ and strain-wedge method⁵. These methods idealize the surrounding soils as a wedge-shaped failure zone, according for continuity and force equilibrium within the soil mass to evaluate lateral capacity. (d) The approximate method^6,7,8. This technique offers a semi-analytical solution for predicting the nonlinear three-dimensional response of the pile-soil system by incorporating elastoplastic constitutive models. These theoretical advances have stimulated extensive research and significantly deepened the understanding of laterally loaded piles among researchers.

The relative pile-soil stiffness plays a critical role in determining the soil-pile interaction mode, as evidenced by existing centrifuge tests^9,10. Wang et al.¹¹ observed distinct soil flow mechanisms around laterally loaded piles with different length-to-diameter ratios. For flexible piles, the soil flow mechanism consists of a wedge-type failure near the ground surface and a full-slow beneath the wedge, with neither deflection nor rotation occurring at the pile base. As the length-to-diameter ratio decreases, piles behave in a semi-rigid or rigid manner, exhibiting not only bending but also rotation with noticeable toe-kick. Consequently, the simpler flow mechanism and well-defined failure boundary make the pile-soil response of flexible piles relatively easier to analyze. From a cross-sectional perspective, the soil wedge near the ground surface around a laterally loaded pile can be approximated by an ellipse shape¹². Based on this assumption of an elliptical influence zone, Pang et al.^13,14 employed the modified Cam-clay (MCC) model to capture the critical state of the soil surrounding the pile and analyzed the pile-soil interaction using cavity expansion theory. The approach offers a relatively rigorous method for predicting the monotonic lateral load-displacement behavior of offshore flexible piles.

Unlike clay soils, the analysis of sands must account for drainage condition during shearing, volumetric strains around the cavity, and the total stress path. These factors complicate the accurate representation of the peak strength and dilatancy characteristics unique to sands¹⁵. Only a limited number of studies have incorporated constitutive models with state-dependent dilatancy to describe sand behavior^16,17,18. However, these models often rely on relatively complex conceptual frameworks and mathematical formulations, limiting their direct applicability in engineering practice. It is widely recognized that the mechanical behavior of sand is state-dependent¹⁹, with both shear resistance and dilation angles being influenced by relative density and stress levels²⁰. Building upon the widely used Mohr-Coulomb failure criterion, a modified version knowns as the state-dependent Mohr-Coulomb (SDMC) model has been shown to provide satisfactory simulations of pile response under monotonic lateral loading^21,22. The key parameters of this model are formulated to reflect state-dependent behavior and are carefully calibrated against high-quality element tests conducted over relevant pressure ranges²³. In present study, the SDMC model will be employed to offer a practical and straightforward approach for capturing state-dependent response of sands.

In the present study, the zone surrounding the pile is idealized as a cylindrical cavity, and the stress-strain response of the adjacent soil is described using the state-dependent Mohr-Coulomb (SDMC) model. By converting the conventional axisymmetric cylindrical cavity expansion problem into a non-axisymmetric problem, a set of equivalent cylindrical cavities is introduced to represent the radial stresses in various directions around the pile. Utilizing the stress equilibrium equation, displacement compatibility condition, boundary conditions, and the SDMC model, a system of governing equations of first-order partial differential equations (PDEs) is established to derive the stress distributions around the pile. Then, monotonic p-y curves, which capture elastic, elastoplastic, and critical processes, are derived and subsequently incorporated into the deflection equilibrium differential equation. This system is solved numerically using the finite-difference method (FDM) to predict the load-displacement behavior of flexible laterally loaded piles. The validity of the proposed cavity expansion approach is verified through finite element method (FEM) simulations employing a user-modified Mohr-Coulomb model that incorporates state-dependent characteristics. Furthermore, the accuracy and predictive capacity of the method for estimating the lateral load capacity of piles are thoroughly evaluated using existing centrifuge test results.

Description of the problem and basic assumptions

The paper addresses the prediction of lateral load capacity and load-displacement behavior of monotonically laterally loaded piles in sandy soils. The soil movement around the pile is conceptualized through i.a. cavity expansion process. Figure 1 shows front and plane views of an investigated pile segment and the surrounding soil disk. Under a lateral load ($\:{H}_{0}$) applied at the pile head, the pile undergoes a horizontal displacement ($\:y$), inducing outward deformation of the adjacent soil. The soil mass is modeled as a semi-infinite sandy medium characterized by an initial void ratio ($\:{e}_{0}$) and an at-rest coefficient of the earth pressure ($\:{K}_{0}$), where $\:{K}_{0}={\sigma\:}_{h0}/{\sigma\:}_{v0}$. Here, $\:{\sigma\:}_{h0}$ and $\:{\sigma\:}_{v0}$ represent the initial horizontal and vertical in-situ stresses, respectively, with $\:{\sigma\:}_{v0}={\gamma\:}_{s}z$ (where $\:{\gamma\:}_{s}$ is the soil weight of the soil and $\:z$ is the depth).

Unlike traditional cavity expansion problems where the cavity center remains fixed and the radius changes, the laterally loaded pile maintains a constant radius ($\:{r}_{a}$) while the centerline displaces, as depicted in Fig. 1. For any direction defined by angle $\:\alpha\:$ with vertex at $\:O{\prime\:}{\prime\:}$, an equivalent cylindrical cavity centered at $\:O{\prime\:}$ can be defined to represent the radial stress around the pile. The radius of this equivalent cavity is given by $\:{r}_{a}^{{\prime\:}}=y\text{c}\text{o}\text{s}\alpha\:+{r}_{a}$. As $\:\alpha\:$ varies from $\:0^\circ\:$ to $\:\pm\:90^\circ\:$, $\:{r}_{a}^{{\prime\:}}$ ranges from $\:y+{r}_{a}$ to $\:{r}_{a}$, enabling the construction of a family of equivalent cavities for different direction and lateral displacements. Initially, the soil responds elastically to pile movement. With increasing load, the soil yields and enters an elastoplastic state, bounded by an equivalent elastoplastic boundary ($\:{r}_{p}^{{\prime\:}}$). The development of radial ($\:{\sigma\:}_{r}$), hoop ($\:{\sigma\:}_{\theta\:}$), and vertical ($\:{\sigma\:}_{z}$) stresses corresponds to radial ($\:{\epsilon\:}_{r}$), hoop ($\:{\epsilon\:}_{\theta\:}$), and vertical ($\:{\epsilon\:}_{z}$) strains. The displacement of an arbitrary soil element from its equivalent initial radius ($\:{r}_{0}^{{\prime\:}}$) to its current radius ($\:r{\prime\:}$) is described by the an equivalent radial displacement $\:{u}_{r}^{{\prime\:}}={r}^{{\prime\:}}-{r}_{0}^{{\prime\:}}$.

Fig. 1

Schematic diagram for the load-displacement mode of a pile element.

The following additional assumptions are adopted in this study:

• 1) The investigated pile is long and flexible. Consequently, mechanisms such as soil rotation and flow failure, which are typically associated with short rigid piles^11,24,25, are not considered. Although previous studies have indicated that pile-soil interface friction can influence the interaction behavior under lateral loading²⁶, the lateral soil resistance remains the primary component of lateral bearing capacity of the pile¹³. Thus, the pile-soil interface friction is neglected in that analysis of long-flexible piles.

• 2) Small-strain behavior is assumed in the elastic zone, with the soil obeying Hooke’s law. In contrast, large-strain assumption is adopted in the elastoplastic zone, where the stress-strain response is described by the state-dependent Mohr-Coulomb (SDMC) model.

The quasi-static stress state of an arbitrary soil element of the equivalent cavity can be given by the stress equilibrium equation in the radial direction as follows:

$$\:\frac{\text{d}{\sigma\:}_{r}}{\text{d}r{\prime\:}}+\frac{{\sigma\:}_{r}-{\sigma\:}_{\theta\:}}{r{\prime\:}}=0$$

(1)

where $\:\text{d}(\bullet\:)$ is the spatial differential of ($\:\bullet\:$) for a given time (i.e. Eulerian description).

The small-strain assumption is assumed in the elastic zone, and radial and hoop strains can be written as:

$$\:{\epsilon\:}_{r}=-\frac{\text{d}{u}_{r}^{{\prime\:}}}{\text{d}r{\prime\:}}$$

(2a)

$$\:{\epsilon\:}_{\theta\:}=-\frac{{u}_{r}^{{\prime\:}}}{r{\prime\:}}$$

(2b)

To accommodate the effect of large deformations during the elastoplastic loading process, the radial and hoop strains adopt logarithmic form, defined as:

$$\:{\epsilon\:}_{r}=-\text{l}\text{n}\left(\frac{\text{d}{r}^{{\prime\:}}}{\text{d}{r}_{0}^{{\prime\:}}}\right)$$

(3a)

$$\:{\epsilon\:}_{\theta\:}=-\text{l}\text{n}\left(\frac{{r}^{{\prime\:}}}{{r}_{0}^{{\prime\:}}}\right)$$

(3b)

$$\:{\epsilon\:}_{v}=-\text{l}\text{n}\left(\frac{v}{{v}_{0}}\right)$$

(3c)

where $\:{v}_{0}$ and $\:v$ are initial and current volumes of the soil element, respectively.

The plane-strain condition and drained condition are adopted in this study, and the volumetric strain increment ($\:\text{d}{\epsilon\:}_{v}$) can be written as:

$$\:\text{d}{\epsilon\:}_{v}=\text{d}{\epsilon\:}_{r}+\text{d}{\epsilon\:}_{\theta\:}+\text{d}{\epsilon\:}_{z}$$

(4)

where $\:\text{d}{\epsilon\:}_{z}=0$ in Eq. (4) is adopted to ensure the plane-strain condition.

Solution to the pile-soil interaction

This section provides a detailed analysis to establish the relationship between stresses and strains in the soil surrounding the pile, and to derive the corresponding p-y curve. The Hooke’s law is used to figure out the elastic behavior of the soil, while the elastoplastic response is characterized by the state-dependent Mohr-Coulomb (SDMC) model.

Elastic response of the soil in the elastic zone

According to the Hooke’s law, the stress-strain relationship of the soil in the elastic zone can be given as:

$$\:\left[\begin{array}{c}\text{D}{\epsilon\:}_{r}^{\text{e}}\\\:\text{D}{\epsilon\:}_{\theta\:}^{\text{e}}\\\:\text{D}{\epsilon\:}_{z}^{\text{e}}\end{array}\right]=\frac{1}{{E}_{s}}\left[\begin{array}{ccc}1&\:-\mu\:&\:-\mu\:\\\:-\mu\:&\:1&\:-\mu\:\\\:-\mu\:&\:-\mu\:&\:1\end{array}\right]\left[\begin{array}{c}\text{D}{\sigma\:}_{r}\\\:\text{D}{\sigma\:}_{\theta\:}\\\:\text{D}{\sigma\:}_{z}\end{array}\right]$$

(5)

where $\:\mu\:$ is the Poisson’s ratio; $\:{E}_{s}$ is the elasticity modulus of the soil; the superscript “e” means that the symbolic variable belongs to the elastic process; $\:\text{D}(\bullet\:)$ is the time differential of ($\:\bullet\:$) for a given soil particle (i.e. Lagrangian description).

The elastic solution of Eqs. (1), (2), and (5) can be expressed following the small-strain assumption as²⁷:

$$\:{\sigma\:}_{r}={\sigma\:}_{h0}+({\sigma\:}_{rp}-{\sigma\:}_{h0})\frac{{r}_{p}^{{\prime\:}2}}{{r{\prime\:}}^{2}}$$

(6a)

$$\:{\sigma\:}_{\theta\:}={\sigma\:}_{h0}-({\sigma\:}_{rp}-{\sigma\:}_{h0})\frac{{r}_{p}^{{\prime\:}2}}{{r{\prime\:}}^{2}}$$

(6b)

$$\:{\sigma\:}_{z}={\sigma\:}_{v0}$$

(6c)

$$\:{u}_{r}^{{\prime\:}}=\frac{{\sigma\:}_{rp}-{\sigma\:}_{h0}}{2{G}_{s}}\frac{{r}_{p}^{{\prime\:}2}}{r{\prime\:}}$$

(6d)

where $\:{G}_{s}$ is the shear modulus of the soil, can be given as $\:{E}_{s}=2{G}_{s}(1+\mu\:)$; $\:{\sigma\:}_{rp}$ is the radial stress at the elastoplastic boundary.

Then, integrating the radial stress ($\:{\sigma\:}_{r}$) from the infinity to a certain equivalent cavity radius ($\:{r}_{a}^{{\prime\:}}$) in the direction of an angle ($\:\alpha\:$), the lateral soil resistance ($\:{p}_{r}^{\text{e}}$) for the elastic process can be given as:

$$\:{p}_{r}^{\text{e}}={\int\:}_{\infty\:}^{{r}_{a}^{{\prime\:}}}{\Delta\:}{\sigma\:}_{r}\text{d}r{\prime\:}={\int\:}_{\infty\:}^{{r}_{a}^{{\prime\:}}}({\sigma\:}_{r}-{\sigma\:}_{h0})\text{d}r{\prime\:}$$

(7)

Finally, the distribution of the lateral soil resistance per unit ($\:{p}^{\text{e}}$) in the elastic zone can be determined as:

$$\:{p}^{\text{e}}(z,y)={\int\:}_{-\frac{\pi\:}{2}}^{\frac{\pi\:}{2}}{p}_{r}^{\text{e}}\text{c}\text{o}\text{s}\alpha\:\text{d}\alpha\:={\int\:}_{-\frac{\pi\:}{2}}^{\frac{\pi\:}{2}}\left[{\int\:}_{\infty\:}^{{r}_{a}^{{\prime\:}}}({\sigma\:}_{r}-{\sigma\:}_{h0})\text{d}r{\prime\:}\right]\text{c}\text{o}\text{s}\alpha\:\text{d}\alpha\:$$

(8)

Substituting $\:{r}_{a}^{{\prime\:}}=y\text{c}\text{o}\text{s}\alpha\:+{r}_{a}$ into Eq. (8), $\:{p}^{\text{e}}$ during elastic loading at a certain depth ($\:z$) can be given in terms of lateral displacement per unit ($\:y$).

Elastoplastic response of the soil in the elastoplastic zone

State-dependent Mohr-Coulomb (SDMC) model

The SDMC model employed in this study is compatible with the critical state theory, which asserts that the soil reaches the critical state when it deforms under extended shearing with constant stress and volume. The critical state void ratio ($\:{e}_{c}$) is a function of the mean stress ($\:p$), and can be expressed as²⁸:

$$\:{e}_{c}={e}_{c}^{0}-{\lambda\:}_{c}{\left(\frac{p}{{p}_{\text{a}\text{t}\text{m}}}\right)}^{\zeta\:}$$

(9)

where $\:{e}_{c}^{0}$, $\:{\lambda\:}_{c}$, and $\:\zeta\:$ are three material constants obtained from conventional triaxial tests, and $\:{p}_{\text{a}\text{t}\text{m}}=101\:\text{k}\text{P}\text{a}$ is the atmospheric pressure. The state parameter ($\:\zeta\:$)¹⁹ is used to capture the variation of both friction and dilation angles with current sand conditions:

$$\:\psi\:=e-{e}_{c}$$

(10)

where $\:\psi\:$, so called the “state parameter”, is a parameter representing the soil state. The state dilatancy ($\:{D}_{s}$) relates to the state parameter ($\:\psi\:$) as follows:

$$\:{D}_{s}={d}_{0}[\text{exp}\left({m}_{0}\psi\:\right)-\frac{\eta\:}{M}]$$

(11)

where $\:{d}_{0}$ and $\:{m}_{0}$ are two model parameters; $\:\eta\:$ and $\:M$ are the current and critical stress ratio, respectively. According to the Mohr-Coulomb model, failure occurs when $\:\eta\:=M$, thus, the dilatancy can be obtained as:

$$\:{D}_{s}={d}_{0}\left[\text{exp}\left(m\psi\:\right)-1\right]={d}_{0}m\psi\:={n}_{d}\psi\:$$

(12)

where $\:{n}_{d}$ is a parameter combining $\:{d}_{0}$ and $\:m$. In this way, the dilatancy angle ($\:{\phi\:}_{d}$) is expressed as a function of the state parameter:

$$\:{\phi\:}_{d}=\text{arctan}\left(-{D}_{s}\right)=\text{a}\text{r}\text{c}\text{t}\text{a}\text{n}(-{n}_{d}\psi\:)$$

(13)

Classical Mohr-Coulomb model cannot simulate the drained strain-softening behavior of dense sands. To overcome this, the state parameter is introduced into the internal friction angle ($\:\phi\:$):

$$\:\phi\:=\text{a}\text{r}\text{c}\text{t}\text{a}\text{n}[\text{t}\text{a}\text{n}({\phi\:}_{c})\text{e}\text{x}\text{p}(-{n}_{b}\psi\:\left)\right]$$

(14)

where $\:{\phi\:}_{c}$ is the critical-state internal friction angle, and $\:{n}_{b}$ is a parameter associated with $\:\phi\:$.

The elastic modulus of the soil in the SDMC model is related to the mean stress ($\:p$) as:

$$\:{E}_{s}={E}_{i}{\left({p}_{\text{a}\text{t}\text{m}}p\right)}^{0.5}$$

(15)

where $\:{E}_{i}$ is the initial elastic modulus of the soil.

Then, with the classical Mohr-Coulomb model, the yield function is given as:

$$\:f={\sigma\:}_{1}-\alpha\:{\sigma\:}_{3}-Y=0$$

(16)

where $\:\alpha\:=(1+\text{s}\text{i}\text{n}\phi\:)/(1-\text{s}\text{i}\text{n}\phi\:)$, $\:Y=2c\text{c}\text{o}\text{s}\phi\:/(1-\text{s}\text{i}\text{n}\phi\:)$; $\:c$ is the cohesion force of the soil; $\:{\sigma\:}_{1}$ and $\:{\sigma\:}_{3}$ are maximum and minimum stresses, respectively.

The plastic potential surface of the SDMC model can be given as:

$$\:g={\sigma\:}_{1}-{\alpha\:}_{d}{\sigma\:}_{3}=0$$

(17)

where $\:{\alpha\:}_{d}=(1+\text{s}\text{i}\text{n}{\phi\:}_{d})/(1-\text{s}\text{i}\text{n}{\phi\:}_{d})$.

Governing equations in the elastoplastic zone

The total strain increment of the soil ($\:\text{D}\varvec{\epsilon\:}$) can be decomposed into its elastic ($\:\text{D}{\varvec{\epsilon\:}}^{\text{e}}$) and plastic ($\:\text{D}{\varvec{\epsilon\:}}^{\text{p}}$) components. The elastic component is represented by the Hooke’s law from Eq. (5), while the plastic component is modelled using the plastic flow rule of the SDMC model, which can be expressed as follows:

$$\:\text{D}{\varvec{\epsilon\:}}^{\text{p}}={\lambda\:}_{S}\frac{\partial\:g}{\partial\:\varvec{\sigma\:}}$$

(18)

where $\:{\lambda\:}_{S}$ is the scalar multiplier and $\:\frac{\partial\:g}{\partial\:\varvec{\sigma\:}}$ is the vector normal to the yield surface. In the general space, $\:\frac{\partial\:g}{\partial\:\varvec{\sigma\:}}$ can be decomposed into maximum, intermediate, and minimum directions, given as:

$$\:\text{D}{g}_{1}=\frac{\partial\:g}{\partial\:{\sigma\:}_{1}}=1-\frac{\partial\:{\alpha\:}_{d}\left(p\right)}{\partial\:{\sigma\:}_{1}}{\sigma\:}_{3}=1-\frac{\partial\:{\alpha\:}_{d}\left(p\right)}{\partial\:\psi\:}\frac{\partial\:\psi\:}{\partial\:{\phi\:}_{d}}\frac{\partial\:{\phi\:}_{d}}{\partial\:p}\frac{\partial\:p}{\partial\:{\sigma\:}_{1}}{\sigma\:}_{3}=1+{g}_{s}\left(p\right){\sigma\:}_{3}$$

(19a)

$$\:\text{D}{g}_{2}=\frac{\partial\:g}{\partial\:{\sigma\:}_{2}}=-\frac{\partial\:{\alpha\:}_{d}\left(p\right)}{\partial\:{\sigma\:}_{2}}{\sigma\:}_{3}=-\frac{\partial\:{\alpha\:}_{d}\left(p\right)}{\partial\:\psi\:}\frac{\partial\:\psi\:}{\partial\:{\phi\:}_{d}}\frac{\partial\:{\phi\:}_{d}}{\partial\:p}\frac{\partial\:p}{\partial\:{\sigma\:}_{2}}{\sigma\:}_{3}={g}_{s}\left(p\right){\sigma\:}_{3}$$

(19b)

$\:\text{D}{g}_{3}=\frac{\partial\:g}{\partial\:{\sigma\:}_{3}}=-\left[\frac{\partial\:{\alpha\:}_{d}\left(p\right)}{\partial\:{\sigma\:}_{3}}{\sigma\:}_{3}+{\alpha\:}_{\psi\:}\left(p\right)\right]$$\:=-\left[\frac{\partial\:{\alpha\:}_{d}\left(p\right)}{\partial\:\psi\:}\frac{\partial\:\psi\:}{\partial\:{\phi\:}_{d}}\frac{\partial\:{\alpha\:}_{d}}{\partial\:p}\frac{\partial\:p}{\partial\:{\sigma\:}_{3}}{\sigma\:}_{3}+{\alpha\:}_{d}\left(p\right)\right]=-{\alpha\:}_{d}\left(p\right)+{g}_{s}\left(p\right){\sigma\:}_{3}$(19c).

where $\:{\sigma\:}_{2}$ is the intermediate sterss; $\:{g}_{s}\left(p\right)$ is given as:

$$\:{g}_{s}\left(p\right)=\frac{2{n}_{d}\lambda\:\zeta\:cos\psi\:}{3{\left(1-\text{s}\text{i}\text{n}\psi\:\right)}^{2}\left(1+{n}_{d}^{2}{\phi\:}_{d}^{2}\right)p}{\left(\frac{p}{{p}_{atm}}\right)}^{\zeta\:}$$

(20)

In this way, combining Eqs. (5), (18) and (19a-c), the elastoplastic relationship of the soil can be given as:

$$\:\left[\begin{array}{c}\text{D}{\epsilon\:}_{1}\\\:\text{D}{\epsilon\:}_{2}\\\:\text{D}{\epsilon\:}_{3}\end{array}\right]=\frac{1}{{E}_{s}}\left[\begin{array}{ccc}1&\:-\mu\:&\:-\mu\:\\\:-\mu\:&\:1&\:-\mu\:\\\:-\mu\:&\:-\mu\:&\:1\end{array}\right]\left[\begin{array}{c}\text{D}{\sigma\:}_{r}\\\:\text{D}{\sigma\:}_{\theta\:}\\\:\text{D}{\sigma\:}_{z}\end{array}\right]+{\lambda\:}_{s}\left[\begin{array}{c}\text{D}{g}_{1}\\\:\text{D}{g}_{2}\\\:\text{D}{g}_{3}\end{array}\right]$$

(21a)

and Eq. (21a) can be inversely expressed in terms of the stiffness matrix and incremental strain components as:

$$\:\left[\begin{array}{c}\text{D}{\sigma\:}_{1}\\\:\text{D}{\sigma\:}_{2}\\\:\text{D}{\sigma\:}_{3}\end{array}\right]=\left[\begin{array}{ccc}{H}_{1}&\:{H}_{2}&\:{H}_{2}\\\:{H}_{2}&\:{H}_{1}&\:{H}_{2}\\\:{H}_{2}&\:{H}_{2}&\:{H}_{1}\end{array}\right]\left[\begin{array}{c}\text{D}{\epsilon\:}_{1}-{\lambda\:}_{S}\text{D}{g}_{1}\\\:\text{D}{\epsilon\:}_{2}-{\lambda\:}_{S}\text{D}{g}_{2}\\\:\text{D}{\epsilon\:}_{3}-{\lambda\:}_{S}\text{D}{g}_{3}\end{array}\right]$$

(21b)

which can be further derived as:

$$\:\left[\begin{array}{c}\text{D}{\sigma\:}_{1}\\\:\text{D}{\sigma\:}_{2}\\\:\text{D}{\sigma\:}_{3}\end{array}\right]=-{\lambda\:}_{S}\left[\begin{array}{c}{{H}_{1}\text{D}g}_{1}+{H}_{2}\left({\text{D}g}_{2}+{\text{D}g}_{3}\right)\\\:{{H}_{1}\text{D}g}_{2}+{H}_{2}\left({\text{D}g}_{1}+{\text{D}g}_{3}\right)\\\:{{H}_{1}\text{D}g}_{3}+{H}_{2}\left({\text{D}g}_{1}+{\text{D}g}_{2}\right)\end{array}\right]$$

(21c)

For the cavity expansion problem, $\:{\sigma\:}_{r}>{\sigma\:}_{z}>{\sigma\:}_{\theta\:}$ can be obtained, thus, $\:{\sigma\:}_{1}={\sigma\:}_{r}$, $\:{\sigma\:}_{2}={\sigma\:}_{z}$, and $\:{\sigma\:}_{3}={\sigma\:}_{\theta\:}$ are assumed in this study.

Equation (21c) is the governing equation system in the elastoplastic zone. The current stress can be given as $\:\varvec{\sigma\:}={\varvec{\sigma\:}}_{i}+\text{D}\varvec{\sigma\:}$, where $\:{\varvec{\sigma\:}}_{i}$ is the initial stress vector, and with the yield function in Eq. (17), the scalar multiplier ($\:{\lambda\:}_{s}$) in Eq. (18) can be expressed as:

$$\:{\lambda\:}_{S}=\frac{{\sigma\:}_{1i}-{\sigma\:}_{3i}\alpha\:+2c\sqrt{\alpha\:}}{\left[{{\beta\:}_{1}Dg}_{1}+{\beta\:}_{2}\left({\text{D}g}_{2}+{\text{D}g}_{3}\right)\right]-\left[{{\beta\:}_{1}\text{D}g}_{3}+{\beta\:}_{2}\left({\text{D}g}_{1}+{\text{D}g}_{2}\right)\right]\alpha\:}$$

(22)

where $\:{\beta\:}_{1}={K}_{s}+\frac{4}{3}{G}_{s}$ and $\:{\beta\:}_{2}={K}_{s}-\frac{2}{3}{G}_{s}$, and $\:{K}_{s}=\frac{2(1+\mu\:)}{3(1-2\mu\:){G}_{s}}$ is the bulk modulus of the soil.

Hybrid Eulerian-Lagrangian (HEL) approach

A number of analytical and semi-analytical cavity expansion solutions have been developed to solve the governing equations for stress-strain relationships in the elastoplastic zone, e.g. Equation (21c). These solution methods for this typical boundary value problem can be broadly classified into the so-called Eulerian approach and Lagrangian approach. The Eulerian approach leverages the self-similar nature of cavity expansion, wherein all soil particles experience identical stress and deformation paths^29,30,31. The Lagrangian approach directly integrates the incremental form of constitutive equations to establish a relationship between stresses and total strains^32,33,34. The present study adopts a novel solving strategy termed the “hybrid Eulerian-Lagrangian approach” (HEL)^35,36,37 to solve the system of governing equation given in Eq. (21c). This approach combines features of both Eulerian and Lagrangian approaches, allowing the solution process to account for temporal and spatial variations in a general form. As a result, the HEL approach offers a rigorous and efficient way for solving partial differential equations (PDEs) in Eq. (21c. It is worth noting that Eq. (21c) can also be expressed in terms of stress and strain increments, which are related to the plastic modulus and stiffness matrix. This incremental form enables the use of the Eulerian approach, as demonstrated by Pang et al.^13,14, which offers both computational convenience and higy accuracy.

The radial, hoop, and vertical strains provided in Eqs. (3a-c) are not independent and should satisfy Eq. (4), thereby giving the compatibility equation in terms of $\:r{\prime\:}$ and $\:v$ as:

$$\:\text{d}r{\prime\:}=\frac{v{r}_{0}^{{\prime\:}}}{{v}_{0}r{\prime\:}}\text{d}{r}_{0}^{{\prime\:}}$$

(23)

Combining Eqs. (1) and (23), the stress equilibrium equation can be transformed to be the expression of $\:{\sigma\:}_{r}$ in the Eulerian description as:

$$\:\text{d}{\sigma\:}_{r}=({\sigma\:}_{\theta\:}-{\sigma\:}_{r})\frac{v{r}_{0}^{{\prime\:}}}{{v}_{0}r{\prime\:}}\text{d}{r}_{0}^{{\prime\:}}$$

(24)

Substituting Eqs. (3b-c), and (4) into Eq. (21c), the hoop stress, vertical stress, and specific volume can be expressed in the incremental form:

$$\:\left[\begin{array}{c}\text{D}{\sigma\:}_{z}\\\:\text{D}{\sigma\:}_{\theta\:}\end{array}\right]={\left[\begin{array}{cc}1&\:-\mu\:\\\:-\mu\:&\:1\end{array}\right]}^{-1}\left[\begin{array}{c}\mu\:\text{D}{\sigma\:}_{r}-{E}_{s}{\lambda\:}_{s}\text{D}{g}_{3}-{E}_{s}\frac{\text{D}r{\prime\:}}{r{\prime\:}}\\\:\mu\:\text{D}{\sigma\:}_{r}-{E}_{s}{\lambda\:}_{s}\text{D}{g}_{2}\end{array}\right]$$

(25)

$$\:\frac{\text{D}v}{v}=-\frac{1-2\mu\:}{E}\left(\text{D}{\sigma\:}_{r}+\text{D}{\sigma\:}_{\theta\:}+\text{D}{\sigma\:}_{z}\right)+{\lambda\:}_{s}(\text{D}{g}_{1}+\text{D}{g}_{2}+\text{D}{g}_{3})$$

(26)

Five PDEs for the analysis in the elastoplastic zone have been obtained by the combination use of Eulerian description from Eqs. (23) and (24) and Lagrangian description from Eqs. (25) and (26). Consequently, stresses and strains in the elastoplastic zone can be calculated with the information at the elastoplastic boundary. However, both the material time derivative (Lagrangian description) and the spatial derivative (Eulerian description) are involved in the governing PDEs for the elastoplastic cavity expansion analysis, and they cannot be transformed into ordinary differential equations (ODEs). With the HEL approach, the cylinder of soil is discretized into ($\:m-1$) concentric annuli, where $\:m$ represents the number of nodes. Simultaneously, the loading process is divided into a number of continuous load steps. Each node is marked by its initial position as $\:{r{\prime\:}}_{\left(i\right)}^{\left(0\right)}$, where the subscript $\:i=1,\:2,\:3,\dots\:,\:m$ denotes the $\:i$th node and the superscript denotes the number of load step. The distribution of nodes is set to follow the nonlinear function in this study:

$$\:{r}_{\left(i\right)}^{{\prime\:}\left(0\right)}={\left(\frac{{r}_{\text{m}\text{a}\text{x}}^{{\prime\:}}}{{r}_{a}^{{\prime\:}}}\right)}^{\frac{1}{m-1}}{r}_{\left(i\right)}^{{\prime\:}\left(0\right)}$$

(27)

An information vector ($\:{\varvec{x}}_{\left(i\right)}^{\left(j\right)}$) can be defined for the $\:i$th node at the $\:j$th load step, including equivalent radial location, stresses, and deformation conditions, as:

$$\:{\varvec{x}}_{\left(i\right)}^{\left(j\right)}={\left[{r}_{\left(i\right)}^{{\prime\:}\left(j\right)},{\sigma\:}_{r\left(i\right)}^{\left(j\right)},{\sigma\:}_{\theta\:\left(i\right)}^{\left(j\right)},{\sigma\:}_{z\left(i\right)}^{\left(j\right)},{v}_{\left(i\right)}^{\left(j\right)}\right]}^{\text{T}}$$

(28)

where the superscript “$\:j$” represents the $\:j$th load step.

At the $\:j$th step, the increment of $\:\varvec{x}$ from node ($\:i+1$) to node ($\:i$) is defined as:

$$\:\text{d}\varvec{x}={\varvec{x}}_{\left(i\right)}^{\left(j\right)}-{\varvec{x}}_{(i+1)}^{\left(j\right)}$$

(29)

Similarly, for the $\:i$th node upon loading from the load step ($\:j-1$) to load step ($\:j$), $\:\text{D}\varvec{x}$ equals to

$$\:\text{D}\varvec{x}={\varvec{x}}_{\left(i\right)}^{\left(j\right)}-{\varvec{x}}_{\left(i\right)}^{(j-1)}$$

(30)

where $\:{\varvec{x}}_{\left(i\right)}^{(j-1)}$ is known from the previous step of loading. Having calculated $\:{r}_{\left(i\right)}^{{\prime\:}\left(j\right)}$ and $\:{\sigma\:}_{r\left(i\right)}^{\left(j\right)}$ from Eqs. (23), (24), and (29), $\:\text{D}{r}_{\left(i\right)}^{{\prime\:}\left(j\right)}$ and $\:\text{D}{\sigma\:}_{r\left(i\right)}^{\left(j\right)}$ can be known by Eq. (30). Then, $\:\text{D}{\sigma\:}_{\theta\:}$, $\:\text{D}{\sigma\:}_{z}$, and $\:\text{D}v$ can be determined from Eqs. (25) and (26) with $\:\varvec{x}={\varvec{x}}_{\left(i\right)}^{(j-1)}$, and $\:{\sigma\:}_{\theta\:\left(i\right)}^{\left(j\right)}$, $\:{\sigma\:}_{z\left(i\right)}^{\left(j\right)}$, and $\:{v}_{\left(i\right)}^{\left(j\right)}$ can also be calculated. The detailed solving procedure can refer to Yang et at^35,36,37.

The initial value for the solving procedure of the HEL approach is the same as the soil state at the elastoplastic boundary, derived from Eqs. (6) and (16) as:

$$\:{\sigma\:}_{rp}={\sigma\:}_{h0}+\frac{Y+(\alpha\:-1){\sigma\:}_{h0}}{\alpha\:+1}$$

(31a)

$$\:{\sigma\:}_{\theta\:p}={\sigma\:}_{h0}-\frac{Y+(\alpha\:-1){\sigma\:}_{h0}}{\alpha\:+1}$$

(31b)

$$\:{\sigma\:}_{zp}={\sigma\:}_{v0}$$

(31c)

$$\:{u}_{rp}^{{\prime\:}}=\frac{{\sigma\:}_{rp}-{\sigma\:}_{h0}}{2{G}_{s}}$$

(31d)

where $\:{\sigma\:}_{\theta\:p}$ and $\:{\sigma\:}_{zp}$ are hoop and vertical stresses at the elastoplastic boundary, respectively; $\:{u}_{rp}^{{\prime\:}}$ is the equivalent radial displacement of the soil at the elastoplastic boundary.

In this way, the lateral soil resistance in the elastoplastic zone can be determined by integrating the radial stress ($\:{\sigma\:}_{r}$) from the equivalent radius of the equivalent elastoplastic boundary ($\:{r}_{p}^{{\prime\:}}$) to a certain equivalent cavity radius ($\:{r}_{a}^{{\prime\:}}$) in the direction of the angle ($\:\alpha\:$). Then, together with the lateral soil resistance in the elastic zone, the lateral soil resistance ($\:{p}_{r}^{\text{p}}$) for the elastoplastic process can be given as:

$$\:{p}_{r}^{\text{p}}={\int\:}_{\infty\:}^{{r}_{a}^{{\prime\:}}}{\Delta\:}{\sigma\:}_{r}\text{d}r{\prime\:}={\int\:}_{{r}_{p}^{{\prime\:}}}^{{r}_{a}^{{\prime\:}}}({\sigma\:}_{r}-{\sigma\:}_{rp})\text{d}r{\prime\:}+{\int\:}_{\infty\:}^{{r}_{p}^{{\prime\:}}}({\sigma\:}_{rp}-{\sigma\:}_{h0})\text{d}r{\prime\:}$$

(32)

The distribution of the lateral soil resistance per unit ($\:{p}^{\text{p}}$) for the elastoplastic process is determined as:

$$\:{p}^{\text{p}}(z,y)={\int\:}_{-\frac{\pi\:}{2}}^{\frac{\pi\:}{2}}{p}_{r}^{p}\text{c}\text{o}\text{s}\alpha\:\text{d}\alpha\:={\int\:}_{-\frac{\pi\:}{2}}^{\frac{\pi\:}{2}}[{\int\:}_{{r}_{p}^{{\prime\:}}}^{{r}_{a}^{{\prime\:}}}({\sigma\:}_{r}-{\sigma\:}_{rp})\text{d}r{\prime\:}+{\int\:}_{\infty\:}^{{r}_{p}^{{\prime\:}}}({\sigma\:}_{rp}-{\sigma\:}_{h0})\text{d}r{\prime\:}]\text{c}\text{o}\text{s}\alpha\:\text{d}\alpha\:$$

(33)

Substituting $\:{r}_{a}^{{\prime\:}}={r}_{a}+y\text{c}\text{o}\text{s}\alpha\:$ into Eq. (33), the lateral soil resistance per unit ($\:{p}^{\text{p}}$) during the elastoplastic process at a certain depth can be determined in terms of lateral displacement per unit ($\:y$).

Monotonic p-y curve

From the above analysis, the soil around the pile deforms elastically at the initial loading, then it occurs elastoplastic deformation with the lateral load increasing. Combining Eqs. (8) and (33), the total monotonic p-y curve for a point on the pile at a specific depth ($\:z$) can be derived as follows:

$$\:{p}_{a}=\left\{\begin{array}{c}{\int\:}_{-\frac{\pi\:}{2}}^{\frac{\pi\:}{2}}\left[{\int\:}_{\infty\:}^{{r}_{a}^{{\prime\:}}}{\Delta\:}{\sigma\:}_{r}\text{d}r{\prime\:}\right]\text{c}\text{o}\text{s}\alpha\:\text{d}\alpha\:,\:\:\:(p<{p}_{a,p})\\\:{\int\:}_{-\frac{\pi\:}{2}}^{\frac{\pi\:}{2}}[{\int\:}_{{r}_{p}^{{\prime\:}}}^{{r}_{a}^{{\prime\:}}}{\Delta\:}{\sigma\:}_{r}\text{d}r{\prime\:}+{\int\:}_{\infty\:}^{{r}_{p}^{{\prime\:}}}{\Delta\:}{\sigma\:}_{r}\text{d}r{\prime\:}]\text{c}\text{o}\text{s}\alpha\:\text{d}\alpha\:,\:\:\:\:(p\ge\:{p}_{a,p})\end{array}\right.$$

(34)

where $\:{p}_{a,p}$ is the lateral soil resistance per unit when the soil at the soil-pile interface yields. As shown in Fig. 2, the monotonic p-y curve encompass elastic, elastoplastic, and critical parts, with a plastic point ($\:{y}_{a,p},{p}_{a,p}$) and a critical point ($\:{y}_{a,u},{p}_{a,u}$).

Fig. 2

Schematic diagram for the p-y curve.

Solving procedure

The lateral load capacity of the pile is reflected by the pile-head load-displacement behavior, i.e. $\:{H}_{0}-{y}_{0}$ curve, which can be determined by the p-y curve proposed in Eq. (34). An equilibrium equation for pile deflection can be derived by considering the balance of lateral forces on any analytical element along the vertical direction of the pile at a certain depth ($\:z$) as follows:

$$\:{E}_{p}{I}_{p}\frac{{\text{d}}^{4}y}{\text{d}{z}^{4}}+{p}_{a}\left(z\right)=0$$

(35)

where $\:{E}_{p}$ is the Young’s modulus of the pile; $\:{I}_{p}=\frac{\pi\:{r}_{a}^{4}}{4}$ is the inertia moment of the pile, $\:{E}_{p}{I}_{p}$ is the flexural rigidity of the pile.

Fig. 3

Schematic diagram for the p-y curve.

To solve the equation system comprising Eqs. (34) and (35), the finite-difference method (FDM) is used. The investigated pile is divided into $\:n$ segments of equal length $\:h=L/n$, where $\:L$ is the pile length, as shown in Fig. 3. In this method, the differential equation in Eq. (35) can be converted into a difference equation as:

$$\:\left({y}_{i-2}-4{y}_{i-1}+6{y}_{i}-4{y}_{i+1}+{y}_{i+2}\right)+\frac{{h}^{4}{p}_{a}(z,{y}_{i})}{{E}_{p}{I}_{p}}=0$$

(36)

where the subscript “$\:i$” means the $\:i$th node of the pile. All ($\:n+1$) points for $\:n$ segments along the pile should satisfy Eq. (36). However, lateral displacements of two virtual nodes (i.e. $\:{y}_{-2}$ and $\:{y}_{-1}$) appear above the node at the pile top when $\:i=0$, and those (i.e. $\:{y}_{n+1}$ and $\:{y}_{n+2}$) appear below the node at the pile tip when $\:i=n$. At this time, four equations of boundary conditions must be added to solve the equation system in Eq. (36), as:

$$\:{E}_{p}{I}_{p}\frac{{\text{d}}^{3}y}{\text{d}{z}^{3}}=\frac{{E}_{p}{I}_{p}}{2{h}^{3}}(-{y}_{-2}+2{y}_{-1}-2{y}_{1}+{y}_{2})={H}_{0}$$

(37a)

$$\:{E}_{p}{I}_{p}\frac{{\text{d}}^{2}y}{\text{d}{z}^{2}}=\frac{{E}_{p}{I}_{p}}{{h}^{2}}({y}_{-1}-2{y}_{0}+{y}_{1})={M}_{a,0}$$

(37b)

$$\:{y}_{n}=0$$

(37c)

$$\:{E}_{p}{I}_{p}\frac{{\text{d}}^{2}y}{\text{d}{z}^{2}}=\frac{{E}_{p}{I}_{p}}{{h}^{2}}\left({y}_{n-1}-2{y}_{n}+{y}_{n+1}\right)={M}_{a,n}=0$$

(37d)

where $\:{M}_{a,0}$ and $\:{M}_{a,n}$ are pile-head and pile-tip bead moments, respectively. Since the investigated pile is long and flexible, the pile-tip is assumed to be fixed with no moment ($\:{M}_{a,n}=0$) and lateral displacement ($\:{y}_{n}=0$).

Finally, combining Eqs. (34), (36), and (37), the lateral load capacity of the pile under the pile-head lateral load ($\:{H}_{0}$) can be determined following standard procedure of the finite-difference method, which is provided in Fig. 4.

Fig. 4

The solving procedure for the proposed approach for predicting the lateral load capacity of flexible piles in sands.

Validation of the proposed approach

Cavity expansion solutions against FEM simulation results

A finite element method (FEM) simulation was performed using the ABAQUS to model the cavity expansion problem in sand, employing a user-modified Mohr-Coulomb model that incorporates state-dependent characteristics (referred to as the SDMC model). Unlike the classical Mohr-Coulomb model, the modified version adopted in this study follows the William Mohr-Coulomb formulation, with detailed equations provided in the Appendix³⁸. The parameters and constants used in the SDMC model correspond to those NE34 sand, as reported by Ye et al.²³ and summarized in Table 1.The state parameter ($\:\psi\:$) was varied within the range of $\:[-\text{0.5,0.5}]$ to establish linear relationships with the friction angle ($\:\phi\:$) and the dilation angle ($\:{\phi\:}_{d}$), given by $\:\phi\:=-55.23\psi\:+34.89$ and $\:{\phi\:}_{d}=-86.31\psi\:$, respectively.

Table 1 Soil characteristics and constants in the SDMC model for the NE34 sand²³.

This study adopts the modeling methodology introduced by Zhou et al.³⁹ and further developed by Pang et al.¹⁵. Given the axisymmetric nature of the cavity expansion problem, the simulation model is established using the axisymmetric shell elements, as shown in Fig. 5. The initial radius of the cavity wall ($\:{r}_{a}$) is set to 1, and the potential boundary effects are minimized by extending the outer boundary to a distance of 400$\:{r}_{a}$. To enhance numerical accuracy, the mesh is refined with smaller elements near the cavity wall. The final model comprises 5360 with eight-node axisymmetric reduced-integration (CAX8R) elements, and the initial void ratio is uniformly defined as $\:{e}_{0}=0.63$.

Fig. 5

FEM model for cylindrical cavity expansion.

The numerical simulation consists of three main steps. First, the initial stress field is applied across the entire soil domain using the “predefined field” option in ABAQUS. For the case with $\:{K}_{0}=0.426$, the vertical and horizontal stresses are set to $\:{\sigma\:}_{v0}=158\:\text{k}\text{P}\text{a}$ and $\:{\sigma\:}_{h0}=67.308\:\text{k}\text{P}\text{a}$, respectively. An additional initial stress condition with $\:{\sigma\:}_{v0}=37.308\:\text{k}\text{P}\text{a}$ and $\:{K}_{0}=1$ is also defined to examine the effect of the in-situ coefficient of the earth pressure ($\:{K}_{0}$). The corresponding initial mean stresses ($\:{p}_{0}$) are calculated as 146.8 kPa and 67.308 kPa for $\:{K}_{0}=0.426$ and $\:1.0$, respectively. Additionally, key parameters are introduced via user-defined keywords: for $\:{K}_{0}=0.426$, $\:{\phi\:}_{0}=-0.1404$ and $\:{E}_{s0}=30443\:\text{k}\text{P}\text{a}$; for $\:{K}_{0}=1$, $\:{\phi\:}_{0}=-0.1635$ and $\:{E}_{s0}=19835\:\text{k}\text{P}\text{a}$, providing initial values for the friction angle ($\:\phi\:$) and elastic modulus ($\:{E}_{s}$). Second, the upper and lower boundaries of the cavity are restricted by fixing displacements in the normal direction. A radial stress equivalent to the initial total radial stress of the soil is applied to the outer boundary to achieve geostatic stress equilibrium. Third, the constraints on the cavity wall are released, and a displacement-controlled boundary with a uniform radial displacement is applied to simulate the cavity expansion process. Comparison between the predictions of this study and the FEM simulations, including stress component distributions and load-displacement response at the cavity wall, are plotted in Figs. 6 and 7, respectively.

Fig. 6

Comparison of stress distributions between results of this study and FEM for $\:{r}_{a}^{{\prime\:}}=1.11{r}_{a}$.

Fig. 7

Comparison of load-displacement curves at cavity wall between results of this study and FEM.

As depicted in Fig. 6, the soil exhibits elastic deformation beyond the equivalent elastoplastic boundary ($\:{r}_{p}^{{\prime\:}}$), which is approximately located at $\:8{r}_{a}^{{\prime\:}}/{r}_{a}$, and elastoplastic deformation within this boundary. The distributions of radial stress ($\:{\sigma\:}_{r}$) and hoop stress ($\:{\sigma\:}_{\theta\:}$) obtained from the present approach align closely with those from FEM simulations. However, a noticeable discrepancy is observed in the distributions of vertical stress ($\:{\sigma\:}_{z}$). This difference arises because the vertical stress-strain relationship is governed by the Hooke’s law in this study, whereas the ABAQUS simulation employs a user-modified William Mohr-Coulomb. As illustrated in Fig. 7, the results for the case with $\:{K}_{0}=0.426$ show good agreement, while a discrepancy is evident for $\:{K}_{0}=1.0$. This deviation can be attributed to differences in the mean stress ($\:p$) predicted by the SDMC model used in this study and user-modified William Mohr-Coulomb model implemented in ABAQUS, as described by Eqs. (9) and (15). Overall, the present approach demonstrates a strong capability to predict both the stress distribution and load-displacement behavior of the soil surrounding a laterally loaded pile.

Solutions of load-displacement response against centrifuge tests

Fontainebleau sand was employed in centrifuge tests conducted by Rosquoët et al.^40,41 to investigate the behavior of piles under cyclic lateral loading. Huang⁴² conducted triaxial tests to provide constants in SDMC model for the Fontainebleau sand, as summarized in Table 2. The Fontainebleau sand was prepared in a dense state with a relative density $\:{D}_{r}=86\%$ and an initial void ratio $\:{e}_{0}=0.63$. The tests were carried out under an acceleration of 40 g, simulating a prototype pile with a radius ($\:{r}_{a}$) of 0.36 m, a total length of 14.6 m, and a bending stiffness ($\:{E}_{p}{I}_{p}$) of 476 MN$\:\bullet\:$m². The prototype’s embedded length ($\:L$) was 12 m, and the cyclic lateral load was applied at a prototype elevation of 1.6 m above the ground surface. The unit weight of the soil ($\:{\gamma\:}_{s}$) was 16 kN/m³, and the in-situ coefficient of the earth pressure ($\:{K}_{0}$) was set to 0.426. Two centrifuge tests were conducted on a single pile subjected to one-way cyclic (test P344) and partial one-way cyclic (test P36) lateral loading, as reported by Rosquoët et al.^40,41. Since the present study focuses on monotonic lateral response, the data from the first loading cycle are used to validate the proposed approach.

Table 2 Soil characteristics and constants in the SDMC model for the Fontainebleau sand⁴².

Figure 8 compares the pile-head load ($\:{H}_{0}$) versus displacement ($\:{y}_{0}$) curves obtained from the present approach and centrifuge test, showing generally good agreement between the two. At the initial loading stage, the lateral displacement ($\:{y}_{0}$) predicted by the present approach is slightly smaller than that measured in the centrifuge test. As the lateral load increases, however, the calculated displacement becomes large than the experimental value. This discrepancy can be attributed to the mobilization of pile-soil interface friction with increasing load, which engages more surrounding soil in resisting lateral deformation and thus enhances the overall bearing capacity. Since the present approach does not account for interface friction, it tends to overestimate the lateral displacement under higher load levels.

Fig. 8

Comparison of $\:{H}_{0}$-$\:{y}_{0}$ curve for the P334 test.

Figure 9 compares the distributions of lateral displacement ($\:y$) and bending moment ($\:{M}_{a}$) along the depth ($\:z$) between the results from the present study and the centrifuge test, showing generally good agreement. It should be noted that the pile-head lateral displacements under different lateral load levels in Fig. 9 differ from those in Fig. 8. This discrepancy arises because the load-displacement curve in Fig. 8 was measured at a prototype elevation of 1.6 m above the ground surface, whereas the curves in Fig. 9 correspond solely to the embedded portion below the ground surface. As the lateral load increases, the depth ($\:z$) at which lateral displacement becomes negligible gradually increases, as shown in Fig. 9(a). Moreover, the location of the maximum bending moment shifts slightly downward along the pile, as shown in Fig. 9(b).

Fig. 9

Comparison of $\:y$-$\:z$ and $\:{M}_{a}$-$\:z$ curves for the P36 test.

Conclusion

In light of the empirical and semi-empirical correlations commonly used in existing research to predict the lateral load capacity of piles, this paper proposes a novel and rigorous approach for predicting load-displacement responses of laterally loaded piles in sands, incorporating the state-dependent Mohr-Coulomb (SDMC) model. The soil region around the pile is idealized as a cylindrical cavity, and the stress-strain response of the surrounding soil during monotonic lateral loading of pile is dealt with the cavity expansion theory. The governing equations for the soil response are derived by integrating the stress equilibrium equation, displacement compatibility condition, and boundary conditions. To solve the partial differential equations (PDEs), the hybrid Eulerian-Lagrangian (HEL) approach is adopted. The lateral soil resistance is determined by integrating the radial stresses around the pile. The pile-soil interaction response is then solved using the pile deflection equilibrium equation in conjunction with the finite-difference method (FDM). To validate the proposed approach, a finite-element simulation was performed using a user-modified Mohr-Coulomb model. The results show good agreements between the radial and hoop stresses predicted by the present method and those from the numerical simulation. Additionally, parameters corresponding to Fontainebleau sand were applied to predict the lateral load capacity of a pile tested in a centrifuge experiment, further confirming the accuracy of the proposed approach.

The method is particularly suited for long flexible piles, such as those used in offshore wind turbines, drilling platforms, and bridge foundations. For such piles, soil deformation under lateral loading is dominated by flow mechanisms consistent with plane strain conditions, making cavity expansion theory applicable. However, for large-diameter piles, interface friction and rotational mechanisms become significant and cannot be neglected. Future research should therefore focus on extending the proposed approach to account for pile-soil interface friction and related effects.