Winkler model for axially loaded piles in inhomogeneous soil

Analytical closed-form solutions are developed for the elastic and elasto-plastic settlement of axially loaded piles in inhomogeneous soil. The soil is modelled by way of a bed of Winkler (‘t–z’) s...


NOTATION
A pile cross-sectional area a stiffness inhomogeneity parameter b non-linearity exponent c u (z) undrained shear strength (variable with depth) D b pile base diameter D s pile shaft diameter E p pile elastic modulus G s (z) soil shear modulus at depth z I ν () modified Bessel function of the first kind of order ν K b pile base spring stiffness (unit of force per length) K el pile elastic section stiffness K ν () modified Bessel function of the second kind of order ν K 0 pile elastic head stiffness 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 z R k(z) Winkler modulus at depth z (unit of force per length 2 ) k 0 Winkler modulus at the ground surface L pile length L p length of plastic region m strength inhomogeneity exponent N c bearing capacity factor n stiffness inhomogeneity exponent P applied head load P b load transferred to pile base P u total ultimate resistance of pile P ub ultimate base resistance P(z) pile axial load at depth z q stiffness inhomogeneity parameter r m mobilised radius S n combination of Bessel functions t u (z) ultimate skin friction per unit length at depth z t u0 ultimate skin friction per unit length at the ground surface w b pile base settlement w by pile base yield settlement w(z) pile settlement at depth z w y (z) soil yield settlement at depth z w 0 pile head settlement z depth below ground level z R reference depth α pile-soil adhesion factor Γ() gamma function γ engineering shear strain γ 50 engineering shear strain at 50% of undrained shear strength λ load transfer parameter λ R load transfer parameter at reference depth z R ν order of Bessel function ν s soil Poisson's ratio τ soil shear stress χ 0 argument of Bessel function at pile head χ L argument of Bessel function at pile base Ω dimensionless base stiffness

