Philip Alkhoury
email: philip.alkhoury@etu.univ-nantes.fr

Abdul-Hamid Soubra

Valentine Rey

Mourad Aït-Ahmed

Dynamic analysis of a monopile-supported offshore wind turbine

Keywords: Offshore wind turbine, Dynamic analysis, Monopile, Finite element analysis, Superstructure, Constitutive models

come

Introduction

In recent years, offshore wind energy received a vast global attention as being one of the most promising renewable energy resources. Nowadays, multi-megawatt offshore wind turbines (OWTs) are widely adopted in practice to efficiently extract the steady offshore wind energy.

In order to reduce their cost, modern multi-megawatt OWTs are designed to have a minimum overall weight (tall slender tower) and a powerful generator atop of the tower (large rotor). This makes the OWT structure dynamically-sensitive to the external loadings (particularly the wind and sea wave loads). The simultaneous action of these loads may lead to excessive vibrations of the wind turbine structure, thus causing fatigue damage of its structural components. Therefore, to ensure an efficient functioning of these multi-megawatt OWTs, it is vital to accurately study their dynamic behavior in the presence of realistic environmental loadings.

The foundation of an OWT is an important part of the whole system as it has to withstand the loads coming from the structure and remain stable during the whole lifetime of the OWT.

Monopile foundations remain the preferred choice of developers with over two-thirds of installation in Europe in 2020 [START_REF]The European offshore wind industry[END_REF]. A typical monopile is a long steel cylindrical tube with a diameter ranging from 4 to 6 m; however, due to increasing turbine sizes and water depths, large diameter monopiles of up to 10 m in diameter are being anticipated for future OWT foundations [START_REF] Murphy | 3D FEM approach for laterally loaded monopile design[END_REF]. For such a foundation supporting a dynamically-sensitive structure like the OWT, the interactions between the superstructure, the monopile and the surrounding soil will significantly influence the dynamic response of the OWT and shall be considered during the design stage [START_REF] Mostafa | Response of fixed offshore platforms to wave and current loading including soil-structure interaction[END_REF]. Further, the well prediction of the foundation response requires a soil model which must be sufficiently accurate to reproduce the soil behavior under complex and irregular dynamic loading over a wide range of strain levels [START_REF] Jostad | Evaluation of soil models for improved design of offshore wind turbine foundations in dense sand[END_REF].

In literature, several studies were devoted to the dynamic analysis of wind turbines under various environmental loadings. To simplify the analysis, the wind turbines blades were usually modelled by a lumped mass located at the tower top and thus, the blades geometrical configuration and their interaction with the tower (via their rotational velocity) were completely neglected [START_REF] Bisoi | Dynamic analysis of offshore wind turbine in clay considering soil-monopile-tower interaction[END_REF][START_REF] Haldar | Probabilistic analysis of monopile-supported offshore wind turbine in clay[END_REF][START_REF] Corciulo | Transient response of offshore wind turbines on monopiles in sand: role of cyclic hydro-mechanical soil behaviour[END_REF][START_REF] Zuo | Using multiple tuned mass dampers to control offshore wind turbine vibrations under multiple hazards[END_REF][START_REF] Lombardi | Dynamic soil-structure interaction of monopile supported wind turbines in cohesive soil[END_REF][START_REF] Kementzetzidis | Geotechnical aspects of offshore wind turbine dynamics from 3D non-linear soil-structure simulations[END_REF]. The simplifications on the geometrical characteristics of the blades as well as their interaction with the turbine tower can directly influence the vibration characteristics of the wind turbine, which can affect its dynamic response [START_REF] Zuo | Dynamic analyses of operating offshore wind turbines including soil-structure interaction[END_REF][START_REF] Ali | Seismic assessment of wind turbines: How crucial is rotor-nacelle-assembly numerical modeling?[END_REF]. To study the influence of the blade modelling on the dynamic response of wind turbines, [START_REF] Zuo | Dynamic analyses of operating offshore wind turbines including soil-structure interaction[END_REF] explicitly considered the blades and developed a three-dimensional (3D) finite element model of the NREL 5 MW wind turbine within Abaqus where the tower and blades were modelled using shell elements and the nacelle/hub was represented by a lumped mass located at the tower top.

Their findings have shown that the explicit modeling of the blades, the operational conditions and the rotor velocity have a significant influence on the dynamic behavior of the OWT. Also [START_REF] Ali | Seismic assessment of wind turbines: How crucial is rotor-nacelle-assembly numerical modeling?[END_REF] recently examined the effects of three key rotor-nacelle assembly (RNA) parameters: (i) rotary inertia of the blades, (ii) rotor eccentricity and (iii) the blades flexibility, on the nonlinear dynamic response of OWTs. Their findings have shown that the conventional simplified lumped mass approach of the RNA can lead to a significant design issue for OWTs.

Several other authors have investigated the influence of the soil-structure interaction (SSI) on the vibration characteristics [START_REF] Zuo | Dynamic analyses of operating offshore wind turbines including soil-structure interaction[END_REF][START_REF] Andersen | Natural frequencies of wind turbines on monopile foundations in clayey soils-a probabilistic approach[END_REF][START_REF] Arany | Closed form solution of eigen frequency of monopile supported offshore wind turbines in deeper waters incorporating stiffness of substructure and SSI[END_REF] and dynamic response of wind turbines [START_REF] Bisoi | Dynamic analysis of offshore wind turbine in clay considering soil-monopile-tower interaction[END_REF][START_REF] Haldar | Probabilistic analysis of monopile-supported offshore wind turbine in clay[END_REF][START_REF] Corciulo | Transient response of offshore wind turbines on monopiles in sand: role of cyclic hydro-mechanical soil behaviour[END_REF][START_REF] Zuo | Using multiple tuned mass dampers to control offshore wind turbine vibrations under multiple hazards[END_REF][START_REF] Kementzetzidis | Geotechnical aspects of offshore wind turbine dynamics from 3D non-linear soil-structure simulations[END_REF][START_REF] Zuo | Dynamic analyses of operating offshore wind turbines including soil-structure interaction[END_REF][START_REF] Ali | Seismic assessment of wind turbines: How crucial is rotor-nacelle-assembly numerical modeling?[END_REF].

The SSI was usually considered by using simplified foundations models (i.e. API p-y curves, coupled/uncoupled springs at mudline) [START_REF] Bisoi | Dynamic analysis of offshore wind turbine in clay considering soil-monopile-tower interaction[END_REF][START_REF] Haldar | Probabilistic analysis of monopile-supported offshore wind turbine in clay[END_REF][START_REF] Zuo | Using multiple tuned mass dampers to control offshore wind turbine vibrations under multiple hazards[END_REF][START_REF] Ali | Seismic assessment of wind turbines: How crucial is rotor-nacelle-assembly numerical modeling?[END_REF][START_REF] Andersen | Natural frequencies of wind turbines on monopile foundations in clayey soils-a probabilistic approach[END_REF][START_REF] Arany | Closed form solution of eigen frequency of monopile supported offshore wind turbines in deeper waters incorporating stiffness of substructure and SSI[END_REF]. The major drawback of the simplified foundation models is related to the determination of suitable dynamic stiffness and damping characteristics of the foundation as they proved to have a considerable impact on the global OWT dynamics [START_REF] Page | Impact of foundation modelling in offshore wind turbines: comparison between simulations and field data[END_REF]. Finally notice that existing design guidelines [START_REF] Germanische | Guideline for the certification of offshore wind turbines[END_REF][START_REF]DNV-OS-J101 offshore standard: design of offshore wind turbine structures[END_REF] recommend a profound understanding of the dynamic soil-monopile interaction under cyclic/dynamic loading conditions.

In light of the above-mentioned premises, several experimental studies [START_REF] Leblanc | Response of stiff piles in sand to long-term cyclic lateral loading[END_REF][START_REF] Cuéllar | On the quasi-static granular convective flow and sand densification around pile foundations under cyclic lateral loading[END_REF][START_REF] Bayton | Centrifuge modelling of long term cyclic lateral loading on monopiles[END_REF] and 3D finite element analyses [START_REF] Corciulo | Transient response of offshore wind turbines on monopiles in sand: role of cyclic hydro-mechanical soil behaviour[END_REF][START_REF] Kementzetzidis | Geotechnical aspects of offshore wind turbine dynamics from 3D non-linear soil-structure simulations[END_REF][START_REF] Kementzetzidis | 3D FE dynamic modelling of offshore wind turbines in sand: natural frequency evolution in the pre-to after-storm transition[END_REF][START_REF] Taborda | Finite-element modelling of laterally loaded piles in a dense marine sand at Dunkirk[END_REF][START_REF] Achmus | Behavior of monopile foundations under cyclic lateral load[END_REF] were conducted to investigate the behavior of monopiles in sand under cyclic/dynamic loading assuming both undrained and fully drained conditions.

The obtained results deemed to strongly depend on the constitutive model used to describe the sand behavior. In this regard, [START_REF] Corciulo | Transient response of offshore wind turbines on monopiles in sand: role of cyclic hydro-mechanical soil behaviour[END_REF] developed a 3D hydro-mechanical analysis to study the dynamic response of a monopile-supported OWT in sand. In their study, a multi-surface plasticity model UCSD08 by [START_REF] Yang | Multi-surface cyclic plasticity sand model with lode angle effect[END_REF] was used to describe the dynamic sand behavior. [START_REF] Kementzetzidis | Geotechnical aspects of offshore wind turbine dynamics from 3D non-linear soil-structure simulations[END_REF][START_REF] Kementzetzidis | 3D FE dynamic modelling of offshore wind turbines in sand: natural frequency evolution in the pre-to after-storm transition[END_REF] have also employed a hydro-mechanical analysis as [START_REF] Corciulo | Transient response of offshore wind turbines on monopiles in sand: role of cyclic hydro-mechanical soil behaviour[END_REF] but adopted for the sand the criticalstate SANISAND model [START_REF] Dafalias | Simple plasticity sand model accounting for fabric change effects[END_REF] to capture the response of an 8 MW OWT under two different loading histories of 10 min [START_REF] Kementzetzidis | Geotechnical aspects of offshore wind turbine dynamics from 3D non-linear soil-structure simulations[END_REF] and 2 h [START_REF] Kementzetzidis | 3D FE dynamic modelling of offshore wind turbines in sand: natural frequency evolution in the pre-to after-storm transition[END_REF]. Although these studies [START_REF] Corciulo | Transient response of offshore wind turbines on monopiles in sand: role of cyclic hydro-mechanical soil behaviour[END_REF][START_REF] Kementzetzidis | Geotechnical aspects of offshore wind turbine dynamics from 3D non-linear soil-structure simulations[END_REF][START_REF] Kementzetzidis | 3D FE dynamic modelling of offshore wind turbines in sand: natural frequency evolution in the pre-to after-storm transition[END_REF] incorporated an advanced soil constitutive model, the superstructure was not explicitly considered. Indeed, the tower was modelled by a beam element and the RNA was replaced by a lumped mass placed at the tower top.

From the previous literature, it seems that the dynamic response of an OWT strongly depends on several factors. Thus, it is essential to have a reliable numerical calculation tool which simultaneously considers the key features of the OWT: (i) an explicit modelling of the superstructure and the foundation system taking into account the different interactions between the soil, the foundation and the superstructure and (ii) a suitable modelling of the soil behavior under complex dynamic loading. For this purpose, a detailed 3D finite element analysis of the entire soil-foundation-superstructure system, considering an advanced constitutive model for the soil and a large deformation analysis for the superstructure and the monopile foundation, was performed in this paper to capture the nonlinear dynamic response of an OWT subjected to realistic environmental loading scenarios.

The aim of this paper is to study the dynamic behavior of a large diameter monopilesupported multi-megawatt 10 MW DTU [START_REF] Bak | Description of the DTU 10 MW reference wind turbine[END_REF] OWT installed in sand under perfectly drained conditions using a detailed 3D finite element (FE) model within Abaqus/Standard. The 3D mechanical model simultaneously considers (i) the real geometrical configuration of the OWT superstructure (tower, blades, transition piece, hub and nacelle) with the monopile foundation and, (ii) the monopile-sand interaction. The present mechanical model follows the recent work of the authors [START_REF] Alkhoury | A full three-dimensional model for the estimation of the natural frequencies of an offshore wind turbine in sand[END_REF] but makes use herein of sophisticated constitutive models for the soil and a large deformation analysis for the superstructure and the monopile foundation.

Nonlinear 10-minutes time-domain simulations are performed under stochastically simulated aerodynamic and hydrodynamic loadings, with emphasis on the OWT structural and geotechnical responses. Two advanced soil constitutive models which proved to be capable of capturing the dynamic behavior of granular soils are used in this work: (i) the critical state elastoplastic SANISAND constitutive model and (ii) the extended hypoplastic model with intergranular strain. The recent calibration done on the Karlsruhe fine sand by [START_REF] Wichtmann | Inspection of three sophisticated constitutive models based on monotonic and cyclic tests on fine sand: Hypoplasticity vs. Sanisand vs. ISA[END_REF] for both constitutive models is used herein to compare the structural and geotechnical OWT dynamic responses as predicted using the two calibrated constitutive models. Also, the relevance of both constitutive models was tested by comparing their results with those obtained using the conventional elastic-perfectly plastic Mohr-Coulomb model. Finally, the accuracy of two existing simplified foundation models (based on the distributed p-y soil springs) as suggested in literature in the case of dynamic/cyclic loading is examined and discussed.

The paper is organized as follows: a description of the 3D mechanical model of the DTU 10 MW OWT is presented in section 2. Section 3 presents the soil constitutive models used in this paper. The stochastically simulated aerodynamic and hydrodynamic loads acting on the OWT structure are described in section 4. Section 5 gives the vibration characteristics and the damping sources of the OWT. Numerical simulations related to the nonlinear dynamic analysis are presented in Section 6. Section 7 discusses the numerical results. Finally, some conclusions are presented in section 8.

Modeling of the soil-monopile-superstructure system

As mentioned before, the OWT system chosen in this study is the reference DTU 10 MW three-bladed OWT defined in [START_REF] Bak | Description of the DTU 10 MW reference wind turbine[END_REF]. It is representative of current industry practice of OWTs being manufactured today. The turbine is characterized by a tower height of 115.63 m and a decreasing diameter and thickness profile from bottom (8.3 m, 0.038 m) to top (5.5 m, 0.02 m). A monopile foundation of 8.3 m diameter (and thickness of 9 cm) and of a length of 80 m is chosen to support the wind turbine. The monopile length covers three parts (45 m embedded in the seabed, 25 m in contact with the sea water and 10 m above the MSL to represent the transition piece). Notice that the monopile thickness is chosen to respect the API recommendation [START_REF]Recommended practice for planning, designing and constructing fixed offshore platforms: Working stress design[END_REF]. A summary of the relevant dimensions and properties of the offshore wind turbine are provided in Figure 1 and Table 1.

Figure 2 shows the whole 3D model of the OWT developed within the FE Abaqus software (version 6.14 Standard) [START_REF] Abaqus | Finite element analysis program[END_REF]. As shown in this figure, the developed FE model consists of seven main parts: (i) 3D soil domain, (ii) monopile, (iii) transition piece, (iv) turbine tower, (v) hub/nacelle assembly, (vi) hinge connector and (vii) blades. In the 3D model, shell elements (S4 in Abaqus) were used to discretize the steel structure above the MSL, while solid elements were used to discretize the steel monopile. The steel was considered as a linear elastic isotropic material. Structural and equipment masses as well as the added mass effect (due to the surrounding water) were included in the 3D model by increasing the effective density of the steel. Concerning the nacelle/hub assembly, a reference point (RP) eccentric to the tower top and which coincides with the nacelle center of mass position was used to model the nacelle/hub assembly as a lumped mass (point mass in Abaqus). The mass and the rotary inertia of the nacelle/hub assembly were defined at this RP.

To accurately take into account the influence of the blades stiffness and geometry on the OWT dynamic response, each blade was divided along its length into 51 segments. A generalized beam cross-section was defined for every segment of the partitioned blade and for each cross-section, its corresponding stiffness and mass properties were assigned. Finally, a hinge connector was used to simulate the rotation of the blades with respect to the tower.

The dynamic response of the OWT under the simultaneous effect of wind and sea wave loading is significantly affected by the monopile-soil interaction. A 3D soil domain was adopted and it was represented with 72,400 linear hexahedral elements of type C3D8. To eliminate the boundary effects, a sufficiently large cylindrical soil domain with a diameter of 20D and a height of 1.7L (D and L being respectively the monopile outer diameter and its embedded depth) was adopted (see Figure 2). Only the horizontal displacement was fixed on the cylindrical lateral boundaries while both the horizontal and vertical displacements were restrained on the bottom boundary. The small sliding, surface-to-surface and master/slave contact pair formulation implemented in Abaqus [START_REF] Abaqus | Finite element analysis program[END_REF] was used to model the contact interaction between the inner/outer surfaces of the monopile (master surfaces) and the surfaces of the soil (slave surfaces) around and inside the monopile. The classical Coulomb friction model was used to describe the frictional behavior between the soil and the monopile.

An interface friction angle of δ= 2 3 φ was used where φ is the peak friction angle of the soil which can be obtained from the soil critical state friction angle φ c and the soil relative density D r by the correlation given in literature by [START_REF] Wichtmann | An experimental data base for the development, calibration and verification of constitutive models for sand with focus to cyclic loading. Part I: Tests with monotonic loading and stress cycles[END_REF] as follows:

φ=φ c exp (0.293 D r 1.76 ) (1)

Soil constitutive models

This paper makes use of the SANISAND version 2004 model given by [START_REF] Dafalias | Simple plasticity sand model accounting for fabric change effects[END_REF] and the hypoplastic model with intergranular strain (IS) developed by [START_REF] Wolffersdorff | A hypoplastic relation for granular materials with a predefined limit state surface[END_REF][START_REF] Niemunis | Hypoplastic model for cohesionless soils with elastic strain range[END_REF] to describe the soil behavior as predicted by the two constitutive models when performing the nonlinear dynamic simulations of the large diameter monopile-supported 10 MW OWT. It should be noted herein that these advanced soil constitutive models were extensively validated in literature based on laboratory model tests that simulate the soil mechanical behavior under complex cyclic loads.

Hence, for the analysis presented in this paper, the two advanced soil constitutive models are expected to provide reliable results when used to simulate the response of a real-OWT structure where the soil experiences complex dynamic and cyclic loadings.

SANISAND

SANISAND is the name generally used in literature for a family of Simple ANIsotropic SAND constitutive models. This model was firstly developed by [START_REF] Manzari | A critical state two-surface plasticity model for sands[END_REF] to realistically simulate the stress-strain behavior of sands under monotonic and cyclic, drained and undrained loading conditions. It makes use of the critical state soil mechanics concept and the bounding surface plasticity model. Later improvements and extensions were carried out by [START_REF] Dafalias | Simple plasticity sand model accounting for fabric change effects[END_REF][START_REF] Papadimitriou | Plasticity model for sand under small and large cyclic strains: a multiaxial formulation[END_REF][START_REF] Papadimitriou | Plasticity model for sand under small and large cyclic strains[END_REF][START_REF] Dafalias | SANISAND-Z: zero elastic range sand plasticity model[END_REF] regarding the dilatancy, the effect of fabric changes during loading and the response to radial stress paths. It should be noted herein that the version 2004 given by [START_REF] Dafalias | Simple plasticity sand model accounting for fabric change effects[END_REF] was employed in this paper due to its simplicity and its reasonable predictive potential [START_REF] Kementzetzidis | Geotechnical aspects of offshore wind turbine dynamics from 3D non-linear soil-structure simulations[END_REF][START_REF] Ramirez | Site response in a layered liquefiable deposit: evaluation of different numerical tools and methodologies with centrifuge experimental results[END_REF][START_REF] Taiebat | Propagation of seismic waves through liquefied soils[END_REF].

The SANISAND 2004 version [START_REF] Dafalias | Simple plasticity sand model accounting for fabric change effects[END_REF] relies on the hypo-elastic small-strain shear (G) and bulk moduli (K ) to characterize the elastic region whose size is defined by a parameter m. The small-strain shear modulus G is pressure-and void ratio-dependent and is defined using a dimensionless parameter G 0 . The bulk modulus K is given by the shear modulus G and the Poisson ratio ν. Thus, the hypo-elastic shear and bulk moduli are function of two material parameters ( G 0 and ν). Concerning the critical-sate surface, it is defined by two parameters (

M c
and M e ) in the ṕ-q space, and by three parameters ( λ c , e 0 and ξ) in the ṕ-e space. The hardening behavior of the SANISAND model in the deviatoric space is described by the plastic modulus which is governed by three parameters ( h 0 , c h and n b ). The dilatancy which is defined as the ratio between the plastic volumetric strain and the plastic deviatoric stain is given using two parameters ( A 0 and n d ¿ . Finally, the amount of soil contraction and dilation is enhanced in this model via an equation describing the evolution of the fabric change and which depends on two user-defined parameters ( z max and c z ).

In summary, the SANISAND 2004 version requires in total fifteen material parameters ( G 0 , ν, m, M c , M e , λ c , e 0 , ξ, h 0 , c h , n b , A 0 , n d , z max and c z ).

Hypoplasticity

In this work, the constitutive equations of hypoplasticity given by [START_REF] Wolffersdorff | A hypoplastic relation for granular materials with a predefined limit state surface[END_REF] and further extended for the intergranular strain (IS) approach as proposed by [START_REF] Niemunis | Hypoplastic model for cohesionless soils with elastic strain range[END_REF] are used for the sake of comparison with the SANISAND model to simulate the sand behavior under the dynamic environmental loadings. The extended hypoplastic constitutive model proves to be capable for the modeling of the nonlinear and anelastic behavior of granular materials as it can capture some of the main interesting sand properties under dynamic loading like contractancy, dilatancy, the dependency of stiffness and strength on the pressure and void ratio, in addition to the accumulative effects and the hysteretic sand behavior under cyclic loading.

Within the hypoplastic model, three curves in the e-p space are used to define the maximum ( e i ), critical ( e c ) and minimum ( e d ¿ void ratios. These curves may be described using five parameters ( h s , n, e d 0 , e c 0 and e i 0 ) as given by [START_REF] Bauer | Calibration of a comprehensive hypoplastic model for granular materials[END_REF]. The parameters h s (granular hardness) and n (compression exponent) control the slope and the curvature of the curves respectively and the parameters ( e d 0 , e c 0 and e i 0 ) control the curves position at p=0 Pa. The dependency of the peak friction angle on the relative void ratio and the change of stiffness with the change of relative density are controlled within the hypoplastic model via the parameters α and β respectively. Finally, notice that the hypoplastic formulation depends on the soil critical state friction angle φ c . As a conclusion, the basic hypoplastic model by [START_REF] Wolffersdorff | A hypoplastic relation for granular materials with a predefined limit state surface[END_REF] has eight material parameters: φ c , h s , n, e d 0 , e c 0 , e i 0 , α and β.

To prevent excessive strain accumulation under cyclic loads, the IS concept [START_REF] Niemunis | Hypoplastic model for cohesionless soils with elastic strain range[END_REF] was proposed where five additional parameters (R, m R , m T , β R and χ) are required. Within the IS approach, there is a certain strain range (whose size is defined via the maximum intergranular strain parameter R) for which the incremental stiffness remains almost constant. Two factors m R and m T are used to characterize the stiffness increase upon load reversals. Finally, the smoothing of stiffness variation is described by the two parameters β R and χ, where β R can also influence the evolution of intergranular strain.

In total, the extended hypoplastic constitutive model [START_REF] Wolffersdorff | A hypoplastic relation for granular materials with a predefined limit state surface[END_REF][START_REF] Niemunis | Hypoplastic model for cohesionless soils with elastic strain range[END_REF] requires thirteen material parameters to describe the granular soil medium.

Numerical implementation

The SANISAND elastoplastic model 2004 version [START_REF] Dafalias | Simple plasticity sand model accounting for fabric change effects[END_REF] and the hypoplastic model with intergranular strain (IS) [START_REF] Wolffersdorff | A hypoplastic relation for granular materials with a predefined limit state surface[END_REF][START_REF] Niemunis | Hypoplastic model for cohesionless soils with elastic strain range[END_REF] are implemented in Abaqus/Standard [41] via user-defined material subroutines (UMATs) available at the SoilModels website [START_REF] Gudehus | The soilmodels.info project[END_REF]. In all of the numerical simulations performed in this paper, a homogeneous soil medium made of Karlsruhe fine sand is selected. The parameters of the Karlsruhe fine sand have been extensively calibrated in literature for simulations on FE applications [START_REF] Wichtmann | Inspection of three sophisticated constitutive models based on monotonic and cyclic tests on fine sand: Hypoplasticity vs. Sanisand vs. ISA[END_REF][START_REF] Fuentes | A p-y Model for Large Diameter Monopiles in Sands Subjected to Lateral Loading under Static and Long-Term Cyclic Conditions[END_REF][START_REF] Poblete | On the simulation of multidimensional cyclic loading with intergranular strain[END_REF][START_REF] Fuentes | Evaluating the performance of an ISAhypoplasticity constitutive model on problems with repetitive loading[END_REF]. In this study, the parameters of the Karlsruhe fine sand based on the previous calibration performed by [START_REF] Wichtmann | Inspection of three sophisticated constitutive models based on monotonic and cyclic tests on fine sand: Hypoplasticity vs. Sanisand vs. ISA[END_REF] (for both the SANISAND and the hypoplastic models) are adopted. It is worth to mention that the calibration performed by [START_REF] Wichtmann | Inspection of three sophisticated constitutive models based on monotonic and cyclic tests on fine sand: Hypoplasticity vs. Sanisand vs. ISA[END_REF] was based on extensive laboratory tests data [START_REF] Wichtmann | [END_REF][START_REF] Wichtmann | An experimental data base for the development, calibration and verification of constitutive models for sand with focus to cyclic loading. Part I: tests with monotonic loading and stress cycles[END_REF][START_REF] Wichtmann | An experimental data base for the development, calibration and verification of constitutive models for sand with focus to cyclic loading. Part II: tests with strain cycles and combined cyclic and monotonic loading[END_REF] and involved three constitutive models (SANISAND, Hypoplasticity and ISA). Table 2 and Table 3 summarize the material parameters adopted in this paper for the SANISAND 2004 version (Table 2) and the hypoplastic model with IS (Table 3).

Loads

The OWT dynamic response was studied in this paper using two different 10 min loading scenarios from the UpWind project [START_REF] Fischer | Upwind Design Basis[END_REF]. This project involves an offshore Dutch site (called K13) located in the North Sea. The K13 site is suitable for monopile foundations in shallow water depths. Load case 10 (LC10) and load case 17 (LC17) from the K13 site were used in this study (cf. Table 4). LC10 (where the mean wind speed is between the cut-in and the cutout speed of the DTU 10 MW) represents the design load case (DLC) 1.2 Power production given by IEC 61400-3 [START_REF]Wind Turbines -Part 3: Design requirements for offshore wind turbines[END_REF], while LC17 (with wind speed above cut-out speed) corresponds to DLC 6.4 Parked (i.e. standing still or idling). The decision behind choosing two different DLCs is to compare the OWT responses as obtained using the two advanced well-calibrated soil models for different strain levels in the soil. In this study, for each LC, the wind load acting along the tower and blades and the wave load acting along the monopile in water were stochastically simulated and are briefly described in the following subsections.

Wind loads

The wind loads acting on the OWT include the drag load along the OWT tower and the aerodynamic loads on the OWT blades.

Dragloadingalongthetower

In this paper, the OWT tower was divided into sixteen segments in the simulation (see Figure 3). The length of the seven bottom segments (T1-T7) is 11.5 m and the nine other segments (T8-T16) is 5 m. The nodal wind drag force F i (t ) exerted on the OWT tower may be represented by a mean and a fluctuating component [START_REF] Bhattacharya | Design of foundations for offshore wind turbines[END_REF]. The mean wind velocity varies with height from the ground surface. In the present work, the logarithmic law recommended by IEC [START_REF]Wind Turbines -Part 3: Design requirements for offshore wind turbines[END_REF] is used to model this variation. Concerning the fluctuating wind velocity, the Kaimal spectrum [START_REF] Jonkman | [END_REF] was employed to model its power spectral density (PSD) function.

A method to simulate the spatially correlated turbulent time histories following the simulation algorithm with random phases as proposed by Shinozuka [START_REF] Shinozuka | Monte Carlo solution of structural dynamics[END_REF] and Deodatis [START_REF] Deodatis | Simulation of ergodic multivariate stochastic processes[END_REF] was implemented in Matlab making use of the Davenport coherency function [START_REF] Davenport | The spectrum of horizontal gustiness near the ground in high winds[END_REF]. The simulated and target Kaimal PSDs of the fluctuating wind velocity were found to be in very good agreement (results not shown for conciseness).

Aerodynamicloadsontheblade

The wind speed at any node along the blade not only reflects the stochastic property of wind speed (fluctuating component of wind speed) but also involves the change in the mean wind speed induced by a periodic variation of the node position in space due to the blade rotation.

In this paper, a Matlab code was developed to generate a 1D wind velocity field that covers the rotor disk domain and which is compatible with the Kaimal spectrum. The Sandia National Laboratory (SNL) method developed by Veers [START_REF] Veers | Three-dimensional wind simulation[END_REF] was used for this purpose. The generated 1D wind velocity field contains time series of the fore-aft component of the wind speed vector at discrete nodes along the rotating blades.

The aerodynamic loads acting on the rotating blades were obtained making use of the Blade Element Momentum (BEM) [START_REF] Hansen | Aerodynamics of wind turbines[END_REF] theory. In the present study, a Matlab code was developed to find the axial and tangential induction factors, within the BEM method, using the iterative algorithm suggested by Hansen [START_REF] Hansen | Aerodynamics of wind turbines[END_REF]. Prandtl tip loss and Gluaert correction [START_REF] Hansen | Aerodynamics of wind turbines[END_REF] were included in the Matlab code to improve the accuracy of the calculated aerodynamic loads.

Notice also that for the operating case (LC10), a 9.6 RPM rotor speed and a 16.968° pitch angle were used for the simulation of the aerodynamic loads whereas for the nonoperating case (LC17), a zero-rotor speed and a 90 degrees blades pitch angle were set to make the blades pitched out of the wind [START_REF] Jonkman | AeroDyn v15 User's Guide and Theory Manual National Renewable Energy Laboratory[END_REF].

Figure 4 represents the total in-plane (tangential) and out-of-plane (normal) aerodynamic loads acting on blade 1 in LC10. In the present study, each blade was divided into 37 segments (see Figure 3) and the aerodynamic loads per unit length were assumed to be the same within each segment. The length of each segment is given in Figure 3 and it is consistent with the data provided in the turbine specifications [START_REF] Bak | Description of the DTU 10 MW reference wind turbine[END_REF].

Hydrodynamic loads on the monopile

In this paper, the time series of the irregular sea water surface elevation were generated using the JONSWAP spectrum [START_REF] Hasselmann | Measurements of windwave growth and swell decay during the Joint North Sea wave project (JONSWAP) Technical Report Deutches Hydrographisches Institut[END_REF] given as:

S η (f )=0.3125 H s 2 T p ( f f p ) -5 exp [ -5 4 ( f f p ) -4 ] (1-0.287 log γ ) γ exp( [ (ω-ωp) 2 2 σ 2 ω p 2 ] ¿)¿ (2)
where T p and H s are respectively the wave period and the wave significant height, f p =1/T p , ω p =2 π f p , σ is a constant defined as σ =0.07 for f ≤ f p and σ =0.09 for f >f p and finally, γ is the peak enhancement factor defined as a function of T p and H s [START_REF] Bhattacharya | Design of foundations for offshore wind turbines[END_REF].

The generation of the time histories was performed by partitioning the spectrum into a finite number of components where each component represents a harmonic wave characterized by its amplitude, angular frequency and a random phase angle. Then, by superimposing the harmonic waves, time series of the irregular wave surface elevation were generated. Once the time series were obtained, the wave loads along the monopile in water were computed making use of the Morison equation (following the specifications given in DNV [START_REF]DNV-OS-J101 offshore standard: design of offshore wind turbine structures[END_REF]) as the sum of a linear inertia component and a nonlinear drag component. It is given as:

F w = π 4 D p 2 ρ w C m ux +0.5 ρ w C d D p u x | u x | (3)
where C m and C d are the inertia and drag coefficients respectively ( C m =2 and C d =1.2 being the values adopted in the simulation), ρ w is the water density (1025 kg /m 3 ¿ and D p is the outer monopile diameter.

In the numerical simulations, the monopile in water was divided into ten segments (see Figure 3) of equal length (2.5 m). A very good agreement was found between the generated and target JONSWAP PSDs of the sea surface elevation for LC10 and LC17 (results not shown in this paper for conciseness).

Vibration characteristics and damping

In this section, a structural modal analysis of the OWT was performed in Abaqus/Standard to calculate the OWT natural frequencies. The obtained natural frequencies were used to model the damping of the OWT in the 3D mechanical model. It should be mentioned herein that the natural frequency calculation is concerned with very small amplitude vibrations [START_REF] Wichtmann | Inspection of three sophisticated constitutive models based on monotonic and cyclic tests on fine sand: Hypoplasticity vs. Sanisand vs. ISA[END_REF][START_REF]Eurocode 1: Actions on structures -general actions -wind actions[END_REF] and thus, the consideration of the elastic properties (initial stiffness) of the soil would suffice for this computation. For the FE model used in this study, the depth variation profile of the hypoelastic small-strain shear modulus G obtained from the SANISAND 2004 constitutive equations [START_REF] Li | Dilatancy for cohesionless soils[END_REF][START_REF] Richart | Vibration of soils and foundations[END_REF] (Equation 4) was used to compute the soil Young modulus profile

Notice that G 0 in Equation 4 is a dimensionless material constant (see Table 2), e is the current void ratio, p is the mean pressure and p atm is the atmospheric pressure. An initial void ratio of 0.865 (corresponding to a relative density of 50%) was considered in the calculation (as will be shown later in this paper). Table 5 provides the results of the modal analysis. As may be seen from this table, the major mode shapes of the OWT are the first bending modes of the tower in the side-to-side and the fore-aft directions.

As the natural frequencies of monopile-supported OWTs are generally very close to the excitation frequencies, damping is a critical parameter and it should be carefully chosen to reliably predict the dynamic behavior of the OWT. The main reason is that damping is the only factor that limits the OWT amplitude of response at resonance. For a monopilesupported OWT, there are four sources of damping: aerodynamic, hydrodynamic, soil and structural [START_REF] Bisoi | Dynamic analysis of offshore wind turbine in clay considering soil-monopile-tower interaction[END_REF][START_REF] Arany | Closed form solution of eigen frequency of monopile supported offshore wind turbines in deeper waters incorporating stiffness of substructure and SSI[END_REF][START_REF] Chen | Modelling damping sources in monopile-supported offshore wind turbines[END_REF]. Aerodynamic damping is generated by the air drag experienced by the rotating blades and it is strongly influenced by the wind speed, the rotation speed and the pitch control system. Hydrodynamic damping results from two different sources: (i) wave radiation and (ii) viscous damping. The wave radiation is proportional to the wave velocity whereas viscous damping is proportional to the square of the relative velocity and it is very low in the low frequency vibration regime of OWTs; thus, the larger contribution in our case results from wave radiation. Soil is another source of damping in OWTs and it includes the energy dissipated due to hysteretic damping (plastic deformation) and radiation damping (wave spreading away from the monopile). Radiation damping is typically negligible for an OWT where the loading frequencies are below 1 Hz. Finally, concerning the structural damping, it depends mainly on the energy dissipated in the material in the form of heat in addition to that dissipated in the different connections (welded and grouted connections, etc.).

It is worth mentioning herein that damping in OWTs has been extensively investigated in literature; however, the variations in published values for most damping sources remain large, which may result in a lot of uncertainty upon the prediction of the dynamic response of the OWT. In this regard, [START_REF] Chen | Modelling damping sources in monopile-supported offshore wind turbines[END_REF] presented a literature review on the existing damping values for OWTs where aerodynamic damping in the fore-aft direction for an OWT in operational condition was found to have a wide variation range (4% -8%) followed by structural damping (0.2% -1.5%), soil damping (0.17% -1.3%) and hydrodynamic damping (0.11% -0.39%).

Also, [START_REF] Arany | Closed form solution of eigen frequency of monopile supported offshore wind turbines in deeper waters incorporating stiffness of substructure and SSI[END_REF] recommended damping ratios that are in conformity with those presented in [START_REF] Chen | Modelling damping sources in monopile-supported offshore wind turbines[END_REF] and highlighted the importance of analyzing the vibrations of the OWT in both directions; the total damping of the first mode of vibration in the fore-aft direction (2% -8%) being typically higher than that in the side-to-side direction (1% -4%).

In this paper, damping was modelled in the case of an operational OWT and by considering only the fore-aft vibration as the wind and wave loads acting on the tower and blades are larger in that direction in comparison to those in the side-to-side direction. The chosen damping ratios are selected to be within the ranges recommended by [START_REF] Arany | Closed form solution of eigen frequency of monopile supported offshore wind turbines in deeper waters incorporating stiffness of substructure and SSI[END_REF][START_REF] Chen | Modelling damping sources in monopile-supported offshore wind turbines[END_REF] and are given as follows:

 3.5% for the aerodynamic damping in the fore-aft direction as given by [START_REF] Zuo | Using multiple tuned mass dampers to control offshore wind turbine vibrations under multiple hazards[END_REF],

 0.43% and 1% for the structural damping for the blades and tower respectively based on the turbine specifications [START_REF] Bak | Description of the DTU 10 MW reference wind turbine[END_REF],

 0.12% for the hydrodynamic damping following [START_REF] Leblanc | Design of offshore wind turbine support structure structures[END_REF],

 no additional soil hysteretic damping was assumed, as it was considered only through the soil constitutive model equations.

Summing all components together, which is an acceptable assumption for lightly damped structures like OWTs, the total damping ratio in the fore-aft direction for the rotating blades is 3.93% (3.5% for aerodynamic damping + 0.43% for structural damping). For the tower, the total damping is 1.12% (1% for structural damping + 0.12% for hydrodynamic damping). The damping of the OWT is modelled in Abaqus by means of material Rayleigh damping [START_REF] Ak | Dynamics of structures[END_REF] and the first in-plane and out-of-plane vibration frequencies of the blades and tower (Table 5)

were used to calculate the mass (α) and stiffness (β) Rayleigh damping coefficients for the blades and tower respectively. α and β are therefore: 0.169 and 8.5 × 10 -3 for the blade and 0.014 and 8.98 × 10 -3 for the tower.

Numerical simulations

In the 3D soil domains, an initial stress field was defined by assuming an at-rest lateral earth pressure condition through the relation σ 11 = σ22 =K 0 σ33 where σ 11 and σ 22 are the effective lateral earth pressures at-rest and σ 33 is the effective vertical overburden pressure. The K 0 lateral earth pressure coefficient was set to K 0 =1-sin φ c according to Jacky's relationship, where φ c is the critical state friction angle of the soil. Note that for the hypoplastic model, φ c =¿ 33.1° (see Table 3), whereas for the SANISAND model, φ c is not a material parameter but it can be determined by calculation using M c =6 sin φ c /(3-sin φ c ) where the ultimate critical state stress ratio M c =¿ 1.34(see Table 2). The initial effective vertical stresses ( σ 33 ) in the 3D soil domain were computed by assuming an effective soil unit weight of γ ' =9 kN /m 3 .

An initial void ratio e 0 of 0.865 for both soil models was considered. It was computed for a relative density D r =50 % as follows:

e 0 =e c 0 -D r (e c 0 -e d 0 ) (5)
where e d 0 and e c 0 are the minimum and critical void ratios for a zero mean effective pressure (Table 3). Notice that e d 0 coincides with the minimum void ratio according to ASTM D4253-16e1 [START_REF]Standard test methods for maximum index density and unit weight of soils using a vibratory table[END_REF] and e c 0 coincides with the maximum void ratio e c 0 = e max according to ASTM D4254-16 [START_REF] Astm | Standard test methods for minimum index density and unit weight of soils and calculation of relative density[END_REF].

To perform the dynamic analysis, three steps were performed. In step 1 named "geostatic step", a gravitational acceleration was applied to the 3D soil domain to accomplish geostatic equilibrium with the initially defined soil stress field. In step 2, the weight of the whole OWT was applied. Finally, in step 3 the nonlinear dynamic analysis was performed where the OWT was subjected to the stochastically simulated 10 min aerodynamic and hydrodynamic loadings. It should be noted that the 10 minutes duration is considered as sufficient for steady environmental conditions.

The dynamic equations were solved within Abaqus/Standard version 6.14 using the Hilber-Hughes-Taylor (HHT) [41] implicit integrator. As for time marching, the step sizes ∆ t were tuned to the loading amplitude to ensure a convenient discretization of the load prior to solving the dynamic equations. Maximum time increments of ∆ t=¿ 0.1 s (for the load case with lower amplitude i.e. LC10) and of ∆ t=¿ 0.075 s (for the load case with higher amplitude i.e. LC17) were found suitable to achieve an accurate prediction of the OWT displacement and rotation. For all simulations in the SANISAND and hypoplastic soil domains, all of the solution-dependent state variables (SDVs) were set to zero except for the initial void ratio where SDV7 = e 0 = 0.865.

Numerical results

This section presents the numerical results obtained based on the 3D detailed FE mechanical model described above. Firstly, the nonlinear structural and geotechnical dynamic responses of the OWT as obtained using the two advanced sand constitutive models (calibrated on the same Karlsruhe fine sand) are compared and discussed. The aim is to highlight some relevant aspects related to the OWT dynamics. The significance of using such advanced soil constitutive models is demonstrated by comparing their results with those provided by the commonly used elastic-perfectly plastic Mohr-Coulomb soil model. In addition, the 3D FE model is exploited (i) to investigate the influence of the sand relative density on the dynamic responses of the OWT and (ii) to emphasis the importance of considering a large deformation analysis for the superstructure and the monopile foundation. Secondly, a comparison of the soil damping as obtained when using the two advanced soil constitutive models is presented.

Finally, the suitability of two simplified foundation models suggested in literature for the dynamic analysis of OWTs is investigated.

1. OWT structural and geotechnical responses

The structural dynamic responses of the OWT (displacement and rotation time histories) vary along the tower height and along the blades length. For conciseness and due to space limitation, only the OWT responses at the tower top and at mudline are presented and discussed in the present study.

Figures 5a and6a show the fore-aft displacement time histories at the tower top for LC10 and LC17 respectively and Figures 7a and8a give the corresponding side-to-side displacement time histories. In these figures, the black curves are the results for the hypoplastic model, and the red curves are those for the SANISAND model. The comparison of Figure 5a with Figure 7a and Figure 6a with Figure 8a indicates that the displacements of the tower top in the sideto-side direction are (as expected) much smaller than those in the fore-aft direction for both load cases and both sand constitutive models. The maximum absolute tower top displacements in the fore-aft and side-to-side directions for LC17 are larger by 67.5% and 20% respectively for the hypoplastic constitutive model (and by 59% and 13% respectively for the SANISAND constitutive model) when compared to LC10. It should be noted that the average computation time for simulating the 10 min nonlinear dynamic response of the OWT was about 2.1 weeks for the hypoplastic constitutive model and 2.9 weeks for the SANISAND constitutive model when using a workstation with a memory of 128 GB RAM, a storage of 128GB SSD/1 TB Sata and a processor 2 x Intel® Xeon® i8.

To compare the OWT time series structural response as obtained using the two soil constitutive models, the peak and root mean square (RMS) values are used in this study. It should be noted that the RMS is very important when analyzing vibration data in time domain as it is related to the energy content of the response profile and thus to the destructive capability of the vibration. The peak and RMS relative differences (i.e. R peak and R RMS ) are calculated in this paper by assuming that the values predicted by the SANISAND constitutive model are the reference values. They are given as follows:

R peak = U hypo Peak -U San Peak U San Peak (6) R RMS = U hypo RMS -U San RMS U San RMS (7)
Notice that (U hypo Peak , U San Peak ) and (U hypo RMS , U San RMS ) in these equations are respectively the peak and RMS responses of the OWT as predicted by the hypoplastic and SANISAND constitutive models.

For LC10 (Figure 5a), one may observe that both constitutive models result in almost the same fore-aft tower top displacement where the peak and RMS relative differences are relatively negligible (2.1% for the peak and 0.2% for the RMS). However, when the OWT is subjected to the higher environmental condition 'LC17' (Figure 6a), the fore-aft displacements are no longer the same and the SANISAND model exhibits a slightly stiffer response compared to the hypoplastic model where the peak and RMS relative differences were found to be 8.7% and 1.2 % respectively.

In addition to the time series responses, Figures 5b and6b show the frequency responses (PSDs) of the fore-aft tower top displacement for LC10 and LC17 respectively and Figures 7b and8b show the corresponding curves for the side-to-side tower top displacement. Notice that the PSDs describe the frequency content of the responses, with peaks representing the vibration modes of the OWT.

As shown in Figure 5b, two peaks (with almost the same amplitude) appear at 0.1972 Hz for the hypoplastic model and at 0.1976 Hz for the SANISAND model, which correspond to the first fore-aft natural frequency of the tower. This means that the first vibration mode of the tower is excited by the simulated external environmental loads. For the PSDs in the case of LC17 (see Figure 6b), higher amplitude peaks are obtained compared to the peaks in the case of LC10, but at slightly lower frequencies due to the possible reduction of the soil stiffness.

Indeed, a peak at 0.195 Hz for the hypoplastic model (instead of 0.1972 Hz for LC10) and a peak at 0.1959 Hz for the SANISAND (instead of 0.1976 Hz for LC10) can be observed.

From the results presented in Figure 6b, it is obvious that the peak amplitude (vibration energy) is somewhat smaller in the case of the SANISAND by around 15.4% and the peak frequency is slightly higher by only 0.46%. This is due to the marginally stiffer response predicted by the SANISAND model at high strain levels.

Concerning the results obtained in Figures 7b and8b in the side-to-side direction, one may observe that the PSDs of the displacement response obtained using both constitutive models possess peaks at the same frequency (≈ 0.197 Hz) and are of almost equal amplitude. The frequency corresponding to the peak (≈ 0.197 Hz) denotes the first side-to-side natural frequency of the tower.

Figures 9a and 9b illustrate the out-of-plane rotation time histories at the tower top for LC10

and LC17 respectively as given by the two constitutive models. The same observations made for Figures 5a and 6a remain valid herein. Indeed, the SANISAND model exhibits a stiffer response with respect to the hypoplastic model for LC17 where the peak and RMS relative differences are around 9% and 0.93% respectively; these relative differences being quasi negligible for LC10. Notice herein that the obtained tower top rotations are in accordance with the design guidelines (maximum rotation is limited to 5° at the tower top).

While the focus has been put so far on the structural dynamic responses of the OWT at the tower top, it is worthwhile to investigate the geotechnical dynamic responses of the OWT.

Figures 10a and10b show the maximal principal logarithmic strain distribution in the hypoplastic and SANISAND soil domains at the end of the 10 min simulations for LC10 and LC17 respectively. From these figures, one may observe that the zone having a strain level beyond 10 -2 (in grey) is much more extended in the case of LC17. Figure 10a also shows that relatively the same level of strain is attained in both the hypoplastic and the SANISAND constitutive models for LC10 where the corresponding values are respectively equal to 3.440 ×10 -2 and 3.325 ×10 -2 . However, for LC17, higher strains are produced by the hypoplastic model ( 8.245 ×10

-2

) compared to those developed by the SANISAND model ( 7.081× 10 -2 ). Overall, the obtained results confirm that the SANISAND constitutive equations are predicting a slightly stiffer response than the hypoplastic model.

To gain further insight into the geotechnical responses of the OWT, Figures 11 and12 illustrate respectively the monopile fore-aft displacement and out-of-plane rotation time histories at mudline using both the hypoplastic and the SANISAND constitutive models under LC10 and LC17. From these figures, one can observe that the SANISAND model provides smaller responses for both the fore-aft displacement and out-of-plane rotation as compared to the hypoplastic model. Comparing Figure 11a and Figure 12a with Figure 5a and Figure 9a respectively, one may observe that for LC10 although the choice of a soil constitutive model has a negligible impact on the structural responses at the tower top, it has a notable effect on the monopile responses at mudline. This effect is more pronounced for LC17 as may be seen when comparing Figure 11b and Figure 12b with Figure 6a and Figure 9b respectively. Hence, the prediction of the geotechnical responses strongly depends (as expected) on the choice of the soil constitutive model.

ComparisonwiththeconventionalMohr-Coulombconstitutivemodel

This section aims at examining the adequacy of using the conventional elastic perfectly plastic soil constitutive model with a Mohr-Coulomb failure criterion when performing the dynamic simulations of OWTs. For this purpose, a 10-minute dynamic simulation was performed under LC17 where the behavior of the 3D soil domain was described by using the elastic perfectly plastic soil model with the Mohr-Coulomb failure criterion based on the soil properties given in Table 6. The obtained structural (fore-aft tower top displacement) and geotechnical (fore-aft monopile displacement at mudline) responses were compared with those obtained when using the more advanced soil constitutive models (Hypoplastic with IS and SANISAND).

Figures 13a and13b show respectively the time history and the frequency response (PSD) of the fore-aft tower top displacement under LC17 and Figure 14 illustrates the monopile fore-aft displacement at mudline. In these figures, the blue curves are the results obtained when using the Mohr-Coulomb constitutive model, and the black and red curves are those obtained when employing the hypoplastic and SANISAND models respectively.

The OWT time series as obtained using the Mohr-Coulomb constitutive model were evaluated by calculating the peak and RMS relative differences with respect to the SANISAND model (reference model). For the structural response (fore-aft tower top displacement, Figure 13a), one may observe that the Mohr-Coulomb model resulted in a stiffer response compared to both advanced constitutive models where the peak and RMS differences were found to be 11.2% and 7.6% respectively. These relative differences (peak and RMS) are much more important when observing the geotechnical response (monopile fore-aft displacement at mudline, Figure 14) obtained by the Mohr-coulomb model where they were found to be 15.1% (peak) and 13.9% (RMS). Overall, the obtained results highlight the importance of using advanced soil constitutive models suitable of describing the soil dynamic behavior.

Also, they confirm the observations obtained before in which the soil constitutive model was found to have a notable effect on the monopile responses compared to the structural responses at the tower top.

EffectofthesandrelativedensityontheOWTdynamicresponses

This section is devoted to examine the impact of the sand relative density on the dynamic response of the OWT. For this purpose, three initial relative densities (Dr=¿ 20%, 50% and 80%) are considered. The hypoplastic model with IS was chosen in the dynamic simulations due to its computational efficiency compared to the SANISAND constitutive model. All simulations are performed under LC17.

Figures 15a and15b show respectively the time histories and the corresponding PSDs responses of the fore-aft tower top displacement for the three relative densities (Dr=¿ 20% in black, Dr=¿ 50% in red and Dr=¿ 80% in blue). The tower top vibration decreases with the increase of the relative density. Indeed, the peak and RMS reduction ratios of the tower top displacement (calculated relative to the reference values at Dr=¿ 20%) are respectively around 8% and 4.7% for Dr=¿ 50% and 16.8% and 7.7% for Dr=¿ 80%. Figure 15b reports the peaks in the frequency spectrum for the three relative densities where the corresponding frequencies represent the first natural frequency of the tower in the fore-aft direction. From this figure, one may observe that the vibration frequency increases with the increase of the sand relative density; the sand with lower density (i.e. Dr=¿ 20%) is simulating a softer response as it exhibits a peak at ≈ 0.1915 Hz with a higher energy. The relative increase of the peak frequencies (calculated relative to the peak frequency 0.1915 Hz for Dr=¿ 20%) are around 1.8% for Dr=¿ 50% and 3% for Dr=¿ 80%.

Figure 16 gives the maximal principal logarithmic strain distribution in the hypoplastic soil domain for the three relative densities at the end of the 10 min nonlinear time domain simulations for LC17. High straining can be clearly observed for the lowest relative density where the sand with Dr=¿ 20% develops a straining that is more than five times greater than that developed in the case with Dr=¿ 80%. This result is due to the lower soil stiffness present in the sand with Dr=¿ 20%.

To evaluate the influence of the sand relative density on the monopile responses, Figures 17a and17b present respectively the monopile fore-aft displacement and out-of-plane rotation time histories at mudline under LC17 for the three relative densities. From both figures, one may observe that the relative density strongly affects the monopile responses at mudline where a significant vibration reduction occurs as the sand relative density increases. Indeed, the peak and RMS reduction ratios are respectively 34.6% and 29.6% for the monopile displacement and 21.8% and 25% for the rotation as the relative density increases from 20%

to 80%. Figures 17a and 17b also show that the average values of the displacement and rotation time histories at mudline decrease with the increase of the sand relative density.

Effectofperformingalargedeformationanalysis

This section aims at examining the effect of considering the large deformation of the OWT superstructure and monopile foundation while performing the nonlinear dynamic simulations.

For this purpose, a nonlinear dynamic analysis was performed using Abaqus without considering the large deformations of the foundation and the superstructure and the obtained results were compared with those obtained when considering the large deformations. Notice herein that in a large deformation analysis, the local material directions rotate with the deformation in each element [41].

The hypoplastic model with IS was used in this section and the numerical dynamic simulations were performed under LC17. Figure 18 and Figure 19 give respectively the foreaft displacement time history at the tower top and the monopile fore-aft displacement at mudline as obtained by considering or by neglecting the superstructure and monopile foundation large deformation. From both figures, one may observe the importance of considering a large deformation analysis when performing nonlinear dynamic simulations as the peak and RMS relative differences were found to be respectively 10.3% and 8.4% for the tower top displacement and 14.4% and 12% for the monopile displacement at mudline.

2.

Estimation and comparison of soil damping

This section aims at estimating and comparing the soil damping as obtained using the two calibrated soil constitutive models (SANISAND and hypoplastic) making use of the developed 3D FE model. For this purpose, two dynamic simulations (one for each soil constitutive model) were performed within Abaqus for LC17 where the OWT was subjected during 20s to forced vibration (the first 20s of LC17 time history) followed by sudden unloading to simulate numerically the OWT free vibration time history response. It should be noted herein that the 20s duration of forced vibration was deemed to be sufficient to develop large straining in the soil (for both soil constitutive models) where a strain level beyond 10 -2 was observed in the vicinity of the monopile at mudline prior to unloading. The total damping ξ tot of the system was computed using the logarithmic decrement method applied to the free vibration time history of the fore-aft tower top displacement (see Figure 20) and it is given as follows:

ξ tot = 1 √1+(2π/δ ) 2 (8)
where δ is the logarithmic decrement calculated as:

δ= 1 n ln ²²² ( A 1 A n ) (9)
in which A 1 and A n are two successive amplitudes n periods apart. The soil damping ξ soil as simulated herein (for a given soil constitutive model) can be estimated as follows:

ξ soil =ξ tot -ξ hydro -ξ struct [START_REF] Kementzetzidis | Geotechnical aspects of offshore wind turbine dynamics from 3D non-linear soil-structure simulations[END_REF] where a structural damping ξ struct of 0.19% based on Eurocode [START_REF]Eurocode 1: Actions on structures -general actions -wind actions[END_REF] and a hydrodynamic damping ξ hydro of 0.12% were considered in the simulations via the Rayleigh damping model in Abaqus. Notice that in this section the adopted structural damping (0.19%) is different from that defined in the turbine specifications (1%) and the aerodynamic damping was neglected.

The main reason behind these considerations is to avoid their influence on the estimation of soil damping.

From Figure 20 and Equations 8-10, a slightly higher soil damping ( ξ soil = 0.58%) was obtained when using the SANISAND model compared to that obtained using the hypoplastic model with IS ( ξ soil = 0.51%). Finally, notice that although the magnitude of the soil damping estimated in this study (mainly by the SANISAND model) is within the range (0.17% -1.3%) calculated by some previous studies [START_REF] Kementzetzidis | Geotechnical aspects of offshore wind turbine dynamics from 3D non-linear soil-structure simulations[END_REF][START_REF] Damgaard | Time-varying dynamic properties of offshore wind turbines evaluated by modal testing[END_REF], the forced vibration amplitude prior to unloading significantly influences the magnitude of soil damping [START_REF] Kementzetzidis | Geotechnical aspects of offshore wind turbine dynamics from 3D non-linear soil-structure simulations[END_REF] and thus, rotor-stop tests after a long history of forced vibration are deemed necessary to properly determine the soil damping.

3. Simplified foundation models

This section is devoted to investigate the suitability of two simplified foundations models suggested in literature for the dynamic analysis of OWTs. The two simplified foundation models are based on the well-known Winkler method (called also p-y method). Within the p-y method, the monopile is substituted by an elastic beam and the soil by a set of independent springs (with nonlinear behavior) distributed along the monopile embedded depth. In this paper, the two p-y relations that are investigated are those provided by the American Petroleum Institute (API) [START_REF]Recommended practice for planning, designing and constructing fixed offshore platforms: Working stress design[END_REF] and by Fuentes et al. [START_REF] Fuentes | A p-y Model for Large Diameter Monopiles in Sands Subjected to Lateral Loading under Static and Long-Term Cyclic Conditions[END_REF].

The API method characterizes the p-y relation as a hyperbolic function which depends only on the peak friction angle or the relative density of sand. It is worth to mention here that the widely employed API p-y relations were originally calibrated on tests performed on flexible small-diameter slender piles and thus, lack accuracy when used for large diameter OWT monopiles with diameters D ≥ 4 m and aspect ratio (i.e. length-to-diameter L D ratios) between 4 and 6.

Fuentes et al. [START_REF] Fuentes | A p-y Model for Large Diameter Monopiles in Sands Subjected to Lateral Loading under Static and Long-Term Cyclic Conditions[END_REF] have recently proposed a nonlinear p-y model for the analysis of large diameter monopiles embedded in cohesionless soils and subjected to static or cyclic lateral loading. This model was calibrated by [START_REF] Fuentes | A p-y Model for Large Diameter Monopiles in Sands Subjected to Lateral Loading under Static and Long-Term Cyclic Conditions[END_REF] against a number of 3D FE simulations incorporating the hypoplastic constitutive model for sand. As the rigid behavior of large diameter monopiles may mobilize important shear stresses at the monopile base, [START_REF] Fuentes | A p-y Model for Large Diameter Monopiles in Sands Subjected to Lateral Loading under Static and Long-Term Cyclic Conditions[END_REF] proposed the incorporation within their model of a base shear force at the tip of the monopile in addition to the p-y curves along the monopile length. The distributed nonlinear lateral p-y curves and the nonlinear base shear-displacement curve (called S B -y B curve) were derived by [START_REF] Fuentes | A p-y Model for Large Diameter Monopiles in Sands Subjected to Lateral Loading under Static and Long-Term Cyclic Conditions[END_REF] based on the parameters of the hypoplastic model for sand { φ c , e d 0 , e c 0 , e i 0 ,β, n B } (see Table 3).

Concerning the implementation of the two distributed spring models within Abaqus, a 3D representation of the monopile is adopted in the dynamic simulations. The soil is replaced in each simplified model by a set of distributed spring elements. Indeed, 45 equally spaced p-y curves at 1m interval are placed in each lateral direction along the monopile length (of 45 m) to represent the lateral resistance of the soil. Two additional springs (i.e. a S B -y B curve for each lateral direction) are considered at the base of the monopile to simulate the shear at the monopile tip only in the case of the model developed by Fuentes et al. [START_REF] Fuentes | A p-y Model for Large Diameter Monopiles in Sands Subjected to Lateral Loading under Static and Long-Term Cyclic Conditions[END_REF]. For the case of the API model, the monopile tip is restrained in all directions.

To examine the accuracy of the two simplified foundations models (i.e. the API p-y curves and the p-y and S B -y B curves proposed by [START_REF] Fuentes | A p-y Model for Large Diameter Monopiles in Sands Subjected to Lateral Loading under Static and Long-Term Cyclic Conditions[END_REF]), dynamic simulations of the OWT are performed under LC17. The obtained results are compared with those of the present 3D model that makes use of the hypoplastic constitutive model.

Figures 21a and21b show respectively the time histories and the corresponding PSDs responses of the fore-aft displacement at the tower top for LC17 as obtained using the two simplified foundation models and the present 3D model using hypoplasticity. In these figures, the black curves denote the results obtained using the hypoplastic model, the red curves correspond to the results of the simplified foundation model by [START_REF] Fuentes | A p-y Model for Large Diameter Monopiles in Sands Subjected to Lateral Loading under Static and Long-Term Cyclic Conditions[END_REF] and the blue curves are those for the API model. As shown in Figure 21a, both simplified foundation models underestimate the tower top displacement, with the model proposed by [START_REF] Fuentes | A p-y Model for Large Diameter Monopiles in Sands Subjected to Lateral Loading under Static and Long-Term Cyclic Conditions[END_REF] being better than that of the API when compared to the present 3D model based on hypoplasticity. The peak and RMS reduction ratios calculated with reference to the values given by the hypoplastic model are found to be respectively 15.5% and 9% for the model given by [START_REF] Fuentes | A p-y Model for Large Diameter Monopiles in Sands Subjected to Lateral Loading under Static and Long-Term Cyclic Conditions[END_REF], and 23.2% and 16.1% for the API. Comparing the PSDs in figure 21b, it is obvious that both simplified foundation models exhibit peaks at a higher frequency compared to the hypoplastic model (0.199 Hz for [START_REF] Fuentes | A p-y Model for Large Diameter Monopiles in Sands Subjected to Lateral Loading under Static and Long-Term Cyclic Conditions[END_REF], 0.2052 Hz for the API and 0.1950 Hz for the hypoplastic model). The response obtained using the API model is stiffer than that obtained using [START_REF] Fuentes | A p-y Model for Large Diameter Monopiles in Sands Subjected to Lateral Loading under Static and Long-Term Cyclic Conditions[END_REF].

To evaluate the suitability of the two simplified foundation models for the computation of the monopile responses at mudline, Figures 22 and23 illustrate respectively the fore-aft displacement and out-of-plane rotation time histories of the monopile at mudline under LC17.

From these figures, it is obvious that the average values of the displacement and rotation at mudline are greatly underestimated when using the two simplified foundation models. Also, the peak and RMS reduction ratios with respect to the present 3D model based on hypoplasticity are respectively (32% and 35% for the displacement and 39% and 34% for the rotation) for the model proposed by [START_REF] Fuentes | A p-y Model for Large Diameter Monopiles in Sands Subjected to Lateral Loading under Static and Long-Term Cyclic Conditions[END_REF] and (44.2% and 45.6% for the displacement and 40.5% and 35.2% for the rotation) for the API.

From Figures 212223, one may conclude that the simplified foundation model proposed by [START_REF] Fuentes | A p-y Model for Large Diameter Monopiles in Sands Subjected to Lateral Loading under Static and Long-Term Cyclic Conditions[END_REF] is more accurate than the API model for the structural response at the tower top where a reduction by around 15.5% and 9% of the peak and RMS responses of the fore-aft displacement at the tower top are obtained; however, it underestimates dramatically the foundation responses at mudline, the peak and RMS reduction ratios of the monopile displacement at mudline being respectively 32% and 35%.

Conclusions

In this paper, an explicit three-dimensional finite element model was developed to investigate the dynamic behavior of a large diameter monopile-supported multi-megawatt OWT installed in sandy soil under perfectly drained conditions. The 3D FE model simultaneously considers It should be noted that the present dynamic 3D model in its current state should be handled with care particularly in design applications as the results presented in this paper are not compared with experimental testing. Also, future research is desirable to improve the model reliability in terms of possible pore pressure build-up. Tables Table 1: Properties of the three-bladed upwind DTU 10 MW [START_REF] Alkhoury | A full three-dimensional model for the estimation of the natural frequencies of an offshore wind turbine in sand[END_REF].

Figures

Table 2: Material parameters of the SANISAND model for fine Karlsruhe sand [START_REF] Wichtmann | [END_REF]. Table 3: Material parameters of the Hypoplastic model with IS for fine Karlsruhe sand [START_REF] Wichtmann | [END_REF].

Basic hypoplastic model [START_REF]Recommended practice for planning, designing and constructing fixed offshore platforms: Working stress design[END_REF] Extended hypoplastic model [START_REF] Abaqus | Finite element analysis program[END_REF] φ c

[°]

e i 0 [-] α [-] β [-] R [-] m R [-] m T [-] β R [-] χ [-] 33



  (i) the real geometrical configuration of the OWT superstructure and the monopile foundation, (ii) the monopile-sand interaction, (iii) an advanced sand constitutive model and (iv) a large deformation analysis for the superstructure and monopile foundation. Nonlinear 10-minutes time domain simulations are performed in Abaqus version 6.14 Standard under two stochastically simulated scenarios of environmental loadings. Two advanced sand constitutive models (the critical state elastoplastic SANISAND constitutive model and the hypoplastic model with intergranular strain), recently calibrated on the same Karlsruhe fine sand by[START_REF] Wichtmann | [END_REF], are used in this paper to compare their predictive ability in simulating the OWT structural and geotechnical responses. The use of these two advanced soil constitutive models is expected to provide a reliable prediction of the OWT responses. Indeed, the dynamic behavior of the soil described by these constitutive models was extensively validated in literature at the element level based on laboratory model tests and it is deemed to provide an accurate prediction of the soil response and consequently an accurate prediction of the OWT structural responses. Based on the obtained results, the following conclusions can be drawn:  The two advanced sand constitutive models have a slight influence on the structural responses (displacement and rotation) at the tower top as compared to their influence on the geotechnical responses (displacement/rotation at mudline and strain level in the soil domain). The SANISAND model generally exhibits a stiffer response relative to the Hypoplastic model; The relevance of considering the two advanced soil constitutive models when performing dynamic simulations for an OWT was proven by a comparison of their results with those obtained using the conventional elastic-perfectly plastic Mohr-Coulomb model which resulted in an underestimation of the structural response (by 11.2% for the peak tower top displacement) and geotechnical response (by 15.1% for the peak monopile displacement at mudline); The hypoplastic constitutive model with IS is computationally more efficient than the SANISAND model when performing dynamic simulations of an OWT;  The SANISAND model provides a slightly higher soil damping compared to that obtained when using the hypoplastic model with IS; The large deformation analysis of OWTs was found to be necessary for an accurate prediction of their structural and geotechnical dynamic responses. This may be explained by the large deformation experienced by this type of slender structures with a big rotating mass at their top; The OWT structural and geotechnical dynamic responses strongly depend on the sand relative density. The effect of the sand relative density on the OWT geotechnical response was shown to be more dramatic. A notable increase in the soil straining (by more than five times) is observed as the relative density passes from 80% to 20%;  The two simplified foundation models based on the p-y formulation provide a stiffer response at the tower top and at mudline when compared to that obtained using the present 3D model based on hypoplasticity. Indeed, both models lead to a notable underestimation of the foundation responses at mudline; the peak reduction of the monopile displacement at mudline being 32% for the model by Fuentes et al. (2021) and 44.2% for the API.

Figure 1 :

 1 Figure 1: DTU 10 MW offshore wind turbine.

Figure 2 :

 2 Figure 2: 3D model of the entire soil-foundation-superstructure system [28].

Figure 3 :Figure 4 :Figure 5 :Figure 6 :

 3456 Figure 3: Points of application of the aerodynamic and hydrodynamic loads

Figure 7 :Figure 8 :

 78 Figure 7: Side-to-side tower top displacement as obtained by two constitutive models under LC10 (a) time history and (b) frequency response spectrum

Figure 9 :Figure 10 :Figure 11 :(Figure 12 :Figure 13 :

 910111213 Figure 9: Tower top rotation time history as obtained by two constitutive models for (a) LC10 and (b) LC17

Figure 14 :Figure 15 :Figure 16 :Figure 17 :

 14151617 Figure 14: Comparison of the fore-aft monopile displacement at mudline in [cm] as obtained by the three soil constitutive models LC17

Figure 18 :Figure 19 :

 1819 Figure 18: Fore-aft tower top displacement as obtained using the hypoplastic constitutive model for a nonlinear dynamic analysis performed under LC17 in the absence and presence of large deformation

Figure 20 :Figure 21 :Figure 22 :Figure 23 :

 20212223 Figure 20: Estimation of the global OWT damping based on the decay of the fore-aft displacement time history at the tower top.

Table 4 :

 4 K13 site characteristic parameters from UpWind project.

	.1	4000	0.27	0.677	1.054	1.212	0.14	2.5	10 -4	2.2	1.1	0.1	5.5
	Load Case	Mean wind speed at hub height, U ref [m/s]	Turbulence intensity, I [%]	Significant wave height, H s [m]	Peak spectral period, T

p

[m]