Predicting Lateral Resistance of Piles in Cohesive Soils

: The ultimate lateral resistance of free- and ﬁxed-headed piles in cohesive soil is examined in this paper using the three-dimensional ﬁnite element limit analysis with upper and lower bound theorems. A special concern, and that is the novelty of this study, is devoted to the combined effect of the three important dimensionless parameters; namely, the overburden stress factor ( n ), the pile length-diameter ratio ( L/D ), and the ratio of eccentric length to diameter ( e/D ). Numerical results are expressed by using Broms’s horizontal load factor, and comparisons are made with several published solutions. In addition, the associated failure mechanisms are investigated with respect to the three parametric effects. The adopted new technique has been successfully used to study a number of different geo-stability problems. It is thus the aim of this paper to produce accurate and practical results with design equations and charts that can be used by practitioners to predict the undrained lateral capacity of ﬁxed- and free-headed piles.


Introduction
High-rise buildings, transmission towers, offshore platforms, bridges, or wind turbines are generally supported by pile foundations [1][2][3][4][5][6][7][8][9]. Owing to several circumstances, such as wind loadings, wave forces, or seismic forces from earthquake actions, piles are designed to resist lateral force induced from those actions. To predict the lateral capacity of piles, reliable design equations or charts are essential for practitioners in the preliminary design stage of pile foundations.
Broms [10] was the first to investigate the lateral capacity of a pile in cohesive soil using the limit equilibrium method, where equilibrium equations were established with simple geometrical earth pressure distributions along pile length. The solution was presented in the form of a normalized horizontal force that is a function of the length-diameter ratio and eccentric-length ratio of a pile. Broms [10] also considered both fixed-headed (restrained) and free-headed (unrestrained) piles in the work. This is widely known as Broms's design charts for evaluating the lateral capacity of a pile. Later, Meyerhof et al. [11] and Georgiadis et al. [12] extended Broms's work by using different earth pressure distributions of soil reaction along the pile length in order to consider the cases of laterally loaded piles in layered soils and sloping ground, respectively.
Izadi and Chenari [25,26] H γLD 2 = f s u γD , L D , α, λ from a design chart where α is adhesion factor and λ is strength gradient ratio To bridge the current research gap, this project considers the combined effect of the soil weight, the no-tension (separation) condition at the interface of piles, the eccentric length of piles, and both the fixed-and free-headed piles on the lateral pile capacity of piles. In this paper, three-dimensional finite element limit analysis (3D FELA) is employed to numerically derive lower bound (LB) and upper bound (UB) solutions of laterally loaded piles in clays. Numerical results are presented in a similar form to Broms's design charts. The derived failure mechanisms of the lateral pile problem are discussed and empirical design equations for predicting the undrained lateral capacity of both fixed-and freeheaded piles are developed. The accurate results presented in this study can be used by practitioners to predict the undrained lateral capacity of fixed-and free-headed piles under horizontal loading considering the effects of overburden stress factors, eccentric length ratios, and pile length ratios.

