Three-dimensional Numerical Simulation and Validation of Load-settlement Behaviour of a Pile Group under Compressive Loading

Settlement of pile foundation is one of the controlling pile design parameters and its numerical simulation is one of the techniques widely used to predict the settlement behaviour of piles. This study was on the settlement behaviour of a pile group located in silty-sand deposits using the finite element (FE) approach which is based on the static pile load test. Three different types of analyses were investigated: (1) a linear elastic (LE) analysis in which the soil was assumed to be linear elastic, (2) a complete nonlinear analysis in which the soil adjacent to the pile shaft as well as the soil between the piles were modelled using the Mohr Coulomb (MC) or hardening soil (HS) model, and (3) a combined analysis in which the soil close to the pile shaft was modelled using the HS model while the soil in the remaining area was modelled using the LE or MC model. Numerical results obtained for the load-settlement behaviour of the pile group were validated using the popular RATZ analytical approach. The results of the FE analysis suggest that incorporating a nonlinear zone of soil close to the pile shaft as an interface and leaving the soil beyond this zone as linear elastic give a more reasonable estimation and a much better prediction of the pile group settlement. It is also suggested that for a typical pile group, a nonlinear interface of thickness equal to the pile diameter and extending from the pile shaft to the edge of the zone would be sufficient to capture the load transfer mechanism. The group settlement ratio predicted in this study is in good agreement with the findings made in the previous studies.


Introduction
Over the years, construction of many tall buildings and heavy structures has involved pile group foundations because of the high bearing capacity of these foundations [1,2,3]. Generally, pile-group foundations involve three-dimensional interaction among the pilecap, soil and piles. In the past, problems that arose because of these interactions had been often solved by making simplifying assumptions regarding the geometry and the properties of the materials used [4]. However, because of the realistic nature of the problem, it has become essential to allow for threedimensional geometry, interface effects and non-linear soil properties. Therefore, the purpose of this paper is to present a threedimensional finite element study that allowed for various interface zones and a wide range of material models such as linear elastic, nonlinear Mohr-Coulomb and nonlinear hardening soil models.
In the recent years, a number of research studies [5 -19] have been carried out using analytical and numerical procedures to study the load-settlement behaviour of piles while making simplifying assumptions during the modelling of the behaviour of the pile-soil system.
Poulos [5] has pioneered an analytical method on pile group displacements and introduced the concept of 'interaction effect' of pile groups. Thereafter, many researches [6][7][8] have been undertaken on various analytical methods to study the vertical deformations of pile groups. Butterfield and Banerjeet [6], Chow [7] and Shen et al. [8] have presented analytical solutions for the vertical deformation of pile groups using theoretical load-transfer curves. Randolph and Wroth [9] have proposed a closed-form solution for the settlement of single piles. They have also used the approach to develop an approximate analytical method for pile groups as well using theoretical loadtransfer curves. Lee [10] has presented a simple discrete layer approach using linear elastic solutions for the settlement analysis of axially loaded pile groups. Randolph [11] has developed a simplified RATZ analytical method that uses parabolic or hyperbolic load transfer curves to describe individual pile settlement behaviours, and to derive the group pile interaction effect using the elastic solutions developed by Randolph and Wroth [9]. Seo and Prezzi [12] and Fan et al. [13] have proposed an analytical solution for vertically loaded piles in multi-layered soil and expansive soil respectively based on the theory of pile-soil interaction.
On the other hand, in the past decade, there had been a rapid increase in the development of finite element methods for use in different types of geotechnical engineering applications, since in principle, these methods can deal with soil inhomogeneity and non-linearity in a consistent manner [14,15]. As such, Zakia et al. [16] have carried out a finite element analysis on the reliable selection the modelling parameters for the settlement predictions of pile foundations. Jun Ju [17] and Fuchun et al. [18] have carried out settlement analysis using PLAXIS 3D finite element package for pilegroups located in sleech strata and soft clayey soil respectively. Alnuiam et al. [19] and Muqtadir [4] have conducted a non-linear three-dimensional finite element analysis of sand to capture the deformation behaviour of pile-groups. Jian-lin et al. [1] have conducted a settlement analysis of pile-groups located in deep clayey soil deposits and have proposed a valuable equation to correct the compression modulus of soil during the simulation of deep soil conditions. Obviously, the 3D finite element method is one of the most precise and time-effective approaches for analysing the behaviour of a pile group [2,17]. Such an approach is applicable to pile groups in layered soils with nonlinear soil properties, and helps to understand better the behaviour of pile groups. Therefore, the aim of this research was to numerically simulate the load-settlement behaviour of a vertically loaded group-pile foundation using the PLAXIS 3D numerical package and based on the case history of a pile foundation located in deep silty-sand deposits.

