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

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]).

Figure 1

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: (2) p(y,z)=Apu tanh(kAPIzApuy) p(y,z) = A{p_u}\,\tanh ({{{k_{{\rm{API}}}}z} \over {A{p_u}}}y) 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: (3) kAPI=(0.008085 ϕ2.45−26.09) 103[kPa/m]for 29°≤ϕ≤45 ° (ϕ in degrees) \matrix{{{k_{{\rm{API}}}} = (0.008085\,{\phi ^{2.45}} - 26.09)\,{{10}^3}[{\rm{kPa}}/{\rm{m}}]} \cr {{\rm{for}}\,29^\circ \le \phi \le 45\,^\circ \,(\phi \,{\rm{in}}\,{\rm{degrees}})} \cr}

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. (4) Epy*=dp(y,z)dy|y=0=ApukAPIzApucos h2(kAPIzyApu)|y=0=kAPIz E_{py}^* = {\left. {{{dp(y,z)} \over {dy}}} \right|_{y = 0}} = A{p_u}{\left. {{{{{{k_{API}}z} \over {A{p_u}}}} \over {\cos \,{h^2}({{{k_{API}}zy} \over {A{p_u}}})}}} \right|_{y = 0}} = {k_{API}}z

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: (5) pu(z)=min{(C1z+C2Dp)γz(C3Dpγz) {p_u}(z) = \min \left\{{\matrix{{({C_1}z + {C_2}{D_p})\gamma z} \cr {({C_3}{D_p}\gamma z)} \cr}} \right.

The first expression corresponds to the wedge mode of failure, whereas the second one corresponds to the flow pattern. The coefficients C₁, C₂ and C₃ were given in both charts and tables in function of ϕ. The coefficient A for static loading is given by: (6) A=(3.0−0.8zDp)≥0.9 A = \left({3.0 - 0.8{z \over {{D_p}}}} \right) \ge 0.9

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, 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.

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.

Figure 3

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 showing a similar behavior. Part figure (b) from each set of figures contains the results of the BNWF model when WILDPOWER 1.0 uses Soensen et al. (2010), Sorensen (2012), and Wiemann et al. (2004) p–y curves. These methods, close to each other, showed a different tendency than that of API. Part figure (c) in Figures 4–7 corresponds to the results which showed an entirely different pattern. These were obtained when WILDPOWER 1.0 was executed with the Kallehave et al. (2012)’s p–y curve.

Figure 4

Figure 5

Figure 6

Figure 7

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, almost a perfect match is observed between the FE displacements and those of Winkler model for the upper part of the monopile when WILDPOWER 1.0 employs Sorensen et al. (2010) p–y curve, albeit a significant deviation is noticed at the monopile tip.

The Winkler model using both Sorensen (2012) and Wiemann et al. (2004) p–y curves reveals an underestimation of pile displacements along the entire monopile length, thus showing a deformation pattern close to that exhibited by the API method. The profile of displacements shown by Kallehave et al. (2012) failed to capture the behavior of a rigid monopile completely.

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 al. (2012) model in 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 Figure 4. In the API’s group (Figure 8a), values of K_L by both Reese et al. (1974) and O’Neill and Murchison (1983) are more than twice that of ABAQUS. The two other stiffness coefficients are also much higher than those of FE analysis. This indicates the failure of the API models in producing accurate soil–monopile stiffness. Two important points are worth noting from the examination of Figure 8b. Firstly, an excellent agreement is observed between the results provided by Sorensen et al. (2010) and those of FE analysis. The computed deviations revolve around 2.6%, 12%, and 13% for K_L, K_R and K_LR, respectively. This confirms the good performance of Sorensen et al. (2010)’s model in analyzing large-diameter monopiles. Secondly, the results provided by Wiemann et al. (2004) and Sorensen (2012) 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 large-diameter monopiles. Figure 8c compares the stiffness coefficients provided by Kallehave et al. (2012)’s model with those given by the FE analysis. The overestimation of soil–monopile stiffness is clearly visible and reaches 74%, 43%, and 63% for K_L, K_R and K_LR, respectively, making Kallehave et al.’s model unsuitable for large-diameter monopiles under lateral loads.

Figure 8

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.

3.2.1

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 soil–monopile 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: (7) fFB=12π3EITS[Mtop+(33/140)mTS]L3 {f_{{\rm{FB}}}} = {1 \over {2\pi}}\sqrt {{{3E{I_{TS}}} \over {[{M_{{\rm{top}}}} + (33/140){m_{TS}}]{L^3}}}} (8) mTS=mT+ms {m_{TS}} = {m_T} + {m_s} 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: (9) EITS=(1−α)EIs+αEItop×λav E{I_{TS}} = (1 - \alpha)E{I_s} + \alpha E{I_{{\rm{top}}}} \times {\lambda _{\rm av}} 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.

Figure 9

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: (10) f1=CRCLfFB {f_1} = {C_R}{C_L}{f_{{\rm{FB}}}} 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.

3.2.3

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₁ (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: (11) {uLθR}=[ILILRIRLIR]{HM} \left\{{\matrix{{{u_L}} \cr {{\theta _R}} \cr}} \right\} = \left[ {\matrix{{{I_L}} & {{I_{LR}}} \cr {{I_{RL}}} & {{I_R}} \cr}} \right]\left\{{\matrix{H \cr M \cr}} \right\} 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.

Figure 10

Table 6

Monopile head flexibility coefficients.

Flexibility coefficients	ABAQUS FE analysis	WILDPOWER 1.0
		Reese et al. (1974)	O’Neill and Murchison (1983)	Kallehave et al. (2012)	Sorensen et al. (2010)	Sorensen (2012)	Wiemann et al. (2004)
IL (m/MN)	0.003356	0.002011	0.001775	0.001052	0.003225	0.001871	0.002375
IR (rad/MN m)	0.000048	0.000041	0.000039	0.000033	0.000047	0.000039	0.000043
ILR (1/MN)	0.000314	0.000225	0.000207	0.000142	0.000298	0.000210	0.000252

The stiffness matrix elements are determined by simply inverting the flexibility matrix of equation (11): (12) KL=IRILIR−ILR2KR=ILILIR−ILR2KLR=ILRILIR−ILR2 \matrix{{{K_L} = {{{I_R}} \over {{I_L}{I_R} - I_{LR}^2}}} & {{K_R} = {{{I_L}} \over {{I_L}{I_R} - I_{LR}^2}}} & {{K_{LR}} = {{{I_{LR}}} \over {{I_L}{I_R} - I_{LR}^2}}} \cr}

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. By taking a close look at the reported numbers, three important features may help describing the performance of each p–y model. Firstly, the reported values of both C_R and C_L suggest different contributions. The values indicate that the monopile–soil interaction is mostly quantified by the former, whereas the latter has a marginal contribution as it is very close to unity regardless of the method of analysis used. Secondly, the 1st NFs of both FEA and WILDPOWER 1.0 when using Sorensen et al. (2010) p–y curve almost coincide and are both very close to the measured 1st NF supplied by the manufacturer. This suggests that this model is the best for designing monopiles when it comes to place an OWT’s 1st NF within a safe interval. Thirdly, the other Winkler models present more pronounced deviations, with the highest shown by Kallehave et al.’s model (8.0%) and the lowest by Wiemann et al.’s model (4.86%).

Table 7

Stiffness coefficients, interaction factors, and the different 1st NFs computed for the OWT at North Hoyle.

Parameters	ABAQUS FE Analysis	WILDPOWER 1.0
		Reese et al. (1974)	O’Neil and Murchison (1983)	Kallehave et al. (2012)	Sorensen et al. (2010)	Sorensen (2012)	Wiemann et al. (2004)
KL (MN/m)	771.64	1293.02	1471.54	2307.26	744.07	1335.67	1101.32
KR (MN m/rad)	54,172.94	63,519.50	66,491.98	74,640.81	50,676.57	63,667.97	60,458.79
KLR (MN)	5065.43	7109.84	7770.27	10,061.80	4689.59	7141.61	6412.90
CR	0.8693	0.8859	0.8900	0.9072	0.8703	0.8901	0.8802
CL	0.9977	0.9986	0.9988	0.9993	0.9978	0.9987	0.9984
Fixed base NF fFB (Hz)	0.417
1st NF f1 (Hz)	0.361	0.369	0.371	0.378	0.362	0.371	0.367
Measured NF f measured (Hz)	0.35
Deviation (%)=|fmeasured−f1fmeasured| 100 {\rm{Deviation}}\,(\%) = \left| {{{{f_{{\rm{measured}}}} - {f_1}} \over {{f_{{\rm{measured}}}}}}} \right|\,100	100 3.14	5.43	6.00	8.00	3.43	6.00	4.86