Evaluation of sand p–y curves by predicting both monopile lateral response and OWT natural frequency

: Extending the use of the p–y curves included in the regulation codes API and DNV to design large-diameter monopiles supporting offshore wind turbines (OWTs) was unsuccessful as it resulted in an inaccurate estimation of the monopile behavior. This had prompted many investigators to propose formulations to enhance the performances of Winkler model. In this paper, two case studies are considered. A case consisting of an OWT at Horns Rev (Denmark) supported by a monopile in a sandy soil was studied first. Taking the FEA using ABAQUS as reference, results of WILDOWER 1.0 (a Winkler computer code) using the recently proposed p–y curves giving design parameters were plotted and evaluated. In order to see the ability of proposed p–y curves to predict the monopile head movements, and consequently the first natural frequency (1st NF), a second case study consisting of a monopile supporting an OWT at North Hoyle (UK) was selected. The monopile head stiffness in terms of lateral, rocking, and cross-coupling stiffness coefficients, necessary for the 1st NF, were computed using both ABAQUS and WILDPOWER 1.0. Comparisons with the measured 1st NF showed that with the exception of one p–y model, none of other proposed Winkler methods is able to predict accurately this parameter.


Introduction
Complicated environmental issues combined with the growing harm of the greenhouse gas emission placed the environmentalists in a challenging situation all over the world. The search for clean energy from renewable resources became an urgent need. This has been confirmed during the United Nations climate change and energy policy framework conference where it was agreed on the target that 27% of the total energy consumed in the EU must be renewable by 2030 (EUCO 169/14) [1]. To reach this prospective goal, more significance is given to the development of offshore wind parks as the wind energy produced by offshore wind turbines (OWTs) is considered more attractive for the availability of large offshore areas, where consistent winds are steadier and stronger than those in onshore lands.
The most widely adopted foundation is the monopile, which is a large-diameter hollow steel driven pile with a diameter D p ranging from 3 to 8 m and a length to diameter (L p ⁄D p ) smaller than 8. This is because this type of foundation is easy and convenient to construct [2,3]. However, monopile diameters increased considerably over the past decades and are expected to grow further in the future due to the increasing rotor diameters required for more efficiency in wind power.
It is well recognized that the Winkler model that uses the nonlinear p-y curves proposed by the standards, currently in use [4,5] for monopile design, has gained broad confidence over many decades. However, as this model has been developed on the basis of field testing investigations performed on small-diameter piles (D p =0.61 m) with an (L p ⁄D p ) ratio of 34.0 [6,7], it does not appear capable of reasonably predicting the lateral response of large-diameter monopiles supporting OWTs. For an in-depth review of the drawbacks of the API p-y model, the reader is referred to [8].
The main objective of this paper is to assess the performances of the currently used p-y models (API and DNV) and the enhanced p-y curves recently proposed to design monopiles under lateral loading. This is achieved by a comparative study between the results supplied by a computer program called WILDPOWER 1.0 [8] and those of the Finite Element (FE) analysis obtained using the commercial package ABAQUS.
Two case studies are considered in this paper to evaluate the Beam on a Nonlinear Winkler Foundation (BNWF) model and to study its performances when applied to large-diameter monopiles. The first, which consists of an OWT embedded at Horns Rev (Denmark), is believed to constitute a severe test to both computer programs involved in the comparative study. The monopile lateral behavior is examined, and the assessment is carried out through the profiles of displacements, bending moments, shear forces, and soil reactions.
It is generally accepted that the dynamic characteristics, especially the first natural frequency (1st NF) of an OWT, are affected by the soil-monopile interaction [2,3]. When the BNWF model fails to adequately capture the lateral behavior of a large-diameter monopile, an inaccurate monopile head stiffness is produced, and consequently the OWT design is unsafe. In this context, a second case study consisting of an OWT supported by a monopile at North Hoyle (UK) is studied. This case is selected for the availability of the OWT tower structural details along with the geotechnical and monopile data required for Winkler and FE modeling, and most importantly for the existence of the OWT-measured 1st NF, which forms a key element in this assessment study.
Three work steps are necessary for reaching the 1st NF and then assessing the performance of each p-y model implemented in WILDPOWER 1.0. Firstly, an analytical expression for the 1st NF encompassing the monopile head stiffness (lateral, rocking, and cross-coupling coefficients) is provided and detailed. Secondly, the monopile head stiffness quantifying the soil-monopile interaction is described. Thirdly, the natural frequencies supplied by both FE and Winkler analyses are compared to the measured ones. Comments are given at the end of the paper.

