Hydro-Servo-Aero-Elastic Analysis of Floating Offshore Wind Turbines

Abstract

A fully coupled hydro-servo-aero-elastic simulator for the analysis of floating offshore wind turbines (FOWTs) is presented. All physical aspects are addressed, and the corresponding equations are concurrently solved within the same computational framework, taking into account the wind and wave excitations, the aerodynamic response of the rotor, the hydrodynamic response of the floater, the structural dynamics of the turbine-floater-mooring lines assembly and finally the control system of the wind turbine. The components of the complex multi-physics system of a FOWT interact with each other in an implicitly coupled manner leading to a holistic type of modeling. Different modeling options, of varying fidelity and computational cost, are made available with respect to rotor aerodynamics, hydrodynamic loading of the floater and mooring system dynamics that allow for timely routine certification simulations, but also for computationally intense simulations of less conventional operating states. Structural dynamics is based on nonlinear multibody analysis that allows reproducing the large rigid body motions undergone by the FOWT, as well as large deflections and rotations of the highly flexible blades. The paper includes the description of the main physical models, of the interaction and solution strategy and representative results. Verification is carried out by comparing with other state-of-art tools that participated in the Offshore Code Comparison Collaboration Continuation (OC4) IEA Annex, while the advanced simulation capabilities are demonstrated in the case of half-wake interaction of floating wind turbines by employing the free-wake aerodynamic method.

1. Introduction

Over the last decade, the advantages of offshore wind have increasingly attracted energy supply developers. Offshore sites exhibit substantially higher mean wind speeds and reduced levels of turbulence implying higher power production and lower fatigue loads, while public annoyance is much less than that for onshore installations since visual impact and noise are substantially lower. Over the last decade, the offshore installed capacity in Europe has increased from 2 GW in 2009 to 22 GW in 2019, which, in terms of the total installed wind capacity, corresponds to a change from 2.5% to 11%. It is thus, understood that meeting EU’s target of a 300 GW overall installed capacity by 2030 will greatly rely on offshore wind energy development.

Today, the vast majority of the offshore wind farms are in shallow waters, utilizing support structures that are fixed in the seabed (bottom-fixed). However, as shallow water sites become scarce, the anticipated overwhelming development is expected to exploit sites of high wind potential in deep waters, with floating systems. Floating offshore wind turbines (FOWTs), whose development nowadays is in transition from demonstration to commercialization, face more challenging wind and wave excitations as compared to their bottom-fixed counterparts. Because of that, but also because of the increased cost of floating substructures, the overall cost of wind energy in deep waters is expected to be higher. Therefore, in order to maintain the competitiveness of wind energy, wind turbines need to be up-scaled leading to more delicate designs, with lightweight, flexible rotors. The designs of very big and highly flexible rotors that are mounted on (moving) floating substructures call for dedicated servo-hydro-aero-elastic tools formulated in a holistic, interdisciplinary context. In such an approach, all the components of a FOWT system are simulated simultaneously and are implicitly coupled within the same modeling environment, while all disciplines involved (i.e., aerodynamics, structural dynamics, elasticity, hydrodynamics and control) are modeled using the highest possible level of fidelity. In particular, holistic design tools for FOWT turbines should be able to account for: The complex, stochastic combined wind and wave excitation, the aerodynamics of a rotor undergoing large motions induced by the floating platform, the highly nonlinear structural dynamics of the flexible components (the blades, the shaft, the tower, the floater and the mooring lines), the hydrodynamic loading on complex platform geometries and finally the controls of the turbines. Moreover, the tools should be able to provide design loads for the subcomponents of the system for all possible design and environmental conditions that the turbine encounters during its lifetime. In this regard, the challenges of such all-inclusive modeling are next discussed in more detail also in reference to the current state-of-the art.

Aerodynamics: The modeling challenges include: The unsteady nature of the inflow due to ambient turbulence, wind shear, wind yaw misalignment and the floater induced motions (in particular the pitching and yawing motions); the onset of stall over the rotor blades, a crucial flow feature especially in standstill conditions; the effect of the 3D wake development and finally the nonlinear aeroelastic coupling. This is a quite demanding mix, far more complicated than in other applications, so that one would expect to find in use only advanced aerodynamic models. However, this is not the case. The vast majority of the existing state-of-the-art design and analysis tools are based on the blade element momentum theory (BEMT) [1]. Other more sophisticated models do exist in two main options; those based on free-wake potential methods [2] (medium fidelity) and those based on computational fluid dynamics (CFD) viscous solvers [3] (high fidelity). Both options have been utilized in FOWT hydro-servo-aeroelastic analyses but only to a limited extent, by only few research groups, while their application is mainly restricted to research rather than design work, aiming at a better understanding of the underlying physics and at assessing or calibrating BEMT-based models. It is noted that, although CFD models have been integrated into holistic modeling framework tools, due to their excessive computational cost they have only been employed in simulations of deterministic inflow conditions [3,4], or simplified turbine configurations in which critical flexibilities of the system (e.g., rotor and tower flexibilities) are suppressed, aiming at moderating the computational effort.

With regard to BEMT-based models, there exist various versions that differ in their details, some related to purely implementation aspects [5]. Although there is no complete consensus, still there is good agreement on the quality of predictions obtained with BEMT models when properly calibrated. Semi-empirical models, although not strictly predictive, are tunable and can indeed become very good design tools. Specifically to FOWT, care is needed in the application of the corrections that account for the platform induced dynamic tilt and yaw motions (combined yaw misalignment and flow inclination effects) [6].

The next alternative to BEMT is the free-wake potential models that apply vortex methods [7]. Vortex theory is quite old and is part of the classical aerodynamic theory. The most well-known examples are the lifting-line and lifting surface theories that were developed in the early 1960s for fixed-wing aircraft. Vortex models are 3D by construction, with strong coupling along the span, which is completely absent in BEMT. In vortex models the tip and root corrections implemented in BEMT are no longer needed, while the assumption of infinite number of blades made by BEMT (in the momentum part) is dropped. Additionally, because of the free-wake formulation, the response to dynamic inflow is inherently accounted for. Although vortex methods address most of BEMT shortcomings, they still rely on 2D aerodynamic polars that a-posteriori correct the loads for viscous effects.

CFD viscous flow solvers are regarded as the best choice in terms of completeness and accuracy. The most frequent option is to solve the Reynolds-averaged Navier–Stokes (RANS) equations supplemented with appropriate turbulence closure and transition modeling. Out of the variety of turbulence closures, the most frequently applied are the eddy viscosity models like the k-ω [8,9]. Although still very expensive for aeroelastic simulations, eddy viscosity modeling is by far less expensive than the more advanced closures, detached eddy simulation (DES) or large eddy simulation (LES) [4,10].

Hydrodynamics: Two methods are widely used for estimating the hydrodynamic loading. The first is based on potential hydrodynamic theory and the second on Morison’s empirical equation. Linear potential theory accounts for wave excitation, diffraction and radiation (added mass and added damping) [11] and is valid in case of small to moderate waves. Radiation concerns the six rigid body motions of the floater and should be extended to also include additional elastic modes in case the flexibility of the floater is considered. Potential theory can be extended to second order, in order to include mean drift forces or second order difference and/or sum frequency loads. These loads are proportional to the square of the wave amplitude A and could be of significant importance in case of increased wave steepness (H/λ) at low or high wave frequencies [12]. Often, only the difference frequency loads are considered using Newman’s approximation [13], which only requires the mean drift forces from the solution of the first order hydrodynamic problem [14]. In order to account for the viscous damping, the quadratic damping term appearing in Morison’s equation is usually added on top.

On the other hand, Morison’s semi-empirical equation [15] applies to slender bodies that may be flexible and is valid in case the ratio λ/D > 5, where D is the characteristic length of the body (i.e., the diameter) and λ the wave length. Two hydrodynamic coefficients (inertia and drag coefficients) should be defined. Morison’s equation accounts for wave excitation, diffraction, hydrodynamic added mass and viscous effects, while the linear radiation hydrodynamic damping is not accounted for. Since Morison’s equation is subjected to calibration, a certain level of uncertainty is introduced. At engineering level, if the empirical parameters are fitted based on measured data, then higher order wave theories such as the stream function [16], can be used. Currently, for bottom mounted support structures, the hydrodynamic modeling is exclusively performed by Morison’s equation, while for floating wind turbines both the potential theory and Morison’s equation are applied.

In recent CFD based hydro-aero-elastic simulations, both rotor aerodynamics and hydrodynamic loads on the floater and the mooring lines are computed either simultaneously using the Volume of Fluid method [3], or via a hybrid approach in which the aerodynamic loads are computed using standard grid based Eulerian CFD analysis while hydrodynamic loads are accounted for using a Lagrangian smoothed particle hydrodynamics (SPHs) [4] approach. For the time being, this type of analysis costs a lot and is limited to simplified configurations and/or flow conditions.

Structural modeling: For wind turbines, beam theory is almost exclusively used [17]. Full 3D FEM structural modeling is only applied to certain areas such as the hub and the nacelle, where beam theory fails to provide design information. All early developments considered linear classical beam modeling. Making ad-hoc additions or corrections, certain nonlinearities were included, as, for example, the centrifugal stiffening due to the rotation of the blades. However, the basic modeling remained linear giving relatively stiff designs. In this context, linear models provided sufficiently good predictions, so there was no need to consider more advanced structural modeling. However, as the size of commercial turbines increased, and the pressure for cost reduction became critical, blades started to become more flexible and the validity of linear beam theory was put in doubt. Larger deflections were anticipated leading to significant nonlinear coupling effects and potentially unfavorable blade loading [18].

Therefore, over the last ten years, a number of advanced nonlinear beam models have been implemented in the new generation of aeroelastic design tools. Options in this respect include nonlinear geometrically exact formulations (generalized Timoshenko methods) [19], multi-body formulations [20] and truncation methodologies [21]. Methods in the first category apply an extended form of Hamilton’s variational principle. However, instead of expressing variations in terms of displacement and rotation variables, they adopt an intrinsic formulation of beam theory, see [19]. The formulation is intrinsic if not tied to a specific choice of displacement and rotation variables. Methods in the second category divide each component into a number of interconnected elements that are either considered as (flexible) beams or as rigid bodies. For the flexible beams, all types of models have been used, ranging from simple linear to geometrically exact and nonlinear ones. Finally, the last category comprises methods in which the nonlinear dynamic equations of the deformed beam are derived by following a consistent ordering scheme, whereby higher order terms beyond a certain degree of accuracy are truncated [22]. In deriving the dynamic equations, the assumption of homogenous isotropic material is usually made. In this way, the necessary input is limited to the standard set of structural properties (bending stiffness in flapwise and edgewise direction, torsional stiffness etc.). Recently a number of methods have been developed in which the full stiffness matrix of the composite lay-up is taken into consideration [23]. Methods of this type need, as input, the detailed inner structure of the composite materials in order to derive appropriate equivalent beam structural properties. Such models become necessary as passive load alleviation based on the structural tailoring of the inner structure of the blades [23,24] is gaining attention.

Amongst the existing design tools, some directly apply the finite element method (FEM) to the full dynamic system of equations [18], while others adopt reduced order modeling based on various order reduction techniques, such as the linear modal expansion [1] or the Craig-Bampton method [25]. In the latter case, the aim is to substantially reduce computational effort in view of the long list of simulations required by the IEC standard [26].

Mooring lines: The modeling is defined either in quasi-static or fully dynamic context and discretized using FEM. In every time step, quasi-static models solve the static catenary equations without considering inertial effects or hydrodynamic loading (see for example [27]). On the contrary, dynamic models solve the dynamic equilibrium equations of the mooring lines, as performed for any other flexible part of the wind turbine. Inertia of the mooring lines and gravity load, as well as hydrodynamic loads (by means of Morison’s equation), are usually included. The interaction of the mooring line with the seabed is important and is usually modeled by non-linear springs and dampers along the seabed that prevent the catenary chain to be further submerged.

Hydro-servo-aero-elastic synthetic modeling tools are classified depending on the methods adopted for the modeling of each module (hydro, servo, aero, elastic), as detailed in the previous paragraphs. The hGAST model that is presented in the paper combines different modeling options of varying fidelity for aerodynamic and hydrodynamic analysis. Rotor aerodynamics are handled by either an engineering BEMT model enhanced with various corrections to account for effects that lie beyond the capabilities of momentum theory or a medium fidelity free-wake vortex model. In a similar manner, hydrodynamic loads can be either calculated through the low fidelity Morison’s equation or by means of medium fidelity potential theory. The low fidelity options are mainly used for routine certification and design simulations of mature concepts at conventional operating conditions. The medium fidelity options are selected in case novel concepts or less conventional operational states are addressed and deeper physical insight is required. Selection of the sub-models in such a holistic context balances between accuracy and cost. The least possible loss of fidelity is introduced by the default option, while special care has been taken in reducing the run time when using costly options (e.g., particle mesh acceleration technique in vortex simulations). The structural dynamics modeling employed in the model is the state-of-the-art for this type of analyses. It allows for large rigid body motions and large deflections and rotations of highly flexible components of the wind turbine (e.g., the blades) or of the mooring lines. Finally, in the context of the generalized modeling environment any aerodynamic or hydrodynamic sub-module can be easily plugged into the structural module without cancelling the nonlinear character of the underlying interactions.

The backbone of the solver is based on the multibody methodology implemented in the General Aeroelastic Structural Tool (GAST), developed in the 1990s [28], which over the years went through several revisions and extensions. In its present form, the hydro-servo-aero-elastic tool offers advanced features lying at or beyond the current state-of-the-art employed in most of the existing software used either for academic or industrial purposes. More specifically, geometric non-linear structural and inertial effects due to large deflections and rotations can be captured within the multibody structural dynamic modeling context [18]. Moreover, prebend and/or swept blades can be consistently modeled in the structural dynamics and the aerodynamic modules, respectively [29,30]. In addition, fully populated stiffness properties can be introduced in the beam modeling that account for the anisotropy of the material which is of particular importance when simulating blades with bend twist coupling (BTC) capabilities [23]. Another advantage of the dynamic modeling approach considered in the model is that aeroelastic equations are derived in their linearized formed which allows solving the linear stability problem for the full configuration at any highly deflected state [18]. Both linear eigenvalue [31] and nonlinear [32] aeroelastic stability analysis can be performed, taking into account the control equations in normal operation (closed loop analysis) [33] or in the idling state in storm conditions [34]. Finally, to the best of the authors’ knowledge, the present tool is the only existing in the wind energy sector that is coupled to an efficient, free-wake aerodynamic solver [2,7,30] combined with a multi-block particle mesh technique [35], based on the algorithm presented in [36]. This combination drastically reduces the computational time of the particle-to-particle calculations needed for the evolution of the wake and renders simulations with millions of particles affordable. Through the above, simulation of certification design cases of ten minutes duration are made possible.

