Settlement of axially loaded pile groups in inhomogeneous soil

Accurate prediction of settlement is key to performance-based design of pile groups. Simple methods based on physically motivated modelling assumptions, in conjunction with wisely chosen soil mater...


Introduction
Extensive research has been carried out to predict the response of single piles to an applied axial head load (summarised in the Rankine lectures by Poulos 1989and Randolph 2003, and recent books by Salgado 2008Viggiani, Mandolini andRusso 2011 andGuo 2012). However, piles are commonly installed as part of groups, where proximity to other piles reduces the stiffness of the foundation (Fleming, Weltman, Randolph and Elson 2009). Additionally, these pile-soil-pile interaction effects can lead to an uneven distribution of load between different piles in the group, as well as change the distribution of load down the pile, transferring more load to the base (Mylonakis and Gazetas 1998). Careful consideration of these effects is required to enable efficient, performance-based design.
A variety of simple models are available to predict pile group response, such as the equivalent raft or equivalent pier methods (Poulos 2006), where the pile group is replaced by a single foundation with empirically determined dimensions to be analysed in its place. The problem can be approached more rigorously using the finite element method (e.g. Ottaviani 1975), however a full three-dimensional model is required. As an alternative, Poulos (1968) introduced the concept of interaction factors. These employ superposition to allow the interaction between only two piles to be considered and applied to the whole group. This concept extended the available methods predicting the settlement of single piles to pile groups (Poulos 1968;Butterfield and Banerjee 1971;Banerjee and Davies 1977). Randolph and Wroth (1979) used a simple free-field displacement function describing soil settlement with radial distance from a single pile, to calculate interaction factors for use with a simpler load-transfer model. However, this does not account for the reinforcing effect of the second pile reducing the overall settlement. Mylonakis and Gazetas (1998) proposed a new model considering the response of an unloaded 'receiver' pile affected by the displacement of a loaded 'source' pile to derive an interaction factor. A closedform solution based on this method has been derived for piles embedded in homogeneous soils. However, soil properties, even in a single geological deposit, often vary with depth (Lumb 1966;Phoon and Kulhawy 1999). It is common practice to select a linearor more complex function of depth to describe this behaviour. This paper investigates piles embedded in inhomogeneous soils with stiffness varying according to a power-law function of depth to allow this simple model to be applied more rigorously to these situations. .
Step 1: Calculate the settlement profile for the loaded source pile, w 11 (z), with depth, z. .
Step 2: Define an attenuation function, ψ(s), describing soil settlement attenuation with radial distance, s, from the source pile. For simplicity the attenuation function is considered to be independent of depth. .
Step 3: Calculate the settlement profile for the receiver pile, w 21 (z), when subjected to the attenuated settlement profile of the source pile.
The problem is analysed using the Winkler model. The piles are represented by elastic columns with Young's modulus, E p , and cross-sectional area, A. The shaft resistance due to the soil is modelled using uniformly distributed Winkler springs with stiffness k(z), and the base resistance is modelled by a single spring with stiffness K b . The behaviour of the source pile is described by the well-known governing differential equation (Scott 1981): The displacement field surrounding the loaded source pile due to the shaft displacement can be approximated using the concentric cylinder model (Cooke 1974;Randolph and Wroth 1978). By assuming the variation with depth of the vertical soil deformations is negligible compared to the variation with distance from the pile (i.e. the response is dominated by shear deformations), the vertical displacement is found to follow a logarithmic function with radial distance from the pile. Randolph and Wroth (1979) proposed using the attenuation function given by equation 2.
where d is the diameter of the source pile and r m is an empirical radius beyond which the displacement due to the load on the source pile is presumed to reach zero. r m can be estimated by matching equation 2 to experimental or numerical results. Randolph and Wroth (1978) developed the following approximate expression for r m by matching equation 2 to finite-element data: where L is the pile length, ν s is the Poisson's ratio of the soil and ρ is an inhomogeneity parameter given by the ratio of the average soil shear modulus to the soil shear modulus at the pile base. This equation predicts a typical range of r m from 5-25 pile diameters (usually about 20 diameters for common pile lengths). The attenuation of settlement due to the displacement of the pile base can be well approximated as following that due to a rigid punch on the surface of a half-space (Randolph and Wroth 1979;Mylonakis and Gazetas 1998), which decreases inversely proportionally with radial distance from the pile. This decreases at a much higher rate than the displacement field due to the shaft displacement and is negligible for almost all practical situations.
Due to the influence of the source pile, the governing differential equation of the receiver pile is given by: The response of the receiver pile at the pile head, w 21 (0), normalised by the response of the source pile at the pile head, w 11 (0), yields the desired interaction factor, α: a = w 21 (0)/w 11 (0) (5)