The BNWF model for designing large-diameter monopiles in sands
In current geotechnical practice, analysis of laterally loaded piles is commonly conducted using an approach in which the soil foundation system is modeled as a BNWF with a series of uncoupled nonlinear springs representing the soil. This model involves the resolution of fourth-order differential equation: In this equation, y is the lateral deflection of the monopile at a point z along the monopile length, p is the soil reaction in force unit per unit length, (EI) p (z) is the monopile flexural rigidity which varies with depth as in the case of segmented or tapered monopiles, and V is the applied vertical load at the pile head.

Reese et al. (1974) and API p-y curves
To perform accurate design of piles supporting platforms in gas and oil offshore industry, Reese et al. [6] carried out a total of seven (two static and five cyclic) load tests on two identical instrumented piles installed at Mustang Island, Texas. In the static tests, both piles had a diameter of 0.61 m and a slenderness ratio of 34.4 and were installed in a medium sand. Based on these field tests, Reese et al. (1974) proposed a semi-empirical p-y curve consisting of three straight lines and a parabola ( Figure 1a). For lack of space, the different parameters involved in Figure 1a are not given here. However, for a detailed explanation, the reader can refer any article from the overwhelming number of papers that have been dedicated to this important subject (see, e.g. [9]). To improve Reese et al.'s formulation which has been characterized by a significant amount of empiricism, O'Neill and Murchison [7] proposed a hyperbolic p-y curve in the form: with k API representing the initial coefficient of subgrade reaction depending on the angle of internal friction or the relative density of the cohesionless soil ( Figure 1b). k API in equation (2) has been fitted by [10] as: The initial stiffness of the API p-y curve E * py recommended in the design can be obtained by evaluating the slope of the p-y curve tangent at y=0.
It is quite clear from equation (4) that the initial stiffness is independent of monopile properties (diameter and bending stiffness) and is linearly dependent on the depth z.
The ultimate resistance p u appearing in equation (2) is the lesser value of the two expressions proposed in [6]. These expressions were simplified by [11] as: The first expression corresponds to the wedge mode of failure, whereas the second one corresponds to the flow pattern. The coefficients C 1 , C 2 and C 3 were given in both charts and tables in function of ϕ. The coefficient A for static loading is given by:

Recent propositions of p-y curves to improve the BNWF model performances in sands
As building OWTs supported by monopiles became a necessity, the monopile designers extrapolated the API methodology to design large-diameter monopiles. However, poor predictions were observed, and therefore, .0 m b=0.6, c=0.5, d=3.6, ϕ in radians

Silica
Sorensen et al. [13] = � � Kallehave et al. [14] = � � Sorensen [15] the above-mentioned p-y curves were deemed no longer appropriate for this kind of foundations. A detailed critical review of the API method can be found in [8]. This situation prompted many researchers to suggest analytical p-y curves formulae in an attempt to improve the analysis by accounting for the monopile diameter and sometimes for monopile bending flexural rigidity. Hence, most of the suggested formulations alter the p-y curves at the level of initial stiffness, whereas the ultimate soil resistance p u in [6] formulation was kept unchanged. These different propositions are listed in Table 1.
The analytical formulations appearing in Table 1 are all for silica sands and the only input soil parameter required is the sand internal friction angle ϕ, excepting the formulation by Sorensen (2012), which additionally necessitates the sand Young's modulus E s .
The multi-segment p-y formulation proposed by Reese et al. [6] (Figure 1a), the hyperbolic form by O'Neill and Murchison [7] (Figure 1b), and the other formulations suggested to enhance the p-y curves described in Table 1 have been encoded in a Fortran program called WILDPOWER 1.0. This has been conceived as a general name which stands for Winkler Idealization of Large-Diameter Piles for Offshore Wind Energy Regenerators version 1.0. The implementation of a new p-y curve [16] has been already performed, and the development of a user graphic interface (GUI) for WILDPOWER 1.0 using VB.net is underway.