In the present work, the fully coupled hydro-servo-aero-elastic simulator for the analysis of floating offshore wind turbines, hGAST, is formulated and the theory of the various interacting modules is presented. Because there are no full-scale measured data available in open literature, proper validation is not possible. This led the Wind Energy community to organize code-to-code comparisons, of which the most widely referenced are in the series of the OC IEA Annexes [37] coordinated by NREL. Though, verification of the code is carried out by comparing predictions against other participating state-of-art tools results. The full configuration is considered with concurrent wind and wave exciting, the controller in closed loop operation and all the structural flexibilities enabled. Then in order to demonstrate the simulation capabilities of the tool, the free-wake vortex aerodynamic method in employed and representative results from the half-wake interaction problem (rotor interacts with the wake of an upstream turbine in partial shading state) for floating wind turbines are discussed. It is worth noting that the simulation of 300 s that considers both rotors was performed on a PC with 6 CPUs in two days.

2. The Coupled Simulation Methodology

Fully coupled hydro-servo-aero-elastic simulations of offshore wind turbines can be formulated within the framework of dynamic systems. The complete multi-physics system is considered as a combination of interacting subcomponents interpreted in a generalized context. The overall concept is illustrated in Figure 1 and includes:

• The turbulent wind inflow that interacts with all parts of the structure above the sea level;

• The incoming sea waves and currents that interact with the submerged fixed or floating support structure;

• The structure of the turbine that interacts on one hand with the air and water flows and on the other with the generator (and eventually the electrical grid) and the control system;

• The structure of the seabed that interacts with the foundation of the mooring lines;

• The control system which specifies the pitch for every blade and the generator torque based on response input and by that interacts with the structure;

• The electrical system which provides the generator’s operating condition.

2.1. Structural (Solid) Dynamics Modeling

The solid parts of a FOWT include the turbine, the floater and the mooring lines as subcomponents. Multibody dynamics is adopted in the formulation; all flexible components are represented as 1D (beam) structures, whereas numerical discretization applies nonlinear FEM. Multibody dynamics can concurrently accommodate rigid body motions and structural flexibility together with the nonlinear inertial and structural effects and geometric couplings that large deflections of highly flexible components may generate (e.g., blades). Rigid motions are imposed by either the control system (e.g., blades’ pitch motion, drive train rotational speed) or by external excitation (e.g., floater motions due to waves and blades rotation due to wind).

2.1.1. Wind Turbine and Floater

The solid components of an offshore wind turbine configuration include the members of the support structure (the buoyant elements and brackets of the floater or struts of the bottom-fixed support), the blades, the drive train and the tower. All the above flexible components, which are usually sufficiently elongated, are approximated as a collection of linear Timoshenko beam structures. Other three dimensional structures, such as the hub, the generator, etc., are considered rigid and modeled as lumped masses to which flexibility may be added in the form of a concentrated spring (e.g., the yaw mechanism).

In multibody dynamics, assembly of the components into the full system is carried out by imposing appropriate kinematic and loading conditions at the connection points. In brief, at any connection point, one of the connected component specifies the position (displacements) and orientation (rotations), while all others contribute their reaction loads [18].

A local coordinate system $[\mathrm{O},\mathrm{x},\mathrm{y},\mathrm{z}]$ (see Figure 2a) is assigned to every component (body) with respect to which local elastic displacements are defined. The local frame of every body (e.g., blades, shaft, tower, floater buoys) is subjected to rigid body motions, each defined by a 6-component vector, including 3 translations and 3 rotations. There are rigid body motions that bodies undergo by themselves and motions that are induced through their connections to other bodies. For the k^th body, let ${\mathsf{q}}_{}^{\mathrm{k}}={\left\{\begin{array}{cc}{\mathsf{q}}_{\hspace{0.17em}\mathrm{t}}^{\mathrm{k}},& {\mathsf{q}}_{\mathrm{r}}^{\mathrm{k}}\end{array}\right\}}^{\mathrm{T}}$ denote the set of time invariant displacements and rotations and time variant body self-motions that define the origin and the orientation of the body local system with respect to the inertial frame in the un-deflected state, and ${\mathsf{q}}_0^{\mathrm{k}}={\left\{\begin{array}{cc}{\mathsf{q}}_{0\mathrm{t}}^{\mathrm{k}},& {\mathsf{q}}_{0\mathrm{r}}^{\mathrm{k}}\end{array}\right\}}^{\mathrm{T}}$ denote the motions due to connections. Rigid translations (denoted by sub-script “t”) will displace the body as a whole to positions denoted by ${\mathsf{R}}_{}^{\mathrm{k}}$ and ${\mathsf{R}}_0^{\mathrm{k}}$, respectively, for the two types of motions, while rigid rotations (denoted by sub-script “r”) lead to rotation matrices ${\mathsf{T}}_{}^{\mathrm{k}}$ and ${\mathsf{T}}_0^{\mathrm{k}}$. Based on the above, the (global) position vector ${\mathsf{r}}_{\mathrm{G}\hspace{0.17em}}^{\mathrm{k}}$ of any arbitrary point P over the deflected k^th body with respect to the global inertia frame $\left[{\mathrm{O}}_{\mathrm{G}},{\mathrm{x}}_{\mathrm{G}},{\mathrm{y}}_{\mathrm{G}},{\mathrm{z}}_{\mathrm{G}}\right]$ (see Figure 2a) is expressed through the (local) position vector of P in the un-deflected state ${\mathsf{r}}_l^{\mathrm{k}}={({\mathrm{x}}_l^{\mathrm{k}},{\mathrm{y}}_l^{\mathrm{k}},{\mathrm{z}}_l^{\mathrm{k}})}^{\mathrm{T}}$ and the vector of the local elastic deflections (displacements and rotations) ${\mathit{u}}_{}^{\mathrm{k}}={\left({\mathrm{u}}^{\mathrm{k}},{\mathrm{v}}^{\mathrm{k}},{\mathrm{w}}^{\mathrm{k}},{\mathsf{\theta }}_{\mathrm{x}}^{\mathrm{k}},{\mathsf{\theta }}_{\mathrm{y}}^{\mathrm{k}},{\mathsf{\theta }}_{\mathrm{z}}^{\mathrm{k}}\right)}^{\mathrm{T}}$ (both defined in body local co-ordinates), and the motions ${\mathsf{q}}_{}^{\mathrm{k}}$ and ${\mathsf{q}}_0^{\mathrm{k}}$ as follows:

(1)${\mathsf{r}}_{\mathrm{G}\hspace{0.17em}}^{\mathrm{k}}={\mathsf{R}}_{}^{\mathrm{k}}({\mathsf{q}}_{}^{\mathrm{k}};\mathrm{t})+{\mathsf{T}}_{}^{\mathrm{k}}({\mathsf{q}}_{\mathrm{r}}^{\mathrm{k}};\mathrm{t})\cdot \left({\mathsf{R}}_0^{\mathrm{k}}({\mathsf{q}}_0^{\mathrm{k}};\mathrm{t})+{\mathsf{T}}_0^{\mathrm{k}}({\mathsf{q}}_{0\mathrm{r}}^{\mathrm{k}};\mathrm{t})\cdot \left({\mathsf{r}}_{l\hspace{0.17em}}^{\mathrm{k}}+\mathbf{S}\cdot {\mathit{u}}^{\mathrm{k}}(\mathrm{t})\right)\right)$

where S matrix is given by:

(2)$\mathit{S}=\left[\begin{array}{cccccc}1& 0& 0& 0& {\mathrm{z}}_l^{\mathrm{k}}& 0\\ 0& 1& 0& -{\mathrm{z}}_l^{\mathrm{k}}& 0& {\mathrm{x}}_l^{\mathrm{k}}\\ 0& 0& 1& 0& -{\mathrm{x}}_l^{\mathrm{k}}& 0\end{array}\right]$

In Equation (1), ${\mathsf{q}}_{}^{\mathrm{k}}$ convert local to global co-ordinates (and vice versa) through ${\mathsf{R}}_{}^{\mathrm{k}}$ and ${\mathsf{T}}_{}^{\mathrm{k}}$ and does not include degrees of freedom of rigid body or elastic motion communicated by other bodies connected with k body. These motions are included in ${\mathsf{q}}_0^{\mathrm{k}}$ and appear in (1) through ${\mathsf{R}}_0^{\mathrm{k}}$ and ${\mathsf{T}}_0^{\mathrm{k}}$. In Figure 2b, ${\mathsf{q}}_0^{\mathrm{k}}$ concerns the motions that are induced on the blade by the end of the drive train in the deflected state, while ${\mathsf{q}}_{}^{\mathrm{k}}$ concerns the fixed offset and rotations of the blade local system with respect to the inertial frame plus the pitch motion of the blade (body self-motion).

Kinematic constraints are applied to the end nodes of connected bodies (P1 and P2 in Figure 2b). For example, node P1 of every blade is rigidly connected to node P2 of the drive train. In this case, the rigid connection kinematic conditions are expressed through the nonlinear velocity/angular velocity equations:

(3)$\begin{array}{l}{\dot{\mathsf{r}}}_{\mathrm{G}\hspace{0.17em}}^{\mathrm{b}\mathrm{l}\mathrm{a}\mathrm{d}\mathrm{e}}(\mathrm{P}1)={\dot{\mathsf{r}}}_{\mathrm{G}\hspace{0.17em}}^{\mathrm{d}\_\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}(\mathrm{P}2)\hspace{1em}\\ {\dot{\mathsf{\theta }}}_{\mathrm{G}}^{\mathrm{b}\mathrm{l}\mathrm{a}\mathrm{d}\mathrm{e}}(\mathrm{P}1)={\dot{\mathsf{\theta }}}_{\mathrm{G}}^{\mathrm{d}\_\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}(\mathrm{P}2)\end{array}$

where ${\dot{\mathsf{\theta }}}_{\mathrm{G}}^{}={\left({\dot{\mathsf{\theta }}}_{\mathrm{G}1}^{},{\dot{\mathsf{\theta }}}_{\mathrm{G}2}^{},{\dot{\mathsf{\theta }}}_{\mathrm{G}3}^{}\right)}^{\mathrm{T}}$ is the angular velocity vector of a body’s end node (P1 or P2) with respect to the inertial frame. While linear velocities ${\dot{\mathsf{r}}}_{\mathrm{G}\hspace{0.17em}}^{}$ are derived through differentiation of (1), angular velocities are obtained through:

(4)${\left({\mathsf{{\rm T}}}^{\mathrm{k}}\cdot {\mathsf{{\rm T}}}_0^{\mathrm{k}}\right)}^{{T}}\cdot \left({\mathsf{{\rm T}}}^{\mathrm{k}}\cdot {\mathsf{{\rm T}}}_0^{\mathrm{k}}\right)=\left(\begin{array}{ccc}0& -{\dot{\mathsf{\theta }}}_{l3}^{\mathrm{k}}& {\dot{\mathsf{\theta }}}_{l2}^{\mathrm{k}}\\ {\dot{\mathsf{\theta }}}_{lz}^{\mathrm{k}}& 0& -{\dot{\mathsf{\theta }}}_{l1}^{\mathrm{k}}\\ -{\dot{\mathsf{\theta }}}_{l2}^{\mathrm{k}}& {\dot{\mathsf{\theta }}}_{l1}^{\mathrm{k}}& 0\end{array}\right)$

where ${\dot{\mathsf{\theta }}}_l^{\mathrm{k}}={\left({\dot{\mathsf{\theta }}}_{l1}^{\mathrm{k}},{\dot{\mathsf{\theta }}}_{l2}^{\mathrm{k}},{\dot{\mathsf{\theta }}}_{l3}^{\mathrm{k}}\right)}^{\mathrm{T}}$ is the angular velocity vector of the k^th body expressed with respect to its local co-ordinate system. Then,

(5)${\dot{\mathsf{\theta }}}_G^{\mathrm{k}}=\left({\mathsf{{\rm T}}}^{\mathrm{k}}\cdot {\mathsf{{\rm T}}}_0^{\mathrm{k}}\right)\cdot {\dot{\mathsf{\theta }}}_l^{\mathrm{k}}$

In addition to the kinematic conditions that are imposed at the connection points, loading conditions must be also satisfied. In particular, at every connection point, one of the connected bodies contributes the displacements and rotations to the others, which in turn contribute their internal (reaction) loads. In the above example, the blade receives displacements and rotations from the drive train end node P2, while the drive train end node receives reaction loads of node P1 of the blade. It is noted that reaction loads of a body can be expressed as functions of the degrees of freedom of the body and; therefore, an implicit coupling between bodies can be readily established.

Of particular importance is that the above multibody formulation is extended to the body level. Highly flexible bodies, such as the blades, are divided into a number of interconnected sub-bodies, each considered as a single or as an assembly of linear beam elements. Large deflections and rotations gradually build up and nonlinear dynamics are introduced by imposing to each sub-body, the deflections and rotations of preceding sub-bodies as rigid body motions. Dynamic coupling of the sub-bodies is again introduced by communicating the reaction loads (three forces and three moments), at the first node of every sub-body to the free node of the previous sub-body, as external load. Illustration of the extension of the multibody framework at the level of a body is illustrated in Figure 3.

Linear Timoshenko beam modeling is applied in order to account for the local deflections vector ${\mathit{u}}_{}^{\mathrm{k}}$ of every flexible body/sub-body. The coordinate system $[{\mathrm{O}}^{\prime },\mathsf{\xi },\mathsf{\eta },\mathsf{\zeta }]$ shown in Figure 4 is the cross section local system. Around the axes of this system, local bending and torsion rotations take place. The set of dynamic equilibrium equations (three force and three moment equations) of the k^th body takes the form:

(6)${\displaystyle \underset{A}{\int }\rho \mathrm{d}\mathrm{A}\hspace{0.17em}\left({\mathsf{S}}^{\mathrm{T}}\cdot \mathsf{S}\right)}\cdot {\left({\mathsf{{\rm T}}}^{\mathrm{k}}\cdot {\mathsf{{\rm T}}}_0^{\mathrm{k}}\right)}^{\mathrm{T}}\cdot {\ddot{\mathsf{r}}}_{\mathrm{G}}^{\mathrm{k}}={\left(\begin{array}{c}{{\mathrm{F}}^{\prime }}_{\mathrm{x}}\\ {{\mathrm{F}}^{\prime }}_{\mathrm{y}}\\ {{\mathrm{F}}^{\prime }}_{\mathrm{z}}\\ {{\mathrm{M}}^{\prime }}_{\mathrm{x}}+{\mathrm{F}}_{\mathrm{z}}-{\mathrm{F}}_{\mathrm{y}}{{\mathrm{w}}^{\prime }}^{\mathrm{k}}\\ {{\mathrm{M}}^{\prime }}_{\mathrm{y}}\\ {{\mathrm{M}}^{\prime }}_{\mathrm{z}}-{\mathrm{F}}_{\mathrm{x}}+{\mathrm{F}}_{\mathrm{y}}{{\mathrm{u}}^{\prime }}^{\mathrm{k}}\end{array}\right)}^{\mathrm{k}}\hspace{1em}+\hspace{1em}\begin{array}{c}\mathrm{e}\mathrm{x}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{a}\mathrm{l}\hspace{0.17em}\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d}\mathrm{s}\hspace{0.17em}\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{l}\mathrm{u}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{g}\\ \mathrm{r}\mathrm{e}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\hspace{0.17em}\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d}\mathrm{s}\mathrm{\hspace{0.17em}}\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{m}\mathrm{u}\mathrm{n}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{d}\\ \mathrm{b}\mathrm{y}\hspace{0.17em}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{n}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{d}\hspace{0.17em}\mathrm{b}\mathrm{o}\mathrm{d}\mathrm{i}\mathrm{e}\mathrm{s}\end{array}\hspace{0.17em}$

