# A Numerical Investigation to Determine the p-y Curves of Laterally Loaded Piles

Kexin Yin, Lianghui Li, Eugenia Di Filippo

<figure><img src="image_1.png" /></figure>

## To cite this version:

Kexin Yin, Lianghui Li, Eugenia Di Filippo. A Numerical Investigation to Determine the p-y Curves of Laterally Loaded Piles. Mathematics, 2021, 9 (21), pp.2783. 10.3390/math9212783. hal-05307653

HAL Id: hal-05307653

https://hal.science/hal-05307653v1

Submitted on 10 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Article

# A Numerical Investigation to Determine the *p*–*y* Curves of Laterally Loaded Piles

Kexin Yin<fnref n="1" />,*, Lianghui Li<fnref n="2" />,*, and Eugenia Di Filippo<fnref n="1" />,<fnref n="3" />

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Yin, K.; Li, L.; Di Filippo, E. A Numerical Investigation to Determine the *p*–*y* Curves of Laterally Loaded Piles. Mathematics **2021**, *9*, 2783. https://doi.org/10.3390/math9212783

Academic Editors: Ioannis K. Argyros and Clemente Cesarano

Received: 10 August 2021
Accepted: 29 October 2021
Published: 2 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** This paper focuses on a numerical approach to finding the *p*–*y* curves for laterally loaded piles. The Drucker–Prager plastic model is employed and implemented within a finite element MATLAB code. The pre- and post-processing code for Gmsh and related numerical tools are established as well. The *p*–*y* curve results from this new approach have been validated and compared to the typical design equations of API (American Petroleum Institute) and Matlock. The validation reveals that the code leads to lower *p*–*y* curves than the API and Matlock equations when the horizontal displacement is less than 0.35 times the diameter of the pile (B). A sensitivity analysis of the number of elements and the interface thickness is presented. The results indicate that the obtained *p*–*y* curves are independent of the two factors. Finally, the influence of clay content on the *p*–*y* behavior is investigated by the implemented MATLAB code. When *y* < 0.15B, the same lateral capacity values are resulted at clay contents of 27.5% and 55%, and they are higher than the ones for 0% clay content. The *p*–*y* curves show a decreasing trend with increasing clay content after *y* > 0.15B.

**Keywords:** laterally loaded pile; *p*–*y* curve; soil-pile interface; Drucker–Prager model

## 1. Introduction

Pile foundation is one of the most commonly used foundation types in complicated site conditions. Piles in port engineering, offshore oil platforms, and wind turbines are highly subjected to lateral loads. At the engineering scale, one of the main aspects that we want to focus on is the response of piles when they are subjected to lateral loading. This problem can be typically investigated by tracing the so-called *p*–*y* curves, which represent the relationship between the lateral displacement (*y*) and the lateral force (*p*) for a generic transversal section of the pile. The *p*–*y* approach has been widely employed to design laterally loaded piles [1–6]. The development of computer technology has made it possible to study complicated engineering problems by using numerical methods [4–11]. It is more flexible and less time consuming to analyze the *p*–*y* characteristics by using the finite element method than the full-scale in-situ tests.

This paper aims to construct *p*–*y* curves that take into account the properties of the interface between the soil and the concrete pile as characterized via the experimental campaign of interface direct shear test in the laboratory. To do this, the Drucker–Prager constitutive model is adopted to describe irreversible (plastic) strains occurring at the interface between the surrounding soil and the concrete pile. The Drucker–Prager constitutive law is then implemented into a two-dimensional finite element (FE) code within a MATLAB environment that was initially developed by Bonnet and Frangi [12] and Bonnet et al. [13]. A 2-D problem describing the effect of a transversal displacement applied to a pile section

<footnote n="1">¹ Institut de Recherche en Génie Civil et Mécanique (GeM), Ecole Centrale de Nantes, UMR 6183 CNRS, 1 Rue de la Noë, CEDEX 3, 44321 Nantes, France; eugenia-difilippo@libero.it</footnote>

<footnote n="2">² School of Energy and Mining Engineering, China University of Mining and Technology-Beijing, Beijing 100083, China</footnote>

<footnote n="3">³ Department of Civil and Environmental Engineering, University of Perugia, via G. Duranti, 93-06125 Perugia, Italy</footnote>

<footnote>* Correspondence: kexin.yin.research@gmail.com (K.Y.); lilianghui@cumt.cn (L.L.)</footnote>
has been simulated. The resultant transversal force acting on the pile is found to be a function of the applied displacement. The obtained p-y curves are validated and compared with those in the literature relative to different design methods. A sensitivity analysis with respect to the number of elements and the thickness of the interface is also presented to check whether they affect the p-y results. Finally, since there is a lack of literature on the p-y response of piles in clay-rich soils, the p-y curves of a concrete pile in soil with different percentages of clay are provided to see the effect of the varying amounts of clay present in the soil.

## 2. p-y Method and Code Development

### 2.1. p-y Method

When a pile is subjected to horizontal loading, the lateral bearing capacity can be expressed as [14–16]:

$$p_u = N_p s_u B \quad (1)$$

where $p_u$ is the lateral bearing capacity, $N_p$ is the lateral bearing capacity factor, $s_u$ is the undrained shear strength of the soil, and $B$ is the pile diameter. The soil shear strength $s_u$ depends on soil properties, and $N_p$ is related to the failure mechanism of the soil-pile interface. Equation (1) can be used to normalize the p-y curves in the next section.

The concept of the p-y curve was firstly proposed by McClelland and Focht [17], and developed by others [18]. The p-y model, also known as the non-linear Winkler spring model, is commonly used in the design of piles under lateral loading because of its simplicity and low computational cost [5,6,19,20]. In the p-y model, the pile is simplified as a beam while the soil-pile interface is treated as a set of 1-D, non-linear springs [2,19–22]. As mentioned above, $p$ is the resistance force of the soil per unit of pile length, and $y$ is the local pile deflection caused by the horizontal loading [14,20,23]. Since the p-y concept was introduced, several p-y curves were proposed based on different influence factors by previous researchers. For example, the Matlock p-y curve, API p-y curve, Jeanjean p-y curve, and Zhang p-y curve [14,16,24,25]. Three typical p-y curves are presented in the following.

Matlock [16] proposed a p-y curve for a pile located in soft clay as a power relationship between the lateral resistance and the normalized lateral displacement:

$$p = 0.5 p_u \left(\frac{y}{y_c}\right)^{0.33} \quad (2)$$

$y_c$ can be computed by:

$$y_c = 2.5 \varepsilon_{50} B \quad (3)$$

where $\varepsilon_{50}$ is the axial strain corresponding to the 50% maximum principal stress difference in an undrained compression test and $B$ is the diameter of the pile.

Jeanjean et al. [14] proposed a p-y curve function of the initial shear $G_{max}$ and the undrained shear strength of the soil $s_u$:

$$p = p_u \tanh \left[ \frac{G_{max}}{100 s_u} \left( \frac{y}{B} \right)^{0.5} \right] \quad (4)$$

The p-y method is also adopted by the American Petroleum Institute (i.e., API) for a lateral loading pile design. According to API [25], the API p-y relationship is given as:

$$p = F p_u \tanh \left( \frac{k H y}{F p_u} \right) \quad (5)$$

where $k$ is the initial modulus of subgrade reaction (kN/m³), $H$ is the depth (m), and $F$ is a factor to account for cyclic or monotonic static loading conditions It can be evaluated by:

$F = 0.9$ for cyclic loading;

$F = (3.0 - 0.8 \frac{H}{B}) \ge 0.9$ for monotonic loading.
In what follows, a finite element-based p-y curve deduction is presented, capable of accounting for the plastic behavior of the soil and of the soil-concrete interface. The results will be compared with the curves from the typical equations.

## 2.2. Drucker-Prager Model

The Drucker-Prager (DP) model [26] for the soil is implemented using the following MATLAB code, which performs the numerical simulation for calculating the p-y curves. The DP yield criterion is defined as [27]:

$$f = \sqrt{\frac{2}{3}}q - (A - Bp_1) = 0 \qquad (6)$$

where $p_1$ is the mean normal stress, $\sigma_1$ is the major principal stress, $\sigma_2$ is the principal stress, and $\sigma_3$ is the minor principal stress. A and B are related to the cohesion c and friction angle $\varphi$ in the Mohr-Coulomb (MC) criterion; q is the second invariant of the deviatoric stress [27]:

$$q = \sqrt{3J_2} = \frac{1}{\sqrt{2}}[(\sigma_1 - \sigma_2)^2 + (\sigma_2 - \sigma_3)^2 + (\sigma_1 - \sigma_3)^2] \qquad (7)$$

As the Drucker-Prager model is represented in the stress plane by a cone with a circular cross section, it is possible to relate it to the Mohr-Coulomb model considering the section of the MC pyramid to be inscribed or circumscribed to the DP cone. If the circle of the DP cone passes through the tensile corners of the MC yield surface (the MC criterion is inscribed):

$$A = \frac{2\sqrt{6c} \cos \varphi}{3 + \sin \varphi}, \quad B = \frac{2\sqrt{6} \sin \varphi}{3 + \sin \varphi} \qquad (8)$$

the cohesion (c) and the friction angle (φ) can be obtained by direct shear tests [28–34]. If the circle of the DP cone passes through the compression corners of the MC yield surface (MC criterion is circumscribed):

$$A = \frac{2\sqrt{6c} \cos \varphi}{3 - \sin \varphi}, \quad B = \frac{2\sqrt{6} \sin \varphi}{3 - \sin \varphi} \qquad (9)$$

Furthermore, we can assume a plastic potential function of the form [27]:

$$g = \sqrt{\frac{2}{3}}q + bp_1 \qquad (10)$$

where b is the constant plastic dilatancy parameter. Using a return-mapping algorithm in the stress invariant space, the discrete plastic multiplier is given as:

$$\Delta\lambda = \frac{\sqrt{2/3}q^{tr} - (A - Bp_1^{tr})}{2\mu + KBb} \qquad (11)$$

where K and μ are the elastic bulk and shear moduli, respectively.

## 2.3. MATLAB Code

The MATLAB FE code was originally developed by Bonnet and Frangi [12]. Further improvements have been implemented in GeM (École Centrale de Nantes). The Drucker-Prager constitutive model has been implemented within the MATLAB FE code that considers the non-associate and associate flow rule. Although the code started from the relatively simple DP model, more complicated constitutive laws can be developed based on this primary work. The pre- and post-processing phases are performed with the freely available code Gmsh [35]. The link between the MATLAB FE code and the pre- and post-processing code for Gmsh, and several specific numerical tools, has been established as well.
### 3. 2-D p-y Curve Modeling

#### 3.1. Geometry, Mesh, and Boundary Conditions

A 2-D circular geometry of the finite element model is presented in Figure 1. Due to the symmetry, only half of the domain is considered; the geometry consists of the soil, the pile, and the soil-pile interface. The pile is a rigid disc with diameter $B = 1$ m. The soil boundary at the bottom is made up with two parts in the horizontal direction, each part having a width $W = 11.5B$ (11.5 m) from the pile surface point; therefore, the radius (the center of the circle is the pile center) of the largest half circle is $H = 12B = 12$ m, see Figure 1. The size of the circular domain is chosen to be large enough to avoid parasitic effects of the lateral boundaries. The outer boundary of the soil is fixed both in the $x$ and $y$ directions (Figure 1), while the bottom is fixed only in the vertical direction (i.e., $y$-axis in Figure 1). A horizontal displacement is applied to the pile section in the $x$ direction in Figure 1.

<figure><img src="image_5.png" /><figcaption>Figure 1. The circular geometry and boundary conditions.</figcaption></figure>

The interface thickness in the geometry can be referred to as the typical experimental value, e.g., 5 ~ 14$d_{50}$. In previous modeling studies [36,37], it was set as 0.2 of the pile diameter $B$. In the following, a sensitivity analysis with respect to the interface thickness is provided, which ranges from 5 mm (which considers 20$d_{50}$ the of Fontainebleau sand) to 5 cm. Constant-strain, three-node triangle elements (T3) are used in the mesh, while the mesh is refined near the interface (see Figure 2).

#### 3.2. Initial Stress State

This starting numerical study works on a monophasic continuum, so it does not consider the presence of water, and this is the reason why effective stress is not introduced in the following work. Obviously as there is no water considered, the earth pressure can be expressed by total stresses rather than effective ones. The initial stress is generated with the $K_0$ condition according to the following relationship (Figure 3):

$$ \begin{aligned} \sigma_v &= \sigma_{zz} = \gamma z \\ \sigma_n &= \sigma_{xx} = \sigma_{yy} = K_0\sigma_{zz} = K_0\gamma z \end{aligned} \qquad (12) $$

in which $\gamma$ is the unit weight of the soil, $z$ is the depth from the pile section to the ground surface (Figure 3), and $K_0$ is the earth pressure coefficient. For the $p-y$ modeling, a depth of 10 m and $K_0= 0.5$ are considered; the initial stress-state parameters are summarized in Table 1.
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Example of the mesh with T3 elements: (a) the whole domain and (b) a zoom near the soil-pile interface zone.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Initial stress state.</figcaption></figure>

<table><caption>**Table 1.** Initial stress-state parameters.</caption><thead><tr><th>γ (kN/m³)</th><th>z (m)</th><th>K₀ (-)</th><th>σ<sub>zz</sub> (kPa)</th><th>σ<sub>xx</sub> (kPa)</th></tr></thead><tbody><tr><td>20</td><td>10</td><td>0.5</td><td>200</td><td>100</td></tr></tbody></table>

### 3.3. Validation

To validate the calculation of the *p-y* curves, interface test results of soil consisting of 55% clay and concrete plate are selected as input parameters and the corresponding *p-y* curves are traced and compared with different design methods (API and Matlock).
A total of 11,390 T3 elements are used, and the interface width is equal to 5 mm. The properties of the interface, the soil, and the pile are given in Table 2. The concrete pile is considered a rigid material with a Young's modulus of 33 GPa (Table 2). The pile is loaded by stepwise loading (horizontal displacement) with increments of 0.002 m in the MATLAB code. The lateral capacity (p) of the pile is calculated by integrating the normal and tangential stresses along the circumference of the pile (or of the soil-pile interface).

<table><caption>**Table 2.** Parameters for validation that come from soil with 55% clay content.</caption><thead><tr><th></th><th>Young's Modulus E (kPa)</th><th>Poisson's Ratio v (-)</th><th>Cohesion c (kPa)</th><th>Friction Angle φ (°)</th></tr></thead><tbody><tr><td>Interface</td><td>6000</td><td>0.3</td><td>3.31</td><td>19.33</td></tr><tr><td>Soil</td><td>6000</td><td>0.3</td><td>10</td><td>21.81</td></tr><tr><td>Pile</td><td>33 × 10<sup>6</sup></td><td>0.2</td><td>-</td><td>-</td></tr></tbody></table>