Case Study used in the Research
The case study adopted for this research is similar to that used by Gowthaman et al., [2], in which they have used the load-settlement behaviour of a pile load test conducted on a single pile in the north western part of Singapore (Woodland). The same stratigraphy was used in this study to simulate the loadsettlement behaviour of a square pile group comprising four piles. Gowthaman et al. [2] have conducted a finite element study of an axially loaded single pile located in a deep siltysand deposit called "old alluvium". In their analysis, settlement behaviour of a single pile located in sandy-silt has been simulated and validated using field test (static pile-load test) results. Silty sand is one of the soils that has an economic importance and has often been encountered in recent major civil engineering projects, although it is still not well understood [20,21]. Also, because of the depositional process, old alluvium deposits consist of many layers and it has been found that even within a single layer, its properties are highly variable [20,22]. Elly [22] has concluded that silty-sand deposit, old alluvium, is an excellent material for construction. During the past decade, many major engineering projects such as those involving the construction of high-rise buildings, expressways, Deep Tunnel Sewage Systems (DTSSs) and Mass Rapid Transit (MRT) lines have been undertaken on these silty-sand deposits [22] and all of these massive and deep constructions have been mostly on pile foundations. Therefore, this research mainly focused on the settlement behaviour of pile group foundations located in deep siltysand deposits. In this study, a hypothetical 2 × 2 square pile group with a symmetric arrangement ( Figure  1a) was considered for the simulation of the settlement behaviour of axially loaded group piles located in deep silty-sand deposits. Typically, the spacing among the piles in a pile group will be in the range between 2D and 4D where D is the width or diameter of a single pile [23]. However, the critical design spacing of the pile group to be located in loose or medium-dense sandy soils is found to be 4D [23], and thus it is this critical spacing that was considered in the numerical analysis of the group pile. Each pile in the pile group was 48 m in length and had a diameter of 1500 mm. They were all driven to a depth of 47.5 m in the ground with a free-standing length of 0.5 m above the ground surface ( Figure 1b). The single pile geometry, soil stratigraphy and ground water table location were similar to those used by Gowthaman et al. [2].

Finite Element Model Parameters
The finite element soil models used in this study were the linear elastic (LE), Mohr Coulomb (MC) and hardening soil (HS) models. The linear elastic model is based on Hookes's Law of isotropic elasticity. It involves two basic elastic parameters i.e. Young's modulus (E) and Poisson's ratio (v). Although the linear elastic model is not suitable to model soil, it may be used to model stiff volumes of soil or stiff formulations in soil [2,17,23].
The Mohr Coulomb model is one of the simple nonlinear models used in this study. It is based on soil parameters that are found in most of the practical situations. It has five input parameters, i.e., elastic modulus (E) and Poisson's ratio (vs) for soil elasticity, friction angle (φ) and cohesion (c) for soil plasticity and dilatancy angle (ψ). However, this model does not include all of the non-linear features of soil behaviour.
The hardening soil model is an advanced nonlinear model used in the simulation of soil behaviour. It is based on the framework of classical theory of plasticity. In this model, total strains are calculated using a stress level dependent stiffness with a hyperbolic stressstrain relationship that defers between virgin loading and unloading/reloading ( Figure 2).

Figure 2 -Hyperbolic stress-strain relationship used in the hardening soil model
In the hardening soil model, stiffness is described much more accurately than in other models through the use of three additional stiffness input parameters: triaxial loading stiffness (E50), triaxial unloading stiffness (Eur) and oedometer loading stiffness [17,23]. All these stiffnesses are specified for a given reference stress (p ref ) and are assumed to vary directly with the horizontal effective stress raised to the power of m. The reference stress (p ref ) used in the stiffness measurement (laboratory level) generally falls in the range of 100-200 kPa. Stiffness power in stiffness laws (m) lies in the range from 0.5 to 1. As input parameters, the hardening soil model also requires normally consolidated coefficient of earth pressure at rest (K0), Poisson's ratio for unloading and reloading (vur) and the failure ratio (Rf). A vur of 0.2 and a Rf of 0.9 are appropriate for the HS model under drained conditions [12].
Model parameters used for different soil models (LE, MC and HS models) are summarized in Table 1.

