Finite-element modelling of laterally loaded piles in a dense marine sand at Dunkirk Open Access

The PISA (pile–soil analysis) project combined ground characterisation, field testing and computational analysis to develop new design models for large-diameter monopiles as offshore wind turbine foundations. Outlined in the paper by Zdravković et al. (2019a) and denoted as the PISA design model, it retains the simplicity of the existing one-dimensional (1D) Winkler-type p–y approach but includes additional soil reactions. The objective of the PISA project was to develop new design models from the results of validated site-specific three-dimensional (3D) finite-element (FE) analyses. Consequently, as set out in the paper by Zdravković et al. (2019a), the project first had to demonstrate the ability of 3D FE analyses to predict with sufficient accuracy the response of test piles installed, as part of the project, in a dense marine sand (McAdam et al., 2019) and in a stiff glacial clay (Byrne et al., 2019) at the Dunkirk and Cowden test sites, respectively. Upon completion of this step, 3D FE analyses of full-scale monopiles would be performed, from which 1D simplified Winkler models for piles in clay and sand could be derived. The current paper presents the development of a 3D FE model for the test piles installed in sand and compares the numerical predictions of pile behaviour, obtained before the testing under lateral loading took place, with field measurements. As a result, the predicted behaviour of the test piles is established solely from the interpretation of ground conditions and available experimental evidence for soil behaviour (i.e. without any back-analysis of the pile tests). A companion study for test piles installed in a glacial clay till at Cowden is reported in the paper by Zdravković et al. (2019b).

The characterisation and modelling of sand behaviour is a complex task, with this type of material exhibiting a marked effect of the confining pressure and void ratio on its dilatancy and peak strength (Been & Jefferies, 1985; Bolton, 1986). It is therefore unsurprising that centrifuge testing on laterally loaded monopiles embedded in sand has shown a strong dependency between the monopile ultimate load capacity and the relative density, D_R, of the material (e.g. Rosquoet et al. (2007) report an increase in load capacity of about 30% when the D_R of Fountainebleau sand increases from 53% to 86%). As a result, when modelling the response of laterally loaded piles in sand deposits, relatively simple constitutive models, which predict constant strength and dilatancy properties (e.g. Mohr–Coulomb type models), require the adoption of different sets of parameters based on the relative density of the material (e.g. Achmus et al., 2009). Consistent results have been obtained following this approach, as shown by Abdel-Rahman & Achmus (2005), Achmus et al. (2009) or Stone et al. (2018).

However, for the PISA project the need to establish an implicit link between the state of the sand and the behaviour simulated by the constitutive model is clear: current monopiles are becoming larger in diameter (up to 10 m) and depth (length-to-diameter ratios L/D up to 6), meaning that the range of confining pressures applied to the sand around the pile is becoming sufficiently large to result in significant variations in the sand's mechanical response. Moreover, given the large areas occupied by offshore wind farms, it is unlikely that a limited number of relative densities characterise all the sand deposits encountered. Therefore, in order to account for site-specific conditions, it is necessary that the constitutive model is capable of reproducing the behaviour of sand under a wide range of relative densities using a single set of parameters. Failure to do so would require the determination of a parameter set for each of the values of relative density, a time-consuming and expensive exercise. Ahmed & Hawlader (2016) carried out analysis of monopiles in sand using a version of the Mohr–Coulomb failure criterion where the strength and dilatancy properties are directly related to the state of the material using the relative dilatancy index, I_R, described by Bolton (1986). Alternatively, constitutive models have been proposed based on the state parameter, ψ, introduced by Been & Jefferies (1985) as a form of predicting the influence of mean effective stress and void ratio on the behaviour of sand within the framework of critical state soil mechanics. These are often based on the bounding surface plasticity model proposed by Manzari & Dafalias (1997).

In the current paper, the application of the latter type of constitutive model to the 3D FE modelling of the laterally loaded PISA test piles installed in Dunkirk sand (McAdam et al., 2019; Zdravković et al., 2019a) is presented, providing a detailed and extensive demonstration of the challenges associated with the incorporation of laboratory and field information into a global 3D FE model of such piles. The resulting level of agreement between field measurements and the 3D FE predictions of the test piles’ response was pivotal to the development of the new 1D PISA design model for laterally loaded piles in a dense sand. It is emphasised that the PISA project, and hence the adopted numerical modelling, was focused on monotonic loading only (see Zdravković et al. (2019a)), as an essential starting point for future extensions of the new 1D PISA design model to cyclic and other loading conditions.

The constitutive model chosen to represent the behaviour of Dunkirk sand is the state parameter-based bounding surface plasticity model described in the paper by Taborda et al. (2014), which is an evolution of the model originally proposed by Manzari & Dafalias (1997). This choice follows the interpretation of the Dunkirk ground model (Zdravković et al., 2019a), which demonstrated the applicability of the critical state framework to describe the mechanical behaviour of a dense, quartzitic, marine sand. A further important requirement was that of the model being capable of reproducing the effects of void ratio, e, and mean effective stress, p′, on the sand behaviour, with a single set of calibrated model parameters.