Prediction of large-diameter monopile response: BNWF modeling against Finite Element Analysis (FEA)
The performances of p-y curve models implemented in WILDPOWER 1.0 are assessed against the FE results provided using ABAQUS [17].
The FE mesh used to model the monopile and its surrounding medium is illustrated in Figure 2. Due to the existence of a plane of symmetry in the considered problem, only half of the domain is meshed as shown in Figure 2. The soil domain is discretized using 4600 twentynoded hexahedral displacement-based isoparametric elements type C3D20R. The vertical boundaries are placed at a distance of 15 D p (D p is the monopile diameter) from the monopile centerline to eliminate any possible boundary effects, while the bottom boundary is placed far enough under the monopile tip at a distance of one L p (L p is the monopile length). Nodes belonging to the vertical boundaries are prevented from lateral movements, whereas those of the bottom boundary are fixed against any movement. The monopile is modeled either as a cylindrical solid or as an open tubular monopile using shell elements. When modeled as a solid, the monopile Young's modulus provided as input data in ABAQUS is calculated based on the effective value of bending stiffness (EI) p , considering the real monopile cross section.
To disregard any installation effects, the monopile was assumed to be wished-in-place. Using the same type of elements as those used to model the soil, 3230 twentynoded quadratic hexahedral elements were employed to model the pile. Furthermore, the pile was treated as a linear elastic material with a Young's modulus of E p =210 GPa and a Poisson's ratio ν=0.30. The elasto-plastic response of sandy soil was simulated using the Mohr-Coulomb model with non-associated flow rule. Five parameters are required in the Mohr-Coulomb model for performing a FEA. The strength parameters, which are the cohesion c, the angle of internal friction ϕ, and the dilation angle ψ, control the plastic behavior of soil, whereas the stiffness parameters, namely, the soil Young's modulus E s and the soil Poisson's ratio v s govern the soil elastic deformation. To model the interaction between the monopile and the soil, surface-tosurface master-slave contact formulation was used. Since the soil in contact with the monopile is much softer than the monopile, it was selected as the slave surface and the monopile was selected as the master surface.

BNWF model against FE analysis for a monopile supporting an OWT at Horns Rev (Denmark): Design parameters estimation
The accuracy of the monopile response to lateral loading is estimated by studying a monopile embedded at a site called Horns Rev, which is situated in the Danish sector of the North Sea. The monopile design parameters, which consist of displacement, bending moment, shear force, and soil reaction profiles, are computed first by both WILDPOWER 1.0 using the six p-y curves and ABAQUS and then compared. The monopile considered in this subsection has been chosen from one of the 80 installed OWTs type Vestas V80. It has an outer diameter of 4 m, a length of 31.6 m, and a wall of variable thickness inducing a variable flexural stiffness (EI) p along the monopile shaft as shown in Figure 3a. The monopile has been driven to 31.8 m below the mean sea level, resulting in an embedded depth of 21.9 m.
The identification of soil strata crossed by the monopile has been carried out by Augustesen et al. [10]. The Cone Penetrometer Testing (CPT) experiments revealed that the lithology was mostly formed of sand with roughly six layers. Details of layer thicknesses and their physical and strength properties are summarized in Table  2. The nature and thickness of each layer are illustrated in Figure 3b.
The monopile in this first case study was subjected to a static horizontal load of H=4.6 MN and an overturning moment M=95 MN m, both acting at the seabed level adopting the same loading considered in [10]. Results of both WILDPOWER 1.0 and ABAQUS are depicted in the four sets of Figures 4-7, which describe the variation of lateral displacements, bending moments, shear forces, and soil reactions, respectively, with depth. A close inspection of the results in their gloabality reveals that the obtained outcomes can be categorized in three main distinct classses. So, part figure (a) from each design parameter corresponds to the results of p-y curves already adopted by the regulation codes (API and DNV) and    Figure 4 shows the variation of the monopile displacements with depth. Augustesen et al. [10] had previously studied the monopile at Horns Rev using FLAC 3D . The obtained displacement profile is included in Figure 4 for comparison. With the exception of small deviations observed at both monopile head and tip, a perfect agreement is noticed between FLAC 3D 's displacement profile and that of ABAQUS for almost the entire length of the monopile. This is somehow a validation of the FE analyses carried out by ABAQUS in this paper.
At the first glance, it easy to notice that the pattern exhibited by the Winkler model, when both Reese et al. and the API p-y curves are used, differs from that of ABAQUS. This confirms that these methods are not able to predict correctly the monopile displacements, especially at the top and bottom of the monopile. However, in Figure  4b The different bending moment profiles are plotted in Figure 5. With the exception of Kallehave et al.'s model which underestimates the bending moment on almost the entire monopile length, it is clearly seen that all methods exhibit an identical tendency with nearly the same value of maximum bending moment, which differs slightly from that of ABAQUS. This gives a confirmation that the bending moment cannot be considered as a performance indicator to assess the BNWF model. The profiles of shear forces are shown in Figure 6. With respect to shear forces, all the methods capture well the monopile behavior over the upper third length, while a significant ovestimation of shear force values is observed in the middle part, where the absolute maximum efforts are expected to take place. Figure 7 shows the distribution with depth of soil reaction for all studied models. Three important key points are worth indicating. Firstly, both FE and Winkler models exhibit a similar variation pattern, where the magnitudes of deviations are identical to those observed in the displacement profiles. This is because the lateral displacements and soil reactions are the direct solutions of the differential equation representing the Winkler model. Secondly, the depth locations at which the soil reaction becomes zero differ from one model to another. For the API's group (Figure 7a), the zero soil reaction is situated slightly above that of ABAQUS, whereas in the Sorensen's group (Figure 7b), the depth of zero soil reaction almost coincides with that of FE analysis. However, Kallehave et    Figure 7c shows a point located well above the depth of zero soil reaction provided by ABAQUS. Thirdly, the soil reaction increases in the part of the soil located beneath the soft layer (formed of organic clay) and extending to the monopile tip. This is observed in all methods with the exception made for both API and Kallehave et al. models which show a decreasing pattern. The soil-monopile interaction is a necessary ingredient for a rigorous assessment of the vibration characteristics of an OWT whose quantification relies fundamentally on monopile head stiffness coefficients. In this regard, the load deformation curves for the monopile studied in this section were established first, fitted by sufficiently higher-degree polynomials, and then the first derivative was evaluated at the origin for the determination of flexibility coefficients. The obtained coefficients were inverted to get the stiffness coefficients. Histograms of Figure 8 show the values of the stiffness coefficients of the three groups previously seen in this case study. As it has been expected, the values of stiffness coefficients confirm the tendency of lateral displacements of are very close to each other, but both significantly deviate from those of ABAQUS. The estimated ranges of deviations are 35%-38%, 22%-26%, and 32%-36% for K L , K R and K LR , respectively. Since the stiffness is highly overestimated, these models are not appropriate for designing largediameter monopiles. Figure 8c