where the prime symbol (‘) denotes derivatives with respect to the beamwise local y direction. The terms ${\mathrm{F}}_{\mathrm{y}}{{\mathrm{w}}^{\prime }}^{\mathrm{k}}$ and ${\mathrm{F}}_y{{\mathrm{u}}^{\prime }}^{\mathrm{k}}$ in the moment x and z equations are the only non-linear terms considered in the analysis. This is because they are expected to contribute significantly, especially in the case of rotating beams in which axial force increases due to the centrifugal effect. The two terms give rise to virtual stiffening of the beam as rotational speed increases.

The constitutive relation between the generalized structural loads over a cross section of a beam and the strains-curvatures is given by:

(7)${\mathsf{F}}^{\mathrm{k}}={\mathsf{K}}^{\mathrm{k}}\cdot {\mathsf{\epsilon }}^{\mathrm{k}}\hspace{0.17em},\hspace{1em}{\mathsf{F}}^{\mathrm{k}}={\left({\mathrm{F}}_{\mathrm{x}}^{\mathrm{k}},{\mathrm{F}}_{\mathrm{y}}^{\mathrm{k}},{\mathrm{F}}_{\mathrm{z}}^{\mathrm{k}},{\mathrm{M}}_{\mathrm{x}}^{\mathrm{k}},{\mathrm{M}}_{\mathrm{y}}^{\mathrm{k}},{\mathrm{M}}_{\mathrm{z}}^{\mathrm{k}}\right)}^{\mathrm{T}},\hspace{1em}{\mathsf{\epsilon }}^{\mathrm{k}}={\left({\mathsf{\gamma }}_{\mathrm{x}}^{\mathrm{k}},{\mathsf{\epsilon }}_{\mathrm{y}}^{\mathrm{k}},{\mathsf{\gamma }}_{\mathrm{z}}^{\mathrm{k}},{\mathrm{k}}_{\mathrm{x}}^{\mathrm{k}},{\mathrm{k}}_{\mathrm{y}}^{\mathrm{k}},{\mathrm{k}}_{\mathrm{z}}^{\mathrm{k}}\right)}^{\mathrm{T}}$

where,

(8)${\mathsf{K}}^{\mathrm{k}}={\left\{\begin{array}{cccccc}{\mathsf{K}}_{\mathrm{x}\mathrm{x}}^{\mathrm{A}}& {\mathsf{K}}_{\mathrm{x}\mathrm{y}}^{\mathrm{A}}& {\mathsf{K}}_{\mathrm{x}\mathrm{z}}^{\mathrm{A}}& {\mathsf{K}}_{\mathrm{x}\mathrm{x}}^{\mathrm{B}}& {\mathsf{K}}_{\mathrm{x}\mathrm{y}}^{\mathrm{B}}& {\mathsf{K}}_{\mathrm{x}\mathrm{z}}^{\mathrm{B}}\\ & {\mathsf{K}}_{\mathrm{y}\mathrm{y}}^{\mathrm{A}}& {\mathsf{K}}_{\mathrm{y}\mathrm{z}}^{\mathrm{A}}& {\mathsf{K}}_{\mathrm{y}\mathrm{x}}^{\mathrm{B}}& {\mathsf{K}}_{\mathrm{y}\mathrm{y}}^{\mathrm{B}}& {\mathsf{K}}_{\mathrm{y}\mathrm{z}}^{\mathrm{B}}\\ & & {\mathsf{K}}_{\mathrm{z}\mathrm{z}}^{\mathrm{A}}& {\mathsf{K}}_{\mathrm{z}\mathrm{x}}^{\mathrm{B}}& {\mathsf{K}}_{\mathrm{z}\mathrm{y}}^{\mathrm{B}}& {\mathsf{K}}_{\mathrm{z}\mathrm{z}}^{\mathrm{B}}\\ & \mathrm{s}\mathrm{y}\mathrm{m}.& & {\mathrm{k}}_{\mathrm{x}\mathrm{x}}^{\mathrm{C}}& {\mathrm{k}}_{\mathrm{x}\mathrm{y}}^{\mathrm{C}}& {\mathrm{k}}_{\mathrm{x}\mathrm{z}}^{\mathrm{C}}\\ & & & & {\mathrm{k}}_{\mathrm{y}\mathrm{y}}^{\mathrm{C}}& {\mathrm{k}}_{\mathrm{y}\mathrm{z}}^{\mathrm{C}}\\ & & & & & {\mathrm{k}}_{\mathrm{z}\mathrm{z}}^{\mathrm{C}}\end{array}\right\}}^{\mathrm{k}}$

and,

(9)${\mathsf{\epsilon }}_{\mathrm{y}}={\mathrm{v}}^{\prime },\hspace{0.17em}\hspace{0.17em}{\mathsf{\gamma }}_{\mathrm{x}}={\mathrm{u}}^{\prime }+{\mathsf{\theta }}_z,\hspace{0.17em}\hspace{0.17em}{\mathsf{\gamma }}_{\mathrm{z}}={\mathrm{w}}^{\prime }-{\mathsf{\theta }}_{\mathrm{x}},\hspace{0.17em}\hspace{0.17em}{\mathrm{k}}_{\mathrm{x}}={{\mathsf{\theta }}^{\prime }}_{\mathrm{x}},\hspace{0.17em}\hspace{0.17em}{\mathrm{k}}_{\mathrm{y}}={{\mathsf{\theta }}^{\prime }}_{\mathrm{y}},\hspace{0.17em}\hspace{0.17em}{\mathrm{k}}_{\mathrm{z}}={{\mathsf{\theta }}^{\prime }}_{\mathrm{z}}\hspace{1em}\mathrm{f}\mathrm{o}\mathrm{r}\hspace{0.17em}\mathrm{t}\mathrm{h}\mathrm{e}\hspace{0.17em}\hspace{0.17em}\mathrm{k}-\mathrm{t}\mathrm{h}\hspace{0.17em}\hspace{0.17em}\mathrm{B}\mathrm{o}\mathrm{d}\mathrm{y}$

In (8), ${\mathsf{K}}^{\mathrm{k}}$ is the sectional Timoshenko full stiffness matrix. With respect to a standard Timoshenko beam approach, the elements ${\mathrm{K}}_{\mathrm{x}\mathrm{x}}^{\mathrm{A}}$ and ${\mathrm{K}}_{\mathrm{z}\mathrm{z}}^{\mathrm{A}}$ represent transverse shear rigidity, ${\mathrm{K}}_{\mathrm{y}\mathrm{y}}^{\mathrm{A}}$ axial stiffness, ${\mathrm{K}}_{\mathrm{x}\mathrm{x}}^{\mathrm{C}}$ and ${\mathrm{K}}_{\mathrm{z}\mathrm{z}}^{\mathrm{C}}$ flexural stiffness in flap- and edge-wise directions, respectively, and ${\mathrm{K}}_{\mathrm{y}\mathrm{y}}^{\mathrm{C}}$ torsional stiffness. The off-diagonal elements ${\mathrm{K}}_{\mathrm{x}\mathrm{y}}^{\mathrm{C}}$ and ${\mathrm{K}}_{\mathrm{y}\mathrm{z}}^{\mathrm{C}}$ are responsible for the activation of bend-twist coupling. By introducing (9) into (7), ${\mathsf{K}}^{\mathrm{k}}$ can be split as follows:

(10)${\mathsf{F}}^{\mathrm{k}}={\mathsf{K}}^{\mathrm{k}}\cdot {\mathsf{\epsilon }}^{\mathrm{k}}\hspace{0.17em}={\mathsf{K}}_1^{\mathrm{k}}{{\mathit{u}}^{\prime }}^{\mathrm{k}}+{\mathsf{K}}_2^{\mathrm{k}}{\mathit{u}}^{\mathrm{k}}$

By substituting (10) into (6), the dynamic equations are expressed with respect to the considered degrees of freedom in the form:

(11)${\displaystyle \underset{\mathrm{A}}{\int }\rho \mathrm{d}\mathrm{A}\hspace{0.17em}\left({\mathsf{S}}^{\mathrm{T}}\cdot \mathsf{S}\right)}\cdot {\left({\mathsf{{\rm T}}}^{\mathrm{k}}\cdot {\mathsf{{\rm T}}}_0^{\mathrm{k}}\right)}^{\mathrm{T}}\cdot {\ddot{\mathsf{r}}}_{\mathrm{G}}^{\mathrm{k}}={\left({\mathsf{K}}_1^{\mathrm{k}}{{\mathit{u}}^{\prime }}^{\mathrm{k}}\right)}^{\prime }+{\left({\mathsf{K}}_2^{\mathrm{k}}{\mathit{u}}^{\mathrm{k}}\right)}^{\prime }+\left({\mathsf{K}}_3^{\mathrm{k}}{{\mathit{u}}^{\prime }}^{\mathrm{k}}\right)+\left({\mathsf{K}}_4^{\mathrm{k}}{\mathit{u}}^{\mathrm{k}}\right)\hspace{1em}+\hspace{1em}\mathrm{e}\mathrm{x}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{a}\mathrm{l}\hspace{0.17em}\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d}\mathrm{s}$

In order to apply the FEM, the above set is written in weak form by applying the principle of virtual work and then discretized using first-order shape functions for extension and torsion, and modified C¹ Hermitian functions are used for the two bending displacements that prevent shear locking effect [38] by satisfying static equilibrium compatibility relations.

The discrete FEM Equation (11) and the constraint Equation (3) are analytically linearized and expressed in perturbed form. The advantage of the proposed formulation, in comparison to other multibody formulations applying the Lagrange multipliers approach, is that linear eigenvalue stability analysis can be performed in the resulting linearized dynamic system with respect to a highly deflected steady or periodic state. Unknowns of the structural dynamic problem are the ${\mathsf{q}}_{}^{\mathrm{k}}$ and ${\mathsf{q}}_0^{\mathrm{k}}$ vectors, defined through (i) kinematic constraint equations (like those given by Equation (3)) or (ii) proper control and dynamic equations (e.g., pitch control equation, drive train equation and floater dynamic equations) and the vector of the body local discrete elastic degrees of freedom ${\mathit{u}}_{}^{\mathrm{k}}={\left({\mathrm{u}}^{\mathrm{k}},{\mathrm{v}}^{\mathrm{k}},{\mathrm{w}}^{\mathrm{k}},{\mathsf{\theta }}_{\mathrm{x}}^{\mathrm{k}},{\mathsf{\theta }}_{\mathrm{y}}^{\mathrm{k}},{\mathsf{\theta }}_z^{\mathrm{k}}\right)}^{\mathrm{T}}$.

2.1.2. Mooring Lines

Modeling of the mooring line is based on 1D nonlinear FEM. The mooring lines are modeled as space truss elements that only transfer tensile forces while their compression stiffness is zero [39]. Each element is subjected to inertial, gravitational and hydrodynamic loading, calculated by means of Morison’s equation. The seabed interaction is modeled by appropriate stiffness and damping terms. Coupling of the mooring lines with the floater is established by imposing as kinematic constraint conditions the motions of the floater to the mooring line node attached to the floater, while the mooring line returns reaction loads to the floater. Both catenary mooring lines and tendons for tension leg platforms (TLP) can be modeled.

In Figure 5, a truss element ${\mathrm{P}}_0{\mathrm{Q}}_0$ of initial undeformed length ${\ell }_0=2{\mathsf{\alpha }}_0$ is shown. The non-dimensional coordinate $\mathsf{\xi }$ (ranging from −1 to 1) is used to define the arbitrary point ${\mathrm{R}}_0$ along the element. The position vector of this point with respect to the inertial frame is initially ${\mathrm{r}}_0$. As a result of rigid body motion and elastic deformation, the truss element will move to a new position ${\mathrm{P}}_{\mathrm{n}}{\mathrm{Q}}_{\mathrm{n}}$. The arbitrary point along the element (point ${\mathrm{R}}_{\mathrm{n}}$) will now have a position vector ${\mathrm{r}}_{\mathrm{n}}$, given by the following expression:

(12)${\mathsf{r}}_{\mathrm{n}}={\mathsf{r}}_0+\mathit{u}$

where vector $\mathit{u}$ represents displacement due to elastic deflection and rigid body motion of the element.

By applying the virtual work principle, the dynamic equations of a truss element in space are obtained. The principle requires that, for an arbitrary virtual displacement field, $\delta \mathit{u}$:

(13)${\displaystyle \underset{{\ell }_0}{\int }\delta {\mathit{u}}^{{\rm T}}{\mathsf{q}}^{\mathrm{i}}}\hspace{0.17em}\hspace{0.17em}\mathrm{d}{\ell }_{}=-{\displaystyle \underset{{\mathrm{V}}_0}{\int }\sigma \hspace{0.17em}\delta \epsilon \hspace{0.17em}}\mathrm{d}{\mathrm{V}}_{}+{\displaystyle \underset{{\ell }_0}{\int }\delta {\mathit{u}}^{{\rm T}}{\mathsf{q}}^{\mathsf{e}}}\hspace{0.17em}\hspace{0.17em}\mathrm{d}{\ell }_{}$

where ${\mathsf{q}}^{\mathrm{i}}$ and ${\mathsf{q}}^{\mathsf{e}}$ are the inertial and external (hydrodynamic, gravitational and buoyancy) per unit length forces acting on the element, $\sigma$ and $\delta \epsilon$ are the axial stress and the incremental variation of the axial strain caused by the relative motion of the element ends $P$ and $Q$, and ${\mathrm{V}}_0$ the volume of the undeformed element. The standard definition of the engineering strain [40],

(14)$\epsilon =\frac{{\ell }_{\mathrm{n}}^{}-{\ell }_0^{}}{{\ell }_0^{}}$