This section summarises the aspects of the adopted constitutive model which are directly related to its calibration for Dunkirk sand and to its application to the analysis of the piles tested at this site. Taborda (2011) and Taborda et al. (2014) provide further details on the model formulation, implementation and calibration. In general, the model utilises four distinct surfaces to describe the mechanical response of sands: the yield surface, as the boundary of the elastic region; the critical state surface, defining the stress conditions at failure; the dilatancy surface, defining the transition from plastic contraction to plastic dilation; and the bounding surface, which controls the sand's peak strength and the plastic strain rate.

The yield surface is a narrow cone in stress space, of a constant opening, as the model utilises only kinematic hardening. The other three surfaces have the shape of an open wedge in stress space, with their deviatoric plane shapes depicted in Fig. 1. The opening of the critical state surface is defined by the stress ratios (q/p′) corresponding to the strength of the sand in triaxial compression and triaxial extension, denoted as M^c_c and M^c_e, respectively. The openings of the bounding and dilatancy surfaces are related to that of the critical state surface through the state parameter, ψ (Been & Jefferies, 1985), defined as the difference between the current void ratio, e, and that at the critical state line (e_CS) for the same value of the mean effective stress, p′

In the equation above, e_CS,ref is the void ratio at critical state for p′ = 0, p′_ref is a reference pressure and λ and ξ are parameters defining the shape of the critical state line (CSL) in the e–lnp′ plane. This CSL shape follows the commonly adopted power law expression of Li & Wang (1998), which was used in the paper by Zdravković et al. (2019a) for the interpretation of laboratory test data. With this, building upon the earlier work by Wood et al. (1994), Manzari & Dafalias (1997) proposed that the positions in triaxial compression of the bounding (M^b_c) and dilatancy (M^d_c) surfaces are determined as

where k^b_c and k^d_c are model parameters and 〈□〉 denotes Macauley brackets, according to which 〈x〉 = x if x > 0 and 〈x〉 = 0 if x < 0. For simplicity, the ratios between the openings in triaxial compression and triaxial extension for the bounding (M^b_e) and dilatancy (M^d_e) surfaces are assumed to be identical to that of the critical state surface, resulting in

For very dense sands, for which ψ is generally a large negative number, the adopted formulation (equation (2a)) implies that a very high stress ratio characterises the bounding surface, therefore simulating a high peak strength. On the other hand, the stress ratio defining the dilatancy surface reduces, implying the onset of plastic dilation early during shearing. Other components of the model formulation that control the plastic response, such as the plastic potential, the mapping rule, the hardening modulus and the possibility to introduce a fabric tensor to represent the evolution of the structure of the sand during cyclic loading, are described in detail in the paper by Taborda et al. (2014).

In terms of elastic behaviour, the model adopts a non-linear elastic Ramberg & Osgood (1943) type stiffness degradation, designed to reproduce an hysteretic soil response, and hence energy dissipation, within the elastic region

where G_tan is the tangent shear modulus; κ and a₁ are model parameters; χ^r_ref is a measure of deviatoric stress with respect to the last reversal point; and η₁ is a normalised stress ratio. Moreover, N is a scaling factor which, in accordance with the Masing rules (Masing, 1926; Kramer, 1996), initialises as 1·0 and changes to 2·0 when the first reversal in shear direction is detected. Lastly, G₀ is the elastic shear modulus, given by the expression of Hardin (1978) 

in which B is a model parameter.

The bulk modulus is determined from the theory of elasticity using the following expression, in which ν is the Poisson's ratio

Calibration of the adopted bounding surface plasticity model was carried out following the hierarchical approach proposed for this type of model, as outlined in the papers by Loukidis & Salgado (2009) and Taborda (2011), applying in the process sound engineering judgement to integrate (often conflicting) laboratory and field evidence of the Dunkirk sand's behaviour. A brief overview of the calibration is given here, focusing on the parameters defining the strength, critical state and stiffness of the Dunkirk sand. The experimental data available for model calibration were summarised and discussed in the paper by Zdravković et al. (2019a). This comprised a set of drained triaxial tests sheared from K₀ stress conditions and with a single void ratio value (Aghakouchak, 2015), and a new PISA-specific set of drained triaxial tests sheared from isotropic conditions (Liu et al., 2017), involving a wider range of void ratios, overconsolidation ratio (OCR) values and initial stresses. Both historic (Chow, 1997) and new field measurements of the elastic shear modulus, G₀, were also available.

While recognising the usual difficulty in establishing the ultimate strength from sand samples sheared under triaxial conditions, the characterisation of the Dunkirk sand deposit in the paper by Zdravković et al. (2019a) was able to utilise the available triaxial test data to establish the critical state strength parameters for this model in triaxial compression and extension, as M^c_c = 1·28 and M^c_e = 0·92, respectively. These ratios correspond to angles of shearing resistance ϕ′_TXC = 32° and ϕ′_TXE = 33°, with such similarity being common to other pluviated sands, as discussed in the paper by Zdravković et al. (2019a). The adopted values lead to a ratio $c=Mec/Mcc=0⋅72$⁠, which is the limit value that ensures the convexity of the model surfaces in the deviatoric plane (see Fig. 1), a condition necessary for guaranteeing numerical stability when integrating the constitutive model (Loukidis & Salgado, 2009).

