Development of lateral capacity-based envelopes of twin-pile group under combined V–H–M loading

This study explores the bearing capacity of a three-dimensional (3D) twin-pile group subjected to combined vertical (V), horizontal (H), and bending moment (M) load using finite element limit analysis coupled with lower and upper bounds and mixed Gauss element formulations. The numerical results based on different element formulations are compared with those of the previous study under single-component loading. The failure envelope in the 3D V– H–M loading space is constructed through a simplified loading strategy. Several typical collapse mechanisms are identified to be associated with different load component values. Moreover, the influence of the pile spacing-to-diameter ratio, pile length-to-diameter ratio, and loading direction angle on the combined bearing capacity is discussed herein. A set of charts is developed as well. Finally, a closed-form design expression of failure envelopes considering the abovementioned influential parameters and the combined V–H–M load are derived based on the obtained numerical results.


Introduction
The evaluation of the lateral bearing capacity of the pile group in clay subjected to external loads is a critical issue for pile foundation designs. Considering the group effect, a number of investigations on pile groups have been performed experimentally (e.g., Matlock et al. [1] and Rollins et al. [2]), numerically (e.g., Zhao et al. [3], Pham et al. [4], and Sheil [5]), and analytically (e.g., Georgiadis et al. [6][7] and Zhao et al. [8][9]). However, most of these works analyzed the pile group for vertical and lateral loads separately, ignoring the external load uncertainty. Existing studies that performed investigations under arbitrary loading (i.e., a combination of vertical forces V, horizontal forces H, and bending moments M) mostly focused on footings, buckets, or individual piles [10][11][12], and a rather limited number of investigations have considered a pile group. This is due to the undrained stability of pile groups under a combined V-H-M load being a complex problem because of the influence of many factors, including the pile diameter, length, spacing, and configuration, loading direction, and coupling effect between the load components. Intense analysis work is required to consider these factors. In view of such a context, the discovery of a robust mathematical tool for efficient computation and the provision of valuable insights on the governing failure mechanism are urged. The finite element limit analysis (FELA) method offers a suitable alternative. Extensive commendable studies on the complex problems of the bearing capacity of foundations and the stability of slopes or trenches were recently conducted using the FELA method [13][14][15][16]. Exact limit analysis solutions are not available for pile groups; hence, this study aims to evaluate the influence of combined V-H-M load interaction on the bearing capacity of a twin-pile foundation, which is an essential group configuration that considers pile-to-pile interaction [6][7], embedded in clay. A three-dimensional (3D) failure envelope in the V-H-M loading space is constructed following a simplified loading strategy. The representative failure mechanisms observed in the soil-pile system are also discussed. Different values of the pile length-to-diameter and pile spacing-to-diameter ratios and different loading plane directions are considered and thoroughly analyzed. Finally, a closed-form expression is presented to fit the failure envelopes obtained herein.

Problem statement
The present study deals with rigid twin-pile groups of length L, diameter D, and center-to-center spacing s subjected to a combined load (Fig. 1). A rectangular pile gap connects the two piles on their head. The pile shafts are embedded into a deep homogeneous undrained clay, where su = undrained shear strength, and = unit weight of clay. In this system, the soil-pile interface is perfectly rough, indicating that the interface can maintain the full strength of the neighboring soil. The external force is conceptualized as a combination of vertical, horizontal, and bending moment loading in the same plane applied on the pile cap tip. The sign for H and M complies with the right-handed axes and clockwise positive convention, while V is highlighted with a superscript '+' for compression (Fig. 1). The bearing capacities for the single-component loading, V, H, or M, are denoted as V0, H0, and M0, respectively. Note that the present study focuses on the undrained behavior of piles embedded in clayey soil. The pile behavior in other soil types with different properties is out of the scope of this study.

Analysis methodology
The combined bearing capacity was computed herein via the well-known FELA method. The FELA method is a powerful computational tool combining the capabilities of the finite element discretization for dealing with complex geotechnical problems and the limit analysis principles for bracketing the exact limit load by upper-and lower-bound solutions. One of the impressive merits of the FELA method is that there is no need to make assumptions about the failure modes in advance, which provides abundant flexibility in the practical application of FELA. The detailed description of the theorems of the limit analysis and the FELA is beyond the scope of this paper. The reader is referred to [17] for information.

