Assessment of natural frequency of installed offshore wind turbines using nonlinear finite element model considering soil-monopile interaction

Djillali Amar Bouzid,Suhamoy Bhattaharya,Lalahoum Otsmane

aDepartment of Civil Engineering,Faculty of Technology,University Saad Dahled of Blida,Route de Soumaa,Blida,09000,Algeria

bDepartment of Civil and Environmental Engineering,Tomas Telford Building,University of Surrey,Surrey,GU2 7HX,UK

cDepartment of Material Engineering,Faculty of Sciences and Technology,University of Médéa,Quartier Ain D’Hab,Médéa,26000,Algeria

1.Introduction

Wind,solar power and geothermal heat are representative of clean and renewable energy which have the potential to become alternatives to current supplement of fossil fuel sources of energy in the future.While these alternative energy sources have their advantages and drawbacks,wind energy is widely accepted as the cheapest and most economically available one based on current technology.Today,wind energy has proven to be a valuable feature for large-scale future investment in the energy industries worldwide,and many countries install their proper wind turbine generators(WTGs),mainly on land.As offshore wind turbines(OWTs)have gained their popularity,many WTG manufacturers believe that offshore wind energy will play an increasingly important role in the future development.This is supported by the fact that all principal wind turbine manufacturers currently are spending huge amount of money and effort on developing larger offshore WTGs for deeper waters where wind speed is generally higher and steadier,resulting in an increase in energy output.

Although there are many OWTs support options which may range from gravity foundations(for shallow depths of 0-15 m)to floating foundations(for very deep waters of 60-200 m)(Achmus et al.,2009;Lombardi et al.,2013;Damgaard et al.,2015;Abed et al.,2016),most OWTs are supported on monopile foundations,as they are simple structures which are easy and convenient to construct.The accumulated experience from limited monitored data from OWTs over the last 15 years showed that the available design procedures(mostly contained in the API(API and ISO,2011)and DNV(DNV-OS-J101,2004)regulation codes suffer limitations.

The existing methods were established/calibrated by testing small-diameter piles used for supporting offshore platforms in gas and oil industry,often with design criteria and loading conditions which are different from those encountered in an OWT.The inappropriateness of these methods comes from the fact that:

(1)The continuum(soil)is replaced by a series of uncoupled springs.However,reliable results necessitate a rigorous method which can properly account for the true deformation mechanism of soil-monopile interaction.

(2)As they rotate freely,monopiles supporting OWT energy converters undergo severe degradation in the upper soil layer resulting from cyclic loading,whereas offshore jacket piles are significantly restrained against rotation at their heads.

(3)Monopiles are relatively shorter and rigid piles with a length to diameter ratio(Lp/Dp)in the range of 2-6 and a diameter(Dp)of up to 8 m envisaged for the next generation of turbines,whereas offshore piled foundations in the offshore oil and gas industry have a length to diameter ratio(Lp/Dp)of over 30 and relevant recommendations have been set on the basis of full-scale loading tests on long,slender and flexible piles with a diameter of 0.61 m(Reese et al.,1974).

(4)The API model is calibrated in response to a small number of cycles for offshore fixed platform applications.However,an OWT over its lifetime of 20-25 years may undergo 107-108cycles of loading.

Due to these complex issues,appropriate determination of the dynamic characteristics of these extremely complex structures through their monopiles head stiffnesses is continuing to challenge designers,as the foundation of an OWT behavior is still not well understood,and also not introduced in the current design guidelines.

Concerning accurate prediction of the monopile head stiffnesses,numerical analysis using the finite element method(FEM)constitutes an excellent alternative to capture the real behavior of this type of foundations and hence to accurately estimate the dynamic characteristics of an OWT.

Otsmane and Amar Bouzid(2018)formulated a nonlinear pseudo three-dimensional(3D)computation method,combining the FEM and the finite difference method(FDM).They wrote a Fortran computer code called NAMPULAL(nonlinear analysis of monopiles under lateral and axial loadings)to study monopiles under axial,lateral and moment loadings in a medium characterized by the hyperbolic model for representing the stressstrain relationships.In this paper,we attempt to apply NAMPULAL to examining the lateral behavior of monopiles supporting OWTs chosen from five different offshore wind farms in Europe.These offshore wind farms include Lely A2(UK),Irene Vorrink(Netherlands),Kentish Flats(UK),Walney 1(UK)and Noth Hoyle(UK).

To accurately estimate the natural frequency of the OWT structure(tower+substructure)which is a function of monopilesub soil interaction,the monopiles head movements(displacements and rotations)and consequently,the lateral stiffnessKL,the rotational stiffnessKRand the cross-coupling stiffnessKLRare obtained and substituted in the analytical expression of natural frequency for comparison.In general,the results of comparison between the computed and measured natural frequencies showed a good agreement.

2.Natural frequency and modal analysis

OWTs are dynamically sensitive structures,in which the dynamic soil-structure interaction is a pivotal aspect of their design process and consequently,they require accurate soil stiffness estimation in order to ensure that the design frequency matches the actual operational frequency when the wind turbines are constructed.

The natural frequency of the hub-tower-foundation system is the key feature on which the response of an OWT to wind and wave loads depends.This is due to the dynamic nature of the loads on the wind turbine structure and the slenderness of the system.Through determination of the natural frequency,designer can assess the strains produced by loading cycles,through which the fatigue failure of the structure can be ascertained.Therefore,an accurate estimation of this parameter is essential to assess the working life of a wind turbine.

Unlike most large-scale civil engineering structures,wind turbines are subjected to millions of periodic excitation cycles during their operating life.The rotor spinning at a given velocity induces mass imbalances(gyroscopic effect),causing a frequency known as 1P.In addition to this,the effect of a standard turbine havingnblades induces a further excitation due to the blades passing the tower.The frequency of this shadowing effect isnP,wheren=3 in most cases.

The modern installed wind turbines are characterized by a range of different velocities in which their rotors are operating.This results in two ranges of operating frequencies around 1Pand 3P.In order to avoid resonance,the natural frequency of the tower cannot be in any of these two ranges and must be far from 1Pand 3P.

The OWT design can be performed in such a way that the first eigenfrequency lies within three possible ranges:soft-soft,soft stiff and stiff-stiff as shown in Fig. 1.

(1)Soft-soft range:the natural frequency is less than the lower bound of 1P.This implies that the structure is too flexible,and moreover,this is a range where the frequency of waves may lie,therefore leading to resonance.

(2)Stiff-stiff range:this is a range where the tower frequency is higher than the upper bound of blade passing frequency(3P).This range is economically unfeasible as it leads to a too rigid(heavy and expensive)structure,making it inappropriate for design.

(3)Soft-stiff range:in this interval,the natural frequency lies between 1Pand 3P.This range is the optimum range for the best possible design.

The system stiffness must be such that the natural frequency of the wind turbine does not lie within the rotor frequency excitation bands,as this may induce resonance which could lower the design life significantly.

In order to satisfy these requirements and to keep the natural frequency of the whole structure in the adequate margin of the soft-stiff range,thus avoiding resonance,a joint effort between foundation designers and turbine manufacturers is performed.Foundation designers need careful site investigations to obtain reliable soil data in order to correctly assess the foundation stiffness.

2.1.Appropriate OWT modeling for dynamic analysis

The natural frequency of a wind turbine is highly dependent on the material properties used in its construction,and is significantly affected by the stiffness of the soil surrounding the monopile.Assessment of foundation stiffness is the key to obtain reliable estimate of system frequency.

In the computation of eigenfrequencyf1,most researchers in the past tried to model this complex system according principally to two concepts(Prendergast et al.,2015;Yi et al.,2015).In the first one,Yi et al.(2015)simply considered the soil as a medium having an infinite stiffness.In this regard,Vught(2000)used a model in

Fig. 1.Forcing frequencies against power spectral density for a three-bladed wind turbine(Hz).

which the wind turbine is considered as an inverted pendulum having a flexural rigidityEI,a tower mass per metermTand a top massmt.Expression for the first eigenfrequency is given by

whereLTis the tower height.The natural frequency expressed by Eq.(1)is based on a uniform tower cross-section.A slight different expression has been proposed by Blevins(2001):

It seems from the first sight that these equations are based on a simple model which ignores the fact that the tower is tubular and conical in shape and generally does not have a constant wall thickness.Additionally,Eqs.(1)and(2)depend only on the tower geometrical and mechanical properties without taking into account the OWT foundation characteristics.This physically does not make sense,as the monopile head movements that occur as a result of the applied loading on the tower can lead to a finite stiffness of the monopile-subsoil system.This obviously has an influence on the value of the first natural frequency and shows that the service limit computations based on Eqs.(1)and(2)are inaccurate.

Bearing in mind that the dynamic analysis of the whole system composed of tower-monopile-soil is hard to perform,Prendergast et al.(2015)tried to find a natural frequency expression considering soil stiffness(Zaaijer,2006;Yu et al.,2014;Prendergast et al.,2015).Alternatively,they replaced the subsoil-monopile system by a set of springs through which the tower is connected to the subsoil,as shown in Fig. 2.This figure illustrates a mechanical model,in which the subsoil-monopile interaction is represented by four springs,i.e.a lateral,a rocking,a cross-coupling and a vertical spring whose stiffnesses areKL,KR,KLRandKV,respectively.Most researchers disregarded the axial vibrations since the wind turbines are very stiff vertically.

The stiffness of these springs which represent the subsoilmonopile interaction may be estimated from the monopile head load-deformation curves,provided that these curves are obtained by means of a rigorous modeling method,such as FEM.

On the basis of a numerical solution of transcendental frequency equation,Adhikari and Bhattacharya(2011,2012)proposed an exact approach where only lateral and rotational stiffnesses have been included.Furthermore,in order to improve the first natural frequency equation,Arany et al.(2014,2015)derived expressions of natural frequency of OWTs on three-spring flexible foundations by means of two beam models:Bernoulli-Euler and Timoshenko.The natural frequencies in both cases have been obtained numerically from the resulting transcendental equations.They proposed a closed-form expression containing,in addition toKLandKR,the cross-coupling stiffnessKLRof the monopile.Their equation for the natural frequency is

wherefFBis called fixed base frequency which can be either Eq.(1)or(2).The factorsCRandCLaccount for the stiffness provided by the monopile,and are functions of tower’s geometrical properties.Their analytical expressions are given by

whereaandbare the empirical constants,which have been obtained by fitting closed-form curves;EIηis the equivalent tower bending stiffness;andηis the soil-foundation interaction coefficient depending on tower’s bending stiffness.The applicability of Eqs.(4)and(5)is conditioned by

Fig. 2.The OWT model used:(a)Principal components and(b)Model considering soil stiffness through monopile head springs.

Although Eq.(3)is mathematically attractive as it contains three simple factors,finding its constituting parameters is not an easy task.This equation is very important in the sense that the two first parameters account for the interaction between the soil and the monopile through the spring stiffnesses.Moreover,they incorporate terms related to the equivalent tower stiffness which should be evaluated properly.

The OWTs,especially those installed in deep waters,differ from onshore wind turbines.This difference comes from the fact that these OWTs are considered as slender and heavy structures which require in general three elements in their construction to bear the heavy masses.These include the tower,the monopile overhang and the transition piece which assembles the first two elements(Fig. 2a).

In Fig. 3d,the length of the towerLTaccounts for the distance from the rotor nacelle assembly to the top of the transition piece.The tower has a varying bending stiffnessEIT,a thicknesstTand top and bottom diameters which are respectivelyDtandDb.Neglecting the f l exibilities of the grouted connection,the monopile overhang and transition piece welded together are assumed to constitute one element,called the substructure.The latter has a length ofLswhich is defined as the distance from the mudline(seabed)to the bottom of the tower.The diameterDsand the thicknesstsof the substructure are assumed to have the same values as those of the monopile on which the substructure is founded.Consequently,the bending stiffness of the substructure is the same as that of the monopileEIp.

As most towers are tapered and tubular,the value of the bending stiffness is used to evaluate the fixed base frequency and consequently the natural frequency is difficult to obtain,although the variation law ofEITwith the increasing tower cross-section is easy to establish.In this context,Bhattacharya(2011)studied a tower as a tapered cantilever beam subjected to a concentrated forcepapplied at its free end(Fig. 3a).Then,by means of beam theory,he computed a parameterfp(m)(termed here as‘tower stiffness coefficient’)as the ratio of the top displacement of a tower having a constant cross-sectionto that of a tower having a linearly varying cross-sectionThis parameter has been determined as

wheremis the ratio of bottom diameter to top diameter(Fig. 3d):

However,it is more likely to consider the tower as a tapered beam subjected to an upward tapered load along the whole length(Fig. 3b).Integrating the beam lateral deflection equation twice and setting the suitable boundary conditions,the tower stiffness coefficient corresponding to this load has been obtained as

Fig. 3.OWT model used to evaluate the tower mass and bending stiffness:(a)Tower subjected to p;(b)Tower subjected to downward tapered load;(c)Tower different masses;and(d)Geometrical properties of the tower.msis the mass of the substructure(transition piece+monopile)and Dpis the monopile diameter.

Fig. 4 illustrates the evolution offp(m)andfqtriang(m)with the increasing values ofm.It is clearly seen that for the interval where the ratiomvaries from 1 to 2.5,both curves yield the same values.This corresponds to the majority of practical applications.Outside this range,a neat discrepancy is observed and we suggest using the average value if it occurs to find a ratiomgreater than 2.5.

Fig. 4.Evolution of tower stiffness coefficient with tower diameters ratio.

The bending stiffnessEITof the tower may be evaluated as

The equivalent bending stiffness of the whole structure(support structure)is given by

Eqs.(1)and(2)for an onshore or offshore wind turbine,where the tower is resting directly on the soil,should be altered to properly find a fixed base natural frequency for an OWT having different masses and different bending stiffnesses,as shown in Fig. 3c.If weadopt Eq.(2),this would have the following form for an OWT composed of both tower and substructure parts:

where γsteelis the steel density usually taken as 7860 kg/m3.

Since the support structure composed of tower and substructure is in contact with soil,the empirical soil-foundation interaction coefficients given in Eq.(6)should be corrected in order to properly take into account the true support structure length and its effective bending stiffness.Thus these expressions are given by

The soil-foundation interaction coefficients given by Eqs.(4)and(5)are evaluated on the basis of Eq.(19).

2.2.Procedures to estimate the monopile head stiffnesses

In dealing with monopiles supporting OWTs,design engineers need to computeKL,KRandKLR.Two ways are often considered to compute these stiffnesses.The fi rst way is to model the monopile using the Winkler concept.This procedure which is alternatively calledp-yapproach is assumed to be sufficiently accurate for monopile diameterDp?2 m,asp-ycurves have been established for small-diameter and slender piles in offshore gas and oil industry.However,several investigations indicated that the pile deflections of large-diameter monopiles are underestimated for service loads and overestimated for small operational loads,which has been confirmed in a separate work(Otsmane and Amar Bouzid,2018).The second way is to directly employ values ofKL,KRandKLRgiven in the existing standards(EC8 for example,where pile-head stiffness of flexible pile is provided).Although the expressions containing these coefficients have been determined for various soil profiles(three profiles in most cases:constant soil stiffness,linear variation of soil stiffness with depth,and variation of soil stiffness with square root of depth),they encompass monopile-soil Young’s modulus ratioEp/Es.However,recent research(Higgins et al.,2013;Abed et al.,2016;Aissa et al.,2017)confirmed that the lateral behavior of large-diameter monopile fundamentally depends on monopile slendernessLp/Dprather than monopile-soil relative stiffnessEp/Es.

A different class of researchers suggested that the stiffness coefficients could be obtained from the elastic behavior of soilmonopile system under lateral loading where both soil and monopile are elastic.Indeed,some researchers(e.g.Carter and Kulhawy,1992;Higgins et al.,2013;Abed et al.,2016;Aissa et al.,2017)performed parametric studies using FEM,in which the ratiosEp/Eswere varied and load-deflection curves were drawn.These studies confirmed that the short monopile head stiffness for monopiles embedded in elastic media depends only on the monopile slenderness rather than the stiffness ratio.Results from the aforementioned references are given in Table 1 for homogeneous soils.The corresponding values for Gibson soil are given in Table 2,and Table 3 provides values for soils where stiffness varies with square root of depth.For simplicity,results presented in these tables are restricted to Poisson’s ratio equal to 0.4.

Using the values ofCRandCL,the natural frequency was found very close to unity,i.e.close to fixed base.Precisely,it may be argued that they do not reflect the soil actual stiffness and they do not bring an actual alteration to the natural frequency which remains very close to the fixed base frequency.

As an alternative procedure,it is appropriate to evaluate the stiffness coefficients on the basis of initial stiffness(tangential values at the origin)of the monopile head-deformations curves,resulting from the study of the nonlinear behavior of the subsoil in which the monopile is embedded.Assuming that the monopile head movements and applied efforts are expressed in function of fl exibility coefficients,this can be given in matrix form as whereHandMare the shear force and overturning moment applied at the monopile head,respectively;uLandθRare the lateral displacement and rotation of the monopile head,respectively;andIL,IR,andILRare the lateral,rotational and cross-coupling flexibility coefficients,respectively.

As the aim is to find the stiffness coefficients,it is easy to reverse the matrix in Eq.(20)to obtain:

The stiffness coefficients are related to flexibility ones by the following terms:

Table 1Short monopile stiffness coefficients proposed by different researchers in homogeneous soils.

Table 2Stiffness coefficients for short monopiles proposed for Gibson soils.

Table 3Stiffness coefficients for short monopiles proposed for soils whose stiffness increases with square root of depth.

Determination of the valuesKL,KRandKLRis not straightforward in the FEM analyses controlled by forces.The flexibility coefficients should be determined first,and then be reversed to obtain the stiffness coefficients of Eq.(22).To do so,an arbitrary pure horizontal load(H≠0 andM=0)is applied at the monopile head at the mudline level,and a plot depicting the increasing values ofHversus the corresponding values of monopile head displacementsuis illustrated(Fig. 5a).The parameter 1/ILis then obtained bysimply computing the slope of the resulting curve at the origin.The parameter 1/ILRis computed from the curve giving the variation ofHin function of rotationθissued from the same analysis(Fig. 5c).

As the rocking flexibility coefficient needs a pure bending,the monopile-soil system is analyzed under an overturning moment(M≠0 andH=0)applied at the topof the pile at the mudline level.From the curve portraying theMincreasing values against the obtained rotationsθ,the reciprocal of flexibility coefficient 1/IRis evaluated by simply computing the slope of curve tangent at the origin(Fig. 5b).

This procedure is followed in this paper,when OWT monopiles of the different wind farms are considered in the next sections.

3.Numerical methodology:the computer program NAMPULAL

A pseudo 3D FEM model has been performed to study soilstructure interaction problems in nonlinear media.This procedure,called nonlinear finite element vertical slices model(NFEVSM),involves the combination of the FEM and the FDM for capturing the behavior of the embedded structure and its surrounding soil being considered to obey the hyperbolic model as proposed by Duncan and Chang(1970).The 3D soil-structure problem plotted in Fig. 6 shows a soil-structure interaction problem example(Fig. 6a)and the vertical slices model where different slices are acted upon by external forces and body forces(Fig. 6b).

The stress and deformation analyses in each slice are conducted by the conventional FEM,using two-dimensional(2D)finite elements.According to the standard formulation in the displacement based FEM,the element stiffness matrix in sliceican be written as

where B and BTare the strain field-nodal displacement matrix and its transpose,respectively;N and NTare the shape function matrix and its transpose,respectively;piis the external force vector to which the sliceiis subjected;and biis the body force vector which has the following compact form:

where

Fig. 5.Monopile head load-movement curves permitting to obtain(a)spring lateral stiffness,(b)spring rocking stiffness,and(c)spring coupling stiffness.

Fig. 6.(a)Real-world soil-structure interaction problem,and(b)The vertical slices model showing the interacting slices subjected to external and body forces.

where the superscripts pc,pr and f lst and for proper contribution of the slice itself,contribution of the preceding slice,and contribution of the subsequent slice,respectively;ai-1,aiand ai+1are the element nodal displacement vectors of slicesi-1,iandi+1,respectively;I is the identity matrix;andGi-1,GiandGi+1are the shear moduli at slicesi-1,iandi+1,respectively.

In Eq.(22),the matrix Dpscorresponds to a problem of plane stresses,which may be given as

where νsis the Poisson’s ratio,andEis the Young’s modulus.

From this fact,and unlike a fully 3D or plane strain problem,the value of Poisson’s ratio equal to 0.5 is no longera singular value.It is clear from Eq.(24)that the fictitious body forces applied to a sliceidepend essentially on its own nodal displacements and on those of slices sandwiching it.The numerical analysis of the vertical slicing model has led to the familiar equations of a pseudo plane stress problem with body forces representing the interaction between the slices,forming the structure and its surrounding medium.Substituting Eqs.(24)and(25)into Eq.(23)gives a more detailed governing equation:

In a more compact form,Eq.(29)becomes

This equation cannot be solved straight-fully,since the right hand terms are not available explicitly at the same time.Consequently,an updating iterative process is needed:

wherejdenotes the iteration number andjmaxis the maximum number of iterations allowed in the numerical process.

The nonlinearity in vertical slices model stems from the implementation of the hyperbolic model proposed by Duncan and Chang(1970)for modeling the soil.In fact,they found out that both tangential modulusEiand ultimate stress difference(σ1- σ3)ultare dependent on the minor principal stressσ3.More precisely,they suggested for the initial tangent modulus the following formula:

whereKis the dimensionless factor termed as ‘modulus number’,nis a dimensionless parameter called ‘modulus exponent’,andpais the atmospheric pressure used to makeKandndimensionless.

The ultimate stress difference(σ1- σ3)ultis defined in terms of the actual failure stress difference by another parameter called‘failure ratio’Rfwhich is given by

Using Mohr-Coulomb failure criterion where the envelope is considered as a straight line,the principal stress difference at failure is related to the confining pressureσ3as

wherecis the cohesion andφis the internal friction angle.The tangent modulusEtis given by

For unloading and reloading cycles,Duncan and Chang(1970)proposed the following expression:

whereEuris the unloading-reloading modulus andKuris the corresponding modulus number.

A thorough literature investigation has been performed by Otsmane and Amar Bouzid(2018)to keep the hyperbolic modeling parameters sufficiently accurate and to make their use practical for solving the soil-structure interaction problems.Theauthorsexamined a large number of well-established correlations between soil physical parameters especially those of sandy deposits whose behaviors are mainly governed by their internal friction angles,and proposed relationships between the sand relative density and the confining pressure.This has been achieved using mainly the recommendations made by well-known researchers who carried out a great number of careful experiments.These parameters are listed in Tables 4 and 5 along with references of their origin.

Table 4Soil stiffness parameters in terms of soil friction angle and confining pressure.

Table 5Parameters governing the hyperbolic model according to correlations and recommendations.

Equations in both Tables 4 and 5 have been implemented in the FEM computer code NAMPULAL which will be described in the next paragraph for evaluating soil model parameters related to the five wind farm sites considered in Section 4.In the Duncan-Chang’s basic model,the Poisson’s ratio νswas assumed to be constant throughout the whole process.

A Fortran computer program called NAMPULAL for the analysis of axially and laterally loaded single monopiles has been written.Although approximate,the computer code NAMPULALisa coherent tool which exhibits many advantages over other numerical codes in dealing with nonlinear soil-structure interaction problems.For further details,the reader can refer to Otsmane and Amar Bouzid(2018)and only the features of this computer program are given here.

Although the computational process in NAMPULAL is naturally iterative to fulfill slices equilibrium,it does not require a significant number of iterations to reach convergence.For the problems analyzed so far,a number of 20 iterations are generally sufficient to reach accurate solutions within acceptable margins.

A number of 20 slices have been implemented in NAMPULAL.This number,which has been set on the basis of parametric study involving many monopile behavior parameters(Amar Bouzid et al.,2005),has been found sufficient to accurately model many soilstructure interaction problems(Amar Bouzid et al.,2005;Otsman and Amar Bouzid,2018).

Unlike most implemented elastoplastic constitutive models,which necessitate a significant number of iterations to subdue the unbalanced forces,the implemented hyperbolic model in NAMPULAL requires only two iterations.This fact alleviates considerably the whole process of solution and makes it easier to find fast solutions even for the most complex soil-structure interaction problems.

In addition to the rectangular cross-sectional monopiles that are automatically considered due to the shape of the vertical slice,the solid circular or tubular cross-sectional monopiles are easily dealt with by prescribing the effective bending stiffness.Hence,an equivalent Young’s modulus is adopted according to the following formula:

This equation has been set on the assumption that the square cross-sectional monopile under consideration in NAMPULAL has the same cross-sectional area as the effective circular crosssectional monopile.

The performances of this computer code have been assessed against analysis of the behavior of a number of OWT monopiles where other commercial packages are used,such as FLAC3D,ABAQUS?and PLAXIS(Otsmane and Amar Bouzid,2018).The results were in excellent agreement with those of the aforementioned powerful numerical tools.This computer code is employed to determine the monopile head stiffnesses for the OWTs examined in this paper.

Table 6List of the five OWTs with soil conditions at the sites.

Table 7Input parameters for the five OWTs chosen for this study.

Table 8Masses and bending stiffnesses of the OWTs constitutive elements.

Table 9Adopted soil deformation and strength parameters as well as hyperbolic model parameters for the OWTs chosen.

Fig. 7.Monopile head displacement against applied horizontal load for different OWTs.

Fig. 8.Monopile head rotation against applied horizontal load for different OWTs.

Fig. 9.Monopile head rotation against applied overturning moment for different OWTs.

4.Computed and measured first natural frequencies for five different offshore wind turbines

In order to assess the performances of the computer code NAMPULAL for a wide range of geotechnical applications,in terms of the accuracy,utilityand potential of the NFEVSM,five OWTs have been selected from five wind farm sites.These are Lely A2(UK),Irene Vorrink(Netherlands),Kentish Flats(UK),Walney 1(UK)and North Hoyle(UK).These wind turbines have been chosen for the full availability of their data,especially the measured first natural frequency.Soil conditions at the site and sources from which the OWT data are adopted are summarized in Table 6.

OWTs structural data along with some soil deformation characteristics and measured natural frequencies are presented in Table 7.Since these data come directly from the OWT manufacturers,it is not possible to check their accuracy,except data relevant to the tower mass,provided that the tower height is correct.

Slight differences between computed values ofmTand those provided in the reference(Arany et al.,2016)are noticed.Thus computed data in Table 8 are used in the coming computations.

Although the OWT structural data were available which enable the users to compute any structural behavior parameter,the parameters relevant to soil behavior were not found.However,only the different strata of each site are given,but nothing about strength and deformation parameters.

The site investigations indicate that almost all the OWTs chosen in this paper are installed through deep layers of dense sand.The pertinent hyperbolic parameters have been computed according to the prescribed values and relationships given in Tables 4 and 5.These are given in Table 9.

A comprehensive mesh study has been performed to find the optimal finite element mesh that captures the behaviors of monopiles under lateral loading in a nonlinear medium characterized by the hyperbolic model as a yield criterion.A mesh of 20 times monopile diameterDpin both sides of the monopile and one monopile lengthLpunder the monopile tip has been adopted for the study of all OWTs considered here.Furthermore,35 finite elements in both sides of the monopile and 36 finite elements in vertical direction as well as 20 slices have been chosen to analyze the pseudo 3D medium under consideration.

As the monopile head stiffness does not depend on the loading level,a horizontal loadHof 1000 kN is applied in 10 increments at the top of each monopile in the five wind farms considered,aiming to compute the monopile head flexibility coefficientsILandILR.Fig. 7 shows the evolution of monopile head displacements as a function of the increasing lateral loadH.

The evolution of rotations is a function of applied lateral load and is plotted in Fig. 8.This figure is used to determine the cross coupling flexibility coefficientILRfor all monopiles considered here.As the flexibility coefficientIRrequires a pure bending,an applied momentMat the top of monopile of 20,000 kN m in magnitude has been considered and the corresponding rotations are plotted in Fig. 9.

Table 10Flexibility coefficients IL,ILRand IRand their corresponding stiffness coefficients KL,KLRand KRrelevant to monopiles in the OWTs chosen.

Table 11Fixed base natural frequencies for different OWTs.

The almost linear relationships betweenHanduandHandθon one hand and that betweenMandθon the other hand somewhat make it easy to computeIL,ILRandIRwhich can be performed by simply inverting the slopes of their corresponding load deformation curves.Then using Eq.(22),the stiffness coefficients are obtained.Flexibility and stiffness coefficients are respectively shown in Tables 10 and 11 for all turbines considered in this paper.

Fig. 10.Values of KL,KLRand KRgiven by NAMPULAL against those developed by Arany et al.(2016).

Table 12

CRandCLcomputed values for the OWTs considered in the current study.

Wind farm nameCRCLLely A2 0.867 0.996 Irene Vorrink 0.867 0.997 Kentish Flats 0.857 0.998 Walney 1 0.837 0.997 North Hoyle 0.849 0.998

Table 13

Predicted and measured natural frequencies of all OWTs.Wind farm namePredicted frequency(fη=CRCLfFB)(Hz)

Measured frequency(Hz)Error(%)

Lely A2 0.621 0.634 2.05

Irene Vorrink 0.570/0.579 0.546-0.563 4.39/2.84 Kentish Flats 0.315 0.339 7.08

Walney 1 0.277 0.35 20.86

North Hoyle 0.342 0.35 2.28

As the monopile head stiffness coefficients play an important role in the correct assessment of the natural frequency,which is in turn a significant parameter in the design of any OWT,it is useful to compare the values listed in Table 10 with those provided by other methods.On the basis of the formulas developed by Poulos and Davis(1980),Randolph(1981)and Carter and Kulhawy(1992).Arany et al.(2016)determined the values of the monopile stiffness coefficients which are added to the histograms of Fig. 10 for comparison.

One important point can emerge from the close examination of the histograms shown in Fig. 10.NAMPULAL’s results are approximately half those given by Arany et al.(2016)for the OWTs whose supporting monopiles are driven in dense sand.We believe that the NFEVSM results are more accurate than those of Arany et al.(2016).This is probably due to the fact that these authors used data from works performed on slender piles using the Winkler model for which many questions had been raised about its applicability to large-diameter monopiles.

Eq.(15),whose different constitutive parts are evaluated using Eqs.(16)and(17),is employed here to give the fixed base natural frequency.This expression,which depends only on the OWT structure properties,gives the values of the fixed base natural frequencies for different turbines shown in Table 11.

The factorsCRandCLdepend on values ofIL,IRandILR.Table 12 shows the values ofCRandCLfor the five OWTs considered in this paper.These values make it quite clear thatCRis the dominant factor that can bring the value of the fixed base frequency to the measured one.However,CLis very close to unity,and hence its influence in changing the value offFBis very small.This has been also noticed by Arany et al.(2016).

The natural frequency which is simply obtained by multiplying the flexibility coefficients by the fixed base frequency for each OWT is given in Table 13.Also shown are errors between the measured and the computed natural frequencies.5.Conclusions

In this article,a nonlinear finite element computer code NAMPULAL developed for soil-pile interaction has been used to analyze five different monopiles from five European wind farms.

As the natural frequency of the whole wind turbine structure is of paramount importance in the design of OWTs,developing reliable methods for its determination is an active area of research.Stiffness of foundation is the key of natural frequency calculation and this work is based on three-spring model where the monopilesoil interaction is modeled by lateral springKL,cross-couplingKLRand rotational springKR.NAMPULAL code has been adopted to find the foundation stiffness.The code has the capability to incorporate nonlinear soil model and in this study,Duncan-Chang hyperbolic model has been used.These stiffnesses in turn were used to obtain natural frequency of the whole wind turbine system and the results obtained were compared with the measurements.Good agreement was noted between prediction and observation.Conflict of interest

The authors wish to confirm that there are no known Conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Abed Y,Amar Bouzid Dj,Bhattacharya S,Aissa MH.Static impedance functions for monopiles supporting offshore wind turbines in non-homogeneous soils emphasis on soil/monopile interface characteristics.Earthquakes and Structures 2016;10(5):1143-79.

Achmus M,Kuo YS,Abdel-Rahman K.Behavior of monopile foundations under cyclic lateral load.Computers and Geotechnics 2009;36(5):725-35.

Adhikari S,Bhattacharya S.Dynamic analysis of wind turbine towers on flexible foundations.Shock and Vibrations 2012;19(1):37-56.

Adhikari S,Bhattacharya S.Vibrations of wind-turbines considering soil-structure interaction.Wind and Structures 2011;14(2):85-112.

Aissa MH,Amar Bouzid Dj,Bhattacharya S.Monopile head stiffness for serviceability limit state calculations in assessing the natural frequency of offshore wind turbines.International Journal of Geotechnical Engineering 2017;35(7).https://doi.org/10.1080/19386362.2016.1270794.

Amar Bouzid Dj.Monopile head stiffness and natural frequency assessment of some installed OWTs using a pseudo 3D nonlinear FE model.In:World congress on Advances in Civil,Environmental and Materials research(ACEM16).Guildford,Surrey,UK:The University of Surrey;2016.p.1-15.http://epubs.surrey.ac.uk/id/eprint/844988.

Amar Bouzid Dj,Vermeer PA,Tiliouine B.Finite element vertical slices model:validation and application to an embedded square footing under combined loading.Computers and Geotechnics 2005;32(2):72-91.

American Petroleum Institute(API),International Organization for Standardization(ISO).ANSI/API specification RP 2GEO:geotechnical considerations and foundation design for offshore structures.Washington,D.C.,USA:API;2011.

Arany L,Bhattachary S,Hogan SJ,Macdonald J.Dynamic soil-structure interaction issues of offshore wind turbines.In:Proceedings of the 9th international conference on structural dynamics,Porto,Portugal;2014.p.3611-7.

Arany L,Bhattacharya S,Adhikari S,Hogan SJ,Mcdonald JHG.An analytical model to predict the natural frequency of offshore wind turbines on three-spring flexible foundations using two different beam models.Soil Dynamics and Foundation Engineering 2015;74:40-5.

Arany L,Bhattacharya S,Adhikari S,Mcdonald JHG,Hogan SJ.Closed form solution of Eigen frequency of monopile supported offshore wind turbines in deeper waters stiffness of substructure and SSI.Soil Dynamics and Foundation Engineering 2016;83:18-32.

Bhattacharya S.SDOWT:simplif i ed dynamics of wind turbines.User’s Manual.Bristol,UK:Bristol Laboratory for Advanced Dynamics Engineering;2011.

Blevins RD.Formulas for frequencies and mode shapes.Krieger Pub Co;2001.

Carter JP,Kulhawy FH.Analysis of laterally loaded shafts in rock.Journal of Geotechnical Engineering 1992;118(6):839-55.

Damgaard M,Bayat M,Andersen LV,Ibsen LB.Dynamic response sensitivity of an offshore wind turbine for varying subsoil conditions.Ocean Engineering 2015;101:227-34.

DNV-OS-J101.Offshore standard:design of offshore wind turbine structures.Hellerup,Denmark:Det Norske Veritas(DNV);2004.

Duncan JM,Chang CY.Nonlinear analysis of stress and strain in soils.Journal of Soil Mechanics and Foundation Division 1970;96(5):1629-53.

Duncan JM,Wong KS.Soil properties manual.User’s manual for SAGE,vol.II.Blacksburg,Virginia,USA:Centre for Geotechnical Practice and Research,Virginia Polytechnic Institute and State University;1999.

Higgins W,Vasquez C,Basu D,Griffiths DV.Elastic solutions for laterally loaded piles.Journal for Geotechnical and Geoenvironmental Engineering 2013;139(7):1096-103.

Jaky JA.Nyugalmi nyomas tényez?je(The coefficient of earth pressure at rest).Agyar Mérnok és Epitész Egylet K?zl?nye(Journal for Society of Hungarian Architects and Engineers)1944:355-8(in Hungarian).

Leblanc C.Design of offshore wind turbine support structures-selected topics in the field of geotechnical engineering.PhD Thesis.Aalborg,Denmark:Aalborg University;2009.

Lombardi D,Bhattacharya S,Wood DM.Dynamic soil-structure interaction of monopile supported wind turbines in cohesive soil.Soil Dynamics and Earthquake Engineering 2013;49:165-80.

Otsmane L,Amar Bouzid D.An efficient FE model for SSI:theoretical background and assessment by predicting the response of large diameter monopiles supporting OWECs.Computers and Geotechnics 2018;97:155-66.

Poulos HG,Davis EH.Pile foundation analysis and design.New York,USA:Wiley;1980.

Prendergast LJ,Gavin K,Doherty P.An investigation into the effect of scour on the natural frequency of an offshore wind turbine.Ocean Engineering 2015;101:1-11.

Randolph MF.The response of flexible piles to lateral loading.Géotechnique 1981;31(2):247-59.

Reese L,Cox WR,Koop FD.Analysis of laterally loaded piles in sand.In:Proceedings of the offshore technology conference,Houston,USA;1974.https://doi.org/10.4043/2080-MS.

Vught JH.Considerations on the dynamics of support structures for an offshore wind energy converter.PhD Thesis.Delft,the Netherlands:Delft University of Technology;2000.

Wong KS,Duncan JM.Hyperbolic stress-strain parameters for non-linear finite element analyses of stresses and movements in soil masses.Report No.TE-74-3.Berkeley,USA:University of California;1974.

Yi JH,Kim SB,Yoon GL,Andersen LV.Natural frequency of bottom-fixed offshore wind turbines considering pile-soil-interaction with material uncertainties and scouring depth.Wind and Structures 2015;21(6):625-39.

Yu L,Bhattacharya S,Li L,Guo Z.Dynamic characteristics of offshore wind turbines on different types of foundations.Electronic Journal of Geotechnical Engineering 2014;19:2917-36.

Zaaijer MB.Design methods for offshore wind turbines at exposed sites(OWTES)-sensitivity analysis for foundations of offshore wind turbines.Technical report.Delft,The Netherlands:Delft University of Technology;2002.

Zaaijer MB.Foundation modeling to assess dynamic behavior of offshore wind turbines.Applied Ocean Research 2006;28(1):45-57.