The position of the CSL is assessed by plotting the triaxial compression data in the e–lnp′ plane, as shown in Fig. 2 and discussed in the paper by Zdravković et al. (2019a). The model adopts the non-linear expression of Li & Wang (1998) to define the void ratio at critical state, e_CS, as given in equation (1). Two possible positions of the CSL were derived from these data, emphasising again the complexity of interpreting sand behaviour. Although the shearing of looser samples, in this case from the pre-shearing void ratio, e₀ ≈ 0·74 (corresponding to a relative density, D_R ≈ 45%), is expected to lead to a more accurate estimate of the CSL due to a smaller volumetric dilation, the adopted location of the CSL is that corresponding to drained shearing from e₀ ≈ 0·64, as this corresponds to D_R ≈ 75%, estimated as the value of relative density of the natural Dunkirk sand (Zdravković et al., 2019a). To define the adopted CSL shape in Fig. 2, the model parameter p′_ref in equation (1) is set to the atmospheric pressure of 101·3 kPa, while the e_CS,ref is set to the maximum void ratio of Dunkirk sand, e_max = 0·91, as proposed by Riemer et al. (1990) and discussed in the paper by Zdravković et al. (2019a). Applying a non-linear regression, the remaining two parameters, λ and ξ in equation (1), are calculated as 0·135 and 0·179, respectively.

According to equation (2a), the opening of the bounding surface, which provides an indication of the stress ratio at peak angle of shearing resistance and is expressed as M^b_c, is a linear function of the state parameter. Consequently, the model parameter k^b_c can be evaluated by plotting, as shown in Fig. 3, the values of $(Mcb−Mcc)$ against the corresponding values of ψ_peak = e_peak − e_CS, using the adopted CSL from Fig. 2 and assuming that for ψ_peak = 0 the value of M^b_c must coincide with the critical state stress ratio M^c_c ( = 1·28). Following this procedure, a preliminary estimate of k^b_c = 3·30 is obtained (see also Zdravković et al. (2019a)), being in the range from 0·5 to 4·0, as suggested by Papadimitriou & Bouckovalas (2002). It should be highlighted that this value is only an estimate as it largely depends on the assessed magnitudes of ψ_peak, the accuracy of which is inherently controlled by the experimentally measured pre-peak volumetric changes.

A similar linear relationship is implied in equation (2b) between the stress ratio at phase transformation (PT), M^d_c, and the corresponding state parameter. As only drained tests were available for model calibration, the PT state is characterised by the transition from plastic contraction to plastic dilation, as discussed in the paper by Zdravković et al. (2019a). To derive the model parameter k^d_c, which controls the opening of the dilatancy surface, the values of (⁠$Mcd−Mcc$⁠) are plotted against the corresponding values of ψ_PT = e_PT − e_CS, similar to the procedure shown in Fig. 3 for parameter k^b_c. However, it should be noted that the identification of parameter k^d_c has a higher degree of uncertainty, being affected by the elastic shear and bulk moduli employed during the integration procedure required to estimate the elastic strains in order to calculate the observed plastic strains. A possible way of increasing the confidence in the value of k^d_c would be to perform undrained triaxial compression tests where the transition from plastic contraction to plastic dilation corresponds, in the case of isotropic materials as it is assumed in this constitutive model, to the minimum value of mean effective stress (i.e. dp′ = 0). Using the adopted CSL in Fig. 2, the elastic parameters detailed in the following section and ensuring that the linear regression satisfies the condition $Mcd=Mcc(=1⋅28)$ when ψ_PT = 0, an estimate of k^d_c = 0·88 is obtained. This value is within the range proposed by Papadimitriou & Bouckovalas (2002) (i.e. from 0·1 to 3·0, see also Zdravković et al. (2019a)).

The interpretation of the elastic shear stiffness of Dunkirk sand, from both laboratory and field experiments, is discussed in the paper by Zdravković et al. (2019a). The available drained triaxial tests of Aghakouchak (2015) and new tests of Liu et al. (2017) were equipped only with local strain instruments (i.e. no bender elements), which were analysed first to establish the vertical Young's modulus profile, E_v, that was subsequently converted to the elastic shear modulus, G₀, by adopting a Poisson's ratio, ν = 0·17, as estimated by Kuwano (1999). The discussion in the paper by Zdravković et al. (2019a) demonstrated that the profile of G₀ interpreted in this way is independent of the OCR and that it could be well represented by the classic expressions of Hardin & Black (1968) or Hardin (1978). The latter expression is adopted by the constitutive model (equation (5)) and the parameter B is calibrated by plotting, in Fig. 4, the triaxial data as G₀ against the modified mean effective stress, p*

As expected, the two sets of tests yield distinct values of B due to the different consolidation states prior to shearing (B = 940 for K₀ and B = 620 for isotropic consolidation). Also depicted in this graph is the behaviour corresponding to B = 875, estimated from the interpreted G₀ profiles from the new seismic cone penetration tests (SCPTs, Fig. 5) in the top 10 m of the deposit, this being the maximum depth of the test piles. The deeper data appear anomalous, indicating a reduction in G₀. The resulting G₀ profiles fitted to the individual three sets of data, using equation (5), are depicted in Fig. 5.