is used in order to define the variation $\delta \epsilon$. By considering that the displacements of the two end nodes, $\mathrm{P}$ and $\mathrm{Q}$, of the element are the degrees of freedom of the problem, the displacement of any intermediate point along the element is written as

(15)$\mathit{u}=\mathbf{{\rm N}}\cdot \widehat{\mathit{u}}\hspace{0.17em},\hspace{1em}\hspace{0.17em}\widehat{\mathit{u}}={\left\{\begin{array}{cc}\underset{\mathrm{n}\mathrm{o}\mathrm{d}\mathrm{e}\hspace{0.17em}\mathrm{P}}{\underbrace{{\mathrm{u}}_1\hspace{0.17em}\hspace{0.17em}{\mathrm{v}}_1\hspace{0.17em}\hspace{0.17em}{\mathrm{w}}_1}}& \underset{\mathrm{n}\mathrm{o}\mathrm{d}\mathrm{e}\hspace{0.17em}\mathrm{Q}}{\underbrace{{\mathrm{u}}_2\hspace{0.17em}\hspace{0.17em}{\mathrm{v}}_2\hspace{0.17em}\hspace{0.17em}{\mathrm{w}}_2}}\end{array}\right\}}^{\mathrm{T}}$

where $\mathbf{{\rm N}}$ is the matrix of the shape functions considered linear in the present analysis.

Using the same linear interpolation functions, the coordinates of the arbitrary point ${\mathrm{R}}_0$ along the undeformed element can be defined in terms of the end points ${\mathrm{P}}_0$ $\left({\mathrm{x}}_1\hspace{0.17em},\hspace{0.17em}{\mathrm{y}}_1\hspace{0.17em},\hspace{0.17em}{\mathrm{z}}_1\right)$ and ${\mathrm{Q}}_0$ $\left({\mathrm{x}}_2\hspace{0.17em},\hspace{0.17em}{\mathrm{y}}_2\hspace{0.17em},\hspace{0.17em}{\mathrm{z}}_2\right)$:

(16)${\mathsf{r}}_0^{}=\mathsf{{\rm N}}\cdot \widehat{\mathsf{x}}\hspace{0.17em},\hspace{1em}\hspace{0.17em}\widehat{\mathsf{x}}={\left\{\begin{array}{cccccc}{\mathrm{x}}_1& {\mathrm{y}}_1& {\mathrm{z}}_1& {\mathrm{x}}_2& {\mathrm{y}}_2& {\mathrm{z}}_2\end{array}\right\}}^{\mathrm{T}}$

Using (14), (15) and (16), the incremental variation of the strain $\delta \epsilon$ can be expressed in terms of the virtual displacement $\delta \widehat{\mathit{u}}$ of the element end nodes:

(17)$\delta \epsilon =\frac{1}{{\ell }_0}{t}^{\mathrm{T}}\cdot A\cdot \mathsf{\delta }\widehat{\mathit{u}}={\mathrm{b}}^T\cdot \mathsf{\delta }\widehat{\mathit{u}}$

where t is the unit vector tangent to the element and $A=[-{\mathrm{I}}_{3\times 3},{\mathrm{I}}_{3\times 3}]$ and ${\mathrm{I}}_{3\times 3}$ the 3 × 3 unit matrix.

The stress field is determined through application of Hooke’s laws by making the assumption of linear elastic material. As already stated, the Young modules ${\rm E}$ of the mooring line will be different in compression and extension. The compression value will be equal to zero but in practical terms is taken very small, for numerical stability.

By introducing (17) into (13), and after integration over the cross section area ${\mathrm{A}}_0$, the following system of non-linear dynamic equations is obtained:

(18)${\displaystyle \underset{-1}{\overset{1}{\int }}{\mathsf{N}}^{\mathrm{T}}{\mathsf{q}}^{\mathsf{e}}\hspace{0.17em}{\alpha }_0\hspace{0.17em}\mathrm{d}\xi }-{\displaystyle \underset{-1}{\overset{1}{\int }}{\mathsf{N}}^{\mathrm{T}}{\left(\rho {\mathrm{A}}_0\right)}_{}\hspace{0.17em}{\ddot{\mathsf{r}}}_{\mathrm{n}}\hspace{0.17em}{\alpha }_0\hspace{0.17em}\mathrm{d}\xi }={\displaystyle \underset{-1}{\overset{1}{\int }}\mathsf{b}\hspace{0.17em}\hspace{0.17em}\sigma \hspace{0.17em}{\mathrm{A}}_0\hspace{0.17em}{\alpha }_0\hspace{0.17em}\mathrm{d}\xi }=2{\alpha }_0{\mathrm{A}}_0\hspace{0.17em}\mathsf{b}\hspace{0.17em}\sigma =\mathsf{f}$

where $\left(\rho {\mathrm{A}}_0\right)$ is the linear mass distribution of the mooring line. The non-linearity of Equation (18) lies in two distinct reasons. The first is the inherent capability of the model to accommodate both large displacements of the truss elements and deformations (geometric non-linearity). The second is due to the different stiffness characteristics of the elements in tension and compression (material non-linearity).

The non-linear Equation (18) that concerns one element are assembled for all mooring elements. Then they are expressed in linear perturbed form (linearization about a reference deflected state) and added together with the rest of the structural equations. This final step is done by imposing connection conditions at the attachment points of the mooring lines on the floater. As shown in Figure 6, the floater defines the motions of the mooring line at the attachment point, while the mooring line gives back the tensile reaction force.

2.2. Aerodynamic Modeling

There are two options available for the modeling of the rotor aerodynamics: a low-cost and semi empirical enhanced model based on blade element momentum theory (BEMT) and a medium-cost and fidelity model based on free-wake vortex (VRX) modeling. The first is mainly used for certification routine simulations, while the second is mainly selected when less conventional conditions with strong 3D effects, that lie out of the modeling capabilities of BEMT, are dealt with (e.g., high yaw misalignment, large floater motions, partial wake operation).

2.2.1. Blade Element Momentum Theory Modeling

The standard BEMT model is enhanced in three aspects: (a) a wake-skewed correction model that accounts for yaw and tilt misalignment [41], (b) a dynamic inflow model based on cylindrical wake modeling [42], (c) unsteady aerodynamics and dynamic stall modeling based on the ONERA [43] or the Beddoes–Leishman [44] engineering models. More details on the enhanced BEMT implemented in the model can be found in [45].

2.2.2. Free-Wake Model

The model considers the unsteady flow of an incompressible and inviscid fluid around a multi-component configuration of bodies that may move independently from each other. Thus, in addition to the isolated turbine problem, the case of more than one wind turbine in a cluster can be dealt with.

As illustrated in Figure 7, aerodynamic bodies (e.g., rotating blades) can be considered as (i) lifting lines, (ii) thin lifting surfaces and (iii) thick panel bodies. All aerodynamic bodies are divided into a number of “cross sections”, which are perpendicular to the body reference line. Furthermore, a number of nodal points are defined over every cross section. The wing element between two consecutive sections forms the “strip”. The strip is the main aerodynamic entity of the body, on which the local aerodynamic loads are calculated. In the lifting line analysis, the building block of the numerical grid is the strip. Each section has only two nodes describing the leading and trailing edge of the airfoil, which; therefore, collapses into a single flat plate element (panel) formed upon every strip, as shown in Figure 7. In the lifting surface and thick panel analysis, the building block of the numerical grid is the surface panel element formed by four grid nodes (two pairs belonging to consecutive sections) in clockwise direction. In this case, a number of surface elements (panels) form the strip, as shown in Figure 7. Within every panel, a single control point (collocation point) is defined, on which the circulation, the pressure jump or the local pressure is calculated for lifting line bodies, thin bodies or thick bodies, respectively.

The theoretical backbone of the method is Helmholtz’s decomposition theorem [46], according to which the flow field $\overrightarrow{\mathrm{u}}$ at a point in space $\overrightarrow{\mathrm{x}}$ can be decomposed as:

(19)$\overrightarrow{\mathrm{u}}(\overrightarrow{\mathrm{x}};\mathrm{t})={\overrightarrow{\mathrm{u}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}(\overrightarrow{\mathrm{x}};\mathrm{t})+{\overrightarrow{\mathrm{u}}}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{i}\mathrm{d}}(\overrightarrow{\mathrm{x}};\mathrm{t})+{\overrightarrow{\mathrm{u}}}_{\mathrm{w}\mathrm{a}\mathrm{k}\mathrm{e}}(\overrightarrow{\mathrm{x}};\mathrm{t})$

where ${\overrightarrow{\mathrm{u}}}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{i}\mathrm{d}}=\mathsf{\nabla }\varphi$ is the potential part of the flow associated with the presence of the solid boundaries (e.g., blades) in the flow, ${\overrightarrow{\mathrm{u}}}_{\mathrm{w}\mathrm{a}\mathrm{k}\mathrm{e}}=\nabla \times \overrightarrow{\psi }$ is the vortical part associated with the wake or the free vorticity given by $\overrightarrow{\omega }$ and ${\overrightarrow{\mathrm{u}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}$ is a given external flow field. In the above $\varphi ,\hspace{0.17em}\overrightarrow{\psi }$ denote the scalar and vector potential of the flow, respectively, that satisfy the following field equations:

(20)${\nabla }^2\varphi =\nabla \cdot \overrightarrow{\mathrm{u}}=0, {\nabla }^2\overrightarrow{\psi }=-\nabla \times \overrightarrow{\mathrm{u}}\equiv -\overrightarrow{\omega }$

According to Green’s theorem, ${\overrightarrow{\mathrm{u}}}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{i}\mathrm{d}}$ can be expressed through surface singularity distributions. In the lifting line approach, the body contribution is considered through open “horseshoe” vortex filaments distributed over the strips. In the lifting surface and thick panel options, the body contribution is considered through surface dipole and surface dipole plus source distributions respectively, over the elements. It is noted that for any dipole surface distribution an equivalent line and surface vorticity distribution can be defined. Therefore, a general expression of ${\overrightarrow{\mathrm{u}}}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{i}\mathrm{d}}$ is given by:

(21)${\overrightarrow{\mathrm{u}}}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{i}\mathrm{d}}(\overrightarrow{\mathrm{x}};\mathrm{t})={\displaystyle \underset{\mathrm{S}(\mathrm{t})}{\int }\frac{(\sigma \cdot +\overrightarrow{\gamma }\times )\overrightarrow{\mathrm{r}}}{4\pi {\overrightarrow{\mathrm{r}}}^3}}\hspace{0.17em}\mathrm{d}\mathrm{S}(\overrightarrow{\mathrm{y}})$

where, $\overrightarrow{\mathrm{r}}=\overrightarrow{\mathrm{x}}-\overrightarrow{\mathrm{y}}$, $\mathrm{S}(\mathrm{t})$ denotes collectively the solid boundaries and $\sigma ,\hspace{0.17em}\overrightarrow{\gamma }$ denote the surface source and vorticity distributions. Similarly according to Green’s theorem, ${\overrightarrow{\mathrm{u}}}_{\mathrm{w}\mathrm{a}\mathrm{k}\mathrm{e}}$ can be expressed in the following integral form,

(22)${\overrightarrow{\mathrm{u}}}_{\mathrm{w}\mathrm{a}\mathrm{k}\mathrm{e}}(\overrightarrow{\mathrm{x}};\mathrm{t})={\displaystyle \underset{{\mathrm{D}}_{\omega }(\mathrm{t})}{\int }\frac{\overrightarrow{\omega }(\overrightarrow{\mathrm{y}};\mathrm{t})\times \overrightarrow{\mathrm{r}}}{4\pi {\overrightarrow{\mathrm{r}}}^3}}\hspace{0.17em}\mathrm{d}\mathrm{D}(\overrightarrow{\mathrm{y}})$

where ${\mathrm{D}}_{\omega }(\mathrm{t})$ denotes the region covered by the wake.

In the lifting line option, the bound line vorticity (equivalent to circulation around filaments), $\overrightarrow{\Gamma }$, of every strip is obtained from the solution of the discretized system of the nonlinear lifting line equations:

(23)$\overrightarrow{\mathrm{L}}({\overrightarrow{\mathrm{x}}}_0;\mathrm{t})=\rho \overrightarrow{u}({\overrightarrow{\mathrm{x}}}_0;\mathrm{t})\times \overrightarrow{\Gamma }({\overrightarrow{\mathrm{x}}}_0;\mathrm{t})$

satisfied on every control point ${\overrightarrow{\mathrm{x}}}_0$ of the lifting line grid (i.e., the single control point per strip located at the quarter chord line in the midspan position of the strip). This solution uses as input the airfoil 2D aerodynamic polars at spanwise sections, while corrections to account for viscous and dynamic stall effects are applied based on the ONERA [43] dynamic stall model.

In the thin lifting surface and thick panel body options, the non-penetration, kinematic boundary condition is satisfied at the control points. In this case, the central point of every surface panel element. Furthermore, the Kutta condition is satisfied along the trailing edge (TE) line of all lifting bodies. In the above two modeling options, airfoil polars are not needed, except when applying viscous corrections. As in the lifting line case, this is done on an a-posteriori basis, using the ONERA dynamic stall model.

In potential theory, wakes are introduced as material surfaces carrying surface vorticity $\overrightarrow{\gamma }$. Existing wake models differ on the choice of elements to describe the wake, but most importantly on whether they assume connectivity among these elements. Connectivity is necessary in order to satisfy the fundamental requirement that vorticity is div free or else that the vorticity lines in the wake are either closed or they start and end on a solid boundary or at infinity. These requirements are, by construction, fulfilled when the wake is formed by means of vortex filaments or is retained as a surface. Connectivity; however, can generate numerical problems when the wake is excessively deformed or interacts with solid bodies, as in the present case. In order to bypass this problem, the present model uses freely moving vortex blobs. They are 3D point vortices equipped with core. Vortex blobs are generated at every time step in a two-step procedure. First the wake is released in the form of surface vorticity. At this step, the no-penetration and Kutta conditions are satisfied. Then the convection step is carried out, during which surface vorticity is integrated and transformed into vortex blobs defined by their intensities ${\overrightarrow{\Omega }}_{\mathrm{P}}$, positions ${\overrightarrow{\mathrm{Z}}}_{\mathrm{P}}$ and core sizes $\epsilon$, as shown in Figure 8. The free vorticity is defined as follows:

(24)$\overrightarrow{\omega }(\overrightarrow{\mathrm{x}};\mathrm{t})\cong {\displaystyle \sum _{\mathrm{P}}{\overrightarrow{\Omega }}_{\mathrm{P}}(\mathrm{t})\hspace{0.17em}{\zeta }_{\epsilon }(\overrightarrow{\mathrm{x}}-{\overrightarrow{\mathrm{Z}}}_{\mathrm{P}}(\mathrm{t}))}$