BNWF model against FE analysis for a monopile supporting an OWT at North Hoyle (UK): Wind tower natural frequency assessment
The most significant element in this subsection is the 1st NF, which is a parameter of importance in the safe design of slender structures such as offshore wind towers. Therefore, and before performing the comparisons based on the 1st NF of each model, three preliminary steps are worth to be described. In the first step, the different terms constituting the analytical expression describing the natural vibration of an OWT are given and explained. In the second step, the structural details of the North Hoyle OWT selected for this study are given and the site geotechnical data are provided as well. The soil-monopile interaction in terms of flexibility coefficients and then the resulting monopile head stiffness coefficients necessary for the natural frequency computation are established by both ABAQUS and WILDPOWER 1.0 using all Winkler models in the third step. The 1st NF computed by the FE analysis and by the Winkler models are compared to the measured natural frequency supplied by the OWT manufacturer to assess the performance of each p-y curve.

Approximate OWT modeling for dynamic analysis
The natural frequency of an OWT is highly dependent on the material properties employed in its fabrication and is significantly affected by the stiffness provided by the soilmonopile system. The structure mounted on the monopile is usually formed by a substructure and tower in an OWT (Figure 9). A rigorous expression for the fixed base natural frequency, taking into account the variation of the tower cross section with elevation, has been proposed in [18] as: where m T and m S are the tower and substructure masses, respectively. For more precise representation of the (substructure + tower) flexural rigidity EI TS , Amar Bouzid et al. [3] and Aissa et al. [2] proposed the following formula: where α=L T ⁄L, EI S is flexural rigidity of the substructure, EI top is the top tower bending stiffness, and λ av is a factor accounting for the tower tapering in terms of m=D b ⁄D t . The expressions for computing m T , m S , EI S , EI top and λ av are all given in Appendix A.
In this kind of structures, the monopile head stiffness quantifying the soil-monopile interaction can be better represented by three springs. This is because the translation and the rotational modes of vibration are cross coupled, which requires an additional spring. In this context, Arany et al. [19] derived the expressions of natural frequency of an OWT on three-spring flexible foundations by means of two beam models: Bernoulli-Euler and Timoshenko. The natural frequency in both cases was obtained numerically from the resulting transcendental equations. They proposed a closed form expression containing a lateral stiffness coefficient K L , a rocking stiffness coefficient K R and a cross-coupling stiffness coefficient K LR . The equation for the 1st NF proposed is: where f FB is the fixed base frequency given by equation (7). The factors C R and C L accounting for the stiffness provided by the monopile (K L , K R and K LR ) are functions of the tower geometrical properties. Their analytical expressions are given in Appendix B.