Homogeneous soil
If the Winkler modulus is taken constant with depth, k(z) = k, the single pile head stiffness, K 1 , is given by (Scott 1981 where λ is the load-transfer parameter with units Length −1 and Ω is a dimensionless base stiffness constant: For the most common case of two identical piles the interaction factor is calculated by solving equation 4 and applying the boundary condition of zero force at the head and a spring with modulus K b at the base of the receiver pile (Mylonakis 1995;Mylonakis and Gazetas 1998): where ζ is a so-called diffraction factor representing the reinforcing effect of the receiver pile on the surrounding soil. The variation of the diffraction factor, ζ, with dimensionless pile length, λL, for different base conditions is shown in Fig. 2. Evidently, as the pile length increases the interaction factor asymptotically approaches a common value independent of the base conditions, given by: Therefore, the interaction factor is only one half of the freefield value when the receiver pile is very long. Equation 8 has been included in foundation design textbooks (e.g. Salgado 2008; Fleming, Weltman, Randolph and Elson 2009) and has been applied in pile group design software such as PIGLET (Randolph 2006). However, it is limited to homogeneous soils.

Inhomogeneous soil
In inhomogeneous soils, exact Winkler solutions for the response of single piles are available for a limited number of soil stiffness variations with depth. Closed form solutions for soil stiffness varying according to a power-law function of depth are provided by Scott (1981), Guo (2012) and Crispin, Leahy and Mylonakis (2018). This paper investigates the following Winkler modulus variation: where k 0 and k R are the Winkler modulus at the ground surface and a reference depth, z R , respectively, a is an inhomogeneity parameter accounting for non-zero stiffness at the ground surface, and n is an inhomogeneity exponent.
The general solution to equation 1, when k(z) is described using equation 10, is given by (Crispin, Leahy and Mylonakis 2018): (11) where I ν (χ) and K ν (χ) are the modified Bessel functions of the first and second kind, respectively, of order ν and argument χ, while λ R is analogous to the parameter λ in equation 7.
The head stiffness of the single pile is given by (Crispin, Leahy and Mylonakis 2018): where the dimensionless factors are given by: By solving equation 4 when k(z) is described using equation 10 and applying the boundary conditions at the head and base of the source and receiver pile, the following closed-form expression for the interaction factor is 2 Variation of diffraction factor, ζ = α/Ψ(s) with dimensionless pile length, λL (adapted from Mylonakis and Gazetas 1998, used with permission from ICE publishing) Similarly to equation 8, with increasing pile length the above equation tends to a common value independent of the base conditions; in this case they approach the result in equation 16, which expresses the behaviour of an infinitely long pile. This reduces to α = ν·ψ(s) for a = 0, which, in turn, reproduces equation 9 when n = 0.
Additionally, if the parameters in equation 10 are chosen to represent a homogeneous soil (n = 0 or a = 1), equation 15 duly reduces to equation 8. Equation 15 is suitable for implementation in a software package. However, due to the presence of the Bessel functions, it is inconvenient for use in hand calculations. As an alternative, equation 8 could be applied using an equivalent homogeneous stiffness profile. An intuitive choice to this end is to take the average Winkler modulus, k av , over the pile length. In this case, if the pile length is chosen as the reference depth for equation 10 (i.e. z R = L, k R = k L , q = 1): Figure 3 shows the interaction factor in equation 15, as a function of dimensionless pile length for inhomogeneous soils. Equation 15 is compared to equation 8 using the parameters in equation 17. Six soil stiffness profiles are considered. The equivalent homogeneous approximation shows good agreement for short piles. However, the prediction diverges as the dimensionless length increases. Naturally, the accuracy of the approximation decreases the further k(z) deviates from a constant value (i.e. by increasing n and decreasing a). Figure 4 shows the percentage error due to using the above approximation to the full solution in equation 15, as a function of dimensionless pile length when Ω = 1. The trends noted above are more evident here: the error increases with pile length and is larger the further the stiffness profile diverges from homogeneous. For a Gibson soil (a = 0, n = 1), the error is greater than 20% when λL = 0.6 (corresponding to L/d ≈ 20 for a pile-soil stiffness ratio of 10 3 ). As the approximation overestimates the interaction factor, it would lead to over-conservative designs. How the magnitude of this error translates to the predicted performance of a specific foundation depends on how much the piles in a foundation interact. The error will increase for groups with more piles and closer spacings.
In order to improve the predictions of the approximate method, a model error correction factor, η, was introduced to correct the behaviour as the pile length increases. The factor can be calculated by taking the ratio of the asymptotic behaviour of the analytical solution in equation 16, to that of the homogeneous approximation using equation 9. An exponent depending on pile length has been included to reduce the effect of this factor when the pile is shorter (and the error is lower). The resulting expression is given in equation 18 and the factor is plotted in Fig. 5.

Illustrative example
A group of 4, 0.6 m diameter, 15 m long concrete piles is installed in a normally consolidated clay. The piles are arranged in a square with 1.8 m centre-to-centre spacing and connected by a rigid cap. Ground investigation data from the site showed that the soil stiffness variation with depth can be approximated using a linear profile (n = 1) having zero surface stiffness (a = 0). Back analysis of a load test on a single pile on the same site showed that at the load level of interest, the shear modulus can be modelled as increasing by 2.5 MPa/meter. The problem is shown in Fig. 6. The Winkler modulus, k(z), can be related to the soil shear modulus, G s (z), using the same concentric cylinder model used to approximate the attenuation function (Randolph 2003): Employing equations 3 and 19 for this example gives k(z) = 4.6 z MPa. Setting the reference depth as the pile length (z R = L, k R = k L , q = 1) and assuming a Young's modulus of concrete of 20 GPa yields λ R z R = 1.65. Modelling the pile base as a rigid punch on the surface of a half-space (Randolph and Wroth 1978; Mylonakis and Gazetas 1998) yields a dimensionless base stiffness, Ω R of 0.14.
The settlement ratio, R s , is the ratio of the average settlement of the pile group to the settlement of a single pile at the same average load as the pile group (Poulos and Davis 1980). Calculating this ratio allows the group settlement to be estimated by multiplying it with the predicted settlement of a single representative pile. As the pile cap is rigid (all piles must settle the same amount) and the group is symmetric (each pile is influenced by the same set of piles, the same distance apart), the settlement ratio of a group of m piles can be calculated using equation 20: where α ij is the interaction factor derived from the response of pile i due to a load on pile j and α ii = 1. In this case there are only two interaction factors that need calculating, α 1 for the piles two piles adjacent to pile i and α 2 for the pile opposite. Table 1 shows the calculated settlement ratio using the analytical inhomogeneous solution, the equivalent homogenous approximation and the corrected approximation. Evidently, the pile group settles over 150% the amount of a single pile under the same average load. Note that the homogenous soil approximation underestimated R s, relative to the more accurate solution in equation 15 by 9% in this case, which could lead to a more 3 Variation of diffraction factor, ζ = α/Ψ(s), with dimensionless pile length, λL, for different degrees of soil inhomogeneity conservative design, particularly for larger groups likely to be encountered in practice. On the other hand, the corrected approximation overestimated R s by less than 2%.
Comparison to experimental data Koizumi and Ito (1967) reported the results of a load test on a group of 9 piles. The test was preceded by a test on an identical individual pile in order to isolate the interaction response of the group. The piles were 0.30 m diameter steel tubes with wall thickness 3.2 mm, embedded 5.55 m into in a soft silty clay. The group is arranged in a 3 × 3 square with a spacing of 0.9 m and connected with a rigid cap. The test arrangement is shown in Fig. 7. The tests have previously been analysed by Poulos (1974) and Randolph and Wroth (1979). Both authors back-analysed the single pile test to derive a soil shear modulus profile with depth for use in analysing the pile group, Poulos (1974) assumed a homogeneous profile and subsequently derived a shear modulus, G s = 5.8 MPa. Randolph and Wroth (1979) assumed a linear variation with depth, resulting in a shear modulus variation of the type: G s (z) = 7·z/L MPa. These were then used to predict the load distribution between the piles in the group.
Employing the linear shear modulus variation derived by Randolph and Wroth (1979), the analytical solution for inhomogeneous soils described in this paper has been used to analyse this test. Setting the reference depth to the pile length, the following parameters describe the problem: a = 0, n = 1, q = 1, λ R z R = 0.83. Using the same base stiffness assumptions as the previous example (resulting in Ω R = 0.09) and equations 2 and 15 (which yields ζ = 0.68), the interaction factors between each pair of piles have been calculated. The settlement of any pile in the group, w i , can then be calculated using equation 21 (Poulos and Davis 1980).
where 1/K 1 is the settlement of a single pile under unit load and P j is the load carried by pile j. As the pile cap is rigid, the settlement of each pile is the same, therefore equation 21 has been employed to generate a set of simultaneous 5 Variation in Inhomogeneity correction factor, η = 2 ζ(L 1) with soil inhomogeneity parameters a and n; z R is chosen such that λ R z R = 1 6 Illustrative example dimensions equations relating the loads in each pile to the group settlement. Due to the symmetry of the foundation, this reduces to three simultaneous equations relating to the loads in the centre pile, the corner piles and the piles in the middle of each side. By considering the total load carried by the group, these equations have been solved and the proportion of the total load carried by each pile calculated. The results are shown in Table 2 with comparison against the measured test results (Koizumi and Ito 1967) and the predictions by Poulos (1974) and Randolph and Wroth (1979). When the load on the group was 910kN, a settlement of 7.1 mm was measured (Koizumi and Ito 1967). Inputting the results in Table 2 back into equation 21 and using equation 13 to calculate the single pile head stiffness yields an estimated settlement of 6.7 mm at this load.

Summary and conclusions
A closed-form analytical solution was derived for the interaction factor between piles in a group embedded in an inhomogeneous soil obeying a power-law variation in stiffness with depth (equation 15). The proposed solution employs the three-step model proposed by Mylonakis and Gazetas (1998), and the Winkler assumption. The method is suitable for implementation in a design-oriented software package and design charts have been provided (Fig. 3) for use in hand calculations. The analytical solution has been compared to an approximation considering an equivalent homogeneous soil and, due to the complexity of the analytical solution, a simple model error correction factor developed to improve the accuracy of the approximate method for routine use (equation 18, Fig. 5). This correction factor reduces the discrepancy between this approximate method and the analytical solution presented here to below 10% for most practical configurations. An illustrative example of the application of interaction factors to a simple design problem has been provided. Finally, the proposed solution has been shown to have good agreement with both experimental results and other analytical methods.

Notation
Latin symbols a stiffness inhomogeneity parameter A pile cross sectional area C n integration constant to be determined from the boundary conditions d pile diameter E p pile elastic modulus G s (z) soil shear modulus variation with depth I ν ( ) modified Bessel function of the first kind K ν ( ) modified Bessel function of the second kind K 1 single pile head stiffness K b pile base spring stiffness k(z) Winkler modulus variation with depth k 0 Winkler modulus at the surface k av average Winkler modulus over the pile length k L Winkler modulus at the pile base k R Winkler modulus at the reference depth L pile length n stiffness inhomogeneity exponent P j load carried by pile j q stiffness inhomogeneity parameter R s settlement ratio r m 'magical' radius S n combination of Bessel functions s pile separation distance w 11 (z) source pile settlement profile due to applied head load (variable with depth) w 21 (z) receiver pile settlement profile due to load on source pile (variable with depth) w i head settlement of pile i z depth below ground level z R reference depth Greek symbols α pile-soil-pile interaction factor α ij interaction factor describing the response of pile i due to a load on pile j ζ diffraction factor 7 Koizumi and Ito (1967) pile group test arrangement Table 2 Measured and predicted load distribution in Koizumi and Ito (1967)  η model error correction factor λ load transfer parameter λ R load transfer parameter at reference depth ν order of Bessel function ν s soil Poisson's ratio ρ inhomogeneity parameter in Wroth (1978, 1979) theory χ 0 argument of Bessel function at pile head χ L argument of Bessel function at pile base ψ(s) free field displacement attenuation function Ω dimensionless base stiffness Ω R dimensionless base stiffness in inhomogeneous soil