where ${\zeta }_{\epsilon }$ is the cut-off or distribution function defined over the core. In the current implementation, a cubic exponential function [47] is used, which results in the following velocity representation:

(25)${\overrightarrow{\mathrm{u}}}_{\mathrm{w}\mathrm{a}\mathrm{k}\mathrm{e}}(\overrightarrow{\mathrm{x}};\mathrm{t})={\displaystyle \sum _{\mathrm{P}}\frac{{\overrightarrow{\Omega }}_{\mathrm{P}}(\mathrm{t})\times {\overrightarrow{\mathrm{r}}}_{\mathrm{P}}}{4\pi {\mathrm{r}}^3}}(1-\exp ({\mathrm{r}}^3/{\epsilon }^3)$

where ${\overrightarrow{\mathrm{r}}}_{\mathrm{P}}=\overrightarrow{\mathrm{x}}-{\overrightarrow{\mathrm{Z}}}_{\mathrm{P}}$. ${\overrightarrow{\Omega }}_{\mathrm{P}}$ and ${\overrightarrow{\mathrm{Z}}}_{\mathrm{P}}$ are determined by integrating in time the following evolution equations:

(26)$\frac{\mathrm{d}{\overrightarrow{\mathrm{Z}}}_{\mathrm{p}}(\mathrm{t})}{\mathrm{d}\mathrm{t}}=\overrightarrow{\mathrm{u}}({\overrightarrow{\mathrm{Z}}}_{\mathrm{p}};\mathrm{t}),\hspace{1em}\hspace{0.17em}\frac{\mathrm{d}{\overrightarrow{\Omega }}_{\mathrm{p}}(\mathrm{t})}{\mathrm{d}\mathrm{t}}=\left({\overrightarrow{\Omega }}_{\mathrm{p}}(\mathrm{t})\nabla \right)\overrightarrow{\mathrm{u}}({\overrightarrow{\mathrm{Z}}}_{\mathrm{p}};\mathrm{t})$

Being unconnected, the vorticity field defined by the vortex blobs will be distorted and numerical errors will accumulate. In order to correct this, remeshing is applied, which, besides the improvement in accuracy, also makes it possible to apply the particle-mesh method that substantially reduces the computational cost. Especially when the implementation makes use of the advantages parallel processing offers, it is possible to perform complete 10 min simulations that will involve several millions of particles. In the particle mesh method, free (space) vorticity is first projected on a Cartesian grid and then the Poisson equation ${\nabla }^2\overrightarrow{\psi }=-\overrightarrow{\omega }$ is solved over the same grid by means of a fast Fourier transformation solver [35]. Then ${\overrightarrow{\mathrm{u}}}_{\mathrm{w}\mathrm{a}\mathrm{k}\mathrm{e}}=\nabla \times \overrightarrow{\psi }$ and its derivatives appearing in (26) are calculated through finite differences.

2.3. Hydrodynamic Modeling

2.3.1. Potential Flow Modeling

In the time domain, the equation of motion of a rigid floater takes the following form:

(27)$(\mathsf{M}+{\mathsf{a}}_{\infty }){\ddot{\mathsf{q}}}_{\mathrm{f}}+{\displaystyle \underset{0}{\overset{\mathrm{t}}{\int }}\mathsf{R}(\mathrm{t}-\tau ){\dot{\mathsf{q}}}_{\mathrm{f}}(\tau )\mathrm{d}\tau +({\mathsf{K}}_{\mathrm{H}}+{\mathsf{K}}_{\mathrm{G}}+{\mathsf{K}}_{\mathrm{M}\mathrm{o}\mathrm{o}\mathrm{r}})}{\mathsf{q}}_{\mathrm{f}}={\mathsf{F}}_{\mathrm{e}\mathrm{x}\mathrm{c}}^{\left(1\right)}+{\mathsf{F}}_{\mathrm{e}\mathrm{x}\mathrm{c}}^{\left(2\right)}+{\mathsf{F}}_{\mathrm{M}\mathrm{o}\mathrm{o}\mathrm{r}}+{\mathsf{F}}_{\mathrm{W}\mathrm{T}}+{\mathsf{F}}_{\mathrm{v}\mathrm{i}\mathrm{s}\mathrm{c}}+{\delta }_{\iota 3}(\mathrm{B}-\mathrm{W})$

where vector ${\mathsf{q}}_{\mathrm{f}}$ contains the three displacements and three rotations of the floater; dots denote time derivatives; M is the 6 × 6 generalized mass matrix of the floater; ${\mathrm{a}}_{\infty }$ denotes the added mass matrix corresponding to infinite wave frequency; R(t) denotes the retardation matrix in the time convolution that captures linear hydrodynamic damping due to memory effects and is defined by,

(28)${\mathrm{R}}_{\mathrm{i}\mathrm{j}}(\mathrm{t})=\frac{2}{\pi }{\displaystyle \underset{0}{\overset{\infty }{\int }}{\mathrm{b}}_{\mathrm{i}\mathrm{j}}(\omega )\cos (\omega t)d\omega }$

where b_ij denote the elements of the added damping matrix b; K_H, K_G and K_Moor denote the stiffness matrices due to hydrostatics, gravity and mooring lines (zero if the dynamic mooring line model is used); $F_{\mathrm{e}\mathrm{x}\mathrm{c}}^{\left(1\right)}$ denotes the first --order wave exciting force accounting for the wave force and diffraction, and is defined as follows:

(29)$\begin{array}{l}{F}_{\mathrm{e}\mathrm{x}\mathrm{c}}^{(1)}(\mathrm{t})=\mathrm{A}\left[\frac{F_{\mathrm{e}\mathrm{x}\mathrm{c}}^{(1)}(\omega )}{\mathrm{A}}\right]\cos (\phi (\omega )-\omega \mathrm{t})\\ F_{\mathrm{e}\mathrm{x}\mathrm{c}}^{(1)}(\mathrm{t})={\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{A}}_{\mathrm{i}}\left[\frac{F_{\mathrm{e}\mathrm{x}\mathrm{c}}^{(1)}({\omega }_{\mathrm{i}})}{\mathrm{A}}\right]\cos ({\phi }_{\mathrm{i}}({\omega }_{\mathrm{i}})-{\omega }_{\mathrm{i}}\mathrm{t}+{\epsilon }_{\mathrm{i}})},{\mathrm{A}}_{\mathrm{i}}=\sqrt{2\mathrm{S}({\omega }_{\mathrm{i}})\mathrm{d}\omega }\end{array}$

where the first expression in (29) is valid for a single Airy wave of height 2A, while the second expression is valid for irregular waves described by a spectrum with spectral density S(ω). The term in square brackets and φ(ω) correspond to the magnitude (per wave amplitude) and the phase angle of the exciting force calculated in the frequency domain, while ε_i is the random phase angle of every ith wave component, that is uniformly distributed over [0, 2π]; $F_{\mathrm{e}\mathrm{x}\mathrm{c}}^{\left(2\right)}$ in (27) denotes the second-order wave-induced force vector that is approximated using Newman’s approximation [13]:

(30)$\begin{array}{l}{F}_{\mathrm{e}\mathrm{x}\mathrm{c}}^{\left(2\right)}(\mathrm{t})={\left[{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{A}}_{\mathrm{i}}\sqrt{2\left[\frac{F_{\mathrm{d}\mathrm{r}\mathrm{i}\mathrm{f}\mathrm{t}}^{}({\omega }_{\mathrm{i}})}{{\mathrm{A}}^2}\right]} \mathrm{c}\mathrm{o}\mathrm{s}(-{\omega }_{\mathrm{i}}\mathrm{t}+{\epsilon }_{\mathrm{i}})}\right]}_{F_{\mathrm{d}\mathrm{r}\mathrm{i}\mathrm{f}\mathrm{t}}^{}({\omega }_{\mathrm{i}})>0}^2\\ -{\left[{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{A}}_{\mathrm{i}}\sqrt{-2\left[\frac{F_{\mathrm{d}\mathrm{r}\mathrm{i}\mathrm{f}\mathrm{t}}^{}\left({\omega }_{\mathrm{i}}\right)}{{\mathrm{A}}^2}\right]} \cos (-{\omega }_{\mathrm{i}}\mathrm{t}+{\epsilon }_{\mathrm{i}})}\right]}_{F_{\mathrm{d}\mathrm{r}\mathrm{i}\mathrm{f}\mathrm{t}}^{}({\omega }_{\mathrm{i}})<0}^2 , {\mathrm{A}}_{\mathrm{i}}=\sqrt{2\mathrm{S}\left({\omega }_{\mathrm{i}}\right)\mathrm{d}\omega }\end{array}$

where F_drift is the mean drift force.

In the above, the added mass and damping matrices, the exciting force and the mean drift force vectors are calculated a priori by solving the first-order hydrodynamic problem in the frequency domain using panel solvers such as WAMIT [48] or freFLOW [49,50].

2.3.2. Morison’s Equation

According to Morison’s equation [15], the hydrodynamic force per unit length L of an inclined, submerged and moving body takes the following form:

(31)$\mathrm{d}\mathrm{F}/\mathrm{d}\mathrm{L}= \underset{\mathrm{F}\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{d}\mathrm{e} \mathrm{k}\mathrm{r}\mathrm{y}\mathrm{l}\mathrm{o}\mathrm{v}}{\underbrace{\rho \mathrm{d}\mathrm{V} {\mathrm{a}}_{\mathrm{n}}}}+\underset{\mathrm{D}\mathrm{i}\mathrm{f}\mathrm{f}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}{\underbrace{{\mathrm{C}}_{\mathrm{a}} \rho \mathrm{d}\mathrm{V} {\mathrm{a}}_{\mathrm{n}}}}-\underset{\mathrm{A}\mathrm{d}\mathrm{d}\mathrm{e}\mathrm{d} \mathrm{m}\mathrm{a}\mathrm{s}\mathrm{s}}{\underbrace{{\mathrm{C}}_{\mathrm{a}} \rho \mathrm{d}\mathrm{V} {\ddot{\mathrm{q}}}_{\mathrm{n}}}}+\underset{\mathrm{D}\mathrm{r}\mathrm{a}\mathrm{g}\mathrm{ }\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{m}}{\underbrace{0.5 {\mathrm{C}}_{\mathrm{d}} \rho \mathrm{d}\mathrm{S} |\mathrm{u}-\dot{\mathrm{q}}{|}_{\mathrm{n}}{(\mathrm{u}-\dot{\mathrm{q}})}_{\mathrm{n}}}}$

where a and u denote the total (wave and current) induced acceleration and velocity, respectively, which are derived from Airy’s or stream function wave theory, while $\ddot{\mathrm{q}}, \dot{\mathrm{q}}$ denote the acceleration and the velocity of the body. Subscript n denotes the normal to the surface component. C_a and C_d are the added mass and drag coefficients, respectively. dV and dS denote the infinitesimal volume and the surface normal to the flow, respectively (i.e., for a cylinder of radius R, dV = πRP²PdL and dS = 2RdL). ρ denotes the density of the water. In the above expression, the total force is the sum of the Froude Krylov force, the diffraction force, the added mass force and the drag force.

The Froude Krylov force, as considered in Morison’s equation, is valid only for fully submerged bodies, as reported, for example, in [11]. In the surface piercing (floating) case, the Froude Krylov force on the end nodes should be expressed in terms of the dynamic pressure p_d over the wetted surface.

2.3.3. Hydrostatics

In case the floater is rigid, the total buoyancy (concentrated load) and the linearized hydrostatic stiffness terms are directly considered in the equation of motion (see Equation (27)). On the contrary, if the substructure is flexible, the hydrostatic pressure is integrated along the flexible members. When conical members are used, pressure integration over their wet surface leads to the following expressions (defined in local coordinates):

(32)$$\begin{gathered} \mathrm{d}{F}_{\mathrm{b}}/\mathrm{d}\mathrm{L}=\left\{\begin{array}{c}\rho \mathrm{g}\mathrm{S}{\mathrm{T}}_{31}\cos {\mathsf{\theta }}_{\mathrm{c}}\\ -\mathrm{P}(\rho \mathrm{g}{\mathrm{z}}_{\mathrm{g}}+{\mathrm{p}}_{\mathrm{d}})\sin {\mathsf{\theta }}_{\mathrm{c}}\\ \rho \mathrm{g}\mathrm{S}{T}_{33}\cos {\mathsf{\theta }}_{\mathrm{c}}\end{array}\right\}\mathrm{d}\mathrm{L}, \mathrm{d}{M}_{\mathrm{b}}/\mathrm{d}\mathrm{L}=\left\{\begin{array}{c}\rho \mathrm{g}\mathrm{S}\mathrm{R}{\mathrm{T}}_{33}\sin {\mathsf{\theta }}_{\mathrm{c}}\\ 0\\ -\rho \mathrm{g}\mathrm{S}\mathrm{R}{\mathrm{T}}_{33}\sin {\mathsf{\theta }}_{\mathrm{c}}\end{array}\right\}\mathrm{d}\mathrm{L} \\ {\mathrm{F}}_{\mathrm{b}\mathrm{n}}=\rho \mathrm{g}{\mathrm{z}}_{\mathrm{g}}\mathrm{S} \end{gathered}$$

where the first two terms (distributed force and moment per unit length) are applied on the side surfaces, while the last term (concentrated force) is applied at the end nodes (e.g., bottom of the floater). In the above equations, S, P and R denote the area, the circumference and the radius of the member at the point where the buoyancy force is evaluated, respectively, z_g the global vertical coordinate of the center of the section, p_d the dynamic pressure, θ_c the cone angle of the member and T_ij the elements of the local to global transformation matrix. In the present convection, the y-axis in the local system coincides with the beam axis.

It is noted that, by linearizing the rotation matrix T and the global position vector z_g in (32), convergence is significantly improved, whereas hydrostatic stiffness terms are inherently introduced in the equation of motion of the floater.

2.4. Wind and Wave Excitation

Wind inflow is specified by defining its three velocity components (u, v, w) over a space grid (wind box). The wind field is decomposed into its deterministic and stochastic parts. The deterministic part may vary in space (e.g., due to shear and topography) and time (e.g., in the case of a deterministic wind gust). For the stochastic part, time-series of the wind fluctuations ${\mathrm{u}}^{\prime }, {\mathrm{v}}^{\prime }, {\mathrm{w}}^{\prime }$, in all three directions in space are determined. These series are obtained by specifying the power spectral density (PSD) function S(ω) of all wind components and a spatial correlation function. The standard choice is the von Karman spectrum, combined with Davenport’s exponential coherence function [51]. Time series of turbulent fluctuations are generated over the inlet cross section of the wind (in the yz plane), by combining inverse Fast Fourier Transform (FFT) and sampling of distributed random phases (method proposed by Veers’ [52]). By applying Taylor’s hypothesis, which postulates that disturbances travel with the average wind velocity ${\mathrm{U}}_{\mathrm{W}}$, time dependence is transformed into x dependence; $\mathrm{x}=\mathrm{t} \hspace{0.17em}{\mathrm{U}}_{\mathrm{W}}$, where x is in the wind prevailing direction and so turbulent velocity fluctuations are defined over the so called “turbulent wind box”.