The non-linear degradation of the elastic stiffness in equation (4) is controlled by parameters a₁, γ₁ and κ. The latter is assumed to have a value of 2·0, as suggested by Papadimitriou & Bouckovalas (2002), with 0·40 prescribed for a₁.

The plastic behaviour of the model is governed by the flow rule and the hardening modulus. The former is characterised by a single parameter, A₀, which can be derived from a linear relationship between the stress ratio q/p′ and the plastic dilatancy rate, $Dep=δεvolplas/δεsplas$⁠, as explained in the papers by Loukidis & Salgado (2009) and Taborda (2011), giving a value of A₀ = 1·3.

The last set of parameters to be calibrated are those used in the calculation of the hardening modulus

where only h₀ is a model parameter which controls the overall magnitude of A and, hence, the plastic strain rate. The remaining components, respectively, denote the effects of the void ratio (h_e), non-linear elastic shear modulus (h_g), distance of the stress state to the bounding surface (h_b) and generated soil fabric (h_f), with d^b being the current distance to the bounding surface. Details of each of these components are given in the paper by Taborda et al. (2014). Given the limited available data, it is assumed here that h_e = 1·0 (i.e. the model parameter γ is set to 0·0) and that h_f = 1·0 (i.e. model fabric parameters H₀ = 0·0 and ζ = 0·0). The current value of the nonlinear elastic tangent shear modulus is introduced in the calculation of A by setting the model parameter α = 1·0, as suggested in the paper by Taborda et al. (2014). Parameters μ and β of the component related to the distance and opening of the bounding surface, h_b, are set to their recommended values of 1·0 and 0·0, respectively (Taborda et al., 2014).

The calculation of h₀ requires a trial-and-error procedure and, unlike in situations where a single source of soil data is used, is performed here in two distinct steps: the first is to determine the value required for the best approximation to the available laboratory triaxial tests, while the second step is designed to provide a consistent way of taking into account, in a calibration based on the behaviour measured in triaxial tests, a considerably higher shear stiffness measured in the field. The first step, commonly carried out in the calibration of models of this type, consists of assessing the value of h₀ required for the best reproduction of each test, and then adopting a representative value, usually its average magnitude, for the entire set of tests. In the present case, the first step involved the use of an elastic stiffness parameter, B = 620, as it corresponds to the best approximation of stiffness based on the triaxial data by Liu et al. (2017), yielding a value of h₀ of 0·023. A second estimate of h₀ of 0·045 was then obtained assuming that the elastic behaviour corresponded to that measured in the field data (i.e. B = 875). However, it should be noted that, in the latter case, where a compromise between the laboratory and field data is attempted, the impact of the higher shear stiffness measured in the field on the model calibration is confined to the maximum elastic stiffness, G₀ (equation (5)). In effect, as pointed out by Tatsuoka & Shibuya (1991) and discussed by Pedro et al. (2017), the effect of sample disturbance, which in this case would be significant as the tested samples were reconstituted to the best estimate value of in situ void ratio and would have lost all in situ fabric, should reflect itself at all strain levels, although to different degrees, and not just at very small deformations where G₀ would be representative of soil behaviour. In order to address this issue, parameter h₀ was increased until a smooth variation with deformation level of the shear stiffness was achieved, without resulting in gains in stiffness larger than the ratio $G0field/G0lab$⁠. A value of h₀ of 0·4 was found to adequately meet this criterion, with this increase having to be accompanied by a reduction in k^b_c from 3·30 to 2·70 (see Fig. 3), in order to prevent the over-prediction of the material's peak strength (see Appendix). The last parameter to be established is p′_YS, which controls the position of the secondary yield surface, a model component introduced to increase the stability of the model by preventing excessively low values of mean effective stress from being calculated (see Taborda et al. (2014) for additional details on this aspect of the model). In the present case, a value of 1·0 kPa was assumed, following the recommendations of Taborda (2011).

The calibration procedure outlined above followed the necessary steps to determine two sets of model parameters: one which reproduces best the measured laboratory data and a second one where the higher stiffness observed in field tests is included. However, the application of any constitutive model in a boundary value problem, where both initial stress and loading conditions may be radically different from those imposed in the tests used for calibration, must be accompanied by additional verifications. In the case of the chosen constitutive model, as the simulated behaviour is intrinsically related to the state parameter (equation (1)), it is fundamental to assess the implications of the initial stress conditions on the modelled response. In Fig. 6, the profile of initial state parameter, established from the ground characterisation outlined in the paper by Zdravković et al. (2019a), is compared to the values of state parameter characterising the laboratory tests performed by Liu et al. (2017). Clearly, with the exception of the two tests carried out with an initial void ratio of 0·58, none of the triaxial tests appears to replicate the in situ state parameter for the top 10 m of the deposit, the length of the longest installed pile. Moreover, no test reaches values of state parameter comparable to those estimated for the hydraulic fill layer in the top 3 m of the deposit. Such a mismatch between the tested samples and the initial ground conditions is not surprising, considering the very high relative density of the in situ material and the difficulties in performing triaxial tests under low mean effective stresses, but requires special treatment when establishing a numerical model.

