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

: This paper focuses on a numerical approach to ﬁnding the p–y curves for laterally loaded piles. The Drucker–Prager plastic model is employed and implemented within a ﬁnite 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 inﬂuence of clay content on the p–y behavior is investigated by the implemented MATLAB code. 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 ones for 0% clay content. The p–y curves show a decreasing trend with increasing clay content after y > 0.15 B .


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][2][3][4][5][6]. The development of computer technology has made it possible to study complicated engineering problems by using numerical methods [4][5][6][7][8][9][10][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 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.

p-y Method
When a pile is subjected to horizontal loading, the lateral bearing capacity can be expressed as [14][15][16]: 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][20][21][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: y c can be computed by: where ε 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 : 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: where k is the initial modulus of subgrade reaction (kN/m 3 ), 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 H B ≥ 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.

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]: where p 1 is the mean normal stress, σ 1 is the major principal stress, σ 2 is the principal stress, and σ 3 is the minor principal stress. A and B are related to the cohesion c and friction angle ϕ in the Mohr-Coulomb (MC) criterion; q is the second invariant of the deviatoric stress [27]: 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): the cohesion (c) and the friction angle (ϕ) can be obtained by direct shear tests [28][29][30][31][32][33][34].
If the circle of the DP cone passes through the compression corners of the MC yield surface (MC criterion is circumscribed): Furthermore, we can assume a plastic potential function of the form [27]: 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: where K and µ are the elastic bulk and shear moduli, respectively.

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.

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. The interface thickness in the geometry can be referred to as the typical experimental value, e.g., 5~14d 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 20d 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).

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): (12) in which γ 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.

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). 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 d ) of the interface and the pile diameter B are used to normalize p so as to obtain the lateral bearing capacity factor p/(s d B), i.e., the N p 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.35B (Figure 4). This indicates that at this horizontal displacement range (y < 0.35B) 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.35B. From 0 to 0.4B, 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.4B, 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.45B, the result of the MATLAB code is equal to, then gradually becomes higher than that of Matlock. Finally, at 0.6B, 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.4B and have a small difference (1.38~2.23%) at 0.4~0.6B.  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.35B, 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.35B (Figure 4). However, at horizontal displacement larger than 0.35B, 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.

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.

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. 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 4. Fitting formulations of the p-y curves for different interface thickness.

Interface Thickness
Fitting Formulation 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.

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. 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.  It follows from Figure 8 that the normalized p after y > 0.15B 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 d 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.4B), 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].

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.35B.

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.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. The normalized p-y curves show a decreasing trend with increasing clay content after y > 0.15B.
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.