Incoming waves are generated in a similar approach. The spectrum concerns the energy of the wave and is associated with the wave height. Pierson–Moskowitz or JONSWAP spectra [11] are usually considered, but also other site-specific ones may be used. The shape of the spectrum is parametrically defined by the peak period T_P and the significant wave height H_S, as well as additional specific shape parameters depending on the spectrum considered (see [11]).

For stochastic wave excitation, wave kinematics is almost exclusively based on linear Airy theory and the wave field is defined as a superposition of regular sinusoidal waves. The surface elevation of an irregular sea state described by a S(ω) is defined as follows:

(33)$\zeta (\mathrm{x};\mathrm{t})={\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{A}}_{\mathrm{i}}\cos ({\mathrm{k}}_{\mathrm{i}}\mathrm{x}-{\omega }_{\mathrm{i}}\mathrm{t}+{\epsilon }_{\mathrm{i}})},{\mathrm{A}}_{\mathrm{i}}=\sqrt{2\mathrm{S}\left({\omega }_{\mathrm{i}}\right)\mathrm{d}\omega }$

where ε_i and k_i denote the random phase angle uniformly distributed over [0, 2π] and the wave number, respectively, of the ith wave component. A is the wave amplitude and x is the space coordinate along which ζ is evaluated. It is noted that, in Airy theory, wave kinematics is, by definition, valid up to the mean sea water level (MSL). When using Morison’s equation, extrapolation methods may apply to extend the application range of Airy theory up to the instantaneous water level (IWL). Amongst the various engineering extrapolation methods, Wheeler stretching [53] is here considered.

2.5. Controller

A built-in baseline generic control module may be considered, which simulates a standard power-speed controller (variable speed/variable pitch), designed to maximize power output in partial load operation and to regulate power/rotor speed in above rated operation, respectively. In the present case, standard, linear Proportional Integral Differential (PID) control elements and linear filters are used. Start-up and shut-down modes are also available in order to fulfill the requirement of the standards for simulating start and stop sequences. Additional control logic is implemented in the baseline controller that accounts for individual pitch and/or flap control for load reduction [54,55,56], thrust peak shaving near rated wind speed for load mitigation (standard option in modern commercial turbines) and soft shut down/storm control option to delay the shut-down of the turbine above the standard cut-out wind speed (i.e., 25 m/s). Alternatively, it is possible to plug-in in the simulation external control libraries, which solve internally the control equations and provide control variables.

The controller interacts with the structural dynamic module at the generator side drive-train node, where input generator torque is applied, and at the root of the blades, where pitch motion is imposed.

3. Numerical Results and Discussion

The numerical results next presented, concern the NREL 5MW reference wind turbine (WT) [57] of the Offshore Code Comparison Collaboration Continuation (OC4) phase II and the Offshore Code Comparison Collaboration (OC3) phase IV IEA Annexes, in which it was mounted on a semi-submersible [58] and a spar-buoy floater [59], respectively. The baseline NREL controller [57] is employed, which implements a variable speed/variable pitch power controller.

In Section 3.1, emphasis is given on the modeling of the floating support structure and the mooring lines, as well as on the verification of the coupled methodology. Predictions are compared against data of state-of-the-art numerical tools that participated in OC4 project in 2015 [37]. Both available hydrodynamic modeling options (i.e., potential theory and Morison’s equation) are employed, alongside with the flexibility of the floater, although the floater is rather stiff. In Section 3.2, a challenging case regarding aerodynamics is presented, in which the operation of a floating wind turbine in half-wake condition is considered. The free-wake vortex particle method, medium fidelity aerodynamic option is employed in order to study one- and two-way wake interactions, as well as to calibrate BEMT.

3.1. Verification of the Hydro-Servo-Aero-Elastic Modeling

In this case the WT was mounted on a steel semi-submersible floater at a depth of 200 m. Four codes were selected for comparisons, out of the list of the participating codes in OC4 phase II, as shown in Table 1. All codes use BEMT aerodynamic modeling. The present (hGAST), HAWC2 and Orcaflex are FEM-based codes and also adopt dynamic mooring line modeling, while BLADED and FAST are modal-based codes with quasi static mooring line modeling. HAWC2 and BLADED use Morison’s equation for the hydrodynamic modeling, while FAST and Orcaflex use potential theory.

As mentioned, both hydrodynamic models are implemented in present code, denoted hereafter as hGASTp and hGASTm, to distinguish potential and Morison-based application, respectively. Wave kinematics in Morison’s equation is calculated in the instantaneous position (IP) of the body and up to the instantaneous water level (IWL) using Wheeler stretching extrapolation method. Moreover, the Froude Krylov force in Morison’s equation at the end nodes is expressed in terms of the dynamic pressure, while hydrostatics (i.e., buoyancy force and restoring moments) is calculated in nonlinear manner using the surface integration method (see Section 2.3). It is noted that, in the OC4 project, the considered hydrodynamic coefficients in Morison’s equation were calibrated on the basis of the linear potential theory to allow consistent comparison between Morison and potential based code predictions, while additional damping was introduced in the potential based tools (i.e., drag term in Morison’s equation or external damping terms) to account for viscous drag effects [58].

At first, an overall verification of the dynamic modeling is performed by comparing the natural frequencies of the coupled semi-submersible FOWT system in Table 2. Gravity and structural damping contributions are considered. In general, the agreement is good, indicating consistent modeling of the WT and of the semi-submersible floater. It also proves that the dynamics of the floating system are well captured. Rigid body modes of the floater are similar. The Morison’s equation version of the present code and HAWC2 predict slightly reduced surge and sway frequencies of 0.0088 Hz, compared to 0.0093 Hz predicted by the other codes, including the potential version of the present code. Both versions of the present code predict slightly bigger frequencies for the roll and pitch rigid body motions of the floater in comparison to the others.

As regards the WT frequencies, Orcaflex predicts higher tower bending and drive train torsional frequencies. First P blade modes are consistent in the majority of the codes, except that HAWC2 predicts lower blade asymmetric pitch and yaw frequencies in the edgewise direction. Differences between modal and FEM codes are found in the second tower fore-aft and side-to-side bending modes and the second blade asymmetric yaw mode in the flapwise direction. Modal codes overestimate the frequency of the aforementioned modes. For example, FAST over predicts the second tower side-to-side bending mode at 5 Hz, while the other codes predict it at ~3.5 Hz. In addition, both modal codes overestimate the second blade asymmetric yaw mode in the flapwise direction, giving ~1.9 Hz, instead of 1.68 Hz predicted by hGAST and Orcaflex and 1.61 Hz predicted by HAWC2.

The effect of the different hydrodynamic modeling on the natural frequencies is identified by comparing the two different approaches of the present model. WT frequencies remain unchained, while small differences are found in the six rigid body modes of the floater. In the present case, application of Morison’s equation gives slightly lower frequencies for the six floater motions, which are in agreement with the HAWC2 results. However, since the Morison’s equation is subjected to calibration, a firm conclusion cannot be drawn.

Next, time domain simulation results are presented with all flexibilities of the structure enabled. Uniform wind speed at 8 m/s and an Airy wave of 6 m height and 10 s period, aligned with wind and x-axis, is considered at 200 m depth. The controller is operating in variable speed mode.

In Figure 9, the rigid body motions of the floater are presented. In general, the agreement is good amongst the codes compared. The wave period of 10 s is well represented in all signals. At the 0° wave angle, surge, heave and pitch motions are directly excited by the wave, while sway, roll and yaw motions are only excited due to couplings and; therefore, attain much smaller amplitudes. FAST and the potential version of the present code, that both apply first P order potential theory, predict almost identical surge motion. Small differences in the surge mean value could be attributed to nonlinear hydrodynamic effects [37] or to different aerodynamic thrust force. In the Morison’s equation approach, drift effects are caused when the hydrodynamic loads are applied at the instantaneous position (IP) of the body and/or at the instantaneous water level (IWL) using Wheeler’s stretching. This is clearly depicted by comparing hGAST potential (solid blue line) and Morison-based predictions (dashed blue line). Morison’s equation predicts a negative drift effect in the surge motion. The mean surge is 0.3 m lower as compared to the results based on potential theory that does not include drift effects.

The pitch signals are very similar and the mean value, which is determined by the aerodynamic thrust, is consistent. The higher amplitude motion (especially in heave and pitch motions, but also on loading signals in fore-aft direction) predicted by some Morison-based codes could be probably reduced if the dynamic pressure modification discussed in Section 2.3.2 was applied. A small phase, lagging between the potential and Morison-based predictions by the present code, appears because Morison’s equation is applied at the instantaneous position, while potential theory at the undisturbed, reference position of the floater. A minor increase in the pitch mean value by ~0.2° is also predicted by Morison’s equation.

Although sway, roll and yaw motions attain small values, the agreement is very good. These motions are less damped, so, in certain cases, the provided results have not completely reached a periodic state. The positive roll angle is set by the torque of the rotor. Differences in the sway sign are not straightforward to explain, because of the very small values attained. Heave motion is well represented. Again some of the codes using Morison’s equation predict higher amplitudes.

In Figure 10a, the tension of the mooring line at fairlead 2 (upwind position) is compared. The mean values are identical, while phase shift and higher frequencies appear in the codes that adopt the dynamic mooring lines modeling, as compared to those that adopt the quasi static approach. These differences do not influence the dynamics of the coupled FOWT system (see Figure 9), but will affect the Damage Equivalent Loads (DELs) of the mooring lines [37], which could be of crucial importance in the design of the mooring lines and the foundation.

Figure 10b presents the local axial force of the upper horizontal bracer of the floater at the connection point between the first vertical column and the upper delta connection. Only Morison-based codes that consider the flexibility of the floater (through 1D beam modeling) provided results for this signal. Good agreement is noted, whereas the wave period of 10 s is again depicted in the loading signal.

In Figure 11, the fore-aft and yawing moments at the tower base are compared. The wave period is depicted in both tower signals. Signals in the fore-aft direction are mainly affected by the wave excitation. Higher amplitude predictions in Morison-based codes are explained from the higher amplitudes of the pitch motion. Expected differences in the mean value of the tower yaw moment between FEM and modal-based codes are masked by the increase in amplitude due to the rigid motions of the floater. Additionally, high frequency excitations are predicted by most of the codes, as a result of the lowly damped yaw motion as compared to the fore-aft motion.

In Figure 12, the out-of-plane moment at the root of the blade and the torsion angle at its tip are presented. The rotor period of about 6 s is dominant in both signals. Bending moment in the out-of-plane direction is consistent, although slight differences exist. The out-of-plane signals are also affected by the presence of the wave, leading to increased amplitudes, as compared to onshore or bottom-fixed supporting structures. Both wave and rotor frequencies are depicted.

Differences are seen in the torsion angle. The present model and HAWC2, that both account for geometric nonlinear effects due to large deflections and rotations, predict similar results with the same phase, but with different mean value and amplitude [18]. In BLADED results, the difference in phase and the smaller amplitude is possibly explained by the use of a simplified beam model.

Out of the above comparisons, it is concluded that the present model (hGAST) is thoroughly verified through code-to-code comparisons against state-of-the-art tools participated in IEA OC4 Annex. Moreover, both potential and Morison-based hydrodynamic modules may provide similar predictions in the case that Morison’s equation is properly tuned. Morison’s equation models can be adopted in cases where the flexibility of the floater is considered. Finally, phase shift and higher frequencies appeared in the mooring line tension in the codes that adopt the dynamic mooring lines modeling, as compared to those that adopt the quasi static approach. These differences are not found to influence the dynamics of the coupled FOWT system but could be critical in the design of the mooring lines and the foundation.

3.2. Half-Wake Effects on Floating Offshore Wind Turbines (FOWT)

The present test case concerns two floating wind turbines placed at a distance of 5 diameters (D) along the prevailing wind direction with an off-set of 1 radius (0.5 D) in the lateral direction, so half of the rotor of the second FOWT is submerged into the wake of the first. A simplification that BEMT aerodynamic cannot avoid, assumes that there is one-way interaction, meaning that only the wake of the preceding WT affects the performance of the second one. In theory the interaction is two ways, in the sense that the incoming wake is modified by the presence of the second WT. A similar two-way interaction also takes place when inflow turbulence impinges a WT. In both cases, the incoming disturbance can be interpreted as a connection of free vortices, which renders vortex methods a natural choice.

An illustration of how vortex methods can be applied in this half-wake interaction problem is shown in Figure 13. The two WTs are represented by the grid defined on the blades (thin lifting surfaces), while the wakes are depicted by an iso-vorticity surface in scatter data format. In the upwind wind turbine, the formation of the tip and root vortices is clearly represented over the early stages of the wake development, while in the downwind wind turbine the mixing of the two wakes is substantial and does not allow distinguishing well shaped vortex structures. In the figure the front view of the wake development is also shown, indicating the way the preceding wake covers the downwind rotor, as well as the expansion of the wake.

Quantitative information is given by the instantaneous and averaged normalized deficit velocity profiles at 5 D and 8 D in the wake of an isolated WT, that are shown in Figure 14. The velocity deficit attains a mean value of about 35% at 5 D and 30% at 8 D, respectively. The dots in the figure represent the raw instantaneous data over the last period obtained from the simulation every time step, while the lines give the corresponding average, which is used as input in the one-way interaction model simulations.

The simulations performed are listed in Table 3 and concern the NREL 5MW Reference WT [57] mounted on the OC3 spar-buoy floater [59]. In the numerical simulations, only the second WT is flexible. Uniform wind inflow (no ambient turbulence, yaw and shear) is chosen at the rated wind speed of 11.4 m/s, while, following the definition of the OC3 project, a regular wave of H = 6 m and T = 10 s is considered at the water depth of 320 m. Rated conditions have been chosen, because for a variable pitch, variable speed WT, maximum loading is obtained leading to highest momentum deficit in the wake. The list of the load cases concerns the isolated case in uniform inflow conditions (1 rotor), the one-way interaction modeling or frozen wake at x = 5 D and the two-way interaction modeling (2 rotors) performed either using the BEMT or the vortex free-wake aerodynamic method. Next, these modeling options are compared in terms of time series and power spectral densities (PSDs) for selected signals recorded on the second (hydro-aero-elastic) WT.

In Figure 15a, time series of 200 s long are presented for the rotor speed. The excitation by the wave of 10 s period is clearly seen. In the two-way interaction simulation, at about t = 55 s the wake of the upwind rotor has covered the 5 D distance. The interaction leads to a slight reduction of the mean value of the rotor speed which becomes almost identical to the vortex based one-way interaction simulation (black solid and dashed light blue lines). In both BEMT sets of results, the mean rotor speed obtained is 3% lower than that vortex models predict.