Finite Element Analysis of the Group Pile
The PLAXIS 3D numerical package was used in this finite element study of settlement behaviour of an axially loaded group pile located in deep silty-sand deposits. The soil was modelled using 15-node wedge elements. These elements composed of 6-node triangular faces in the work planes produced by 2D mesh generation, and 8-node quadrilateral faces in the y-direction. Lateral sides of the computational domain were taken sufficiently away from the pile to avoid the boundary effect. Borehole option in PLAXIS 3D was used to define soil stratigraphy, ground surface level and pore pressure distribution in the sub-soil.

Figure 3 -Meshed pile foundation
Geometry lines were used as and when necessary to separate or isolate soil types from one another. Separate boreholes were used to differentiate characteristic properties and soil types. Fine mesh analysis was used throughout the numerical study of the settlement behaviour, to improve the convergent result. The meshed geometry and the meshed pile foundation are shown in Figure 3. Furthermore, all group pile models used in the different analyses were developed using a working area of 70 m × 25 m × 25 m. All movements at the bottom of the models were restrained along with all lateral movements perpendicular to the boundary and at the lateral sides.
A linear elastic non-porous and isotropic material model was used to represent the piles and a linear floor element was used to model the pile cap. The material properties of the piles and the pile cap used in the model are given in Table 2. During the 3D analysis, the PLAXIS 3D programme automatically generated interface elements around the piles. The following three types of FE analyses were performed during the numerical study:  Figure 4) was modelled using the HS model, while the soil in the remaining area (Zone B in Figure 4) was modelled as LE or MC material.
In the combined analysis, three different sizes for the Zone A interface were selected: (i) Zone extending to a distance (d) of D from the pile shaft, (ii) Zone extending to a distance (d) of D/2 from the shaft and (iii) Zone extending to a distance (d) of D/4 from the shaft where D is the pile diameter.

Figure 4 -Pile group layout of the combined analysis (NL-LE/ NL-NL)
A study on the effect of pile spacing on the group load-settlement behaviour was also undertaken in this research. After comparing the different types of analyses chosen with the RATZ approach, the finite element model that was found to be the best among the three models was chosen to predict the actual pile group settlement of the group pile located in sandy soil. The model that agreed most with the RATZ approach was used to analyse the group settlement behaviour with piles in the group placed at different spacings. In this study, the settlement behaviour analysis was conducted for the following three spacings (centre to centre): (i) S = 2D, (ii) S = 4D and (iii) S = 8D where D is the pile diameter.
Finally, a comparison of the load-settlement behaviours of the group pile for the three spacings was done, and the results obtained were used to determine the stiffness of the pile group for each of the spacings. The study was extended to determine the pile group settlement ratio as well. The group settlement ratios obtained from the numerical study for different pile spacings were compared with the findings of previous researches. Moreover, the pile group efficiency of the pile group modelled for each pile spacing was determined using the analytical methods available and the results obtained were compared with the findings of previous researches.

Pile Group Settlement Prediction Using RATZ Analytical Approach
The RATZ approach is one of the best methods proposed by Randolph [11] to predict the loadsettlement curve of a group pile by using a single pile load-settlement curve. It has been said that the results obtained from the RATZ analytical method could provide very good agreement with both numerical and experimental analyses [17]. This simplified settlement prediction method can be applied only when both the single pile and the group pile have been located in the same stratigraphy and with the same average pile head load. In this approach, by specifying a group settlement ratio (Rs), as illustrated in Figure 5, the elastic part of the load transfer curve was first factored to allow for interaction among the piles in the group and the nonlinear component of the load transfer curve of the single pile was added thereafter to the factored elastic component to obtain the group load transfer curve.

Figure 5 -Modeling group effects by factoring the load transfer curve (after Randolph [11])
The results of the static load test performed on the reference single pile were used to obtain the settlement behaviour of the group pile using the RATZ analytical approach. The group settlement ratio for the pile group was obtained from a previous study [24] based on the geometry of the pile group. The group settlements predicted were compared with the results obtained from the finite element analyses.