The impact of the very large (in absolute terms) values of state parameter estimated for the top 5–10 m is most pronounced on the modelled peak strength. It is worth adding a reminder that this type of model does not allow for a peak strength to be prescribed directly since, unlike other bounding surface plasticity models (e.g. Grammatikopoulou et al., 2006), the adopted formulation allows for the bounding surface to be crossed, this being the mechanism that triggers strain softening (d^b < 0 in the hardening modulus, equation (8)). Therefore, it is important to highlight that the opening of the bounding surface may only be used as a proxy when controlling the simulated peak strength. For the adopted parameters, the opening of the bounding surface for various values of the state parameter is shown in Fig. 7, measured both in terms of the mobilised angle of shearing resistance, ϕ_mob, and normalised by the value under triaxial compression loading, ϕ_mob/ϕ_mob^TXC. It can be seen clearly that, for the magnitudes of state parameter estimated in the field, large values (up to about 55°) of mobilised shear resistance in triaxial compression conditions can be potentially achieved by the model. While this is naturally a direct consequence of the state parameter framework and not of the specific model formulation, the same cannot be said for loading conditions other than triaxial compression, where the simulated peak strength is mostly a function of the adopted deviatoric shape for the model, which in the present case is

where, as previously introduced, $c=Mec/Mcc$ and θ is the Lode's angle. For the values of state parameter estimated for the field conditions, values of peak strength of up to 80° are possible due to the change in shape of the bounding surface for non-triaxial compression loading conditions, which is quantified in Fig. 7(b) by the ratio ϕ_mob/ϕ_mob^TXC. Such observations highlight potential pitfalls of extrapolating ground behaviour from the knowledge gathered through a number of laboratory tests in which the material is subjected to limited loading conditions. To address this potential issue, a limit is introduced to the opening of the bounding surface, corresponding to ϕ_mob = 50° under plane-strain conditions (approximately Lode's angle of 0°), which corresponds to a maximum stress ratio under triaxial compression loading conditions of M^b_c,max = 1·631, equivalent to ϕ^peak_TXC,mob = 40°. This value is only slightly below the maximum value of about 42°, estimated from the laboratory triaxial tests of Liu et al. (2017).

Table 1 summarises the model parameters derived considering only the triaxial test data of Liu et al. (2017) and the adjusted set of parameters obtained from the above second step of calibration, which additionally incorporates the higher elastic stiffness measured in the field. The performance of the model with the two sets of parameters is briefly illustrated in the Appendix, demonstrating that model adjustments have ensured the reproduction of the same peak strengths and volumetric behaviour as the laboratory-based model calibration.

In total, 14 open-ended steel piles were installed at the Dunkirk site and subjected to lateral loading, with their geometric characteristics and loading steps summarised in the paper by McAdam et al. (2019). The smallest-diameter piles (D = 0·273 m) were not considered in the numerical analyses as they were principally used to check the bespoke loading and monitoring system on site. The medium- (D = 0·762 m) and large- (D = 2·0 m) diameter piles comprised four different geometries, summarised in Table 2, which were subjected to detailed numerical analyses and the predictions compared with field measurements. Multiple piles with the same geometries were tested to demonstrate the repeatability of the field measurements. All analyses presented here are performed with the Imperial College Finite Element Program (ICFEP) (Potts & Zdravković, 1999).

The geometry contains a plane of symmetry and only half of the problem is discretised in the finite-element mesh, as shown in Fig. 8. The soil domain is represented with 6464 20-noded displacement-based hexahedral elements, with its bottom boundary at 40 m depth. The vertical cylindrical boundary is at a radial distance, R, of 30 m for the medium-diameter piles (D = 0·762 m, R/D ≈ 40) and at R of 40 m for the large-diameter (D = 2·0 m, R/D = 20) piles. The pile is discretised with 280 eight-noded shell elements (Schroeder et al., 2007), with the behaviour of steel assumed linear-elastic and represented with a Young's modulus, E = 200 GPa and a Poisson's ratio, ν = 0·30. Special interface elements are introduced around the outside of the pile to allow appropriate constitutive modelling of the pile–soil interface. These are 16-noded zero-thickness elements (Day & Potts, 1994), which allow for the separation between the pile and soil to occur, if triggered, and are simulated with an elasto-plastic Mohr–Coulomb model. The elastic parameters of this model are the normal (K_N) and shear (K_S) stiffness, which are both set to 10⁵ kN/m³, ensuring a minimal compression deformation across the element and a high shear stiffness prior to the onset of plastic behaviour, respectively. The plastic part is characterised with zero cohesion (c′ = 0) and an angle of shearing resistance equal to that of the soil at failure (i.e. critical state) in triaxial compression (ϕ′ = 32°).