In Figure 15b, time series of 200 s long are presented for the yaw angle of the floater. As previously, transition from uniform to half-wake inflow conditions starts at about t = 55 s in the two-way results. Before, the yaw motion is the same as in the isolated case, whereas, when the transient dies out at about t = 100 s, the motion converges to that of the one-way interaction model with minor deviation in amplitude. In general, BEMT and vortex aerodynamic predictions agree well. It is noted that the yaw motion is chosen here because the wave dominates the rest of the motions, so a clear effect of the half-wake interaction is only visible in the yaw motion due to the loading asymmetry (wave excitation in yaw motion of the spar-buoy floater is zero).

Next, focusing on WT loading, time series of the fore-aft bending moment at the tower base are presented in Figure 16a. Since the half-wake effect is masked by the wave loading, the mean value and the amplitude remain almost unchanged, although high excitation frequencies appear. In order to identify them, in Figure 16b the PSD of the tower base fore-aft bending moment is presented. Clearly the wave frequency at 0.1 Hz attains the highest value. However, in both the one- and two-way interaction models, the 3P (three per revolution) frequency at 0.6 Hz is also excited, as well as the tower natural frequency at 0.45 Hz. The latter is found to be more excited in the two-way model, as shown in both plots.

In Figure 17, the tower yaw moment signal is considered. Time averaged values and amplitudes shown in Figure 17a increase when there is half-wake effect, which is responsible for the loading asymmetry on the rotor. In quantitative terms, there is fair agreement. In the interaction case, the BEMT amplitudes are lower (almost half). In Figure 17b, the PSD of the yaw moment is presented. High excitation takes place at 0.1 Hz, which is the wave frequency, but also at 0.6 and 1.2 Hz, which are the 3P and 6P frequencies of the rotor, respectively. In terms of peaks, the agreement is good. However, the excitation in the two-way interaction spectrum is broader, especially at low frequencies. Because the wave excites the yaw direction only indirectly through the couplings with the other floater motions, the most excited frequency is the 3P at 0.6 Hz, also visible in the time series.

Finally, in Figure 18 the signal of the out-of-plane bending moment at the blade root is considered. In the time series, the mean level and the min-to-max amplitude are in fair agreement. However, in their details the variations of the specific load are different. Deviations of this kind are reflected on the fatigue loading. Additionally, phase differences are noted. Because it is not possible to synchronize the rotor speed with the wave elevation, such differences are expected. The PSDs shown on the right indicate the dominance of the wave seen at 0.1 Hz, but also the significance of the blade passage through the incoming wake flow at 0.2 Hz which corresponds to 1P. It is worth noticing that 2P and 3P are also excited in both modeling options of the half-wake interaction case. The agreement between the two interaction models is good.

Out of the above comparisons, it follows that one- and two-way interaction modeling lead to similar results. It is also found that BEMT, which can only model the one-way interaction, is consistent as compared to vortex modeling. This point requires a deeper investigation that should cover a variety of designs and load cases over the full operation span. Deviations appeared in the load ranges which could affect fatigue. However, in order to support this point, proper fatigue simulations are needed.

4. Concluding Remarks

An all-inclusive hydro-servo-aero-elastic simulation tool for FOWTs was presented together with indicative results. By following the generalized formulation of analytic dynamics, a simulation environment was defined, in which interfaces realizing kinematic and loading compatibility allow to plug-in different sub-models without canceling any of the couplings or changing the non-linear character of the underlying interactions. In this respect, sub-models can be any of the aerodynamic, hydrodynamic, elasto-dynamic and controlling aspects of the full problem. One advantage of such an approach is that formal linearization reformulates the equations in incremental form, which facilitates the iterative solution of the fully non-linear problem and at the same time formulates the linear stability problem for the full configuration [18].

Targeting holistic analysis of complex systems involves compromises that are linked to the available computing capacity, the time frame and the amount of simulations needed. In the present simulation tool, the main compromise concerns the aerodynamic modeling and the choice amongst the BEMT and the vortex options. The BEMT version was verified via the code-to-code comparisons with other state-of-art BEMT tools for the semi-submersible 5 MW FOWT in Annex OC4 (results in Section 3.1). From the comparisons, it was concluded that both potential and Morison-based hydrodynamic modules may provide similar predictions in the case that Morison’s equation is properly calibrated. Additionally, phase shift and higher frequencies appeared in the mooring line tension in the codes that adopt the dynamic mooring lines modeling, as compared to those that adopt the quasi static approach. These differences are not found to influence the dynamics of the coupled FOWT system but could be critical in the design of the mooring lines and the foundation. This justifies the choice of engaging higher order nonlinear submodels in the analysis of the complex FOWT system.

By implementing cost reduction techniques, vortex-based simulations of the certification type were made feasible. This allowed a further verification of BEMT predictions in [2], which concluded that, in terms of equivalent loads, BEMT is conservative in comparison to vortex-based modeling. Even in the complex case of half-wake interaction problem, discussed in Section 3.2, the low-cost BEMT performs well, while the one- and two-way interaction modeling in the free-wake aerodynamic method were found to lead to similar results.

Herein, as well as in previous works, the verification campaigns were limited to specific designs and certain loading conditions, in that further confirmation of the findings obtained so far is needed by comparing vortex and CFD-based predictions with BEMT. Important steps in this direction are to actively include inflow turbulence, as indicated in [60], and to extend the comparisons in other operational conditions that are currently poorly addressed, as, for example, in the occurrence of stall- and vortex-induced vibrations.

Author Contributions

Conceptualization, D.I.M., V.A.R. and S.G.V.; methodology, all; software, D.I.M., G.P.P. and V.A.R.; validation, D.I.M. and S.G.V.; formal analysis, all; writing-original draft, D.I.M. and V.A.R.; writing-review and editing, V.A.R. and S.G.V.; supervision, S.G.V. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the European projects INNWIND.EU (grant number: FP7-ENERGY-2012-1-2STAGE-308974), AVATAR (grant number: FP7-ENERGY.2013.2.3.1/no. 608396), JABACO (grant number: RFSR-CT-2015-00024) and REFOS (grant number: RFCS-02-2015-709526), and the Greek national projects SYNERGASIA-AVRA and ARISTEIA-POSEIDON. The authors would like to acknowledge the European Union and the Greek secretary for the funding.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Sub-components of the multi-physics floating offshore wind turbine (FOWT) system.

Figure 2. (a) Local co-ordinate systems of the different flexible bodies forming the wind turbine system, (b) application of multibody kinematics at the turbine level (connected bodies).

Figure 3. Application of multibody kinematics at the body level (connected sub-bodies).

Figure 4. Definition of body local deflections.

Figure 5. Truss element in undeformed and deformed states.

Figure 6. Kinematic and dynamic coupling of the turbine and the floater with the mooring line.

Figure 7. Aerodynamic modeling options of lifting bodies.

Figure 8. Wake representation by free vortex particles.

Figure 9. Comparison of platform motions: (a) Platform surge motion, (b) platform sway motion, (c) platform heave motion, (d) platform pitch motion, (e) platform roll motion and (f) platform yaw motion (Airy wave H = 6 m, T = 10 s, uniform inflow at 8 m/s).

Figure 10. Comparison of (a) mooring line tension at fairleads 2 (upwind) position and (b) local axial force of the upper horizontal member of the floater at the connection point between column 1 and delta connection (Airy Wave H = 6 m, T = 10 s, uniform inflow at 8 m/s).

Figure 11. Comparison of the tower base (a) fore-aft moment and (b) yawing moment (Airy wave H = 6 m, T = 10 s, uniform inflow at 8 m/s).

Figure 12. Comparison of the (a) blade root out-of-plane bending moment and (b) blade tip torsion angle (Airy wave H = 6 m, T = 10 s, uniform inflow at 8 m/s).

Figure 13. Snapshot of the wake formation in the case of two interacting wind turbines (WTs). An iso-vorticity surface is used to mark the formed wake structure.

Figure 14. Instantaneous (dots) and time averaged (dash lines) axial velocity deficits at x = 5 D and 8 D. The green and red dots correspond to downwind positions of 5 D and 8 D respectively.

Figure 15. Time-series of the (a) rotational speed and (b) yaw motion of the floater. Comparisons between the isolated inflow case and the one- and two-way interaction models for half-wake modeling at x = 5 D using blade element momentum theory (BEMT) and Vortex (VRX) aerodynamic models for the OC3 spar-buoy FOWT at 11.4 m/s.

Figure 16. (a) Time-series and (b) Power spectral densities (PSD) of the fore-aft bending moment at the tower base. Comparisons between the isolated inflow case and the one- and two-way interaction models for half-wake modeling at x = 5 D using BEMT and Vortex (VRX) aerodynamic models for the OC3 spar-buoy FOWT at 11.4 m/s.

Figure 17. (a) Time series and (b) PSD of the yawing moment at the tower base. Comparisons between the isolated inflow case and the one- and two-way interaction models for half-wake modeling at x = 5 D using BEMT and Vortex (VRX) aerodynamic models for the OC3 spar-buoy FOWT at 11.4 m/s.

Figure 18. (a) Time series and (b) PSD of the out-of-plane bending moment at the blade root. Comparisons between the isolated inflow case and the one- and two-way interaction models for half-wake modeling at x = 5 D using BEMT and Vortex (VRX) aerodynamic models for the OC3 spar-buoy FOWT at 11.4 m/s.

Table 1

Description of the codes selected for comparison out of the list of the OC4 phase II.

Participant	4Subsea, Hvalstad, Norway	Technical University of Denmark (DTU), Roskilde, Denmark	DNV GL, Høvik, Norway	National Renewable Energy Laboratory (NREL), Golden, Colorado, USA	National Technical University of Athens (NTUA), Athens, Greece
Code Name	Orcaflex	HAWC2	BLADED	FAST	hGAST
Structural Dynamics	FEM	FEM	Modal	Modal	FEM
Aerodynamics	BEMT	BEMT	BEMT	BEMT	BEMT
Hydrodynamics	Potential	Morison	Morison	Potential	PotentialMorison
Mooring Lines	Dynamic	Dynamic	Quasi static	Quasi static	Dynamic
Control	NREL 5MW baseline controller

Table 2

Natural frequencies [Hz] comparison of the semi-submersible coupled FOWT of the OC4 phase II.

Mode Description	hGASTp	hGASTm	HAWC2	BLADED	FAST	Orcaflex
Platform surge	0.0093	0.0088	0.0086	0.0094	0.0093	0.0093
Platform sway	0.0093	0.0088	0.0088	0.0092	0.0093	0.0093
Platform heave	0.0583	0.0574	0.0573	0.0581	0.0581	0.0556
Platform roll	0.0413	0.0404	0.0384	0.0397	0.0392	0.0385
Platform pitch	0.0413	0.0404	0.0384	0.0397	0.0392	0.0385
Platform yaw	0.0131	0.0126	0.0132	0.0136	0.0132	0.0127
First tower fore-aft	0.424	0.423	0.424	0.425	0.426	0.465
First tower side-to-side	0.415	0.414	0.415	0.417	0.418	0.458
First drivetrain torsion	0.622	0.622	0.608	0.623	0.628	0.672
First blade collective flap	0.717	0.717	0.686	0.706	0.704	0.692
First asymmetric flapwise pitch	0.677	0.677	0.618	0.641	0.670	0.664
First asymmetric flapwise yaw	0.639	0.639	0.648	0.670	0.667	0.635
First asymmetric edgewise pitch	1.079	1.079	1.002	1.080	1.079	1.095
First asymmetric edgewise yaw	1.092	1.092	1.015	1.091	1.092	1.103
Second tower fore-aft	3.417	3.415	3.314	3.864	3.898	3.405
Second tower side-to-side	3.540	3.537	3.494	3.437	5.012	3.875
Second collective flap	2.000	2.000	1.840	1.972	2.023	1.928
Second asymmetric flapwise pitch	1.876	1.876	1.739	1.718	1.914	1.829
Second asymmetric flapwise yaw	1.681	1.681	1.609	1.870	1.934	1.672

Table 3

Load cases at rated wind speed of 11.4 m/s.

Description of Load Case	Vortex	BEMT
Isolated WT: uniform inflow, no wave excitation, pure aerodynamic simulation	Yes	No
Isolated WT: uniform inflow, wave excitation and flexibility included	Yes	Yes
WT in half-wake operation: one-way interaction wake model (the impacted WT is only simulated)	Yes	Yes
WT in half-wake operation: two-way interaction wake model (both WTs are included in the simulation)	Yes	No

References

