Nonlinear finite element Analysis of laterally loaded piles in Layered Soils

This work is based on Winkler’s theory to analyze laterally loaded single piles using the finite element method. Soil nonlinearity is considered via nonlinear springs “p - y” curves. Exact element displacement shape functions and stiffness matrix are used f or the element in the case of linear Winkler’s modulus of subgrade reaction. In the nonlinear stage, an averaging technique for the element secant Winkler modulus is used to calculate the shape functions and stiffness matrix. An iterative technique is used to consider the nonlinearity of the soil. Few elements are required to simulate the pile efficiently. Unlike other analytical methods, the current method can be used to analyze piles with any load-transfer curves with arbitrary variation with depth.


Introduction
Piles are increasingly used to support axial and lateral loads in many structures such as tall buildings, earth retaining structures, and bridge substructures. The important aspects in the design of laterally loaded piles are the amount of deflection and the moment in the pile [Poulos and Davis 1980]. Laterally loaded piles can be analyzed by different methods [Liao and Lin, 2003;Zhang et al., 2013;Gupta and Basu, 2019;Franklin and Scott, 1979;Basu et al., 2009;Budhu and Davies, 1987;Sun, 2011]. Some methods treat the soil as an elastic continuum [Basu and Prezzi, 2009;Budhu and Davies, 1987;Sun, 2011]; which requires significant modeling and computation effort. These methods are used primarily for research and not for design purposes [Sun, 2011]. One of the most suitable and popular methods is known as the Winkler foundation method in which the supporting soil is represented by closely spaced springs. For the nonlinear analysis, these springs are assumed to have a nonlinear force-displacement relationship [Reese and Van, 2011]. These relationships are sometimes called p-y curves. The p-y method is very versatile and can be used to analyze and design a variety of cases encountered in practice. The main advantage of the p-y analysis method of laterally loaded piles is that it does not require discretizing the supporting soil into elements. Therefore, the API standard and many commercial software packages such as SAP, pileLAT, GEO5, LAP and many others use this method to analyze and design laterally loaded piles. Therefore, many researchers [Amirmojahedi et al., 2023;Wei et al., 2023;Franke, and Rollins, 2013;Lianyang, 1992;Xiaoling et al, 2020] focused on developing p-y curves for different types of soils. Psaroudakis et al. in 2021 developed a semi-analytical method based on Winkler's theory to analyze laterally loaded single piles. Shape functions that describe the elastic displacements of the pile under loading were used in combination with the principle of virtual work in the analysis. Soil nonlinearity was simulated via nonlinear ''p-y'' springs located along the pile axis. Using the finite difference method, Yin et al. in 2018 proposed a simplified iterative method based on p-y curves to analyze laterally loaded single piles in or near a sloping ground. Basu et al. in 2011 developed an analytical method to analyze piles subjected to horizontal loads and moments at the pile head. The soil supporting the pile was considered a multi-layered and elastic continuum. The energy principles were used in deriving the governing differential equations. An iterative method based on the finite difference method was used to solve the differential equations and determine the displacements and forces along the pile. This study makes use of the p-y curves with arbitrary shapes to efficiently analyze piles using the finite element method with very few elements.

Problem Description and Modeling of the Pile-soil System
The pile model under consideration has length L with flexural rigidity EI. The pile material is considered linear elastic. The pile is considered embedded in multi-layered soil of random mechanical properties. Each soil layer may have a nonlinear reaction modulus (k) and varying ultimate lateral resistance along with its depth. The reaction modulus (k) is defined as the soil reaction per unit length of the pile due to unit lateral displacement of the pile at that location. A force P0 and moment M0 are applied at the pile head.

Homogeneous elastic soil
For the case of homogeneous elastic soil, the reaction modulus k, sometimes referred to as the modulus of subgrade reaction is constant. Using Winkler's hypotheses, the governing differential equation for bending of laterally loaded piles in homogeneous elastic soil is the same as the equation of beam on elastic foundation [Reese and Van, 2011]:

+ − sin
(4) Where A, B, C and D are arbitrary constants and can be determined from the boundary conditions of the pile element. The slope at any point within the pile element is the derivative of the general solution of the differential equation (Eq.4). Therefore, To derive the element shape functions, the following four boundary conditions are set for the element: 1: Where 1 and 1 represent displacement and rotation respectively at node 1 of the element. Similarly, 2 and 2 represent displacement and rotation respectively at node 2 of the element.
Substituting these four boundary conditions into equations 4 and 5, the arbitrary constants 1 − 4 can be determined by solving the resulting four simultaneous equations. And Eq. 4 can be written in terms of nodal displacements as: Where N1, N2, N3 and N4 are the element displacement shape functions and are given by: Where = and = cosh(2 ) − 2 2 − 1.
These shape functions depend on the value of λ and reduce to the Lagrange interpolation polynomials when λ approaches zero. Figure 2 shows the shape functions for different λ values.

