Finite-element modelling of laterally loaded piles in a stiff glacial clay till at Cowden Open Access

The design of laterally loaded piles in clay is a complex problem, involving aspects of soil–structure interaction which are difficult to consider in most methods of analysis, such as the possibility of a gap opening on the active soil side. However, despite these challenges, there are various published studies investigating the behaviour of laterally loaded piles, with solution strategies including elastic methods (e.g. Poulos & Davis, 1980; Basu et al., 2009) and numerical approaches, such as those based on the finite-element (FE) method (e.g. Brown & Shie, 1990; Yang & Jeremić, 2002; Jung et al., 2015; Haiderali & Madabhushi, 2016). In these studies, the response of clay is typically simulated using a total stress approach with a Tresca (e.g. Jung et al., 2015) or a von Mises failure criterion (e.g. Brown & Shie, 1990; Yang & Jeremić, 2002). The undrained shear strength, $Su$⁠, is often assumed to be either uniform (e.g. Yang & Jeremić, 2002) or increasing with depth (e.g. Jung et al., 2015), while the associated Young's modulus is specified as being uniform (e.g. Brown & Shie, 1990), a function of $Su$ (e.g. Jung et al., 2015), or to depend directly on the confining pressure (e.g. Yang & Jeremić, 2002). Although the abovementioned studies provide significant insight into the response of laterally loaded piles in clay, particularly in terms of ultimate limit states, the accurate prediction of displacements requires the use of more comprehensive modelling approaches, capable of simulating accurately the non-linear pre-yield response of the soil, including the dependency of stiffness on the stress and strain levels (see e.g. Jardine et al. (1986) for a discussion on the impact of this aspect of soil response on the behaviour of geotechnical structures). In this paper, the procedure adopted in the advanced three-dimensional (3D) FE analyses carried out in the context of the PISA project, which include the abovementioned features, is presented and discussed.

As outlined in the paper by Zdravković et al. (2019), the PISA (Pile–Soil Analysis) project proposed a new design model for laterally loaded piles. Referred to as the PISA design model, it is consistent with the existing one-dimensional (1D) Winkler-type p–y approach, but extended to include soil reactions in addition to that representing just a distributed lateral load. The principal premise of this development was that the new 1D PISA design model could be derived from the results of site-specific 3D FE modelling. Consequently, the 3D FE analyses discussed in the current paper were not back-analyses of pile tests. Instead, they were performed before the medium-scale piles, installed in a glacial clay till at the Cowden test site, were tested under lateral loading, thus providing Class A predictions according to Lambe (1973). The results of these analyses were used to: (a) initially aid the design of the field testing programme (Byrne et al., 2019); (b) subsequently validate the established numerical model by comparison with field testing results; and (c) provide insight into the main load-transfer mechanisms of laterally loaded monopiles in glacial clay, which would form the basis for formulating soil-reaction curves for new Winkler-type design models. Similar field testing (McAdam et al., 2019) and numerical (Taborda et al., 2019) studies were performed on piles installed in a dense marine sand at the Dunkirk test site.

The central focus of the PISA project was to develop a design model for monotonic loading only (see Zdravković et al., 2019), as a key step for future extensions of the model to cyclic and other loading conditions. Consequently, the 3D FE model and most of the field testing programme were developed with this objective. Some field testing explored the effects of different loading rates and of load cycles, to maximise the output from the field tests, but these aspects were not considered at the current stage of the new design model development.

The paper starts with elaborating a consistent 3D numerical model for laterally loaded piles installed at the Cowden test site, based on available ground information interpreted and discussed in the paper by Zdravković et al. (2019). This is followed by demonstrating the capabilities of the developed numerical model through comparison of the predicted response with the field measurements presented in the paper by Byrne et al. (2019).

The interpretation of the Cowden ground model discussed in the paper by Zdravković et al. (2019) demonstrated the applicability of the critical state framework to describe the mechanical behaviour of the low-plasticity, overconsolidated clay till. The constitutive model employed for the numerical analyses presented here is selected to (a) satisfy design requirements, as outlined in the paper by Zdravković et al. (2019), for realistic representation of initial and ultimate loading behaviour of the test piles, and (b) simulate accurately aspects of the mechanical behaviour of Cowden till, observed from the available historic and new field and laboratory investigations. The model adopts elements of the original modified Cam Clay model formulation (Roscoe & Burland, 1968), but is enhanced with: (a) a non-linear Hvorslev-type surface (Tsiampousi et al., 2013) to capture accurately the undrained strength of Cowden till dry of critical; (b) a generalised shape for the yield and plastic potential surfaces in the deviatoric plane (Van Eekelen, 1980) to account for the effect of the intermediate principal stress, as evidenced by the different strength of Cowden till in triaxial compression and extension; and (c) a modified hyperbolic equation for simulating the non-linear variation of the elastic shear modulus with mean effective stress and deformation level (Taborda & Zdravkovic, 2012; Taborda et al., 2016). The latter component of the numerical model can be coupled with a reversal-detection algorithm and non-linear scaling laws (see Taborda et al., 2016) in order to replicate the hysteretic response exhibited by geomaterials when subjected to cyclic loading. However, given that the objective of the PISA project was to establish the response of piles subjected to monotonic lateral load, this aspect of the constitutive model formulation was deemed unimportant to the analyses presented herein and was deactivated.