Finite element limit analysis model
The 3D finite element models of the soil-pile system were built in OptumG3 [18] (Fig. 2), a novel, but robust commercially available FELA software with advanced optimization algorithms. The soil domain was considered as a 3D rectangular box. The two piles were regarded as multifaceted prisms with 36 facets. The pile cap on the top was modeled by a rectangular plate element (without thickness). The soil material behavior was described using an elastoplastic Tresca constitutive model with Young's modulus Eu = 30 Mpa, undrained shear strength su = 50 kPa, and unit weight = 20 kN/m 3 . The piles and the cap were assumed weightless and rigid to omit the structure's weight contribution to the overall capacity. The contact surfaces between the piles and the soil mass were modeled using shear joint elements with strength similar to that of soil. The model size was determined via a sensitivity analysis. The model length varied from 20D to 40D, whereas the distance from the bottom boundary to the pile tip was from 2D to 8D. As shown in Table 1, the results became insensitive to the model size variation when the dimension expanded beyond 30D (length) × (L + 5D) (depth), which was the dimension of the models used in all the analyses (Fig. 2). The standard boundary conditions of the roller in all sides and the fixed bottom were applied to the model. Due to the soil-pile system symmetry, only half of the model was considered in all the analyses. Table 1. Influence of the model size on pure ultimate loads (L/D = 9 and s/D = 6). Three-dimensional tetrahedron elements with upper bound (UB), lower bound (LB), and mixed Gauss (MG) element formulations were used to discretize both the soil mass and the piles. The UB and LB elements were constructed directly based on the upper-and lower-bound principles and expectedly led to rigorous upper and lower bounds on the exact solution. Meanwhile, the MG elements, in which the elastoplastic Hellinger-Reissner principle is considered, are often more accurate, but do not lead to rigorously bounded solutions.
Exploiting the unique technique of mesh adaptivity embedded in Optum G3, acceptable results can be obtained by applying an appropriate number of elements and adaptive steps. In the simulations, three rounds of mesh adaptivity were implemented for the refinement procedure with the number of elements increasing from 1000 to 15,000.

Figure 2.
Dimensions of the symmetric finite element model.

Load paths and parameters
A specific loading strategy in the V-H-M space must be determined to construct the failure surface of a twin-pile group under the combined V-H-M loading. Among the various loading strategies available, the most common approach is to consider a single-load multiplier that affects all the loading components in spherical coordinates. The strategy adopted herein takes more account of the load feature in reality. In view of the fact that the vertical load imposed on the pile cap is generally from the self-weight of the superstructure and does not vary much, the soil-pile system collapse is mostly caused by the increase of H or M. Therefore, a more convenient method of constructing the capacity envelope is to keep the V value constant and apply H and M under a fixed ratio (M/H = tan) [10][11][12] (Fig. 3). Each radial loading path in the H-M plane is identified by the angle with regard to the horizontal axis. Due to the symmetry in the horizontal load application on the H-M plane, only the right half of the coordinate system was investigated (= −90°to 90°). The eight V loading levels considered were V = 0V0, 0.25V0, 0.5V0, 0.75V0, 0.8V0, 0.85V0, 0.9V0, and 0.95V0. In the 3D V-H-M space, at least 120 simulations were required to complete a failure envelope for a certain twin-pile model. The effects of the varying pile length-to-diameter ratio (L/D = 3, 6, and 9), pile spacing-to-diameter ratio (s/D = 2, 3, 4, 5, and 6), and loading plane direction (along the x-axis (Loading Condition A) or the y-axis (Loading Condition B)) on the failure envelope were further investigated. More than 2100 FELA models were established to perform such a comprehensive parametric analysis.

Development of the failure envelope
The ultimate bearing capacity of the twin-pile groups under pure loading (V0, H0, or M0) is investigated in this section. The failure envelope projection on the H-M plane is presented before the 3D envelope construction in the V-H-M space.