The applied boundary conditions prevent movements in all coordinate directions (X, Y, Z) at the base of the mesh (at Z = −40 m) and in the direction normal to the vertical cylindrical boundary. The plane of symmetry at Y = 0 requires that the displacements in the Y-direction must be set to zero over this plane. Equally, the rotational degrees of freedom about the X- and Z-axes must be set to zero at the nodes of shell elements contained in the Y = 0 plane. The loading of the pile (at Z = h = +10 m) is simulated by applying uniform increments of horizontal displacement in the X-direction around the perimeter of the pile. The resulting reactions at all of these nodes are summed to give the total load, H/2.

Overall characterisation of the ground conditions at the Dunkirk site was discussed in the paper by Zdravković et al. (2019a). Consequently, only information relevant for the numerical input is presented here. This requires initialisation of stresses in the ground, using the saturated bulk unit weight γ_sat = 17·1 kN/m³ above the water table and 19·9 kN/m³ below it, as reported by Chow (1997), to calculate the total vertical stress. The pore water pressure profile on site is estimated from the new cone penetration test using piezocone (CPTu) data shown in Fig. 9, indicating the groundwater table at 5·4 m depth and a hydrostatic distribution below it. Above the water table the conditions are uncertain, as piezocones are not suited for measuring large negative pore water pressures. However, from the work on site it was evident that the ground surface was able to sustain significant effective stresses. As discussed in the paper by Zdravković et al. (2019a), no high-quality sampling was conducted at Dunkirk, in particular for the fill material in the first 3·0 m, and therefore the source of strength and stability of the fill is unclear. Consequently, the indication of some suction above the groundwater table in Fig. 9 was considered a reasonable option for generating additional strength in the fill in the numerical analyses. It was not expected that full hydrostatic suction would develop in the sand over the 5·4 m height from the groundwater table. Therefore, the profile of suction depicted in Fig. 9 is adopted: hydrostatic suction from the groundwater table to 4·0 m depth, constant suction of 13·5 kPa to 2·0 m depth, then reducing to 10 kPa at the ground surface. The horizontal effective stress is then initiated by applying the coefficient of earth pressure at rest, K₀ = 0·4, as recommended by Chow (1997) for a normally consolidated Dunkirk sand.

For the chosen constitutive model it is also necessary to initialise the void ratio, e₀, in the ground. Following the interpretation of the initial density profile in the paper by Zdravković et al. (2019a), e₀ = 0·54 is prescribed in the 3·0 m thick hydraulic fill, estimated to be at 100% relative density, while e₀ = 0·628 is specified for the natural Dunkirk sand at 75% relative density.

The same initial pore water pressure profile is prescribed in the interface elements located at the pile–soil interface, with the interface effective normal stress prescribed as equal to the horizontal effective stress in the adjacent soil elements. This ensures the continuity and equilibrium of the stress field at the start of analysis. The pile is modelled as wished in place, hence installation effects are not considered as they would need to be accurately quantified for a reasonable inclusion in the 3D FE model. Otherwise, this brings additional uncertainty in the numerical model. This modelling decision is further supported by the knowledge that, under lateral loading, a large volume of soil is mobilised around the pile, exceeding the boundary of the interface zone that may be disturbed by pile installation. The implication is that the disturbance of the pile–soil interface is considered less significant when determining pile capacity for this type of loading, compared to an axially loaded pile whose capacity strongly depends on the interface conditions. All analyses are performed under drained conditions, according to the envisaged field loading rates, to establish monotonic backbone load–displacement curves.

The normalised deflected shapes of all four test piles, shown in Fig. 10, indicate a predominant rigid-body rotation mode of deformation, even for the longest medium-diameter pile (DM3, L/D = 8). The graphs are shown for the ground-level horizontal displacements of 0·01D and 0·1D, the latter representing the adopted ultimate condition of the piles. Although some flexibility is retained in the top part of the DM3 pile to the end of loading, the movement at the toe is comparable to that of the shorter piles and smaller L/D ratios. This behaviour may be explained by interrogating the computed changes in the mean effective stress, p′, in the soil around the pile (Fig. 11) and the resulting void ratio, e (Fig. 12). The data are shown for the large-diameter DL2 pile at ultimate conditions, but are representative of all other analysed piles. As would be expected, the p′ increases in front of the pile (in relation to the direction of loading) near the top, and at the back of the pile at its toe, as shown by the dark/black contours in Fig. 11. On the contrary, areas at the back of the pile near the top and in front of the pile near the toe indicate a reduction in p′ due to unloading and, consequently, an increase in void ratio and loosening of the sand (Fig. 12) around the pile, leading to its more rigid response. This demonstrates the importance of employing a sand model that is capable of reproducing the effects of void ratio and stress level changes in the soil.