The results are presented in terms of normalized *p-y* curves (Figure 4) using the shear strength of the soil and the pile diameter (see also [20,37]). More specifically, in this study, the soil is isotropic normally consolidated and the interface behavior is characterized by a direct shear test in drained condition, hence the critical drained interface shear strength (*s*<sub>d</sub>) of the interface and the pile diameter *B* are used to normalize *p* so as to obtain the lateral bearing capacity factor *p*/*(s*<sub>d</sub>*B*), i.e., the *N*<sub>p</sub> in Equation (1). The horizontal displacement is normalized by the pile diameter *B*.

The *p-y* curves calculated by the MATLAB code with DP model on the circular geometry are compared to the ones computed from the API and Matlock equations in Figure 4. The *p-y* curves from the DP model are lower than the API curve when the horizontal displacement is less than 0.35*<i>B</i>* (Figure 4). This indicates that at this horizontal displacement range (*y* < 0.35*<i>B</i>*) a smaller reaction force is mobilized by the DP model than the API model. The curve from the DP model is higher than that from the API when horizontal displacement is larger than 0.35*<i>B</i>. From 0 to 0.4*<i>B</i>, the MATLAB curve presents a similar shape to the Matlock curve, but with lower normalized *p* (Figure 4). The two curves become closer with increasing horizontal displacement. At *y* = 0.4*<i>B</i>, the normalized *p* from the DP model in the MATLAB code stays between the API and Matlock values. The normalized *p* calculated from the MATLAB code is 8.58, which is 0.5 higher than the one of API (8.08) and 0.4 lower than the one from the Matlock empirical calculation (8.98), see Figure 4. As the horizontal displacement goes up to 0.45*<i>B</i>, the result of the MATLAB code is equal to, then gradually becomes higher than that of Matlock. Finally, at 0.6*<i>B</i>, the normalized lateral capacity of the MATLAB code is 10.57, which is 104.6% of that of Matlock (10.28).

In the MATLAB code, *p* is calculated by integrating the normal and tangential stresses along the circumference of the pile or the soil-pile interface. The two results are compared in Figure 5, which reveals that the *p-y* curves are the same before 0.4*<i>B</i> and have a small difference (1.38 ~ 2.23%) at 0.4 ~ 0.6*<i>B</i>.
<figure><img src="image_8.png" /><figcaption>**Figure 4.** *p*–*y* curves from the MATLAB code and the design methods.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** *p*–*y* curves integrating the stresses along the pile and the soil–pile interface.</figcaption></figure>

Even if the *p*–*y* curve obtained from the MATLAB code presents a shape which is not identical to those from the API and Matlock empirical equations, the normalized horizontal reaction force does not exceed the empirical values when *y* < 0.35*B*, then it is characterized with smaller difference to the Matlock curve up until *y*/*B* = 0.6. Furthermore, the validation procedure reveals that in engineering design, the DP model leads to lower *p*–*y* curves than the API and Matlock equations at *y* < 0.35*B* (Figure 4). However, at horizontal displacement larger than 0.35*B*, the MATLAB code generates lateral reaction force that is close to that of the Matlock formulation. The integrated results from the pile and soil–pile interface are the same. Both confirm that the *p*–*y* curves from the MATLAB code are reliable and consistent with the Matlock design method.

**4. Sensitivity Analysis**

4.1. *Number of Elements*