Method of Analysis
To derive numerical solutions for the undrained lateral capacity of 3D rigid circular piles, the commercial software OptumG3 (Copenhagen, Denmark) [44] is employed. The historical development of 3D FELA can be found in Lyamin and Sloan [45,46], Kranbbenhoft et al. [47], Sloan [48], Keawsawasvong and Ukritchon [49], Ukritchon and Keawsawasvong [50], and Ukritchon et al. [51,52]. Figure 1 shows the problem statement of 3D rigid circular piles using 3D meshes from OptumG3. A free-headed (unrestrained) pile is shown in Figure 1a, whilst Figure 1b shows a fixed-headed (restrained) pile. The pile is subjected to a horizontal load (H). Since the problem of a monotonically and laterally loaded pile is symmetrical, only half of the model is used in all analyses. The rigid circular pile has a diameter (D) and a length (L). To model the interface between piles and soils, zero-thickness interface elements are used. The undrained shear strength at the interface is equal to that of the surrounding soil (i.e., s ui = s u ) which means that the adhesion factor at the interface is equal to one (α = 1). In addition, the tension cut-off condition at the soil-pile interface is adopted, which allows the soil to separate from the back of the circular pile whenever "necessary" during the solution process (see [17,[53][54][55]). Additionally, the tension cut-off condition is implemented at the soil-pile interface, allowing the soil to separate from the rectangular pile's backside since it was required throughout the solution procedure. By using OptumG3, the interface elements are added at the soil-pile interface, where the properties of this interface are the same as the surrounding clay, except the tension stresses are not allowed to take place at the interface. The surrounding clay is defined as homogeneous and isotropic clays with an elastic-perfectly plastic material obeying the Tresca failure criterion. The soil parameters required for the model are the undrained shear strength (s u ) and the soil unit weight (γ). Note the eccentric length (e) in a free-headed pile (see Figure 1a). On the other hand, in Figure 1b, there is no eccentric length for a fixed-headed pile since the top of the pile is restrained against vertical movement and thus pile rotation. Sloan [48] stated that the stability analysis using the FELA based on the limit analysis theory requires only the conventional strength parameters, such as undrained shear strength, but does not use the deformation parameters, such as Poisson's ratio and Young's modulus, which is different from the conventional displacement-based FEM. Thus, the Poisson's ratio and Young's modulus are not considered in this study using FELA. Hence, the displacements of the pile cannot be investigated by using FELA. In all numerical analyses of the paper, the model domains are built to be sufficiently large so that collapsible areas are fully contained, and the failure zone would not reach or intersect the far boundaries of the model. The boundary conditions shown in Figure 2 are described next. At the domain bottom, zero velocity is applied in x-y-z directions so that  In all numerical analyses of the paper, the model domains are built to be sufficiently large so that collapsible areas are fully contained, and the failure zone would not reach or intersect the far boundaries of the model. The boundary conditions shown in Figure 2 are described next. At the domain bottom, zero velocity is applied in x-y-z directions so that there is no movement in the bottom plane in any direction. On all three side-planes, only the velocity in normal directions is fixed so that the roller supports are applied along all side planes. This condition is the same as the symmetrical plane where all nodes are fixed in the normal direction. The other two tangential directions are free to move. On the top plane of soils is a free surface, and the nodes are free to move in all directions. Note that a special boundary condition is given to the fixed-headed pile, where movement in the vertical direction is not permitted at the pile cap. In this way, no rotation can occur to the pile. The sizes of boundaries are selected to be large enough to confirm the failure zones are captured within the boundary domain so that it cannot cause an insignificant effect on the computed FELA solutions. However, the extremely large size of the domain can result in an increase in the number of elements directly reducing the computational performance of numerical simulations (i.e., increasing calculation time). Based on several trial-and-error checks, the optimal sizes of the model are 5D × 10D × 1.4 L for the width, length, and depth of the domain (see Figure 2). In the proposed study, six dimensional parameters of laterally loaded piles are considered; they are (H, n, D, L, Su, and ). These six parameters can be further reduced to the following dimensionless parameters as: For fixed-headed piles: For free-headed piles: where H/(suLD) is the normalized horizontal load factor and it is a function of (L/D, n, and e/D), depending on whether it is a free-or a fixed-headed pile. (L/D) is the ratio of pile length to diameter, (n = γL/su) is the overburden stress factor, and (e/D) is the ratio of eccentric length to pile diameter. The n parameter may be interpreted as the degree of the total overburden pressure of soil at the pile tip as compared to its average undrained shear strength. Conversely, the reciprocal of n may be interpreted as the ratio of the average undrained shear strength of soil normalized by the total overburden pressure. The special condition of n = 0 corresponds to the ideal case of weightless soil (i.e., γ = 0). The ranges of the overburden stress factor n are chosen to be n = 0 to 80 in the study. The selected pile length-diameter ratios are L/D = 5 to 60 and the ratio of eccentric length to pile diameter are e/D = 0 to 16. These dimensionless parameters can be simply used in practice as they were employed by Keawsawasvong [Sauer] to investigate the undrained capacity of loaded circular piles One of OptumG3's advanced features is the automatic adaptive mesh refinement scheme, which is based on the previous development in Ciria et al. [56]. The numbers of meshes in sensitive zones with very high shear power dissipation, which is the control variable for error estimation in all analyses, are automatically increased through successive iterations using the adaptive technique. The original and targeted number of elements used in all numerical runs are 5000 and 10,000 elements with five adaptive iterations [57][58][59][60][61][62]. Utilizing the adaptive mesh techniques, meshes will automatically enlarge in sensitive zones with considerable plastic shearing strain. All numerical simulations employ 5000 to 10,000 elements as the initial and goal number of elements, respectively, with five adaptive iterations. Examples of the final mesh refinements for a free-headed pile and a fixed-headed pile are shown in Figure 2a,b, respectively.
In the proposed study, six dimensional parameters of laterally loaded piles are considered; they are (H, n, D, L, S u , and γ). These six parameters can be further reduced to the following dimensionless parameters as: For fixed-headed piles: For free-headed piles: Sustainability 2022, 14, 12940 where H/(s u LD) is the normalized horizontal load factor and it is a function of (L/D, n, and e/D), depending on whether it is a free-or a fixed-headed pile. (L/D) is the ratio of pile length to diameter, (n = γL/s u ) is the overburden stress factor, and (e/D) is the ratio of eccentric length to pile diameter. The n parameter may be interpreted as the degree of the total overburden pressure of soil at the pile tip as compared to its average undrained shear strength. Conversely, the reciprocal of n may be interpreted as the ratio of the average undrained shear strength of soil normalized by the total overburden pressure. The special condition of n = 0 corresponds to the ideal case of weightless soil (i.e., γ = 0). The ranges of the overburden stress factor n are chosen to be n = 0 to 80 in the study. The selected pile length-diameter ratios are L/D = 5 to 60 and the ratio of eccentric length to pile diameter are e/D = 0 to 16. These dimensionless parameters can be simply used in practice as they were employed by Keawsawasvong [Sauer] to investigate the undrained capacity of loaded circular piles under combined horizontal load and moment. Furthermore, it may be used to validate the accuracy of the current solutions to those suggested by Keawsawasvong [63], Izadi and Chenari [25], and Brom [10]. For the range of e/D = 0-16, we follow the same length as proposed by Brom [10].
It is important to note that this study uses the rigorous LB and UB FELA to derive rigorous solutions of laterally loaded piles by also considering the soil weight, whereas the solutions in the original work by Broms [10] were obtained from the limit equilibrium method with the assumptions of weightless soils.

Comparisons with Published Results
To verify the numerical results presented in the paper, Figure 3 compares the freeheaded results of the present study with those from the displacement finite element method (FEM) by Keawsawasvong and Ukritchon [17] as well as the lower bound (LB) finite element limit analysis (FELA) by Izadi and Chenari [25]. The presented results are the average solutions of H/(s u LD) from the UB and LB solutions, and they are for cases of e/D = 0 and n = γL/s u = (0, 5, and 10). In general, the present average solutions are in good agreement with those published solutions. All solutions indicate that an increase in L/D results in a nonlinear rise of H/(s u LD). The larger the overburden stress factor n, the greater the H/(s u LD). For n = 0 (i.e., weightless soil), the LB solutions by Izadi and Chenari [25] are over-conservative as they are lower than the present solutions by 15-20%.
For the verification of fixed-headed piles, Figure 4 shows a comparison between the averaged LB and UB results with those from the displacement-based finite element method (FEM) by Keawsawasvong [63] and the limit equilibrium method (LEM) from Brom [10]. Numerical results have shown that the averaged LB and UB solutions agree very well with those solutions produced by Keawsawasvong [63]. Nevertheless, the LEM solutions by Brom [10] are much lower than the present solutions. This is due to the assumption of the failure mechanisms used in the LEM method by Brom [10] were not quite correct since the wedge-shaped failure was used in his analysis. It is necessary to re-examine the failure mechanisms assumed in their analytical UB and LEM studies. This has also highlighted the advantages of using the current FELA technique, in which a prior assumption of the failure mechanism is not needed [48].
A comprehensive comparison with the classic LEM solutions in Broms [10] is presented in Figure 5 for both fixed-headed piles and free-headed piles with e/D = 0-16. Since the work presented by the author did not consider the soil unit weight, therefore, the comparison is for the case of n = 0. Numerical comparisons have shown that Broms's LEM solutions are conservative, as they are consistently lower than those in the present study for all cases of both free-and fixed-headed piles.
The above three comparisons have improved the confidence in using the numerical results produced in the paper. Consequently, a series of parametric studies are presented in the next section. the average solutions of H/(suLD) from the UB and LB solutions, and they are for cases of e/D = 0 and n = γL/su = (0, 5, and 10). In general, the present average solutions are in good agreement with those published solutions. All solutions indicate that an increase in L/D results in a nonlinear rise of H/(suLD). The larger the overburden stress factor n, the greater the H/(suLD). For n = 0 (i.e., weightless soil), the LB solutions by Izadi and Chenari [25] are over-conservative as they are lower than the present solutions by 15-20%. For the verification of fixed-headed piles, Figure 4 shows a comparison between the averaged LB and UB results with those from the displacement-based finite element method (FEM) by Keawsawasvong [63] and the limit equilibrium method (LEM) from Brom [10]. Numerical results have shown that the averaged LB and UB solutions agree very well with those solutions produced by Keawsawasvong [63]. Nevertheless, the LEM solutions by Brom [10] are much lower than the present solutions. This is due to the assumption of the failure mechanisms used in the LEM method by Brom [10] were not quite correct since the wedge-shaped failure was used in his analysis. It is necessary to reexamine the failure mechanisms assumed in their analytical UB and LEM studies. This has also highlighted the advantages of using the current FELA technique, in which a prior assumption of the failure mechanism is not needed [48].   A comprehensive comparison with the classic LEM solutions in Broms [10] is presented in Figure 5 for both fixed-headed piles and free-headed piles with e/D = 0-16. Since the work presented by the author did not consider the soil unit weight, therefore, the comparison is for the case of n = 0. Numerical comparisons have shown that Broms's LEM solutions are conservative, as they are consistently lower than those in the present study for all cases of both free-and fixed-headed piles.  A comprehensive comparison with the classic LEM solutions in Broms [10] is presented in Figure 5 for both fixed-headed piles and free-headed piles with e/D = 0-16. Since the work presented by the author did not consider the soil unit weight, therefore, the comparison is for the case of n = 0. Numerical comparisons have shown that Broms's LEM solutions are conservative, as they are consistently lower than those in the present study for all cases of both free-and fixed-headed piles.

Results and Discussion
This section presents parametric studies for the following three dimensionless parameters, namely, the overburden stress factor (n), the pile length-diameter ratio (L/D), and the ratio of eccentric length to diameter (e/D). The effects of the individual parameter on the undrained lateral capacity of piles are examined using the normalized horizontal load factor H/(s u LD).
The influences of L/D on H/(s u LD) are presented in Figure 6 for the various values of (e/D = 0-16) and (n = 0-80). In Figure 6a, where n = 0, an increase in L/D results in a nonlinear increase in H/(s u LD). The larger the e/D, the smaller the H/(s u LD). Numerical results have shown that the fixed-headed piles yield greater values of H/(s u LD) than those of the free-headed piles. Interestingly, the difference between the two piles is not small which is about three times larger than the free-headed pile with e/D = 0. The same observation applies to other values of n, which are also presented in Figure 6b-f. The maximum and minimum values of each case (different n and e/D) can be found in Table 2. results have shown that the fixed-headed piles yield greater values of H/(suLD) than those of the free-headed piles. Interestingly, the difference between the two piles is not small which is about three times larger than the free-headed pile with e/D = 0. The same observation applies to other values of n, which are also presented in Figure 6b-f. The maximum and minimum values of each case (different n and e/D) can be found in Table  2.     shown in Figure 7b with L/D = 60, a reduction of H/(s u LD) as e/D increases is reported, where the function is almost linear when e/D is approximately larger than 4. One possible reason for this could be due to the better and more uniform stress distribution for piles with large L/D so that the impact of the eccentricity on the lateral force becomes less.   As for fixed-headed piles, Figure 9 presents the effects of L/D on H/(suLD) for the various values of n. Similar to the free-head piles, H/(suLD) increases nonlinearly with the increasing L/D. It is interesting to see that the gradient of the nonlinear curve decreases (the curve becomes flattered) as the value of n increases. Note that the smallest H/(suLD) corresponds to the cases of n = 0. Moreover, see the case of n = 80 with a near linear line. This is mostly due to the effect of the large overburden stress factor n = γL/su. It can therefore be concluded that the effect of n on H/(suLD) cannot be neglected in the design of piles subject to lateral loading. A similar observation can be made with respect to Figure  10     As for fixed-headed piles, Figure 9 presents the effects of L/D on H/(suLD) for the various values of n. Similar to the free-head piles, H/(suLD) increases nonlinearly with the increasing L/D. It is interesting to see that the gradient of the nonlinear curve decreases (the curve becomes flattered) as the value of n increases. Note that the smallest H/(suLD) corresponds to the cases of n = 0. Moreover, see the case of n = 80 with a near linear line. This is mostly due to the effect of the large overburden stress factor n = γL/su. It can therefore be concluded that the effect of n on H/(suLD) cannot be neglected in the design of piles subject to lateral loading. A similar observation can be made with respect to Figure  10, where the relationship between n and H/(suLD) is shown for the various values of L/D. A nonlinear increase is depicted as n increases. The larger the L/D, the greater the H/(suLD). The rate of increase (gradient) is the smallest when L/D is the largest. As for fixed-headed piles, Figure 9 presents the effects of L/D on H/(s u LD) for the various values of n. Similar to the free-head piles, H/(s u LD) increases nonlinearly with the increasing L/D. It is interesting to see that the gradient of the nonlinear curve decreases (the curve becomes flattered) as the value of n increases. Note that the smallest H/(s u LD) corresponds to the cases of n = 0. Moreover, see the case of n = 80 with a near linear line. This is mostly due to the effect of the large overburden stress factor n = γL/s u . It can therefore be concluded that the effect of n on H/(s u LD) cannot be neglected in the design of piles subject to lateral loading. A similar observation can be made with respect to Figure 10 It is common to use the contour plots of shear dissipation to depict possible failure mechanisms of a soil structure, as it provides a good indicator of the intensity of non-zero plastic strains. Technically speaking, the actual values of the contour are not important in such a perfectly plastic material model, and therefore the contour bars for these plots are not normally shown in technical documents. Figure 11 shows the failure mechanisms of a free-headed pile with six different overburden stress factors (n = γL/su). The chosen freeheaded pile is for (e/D = 0, L/D = 20). As n increases (i.e., the overburden stress increases), the size of the passive failure zone (see the area near the ground surface) decreases. Accompanying this is an increase in the size in the active failure zone (see the area near the pile tip, especially for n = 50 and 80). It is common to use the contour plots of shear dissipation to depict possible failure mechanisms of a soil structure, as it provides a good indicator of the intensity of non-zero plastic strains. Technically speaking, the actual values of the contour are not important in such a perfectly plastic material model, and therefore the contour bars for these plots are not normally shown in technical documents. Figure 11 shows the failure mechanisms of a free-headed pile with six different overburden stress factors (n = γL/su). The chosen freeheaded pile is for (e/D = 0, L/D = 20). As n increases (i.e., the overburden stress increases), the size of the passive failure zone (see the area near the ground surface) decreases. Accompanying this is an increase in the size in the active failure zone (see the area near the pile tip, especially for n = 50 and 80). It is common to use the contour plots of shear dissipation to depict possible failure mechanisms of a soil structure, as it provides a good indicator of the intensity of non-zero plastic strains. Technically speaking, the actual values of the contour are not important in such a perfectly plastic material model, and therefore the contour bars for these plots are not normally shown in technical documents. Figure 11 shows the failure mechanisms of a freeheaded pile with six different overburden stress factors (n = γL/s u ). The chosen free-headed pile is for (e/D = 0, L/D = 20). As n increases (i.e., the overburden stress increases), the size of the passive failure zone (see the area near the ground surface) decreases. Accompanying this is an increase in the size in the active failure zone (see the area near the pile tip, especially for n = 50 and 80).
In Figure 12, the effects of the ratio of eccentric length to pile diameter e/D on the failure mechanisms of a free-headed pile are presented. The chosen free-headed pile is for (n = 10, L/D = 30). As e/D increases (i.e., the pile length above the ground surface increases), the size of the passive failure zone (see the colored area near the ground surface) decreases. Nevertheless, unlike in Figure 11, the developed active zone near the pile tip is not pronounced. One of the possible reasons for this could be attributed to a location change of the inflection point of the pile as e/D increases. In Figure 12, the effects of the ratio of eccentric length to pile diameter e/D on the failure mechanisms of a free-headed pile are presented. The chosen free-headed pile is for (n = 10, L/D = 30). As e/D increases (i.e., the pile length above the ground surface increases), the size of the passive failure zone (see the colored area near the ground surface) decreases. Nevertheless, unlike in Figure 11, the developed active zone near the pile tip is  For the chosen fixed-headed pile study (e/D = 0, L/D = 30), Figure 13 shows the failure mechanisms of six different overburden stress factors (n = γL/s u ). By fixing the vertical movement of a fixed-headed pile, only horizontal translations are allowed. This is similar to the classical passive earth pressure problem. We can therefore expect that, for the weightless case with n = 0, a full passive failure zone can be developed. This is seen in Figure 13a. As n increases, the overburden stress increases, and the overall size of the passive zone reduces, resulting in partially mobilized passive zones. It can therefore be concluded that the overburden stress factor n plays an important role in the determination of piles subject to lateral loading. For the chosen fixed-headed pile study (e/D = 0, L/D = 30), Figure 13 shows the failure mechanisms of six different overburden stress factors (n = γL/su). By fixing the vertical movement of a fixed-headed pile, only horizontal translations are allowed. This is similar to the classical passive earth pressure problem. We can therefore expect that, for the weightless case with n = 0, a full passive failure zone can be developed. This is seen in Figure 13a. As n increases, the overburden stress increases, and the overall size of the passive zone reduces, resulting in partially mobilized passive zones. It can therefore be concluded that the overburden stress factor n plays an important role in the determination of piles subject to lateral loading.

Empirical Design Equations
For research to be practical and used by engineers, it is literally important to develop accurate design equations. By employing a nonlinear regression analysis to the computed solutions, a mathematical expression of H/(suLD) is proposed using the combined

Empirical Design Equations
For research to be practical and used by engineers, it is literally important to develop accurate design equations. By employing a nonlinear regression analysis to the computed solutions, a mathematical expression of H/(s u LD) is proposed using the combined relationship between power functions of L/D and n. This is shown in Equation (3) as follows: where (a 1 to a 3 ), (b 1 to b 3 ), and (c 1 to c 3 ) are constant coefficients and they are presented in Table 3 for practical uses. The least-square method proposed by Sauer [64] was employed to determine the optimal values of these coefficients. Even though both sides of Equation (3) have the parameters L and D, they are in fact dimensionless parameters. To use the design Equation (3), the values of e, L and D should be defined first. Then, the values of n, e/D, and L/D are computed based on the soil investigation results (e.g., γ and s u ). Substituting n and L/D into Equation (3), the value of H/(s u LD) is then obtained. By multiplying the value of H/(s u LD) with s u , L, and D, the result of H can be acquired which can be used as a maximum lateral force that a pile can resist. A comparison of the H/(s u LD) values between the FELA solutions (Avg) and the proposed design equation is shown in Figure 14 for the purpose of verifying the equation uses for the fixed-headed piles. For the free-headed piles, the comparisons are shown in Figure 15a-d, respectively, for the four different values of e/D. The comparisons have shown that the coefficients of determination (R 2 ) are all greater than 99.44%. The proposed Equation (3) is considered as highly accurate with the coefficients provided in Table 3, and it can be used with confidence in practice to estimate the lateral capacity of a pile in cohesive soil.
proposed design equation is shown in Figure 14 for the purpose of verifying the equation uses for the fixed-headed piles. For the free-headed piles, the comparisons are shown in Figure 15a-d, respectively, for the four different values of e/D. The comparisons have shown that the coefficients of determination (R 2 ) are all greater than 99.44%. The proposed Equation (3) is considered as highly accurate with the coefficients provided in Table 3, and it can be used with confidence in practice to estimate the lateral capacity of a pile in cohesive soil.

Design Charts Limitations
The undrained lateral resistance of circular piles is investigated using the UB and LB limit analysis methodologies. Note that the following are some limitations that need to be investigated further: 1. All computed design charts are based on the assumption of undrained clay and Data Availability Statement: All data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.