Estimation of the ultimate loads for pure loading
The dimensionless bearing capacities V0/suD 2 , H0/suD 2 , and M0/suD 3 computed using the UB, LB, and MG element formulations, respectively, are shown and compared with those in an available study (Leng et al. [12]) ( Table 2). The results based on the MG element formulation fell within the bracket of the numerical LB and UB solutions for all pile lengths and spacings. The difference among the three types of results was acceptably below 5%. Moreover, the MG element formulation predicted the dimensionless capacities, especially (M0/suD 3 )x, to be closer to the results obtained by Leng et al. [12], although such comparisons were only within limited cases. In view of the foregoing, the FELA results based on the MG element formulation offered the most reliable estimates and were used for further analyses in this study. Table 2 also lists some conventional rules. Increases in the pile length (L/D) and the pile spacing (s/D) generally led to expected increases in the dimensionless capacities, although the (H0/suD 2 )y and (M0/suD 3 )y increments were rather limited. The effect of the loading direction on the lateral bearing capacity was also predominant, which is exactly as Georgiadis et al. [7] emphasized. For the worst case of L/D = 3 and s/D = 6, a loading angle rotation of 90°nearly halved the magnitude of the dimensionless horizontal bearing capacity. The investigations on the ultimate loads for pure loading, which will be elaborated in the subsequent subsections, aim to determine the edges of the failure envelope in the V-H-M space and better assess the coupling effects of the multiple loads.   Fig. 4, the increase of the bending moment linearly reduced the value of the ultimate horizontal load. This was caused by H and M acting in the same direction and the presence of M that impelled the soil-pile system vulnerable to failure. In the fourth quadrant (Fig. 4) corresponding to the case of H and M acting in the opposite direction, the ultimate horizontal load first increased, and then decreased as the M/H load path varied. Note that the maximum H the piles can mobilize does not occur under pure loading, but rather for a certain combination of positive H and negative M. This turning point corresponding to P1 for Loading Condition A and P2 for Loading Condition B roughly divided the curve into two parts in light of the motion of the twin-pile group. Before this point, the piles moved in the horizontal load direction, and more shear dissipations were produced in front of the moving piles (Figs. 5a and c); thus, an increasing H was required to overcome the opposite M. After the turning point, the limit state of the twin-pile group was achieved for the radial path with a high negative M/H ratio. The moment load dictated the pile motion direction (Figs. 5b and d), and energy dissipations occurred around the twinpile group. Comparing the failure surfaces for different loading conditions, the loading plane rotation of 90°did not change the curve shape, but shrank the failure surface in the H-M space. This result indicates the reduction in the overall bearing capacity due to the variation of the loading direction. More discussions on this point can be found in the "Parametric analyses" section.

Failure envelope in the V-H-M loading space
The effect of the vertical load was then considered, complying with the loading strategy used in this study. Eight levels of constant V were involved in the construction of the failure surfaces on the normalized H/H0-M/M0 plane (Fig. 6). The coupling capacity of the horizontal load and the bending moment generally reduced with the increasing vertical load. The outermost failure surface in Fig. 6 corresponded to the case without V. The vertical load variations within the range of 0V0 to 0.5V0 produced a limited impact on the size of the failure surfaces. The failure surface contraction became increasingly obvious only when the vertical load increased beyond 0.75V0. Hence, the number of samples with a vertical load between 0.75V0 and 0.95V0 was increased to improve the accuracy of the failure surface prediction (Fig. 6). Although such an evolution law of the failure surface with vertical load has been observed for the twin-pile group with a certain geometry (i.e., L/D = 6 and s/D = 4) and under a certain planar load direction (i.e., along the y-axis), it has been discovered suitable for all cases investigated with various combinations of the influential parameters. In addition, the presence of V may change the shear dissipation distribution in the soil mass. Fig. 7 presents the collapse modes of the pile-soil system under different vertical component loads. For both cases, the failure mechanism is a mixture of the pure vertical-loaded mode and the pure lateral-loaded mode. The vertical-loaded feature of the mixed mechanism increasingly dominated with the increase of the V value following the arrow direction from P3 to P4 in Fig. 6. Similar mechanisms can be broadly observed for other influential parameters. All data points for the loading paths concerning various V levels were subsequently integrated into the 3D loading space of V/V0-H/H0-M/M0. Fig. 8 depicts the overall failure envelope.

Parametric analyses
The effects of the vertical load on the combined capacity were extremely similar for various influential parameters (e.g., pile length-to-diameter ratio L/D, pile length-to-diameter ratio s/D, and loading plane direction); hence, the parametric analysis in this section can be simplified and considered in the H-M loading space. Finally, a closed-form expression is proposed to predict the combined bearing capacity.

Effect of L/D
Taking s/D = 4 and Loading Condition A as an example, Fig. 9 depicts the influence of the pile length on the failure envelope. As expected, the overall bearing capacity of the twin-pile group sharply increased when the pile length increased. Fig. 9a also shows that for the larger L/D, the maximum horizontal load the piles can sustain occurred under the radial path with a higher negative value of M/H. The L/D increase improved the coupling effect of the combined H-M load. A larger value of the maximum normalized load (H/H0 or M/M0) can be obtained at a certain radial path when L/D increases (Fig. 9b). A similar trend was observed for the other cases investigated herein.

Effect of s/D
Considering L/D = 6 and Loading Condition A first, Fig. 10a illustrates the influence of the pile spacing on the failure envelope. The combined bearing capacity of the twin-pile group gradually increased when the pile spacing increased. Fig. 10b also shows that the maximum value of H/H0 that can be obtained at a certain path was the lowest for the medium values of s/D (e.g., s/D = 4 and 5) compared with the cases of large or small s/D (e.g., s/D = 2 or 6). The evolution trend of the failure surfaces with s/D was totally changed when the planar load direction varied. For the case of L/D = 6 and Loading Condition B (Fig. 11), the effect of s/D on the bearing capacity of the twin-pile group was limited, especially when the curves were plotted on the normalized H/H0-M/M0 plane. In addition, the L/D variations will not change the abovementioned trends.

Effect of the loading plane direction
The bearing capacity of the twin-pile group not only depends on the pile geometry and configuration, but also on the loading plane direction. Such is depicted in Fig. 12 for the case of L/D = 3 and s/D = 6. When the loading plane rotated from Loading Condition A to Loading Condition B, the failure surface on the dimensionless H/suD 2 -M/suD 3 plane sharply contracted, and the ultimate bearing capacity under pure loading (for H = 0 or M = 0) was reduced by nearly 50% compared with its original value (Fig.  12a). The maximum value of H/H0 investigated for Loading Condition B was greatly larger than that for Loading Condition A (Fig. 12b) when the curves were plotted on the normalized H/H0-M/M0 plane.
where A = A1 and q = q1 when the loading plane direction is along the x-axis and A = A2 and q = q2 when it is along the y-axis. Fig. 13 presents comparisons of the proposed expressions and the numerical results obtained in the present study for the case of L/D = 6, s/D = 4, and Loading Condition A. For the special case of V = 0 (Fig. 13a), the results from Leng et al. [12] were overlapped to further validate the proposed formulas. A reasonable agreement can be observed generally.

Conclusions
In this study, a series of FELAs using the LB, UB, and MG element formulations were performed to explore the bearing capacity of a 3D twin-pile group under the combined vertical, horizontal, and bending moment load. First, under single-component loading, the numerical results based on different element types were presented and compared with data from other published papers. Subsequently, the failure envelopes in the 3D V-H-M loading space were constructed by varying the ratio of H/M on several H-M planes for different V levels. The effects of the L/D, s/D, and the loading plane direction on the failure surface of the twin-pile group were discussed. The empirical expression was additionally given to predict the combined bearing capacity. Some detailed conclusions have been summarized as follows: (1) The failure surface on the H-M plane presents an inclined elliptical shape, and the pile motion direction is judged through the turning point on the curve. Lastly, note that the conclusions drawn in this study are only applicable for clayey soil foundations. The behaviors of twin-pile groups in soil regions with other property types will be explored in future works.