In order to check the spatial discretization, the influence of the number of elements is
studied hereafter. The circular geometry with a 5 mm interface thickness is used for the
FE number sensitivity analysis. More specifically 3600, 6300, and 11,000 T3 elements are
considered. The input parameters of the simulations are listed in Table 2.
The normalized *p*–*y* curves are similar when the number of elements is 3600, 6300, and 11,000, see Figure 6. At 0.4B horizontal displacement, the normalized *p* values rank as 8.71, 8.69, and 8.58 for the three element numbers studied (3600, 6300, and 11,000). The *p*–*y* curves are fitted and presented in Table 3, indicating that no significant difference exists for the three element numbers. From a coarse mesh with 3600 elements to a dense mesh with 11,000 elements, the final normalized *p* (at 40% of *y*/*B*) decreases by a small percentage of 1.57%. Considering the necessary computational time, the element number for the *p*–*y* modeling can be chosen as 6000 ~ 10,000, to insure both effective calculation and enough elements on the interface zone.

<figure><img src="image_10.png" /><figcaption>Figure 6. Influence of the number of elements.</figcaption></figure>

<table><caption>**Table 3.** Fitting formulations of the *p*–*y* curves for different numbers of elements.</caption><thead><tr><td>Element Number (-)</td><td>Fitting Formulation</td><td>R<sup>2</sup></td></tr></thead><tbody><tr><td>3600</td><td><sup>p</sup>/<sub>s<sub>d</sub>B</sub> = -0.0058(<sup>y</sup>/<sub>B</sub>)<sup>2</sup> + 0.4337(<sup>y</sup>/<sub>B</sub>)</td><td>0.9524</td></tr><tr><td>6300</td><td><sup>p</sup>/<sub>s<sub>d</sub>B</sub> = -0.0058(<sup>y</sup>/<sub>B</sub>)<sup>2</sup> + 0.4335(<sup>y</sup>/<sub>B</sub>)</td><td>0.9520</td></tr><tr><td>11,000</td><td><sup>p</sup>/<sub>s<sub>d</sub>B</sub> = -0.0058(<sup>y</sup>/<sub>B</sub>)<sup>2</sup> + 0.4280(<sup>y</sup>/<sub>B</sub>)</td><td>0.9510</td></tr></tbody></table>

## 4.2. Interface Thickness

To check the influence of the interface thickness on the numerical results, the following different interface-thickness values are considered: 5 mm, 7 mm, 1 cm, 2.5 cm, 5 cm, and 10 cm. The input parameters of the simulations are listed in Table 2.

Figure 7 presents the effect of interface thickness on the soil reaction curves. All the *p*–*y* curves are independent of the interface thickness at the horizontal displacement considered 0 ~ 0.4B (Figure 7). From 5 mm to 10 cm interface width, the resulting normalized *p*–*y* curves have a close shape and almost the same value of the lateral bearing capacity factor at *y* = 0.4B: 8.58, 8.69, 8.55, 8.75, 8.81, and 8.67. The difference between the final normalized *p* of the six curves from the MATLAB code is 0.02 ~ 0.26. These values are higher than the one calculated by the API method (8.08) and lower than the one calculated by the Matlock equation (8.98), as shown in Figure 7.
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Influence of the interface thickness.</figcaption></figure>

The obtained *p*–*y* curves are therefore independent of the interface thickness, from thickness ranges from 5 mm to 10 cm. This is also confirmed by the fitting formulations of the *p*–*y* curves as presented in Table 4.

<table><caption>**Table 4.** Fitting formulations of the *p*–*y* curves for different interface thickness.</caption><thead><tr><td>Interface Thickness</td><td>Fitting Formulation</td><td>R<sup>2</sup></td></tr></thead><tbody><tr><td>5 mm</td><td><sup>p</sup>/<sub>s<sub>d</sub>B</sub> = -0.0058(<sup>y</sup>/<sub>B</sub>)<sup>2</sup> + 0.4280(<sup>y</sup>/<sub>B</sub>)</td><td>0.9510</td></tr><tr><td>7 mm</td><td><sup>p</sup>/<sub>s<sub>d</sub>B</sub> = -0.0058(<sup>y</sup>/<sub>B</sub>)<sup>2</sup> + 0.4338(<sup>y</sup>/<sub>B</sub>)</td><td>0.9525</td></tr><tr><td>1 cm</td><td><sup>p</sup>/<sub>s<sub>d</sub>B</sub> = -0.0056(<sup>y</sup>/<sub>B</sub>)<sup>2</sup> + 0.4227(<sup>y</sup>/<sub>B</sub>)</td><td>0.9508</td></tr><tr><td>2.5 cm</td><td><sup>p</sup>/<sub>s<sub>d</sub>B</sub> = -0.0059(<sup>y</sup>/<sub>B</sub>)<sup>2</sup> + 0.4393(<sup>y</sup>/<sub>B</sub>)</td><td>0.9360</td></tr><tr><td>5 cm</td><td><sup>p</sup>/<sub>s<sub>d</sub>B</sub> = -0.0057(<sup>y</sup>/<sub>B</sub>)<sup>2</sup> + 0.4335(<sup>y</sup>/<sub>B</sub>)</td><td>0.9528</td></tr><tr><td>10 cm</td><td><sup>p</sup>/<sub>s<sub>d</sub>B</sub> = -0.0058(<sup>y</sup>/<sub>B</sub>)<sup>2</sup> + 0.4316(<sup>y</sup>/<sub>B</sub>)</td><td>0.9481</td></tr></tbody></table>

In the following, a circular mesh with 6000 T3 elements and an interface thickness of 5 mm is considered. This thickness also corresponds to the experimental laboratory results.

### 5. Effect of Clay Fraction on *p*–*y* Curves

The adhesion and friction angles from the interface direct shear results are used as input parameters for the *p*–*y* modeling to investigate how the clay fraction affects the curves. The influence of the clay content on the *p*–*y* behavior of pile with lateral loading is presented and discussed below.

The adhesion and friction angles from three clay fractions of 0%, 27.5%, and 55% are chosen as the input parameters for the *p*–*y* calculation (Tables 5 and 6). The parameters of the pile used are the same, as shown in Table 2. The horizontal loading is applied on the pile section considering 200 steps of 0.002 m each (i.e., in total 0.4B horizontal displacement). The residual drained shear strength, $s_d$, from the interface tests under 100 kPa normal stress and the pile diameter B are used to do the normalization.
<table><caption>**Table 5.** Parameters of the interface.</caption><thead><tr><th>Clay Fraction (%)</th><th>Young's Modulus E (kPa)</th><th>Poisson's Ratio v (-)</th><th>Cohesion c (kPa)</th><th>Friction Angle δ (°)</th></tr></thead><tbody><tr><td>0</td><td>8000</td><td>0.3</td><td>0</td><td>26.40</td></tr><tr><td>27.5</td><td>7500</td><td>0.3</td><td>2.22</td><td>22.91</td></tr><tr><td>55</td><td>6000</td><td>0.3</td><td>3.31</td><td>19.33</td></tr></tbody></table>

<table><caption>**Table 6.** Parameters of the soil.</caption><thead><tr><th>Clay Fraction (%)</th><th>Young's Modulus E (kPa)</th><th>Poisson's Ratio v (-)</th><th>Cohesion c (kPa)</th><th>Friction Angle φ (°)</th></tr></thead><tbody><tr><td>0</td><td>8000</td><td>0.3</td><td>0</td><td>34.11</td></tr><tr><td>27.5</td><td>7500</td><td>0.3</td><td>2.22</td><td>25</td></tr><tr><td>55</td><td>6000</td><td>0.3</td><td>10</td><td>21.81</td></tr></tbody></table>

The effect of clay content in the soil on the *p*-y curves is illustrated in Figure 8 and the reformulated functions are provided in Table 7. The curves show non-linearity shapes due to the non-linearity characteristics of the soil-pile interface. The overall shapes of the *p*-y curves in Figure 8 do not exhibit a strong asymptotic behavior in the horizontal displacement range considered (0 ~ 40 cm), which agrees with the results calculated using the Mohr-Coulomb model in [2] and [36]. When *y* < 0.15B, the same lateral capacity values are resulted at clay contents of 27.5% and 55%, and they are higher than the curve of the sand (Figure 8), but the differences are not significant. At the 0.15B, the *p* is about 5.22 for the three clay contents presented. To conclude, clay content has nearly no effect on the *p*-y results when *y* is smaller than 0.15B.

<figure><img src="image_12.png" /><figcaption>**Figure 8.** *p*-y curves as a function of clay content.</figcaption></figure>

<table><caption>**Table 7.** Fitting formulations of the *p*-y curves at different clay contents.</caption><thead><tr><th>Clay Content (%)</th><th>Fitting Formulation</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td>0</td><td>$\frac{p}{s_d B} = -0.0042(\frac{y}{B})^2 + 0.4109(\frac{y}{B})$</td><td>0.9839</td></tr><tr><td>27.5</td><td>$\frac{p}{s_d B} = -0.0053(\frac{y}{B})^2 + 0.4249(\frac{y}{B})$</td><td>0.9668</td></tr><tr><td>55</td><td>$\frac{p}{s_d B} = -0.0058(\frac{y}{B})^2 + 0.4280(\frac{y}{B})$</td><td>0.9510</td></tr></tbody></table>
It follows from Figure 8 that the normalized *p* after *y* > 0.15*B* is influenced by the clay content. The curve for 0% clay content (i.e., sand) becomes higher than the one for 27.5% clay content at each load step, see Figure 8. The curve of 27.5% clay content is also higher than the one for 55% clay content but with a smaller *p*/(*s*<sub>d</sub>*B*) gap (0.35) at each loading step compared to the one (0.46) for between 0% and 27.5% clay content. At a horizontal displacement of 40 cm (0.4*B*), the normalized *p* attains 10.18, 9.23, 8.58 for 0%, 27.5%, and 55% clay contents, respectively (Figure 8); this is because the adhesion, friction angles, and Young's modulus vary with the increasing clay content (Tables 5 and 6). The *p*-*y* curves are affected by the cohesion and friction angles, which are consistent with the literature [5,6,37].

**6. Conclusions**

This paper focuses on a numerical approach to finding the *p*-*y* curves for laterally
loaded piles. The Drucker-Prager plastic model has been employed and the approach
has been validated and compared to the API and Matlock design equations. A sensitivity
analysis in terms of the number of elements and interface thickness has been presented.
Finally, the influence of the clay content on the *p*-*y* behavior is presented and discussed.
The main conclusions of this paper are summarized as follows:

1. The *p*-*y* results agree with the empirical results of Matlock. The validation reveals that the DP model leads to lower *p*-*y* curves with respect to those from Matlock and API when the horizontal displacement is less than 0.35*B*.

2. The number of elements has no important effect on the *p*-*y* curves. Considering the necessary computational time, the *p*-*y* modeling with the MATLAB code is effective when 6000 ~ 10,000 elements are adopted for the spatial discretization.

3. The *p*-*y* curves are independent of interface thicknesses from 5 mm to 10 cm.

4. Clay content influences the *p*-*y* curve results. When *y* < 0.15*B*, the same lateral capacity values are resulted at clay contents of 27.5% and 55%, and they are higher than the curve of the sand. The normalized *p*-*y* curves show a decreasing trend with increasing clay content after *y* > 0.15*B*.

The primary results in this paper shed light on the relationship between *p*-*y* curves and interface parameters. However, the FE modeling of *p*-*y* curves in this study is performed by an implemented code using only the DP model, and the input parameters concern just three clay contents. Therefore, further numerical studies should involve more advanced constitutive laws, water, 3-D conditions, and the possibility to have gaps between the pile and the soil. A new function relating the lateral capacity and the horizontal displacement, as well as the clay content, should be proposed.

**Author Contributions:** Conceptualization, methodology, data curation, writing-original draft preparation, K.Y.; formal analysis, K.Y.; writing-review and editing, K.Y., L.L. and E.D.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the
corresponding author.

**Acknowledgments:** Thanks for the editors' work and reviewers' comments.

**Conflicts of Interest:** The authors declare that they have no conflict of interest.

**References**

1. Zhang, Y.; Andersen, K.H.; Jeanjean, P. Cyclic py curves in clays for offshore structures. In Proceedings of the Offshore Technology Conference, Houston, TX, USA, 6–9 May 2019.

2. Wolf, T.K.; Rasmussen, K.L.; Hansen, M.; Ibsen, L.; Roesen, H. *Assessment of Py Curves from Numerical Methods for a Non-Slender Monopile in Cohesionless Soil*; Department of Civil Engineering, Aalborg University: Aalborg, Denmark, 2013.
3. Dash, S.; Rouholamin, M.; Lombardi, D.; Bhattacharya, S. A practical method for construction of py curves for liquefiable soils. *Soil Dyn. Earthq. Eng.* **2017**, *97*, 478–481. [CrossRef]

4. Liu, X.; Cai, G.; Liu, L.; Liu, S.; Duan, W.; Puppala, A.J. Improved py curve models for large diameter and super-long cast-in-place piles using piezocone penetration test data. *Comput. Geotech.* **2021**, *130*, 103911. [CrossRef]

5. Yang, Z.; Jeremić, B. Numerical analysis of pile behaviour under lateral loads in layered elastic-plastic soils. *Int. J. Numer. Anal. Methods Geomech.* **2002**, *26*, 1385–1406. [CrossRef]

6. Yang, Z.; Jeremić, B. Study of soil layering effects on lateral loading behavior of piles. *J. Geotech. Geoenvironmental Eng.* **2005**, *131*, 762–770. [CrossRef]

7. Wang, S.-F.; Tang, Y.; Wang, S.-Y. Influence of brittleness and confining stress on rock cuttability based on rock indentation tests. *J. Cent. South Univ.* **2021**, *28*, 2786–2800. [CrossRef]

8. Yu, Q.; Yin, K.; Ma, J.; Shimada, H. Vertical Shaft Support Improvement Studies by Strata Grouting at Aquifer Zone. *Adv. Civ. Eng.* **2018**, *2018*, 5365987. [CrossRef]

9. Wang, J.; Yang, S.; Wei, W.; Zhang, J.; Song, Z. Drawing mechanisms for top coal in longwall top coal caving (LTCC): A review of two decades of literature. *Int. J. Coal Sci. Technol.* **2021**, *1*, 26. [CrossRef]

10. Mikkelsen, A.T.; Nielsen, S.D.; Østergaard, M.U. Finite Element Modelling of p-y curves for Monopiles in Liquefied Soil. In Proceedings of the 29th International Ocean and Polar Engineering Conference, Honolulu, HI, USA, 16–21 June 2019.

11. Wang, S.-F.; Yu, T.; Li, X.-B.; Kun, D. Analyses and predictions of rock cuttabilities under different confining stresses and rock properties based on rock indentation tests by conical pick. *Trans. Nonferrous Met. Soc. China* **2021**, *31*, 31–1766. [CrossRef]

12. Bonnet, M.; Frangi, A. Analyse des structures mécaniques par la méthode des éléments finis. *Notes de Cours, Ecole Polytechnique*. 2005. Available online: http://nguyen.hong.hai.free.fr/EBOOKS/SCIENCE%20AND%20ENGINEERING/MATHEMATIQUE/PROGRAMMING/FINITE%20ELEMENTS%20METHODE/Analyse%20m%E9canique%20par%20FEM.Bonnet.pdf (accessed on 23 May 2021).

13. Bonnet, M.; Frangi, A.; Rey, C. *The Finite Element Method in Solid Mechanics*; McGraw Hill Education: New York, NY, USA, 2014; p. 365.

14. Jeanjean, P.; Zhang, Y.; Zakeri, A.; Andersen, K.; Gilbert, R.; Senanayake, A. A framework for monotonic py curves in clays. In Proceedings of the Offshore Site Investigation Geotechnics 8th International Conference Proceeding, London, UK, 12–14 September 2017; pp. 108–141.

15. Martin, C.; Randolph, M. Upper-bound analysis of lateral pile capacity in cohesive soil. *Géotechnique* **2006**, *56*, 141–145. [CrossRef]

16. Matlock, H. Correlations for design of laterally loaded piles in soft clay. In Proceedings of the Offshore Technology in Civil Engineering's Hall of Fame Papers from the Early Years, Houston, TX, USA, 1 April 1970; pp. 77–94.

17. McClelland, B.; Focht, J. Soil modulus for laterally loaded piles. *J. Soil Mech. Found. Div.* **1956**, *82*, 1081-1. [CrossRef]

18. Yang, M.; Ge, B.; Li, W.; Zhu, B. Dimension effect on py model used for design of laterally loaded piles. *Procedia Eng.* **2016**, *143*, 598–606. [CrossRef]

19. Nogami, T.; Otani, J.; Konagai, K.; Chen, H.-L. Nonlinear soil-pile interaction model for dynamic lateral motion. *J. Geotech. Eng.* **1992**, *118*, 89–106. [CrossRef]

20. Ahayan, S. *A constitutive Model for natural Clays: From Laboratory Testing to Modelling of Offshore Monopiles*; École centrale de Nantes, Université de Liège, Faculté des Sciences: Liège, Belgium, 2019.

21. Sjrensen, S.P.H.; Brjdbøk, K.T.; Mjller, M.; Augustesen, A.H. *Review of Laterally Loaded Mono-Piles Employed as the Foundation for Offshore Wind Turbines*; Department of Civil Engineering, Aalborg University: Aalborg, Denmark, 2012.

22. Ahayan, S.; Cerfontaine, B.; Collin, F.; Kotronis, P. Behaviour of laterally loaded pile. In Proceedings of the 9th European Conference on Numerical Methods in Geotechnical Engineering, NUMGE 2018, Porto, Portugal, 25–27 June 2018; pp. 1511–1518.

23. Suryasentana, S.K.; Lehane, B.M. Numerical derivation of CPT-based p-y curves for piles in sand. *Géotechnique* **2014**, *64*, 186–194. [CrossRef]

24. Zhang, Y.; Andersen, K.H. Scaling of lateral pile py response in clay from laboratory stress-strain curves. *Mar. Struct.* **2017**, *53*, 124–135. [CrossRef]

25. API. American petroleum institute recommended practice for planning, designing and constructing fixed offshore platforms-working stress design. In *API RP2A-WSD*; American Petroleum Institute: Washington, DC, USA, 2007.

26. Drucker, D.C.; Prager, W. Soil mechanics and plastic analysis or limit design. *Q. Appl. Math.* **1952**, *10*, 157–165. [CrossRef]

27. Borja, R.I. *Plasticity: Modeling & Computation*; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2013.

28. Yin, K.; Liu, J.; Lin, J.; Vasilescu, A.-R.; Othmani, K.; Di Filippo, E. Interface Direct Shear Tests on JEZ-1 Mars Regolith Simulant. *Appl. Sci.* **2021**, *11*, 7052. [CrossRef]

29. Vasilescu, A.-R. Design and Execution of Energy Piles: Validation by In-Situ and Laboratory Experiments. Ph.D. Theis, École Centrale de Nantes, Nantes, France, 2019.

30. Dafalla, M.A. Effects of clay and moisture content on direct shear tests for clay-sand mixtures. *Adv. Mater. Sci. Eng.* **2013**, *2013*, 562726. [CrossRef]

31. Maghsoodi, S. Thermo-Mechanical Behavior of Soil-Structure Interface under Monotonic and Cyclic Loads in the Context of Energy Geostructures. Ph.D. Theis, Université de Lorraine, Nancy, France, 2020.

32. Maghsoodi, S.; Cuisinier, O.; Masrouri, F. Effect of Temperature on the Cyclic Behavior of Clay-Structure Interface. *J. Geotech. Geoenvironmental Eng.* **2020**, *146*, 04020103. [CrossRef]
33. Konkol, J.; Mikina, K. Some Aspects of Shear Behavior of Soft Soil-Concrete Interfaces and Its Consequences in Pile Shaft Friction Modeling. *Materials* **2021**, *14*, 2578. [CrossRef] [PubMed]

34. Yin, K.; Liu, J.; Vasilescu, A.-R.; Di Filippo, E.; Othmani, K. A Procedure to Prepare Sand-Clay Mixture Samples for Soil-Structure Interface Direct Shear Tests. *Appl. Sci.* **2021**, *11*, 5337. [CrossRef]

35. Geuzaine, C.; Remacle, J.F. Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities. *Int. J. Numer. Methods Eng.* **2009**, *79*, 1309–1331. [CrossRef]

36. Østergaard, M.U.; Knudsen, B.S.; Ibsen, L.B. P-y curves for bucket foundations in sand using finite element modeling. In Proceedings of the 3rd Internationl symposium on Frontiers in offshore Geotechnics, Oslo, Norway, 10–12 June 2015; pp. 343–348.

37. Vethanayagam, V.; Ibsen, L.B. *Determination of py Curves for Bucket Foundations in Silt and Sand Using Finite Element Modelling*; Department of Civil Engineering, Aalborg University: Aalborg, Denmark, 2017.