Simulation of the Load-settlement
Behaviour using Different Material Models and Interface Thicknesses The suitability of the different material models was analysed using a simple 2 × 2 symmetric square pile group with spacing of 4D (critical design spacing of a pile group located in loose or medium-dense sandy soil). Load-settlement behaviours predicted from the LE and NL analyses and the group pile load-settlement curve obtained from the RATZ analytical approach are shown in Figure 06. Based on the results, the settlement of the group pile measured using the RATZ approach is 8 mm for an average pile head load of 7,000 kN. For the same working load, the settlement predicted from the HS model is 16 mm which is about 2 times the RATZ prediction. Settlement values predicted from the linear elastic (LE) model and the Mohr Coulomb (MC) model for the same working load are 12 mm and 12.5 mm respectively. Thus, it can be seen that the settlements derived from both LE and MC analyses are higher than the field measurement (using the RATZ approach) obtained for a typical working load, by about 1.5 times.

Figure 7 -Comparison between RATZ prediction and finite element analysis results with corrected Moduli
It can be clearly observed that there is a significant difference between the values calculated for the settlement of the group pile and the corresponding values obtained from the RATZ prediction curve which indicates that the finite element analysis either significantly over predict the pile head settlement or underestimate the pile head stiffness. The differences between the RATZ predictions and the results of the finite element analysis are due to the ignorance of the deep in-situ effect in the modulus of soils. In a real field situation, the modulus that is obtained from laboratory tests can significantly differ from the in-situ modulus of deep soil and at times the difference can be very high [1,2]. Therefore, the equation proposed by Jian-lin [1] was used to increase the accuracy of the settlement prediction. The elastic modulus obtained from laboratory tests was modified as follows using Equation 1: where z is the depth of the soil layer (m), h0 is the reference depth (generally 1 m), Es,0.1−0.2 is the compression modulus obtained in the laboratory under a pressure in the range of 100-200 kPa, and β is the plasticity of the soil. Using the specifications give in the BS code, the value for β can be obtained based on the liquid limit and the plasticity index (IP) values obtained from the data. For silty sandy soils, β would be between 3.5 and 5 [1].
Finally, all three types of FE analyses were repeated using the corrected modulus (Es,z) obtained using the equation proposed by Jianlin [1]. Figure 7 shows the settlement results obtained from the three types of models and their comparisons with the RATZ predictions. From the results of the comparison between RATZ predictions and each of the finite element analysis results, it can be seen that the prediction of the pile head settlement has improved. At an average pile head load of 7,000 kN, the settlement calculated from either LE or MC analysis is 8 mm, whereas the settlement obtained from the RATZ prediction is also exactly 8 mm. At the same time, the settlement predicted from the HS model at the same working load is 11.2 mm which is 40 % higher than the RATZ analytical prediction. Although LE and MC models showed better agreement with the RATZ prediction at lower working loads (<12,000 kN), it was not so at higher working loads (>12,000 kN). The Mohr Coulomb (MC) model is a simple nonlinear model which does not accommodate all of the nonlinear parameters of soil and therefore it can be used to accurately predict only a limited level of settlement behaviour (up to 12,000 kN). The MC model is the most suitable model for settlement prediction of group piles at small working loads. Even though the HS model overpredicts the settlement or underestimates the soil stiffness, it is the only model that will capture the actual nonlinear behaviour of soils. The settlement contours obtained from the finite element analysis using LE, MC and HS models are shown in Figures 8 (a), (b) and (c) respectively. As the nonlinearity affects the interactions within the pile group, LE and MC models will fail to predict the settlement behaviour of a group pile. The soil nonlinearity and the interaction effects are clearly captured by the HS model as shown in Figure 08c and the incorporation of the two nonlinear effects provide for the accurate prediction of the settlement behaviour of group piles. However, the finite element results obtained from the calculations provide only a limited level of agreement with the analytical RATZ results, and it can thus be concluded that combined models will be able to predict better the settlement of group piles. It is well-established that the modulus of a soil mass will decrease with increasing strain levels. In a group of piles, the strain level can increase as the pile shaft is approached and hence the stiffness of the soil in this narrow zone close to the pile shaft will be smaller than that in the space between the piles at some distance away from the pile shaft. Therefore, to account for this stiffness variation, two types of combined analyses were considered: (i) HS-LE and (ii) HS-MC. Results obtained for the settlement of the pile group using HS-LE and HS-MC analyses for different interface thicknesses are given in Figures 9 and 10 respectively. From these findings, it can be seen that the pile group load-settlement behaviour can be better predicted by having a nonlinear interface of thickness equal to pile diameter, while keeping the remaining zone as either LE or MC. However, the HS-MC (d = D) model is unable to predict the actual behaviour between 5,000 kN and 17,000 kN, as it over predicts the settlement in that range. When the completely nonlinear Mohr Coulomb model was used for the remaining zone, the actual stiffness of the strata got reduced. Therefore, by considering the interface with a thickness equal to the pile diameter and extending from the pile shaft an advanced nonlinear HS model while considering the remaining zone as a LE model, it will be possible to predict the actual settlement behaviour of the pile group. Jian-lin [1] and Lee and Poulos [25] also have concluded that as compared to a single model, the combined model which combines a nonlinear model with a linear model (NL-LE) would provide more realistic predictions of the settlement behaviour of pile groups. However, they have suggested that a nonlinear interface of thickness equal to half the pile width and extending for a distance from the pile shaft (d = D/2) would be sufficient to capture the load transfer mechanism, whereas the findings of this study reveal that the interface thickness should be equal to the pile diameter.