Structural OWT and geotechnical site details at North Hoyle (UK)
The geotechnical investigations at the North Hoyle site revealed a lithology formed mostly of sand layers. The adopted soil strength and deformation parameters are summarized in Table 3. The North Hoyle OWT structural data and measured 1st NF supplied by the manufacturer are presented in Table 4. Based on the data provided in Table 4, other parameters necessary for the 1st NF computation are given in Table 5 in detail. Here, an important point is worth noting. The value of soil Young's modulus reported in Table 5 was deemed extremely high for such a soil and a value of E s =110 MPa was adopted in the present computations to mimic the sandy layer at Horns Rev site which has a close angle of friction in Table 2 of the first case study in Section 3.1.

Approach to compute the OWT natural frequency and comparison
The last step leading to the comparative study is to evaluate the OWT natural frequency using FE by ABAQUS on one hand and the 1st NF using WILDPOWER 1.0 for each p-y curve model on the other hand. FE computations using ABAQUS are based on the same FE mesh used for the first case study in Section 3.1. Although finding f FB is easy to perform since equation (7) depends only on the OWT structural parameters, finding f 1 (in equation 10) is not a straightforward task. The latter depends on the values of stiffness coefficients K L , K R and K LR . These coefficients cannot be evaluated by a direct process in boundary value problems controlled by forces (FE or    Winkler approach). However, the flexibility coefficients at the monopile head are the direct results of any structural method controlled by forces because the monopile head movements (displacement and rotation) are related to the applied loading by the following equation: where H and M are the shear force and the overturning moment applied at the monopile head, respectively, and u L and θ R are the lateral displacement and the rotation of the monopile head, respectively. The curves H-u L , M-θ R and H-θ R are respectively plotted in Figure 10a-c. The flexibility coefficients determined by both ABAQUS and the p-y models included in WILDPOWER 1.0 are presented in Table 6. Once again, the most flexible behavior is exhibited by the FE analysis, which is seen either in translational or rocking mode. However, when the p-y curve of Kallehave et al. (2012) is considered, the flexibility coefficients have the lesser values, indicating thus the stiffest behavior.
The stiffness matrix elements are determined by simply inverting the flexibility matrix of equation (11): The monopile head stiffness coefficients, C R and C L and the values of 1st NF generated by the FEA and the BNWF model are grouped in Table 7

Conclusions
Although successful in designing small-diameter piles, the extension of the p-y method already included in the piling regulation codes to the design of large-diameter monopiles supporting OWTs resulted in unsatisfactory predictions and showed that the API p-y curves are not able to predict the correct deformation pattern of this kind of stiff foundations. Therefore, many researchers have been prompted to propose other formulations of p-y curves which have the ability to produce better performances.
To assess the performances of recently proposed p-y curves, two case studies were considered in this paper. The first case study, which consisted of a monopile embedded in multi-layered sandy soil at Horns Rev (Denmark), was considered to illustrate the evolution of all the monopile design parameters with depth. A monopile supporting an OWT at North Hoyle (UK) was chosen as a second case to evaluate the Winkler model performances by computing      The analytical expressions for both substructure and tower masses are given by: where γ steel is the steel density usually taken as 7860 kg⁄m 3 , while the other parameters are detailed in Figure 9c. The flexural rigidities for both substructure and tower top are given by: where E is the steel Young's modulus, which is equal to both E T and E p . The factor λ av , which depends on the ratio m, has been determined by Aissa

Appendix B: Factors quantifying the soil-monopile interaction
The factors C L , C R appearing in equation (10) are given by: where a=0.6, b=0.5, and η L , η R and η LR are given by: