Non-linear hydroelastic response of a monopile foundation in regular waves

This paper presents an original contribution to the coupled hydroelastic simulation of offshore wind turbines foundations. A non-linear coupled simulation tool has been developed to calculate the response of these systems. The hydrodynamic loads are computed by the non-linear potential flow solver WSCN, developed at Centrale Nantes and based on the weak-scatterer theory. The response of the structure is calculated with a finite element analysis tool using an Euler-Bernoulli beam model. The two solvers are tightly coupled. A first case study is presented here. It consists of a uniform monopile foundation in regular waves. Numerical simulations are compared to results obtained with the SIMA software, which uses Morison equation and non-linear wave kinematics. The two simulation tools show a very good agreement in a series of regular waves up to the third-order load harmonic. Higher-order harmonics appear with WSCN in steep waves, which are not computed by SIMA.


Introduction
Offshore wind turbines can be submitted to violent sea-states and hence to very large hydrodynamic loads. The response of these structures should be simulated accurately to optimise their design. The high-order hydrodynamics of bottom-fixed vertical cylinders in waves have been extensively studied, including notably with the FNV model [1] that has recently been adapted for finite water depth [2]. These models can be used to design monopile foundations, for example to study springing and ringing responses. The well-known Morison equation [3] can also be used for slender structures, potentially including non-linear wave kinematics in severe sea-states [4].
This study presents a new non-linear hydro-elastic numerical model WSCN-FEM, including the weak-scatterer theory-based hydrodynamic solver WSCN [5], developed at Centrale Nantes, and a finite element model to compute the response of the structure. The latter uses the Euler-Bernoulli theory to compute the response of beams. Such coupling can account for fully non-linear incident waves and compute the hydrodynamic loads at the exact position of the body. Similar couplings have for example been developed to compute the hydro-elastic response of ships [6,7], showing good agreement with experimental measurements made on a flexible ship-like barge.
A first case study is presented in this paper. It consists of a monopile foundation in regular waves of various wave periods and steepnesses. The simulations with WSCN-FEM are compared with the reference design tool SIMA [8] which uses Morison's equation and non-linear wave kinematics to compute the hydrodynamic loads.

. Hydrodynamics
The hydrodynamic solver used in this study is based on the weak-scatterer hypothesis. It has been developed at Centrale Nantes since 2011 and is called WSCN. It was initially developed for submerged wave energy converters [5], and was successfully improved to account for free-surface piercing and floating bodies [9,10]. In this solver, the exact position of the body can be considered, allowing for integration of the hydrodynamic pressure on the exact wet body surface below the incident and undisturbed free surface elevation. Non-linear wave fields can be explicitly accounted for in this theory. WSCN can use either linear Airy waves or fully non-linear waves given by the stream-function theory of Rienecker-Fenton [11].
Let us consider a body in waves (fixed or floating) as shown in Figure 1. The wet surface of the body is called S B , the free surface elevation is S FS and the submerged external boundary of the domain (i.e. the tank walls and seabed) is called S D . S is the boundary of the domain, defined as S = S B ∪ S FS ∪ S D . The weak-scatterer theory assumes a potential flow (inviscid, incompressible and irrotational), leading to Laplace equation on the velocity potential φ in the whole fluid domain, as written in equation (1).
φ and the free surface elevation η are decomposed into an incident and a perturbed (or scattered) quantity. The weak-scatterer hypothesis is that the perturbed quantities are small compared to the incident velocity potential and free-surface elevation.
This theory is hence weakly non-linear. The wet surface S B , on which the hydrodynamic pressures are integrated, is then the surface of the body that is below η I . The radiation condition implies that the perturbed velocity potential and free surface elevation are zero far from the body, as written in equation (3).
Boundary conditions are then written for the velocity potential. Using the weak-scatterer approximation given in equation (2), the kinematic and dynamic free surface boundary conditions are written at the undisturbed elevation at z = η I using a Taylor expansion and ignoring the terms O(φ P 2 )  3 and O(η P 2 ). The detailed equations are given in [10]. The slip conditions are written on the other boundaries in equation (4), where v is the velocity of the body mesh.
The second Green's identity can be written for the two twice continuously differentiable functions: φ P and the Green function G(M, P) = 1 MP , M and P being two points of the fluid volume. As detailed by Wuillaume [10], the second Green's identity yields where Ω(M) is interpreted as the solid angle seen from the point M. It is given in equation (6).
In equations (5) and (6), the normal vectors on the fluid boundary are oriented inwards. The value of the potential φ P (M) is known on the free surface, and its normal gradient ∂ φ P ∂ n (M) is known on S B , on the seabed and on the numerical tank walls.
Equation (5) can be written as a linear system in the matrix form as given in equation (7), where φ P ,n is the vector containing the values of ∂ φ P ∂ n at each computation node. This equation is called the first Boundary Value Problem (BVP).
G and H are called influence coefficients matrices and represent the influence of every panel on every node. Their values do not depend on the flow state but on the geometry of the body. This equation is solved together with the flow boundary conditions to compute the flow state. The calculation of the loads implies a fluid-structure coupling, as detailed in section 2.2.

Structure dynamics and hydroelastic coupling
A python Finite Element Method (FEM) solver called Beampy has been developed at Centrale Nantes in partnership with NTNU. The solver is based on Euler-Bernoulli beam theory, and has been validated with comparison to analytic results for static and dynamic cases. Beampy has been tightly coupled to WSCN, which means that the structural and hydrodynamic equations are solved at the same time. In the following, the coupling is called WSCN-FEM.
WSCN-FEM has been implemented following an explicit 4th order Runge-Kutta integration scheme. A direct integration of the dynamics in the FEM model would need very small time-steps in an explicit scheme. For sake of stability and efficiency of the coupling, it has hence been chosen to solve the equation of motion using modal superposition.
The structure modes are computed without accounting for the hydrodynamic coupling. The nodal positions x can be related to the modes with the relation between the modal matrix ψ and the modal amplitudes vector y according to: x = ψy.
The equation of motion in time domain is written in equation (8). In equation (8), F ext refers to all external forces applied to the body (weight, mooring), except for hydrodynamic forces that are included in F W SC .
To compute the hydrodynamic loads, the pressure p W SC calculated by WSCN needs to be integrated over the wet surface of the structure, following where n is the normal vector of the body surface, oriented towards the fluid. Following Bernoulli-Lagrange equation, the pressure computed by the WSCN solver takes the form given in equation (10).
At the start of a given time-step, ∂ φ P ∂t is unknown and needs to be computed over the wet surface. All other terms of equation (10) are known from calculations presented in section 2.1. The hydrodynamic force is then written as where F W SC 0 contains the known part of the hydrodynamic force, i.e. the force induced by the distribution of the known part of the pressure p W SC 0 over the wet surface of the platform, as defined in equation (12).
Solving the equation of motion (8) hence involves two unknowns: the accelerations of all nodesü, or the modal accelerationÿ, and the time-differentiation of the scattered velocity potential ∂ φ P ∂t . Instead of implementing an iterative process which can be expensive with respect to CPU time, a strong fluid-structure coupling is implemented to solve both the structural and the hydrodynamic problems at the same time. This approach is called the implicit boundary method [10].
To calculate ∂ φ P ∂t , a second BVP is introduced. The second Green identity, used in section 2.1, is written for the time-differentiation of the velocity potential and its normal gradient at a point M as with the boundary conditions written in equation (14).
In the Neumann boundary condition on the body written in equation (14), both the nodal acceleration x(M) and ∂ 2 φ P ∂ n∂t are unknown. The term q appears when the slip condition on the body of equation (4) is differentiated in time [12]. It is a convection term calculated from known quantities.
To solve the fluid problem together with the structural problem, the equation of motion (8) has to be solved together with equations (13) and (14). The linear system of equations is then constructed and written in equation (15). It includes the Green equation (solved on the whole mesh), the equation of motion and the slip condition on the body (written on the body mesh).
The matrices G and H are the same as in the first BVP, defined in section 2.1. L is the matrix that projects the hydrodynamic loads computed on the hydrodynamic mesh onto the beam structural discretisation. D is the matrix that relates the accelerations of the hydrodynamic nodes on the body to the accelerations of the structural mesh. B is calculated when differentiating the slip condition on the body. Q is the convection terms vector.
It should be noted that the hydrodynamic domain may change at each time step, following the intersection between the incident free-surface elevation around the body geometry. Thus, the influence coefficient matrices G and H, and the projection matrices L and D must be calculated at each time-step. Computing G and H induces a large CPU cost, they are hence kept constant throughout the 4 RK4 substeps over one time-step. A GMRES (generalised minimal residual) iterative scheme is used to solve the linear system of equations.
On the one hand, keeping stable the integration of high-frequency modes of a structure in an explicit time scheme requires a very small time-step. The time-step depends on the number of selected modes, and on the shortest natural frequency. On the other hand, the resolution of the hydrodynamic problem at each time-step needs a high CPU time, and a good convergence of the hydrodynamics can be reached with a much larger time-step (about 30 to 50 time-steps per wave period is enough for regular waves). To improve the efficiency of the tight coupling and the synchronisation of the two solvers at each hydrodynamic resolution, it is then chosen to integrate the dynamics of the structure with a sufficiently small time-step over one hydrodynamic step. This consists in the resolution and integration of the uncoupled equation of motion alone, keeping the hydrodynamic terms constant. When the time-schemes of the two solvers are synchronised, the whole system of coupled equations is solved.
The hydrodynamic mesh does not deform with the flexible motion of the body. Only the rigid motion is followed by the hydrodynamic mesh. The flexible modes may however induce velocities and accelerations to the nodes. Hydrodynamic loads induced by elastic deformations are hence accounted for and tightly coupled with the structural response.

Geometry and fluid domain
This study focuses on the simulation of a monopile-like foundation in regular waves. The monopile is a homogeneous cylinder of L 0 = 100 m length, D = 6 m diameter and ε = 7.5 cm thickness, sketched in figure 2. The steel density is 7800 kg.m -3 , its Young modulus is 210 GPa and its Poisson coefficient is 0.3. The cylinder is assumed fully fixed in the seabed, at 30 m depth. The convergence of the beam discretisation is ensured and the beam is eventually described by 50 beam elements of equal length. The 12 first modes are selected to accurately describe the response of the monopile. The frequency of the first dry (i.e. without hydrodynamic added mass) bending mode is 3.82 rad.s -1 .
The convergence of the hydrodynamic mesh has been ensured, focusing on the convergence of the three first hydrodynamic load harmonics. The selected reference panel size on the body is 0.6 m, yielding approximately 2000 vertices on the body. A symmetry condition is applied at the seabed to ensure the slip condition. The domain radius is taken equal to 3 wave lengths. The total mesh is then composed of approximately 5000 nodes, see figure 3. To avoid wave reflections at the domain boundary, an absorption zone has been implemented on the periphery of the domain. Its length is equal to a wavelength.

Reference model
The WSCN-FEM coupling is here compared to results from SIMA [8], developed at SINTEF.  the hydrodynamic loads. The drag loads are neglected in the Morison equation, and the added mass coefficient is calculated using the MacCamy-Fuchs's flow solution [13]. The wave kinematics follow Stokes 2 nd order model, and the hydrodynamic forces are integrated up to z = η I . Rayleigh stiffness-proportional damping is used in both models, with a 2% damping ratio on the first mode. This may induce more damping on the higher modes, but it should not have an impact on the results presented in this paper.
A semi-analytical solution is also used as reference. The response of the structure is calculated using the 12 first analytical modes shapes. The hydrodynamic loads are calculated using the Morison formula but neglecting the structure velocity and acceleration, and using Airy's linear wave kinematics. This solution is referenced as "semi-analytic" in the following.

Load cases
The first verifications of the coupling and comparisons with SIMA consider linear waves with very small wave amplitudes. These simulations showed a very good agreement between the different models but are not presented in this paper because of lack of space.
Here, a series of regular waves, with different wave steepness, are considered. In WSCN-FEM, simulations are run with a hydrodynamic time-step of 0.1 s and structural sub-time-steps are performed with a step of 0.002 s. In SIMA, the simulations are made with an implicit Newmark-β scheme with a time-step of 0.01 s.
The water depth is 30 m. In WSCN-FEM, the waves are generated with the stream-function solution of Rienecker-Fenton [11] and are described in Table 1. Three wave lengths are chosen. The larger waves, with height H 2 , are chosen so that they approach the limit of validity of the Stokes 2 nd order wave model. The mass coefficient C M used in the SIMA model calculated with the MacCamy and Fuchs theory for each wave period (with a lower limit of C M = 1 applied) is given in the fourth column of table 1. The wave steepness, given as the ratio between wave height and wave length H λ is also given. For T = 3 s, the Keulegan-Carpenter number is below 0.5 for both waves which makes the use of Morison equation for slender bodies questionable.
Simulations are run for 150 s, and only the steady state is considered for the comparison. We can then study the calculated hydrodynamic loads and structure response (tower-top motion and mud-line bending moment).

Results
The results are first presented in time-series. Then, harmonic analysis is performed to better understand the differences between the models. Only selected, representative results are presented due to lack of space.

Time-domain response
The simulated time-series are first compared over a few wave periods. Results are shown for the bending moment around the y axis at the mudline. The results from the coupling WSCN-FEM for two wave periods (T = 5 s and T = 8 s) are compared with those from SIMA and the semi-analytic model and plotted in figures 4 and 5. The agreement is rather good with the least steep waves (figures on the lefthand side): responses are in phase with less than 1.5% relative difference in peak-to-peak variation. In these cases, the difference between SIMA and the semi-analytic model come mostly from non-linearities in the considered wave kinematics. Differences grow when the waves are steeper. The bending moment computed by WSCN-FEM shows higher frequency oscillations. Similar observations have been made on the other calculated quantities, as the tower-top deflection for instance.  frequency is close to 3ω. For T = 8 s, it is close to the 5 th wave harmonic. However, it can be expected that the wet natural frequency of the monopile is slightly lower due to the hydrodynamic added mass. Despite only having 2 nd order wave kinematics, partial 3 rd harmonic effects can also be obtained with the SIMA model. These effects are captured notably because the moment is calculated by integration of the wave loads up to the incident 2 nd order free-surface elevation in the Morison loads calculation. However not all the 3 rd order effects are captured and the harmonic amplitude may be inaccurate.
In these figures, we can see that the mean bending moment (harmonic 0) is much smaller with SIMA than with WSCN-FEM. For example for T = 5 s and H = 1 m, it is equal to 58 kN.m with WSCN-FEM and to 15 N.m with SIMA. The mean value computed by WSCN-FEM corresponds closely with the potential flow theory 2 nd order drift load [3]. As explained in section 3.2, drag loads are not included in SIMA. The viscous drift is hence not computed. This harmonic amplitude is however small compared to the 1 st , 2 nd and 3 rd order harmonics so it has a limited impact on the analysis.
For a clearer understanding, the relative differences between the 1 st to 3 rd harmonics amplitudes are given in table 2. This relative difference is defined for the harmonic amplitude X as δ = X SIMA −X W SCN X SIMA . The agreement on the 1 st harmonic amplitudes is always good (below 4% relative difference). However, the 2 nd harmonic amplitude is larger with SIMA when T = 5 s and lower when T = 8 s. The  agreement on the 2 nd harmonic amplitude between the two solvers is better at T = 5 s. The differences between SIMA and WSCN-FEM grow with the wave steepness. In steep waves, the wave diffraction calculated by MacCamy-Fuchs and used in the Morison equation of SIMA may become inaccurate. As previously mentioned, the accuracy of the 3 rd order effects in the SIMA model is questionable as only part of the third order effects are taken into account. A good agreement on this harmonic between WSCN-FEM and SIMA is hence not expected. It is also observed that the agreement for this harmonic is better in the small wave with T = 8 s, when the bending mode does not match with its frequency.
For the steepest waves, the Stokes 2 nd order wave model almost reaches its limit of validity [14]. The differences induced by the high-order harmonics on the kinematics become noticeable and could increase the differences in the responses of the monopile.
The SIMA model does not capture harmonics above the 3 rd order harmonic. It is interesting to note that higher-order effects are captured by WSCN-FEM and that they impact the response of the monopile. The 5 th harmonic in the response computed by WSCN-FEM is particularly strong with the steepest wave (see figure 7). Here the 5 th harmonics match with the first natural frequency. Note that these results should be compared to experimental measurements or other high-order simulation tools in order to evaluate their accuracy.
Additionally, with the largest wave (T = 8 s and H = 5.02 m), the Keulegan-Carpenter number for the linear Airy wave reaches KC = 2.6. With steep and non-linear waves, the local KC number calculated at the wave crests would increase and flow separation could then occur [2]. Both models, which neglect the viscous effects, would then be inaccurate. For smaller waves however, the wave scattering is expected to be dominant.  10

Conclusion
The new simulation tool WSCN-FEM, coupling non-linear potential flow theory with a finite element method solver, has been developed to simulate the hydro-elastic response of bottom-fixed offshore wind turbines support structures. The hydrodynamic solver is based on the weak-scatterer theory, which can include fully non-linear wave kinematics and computes the hydrodynamic loads at the exact position of the body, accounting for the wet surface up to the incident and undisturbed free surface elevation. The boundary value problem is then solved in time-domain, which induces a high CPU cost, but the solver can then include non-linear hydrodynamic loads on arbitrary geometries.
A first verification test is presented in this paper. It is a monopile foundation in a series of regular waves of various periods and steepnesses. The results from WSCN-FEM are compared with those from the reference design tool SIMA that uses the Morison equation with Stokes second-order wave kinematics to compute the response of the structure. The paper focuses on the mudline bending moment calculated by the solvers. Good agreement is obtained, particularly in small steepness waves. A harmonic analysis shows that the agreement on the first load harmonics are very good, and relatively good with the second harmonics. Differences between the second and higher-order harmonics grow with increasing steepness.
Higher than third order effects cannot be computed accurately using Morison's equation and secondorder wave kinematics. It would thus be interesting to use SIMA with fifth-order wave kinematics within the Morison formula for an additional comparison. Also it would be particularly interesting to compare the results of WSCN-FEM with experimental measurements on flexible monopile foundations.
The present work corresponds to a first step in the development of a tool dedicated to the hydroelastic analysis of floating wind turbines. Further developments are intended to reach this goal.