A further characteristic of material modelling is shown in Fig. 13, in terms of the mobilised loading conditions around the pile that are different from triaxial compression. The figure shows the spatial distribution of the Lode's angle, θ, for pile DL2 at ultimate conditions. As Dunkirk sand is normally consolidated, with K₀ < 1, the initial value of θ = −30° corresponds to triaxial compression. A rapid variation of θ is seen in the soil in front of the pile as it is loaded, engaging again the value of θ = −30° several diameters away. Different operational strengths, in terms of ϕ′, will be mobilised with respect to θ due to the shape of the yield surface in the deviatoric plane the model adopts (equation (9)), which justifies the above imposed limit on the opening of the bounding surface. Immediately behind the pile, the soil state is in triaxial extension (θ = +30°) due to the opening of a gap. In the analyses, the tensile capacity of the sand interface is enabled only by the prescribed initial suction in the interface elements above the groundwater table, as explained earlier and shown in Fig. 9. Consequently, the maximum depth of a gap is limited to 5·4 m, corresponding to the depth of the suction profile, and this is fully mobilised only for the largest, DL2, pile. For the medium piles, the depths of the gap are 1·45 m (DM7, L = 2·3 m), 2·5 m (DM4, L = 4·0 m) and 3·5 m (DM3, L = 6·1 m). It should be noted that gapping due to superficial suctions in the soil is irrelevant for offshore conditions, as the seabed soil is likely to be fully saturated.

Figure 14 compares the load–displacement curves predicted for all four test piles from the 3D FE analyses, to those obtained from field testing. These are presented for ground level horizontal displacements, v_G, up to 0·1D, where v_G is calculated as the average between the displacements of the leading and trailing edges of the pile, thus mitigating the impact of any possible ovalisation of the section due to loading. The analyses were not designed to account for strain-rate effects, therefore the predicted pile capacities should be compared with the end points of experimental holding stages, as these correspond to negligible strain rates. Apart from the shortest pile DM7, the 3D FE predictions for the remaining test piles are consistent and in good overall agreement with pile responses measured in the field, particularly when considering that these result directly from the interpretation of field conditions and laboratory soil behaviour (i.e. a back-analysis procedure of adjusting model parameters, by using pile test data, was not carried out to improve the numerical modelling results). The 3D FE analysis over-predicts the capacity of the shortest pile DM7 which, with its 2·3 m length, is fully embedded in the top 3·0 m thick layer of hydraulic fill. As discussed earlier, the behaviour of the fill material has not been sufficiently characterised during the Dunkirk testing campaign and, consequently, the reasons for this over-prediction are unclear and could not be quantified.

To demonstrate the appropriateness of the approach to integrate, in the constitutive model calibration, the field measurements of the maximum shear modulus, G₀, with triaxial test data, all four test piles were additionally analysed by adopting the constitutive model calibration from triaxial tests only (see Table 1) and the results are added to Fig. 14 (as ‘3D FE – lab only’). The impact of ignoring the field measurement of G₀ is obvious.

The early-stage loading (i.e. to a v_G = 0·01D) presented in Fig. 15 demonstrates again consistent predictions from FE analyses for all piles, albeit with a slightly softer response compared to the field data. Consistent with Fig. 14, the laboratory-only-based constitutive model calibration shows softer predictions of pile response.

The observed differences are reflected in the accuracy metric, η, used to assess the quality with which the measured load–displacement response is reproduced by the 3D FE model (considering here only the results with the adjusted set of constitutive model parameters)

where A_ref is the area below the reference load–displacement curve (i.e. the test data) and A_diff is the area of the region bounded by the reference and simulated curves, as shown in Fig. 16. The procedure followed for determining the reference curves is explained in the paper by McAdam et al. (2019). The accuracy metric is assessed at two distinct stages: for the overall load–displacement curve, η_0·1D, corresponding to ground-level displacements up to v_G = 0·1D; and for the initial loading, η_sd, which corresponds to small displacements of up to 0·001D. The results in Fig. 17 show that the 3D FE analyses achieve an average accuracy for the ultimate response, η_0·1D, of 81%, while that at small displacements, η_sd, is 72%.

The prediction of the embedded pile response compared to the field measurements is discussed here only for the largest DL2 pile as this is the longest of the tested piles (L = 10·5 m) and is therefore the least affected by the uncertainties in the behaviour of the top fill layer. The graphs in Fig. 18 are obtained by selecting stages of the analysis where either the ground-level displacements (Figs. 18(a) and 18(c)) or the ground-level moments (Figs. 18(b) and 18(d)) are the same as those in the field tests. Consequently, despite the good agreement between the computed and measured load–displacement curves, it should be recognised that the deformed shapes may not correspond to the same load, and equally the bending moment distributions may not correspond to the same ground-level displacement. The interpretation of the embedded pile response from field measurements is explained in the papers by Burd et al. (2019) and McAdam et al. (2019). For the 3D FE analyses the displacements were obtained directly from the nodal values, while the bending moments were obtained by summing the contributions from nodal forces and bending moment components of the shell elements at any given level.

One set of comparisons is shown in Figs 18(a) and 18(b) at early loading (v_G = 10·1 mm or M_G = 10·8 MNm). The 3D FE analysis predicts with good accuracy the measured deflected shape of the pile, although it overestimates slightly the depth of the point of rotation and the displacement at the pile base. The comparison of bending moments in Fig. 18(b) shows a good prediction from the 3D FE analysis. At significantly larger loads (v_G = 110 mm or M_G = 34·1 MNm), the deflected shape in Fig. 18(c) continues to be well predicted by the 3D FE calculations for most of the pile length, with a slight over-prediction of the toe-kick displacements around the pile base. The comparison of the bending moment in Fig. 18(d) is reasonable, with the 3D FE analysis under-predicting the moment profile in the top quarter of the pile when compared to the optimised structural model. The latter model, which results from the best fit between a structural model of the embedded pile and the measurements from strain gauges and inclinometers, is explained in the paper by Burd et al. (2019).