Fig. 2 Shape functions for beam on elastic foundation.
The element stiffness matrix is derived by applying the principle of minimum potential energy and may be written as follows.
The terms in the stiffness matrix [K] corresponding to the degrees of freedom shown in Fig. 1  When the value of λ approaches zero, the stiffness matrix reduces to the ordinary beam element stiffness matrix.

Nonlinear layered soil
Assuming the supporting soil as a linear elastic material is an oversimplification of reality. Soil behaves non-linearly in the range of loading of practical interest [Reese, 1977;Davies and Budhu 1986]. Moreover, the hardness of soil profiles changes with depth due to soil stratification. Therefore, the value of soil reaction ( ) is not constant along the element length. Therefore, the differential equation of the laterally loaded pile element will be as follows: Where ( , ) represents the variation of the reaction modulus with depth and displacements. The exact solution for Eq.10 is available for special cases in which the reaction modulus is elastic and varies linearly with depth. Psaroudakis et al. in 2021 proposed a solution for an elastic and parabolic profile of reaction modulus. However, for nonlinear analysis and randomly varying soil profiles, some approximation must be done to solve the governing differential equation.
In this study, an average reaction modulus is calculated for each element as follows: Where ( , ) is the secant soil reaction modulus which is presented graphically in Figure 3.

Fig. 3 Secant reaction modulus
Therefore, the average secant reaction modulus can be written as: In which represents the total lateral soil reaction within the element divided by the total element lateral displacements. The absolute value notation is to take account of displacements in the negative direction. The value of ks depends on the p-y curve of the soil layers, and explicit integration of the numerator may not be possible for many p-y models. Therefore, numerical integration is used to evaluate integrals in equation 13. Using Gaussian quadrature rule, Eq. 13 is evaluated as: where , are lateral displacement, secant soil reaction modulus, and Gaussian weight at Gauss point on the element respectively. The number of Gaussian points (n) used in this study is five. Due to soil nonlinearity, an iterative procedure is used to determine and the related lateral displacements. A flow chart of the solution procedure is shown in Figure 4.

Fig. 4 Flow chart of the proposed solution algorithm
The iterative solution algorithm is shown graphically in Fig. 5.

Fig. 5 Schematic diagram of iterative secant stiffness method.
The proposed method can analyze piles in multilayered soils with different p-y curves. A computer program was developed to implement the proposed analysis method. Only a few elements are required for good accuracy.

Validation of the of the proposed method
To check the feasibility of the proposed analysis method, two examples were solved and compared with others well documented in the literature.

Example 1
A case history analyzed by Reese [Reese 1977] is reanalyzed and the results are compared with those of Reese. The geometrical properties and loading of the pile are shown in Fig. 6. The load transfer "p-y "curves of the soil are shown in Fig. 7 in which the soil properties vary with depth. Linear variation is assumed between any two curves. The soil properties are constant for a depth of more than 3.91 m. This example was reanalyzed using three elements only with five Gauss points along each element. The lateral deflection of the pile is calculated using displacement shape functions of each element and is shown in Fig. 8. On the same figure is shown the results of [Reese, 1977] using the finite difference method obtained by dividing the pile into 72 segments. The bending moment diagram of the pile is shown in Fig. 9. The results are very close to those of Reese by using only three elements.

Example 2
Reese and William in 2011 analyzed a laterally loaded steel pipe pile driven into stiff clay. The p-y curves for the clay layers are shown in Fig.  10. The experimental p-y curves in Fig. 10 were derived from data from the tests at Manor-Texas [Reese, and Robert 1975]. The soil properties were assumed constant below a depth of 3.05m below the ground surface. Linear variation is assumed between any two curves. The pile consists of two pieces with a total length of 15.24 m. The length of the upper piece is 7.01 mm, and its outside diameter is 641 mm. The bottom piece has a length of 8.23 m and an outer diameter of 610 mm. The flexural rigidity (EI) of the two pieces is 493,700 and 168,400 kN.m 2 respectively. The load was applied 0.305m above the ground surface.

Fig. 10 p-y curves of supporting soil of example 2 [Reese and William 2011]
This example was reanalyzed using four elements with only five Gauss points along each element. The load-deflection relationship of the pile head is shown in Fig. 11. In the same Figure, the experimental results and analytical solution of [Reese and William, 2011] are shown for comparison.

Fig. 11 Load-Deflection curve of the pile head
The pile was instrumented with strain gauges along its shaft for measurement of the bending moments. The maximum bending moment of the pile using the present study in addition to the experimental and analytical results of [Reese and William, 2011] are also shown in Figure  12. The agreement is excellent among the experimental, analytical [Reese and William, 2011] and the present study for both the loaddeflection and the load-maximum moment on the pile.

Conclusion
A simple and efficient technique is presented for the nonlinear analysis of laterally loaded piles in multilayered soils using p-y curves. Averaging technique is used to calculate the soil reaction modulus for each element. The corresponding exact shape functions are used to calculate displacements, soil reaction and stiffness matrix for the elements. An iterative method using the secant stiffness method is used to track the soil nonlinearity. The present method does not necessitate discretizing the pile into many elements. Only two to four elements are sufficient for good accuracy.