INTRODUCTION
The Winkler model and related load-transfer analyses have been used extensively to model piles under axial load (e.g. Randolph & Wroth, 1978;Poulos & Davis, 1980;Scott, 1981;Mylonakis & Gazetas, 1998;Salgado, 2008;Guo, 2012). With careful selection of the Winkler modulus, the model has been shown to agree well with more rigorous continuum analyses (Mylonakis, 2001a(Mylonakis, , 2001bSyngros, 2004;Guo, 2012;Anoyatis, 2013;Anoyatis et al., 2018). Scott (1981) and Guo (2012) extended the Winkler model to inhomogeneous soil and provided closed-form solutions for pile head stiffness when soil stiffness varies as a power function of depth. This paper considers a similar stiffness variation to investigate the necessity of considering soil inhomogeneity when predicting pile settlement. New solutions are provided for the limiting cases of pile length, base conditions and soil inhomogeneity. The method is extended to the non-linear range by employing elastic-perfectly plastic Winkler springs, and applied to predict complete load-settlement curves for a series of tests on a site in London.

LINEAR RESPONSE
The pile model considered is depicted in Fig. 1. The base resistance is represented by a spring with stiffness K b and shaft resistance by uniformly distributed Winkler springs of stiffness k(z). The equilibrium of an arbitrary section of the pile yields the familiar governing differential equation (Scott, 1981) where E p and A denote the pile elastic modulus and cross-sectional area, respectively, and w is the pile settlement at depth z.
For the simple case of homogeneous soil, the Winkler modulus is constant with depth and the head stiffness, K 0 , of a pile of length, L, is given by the familiar solution (Mylonakis, 1995;Mylonakis & Gazetas, 1998;Randolph, 2003;Fleming et al., 2008;Salgado, 2008) where λ is a load transfer parameter (units of length −1 ) and Ω a dimensionless base stiffness constant This paper considers a variable Winkler modulus following the power-law function where k 0 is the Winkler modulus at the soil surface, k R the Winkler modulus at a reference depth, z R and a a dimensionless coefficient accounting for non-zero surface stiffness and n a positive inhomogeneity exponent. A simple approach to account for soil inhomogeneity is to calculate an average Winkler modulus, k = k av , over the whole pile length and use equation (2) to obtain the pile stiffness treating the soil as a homogeneous medium. To assess the validity of this approximation, one may compare against the exact solution to equations (1) and (4), given by equation 5 where I ν () is the modified Bessel functions of the first kind, of order ν. In the above equations, λ R , Ω R , q, ν, χ 0 and χ L denote the parameters ð6d; e; fÞ The full derivation of this expression is shown in section A1 of the Supplementary Appendix. Alternative forms of equation (5) have been presented by Scott (1981) and Guo (2012).
Equation (5) attains several different special forms which occur at the limits of the parameters a, n, L and Ω R , which are summarised in Table 1. Additionally, Tables 2 and 3 present solutions for normalised base displacement and base load, respectively, for different base conditions. The derivation of these solutions and alternative representations are provided in sections A2 and A3 of the Supplementary Appendix. Figure 2 shows normalised plots of pile head stiffness as a function of normalised pile length in which equation (2) is compared with equation (5). Seven Winkler modulus variations are considered by varying a and n and setting the reference depth to the pile length. Each variation is plotted for four dimensionless pile base stiffness values. As expected, when a = 1 and/or n = 0, k(z) is constant with depth and both solutions match. Naturally, the accuracy of equation (2) drops as k(z) deviates from a constant distribution. Figure 3 shows the percentage error of equation (2) compared with equation (5) as a function of dimensionless pile length for different inhomogeneous soils. If k 0 /k L > 0·5 and λL < 2, which covers a wide range of practical configurations, the discrepancy is less than 10% which indicates that the homogeneous soil approximation is acceptable. However, for k 0 /k L < 0·2, which applies to many normally consolidated soils, this approximation ceases to be acceptable as the error exceeds 30% beyond λL = 1·5 or so.

NON-LINEAR RESPONSE
As pile head load increases, the non-linear behaviour of soil gradually dominates response and the above solutions are no longer applicable. Scott (1981) and Guo (2012) propose Table 2. Summary of closed-form solutions for normalised pile base displacement in an inhomogeneous soil following a power-law stiffness variation with depth Perfectly end-bearing pile, Ω R ! ∞ 0 0 0 *For n = 0 these solutions reproduce those for homogeneous soil (in the first column). Dimensionless parameters λ R , Ω R , a, q, ν, χ 0 and χ L are provided in equations (4) and (6), Γ() = gamma function, S 2 = q n/2 I ν−1 (χ L ) + Ω R I +ν (χ L ), Table 1. Summary of closed-form solutions for pile stiffness in an inhomogeneous soil following a power-law stiffness variation with depth Mylonakis (1995) and Mylonakis & Gazetas (1998). †For n = 0 these solutions reproduce those for homogeneous soil (in the first column). Dimensionless parameters λ R , Ω R , a, q, ν, χ 0 and χ L are provided in equations (4) and (6), Γ() = gamma function, Table 3. Summary of closed-form solutions for normalised pile base load in an inhomogeneous soil following a power-law stiffness variation with depth *For n = 0 these solutions reproduce those for homogeneous soil (in the first column). Dimensionless parameters λ R , Ω R , a, q, ν, χ 0 and χ L are provided in equations (4) and (6), Γ() = gamma function, S 1 = q n/2 I 1−ν (χ L ) + Ω R I −ν (χ L ), a simple family of models to describe the non-linear shaft and base resistance of the pile, shown in Fig. 4. In this approach, linear elastic Winkler springs are replaced by linear elastic-perfectly plastic 't-z' curves, with elastic stiffness equal to k(z) and ultimate skin friction per unit pile length, t u (z), that can be established using traditional means (e.g. the α-method, Skempton, 1959 or the β-method, Burland, 1973). To derive a simple solution, the additional assumption is made that yielding originates at the pile top and propagates downwards. A sufficient condition to this end is that the displacement at yield at depth z, w y (z) = t u (z)/k(z), must be a non-strictly increasing function of depth. Finally, the base spring is assumed to behave linearly until the whole shaft friction is mobilised, transitioning to perfectly plastic behaviour beyond a base settlement, w by . The corresponding ultimate load can then be calculated using traditional bearing capacity approaches (e.g. Skempton, 1959;Salgado, 2008;Viggiani et al., 2011). Figure 4 shows the general case where the head load, P, is such that the shaft resistance has reached its limiting value, t u (z), over a certain length, L p . To calculate the corresponding settlement for a given applied load, L p must be determined, but this is not trivial as the associated equations are transcendental. Instead, L p can be taken as a known parameter and the rest of the solution, including forces and displacements, calculated as a function of L p . Repeating this for a series of L p values provides a complete load-settlement curve. In the realm of the above approach, the load-settlement response of the pile can be split into four stages, shown in Fig. 5 (a) Linear elastic response, L p = 0. Before the yield displacement at the top of the pile, w y (0), is reached, equation (5) can be used directly. (b) Non-linear shaft behaviour, 0 < L p < L. For an assumed plastic length, L p , the pile axial force at the plastic-elastic interface can be calculated by multiplying the yield displacement at that depth, w y (L p ), by the stiffness of the elastic pile section below, K el , from equation (5) with L p -L as the pile length. The plastic section can then be treated as a rod under known base and distributed shaft loads to obtain both the load and settlement at the pile head.
(c) Shaft resistance exhausted, L p = L. The behaviour of the pile is now governed by base response. A pile base settlement, w b , can be input into the load-settlement curve of the base to get a reaction. The shaft can then be treated as a rod with known distributed load. (d ) Base resistance exhausted, w b = w by . There is no additional resistance to load, therefore the applied head load, P, is equal to the ultimate resistance, P u , and head displacement increases without limit.
Considering a power-law variation in ultimate shaft resistance per unit length, t u (z), with exponent, m, and surface value t u,0 ; the load, P, and displacement, w 0 , at the pile head are given by the following parametric expressions where w(L p ), P(L p ) and t u (L p ) are the axial displacement, load and ultimate shaft resistance per unit length at the interface depth, z = L p , respectively. w(L p ) and P(L p ) depend on the response stage being considered, for stage 2, w(L p ) = w y (L p ) and P(L p ) = w y (L p )K el , where K el is the head stiffness of the elastic region of the pile, obtained from equation (5). For stages 3 and 4, the shaft resistance is exhausted, therefore, the interface depth is at the pile base, L p = L, and the load and displacement at this depth are given by the base load-settlement curve, w(L p ) = w b and P(L p ) = P b . Note that as the base spring is modelled as linear elastic, the pile head load, P, varies linearly with pile head settlement, w 0 , in stage 3. The derivation of these expressions is shown in section A4 of the Supplementary Appendix. Guo (2012) presents alternative forms of equations (7) and (8).
The above solution provides a closed-form alternative to the popular method by Coyle & Reese (1966), which, unlike the proposed one, is iterative. Additionally, the proposed solution does not require discretising the pile into elements, the number of which can affect the accuracy of the calculation in the iterative method. In fact, choosing different L p values only increases the number of points  (2) and the exact Winkler solution in equation (5). As equation (2) overestimates stiffness, the values in the graph are always positive used to plot the full load-settlement curve with each point being an exact result. However, discretising the pile into a sufficient number of elements and/or performing enough iterations will yield identical results to this solution at any required precision. It is therefore up to the designer/modeller to select their preferred method. Whitaker & Cooke (1966) conducted 11 maintained load tests on straight-shafted and under-reamed piles embedded in London Clay in the Wembley area. The pile dimensions are summarised in Table 4. The method described in this paper is applied here to predict the load-settlement response of these piles. To this end (a) Ultimate skin friction, t u (z), is estimated using the α-method (Skempton, 1959;Patel, 1992): t u (z) = αc u (z) πD s where D s is the shaft diameter, c u (z) the undrained shear strength, derived from triaxial tests on 38 mm samples by Whitaker & Cooke (1966) and reduced by a factor of 1·3 to account for the underestimation of fissured strength by tests on 38 mm samples (Patel, 1992) and α the pile-soil adhesion factor, taken here as 0·5 for London Clay (Patel, 1992;LDSA, 2017). (b) The ultimate load at the pile base, P ub , is obtained using Skempton (1951): P ub = N c c u (L) πD b 2 /4, where N c is the bearing capacity factor, taken here as 9 and D b the base diameter. (c) The Winkler modulus, k(z), is related to soil shear modulus, G s = G s (z), using the concentric cylinder theory of Cooke (1974) and Randolph & Wroth (1978): k(z) ≃ 2πG s (z) /ln(2r m /D s ), where ν s is the soil Poisson's ratio and r m an empirical radius beyond which displacement due to the load on the pile is presumed to reach zero. r m was estimated by Randolph & Wroth (1978) by matching Winkler analyses to numerical solutions of the full continuum problem, to get: r m ≈ 2·5L(1 − ν s )G s (L/2)/G s (L). Guo (2012) developed enhanced formulations for r m and identified some dependence on soil inhomogeneity. Mylonakis & Gazetas (1998) showed that the resulting Winkler modulus is rather insensitive to r m for piles of practical dimensions so the above expression is used here. (d ) The stiffness of the pile base spring, K b , is approximated as a rigid punch on a half-space (Randolph & Wroth, 1978):

APPLICATION TO PILE LOAD TESTS IN LONDON CLAY
where η = a factor to account for depth below the ground surface, taken as unity. (e) A relationship between maximum shear modulus, G 0 , and c u is obtained from Vardanega & Bolton (2011a), who analysed high-quality triaxial tests in London Clay and derived the correlation: G 0 ≈ 320c u , which is in meaningful agreement with the recommendations by Poulos (1989). ( f ) For shaft resistance and the early stages of base resistance, G 0 will adequately predict settlements, therefore it is used to derive the Winkler modulus and pile base stiffness for stages 1 and 2. However, as more base resistance is mobilised, G 0 will under-predict settlements, therefore a different model is required for stage 3 of the response. Vardanega & Bolton (2011b) proposed a simple power-law equation to describe soil shear stress-strain (τ-γ) response where γ 50 is the shear strain at 50% of the undrained shear strength and b is an exponent. Vardanega & Bolton (2011a) found that the tests on London Clay had an average b value of 0·58 and a γ 50 of about 7 × 10 −3 . A secant shear modulus between the estimated mobilisation shear strength of the base at the beginning of stage 3 and at failure interpreted from this model is used to predict the pile base stiffness for stage 3. The elastic portion of the base spring is, therefore, modelled as bilinear.
Predicted and measured load-settlement curves for the straight-shafted piles are shown in Fig. 6 and for the under-reamed piles in Fig. 7. The good agreement for small settlements indicates that the adopted shear modulus variation is realistic. Although the linear elastic-perfectly plastic model provides a sharper transition to yielding, the non-linear response is also approximated satisfactorily. In particular, the bilinear model for the pile base satisfactorily reproduces the response after shaft resistance is exhausted.

CONCLUSIONS
A simple analytical solution was derived for the analysis of elastic and inelastic settlement of an axially loaded pile in inhomogeneous soil obeying a power-law variation in stiffness and strength with depth. The solution makes use of the Winkler assumption and models soil-pile interaction using elastic and elastic-perfectly plastic springs distributed along the pile, and a concentrated spring at the pile base. The main findings of the study are as follows.
(c) Under the assumption of soil yielding propagating downwards, the solution was extended to the inelastic regime considering elastic-perfectly plastic 't-z' and base springs, expressed as a function of depth of the yield zone, L p . (d ) Using common site investigation data, such as c u (z), a database of high-quality tests in London Clay, and simple formulations to relate these properties to model parameters such as k(z), a successful comparison was presented against the classical measurements by Whitaker & Cooke (1966). The application of this method to a database of pile tests in London Clay can be found in Voyagaki et al. (2018).  Whitaker & Cooke (1966) As a final remark, the solution presented here can be applied to more general classes of elastic and elastoplastic Winkler beds (e.g. exponentially varying stiffness with depth, hyperbolic 't-z' and base resistance curves). It can also be extended to calculate interaction factors between piles in a group using the three-step method proposed in Mylonakis & Gazetas (1998). However, this lies beyond the scope of this study.