The paper presents a comprehensive computational study on the soil–structure interaction of laterally loaded piles installed at the PISA test site in Dunkirk (McAdam et al., 2019; Zdravković et al., 2019a). The 3D FE model, created in the FE code ICFEP (Potts & Zdravković, 1999), has demonstrated the following points.

• (a)

  Rigorous integration of soil data from laboratory and field investigations, applying sound engineering judgement, is a key to (i) establishing realistic initial ground conditions on the site and (ii) enabling detailed calibration of the chosen constitutive model. The latter is an enhanced version of the state-parameter-based bounding surface plasticity model developed by Taborda et al. (2014).

• (b)

  The constitutive model's advanced features, reflected in its ability to reproduce the dependency of sand behaviour on stress level and void ratio, are critical for the correct simulation of the observed behaviour of test piles, in terms of load–displacement curves as well as embedded deflected shapes and bending moments. In particular, the 3D FE predictions are systematically consistent at both small displacements and ultimate loads, for the selected range of pile geometries which include two diameters (0·762 m and 2·0 m) and three slenderness ratios (L/D = 3·0, 5·25 and 10·0).

• (c)

  The presented numerical approach establishes a systematic and comprehensive procedure for advanced modelling of laterally loaded piles in natural marine sands, aiming at providing a template for the use of complex constitutive models in boundary value problems. In effect, rather than performing back-analyses, this study focuses on the importance of void ratio, stress level and intermediate principal stress effects on the sand's evolving strength and stiffness.

The application of this FE model was key in the development of new 1D PISA design methodology for monopiles in sands, as 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.

To demonstrate the performance of the model using the two sets of calibrated model parameters listed in Table 1, five drained tests in triaxial compression are reproduced (see Zdravković et al. (2019a) for the laboratory testing programme), with varying initial mean effective stress and void ratio. The tests are summarised in Table 3 and comparisons shown in terms of ε_a–q and ε_a–ε_vol curves. As expected, the model calibrated only on triaxial data reproduces very closely the experimental behaviour in terms of initial stiffness and both the peak and ultimate strengths, as shown in Fig. 19. When model adjustments are applied in order to reproduce the shear stiffness measured in the field, the simulated stress–strain curves in Fig. 20, while showing initially stiffer response compared to experimental curves, as intended, retain a good overall reproduction of peak and ultimate strengths compared to experimental curves. The overall volumetric response is reasonably well reproduced in both cases, although the model predicts a larger effect of stress level on the volumetric behaviour for samples with the initial void ratio of about 0·64. The experimental data for these specimens show an unusually narrow range of the volumetric strain variation and therefore appear not to adhere entirely to well-established features of sand behaviour.

A

hardening modulus in equation (8)

A_diff

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

A_ref

area below reference curve in equation (10)

a₁, κ

model parameters in equation (4)

B

model fitting parameter in equation (5)

D

pile diameter

D^ep

plastic dilatancy ratio

D_R

relative density

e

void ratio

e₀

initial void ratio

e_0,ref

reference void ratio in equation (1)

e_CS

void ratio at critical state

e_min

minimum void ratio

e_max

maximum void ratio

e_peak

void ratio at peak

e_PT

void ratio at phase transformation

G₀

elastic shear modulus

G_tan

tangent shear modulus

g(θ, c)

deviatoric shape of model surfaces in equation (9)

H

applied lateral load

h

height of point of application of lateral load above ground level

h₀, h_b, h_e, h_f, h_g

parameters of the hardening modulus in equation (8)

K₀

earth pressure coefficient at rest

K_tan

tangent bulk modulus

k^b,d_c

model parameters for bounding (b) and dilatancy (d) surfaces in compression (c) in equation (2)

L

embedded pile length

M

bending moment in pile

m

size of the yield surface

p′

mean effective stress

p′_ref

reference pressure

p′_YS

size of the secondary yield surface

t

pile wall thickness

v

horizontal displacement of pile

v_G

horizontal displacement of pile at ground level

X, Y, Z

finite-element mesh coordinate system

z

depth below ground surface

η

accuracy parameter in equation (10)

θ

Lode's angle

λ

fitting parameter in equation (1)

M^b_c,e

gradient of the bounding surface (b) in triaxial compression (c) or extension (e)

M^c_c,e

gradient of the critical state surface (c) in triaxial compression (c) or extension (e)

M^d_c,e

gradient of the dilatancy state surface (d) in triaxial compression (c) or extension (e)

ν

Poisson's ratio

ξ

fitting parameter in equation (1)

ϕ′_mob

mobilised angle of shearing resistance

ϕ′_TXC

angle of shearing resistance in triaxial compression

ϕ′_TXE

angle of shearing resistance in triaxial extension

ψ

state parameter

ψ_peak

state parameter at peak

ψ_PT

state parameter at phase transformation