It is noted that the selected constitutive model does not reproduce cyclic and strain-rate effects, as this was not the objective of the project, nor does it reproduce strength anisotropy of the soil, as characterisation of these aspects of Cowden till behaviour was not possible from the available experimental information.

The Hvorslev surface, presented in detail by Tsiampousi et al. (2013), introduces alterations to the shape of the yield and plastic potential surfaces in the $p′$–J plane (where $p′$ is mean effective stress and $J$ is generalised deviatoric stress, see Appendix 1) on the dry (or supercritical) side of the critical state, thus defining the elastic region for overconsolidated clays and controlling their plastic response. The adopted yield surface, $F(σ′,p′0)$⁠, is defined by

where $θ$ is the Lode's angle (see Appendix 1); $g(θ)$ defines the inclination of the critical state line (CSL) to the p′ axis for the current value of $θ$⁠; $α$ and $n$ are parameters defining the non-linearity of the yield surface on the dry side (i.e. $p′<p′0/2$⁠); and $p′0$ is the pre-consolidation mean effective stress and the state/hardening parameter for this model. Fig. 1 illustrates the shape of the yield surface defined by equation (1) in the $p′$–J plane. The plastic potential, $P(σ′,p′0)$⁠, is modified accordingly

where $β$ and $m$ are parameters defining the shape of the plastic potential and therefore control the rate at which the material reaches the critical state, while $p′c$ and $θc$ denote the current mean effective stress and Lode's angle, respectively (see Tsiampousi et al., 2013). As shown in the expressions above, on the wet side of critical (i.e. $p′≥p′0/2$⁠) the yield surface and plastic potential are similar to the original proposal by Roscoe & Burland (1968), where associated plasticity was assumed.

As shown in Fig. 1, the CSL inclination is controlled by the function $g(θ)$⁠, which defines the shape of the yield surface in the deviatoric plane (see Potts & Zdravković (1999) for a detailed discussion on the importance of this function). In the present case, the generalised shape proposed by Van Eekelen (1980) was employed

where $X$⁠, $Y$ and $Z$ are model parameters, the values of which need to satisfy a number of conditions to ensure the surface is convex. The use of this shape takes account of the effect of the intermediate principal stress by allowing the simulation of the shear strength variation in the deviatoric plane, such as that observed in Cowden till where the strength in triaxial compression is different from that in triaxial extension loading (Zdravković et al., 2019). The Van Eekelen (1980) function, with the appropriate choice of parameters, can reproduce a variety of deviatoric shapes, including the Matsuoka & Nakai (1974), or Lade & Duncan (1975) surfaces.

The final modification to the formulation is the introduction of a non-linear variation for the tangent shear modulus, which is based on the modified hyperbolic expression used in the ICG3S model described in the papers by Taborda & Zdravkovic (2012) and Taborda et al. (2016). According to this, the tangent non-linear elastic shear and bulk moduli are given by

where $G0∗$ and $K0∗$ are the basic stiffness properties of the material (i.e. without taking into account the effects of void ratio, mean effective stress and strain level); $p′ref$ is a reference pressure; $fG(e)$ and $fK(e)$ are user-defined functions expressing the influence of void ratio on the values of the shear and bulk moduli, respectively; $mG$ and $mK$ are parameters controlling the non-linearity of the variation of stiffness with mean effective stress; $Ed$ is the generalised deviatoric strain (i.e. second invariant of the strain tensor, see Appendix 1); $εvol$ is the volumetric strain (i.e. first invariant of the strain tensor, see Appendix 1); and $a,b,RG,min$⁠, $r$⁠, $s$ and $RK,min$ are parameters defining the reduction in stiffness with deformation level. Clearly, there is an incompatibility between the ICG3S model and the original formulation of the modified Cam Clay model, since the assumption that the swelling lines have a constant slope (⁠$κ$⁠) in the plane plotting the specific volume (⁠$v$⁠) against $ln⁡p′$ implies that the value of the tangent elastic bulk modulus is

Therefore, the use of this small-strain stiffness overlay is only possible if the adopted function of void ratio is defined as $fKe=1+e=v$ and $mK$ and $RK,min$ are set to 1·0 (the latter parameter renders the values of $r$ and $s$ irrelevant). In this situation, the ICG3S model will reproduce swelling lines of constant slope, and hence be compatible with the original formulation of the modified Cam Clay, providing the value of $K0∗$ is chosen to be

As a result, the integrity of the assumptions central to the modified Cam Clay model is guaranteed, meaning that, for example, the closed-form solutions for the undrained shear strength (Potts & Zdravković, 1999) remain unchanged. Concurrently, the use of the ICG3S model enables a more accurate representation of the considerably non-linear variation of the shear modulus exhibited by Cowden till (Zdravković et al., 2019), which would not be captured by the alternative formulations for the elastic shear modulus typically employed when using the modified Cam Clay model (constant Poisson ratio, $ν$⁠, constant shear modulus, $G$⁠, or constant ratio $G/p′0$⁠).

The interpretation of the historic and newly performed laboratory tests on Cowden till in terms of strength at critical state under triaxial compression and triaxial extension loading conditions is discussed in detail in the paper by Zdravković et al. (2019). These data indicate angles of shearing resistance of $ϕ′TXC=27deg$ and $ϕ′TXE=32deg$⁠, in triaxial compression and extension, respectively. Since no information on the strength of the material is available for intermediate values of the Lode's angle, a number of combinations are possible for the parameters $X$⁠, $Y$ and $Z$ in equation (3), ensuring the convexity of the yield surface. To achieve the experimentally observed difference of up to 6° in the angle of shearing resistance between triaxial and plane-strain loading conditions for most soils (e.g. Bishop, 1966; Gens, 1982), an iterative procedure suggested that it was appropriate to assume a value of 0·1 for $Z$⁠, thus allowing the calculation of the remaining quantities using

In the expressions above, $MJ,TXC$ and $MJ,TXE$ designate the ratios $J/p′$ at critical state for triaxial compression (⁠$θ=−30o$⁠) and triaxial extension (⁠$θ=30o$⁠) loading conditions, respectively. These can be calculated using (Potts & Zdravković, 1999)

For the calibrated parameters, the variation of the angle of shearing resistance with the value of Lode's angle is illustrated in Fig. 2(a). Note that higher values of $Z$ would have increased substantially the simulated strength of the material for non-triaxial loading conditions, which would be unrealistic.

The characterisation of the strength of the material requires the calibration of the parameters controlling the shape of the non-linear Hvorslev surface. These were determined by plotting the measured effective stress paths, of highly overconsolidated samples sheared under undrained triaxial compression, in the plane $p′/p′0$–$J/p′$ (Fig. 2(b)) and performing a non-linear regression, which gave $α=0⋅25$ and $n=0⋅40$⁠. The calibrated non-linear Hvorslev surface is plotted together with the normalised stress paths in Fig. 2(b), demonstrating good agreement between the measured and simulated yield loci.

The results of the oedometer tests, described in the paper by Zdravković et al. (2019), indicate that the consolidation behaviour of Cowden till, particularly that of the more superficial layers which are dominant for shorter piles, can be represented with good accuracy by $ν1=2⋅20$ and $λ=0⋅115$⁠. The swelling response, as measured in the same tests, is characterised by $κ=0⋅021$⁠, thus determining the tangent bulk modulus (equation (6)). Since the latter is linearly dependent on the value of mean effective stress, a similar variation with p′ was assumed for the elastic shear modulus, meaning that $mG=1⋅0$ in equation (4).

A summary of field and laboratory measurements of the Cowden till stiffness is presented in the paper by Zdravković et al. (2019), including results from the new seismic cone penetration tests (SCPTs), bender element (BE) tests and triaxial compression (TXC) and triaxial extension (TXE) tests using high-resolution linear variable differential transducers (LVDTs), as well as from historic cross-hole (CH) and down-hole (DH) geophysics. As discussed in the paper by Zdravković et al. (2019) and shown in Fig. 3(a), the dynamic measurements SCPT, BE and DH indicate considerably higher stiffness values than those obtained using local transducers (TXC, TXE), which may be attributed to insufficient resolution of the latter to capture elastic stiffness at very small strains. Given the consistency exhibited by the dynamic measurements, a value of $G0/p′=1100⋅0$ was assumed representative of the small-strain shear modulus. As discussed in the paper by Zdravković et al. (2019), a much higher stiffness profile from the historic CH data might be attributed to soil variability at the previous Cowden test area (described in the paper by Powell & Butcher (2003)), which was substantially larger than the PISA test area and located in a different part of the Cowden site (see Zdravković et al., 2019). Adopting a reference pressure, $p′ref$⁠, of 100 kPa and discarding the effect of void ratio on the value of the elastic shear modulus due to lack of clear experimental evidence (i.e. $fGe=1⋅0$⁠), implies that parameter $G0∗$ in equation (4) must assume a value of 110 MPa.

The measured variations of the normalised secant shear modulus, $Gsec/p′$⁠, with strain level, $Ed$⁠, for individual triaxial tests are shown in Fig. 3(b), together with the simulated behaviour obtained using equation (4), with parameters $a=9⋅78×10−5$⁠, $b=0⋅987$ and $RG,min=0⋅05$⁠. As seen in Fig. 3(a), there are a number of triaxial tests for which the maximum stiffness is below that given by the assumed relationship $G0=1100p′$⁠. Therefore, the values chosen for $a$⁠, $b$ and $RG,min$ in most cases imply higher-than-measured stiffness for relatively small strains. To mitigate the impact of this overestimation, the non-linear elastic part of the constitutive model was calibrated to predict a considerably sharper reduction in stiffness with increasing deviatoric strain than observed in most tests, such that at moderate strains the simulated behaviour is closer to the average stiffness measured in all the triaxial tests. This is shown more clearly in the normalised $Gsec/G0$ degradation with deviatoric strain in Fig. 3(c), further highlighting the need to consider globally the stiffness during calibration, rather than focusing on producing, in isolation, the maximum stiffness (Fig. 3(a)) and its normalised reduction curve (Fig. 3(c)).

The final parameters to be determined are those controlling the shape of the plastic potential function associated with the non-linear Hvorslev surface, $β$ and $m$⁠. These were evaluated using a trial-and-error approach, focusing on the rate at which overconsolidated samples approach the critical state, as discussed in Tsiampousi et al. (2013). In the present case, a simple linear form for the plastic potential function was assumed (i.e. $m=1⋅00$⁠), with a value of $β=0⋅20$ appearing to reproduce the results of triaxial tests with greater accuracy, as shown in Appendix 2. The resulting set of calibrated parameters is listed in Table 1, while the performance of the model is briefly illustrated by comparisons with experimental data in Appendix 2.

All analyses reported herein were carried out with the Imperial College Finite Element Program (ICFEP) (Potts & Zdravković, 1999). Given the existence of a plane of symmetry in the problem of a laterally loaded pile, only half of the domain was discretised, as shown in an example FE mesh in Fig. 4. In total, 14 test piles were installed at Cowden (Byrne et al., 2019), with four distinct geometries, as summarised in Table 2, representing medium- (⁠$D=0⋅762m$⁠) and large (⁠$D=2⋅0m$⁠)-diameter piles that were the focus of numerical predictions in the current paper. Multiple tests were conducted on the same pile geometries to confirm the repeatability of the field measurements.

The soil domain was discretised using 6464 20-noded hexahedral displacement-based isoparametric solid elements. The far vertical boundary was positioned at a radial distance of 30 m for analyses CM2, CM3 and CM9 (⁠$D=0⋅762m$⁠), while for analysis CL2, where a pile of larger diameter (2·0 m) was simulated, this was increased to 80 m, mitigating any possible boundary effects. In all analyses the depth to the bottom boundary was 40 m, at the interface of the Cowden till deposit with chalk. The tubular open-ended pile is discretised with 280 eight-noded shell elements (Schroeder et al., 2007), of which 18 rows of elements in the $Z$-direction discretise the embedded part of the pile and ten rows the stickup. The steel was assumed to be linear elastic, with a Young's modulus of $E$ = 200 GPa and a Poisson ratio, $ν$ = 0·30. Furthermore, the pile was assumed to be wished in place, hence disregarding any installation effects. This was considered acceptable as a laterally loaded pile derives its capacity from the mobilisation of a significant volume of soil around the pile, well beyond the pile–soil interface zone that may be disturbed by installation. Consequently, the effect of the disturbed zone is thought to be less significant compared to an axially loaded pile, which relies on interface conditions for its capacity. As shown later, for lateral loading it is much more important that the interface can simulate opening of a gap around the pile, as this would affect its capacity. Consequently, the outer pile–soil interface is discretised using 180 16-noded zero-thickness interface elements (Day & Potts, 1994), the behaviour of which is represented by an elasto-plastic Tresca model with elastic normal (⁠$KN$⁠) and shear stiffness (⁠$KS$⁠) values of $1⋅0×105kN/m3$ and a limited tensile capacity (in the present case, assumed to be zero). The interface shear strength in compression, denoted as $Su,int$⁠, is defined using an innovative approach where the interface element adopts the initial undrained shear strength of the adjacent soil element, thus seamlessly reflecting the contributions of the spatial variations in stress conditions, overconsolidation ratio (⁠$OCR$⁠) and loading direction (⁠$θ$⁠). For the modified Cam Clay model, an analytical solution for the undrained shear strength is utilised (Potts & Zdravković, 1999)

where $OCR$ is determined using $σ′vm/σ′v0$⁠, with $σ′vm$ being the maximum vertical effective stress that the soil has been subjected to and $σ′v0$ being the initial vertical effective stress, $K0OC$ is the current value of the coefficient of earth pressure at rest, while $K0NC$ is its value associated to normally consolidated conditions and calculated as $(1−sin⁡ϕ′cs)$⁠. Lastly, $B$ is given by

To avoid rigid body movements of the mesh, the displacements along all directions (⁠$X$⁠, $Y$ and $Z$⁠) are set to zero over the bottom boundary of the mesh (i.e. at $Z=−40⋅0m$⁠). Moreover, the displacements normal to the vertical cylindrical boundary are prescribed as zero. To ensure that the X–Z plane at $Y=0$ is a plane of symmetry, the displacements in the Y-direction over this plane are set to zero. Similarly, at the nodes of shell elements, defining the edges of the pile contained in this plane, the rotational degrees of freedom about the X- and Z-axes are also set to zero. The load applied at the top of the pile (i.e. $Z=h=+10m$⁠) is simulated by imposing an incremental uniform horizontal displacement in the X-direction to all the nodes defining the perimeter of the section. The calculated reactions at each node are then summed up to give the total load $H/2$⁠.

A detailed description of the ground conditions is provided in the paper by Zdravković et al. (2019), while herein only information relevant to the numerical analyses is presented. In terms of initial conditions, the total vertical stress is obtained by assuming a saturated bulk unit weight, $γsat$⁠, of $21⋅19kN/m3$⁠, which is consistent with the values of density reported by Powell & Butcher (2003) and Ushev et al. (2015). The adopted variation with depth of the at-rest coefficient of earth pressure, $K0$⁠, is shown in Fig. 5(a), with a profile that follows the Powell & Butcher (2003) data and imposes the shallow cut-off of 1·5, as discussed in the paper by Zdravković et al. (2019). The pore water pressures measured in situ revealed an under-drained distribution in the top 12 m, Fig. 5(b), as a result of the soil permeability, $k$⁠, reducing considerably with depth and the presence of a lower aquifer with a phreatic surface below the groundwater table (Zdravković et al., 2019). The profile adopted in the analyses (solid line) is derived at the initiation of the ground conditions, using an exponential law for the permeability variation with mean effective stress (see Vaughan, 1994; Potts & Zdravković, 1999), fitted to reproduce the reported reduction in permeability while maintaining pore water pressures in equilibrium with the water table at 1 m depth and the values estimated in the sand layer. The resulting non-linear under-drained pore water pressure profile exhibits good agreement with field measurements. Below the sand layer the available data on pore pressures are limited and a hydrostatic variation is adopted. Clearly, as the value of permeability simultaneously depends on the mean effective stress and controls the pore water pressure profile, an iterative approach is required to establish the correct initial pore water pressure profile to be used in the analysis. However, it should be noted that the procedure followed was needed solely to ensure that the under-drained pore water pressure profile would represent the steady-state conditions and that the correct initial effective stress state was reproduced, as it controls not only the simulated undrained shear strength profile (equation (11)), but also the soil stiffness (equations (4) and (5)).

Since all the piles were wholly embedded in Cowden till (maximum length was 10·5 m, see Table 2), the soil profile adopted in all analyses is simplified to a single layer of this material, modelled in its entirety with the parameters listed in Table 1. The existence of the two deeper sand layers is ignored. Moreover, given the low permeability of Cowden till, the numerical analyses were performed as undrained, a condition achieved by specifying that the bulk modulus of the pore fluid was three orders of magnitude larger than that of the soil (see Potts & Zdravković (1999) for further details of this procedure).

As indicated by equation (11), the undrained shear strength simulated by the adopted constitutive model is a function of the stress state, the model parameters and the $OCR$⁠. Therefore, based on the measured variation with depth of the undrained shear strength in triaxial compression, $Su,TXC$ (Fig. 5(d); see Zdravković et al. (2019) for additional details), the $OCR$ profile is derived, as shown in Fig. 5(c). The fact that the adopted $OCR$ profile agrees very well with measured data is an indication of a high level of consistency between the calibrated model parameters and adopted ground conditions. This is further demonstrated in Fig. 5(e) by the very good agreement between the measured and predicted profiles of undrained strength in triaxial extension, $Su,TXE$⁠, in the top 8 m of the deposit. As $K0>1$ in this part of the ground, the initial soil state is in triaxial extension (TXE) and the predicted profile (grey line) is the result of the adopted deviatoric shape of the yield surface (Van Eekelen, 1980). The experimental data show a slightly larger scatter compared to the TXC tests in Fig. 5(d), but this is attributable to the usual issue of ‘necking’ failure, typical for triaxial testing in extension. The undrained strength profile at the start of 3D FE analysis is therefore a combination of the $Su,TXE$ profile in the top 8 m of the ground and $Su,TXC$ below 8 m.

The influence of the geometric characteristics of the monopile on its response is illustrated in Fig. 6, where the normalised deformed shapes of piles CM2 and CM3, characterised by the two extreme values of $L/D$⁠, are shown for a ground-level displacement, $vG$⁠, in the direction of loading (⁠$X$⁠) corresponding to 10% of their diameter (⁠$0⋅1D=76⋅2mm$⁠). Note that $vG$ corresponds to the displacement at the leading (front) edge of the pile. Throughout the PISA project, $vG=0⋅1D$ was assumed to correspond to the ultimate limit state to allow direct comparisons between the different piles. Although the two piles have the same diameter, they have different lengths, with values of $L/D$ of 3·0 and 10·0, respectively. Pile CM2 exhibits a rigid-body rotation about a point located at a depth corresponding approximately to 70% of the pile length, while the response of pile CM3 is significantly more flexible. Also shown in Fig. 6 are the maximum depths of the gaps formed around the two piles as a consequence of the applied loading. The depth of the gap varies considerably around the pile circumference, with its maximum occurring in the plane of symmetry, on the active side of the pile. The rigid-body rotation of pile CM2 results in the formation of a deeper gap in normalised terms, $LG/L=0⋅67$⁠, than is the case for pile CM3, for which $LG/L=0⋅53$⁠, $LG$ being depth of the gap. The significant extent of the gaps determined from the numerical model closely replicates the observed field behaviour (Byrne et al., 2019), with measured $LG/L$ of 0·72 and 0·48 for piles CM2 and CM3, respectively. This highlights the need to employ non-linear elasto-plastic interface elements, incorporating the possibility of gapping, in the simulation of laterally loaded piles, particularly for relatively low values of $L/D$⁠. In the analyses discussed here, it is assumed that the gap remained dry after formation and therefore the tension cut-off was set to zero. However, this cut-off value can be set as depth-dependent, corresponding to a possible hydrostatic pore water pressure distribution, in order to simulate water-filled gaps, such as the ones expected in piles installed offshore.

The importance of adopting the Van Eekelen (1980) deviatoric shape for the yield surface, which is capable of reproducing the effect of loading conditions on the strength of the material, is demonstrated in Fig. 7(a). This figure shows the spatial distribution of the value of the Lode's angle, $θ$⁠, obtained for pile CM2 at a ground-level displacement of $vG=0⋅1D$⁠. According to the adopted definition for $θ$ (see Appendix 1), $θ=−30deg$ corresponds to triaxial compression loading, whereas $θ=30deg$ characterises that in triaxial extension. As expected, the value of $K0>1$ adopted in the top 8 m (Fig. 5(a)) implies that the material surrounding the pile is initially characterised by $θ=30deg$⁠. As the pile is pushed forward in the $X$-direction, the loading conditions applied to the soil evolve rapidly and substantially. In effect, at the chosen stage of $vG=0⋅1D$⁠, the superficial zone in front of the pile, where relatively large plastic strains occur, assumes values close to $θ=−30deg$ in the immediate vicinity of the pile and therefore an operational strength of $ϕ′=27deg$⁠. The magnitudes of $θ$ increase with increasing radial distance in front of the pile, engaging again the initial $θ=30deg$ in the distant parts of the superficial soil domain that are unaffected by the pile loading. Furthermore, a considerable variation in the circumferential direction can be seen in Fig. 7(a), with the region immediately behind the pile being characterised by high values of $θ$⁠, for which the angle of shearing resistance is 32°, caused by the lateral unloading of the material which takes place until the gap is formed.

Contours of mobilised deviatoric strains, $Ed$⁠, are illustrated for the $vG=0⋅1D$ stage in Fig. 7(b). Note that, in this case, a logarithmic scale is adopted for the contour lines. Clearly, in the soil surrounding the pile, the strain level, and hence the shear stiffness of the material (equation (4)), changes rapidly, further reinforcing the need for a robust modelling approach which is capable of simulating soil behaviour over the entire range of deformation levels that occur in the problem.

Figure 8 compares the simulated embedded response of piles CM2 (⁠$L/D$ = 3) and CM3 (⁠$L/D$ = 10) with that interpreted from the field tests (Burd et al., 2019; Byrne et al., 2019). Note that, for a meaningful comparison, the plots in Fig. 8 are obtained by selecting stages of the analyses where either the ground-level displacement (Figs 8(a) and 8(c)) or the ground-level moment (Figs 8(b) and 8(d)) match those of the field test. Therefore, it is important to highlight that the deformed shapes do not necessarily correspond to the same loading level, while the bending moment distributions are not calculated for an identical value of ground-level displacement. Also shown in Fig. 8 are data, labelled ‘optimised structural model’, for the embedded pile response. This model was determined by obtaining a best fit between a structural model of the embedded pile and the measured strain gauge and inclinometer data. The process of determining this model is described in the paper by Burd et al. (2019).

The results in Fig. 8 demonstrate the excellent performance of the FE model, which is capable of reproducing with equal accuracy the measured deformations (Figs 8(a) and 8(c)) and bending moments (Figs 8(b) and 8(d)) for the two selected pile geometries, despite their significantly different values of slenderness ratio. More importantly, both the measurements and the model demonstrate a significant moment reaction at the toe of the rigid CM2 pile (Fig. 8(b)) and practically a zero value of this component at the toe of the flexible CM3 pile (Fig. 8(d)). The agreement between the numerical and field data confirms that a successful reproduction of a complex boundary value problem of this sort requires the accurate simulation of soil behaviour under the full range of strain levels, since ground movements appear to be concentrated in the top half of the pile, with the lower half imposing significantly smaller deformations on the soil.

The analyses results for the four chosen pile geometries in terms of load–displacement curves are compared in Fig. 9 with the observed field behaviour (Byrne et al., 2019). As the analyses were not designed to simulate strain-rate effects, the predicted pile capacities should be compared with the end points of experimental holding stages as these correspond to negligible strain rates. Clearly, up to ground-level displacements of $0⋅1D$⁠, there is a very good agreement between the FE simulations and the field results, with the loading capacity at $vG=0⋅1D$ being predicted with considerable accuracy for all pile geometries, thus suggesting that the predictive ability of the numerical model is practically unaffected by changes in pile slenderness ratio or diameter.

The importance of adopting a small-strain stiffness model within the elastic region of the modified Cam Clay model is further reinforced in Fig. 10, where the initial stages of the loading tests and associated predictions are shown in greater detail, by focusing on ground-level displacements of up to $vG=0⋅01D$⁠. In addition to reproducing the behaviour at large displacements where the stiffness at moderate to large strains is dominant (Fig. 9), the FE model successfully captures the initial stiffness of the pile–soil system for all geometries. This excellent agreement between the FE model and the measured response demonstrates the suitability of the chosen modelling approach for the shear modulus at very small to moderate strains, which in turn provides a systematic connection with the results of the field and laboratory testing (Fig. 3).

To assess the quality with which the observed load–displacement response is reproduced by the FE model, the following accuracy metric is adopted

where $Aref$ is the area below the reference load–displacement curve (in the present case, that obtained in the field tests) and $Adiff$ is the area of the region delimited by the reference and simulated curves. The procedure followed for determining the reference curves is explained in the paper by Byrne et al. (2019).

The geometric representations of $Aref$ and $Adiff$ are shown in Fig. 11 for two distinct cases: in the first (Fig. 11(a)), which is designated $η0⋅1D$⁠, the overall performance of the developed numerical model is assessed by considering ground-level displacements up to 10% of the pile diameter, while in the second (Fig. 11(b)), the range of displacements is limited to 0·1% of the pile diameter, to evaluate the ability of the numerical model to reproduce the observed initial stiffness of the pile. The latter metric is labelled $ηsd$⁠, as it focuses on the small-displacement range and its introduction is justified by the need to assess the adequacy of the numerical model to predict the response of the soil–pile system under operational conditions, which tend to be characterised by a large number of small loading cycles. A value of $η=1$ corresponds to a perfect reproduction of the reference curve, with any deviations from that behaviour leading to a reduction of its value (⁠$Adiff$ is always positive and therefore $η$ assumes values always equal to or below $η=1$⁠.

The interpreted values of $η0⋅1D$ and $ηsd$ are illustrated in Figs 12(a) and 12(b), respectively, for the four considered pile geometries. The results of the FE analyses are seen to replicate with reasonably high accuracy both the overall observed pile response (average $η0⋅1D>88%$⁠) and the behaviour at small displacements (average $ηsd>85%$⁠). It is also noted that the values of $η0⋅1D$ are consistently high for all considered geometries, which include two distinct values of diameter and three values of the pile slenderness ratio, $L/D$⁠, ranging from 3 to 10.

This paper presents a systematic and comprehensive procedure for developing a numerical model for the simulation of laterally loaded piles in overconsolidated glacial clays, which is applied in the study of laterally loaded piles installed at the PISA test site at Cowden. The 3D FE analyses were conducted before the pile tests were performed and the following were the key steps of the numerical procedure.

• (a)

  The choice of the soil constitutive model and its calibration for Cowden till were based on the project objectives (monotonic loading, no cyclic or strain-rate considerations) and on available experimental evidence for the foundation soil. Particular emphasis was placed on the combined interpretation of field and laboratory testing results, as well as on the significant engineering judgement needed to ensure the consistency of calibration and realistic soil response under a wide range of strains. Based on experimental data, the adopted modified Cam Clay framework (Roscoe & Burland, 1968) was extended crucially with a non-linear Hvorslev surface (Tsiampousi et al., 2013), a generalised shape of the yield and plastic potential surfaces in the deviatoric plane (Van Eekelen, 1980), and a non-linear variation of the elastic shear modulus, dependent on both stress and strain levels (Taborda et al., 2016).

• (b)

  The use of advanced features such as limited-tension interface elements, was shown to be critical for the precise prediction of the observed field behaviour in terms of load–displacement curves, pile deflections and bending moments.

• (c)

  The FE model developed in ICFEP (Potts & Zdravković, 1999) was validated by predicting accurately the response to lateral loading of four distinct PISA test piles installed at Cowden. Their geometries include two different diameters (0·76 m and 2·0 m) and three values of slenderness ratio, $L/D$ (3·0, 5·25 and 10·0).

The close agreement of the predicted and measured responses of the four piles proved the premise of the PISA project that the new Winkler-type design method for laterally loaded piles could be derived directly from advanced, site-specific, FE modelling. As a result, the application of this FE model was key in the development of the new PISA monopile design methodology, which is outlined in the paper by Byrne et al. (2017).

The PISA project was funded by the UK Department for Energy and Climate Change (DECC) and the PISA industry partners under the umbrella of the offshore wind accelerator (OWA) programme, which was designed and is led by the Carbon Trust. The authors acknowledge the provision of financial and technical support by the following project partners: Ørsted Wind Power (formerly DONG Energy), Alstom Wind, E.ON, EDF, Equinor (formerly Statoil), innogy, SPR, SSE, Vattenfall and Van Oord. The authors acknowledge gratefully the work of Socotec UK Ltd (formerly ESG) as the main contractor for the design and execution of the field testing programme.

Given the effective stress tensor, $σ′$⁠, defined as

and a compression-positive sign convention, the invariants p′ (mean effective stress), $J$ (generalised deviatoric stress) and $θ$ (Lode's angle) are calculated using

where

Similarly, for the strain tensor, $ε$

the invariants $εvol$ (volumetric strain) and $Ed$ (generalised deviatoric strain) are calculated using

The calibration process of the adopted constitutive model was concluded by assessing parameter $β$⁠, which controls the shape of the plastic potential on the dry side. This was carried out using a trial-and-error procedure, where the simulated response for various values of $β$ was compared to that measured experimentally. Although this was carried out for multiple tests, in Fig. 13 only the results obtained for test CR100IUC8.2 – an undrained triaxial compression test on a 100 mm sample collected at 8·2 m depth (see Table 3 and Zdravković et al. (2019) for test details) – are shown in terms of the axial strain–deviatoric stress (⁠$εa$–q) and axial strain–excess pore pressure (⁠$εa$–$Δu$⁠) diagrams. Based on the observed behaviour, a value of $β$ = 0·20 was selected.

To demonstrate the model's performance using the final set of model parameters, listed in Table 1, two other triaxial tests carried out on 100 mm samples, collected at two depths – 2·50 m and 5·00 m – were simulated. The initial conditions of the tests are listed in Table 3 (also in the paper by Zdravković et al. (2019)) and the results of the simulations, in terms of $εa$–q and $εa$–$Δu$ plots, are compared with the test data in Fig. 14. Clearly, there is a good agreement between observed and simulated responses, confirming the suitability of the chosen modelling approach and respective calibration.

$Adiff$

area difference between reference and simulated curve in equation (13)

$Aref$

area below reference curve in equation (13)

$a,b,mG,RG,min$

fitting parameters for shear modulus in equation (4)

$D$

pile diameter

$detssfs$

determinant of stress tensor

$fG(e)$

void ratio function for shear modulus in equation (4)

$fK(e)$

void ratio function for bulk modulus in equation (5)

$G0$

elastic shear modulus

$G0∗$

basic shear modulus in equation (4)

$Gsec$

secant shear modulus

$Gtan$

tangent shear modulus

$gθ$

inclination of the critical state line

$H$

applied lateral load

$h$

height of point of application of lateral load above ground level

$J$

generalised deviatoric stress

$K0$

earth pressure coefficient at rest

$K0∗$

basic bulk modulus in equation (5)

$Ktan$

tangent bulk modulus

$k$

soil permeability

$L$

pile embedded length

$LG$

depth of the gap

$M$

bending moment in pile

$m$

fitting parameter for the shape of plastic potential on dry side

$n$

fitting parameter for the shape of yield surface on dry side

$p′$

mean effective stress

$p′0$

hardening parameter for the modified Cam Clay model

$p′c$

current mean effective stress in equation (2)

$p′ref$

reference pressure

$r,s,mK,RK,min$

fitting parameters for bulk stiffness in equation (5)

$Su$

undrained shear strength

$Su,TXC$

undrained shear strength in triaxial compression

$Su,TXE$

undrained shear strength in triaxial extension

$t$

pile wall thickness

$v$

specific volume

v

horizontal displacement of pile

$vG$

horizontal displacement of pile at ground level

$X,Y,Z$

parameters of the yield surface in deviatoric plane in equation (3)

$X,Y,Z$

finite-element mesh coordinate system

z

depth below ground surface

$α$

fitting parameter for the shape of yield surface on dry side

$β$

fitting parameter for the shape of plastic potential on dry side

$γxy,γxz,γyz$

shear strain components of the strain tensor

$Ed$

generalised deviatoric strain

$ε$

strain tensor

$εvol$

volumetric strain

$εx,εy,εz$

direct strain components of the strain tensor

$η$

accuracy parameter in equation (13)

$θ$

Lode's angle

$θc$

current value of Lode's angle in equation (2)

$κ$

inclination of isotropic swelling line in ν–$ln⁡p′$ plane

$λ$

inclination of isotropic normal compression line in ν–$ln⁡p′$ plane

$MJ,TXC$

generalised stress ratio in compression (⁠$J/p′$⁠)

$MJ,TXE$

generalised stress ratio in extension (⁠$J/p′$⁠)

$MJ,θ$

generalised stress ratio as function of Lode's angle, $θ$⁠, in equation (10)

$ν$

Poisson's ratio

$σ′$

effective stress tensor

$σ′x,σ′y,σ′z$

direct stress components of the effective stress tensor

$τxy,τxz,τyz$

shear stress components of the effective stress tensor

$ϕ′TXC$

angle of shearing resistance in triaxial compression

$ϕ′TXE$

angle of shearing resistance in triaxial extension

$ϕ′θ$

angle of shearing resistance as function of Lode's angle, $θ$⁠, in equation (10)