5.2
Effect of Spacing A study of the effect of pile spacing on the group load-settlement behaviour was done by using the model identified as best and based on the outcome provided in Section 5.1. Hence, the combined model that uses a thickness equal to the pile diameter for the interface that extends from the pile shaft and models that interface using the HS model with the remaining zone modelled as MC was used for the analysis. The results which indicate the effect of the spacing are given in Figure 11. The settlements predicted for an average working load of 7,000 KN, at S = 2D, S = 4D and S = 8D are 13 mm, 8 mm and 4 mm respectively. From Figure 11, it can be seen that the settlement of the pile group decreases as the spacing is increased. When the spacing is increased by four times, the settlement decreases to approximately one third of its original value.
From Figure 11, it can also be seen that the pile group stiffness decreases as the pile spacing decreases. The stiffness of the pile group at a 4D spacing is higher than its stiffness at a spacing of 2D by 1.5 times. Also, the stiffness of the pile group at a spacing of 8D is more than 2 times its value at a spacing of 2D. This result is also supported by Pressley and Poulos [26] who have shown that the block failure mechanism occurs at closer spacings (Figure 12a) with significant plastic zones developed below the group and only the full pile-slip developed along the outer piles. From Figures 11 and 12, it can be seen that when the spacing is eight times the diameter (8D), the group settlement behaviour gets closer to the settlement behaviour of a single pile. As the pile spacing increases, the failure mechanism gradually changes to the "single pile" mode ( Figure 12b), in which full pile-soil slip occurs along the full length of each individual pile.

5.3
Group Settlement Ratio The group settlement ratios obtained from the finite element analysis for three different spacings and under average ultimate load conditions were compared with the results obtained in previous researches. Pressley and Poulos [26], using the finite element analysis and the elastic method have conducted a study of the group settlement behaviour of various sizes of pile groups. They have also calculated the group settlement ratios of pile groups of various geometries. The group settlement ratio (Sg/S) can be defined as the ratio of the settlement of the pile group (Sg) to the settlement of the single pile (S). The comparisons of the results are shown in Figure  13a from which it can be observed that the group settlement ratio predicted in this study is in good agreement with the prediction made by Pressley and Poulos [26] using finite element and elastic methods.

Pile Group Efficiency
Pile group efficiency (η) was calculated for the pile group for the three different spacings using the formulae proposed by Castelli and Maugeri [27]; McCabe and Lehane [28]; and Chellis [29] and the Converse-Labarre Formula. The pile group efficiencies calculated using the above methods and their comparison with the values given in the previous studies [26,30] are shown in Figure 13b. From the results (Figure 13b), it can be seen that the pile group efficiencies calculated using all three methods are in very good agreement with the results of the previous studies. It can also be seen that the pile group efficiency depends on the geometry of the pile group and that it increases with the spacing/diameter ratio.

Conclusions
The numerical simulation of the settlement behaviour of an axially loaded group pile located in deep silty-sand deposits was conducted using the PLAXIS 3D finite element package. From the findings of this study, the following conclusions can be drawn: 1. Linear elastic (LE) and Mohr Coulomb (MC) analyses underestimate the settlement of the pile group due to the ignorance of soil nonlinearity and simple nonlinearity respectively, while the hardening soil (HS) analysis overestimates the settlement of the