1. Jonkman, J.M.; Buhl, M.L., Jr. FAST User’s Guide. Technical Report No. NREL/EL-500-38230; NREL: Golden, CO, USA, 2005; [DOI: https://dx.doi.org/10.2172/15020796]

2. Manolas, D.; Riziotis, V.; Voutsinas, S. Assessment of 3D aerodynamic effects on the behaviour of floating wind turbines. J. Phys. Conf. Ser.; 2014; 555, 012067. [DOI: https://dx.doi.org/10.1088/1742-6596/555/1/012067]

3. Liu, Y.; Qing Xiao, Q.; Incecik, A.; Peyrard, C.; Wan, D. Establishing a fully coupled CFD analysis tool for floating offshore wind turbines. Renew. Energy; 2017; 112, pp. 280-301. [DOI: https://dx.doi.org/10.1016/j.renene.2017.04.052]

4. Leble, V.; Barakos, G.N. A coupled floating offshore wind turbine analysis with high-fidelity methods. Proceedings of the 13th Deep Sea Offshore Wind R&D Conference, EERA DeepWind; Trondheim, Norway, 20–22 January 2016.

5. Madsen, H.A.; Riziotis, V.; Zahle, F.; Hansen, M.O.L.; Snel, H.; Grasso, F.; Larsen, T.J.; Politis, E.; Rasmussen, F. BEM modelling of inflow with shear in comparison with advanced model results. Wind Energy; 2012; 15, pp. 63-81. [DOI: https://dx.doi.org/10.1002/we.493]

6. Schepers, J.G.; Boorsma, K.; Soerensen, N.; Voutsinas, S.G.; Sieros, G.; Rahimi, H.; Heisselmann, H.; Jost, E.; Lutz, T.; Maeder, T. et al. Final results from the EU project AVATAR: Aerodynamic modelling of 10 MW wind turbines. J. Phys. Conf. Ser.; 2018; 1037, 022013. [DOI: https://dx.doi.org/10.1088/1742-6596/1037/2/022013]

7. Voutsinas, S.G. Vortex Methods in Aeronautics: How to make things work. Int. J. Comput. Fluid Dyn.; 2006; 20, pp. 3-18. [DOI: https://dx.doi.org/10.1080/10618560600566059]

8. Menter, F.R. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA J.; 1994; 32, pp. 1598-1605. [DOI: https://dx.doi.org/10.2514/3.12149]

9. Wilcox, D.C. Turbulence Modeling for CFD; DCW Industries: La Canada, CA, USA, 1998; Volume 2.

10. Heinz, J.C.; Sørensen, N.N.; Zahle, F.; Skrypinski, W. Vortex-induced vibrations on a modern wind turbine blade. Wind Energy; 2016; 19, pp. 2041-2051. [DOI: https://dx.doi.org/10.1002/we.1967]

11. Faltinsen, O. Sea Loads on Ships and Offshore Structures; Cambridge University Press: Cambridge, UK, 1993.

12. Duarte, T.; Sarmento, A.; Jonkman, J. Effects of Second-Order Hydrodynamic Forces on Floating Offshore Wind Turbines. Proceedings of the 32nd ASME Wind Energy Symposium; National Harbor, MD, USA, 13–17 January 2014; [DOI: https://dx.doi.org/10.2514/6.2014-0361]

13. Newman, J.N. Second-order, slowly-varying forces on vessels in irregular waves. Proceedings of the International Symposium on Dynamics of Marine Vehicles and Structures in Waves; London, UK, 1–5 April 1974.

14. Mazarakos, T.; Konispoliatis, D.; Manolas, D.; Voutsinas, S.; Mavrakos, S. Modelling of an Offshore Multi–Purpose Floating Structure Supporting a Wind Turbine Including Second–Order Wave Loads. Proceedings of the 11th European Wave and Tidal Energy Conference EWTEC2015; Nantes, France, 6–11 September 2015.

15. Morison, J.; Johnson, J.; Schaaf, S. The force exerted by surface waves on piles. J. Pet. Technol.; 1950; 2, pp. 149-154. [DOI: https://dx.doi.org/10.2118/950149-G]

16. Dean, R. Stream function representation of nonlinear ocean waves. J. Geophys. Res.; 1965; 70, pp. 4561-4572. [DOI: https://dx.doi.org/10.1029/JZ070i018p04561]

17. Hansen, M.O.L.; Sørensen, J.N.; Voutsinas, S.; Sørensen, N.; Madsen, H.A. State of the art in wind turbine aerodynamics and aeroelasticity. Prog. Aerosp. Sci.; 2006; 42, pp. 285-330. [DOI: https://dx.doi.org/10.1016/j.paerosci.2006.10.002]

18. Manolas, D.; Riziotis, V.; Voutsinas, S. Assessing the importance of geometric non-linear effects in the prediction of wind turbine blade loads. J. Comput. Nonlinear Dyn.; 2015; 10, 041008. [DOI: https://dx.doi.org/10.1115/1.4027684]

19. Hodges, D.H. A mixed variational formulation based on exact intrinsic equations for dynamics of moving beams. Int. J. Solids Struct.; 1990; 26, pp. 1253-1273. [DOI: https://dx.doi.org/10.1016/0020-7683(90)90060-9]

20. Bottasso, C.L.; Croce, A.; Savini, B.; Sirchi, W.; Trainelli, L. Aero-servo-elastic modeling and control of wind turbines using finite-element multibody procedures. Multibody Syst. Dyn.; 2006; 16, pp. 291-308. [DOI: https://dx.doi.org/10.1007/s11044-006-9027-1]

21. Kallesøe, B.S. Equations of motion for a rotor blade, including gravity, pitch action and rotor speed variations. Wind Energy; 2007; 10, pp. 209-230. [DOI: https://dx.doi.org/10.1002/we.217]

22. Hodges, D.H.; Dowell, E. Nonlinear Equations of Motion for the Elastic Bending and Torsion of Twisted Nonuniform Rotor Blades; NASA Technical Note D-7818 NASA: Washington, DC, USA, 1 December 1974.

23. Bagherpour, T.; Li, X.M.; Manolas, D.I.; Riziotis, V.A. Modeling of Material Bend-Twist Coupling on Wind Turbine blades. Compos. Struct.; 2018; 193, pp. 237-246. [DOI: https://dx.doi.org/10.1016/j.compstruct.2018.03.071]

24. Bottasso, C.; Campagnolo, F.; Croce, A.; Tibaldi, C. Integrating active and passive load control in wind turbines. Proceedings of the EWEA Annual Event; Copenhagen, Denmark, 16–19 April 2012; pp. 516-521.

25. Øye, S. FLEX4 simulation of wind turbine dynamics. Proceedings of the 28th IEA Meeting of Experts Concerning State of the Art of Aeroelastic Codes for Wind Turbine Calculations; Lyngby, Denmark, 11–12 April 1996.

26. IEC 61400-3-1:2019. Wind Energy Generation Systems—Part 3-1: Design Requirements for Fixed Offshore Wind Turbines; International Electrotechnical Commission (IEC): London, UK, 5 April 2019.

27. Jonkman, J.M. Dynamics Modeling and Loads Analysis of an Offshore Floating Wind Turbine; Technical Report No NREL/TP-500-41958 NREL: Golden, CO, USA, 2007; [DOI: https://dx.doi.org/10.2172/921803]

28. Riziotis, V.A.; Voutsinas, S.G. GAST: A general aerodynamic and structural prediction tool for wind turbines. Proceedings of the EWEC Conference; Dublin Castle, Ireland, 6–9 October 1997; pp. 448-452.

29. Riziotis, V.A.; Voutsinas, S.G.; Manolas, D.I.; Politis, E.S.; Chaviaropoulos, P.K. Aeroelastic analysis of pre-curved rotor blades. Proceedings of the EWEC 2010 Conference; Warsaw, Poland, 20–23 April 2010.

30. Riziotis, V.A.; Manolas, D.I.; Voutsinas, S.G. Free-wake Aeroelastic Modelling of Swept Rotor Blades. Proceedings of the EWEA 2011 Conference; Brussels, Belgium, 14–17 March 2011.

31. Riziotis, V.A.; Voutsinas, S.G.; Politis, E.S.; Chaviaropoulos, P.K. Aeroelastic stability of wind turbines: The problem, the methods and the issues. Wind Energy.; 2004; 7, pp. 373-392. [DOI: https://dx.doi.org/10.1002/we.133]

32. Riziotis, V.A.; Voutsinas, S.G.; Politis, E.S.; Chaviaropoulos, P.K. Assessment of passive instability suppression means on pitch-regulated wind turbines. Wind Energy; 2008; 11, pp. 171-192. [DOI: https://dx.doi.org/10.1002/we.241]

33. Riziotis, V.A.; Voutsinas, S.G.; Politis, E.S.; Chaviaropoulos, P.K. Stability analysis of pitch-regulated, variable-speed wind turbines in closed loop operation using a linear eigenvalue approach. Wind Energy; 2008; 11, pp. 517-535. [DOI: https://dx.doi.org/10.1002/we.276]

34. Wang, K.; Riziotis, V.A.; Voutsinas, S.G. Aeroelastic stability of idling wind turbines. Wind Energy Sci.; 2017; 2, pp. 415-437. [DOI: https://dx.doi.org/10.5194/wes-2-415-2017]

35. Papadakis, G.; Voutsinas, S.G. In view of accelerating CFD simulations through coupling with vortex particle approximations. J. Phys. Conf. Ser.; 2014; 524, 012126. [DOI: https://dx.doi.org/10.1088/1742-6596/524/1/012126]

36. Balls, G.T.; Colella, P. A finite difference domain decomposition method using local corrections for the solution of Poisson’s equation. J. Comput. Phys.; 2002; 180, pp. 25-53. [DOI: https://dx.doi.org/10.1006/jcph.2002.7068]

37. Robertson, A.; Jonkman, J.; Vorpahl, F.; Popko, W.; Qvist, J.; Frøyd, L.; Chen, X.; Azcona, J.; Uzunoglu, E.; Guedes Soares, C. et al. Offshore code comparison collaboration continuation within IEA Wind Task 30: Phase II results regarding a floating semisubmersible wind system. Proceedings of the ASME 2014 33rd International Conference on Ocean, Offshore and Arctic Engineering OMAE; San Francisco, CA, USA, 8–13 June 2014; [DOI: https://dx.doi.org/10.1115/OMAE2014-24040]

38. Stäblein, A.; Hansen, M.H. Timoshenko beam element with anisotropic cross-sectional properties. Proceedings of the VII ECCOMAS; Crete Island, Greece, 5–10 June 2016.

39. Riziotis, V.A.; Katsaounis, G.M.; Papadakis, G.; Voutsinas, S.G.; Bergeles, G.; Tzabiras, G.D. Numerical and experimental analysis of the hydroelastic behavior of purse seine nets. Ocean Eng.; 2013; 58, pp. 88-105. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2012.09.022]

40. Crisfield, M.A. Non-Linear Finite Element Analysis of Solids and Structures; John Wiley and Sons Inc.: New York, NY, USA, 1991; Volume 1.

41. Rahimi, H.; Hartvelt, M.; Peinke, J.; Schepers, J.G. Investigation of the current yaw engineering models for simulation of wind turbines in BEM and comparison with CFD and experiment. J. Phys. Conf. Ser.; 2016; 753, 022016. [DOI: https://dx.doi.org/10.1088/1742-6596/753/2/022016]

42. Schepers, J.G.; Snel, H. Joint Investigation of Dynamic Inflow Effects and Implementation of an Engineering Method for Response Analysis; ECN Report ECN-C-94-107 ECN: Petten, The Netherlands, 1994.

43. Petot, D. Differential Equation Modeling of Dynamic Stall. La Recherche Aerospatiale; 1989; 5, pp. 59-72.

44. Leishman, J.G.; Crouse, G.L. State-Space Model for Unsteady Airfoil Behavior and Dynamic Stall. Proceedings of the 30th Structures, Structural Dynamics and Materials Conference; Mobile, AL, USA, 3–5 April 1989; [DOI: https://dx.doi.org/10.2514/6.1989-1319]

45. Voutsinas, S.G.; Riziotis, V.A.; Chaviaropoulos, P. Non-Linear Aerodynamics and Fatigue Loading on Wind Turbines Operating at Extreme Sites. Proceedings of the 35th AIAA Annual Meeting; Reno, NV, USA, 6–9 January 1997; AIAA Paper 97-0935 [DOI: https://dx.doi.org/10.2514/6.1997-935]

46. Richardson, S.M.; Cornish, A.R.H. Solution of three dimensional incompressible flow problems. J. Fluid Mech.; 1977; 82, pp. 309-319. [DOI: https://dx.doi.org/10.1017/S0022112077000688]

47. Beale, J.T.; Majda, A. Higher order accurate vortex methods with explicit velocity kernels. J. Comput. Phys.; 1985; 58, pp. 188-208. [DOI: https://dx.doi.org/10.1016/0021-9991(85)90176-7]

48. Lee, C. WAMIT—A Radiation–Diffraction Panel Program for Wave-Body Interactions, User Manual; Department of Ocean Engineering: Massachusetts Institute of Technology: Cambridge, MA, USA, 1998.

49. Manolas, D. Hydro-Aero-Elastic Analysis of Offshore Wind Turbines. PhD Thesis; NTUA: Athens, Greece, 2015.

50. Manolas, D.I.; Karvelas, C.G.; Kapogiannis, I.A.; Riziotis, V.A.; Spiliopoulos, K.V.; Voutsinas, S.G. A comprehensive method for the structural design and verification of the INNWIND 10MW tri-spar floater. J. Phys. Conf. Ser.; 2018; 1104, 012025. [DOI: https://dx.doi.org/10.1088/1742-6596/1104/1/012025]

51. Breard, C.; Riziotis, V.; Voutsinas, S. Fatigue loading on wind turbines operating at extreme sites. Proceedings of the EWEC Conference; Dublin Castle, Ireland, 6–9 October 1997; pp. 642-646.

52. Veers, P.S. Three-Dimensional Wind Simulation; Sandia Report SAND 88-0152 Sandia Laboratories: Albuquerque, NM, USA, 1988.

53. Wheeler, J. Methods for calculating forces produced on piles in irregular waves. J. Pet. Technol.; 1970; 22, pp. 359-367. [DOI: https://dx.doi.org/10.2118/2712-PA]

54. Karakalas, A.A.; Manolas, D.I.; Machairas, T.T.; Riziotis, V.A.; Saravanos, D.A. Active Load Alleviation Potential of Adaptive Wind Turbine Blades using Shape Memory Alloy Actuators. Wind Energy; 2019; 22, pp. 620-637. [DOI: https://dx.doi.org/10.1002/we.2311]

55. Manolas, D.I.; Serafeim, G.P.; Chaviaropoulos, P.K.; Riziotis, V.A.; Voutsinas, S.G. Assessment of load reduction capabilities using passive and active control methods on a 10MW-scale wind turbine. J. Phys. Conf. Ser.; 2018; 1037, 032042. [DOI: https://dx.doi.org/10.1088/1742-6596/1037/3/032042]

56. Manolas, D.I.; Spyropoulos, N.A.; Serafeim, G.P.; Riziotis, V.A.; Chaviaropoulos, P.K.; Voutsinas, S.G. Inflow-based flap control on a 10MW-scale wind turbine using a spinner anemometer. J. Phys. Conf. Ser.; 2018; 1037, 032045. [DOI: https://dx.doi.org/10.1088/1742-6596/1037/3/032045]

57. Jonkman, J.M.; Butterfield, S.; Musial, W.; Scott, G. Definition of a 5-MW Reference Wind Turbine for Offshore System Development; Technical Report No NREL/TP-500-38060 NREL: Golden, CO, USA, 2009; [DOI: https://dx.doi.org/10.2172/947422]

58. Robertson, A.; Jonkman, J.; Masciola, M.; Song, H.; Goupee, A.; Coulling, A.; Luan, C. Definition of the Semisubmersible Floating System for Phase II of OC4; Technical Report No NREL/TP-5000-60601 NREL: Golden, CO, USA, 2014; [DOI: https://dx.doi.org/10.2172/1155123]

59. Jonkman, J.M. Definition of the Floating System for Phase IV of OC3; Technical Report No NREL/TP-500-47535 NREL: Golden, CO, USA, 2010; [DOI: https://dx.doi.org/10.2172/979456]

60. Kim, Y.; Lutz, T.; Jost, E.; Gomez-Iradi, S.; Muñoz, A.; Méndez, B.; Lampropoulos, N.; Stefanatos, N.; Sørensen, N.N.; Madsen, H. et al. AVATAR Deliverable D2. 5: Effects of Inflow Turbulence on Large Wind Turbines. 31 December 2016. Technical Report Available online: http://www.eera-avatar.eu/fileadmin/avatar/user/avatar_D2p5_revised_20161231.pdf (accessed on 2 November 2020).