An efficient frequency-domain model for quick load analysis of floating offshore wind turbines

A model for Quick Load Analysis of Floating wind turbines (QuLAF) is presented and validated here. The model is a linear, frequency-domain, efficient tool with four planar degrees of freedom: floater surge, heave, pitch and first tower modal deflection. The model relies on state-of-the-art tools from which hydrodynamic, aerodynamic and mooring loads are extracted and cascaded into QuLAF. Hydrodynamic and aerodynamic loads are pre-computed in WAMIT and FAST, respectively, while the mooring system is linearized around the equilibrium position for each wind speed using MoorDyn. An approximate approach to viscous hydrodynamic damping is developed, and the aerodynamic damping is extracted from decay tests specific for each degree of freedom. Without any calibration, the model predicts the motions of the system in stochastic wind and waves with good accuracy when compared to FAST. The damage-equivalent bending moment at the tower base is estimated with errors between 0.2 % and 11.3 % for all the load cases considered. The largest errors are associated with the most severe wave climates for wave-only conditions and with turbine operation around rated wind speed for combined wind and waves. The computational speed of the model is between 1300 and 2700 times faster than real time.

1 Introduction: the need for an efficient, frequency-domain tool

Offshore wind energy is a key contributor to a carbon-free energy supply. Most of today's offshore wind farms are bottom-fixed, meaning their feasibility is limited to shallow and intermediate water depths. On the other hand, the wind resource in deep water represents an enormous potential that can be unlocked with the deployment of floating wind farms. An important step in making floating wind turbines economically feasible is the application of larger wind turbines and the ability to design the floater to a minimum cost. The design of a floating substructure for offshore wind deployment depends on many design variables, and each possible combination is a potential design. In the design process, the candidate designs need to be simulated in different environmental conditions in order to assess the magnitude of the motions and loads in the system. These simulations are typically carried out with time-domain numerical tools, which allow a representative modelling of the physical phenomena involved and can simulate at about real-time CPU speed. However, this approach can be computationally expensive, especially if one needs to evaluate different floater designs under several environmental conditions. For an improved design process, faster tools are needed to allow optimization in the initial design stage, where the design space has to be thoroughly explored and a broad overview of the system response is desirable.

A few studies of simplified design models for offshore wind turbine floaters exist in the literature. Lupton (2014) presented a frequency-domain numerical tool for the analysis of the OC3-Hywind spar floating wind turbine (Jonkman, 2010), with eight degrees of freedom (DoFs): one normal mode per blade, two tower fore-aft modes, and floater surge, heave and pitch. The model included linear hydrodynamics computed with a potential-flow panel code. The aerodynamic forces were included through harmonic linearization, and the mooring lines were represented by a stiffness matrix. The frequency-domain code was benchmarked against an equivalent Bladed (DNV-GL AS, 2016) model with Morison-based hydrodynamics, and with a stiffness mooring matrix. Neither the frequency-domain model nor the Bladed model included viscous drag. Results were shown for regular waves and uniform, harmonic wind, and the frequency-domain code was reported to be up to 37 times faster than the Bladed model. In Lemmer et al. (2016), simplified time-domain models of the OC3-Hywind spar (Jonkman, 2010) and OC4-DeepCwind semi-submersible (Robertson et al., 2014) floating wind turbines were introduced. The models had four DoFs: floater surge and pitch, tower first fore-aft mode, and rotor azimuthal position. Linear hydrodynamics from a radiation-diffraction panel code were included in the time-domain model through the Cummins equation (Cummins, 1962). Aerodynamics were computed by coupling the code to AeroDyn. Quasi-static mooring forces were computed by solving the catenary mooring equations at each time step. A linearized version of the code was also presented. In the results, the linearized frequency-domain version was successfully benchmarked against the nonlinear time-domain version, by comparing the linear transfer function from wave height to tower-top displacement with its nonlinear equivalent. The work of Wang et al. (2017) involved a frequency-domain model of the DeepCwind semi-submersible (Robertson et al., 2014) with two rigid-body DoFs: floater surge and pitch. Linear hydrodynamics, linearized drag and drift forces were computed with the commercial software AQWA. The aerodynamic loads were included through a linearized version of the actuator point equation, where the aerodynamic contribution was divided into a constant force and a damping term – thus neglecting stochastic wind forcing. The mooring loads were included through a stiffness matrix, obtained from both quasi-static and dynamic mooring models. The model was validated against DeepCwind test data in terms of natural frequencies, response-amplitude operators and power spectral density (PSD) plots of surge and pitch response, generally obtaining a good agreement. However, a frequency-domain model for floating wind turbines able to incorporate realistic aerodynamic loads is still needed.

For bottom-fixed offshore wind turbines, Schløer et al. (2018) recently developed an efficient, frequency-domain model named QuLA (Quick Load Analysis), considering the DTU 10 MW Reference Wind Turbine (RWT; Bak et al., 2013). The mono-pile foundation and the wind turbine tower were defined as an Euler beam, and the first fore-aft modal deflection of this beam was the only DoF. Inspired by the work of van der Tempel (2006), the rotor and nacelle were represented by a point mass at the tower top, and aerodynamic loads and damping coefficients were pre-computed in the time-domain aero-elastic tool Flex5 (Øye, 1996). Compared to the work of van der Tempel (2006), the aerodynamic damping in QuLA was considered as dependent on mean wind speed. Hydrodynamic forcing was included through the Morison equation (Morison et al., 1950), where the structure velocity and acceleration were neglected. The code was validated against Flex5 in terms of time series, PSD, exceedance probability curves and fatigue damage-equivalent load (DEL). The bending moment at the seabed was estimated by QuLA within a 5 % error, and the code was reported to be approximately 40 times faster than its Flex5 equivalent.

This study presents the extension of QuLA to floating offshore wind turbines. The resulting model, QuLAF (Quick Load Analysis of Floating wind turbines), was first presented in Pegalajar-Jurado et al. (2016), with only two DoFs: floater surge and tower first fore-aft bending mode. Here we present an improved version of the model, a frequency-domain code that captures the four dominant DoFs in the in-plane global motion: floater surge, heave and pitch, and tower first fore-aft modal deflection. The model, which is here adapted to the DTU 10 MW RWT mounted on the OO-Star Wind Floater Semi 10 MW (Yu et al., 2018) was set up through cascading techniques. In the cascading process, information is pre-computed or extracted from more advanced models (parent models) to enhance the simplified models (children models). In this case, the hydrodynamic loads are extracted from the radiation-diffraction potential-flow solver WAMIT (Lee and Newman, 2016). The aerodynamic loads and aerodynamic damping coefficients are pre-computed in the numerical tool FAST v8 (Jonkman and Jonkman, 2016), and the mooring module MoorDyn (Hall, 2017) is employed to extract a mooring stiffness matrix for different operating positions. This way, the model includes standard radiation-diffraction theory and realistic rotor loads through pre-computed aero-elastic simulations. In the model, the system response is obtained by solving the linear equations of motion (EoMs) in the frequency domain, leading to a very efficient tool. While the radiation-diffraction results allow for a full linear response evaluation for rigid structure motion in waves, the ambition of this model is to extend them with the flexible tower and realistic stochastic rotor loads, thus going one step further than other simplified models in the literature. The results from QuLAF are here benchmarked against results from its time-domain, state-of-the-art (SoA) parent model in terms of time series, PSD, exceedance probability and fatigue DEL. We assess the strengths and weaknesses of the cascading process by comparison with the original time-domain model, and further develop techniques to improve the accuracy of the simplified model with respect to planar motion and tower-base loads. In this way, the potential of the model as a reliable tool for pre-design and as a complement to SoA models is demonstrated. The idea is that once the conceptual floater design is established with the efficient pre-design model, more advanced SoA models can be used for further design verification with a full design load basis that includes extreme and transient events.

Figure 1 The OO-Star Wind Floater Semi 10 MW concept (http://www.olavolsen.no; last access: 24 November 2017).

Figure 2 (a) Surge response of the OC3-Hywind floating wind turbine (Jonkman, 2010) to stochastic hydrodynamic linear forcing in a one-DoF linear model, computed in both the time and frequency domains. (b) Relative error between the time- and frequency-domain solutions. Only the first 1000 s are shown.

2 The case study

The floating wind turbine chosen for the present study is the DTU 10 MW RWT (Bak et al., 2013) mounted on the OO-Star Wind Floater Semi 10 MW (Yu et al., 2018). The main properties of the DTU 10 MW RWT are given in Table 1, and further information can be found in Bak et al. (2013). The basic DTU Wind Energy controller (Hansen and Henriksen, 2013) is utilized, and tuned to avoid the floater pitch instability (commonly known as the “negative damping problem”) reported in, for example, Larsen and Hanson (2007).

The floating substructure (see Fig. 1), developed by Dr.techn. Olav Olsen AS (http://www.olavolsen.no; last access: 24 November 2017), is a semi-submersible floater made of post-tensioned concrete. It has a central column and three outer columns mounted on a star-shaped pontoon with three legs. Each outer column is connected to the seabed by a catenary mooring line with a suspended clump weight. The main properties of the floating substructure are collected in Table 2, and further information can be found in Yu et al. (2018).

Table 1 Key figures for the DTU 10 MW RWT.

Table 2 Key figures for the OO-Star Wind Floater Semi 10 MW anchored at the selected site.

3 Time vs. frequency domain: advantages and disadvantages

Floating wind turbines can be considered harmonic oscillators with multiple coupled DoFs. To illustrate the strengths and weaknesses of solving the relevant EoM in the time or the frequency domain, a simple one-DoF mass-spring-damper system is considered,

$$\begin{array}{}\text{(1)}& m\ddot{\mathit{\xi }}\left(t\right)+b\dot{\mathit{\xi }}\left(t\right)+c\mathit{\xi }\left(t\right)=F\left(t\right),\end{array}$$

where m is the system mass, b is the damping coefficient, c is the restoring coefficient, ξ(t) is the system displacement from its equilibrium position, and F(t) is a harmonic excitation force. Equation (1) can be also written in complex notation, by expressing the excitation force at the frequency ω as $F\left(t\right)=\mathrm{\Re }\mathit{\left\{}\widehat{F}\right(\mathit{\omega }\left)e^{i\mathit{\omega }t}\mathit{\right\}}$, where ℜ indicates the real part, $\widehat{F}\left(\mathit{\omega }\right)$ is the Fourier transform of F(t) and i is the imaginary unit. If the initial transient part of the response is neglected, the steady-state system response at the given frequency can also be written as $\mathit{\xi }\left(t\right)=\mathrm{\Re }\mathit{\left\{}\widehat{\mathit{\xi }}\right(\mathit{\omega }\left)e^{i\mathit{\omega }t}\mathit{\right\}}$, leading to the EoM in the frequency domain,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }[-{\mathit{\omega }}^{\mathrm{2}}m+i\mathit{\omega }b+c]\widehat{\mathit{\xi }}\left(\mathit{\omega }\right)=\widehat{F}\left(\mathit{\omega }\right),\\ \text{(2)}& {\displaystyle }& {\displaystyle }\Rightarrow \widehat{\mathit{\xi }}\left(\mathit{\omega }\right)={\displaystyle \frac{\widehat{F}\left(\mathit{\omega }\right)}{-{\mathit{\omega }}^{\mathrm{2}}m+i\mathit{\omega }b+c}}\equiv H\left(\mathit{\omega }\right)\widehat{F}\left(\mathit{\omega }\right).\end{array}$$

The frequency-domain response $\widehat{\mathit{\xi }}\left(\mathit{\omega }\right)$ may be obtained by simply multiplying the frequency-domain excitation force $\widehat{F}\left(\mathit{\omega }\right)$ by the transfer function H(ω). This can be done at all frequencies and, due to the linearity, one can add the results at each frequency to obtain the total solution. Thus, once $\widehat{\mathit{\xi }}\left(\mathit{\omega }\right)$ has been determined for all frequencies, the time-domain response ξ(t) is obtained through an inverse Fourier transform of $\widehat{\mathit{\xi }}\left(\mathit{\omega }\right)$. If fast Fourier transform (FFT) and inverse fast Fourier transform (iFFT) are used, the solution can be obtained very quickly.

Figure 2 shows the surge response ξ(t) of a one-DoF model of the OC3-Hywind spar floating wind turbine (Jonkman, 2010) subjected to stochastic hydrodynamic linear forcing. The response labelled as “Time domain” was obtained by time-stepping of Eq. (1) with the classical fourth-order Runge–Kutta method and initial conditions ξ(0)=0 and $\dot{\mathit{\xi }}\left(\mathrm{0}\right)=\mathrm{0}$. The response labelled as “Frequency domain” was computed by first obtaining the frequency-domain excitation force $\widehat{F}\left(\mathit{\omega }\right)=\mathrm{FFT}\left(F\right(t\left)\right)$, calculating the frequency-domain response using Eq. (2) and finally going back to the time-domain response $\mathit{\xi }\left(t\right)=\mathrm{\Re }\mathit{\left\{}\mathrm{iFFT}\right(\widehat{\mathit{\xi }}\left(\mathit{\omega }\right)\left)\mathit{\right\}}$. The simulation time step was 0.025 s and the total simulated time was 5400 s, although only the first 1000 s are shown here. The time-domain solution took 69.41 s to run, while the frequency-domain solution was done in 0.03 s, or 2344 times faster. The two responses diverge at the beginning, where the time-domain solution is dominated by the transient response to the initial conditions, which is not present in the frequency-domain solution. However, after approximately 800 s (or six natural periods) and until the end of the simulation, the two solutions are practically identical, with errors between 0.2 % and 0.5 %, likely due to the time and frequency discretizations.

In addition to the gain in CPU speed, solving the EoM in the frequency domain allows for an easier handling of frequency-dependent properties, such as hydrodynamic added mass and radiation damping. On the other hand, it has also been shown that transient response due to initial conditions is only captured by time-domain models. However, as in the above example, the transient response due to initial conditions is an artifact of the time-domain formulation and is often discarded in the analysis. Perhaps the most clear disadvantage of frequency-domain models is that they can only accommodate loads that depend linearly on the response and its time derivatives, such as hydrodynamic added mass loads or hydrostatic loads. They cannot directly accommodate loads that depend on the response in a nonlinear manner, such as viscous drag from relative structural motion or catenary mooring loads. In those cases, simplified or linearized formulations have to be implemented instead.

Figure 3 Sketch of the floating wind turbine as seen by the QuLAF model.

4 The time-domain, state-of-the-art numerical model

In SoA models the nacelle, hub and floater are often considered rigid, whereas the tower and blades are flexible. The floater motion typically has six DoFs: surge, sway, heave, roll, pitch and yaw. Aerodynamics are normally computed using unsteady blade element momentum theory (Hansen, 2008). Hydrodynamics are typically represented by radiation-diffraction theory (Newman, 1980), the Morison equation or a combination of both. The mooring lines can be modelled with either quasi-static or dynamic approaches. In general, SoA models are more accurate than simplified models, but they also have a higher CPU cost.

A SoA, time-domain numerical model of the OO-Star Semi + DTU 10 MW floating wind turbine was used in this study as a parent model to QuLAF. The SoA model was implemented in FAST v8.16.00a-bjj (Jonkman and Jonkman, 2016) with active control and 15 DoFs for turbine and floater: first and second flap-wise blade modal deflections; first edge-wise blade modal deflection; drivetrain rotational flexibility; drivetrain speed; first and second fore-aft and side-side tower modal deflections; and floater surge, sway, heave, roll, pitch and yaw. The turbulent wind fields were computed in TurbSim, and the aerodynamic loads were modelled with AeroDyn v14. The basic DTU Wind Energy controller was applied through a dynamic-link library (DLL). The mooring loads, calculated by MoorDyn (Hall, 2017), included buoyancy, mass inertia and hydrodynamic loads resulting from the motion of the mooring lines in calm water. Hydrodynamic loads on the floater were first computed in WAMIT (Lee and Newman, 2016) and coupled to FAST through the Cummins equation. Viscous effects were modelled internally by the Morison drag term. Further details on the modelling of floating wind turbines in FAST can be found in Jonkman (2009), while a thorough description of the FAST model used in this study is presented in Pegalajar-Jurado et al. (2018b) and Pegalajar-Jurado et al. (2018a).

5 The frequency-domain, cascaded numerical model

The simplest model for the dynamic analysis of floating wind turbines would only have a few DoFs, typically rigid-body motion of the floater in surge and pitch. Aerodynamic loads would be represented by a point force at the rotor hub and defined by an actuator point model. If the floating substructure is slender compared to the incident waves, a strip-theory approach may be applied to compute the hydrodynamic loads from the Morison equation. The forces exerted by the mooring system can be included through a stiffness matrix in the linear EoMs. Simplified, low-order models are very CPU efficient but their accuracy is often limited. In the following we present a simplified model that combines elements extracted from a SoA model into a very efficient tool, which aims at getting close to the accuracy of the SoA model while still retaining the CPU efficiency of low-order models.

QuLAF represents the floating wind turbine as two lumped masses – floater and rotor-nacelle assembly – connected by a flexible tower. The model captures four planar DoFs – floater surge, heave, pitch and first tower fore-aft modal deflection – and is thus applicable to aligned wind and wave situations. The floating wind turbine is represented as depicted in Fig. 3. The EoM is a matrix version of Eq. (2),

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\left[-{\mathit{\omega }}^{\mathrm{2}}(\mathbf{M}+\mathbf{A}(\mathit{\omega }\left)\right)+i\mathit{\omega }\mathbf{B}\left(\mathit{\omega }\right)+\mathbf{C}\right]\widehat{\mathit{\xi }}\left(\mathit{\omega }\right)=\widehat{\mathit{F}}\left(\mathit{\omega }\right),\\ \text{(3)}& {\displaystyle }& {\displaystyle }\Rightarrow \widehat{\mathit{\xi }}\left(\mathit{\omega }\right)=\mathbf{H}\left(\mathit{\omega }\right)\widehat{\mathit{F}}\left(\mathit{\omega }\right),\end{array}$$

where M is the structural mass and inertia matrix, A(ω) is the frequency-dependent, hydrodynamic added mass and inertia matrix, B(ω) is the frequency-dependent damping matrix, and C is the restoring matrix. The vector $\widehat{\mathit{\xi }}\left(\mathit{\omega }\right)$ is the dynamic response in the frequency domain for the four DoFs and $\widehat{\mathit{F}}\left(\mathit{\omega }\right)$ is the dynamic vector of excitation forces and moments in the frequency domain. The system transfer function is given by H(ω). The different elements in Eq. (3) are described in detail below.

5.1 Dynamic response vector

The dynamic response vector,

$$\begin{array}{}\text{(4)}& \widehat{\mathit{\xi }}\left(\mathit{\omega }\right)=\left[\begin{array}{c}{\widehat{\mathit{\xi }}}_{\mathrm{1}}\left(\mathit{\omega }\right)\\ {\widehat{\mathit{\xi }}}_{\mathrm{3}}\left(\mathit{\omega }\right)\\ {\widehat{\mathit{\xi }}}_{\mathrm{5}}\left(\mathit{\omega }\right)\\ \widehat{\mathit{\alpha }}\left(\mathit{\omega }\right)\end{array}\right],\end{array}$$

has one element for each DoF: floater surge, heave, pitch and first tower fore-aft modal deflection. The sign convention is that shown in Fig. 3, with positive surge in the downwind direction, positive heave upwards, positive pitch (about flotation point O) clockwise and positive tower deflection in the downwind direction. The physical tower deflection at any height z can be obtained by multiplying the mode shape ϕ(z) and the modal deflection α(t). The tower deflection at the hub height h_hub is therefore given by δ(t)=ϕ_hubα(t). If the absolute nacelle displacement is sought, the contributions from floater surge and pitch motions must be added to the tower deflection, and the global response vector ${\widehat{\mathit{\xi }}}_{\mathrm{glob}}\left(\mathit{\omega }\right)$ is found by introducing a transformation matrix T_glob,

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }{\widehat{\mathit{\xi }}}_{\mathrm{glob}}\left(\mathit{\omega }\right)=\left[\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& h_{\mathrm{hub}}& {\mathit{\varphi }}_{\mathrm{hub}}\end{array}\right]\left[\begin{array}{c}{\widehat{\mathit{\xi }}}_{\mathrm{1}}\left(\mathit{\omega }\right)\\ {\widehat{\mathit{\xi }}}_{\mathrm{3}}\left(\mathit{\omega }\right)\\ {\widehat{\mathit{\xi }}}_{\mathrm{5}}\left(\mathit{\omega }\right)\\ \widehat{\mathit{\alpha }}\left(\mathit{\omega }\right)\end{array}\right]={\mathbf{T}}_{\mathrm{glob}}\widehat{\mathit{\xi }}\left(\mathit{\omega }\right).\end{array}$$

5.2 Dynamic load vector

The dynamic load vector,

$$\begin{array}{}\text{(6)}& \widehat{\mathit{F}}\left(\mathit{\omega }\right)={\widehat{\mathit{F}}}_{\mathrm{hydro}}\left(\mathit{\omega }\right)+{\widehat{\mathit{F}}}_{\mathrm{aero}}\left(\mathit{\omega }\right),\end{array}$$

contains hydrodynamic loads ${\widehat{\mathit{F}}}_{\mathrm{hydro}}\left(\mathit{\omega }\right)$ and aerodynamic loads ${\widehat{\mathit{F}}}_{\mathrm{aero}}\left(\mathit{\omega }\right)$. Hydrodynamic loads are extracted from the solution to the diffraction problem, which provides a vector of wave excitation forces and moments in all six DoFs, namely $\widehat{\mathit{X}}\left(\mathit{\omega }\right)$. These excitation forces are normalized to waves of unit amplitude, therefore the wave loads for a specific time series of free-surface elevation η(t) are obtained by the product $\widehat{\mathit{X}}\left(\mathit{\omega }\right)\widehat{\mathit{\eta }}\left(\mathit{\omega }\right)$. The vector of wave excitation forces and moments is reduced to adapt it to the simplified model, and a zero is added in the fourth element for the tower DoF,

$$\begin{array}{}\text{(7)}& {\widehat{\mathit{F}}}_{\mathrm{hydro}}\left(\mathit{\omega }\right)=\widehat{\mathit{X}}\left(\mathit{\omega }\right)\widehat{\mathit{\eta }}\left(\mathit{\omega }\right)\equiv \left[\begin{array}{c}{\widehat{X}}_{\mathrm{1}}\left(\mathit{\omega }\right)\\ {\widehat{X}}_{\mathrm{3}}\left(\mathit{\omega }\right)\\ {\widehat{X}}_{\mathrm{5}}\left(\mathit{\omega }\right)\\ \mathrm{0}\end{array}\right]\widehat{\mathit{\eta }}\left(\mathit{\omega }\right),\end{array}$$

where $\widehat{\mathit{\eta }}\left(\mathit{\omega }\right)$ can be computed from an input time series η(t) or from a theoretical wave spectrum. The only viscous effect considered in the model is viscous damping (see Sect. 5.8.1), but viscous forcing is neglected to keep the model computationally efficient. This simplification, however, is considered reasonable because hydrodynamics for this floater are dominated by inertia loads, and viscous forcing is expected to be relevant mainly for severe sea states, which lie on the border of the model's applicability. The vector of aerodynamic loads only contains the dynamic part of the wind loads and has the format

$$\begin{array}{}\text{(8)}& {\widehat{\mathit{F}}}_{\mathrm{aero}}\left(\mathit{\omega }\right)=\left[\begin{array}{c}{\widehat{F}}_{\mathrm{aero},\mathrm{1}}\left(\mathit{\omega }\right)\\ {\widehat{F}}_{\mathrm{aero},\mathrm{3}}\left(\mathit{\omega }\right)\\ {\widehat{F}}_{\mathrm{aero},\mathrm{1}}\left(\mathit{\omega }\right)h_{\mathrm{hub}}+{\widehat{\mathit{\tau }}}_{\mathrm{aero}}\left(\mathit{\omega }\right)\\ {\widehat{F}}_{\mathrm{aero},\mathrm{1}}\left(\mathit{\omega }\right){\mathit{\varphi }}_{\mathrm{hub}}+{\widehat{\mathit{\tau }}}_{\mathrm{aero}}\left(\mathit{\omega }\right){\mathit{\varphi }}_{z,\mathrm{hub}}\end{array}\right],\end{array}$$

where ${\widehat{F}}_{\mathrm{aero},\mathrm{1}}\left(\mathit{\omega }\right)$ and ${\widehat{F}}_{\mathrm{aero},\mathrm{3}}\left(\mathit{\omega }\right)$ represent the horizontal and vertical components of the aerodynamic loads on the rotor, respectively. The aerodynamic tilt torque on the rotor is given by ${\widehat{\mathit{\tau }}}_{\mathrm{aero}}\left(\mathit{\omega }\right)$. The fourth element of ${\widehat{\mathit{F}}}_{\mathrm{aero}}$ represents the effect of the aerodynamic loads on the tower modal deflection, hence the mode shape deflection ϕ_hub and its slope ϕ_z,hub evaluated at the hub are involved. The time-domain aerodynamic loads for each mean wind speed W are pre-computed in the SoA model, as detailed in Sect. 5.8.2.

5.3 Structural mass and inertia matrix

The symmetric matrix of structural mass and inertia, obtained by looking at the forces needed to produce unit accelerations in the different DoFs, is defined as

$$\begin{array}{}\text{(9)}& \mathbf{M}=\left[\begin{array}{cccc}{m}_{\mathrm{tot}}& \mathrm{0}& m_{\mathrm{tot}}{z}_{\mathrm{tot}}^{\mathrm{CM}}& m_{\mathrm{rn}}{\mathit{\varphi }}_{\mathrm{hub}}+\sum _{i=\mathrm{1}}^{N_t}{\stackrel{\mathrm{̃}}{\mathit{\rho }}}_i{\mathit{\varphi }}_i\mathrm{\Delta }{z}_i\\ & m_{\mathrm{tot}}& \mathrm{0}& \mathrm{0}\\ & & I_{\mathrm{tot}}^O& m_{\mathrm{rn}}{\mathit{\varphi }}_{\mathrm{hub}}{h}_{\mathrm{hub}}+I_{\mathrm{rn}}^{\mathrm{TT}}{\mathit{\varphi }}_{z,\mathrm{hub}}\\ & & & +\sum _{i=\mathrm{1}}^{N_t}{\stackrel{\mathrm{̃}}{\mathit{\rho }}}_i{\mathit{\varphi }}_i{z}_i\mathrm{\Delta }{z}_i\\ & & & m_{\mathrm{rn}}{\mathit{\varphi }}_{\mathrm{hub}}^{\mathrm{2}}+I_{\mathrm{rn}}^{\mathrm{TT}}{\mathit{\varphi }}_{z,\mathrm{hub}}^{\mathrm{2}}\\ & & & +\sum _{i=\mathrm{1}}^{N_t}{\stackrel{\mathrm{̃}}{\mathit{\rho }}}_i{\mathit{\varphi }}_i^{\mathrm{2}}\mathrm{\Delta }{z}_i\end{array}\right],\end{array}$$

where m_tot is the total mass of the floating wind turbine, $m_{\mathrm{tot}}=m_{\mathrm{f}}+m_{\mathrm{rn}}+{\sum }_{i=\mathrm{1}}^{N_t}{\stackrel{\mathrm{̃}}{\mathit{\rho }}}_i\mathrm{\Delta }{z}_i$, which includes the mass of the floater m_f, the rotor-nacelle mass m_rn and the mass sum of all the N_t elements that compose the flexible tower, each with a mass per length ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}_i$ and a height Δz_i. The total mass inertia of the system about the y axis at the flotation point O is given by $I_{\mathrm{tot}}^O=I_{\mathrm{f}}^O+I_{\mathrm{rn}}^O+\sum _{i=\mathrm{1}}^{N_t}{\stackrel{\mathrm{̃}}{\mathit{\rho }}}_i{z}_i^{\mathrm{2}}\mathrm{\Delta }{z}_i$, including the floater inertia $I_{\mathrm{f}}^O$, the rotor-nacelle inertia $I_{\mathrm{rn}}^O$ and the inertia of each of the tower elements, located at an absolute height $z_i=(z_{t,i}+h_t)$, where z_t,i is the vertical position of the element i with respect to the tower base, located at a height h_t. The centre of mass (CM) of the whole structure is located at $z_{\mathrm{tot}}^{\mathrm{CM}}=(m_{\mathrm{f}}{z}_{\mathrm{f}}^{\mathrm{CM}}+m_{\mathrm{rn}}{h}_{\mathrm{hub}}+{\sum }_{i=\mathrm{1}}^{N_t}{\stackrel{\mathrm{̃}}{\mathit{\rho }}}_i{z}_i\mathrm{\Delta }{z}_i)/m_{\mathrm{tot}}$, with contributions from the floater CM at $z_{\mathrm{f}}^{\mathrm{CM}}$, the rotor-nacelle CM at the hub height h_hub and the CM of each of the tower elements. The mode shape deflection of the tower evaluated at a generic tower element i is ϕ_i, while ϕ_hub and ϕ_z,hub are the mode shape deflection and its slope evaluated at the hub. Finally, $I_{\mathrm{rn}}^{\mathrm{TT}}$ represents the mass inertia of the rotor-nacelle assembly referred to the tower top. The tower structural properties and first mode shape are the same as the ones given as an input to the SoA model.

5.4 Hydrodynamic added mass matrix and damping matrix

The frequency-dependent, hydrodynamic added-mass and radiation-damping matrices, A(ω) and B_rad(ω), can be pre-computed in a radiation-diffraction solver. Here, the same WAMIT output files used for the SoA model are loaded into QuLAF. However, the original 6×6 matrices are reduced by removing the rows and columns corresponding to the DoFs not included in the simplified model (sway, roll, yaw), and a row and column of zeros is added for compatibility with the tower DoF,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathbf{A}\left(\mathit{\omega }\right)=\left[\begin{array}{cccc}{a}_{\mathrm{11}}\left(\mathit{\omega }\right)& a_{\mathrm{13}}\left(\mathit{\omega }\right)& a_{\mathrm{15}}\left(\mathit{\omega }\right)& \mathrm{0}\\ a_{\mathrm{31}}\left(\mathit{\omega }\right)& a_{\mathrm{33}}\left(\mathit{\omega }\right)& a_{\mathrm{35}}\left(\mathit{\omega }\right)& \mathrm{0}\\ a_{\mathrm{51}}\left(\mathit{\omega }\right)& a_{\mathrm{53}}\left(\mathit{\omega }\right)& a_{\mathrm{55}}\left(\mathit{\omega }\right)& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right],\\ \text{(10)}& {\displaystyle }& {\displaystyle }{\mathbf{B}}_{\mathrm{rad}}\left(\mathit{\omega }\right)=\left[\begin{array}{cccc}{b}_{\mathrm{11}}\left(\mathit{\omega }\right)& b_{\mathrm{13}}\left(\mathit{\omega }\right)& b_{\mathrm{15}}\left(\mathit{\omega }\right)& \mathrm{0}\\ b_{\mathrm{31}}\left(\mathit{\omega }\right)& b_{\mathrm{33}}\left(\mathit{\omega }\right)& b_{\mathrm{35}}\left(\mathit{\omega }\right)& \mathrm{0}\\ b_{\mathrm{51}}\left(\mathit{\omega }\right)& b_{\mathrm{53}}\left(\mathit{\omega }\right)& b_{\mathrm{55}}\left(\mathit{\omega }\right)& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right].\end{array}$$

The global damping matrix includes contributions from the hydrodynamic radiation damping B_rad(ω), the hydrodynamic viscous damping B_vis, the aerodynamic damping B_aero(ω) and the tower structural damping B_struc:

$$\begin{array}{}\text{(11)}& \mathbf{B}\left(\mathit{\omega }\right)={\mathbf{B}}_{\mathrm{rad}}\left(\mathit{\omega }\right)+{\mathbf{B}}_{\mathrm{vis}}+{\mathbf{B}}_{\mathrm{aero}}\left(\mathit{\omega }\right)+{\mathbf{B}}_{\mathrm{struc}}.\end{array}$$

The hydrodynamic viscous damping matrix B_vis is analytically extracted from the Morison equation, as shown in Section 5.8.1. The diagonal matrix of aerodynamic damping,

$$\begin{array}{}\text{(12)}& {\mathbf{B}}_{\mathrm{aero}}\left(\mathit{\omega }\right)=\left[\begin{array}{cccc}{b}_{\mathrm{aero},\mathrm{11}}\left(\mathit{\omega }\right)& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ & \mathrm{0}& \mathrm{0}& \mathrm{0}\\ & & b_{\mathrm{aero},\mathrm{55}}\left(\mathit{\omega }\right)& \mathrm{0}\\ & & & b_{\mathrm{aero},\mathrm{tow}}\end{array}\right],\end{array}$$

is extracted from the SoA model for each mean wind speed W, as detailed in Section 5.8.2. The matrix of structural damping only concerns the tower and is given by

$$\begin{array}{}\text{(13)}& {\mathbf{B}}_{\mathrm{struc}}=\left[\begin{array}{cccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ & \mathrm{0}& \mathrm{0}& \mathrm{0}\\ & & \mathrm{0}& \mathrm{0}\\ & & & \mathrm{2}{\mathit{\zeta }}_{\mathrm{struc},\mathrm{tow}}\sqrt{C_{\mathrm{tow}}{M}_{\mathrm{tow}}}\end{array}\right],\end{array}$$

where the structural damping ratio for the first fore-aft tower mode, ζ_struc,tow, is directly taken from the input to the SoA model, and C_tow and M_tow are the last diagonal elements of the system restoring matrix C and the mass inertia matrix M, respectively.

5.5 Restoring matrix

The restoring matrix includes hydrostatic stiffness C_hst, structural stiffness C_struc and mooring stiffness C_moor:

$$\begin{array}{}\text{(14)}& \mathbf{C}={\mathbf{C}}_{\mathrm{hst}}+{\mathbf{C}}_{\mathrm{struc}}+{\mathbf{C}}_{\mathrm{moor}}.\end{array}$$

The hydrostatic matrix should only include the contributions from the centre of buoyancy (CB) and waterplane area. It is computed as part of the radiation-diffraction solution, and is reduced following the same procedure as for the added mass and radiation damping matrices. The symmetric matrix of structural stiffness is given by

$$\begin{array}{}\text{(15)}& {\mathbf{C}}_{\mathrm{struc}}=\left[\begin{array}{cccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ & \mathrm{0}& \mathrm{0}& \mathrm{0}\\ & & -m_{\mathrm{tot}}g{z}_{\mathrm{tot}}^{\mathrm{CM}}& -m_{\mathrm{rn}}g{\mathit{\varphi }}_{\mathrm{hub}}-\sum _{i=\mathrm{1}}^{N_t}{\stackrel{\mathrm{̃}}{\mathit{\rho }}}_{i}g{\mathit{\varphi }}_i\mathrm{\Delta }{z}_i\\ & & & \sum _{i=\mathrm{1}}^{N_t}{\mathrm{EI}}_i{\mathit{\varphi }}_{zz,i}^{\mathrm{2}}\mathrm{\Delta }{z}_i\end{array}\right],\end{array}$$

where g is the acceleration of gravity, and EI_i and ϕ_zz,i are the bending stiffness and the curvature of the mode shape for the tower element i, respectively. The off-diagonal term represents the negative restoring effect of the tower and rotor-nacelle mass on the tower DoF when the floater pitches. The mooring restoring matrix C_moor is position-dependent and therefore extracted from the SoA model for each mean wind speed W, as detailed in Sect. 5.8.3. Although in this study wind is the only effect considered to affect the mean position of the floating wind turbine, other effects such as mean drift forces and current can be taken into account in the SoA model when linearizing the mooring system.

5.6 Static load and response

Static loads are related to the equilibrium of the structure. In the model, the static part of the response, ξ_st, is added to the dynamic part $\widehat{\mathit{\xi }}\left(\mathit{\omega }\right)$ when it is converted from the frequency to the time domain via iFFT. The static loads applied are

$$\begin{array}{}\text{(16)}& {\mathit{F}}_{\mathrm{st}}={\mathit{F}}_{\mathrm{aero},\mathrm{st}}+{\mathit{F}}_{\mathrm{grav}}+{\mathit{F}}_{\mathrm{buoy}},\end{array}$$

which include the static part of the aerodynamic loads F_aero,st, the gravity loads F_grav and the buoyancy loads F_buoy. The gravity load vector is given by

$$\begin{array}{}\text{(17)}& {\mathit{F}}_{\mathrm{grav}}=\left[\begin{array}{c}\mathrm{0}\\ -m_{\mathrm{tot}}g-F_{\mathrm{moor},z}\\ m_{\mathrm{rn}}g{x}_{\mathrm{rn}}^{\mathrm{CM}}\\ m_{\mathrm{rn}}g{x}_{\mathrm{rn}}^{\mathrm{CM}}{\mathit{\varphi }}_{z,\mathrm{hub}}\end{array}\right],\end{array}$$

where F_moor,z is the vertical force exerted by the mooring lines in equilibrium, and $x_{\mathrm{rn}}^{\mathrm{CM}}$ is the horizontal coordinate of the rotor-nacelle CM.

The buoyancy load vector is

$$\begin{array}{}\text{(18)}& {\mathit{F}}_{\mathrm{buoy}}=\left[\begin{array}{c}\mathrm{0}\\ {\mathit{\rho }}_{\mathrm{w}}g{V}_{\mathrm{f}}\\ -{\mathit{\rho }}_{\mathrm{w}}g{V}_{\mathrm{f}}{x}_{\mathrm{f}}^{\mathrm{CB}}\\ \mathrm{0}\end{array}\right],\end{array}$$

where ρ_w is the water density, V_f is the volume displaced by the floater and $x_f^{\mathrm{CB}}$ is the horizontal coordinate of the floater CB. With no wind and only linear wave forcing, the floating wind turbine operates around its equilibrium position with a stiffness matrix C₀. If wind (or any other mean force) is introduced, the floating wind turbine is moved to a new equilibrium position, where the stiffness matrix is C_W. The static response ξ_st is therefore obtained from the static loads by considering a mean stiffness matrix C_st,

$$\begin{array}{}\text{(19)}& {\mathbf{C}}_{\mathrm{st}}={\displaystyle \frac{{\mathbf{C}}_{\mathrm{0}}+{\mathbf{C}}_{\mathrm{W}}}{\mathrm{2}}}\Rightarrow {\mathbf{C}}_{\mathrm{st}}{\mathit{\xi }}_{\mathrm{st}}={\mathit{F}}_{\mathrm{st}}.\end{array}$$

This approximation is accurate to second order.

5.7 System natural frequencies

The vector of natural frequencies ω₀ is found by solving the undamped eigenvalue problem given by

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\left[-{{\mathit{\omega }}_{\mathbf{0}}}^{\mathrm{2}}(\mathbf{M}+\mathbf{A}({\mathit{\omega }}_{\mathbf{0}}\left)\right)+\mathbf{C}\right]\widehat{\mathit{\xi }}\left({\mathit{\omega }}_{\mathbf{0}}\right)=\mathrm{0},\\ \text{(20)}& {\displaystyle }& {\displaystyle }\Rightarrow {\mathit{\omega }}_{\mathbf{0}}^{\mathbf{2}}\widehat{\mathit{\xi }}\left({\mathit{\omega }}_{\mathbf{0}}\right)=(\mathbf{M}+\mathbf{A}({\mathit{\omega }}_{\mathbf{0}}\left){)}^{-\mathrm{1}}\mathbf{C}\widehat{\mathit{\xi }}\right({\mathit{\omega }}_{\mathbf{0}}).\end{array}$$

Since the matrix of added mass depends on frequency, the eigenvalue problem is solved in a frequency loop. For each frequency ω, the four possible natural frequencies are computed. When one of the four possible frequencies obtained is equal to the frequency of that particular iteration in the loop, then a system natural frequency has been found. The system natural frequencies computed in QuLAF are compared to those obtained with the SoA model in Sect. 6.1.

5.8 Cascading techniques applied to the simplified model

In Sect. 3 it was stated that one disadvantage of frequency-domain models is their inability to directly capture loads that depend on the response in a nonlinear way. Some relevant examples are viscous drag, aerodynamic loads and catenary mooring loads. This section gives a description of the cascading methods employed to incorporate such nonlinear loads into the simplified model.

5.8.1 Hydrodynamic viscous loads

Viscous effects on submerged bodies depend nonlinearly on the relative velocity between the wave particles and the structure, hence they can only be directly incorporated in time-domain models. In the offshore community this is normally done through the drag term of the Morison equation, which provides the transversal drag force F_d on a cylindrical member section of diameter D and length dl as

$$\begin{array}{}\text{(21)}& F_{\mathrm{d}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }{C}_{\mathrm{D}}D|v_{\mathrm{f}}-v_{\mathrm{s}}|(v_{\mathrm{f}}-v_{\mathrm{s}})\mathrm{d}l,\end{array}$$

where ρ is the fluid density, C_D is a drag coefficient, and v_f and v_s are the local fluid and structure velocities perpendicular to the member axis. The equation can be also written as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{F}_{\mathrm{d}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }{C}_{\mathrm{D}}D\mathrm{sgn}\left(v_{\mathrm{f}}-v_{\mathrm{s}}\right){\left(v_{\mathrm{f}}-v_{\mathrm{s}}\right)}^{\mathrm{2}}\mathrm{d}l\\ \text{(22)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }{C}_{\mathrm{D}}D\mathrm{sgn}\left(v_{\mathrm{f}}-v_{\mathrm{s}}\right)\left(v_{\mathrm{f}}^{\mathrm{2}}+v_{\mathrm{s}}^{\mathrm{2}}-\mathrm{2}{v}_{\mathrm{f}}{v}_{\mathrm{s}}\right)\mathrm{d}l,\end{array}$$

which shows that the drag effects can be separated into a pure forcing term, a nonlinear damping term and a linear damping term. Since the hydrodynamics on the given floating substructure are inertia-dominated and under the assumption of small displacements around the equilibrium position, the two first terms are neglected and only the linear damping term is retained in the QuLAF model. Invoking further the assumption of small displacements and velocities relative to the fluid velocity, we have $\mathrm{sgn}(v_{\mathrm{f}}-v_{\mathrm{s}})\approx \mathrm{sgn}\left(v_{\mathrm{f}}\right)$. With this assumption, the linear damping term of the viscous force becomes

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{F}_{{\mathrm{d}}_{\mathrm{l}}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }{C}_{\mathrm{D}}D\mathrm{sgn}\left(v_{\mathrm{f}}-v_{\mathrm{s}}\right)\left(-\mathrm{2}{v}_{\mathrm{f}}{v}_{\mathrm{s}}\right)\mathrm{d}l\\ \text{(23)}& {\displaystyle }& {\displaystyle }\approx -\mathit{\rho }{C}_{\mathrm{D}}D\left|v_{\mathrm{f}}\right|v_{\mathrm{s}}\mathrm{d}l.\end{array}$$

A symmetric viscous damping matrix B_vis is now derived by applying Eq. (23) to the different DoFs. For the surge motion, integration over the submerged body gives the total viscous force in the x direction as

$$\begin{array}{}\text{(24)}& F_{\mathrm{1}}=-\underset{z_{\min }}{\overset{\mathrm{0}}{\int }}\mathit{\rho }{C}_{\mathrm{D}}D\left|u\right|{\dot{\mathit{\xi }}}_{\mathrm{1}}\mathrm{d}z,\end{array}$$

where z_min is the structure's deepest submerged point, u is the horizontal wave particle velocity and ${\dot{\mathit{\xi }}}_{\mathrm{1}}$ is the surge velocity. The integral in Eq. (24) requires the estimation of drag coefficients and the computation of wave kinematics at several locations on the submerged structure, which can be involved for complex geometries. These computations would reduce the CPU efficiency relative to the radiation-diffraction terms, so instead the local drag coefficient and wave velocity inside the integral are replaced by global, representative values outside the integral, C_Dx and u_rep. Hereby the force becomes

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{F}_{\mathrm{1}}=-\mathit{\rho }{\dot{\mathit{\xi }}}_{\mathrm{1}}\underset{z_{\min }}{\overset{\mathrm{0}}{\int }}{C}_{\mathrm{D}}D\left|u\right|\mathrm{d}z\approx -\mathit{\rho }{C}_{\mathrm{D}x}{u}_{\mathrm{rep}}{\dot{\mathit{\xi }}}_{\mathrm{1}}\underset{z_{\min }}{\overset{\mathrm{0}}{\int }}D\mathrm{d}z\\ \text{(25)}& {\displaystyle }& {\displaystyle }=-\mathit{\rho }{C}_{\mathrm{D}x}{A}_x{u}_{\mathrm{rep}}{\dot{\mathit{\xi }}}_{\mathrm{1}}\equiv -b_{\mathrm{11}}{\dot{\mathit{\xi }}}_{\mathrm{1}},\end{array}$$

where A_x is the integral of the local diameter D over depth, or the floater's area projected on the yz plane. This defines the surge-surge element of the viscous damping matrix B_vis. Further the b₅₁ element of the matrix is obtained by consideration of the moment from F₁ around the point of flotation,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\mathrm{1}}=-\mathit{\rho }{\dot{\mathit{\xi }}}_{\mathrm{1}}\underset{z_{\min }}{\overset{\mathrm{0}}{\int }}{C}_{\mathrm{D}}D\left|u\right|z\mathrm{d}z\approx -\mathit{\rho }{C}_{\mathrm{D}x}{u}_{\mathrm{rep}}{\dot{\mathit{\xi }}}_{\mathrm{1}}\underset{z_{\min }}{\overset{\mathrm{0}}{\int }}Dz\mathrm{d}z\\ \text{(26)}& {\displaystyle }& {\displaystyle }=-\mathit{\rho }{C}_{\mathrm{D}x}{S}_{y,Ax}{u}_{\mathrm{rep}}{\dot{\mathit{\xi }}}_{\mathrm{1}}\equiv -b_{\mathrm{51}}{\dot{\mathit{\xi }}}_{\mathrm{1}},\end{array}$$

where $S_{y,A_x}$ is the first moment of area of A_x about the y axis (negative due to z≤0) and b₅₁ is the surge-pitch element of the viscous damping matrix. In a similar way, the heave-heave and heave-pitch coefficients of B_vis are obtained by applying Eq. (23) to the heave motion,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{F}_{\mathrm{3}}=-\mathit{\rho }{\dot{\mathit{\xi }}}_{\mathrm{3}}\underset{x_{\min }}{\overset{x_{\max }}{\int }}{C}_{\mathrm{D}}D\left|w\right|\mathrm{d}x\approx -\mathit{\rho }{C}_{Dz}{w}_{\mathrm{rep}}{\dot{\mathit{\xi }}}_{\mathrm{3}}\underset{x_{\min }}{\overset{x_{\max }}{\int }}D\mathrm{d}x\\ \text{(27)}& {\displaystyle }& {\displaystyle }=-\mathit{\rho }{C}_{Dz}{A}_z{w}_{\mathrm{rep}}{\dot{\mathit{\xi }}}_{\mathrm{3}}\equiv -b_{\mathrm{33}}{\dot{\mathit{\xi }}}_{\mathrm{3}},{\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\mathrm{3}}=\mathit{\rho }{\dot{\mathit{\xi }}}_{\mathrm{3}}\underset{x_{\min }}{\overset{x_{\max }}{\int }}{C}_{\mathrm{D}}D\left|w\right|x\mathrm{d}x\approx \mathit{\rho }{C}_{Dz}{w}_{\mathrm{rep}}{\dot{\mathit{\xi }}}_{\mathrm{3}}\underset{x_{\min }}{\overset{x_{\max }}{\int }}Dx\mathrm{d}x\\ \text{(28)}& {\displaystyle }& {\displaystyle }=\mathit{\rho }{C}_{\mathrm{D}z}{S}_{y,Az}{w}_{\mathrm{rep}}{\dot{\mathit{\xi }}}_{\mathrm{3}}\equiv -b_{\mathrm{53}}{\dot{\mathit{\xi }}}_{\mathrm{3}}.\end{array}$$

Here ${\dot{\mathit{\xi }}}_{\mathrm{3}}$ is the heave velocity, w is the wave particle vertical velocity, A_z is the floater's bottom area projected on the xy plane and $S_{y,A_z}$ is the first moment of area of A_z about the y axis, which is zero for the present floating substructure due to symmetry. Finally, by applying Eq. (23) to the pitch motion, the pitch-pitch element of the viscous damping matrix, b₅₅, is found. When the floater pitches with a velocity ${\dot{\mathit{\xi }}}_{\mathrm{5}}$, an arbitrary point on the floater with coordinates (x,z) moves with a velocity $(z{\dot{\mathit{\xi }}}_{\mathrm{5}},-x{\dot{\mathit{\xi }}}_{\mathrm{5}})$. The motion creates a moment due to viscous effects given by

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\mathrm{5}}=-\mathit{\rho }{\dot{\mathit{\xi }}}_{\mathrm{5}}\underset{z_{\min }}{\overset{\mathrm{0}}{\int }}{C}_{\mathrm{D}}D\left|u\right|z^{\mathrm{2}}\mathrm{d}z-\mathit{\rho }{\dot{\mathit{\xi }}}_{\mathrm{5}}\underset{x_{\min }}{\overset{x_{\max }}{\int }}{C}_{\mathrm{D}}D\left|w\right|x^{\mathrm{2}}\mathrm{d}x\\ \text{(29)}& {\displaystyle }& {\displaystyle }\approx -\mathit{\rho }(C_{\mathrm{D}x}{I}_{y,Ax}{u}_{\mathrm{rep}}+C_{\mathrm{D}z}{I}_{y,Az}{w}_{\mathrm{rep}}){\dot{\mathit{\xi }}}_{\mathrm{5}}\equiv -b_{\mathrm{55}}{\dot{\mathit{\xi }}}_{\mathrm{5}},\end{array}$$

where I_y,Ax and I_y,Az are the second moments of area of A_x and A_z about the y axis, respectively. The complete symmetric matrix of viscous damping is therefore

$$\begin{array}{}\text{(30)}& {\mathbf{B}}_{\mathrm{vis}}=\left[\begin{array}{cccc}\mathit{\rho }{C}_{\mathrm{D}x}{A}_x{u}_{\mathrm{rep}}& \mathrm{0}& \mathit{\rho }{C}_{\mathrm{D}x}{S}_{y,Ax}{u}_{\mathrm{rep}}& \mathrm{0}\\ & \mathit{\rho }{C}_{Dz}{A}_z{w}_{\mathrm{rep}}& \mathrm{0}& \mathrm{0}\\ & & \mathit{\rho }\left(C_{\mathrm{D}x}{I}_{y,Ax}{u}_{\mathrm{rep}}\right.& \mathrm{0}\\ & & \left.+C_{Dz}{I}_{y,Az}{w}_{\mathrm{rep}}\right)& \\ & & & \mathrm{0}\end{array}\right].\end{array}$$

The global drag coefficients above have been chosen as C_Dx=1 and C_Dz=2, given that the bottom slab of the floater under consideration has sharp corners and is expected to oppose a greater resistance to the flow than the smooth vertical columns (see Fig. 1). To obtain the representative velocity u_rep, the time- and depth-dependent horizontal wave velocity at the floater's centreline $u(\mathrm{0},z,t)$ is first averaged over depth and then over time,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{u}_{\mathrm{avg}}\left(t\right)={\displaystyle \frac{\mathrm{1}}{\left|z_{\min }\right|}}\underset{z_{\min }}{\overset{\mathrm{0}}{\int }}u(\mathrm{0},z,t)\mathrm{d}z\equiv \\ {\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathrm{1}}{\left|z_{\min }\right|}}\mathrm{\Re }\left\{\mathrm{iFFT}\left({\displaystyle \frac{\mathit{\omega }\widehat{\mathit{\eta }}\left(\mathit{\omega }\right)}{k}}\left(\mathrm{1}-{\displaystyle \frac{\sinh \left(k\right(z_{\min }+h\left)\right)}{\sinh \left(kh\right)}}\right)\right)\right\},\\ \text{(31)}& {\displaystyle }& {\displaystyle }\Rightarrow u_{\mathrm{rep}}=\stackrel{\mathrm{‾}}{\left|u_{\mathrm{avg}}\right|}.\end{array}$$

Here k is the wave number for the angular frequency ω and h is the water depth. The representative velocity w_rep is chosen as the time average of the vertical wave velocity at the centre of the bottom plate,

$$\begin{array}{}\text{(32)}& w_{\mathrm{avg}}\left(t\right)=w(\mathrm{0},z_{\min },t)\Rightarrow w_{\mathrm{rep}}=\stackrel{\mathrm{‾}}{\left|w_{\mathrm{avg}}\right|}.\end{array}$$

This simplification of the wave kinematics history, although drastic, allows for the characterization of the viscous damping for each sea state and avoids the need to compute wave kinematics locally and integrate the drag loads.

Figure 4 Example of time series of hub position and selected peaks for the extraction of aerodynamic damping. From top to bottom: surge, pitch and clamped tower DoFs.

5.8.2 Aerodynamic loads

Aerodynamic loads depend on the square of the relative wind speed seen by the blades. The relative wind speed includes contributions from the rotor speed, the blade deflection, the tower deflection and the motion of the floater. The fact that the aerodynamic thrust depends on the blade relative velocity produces the well-known aerodynamic damping. State-of-the-art numerical models incorporate aerodynamic loads based on relative velocity, because both the wind speed and the blade structural velocity are known at each time step. However, this cannot be done in a frequency-domain model. In the approach implemented in QuLAF, the aerodynamic loads considering the motion of the blades are simplified and approximated by loads considering a fixed hub with rigid blades and linear damping terms. The time series of fixed-hub loads and the aerodynamic damping coefficients are extracted from the SoA model for each mean wind speed.

The aerodynamic loads are obtained at each wind speed W by a SoA simulation with turbulent wind and no waves, where all DoFs except shaft rotation and blade pitch are disabled and where the wind turbine controller is enabled. The time series of fixed-hub, pure aerodynamic loads F_aero,1(t), F_aero,3(t) and τ_aero(t) are extracted from the results and stored in a data file which is loaded into the model. Hence, these SoA simulations need to be as long as the maximum simulation time needed in the simplified model (5400 s in this case).

For a given rotor, the work carried out by the aerodynamic damping is a function of wind speed, rotational speed, turbulence intensity, motion frequency and oscillation amplitude. Here, we define an equivalent linear damping which delivers the same work over one oscillation cycle and can be extracted from a decay test. Schløer et al. (2018) used this principle for the tower fore-aft mode of a bottom-fixed offshore turbine and found that the damping was only slightly dependent of the motion amplitude. We make a further simplification and carry out the decay tests in steady wind. Since the mass and stiffness of floater and tower only affect the aerodynamic damping through the motion frequency, we transfer the damping coefficients b from the decay tests in FAST to the QuLAF model. On the contrary, if the damping ratio ζ was transferred, changes in mass or stiffness properties would imply a change in the aerodynamic forcing, which is not physically correct. With the transfer of damping coefficients b, recalculation of the decay tests is only necessary in the event that the change of natural frequencies should affect the damping values significantly. Here, the decay tests from which aerodynamic damping ratios were extracted were carried out at representative natural frequencies equal to those of the present floater. These decay tests in calm water and with the wind turbine controller active were carried out for each DoF with all the other DoFs locked. This way, the floating wind turbine was a one-DoF spring-mass-damper system in each case, where the horizontal position of the hub x_hub was of interest. The decay tests were carried out as a step test in steady wind where the wind speed goes from the minimum to the maximum value with step changes every 600 s. With every step change of wind speed, the structure moves to a new equilibrium position. If all sources of hydrodynamic and structural damping are disabled, the aerodynamic damping is the only one responsible for the decay of the hub motion, and it can be extracted from the time series of x_hub. The n peaks extracted from the signal are used in pairs to estimate each local logarithmic decrement d_i and, from it, a local damping ratio ζ_i, which is then averaged to obtain the aerodynamic damping ratio ζ_aero for the given DoF and W:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{d}_i=\log {\displaystyle \frac{x_{\mathrm{hub},i}}{x_{\mathrm{hub},i+\mathrm{1}}}},\Rightarrow {\mathit{\zeta }}_i={\displaystyle \frac{d_i}{\sqrt{\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}+d_i^{\mathrm{2}}}}},\\ \text{(33)}& {\displaystyle }& {\displaystyle }\Rightarrow {\mathit{\zeta }}_{\mathrm{aero}}={\displaystyle \frac{\mathrm{1}}{n-\mathrm{1}}}\sum _{\mathrm{1}}^{n-\mathrm{1}}{\mathit{\zeta }}_i.\end{array}$$

Figure 5 Aerodynamic damping ratios for different DoFs as a function of wind speed.

Figure 4 shows examples of x_hub(t) and selected peaks for surge, pitch and clamped tower DoFs for a wind speed of 13 m s⁻¹. The wind changed from 12 to 13 m s⁻¹ at t=0, and the mean of the signals has been subtracted. For surge and pitch, peaks within the first 40 s are neglected to allow the unsteady aerodynamic effects to disappear. For the tower DoF, however, the frequency is much higher and the signal has died out by the time the aerodynamics are steady. For that reason, the tower decay peaks are extracted after 300 s, and a sudden impulse in wind speed is introduced at t=300 s to excite the tower. This method was chosen since the standard version of FAST does not allow an instantaneous force to be applied.

In Fig. 5 the aerodynamic damping ratio is shown for all DoFs as a function of W. It is observed that the aerodynamic damping in surge is negative for wind speeds between 11.4 and 16 m s⁻¹, due to the wind turbine controller. However, in real environmental conditions with wind and waves, it has been observed that the hydrodynamic damping contributes to a positive global damping of the surge motion. This controller effect is similar to the “negative damping problem” reported in, for example, Larsen and Hanson (2007). The negative aerodynamic damping in surge may be eliminated if one tunes the controller natural frequency so it lies sufficiently below the surge natural frequency of the floating wind turbine, as it was done in Larsen and Hanson (2007) for the floater pitch motion. This solution, however, would make the controller too slow and would affect power production, thus it was not adopted here because the global damping in surge has been observed to be positive when all other damping contributions are taken into account.

The damping ratio for the ith DoF and each wind speed, ζ_aero,i(W), is next converted to a damping coefficient by

$$\begin{array}{}\text{(34)}& b_{\mathrm{aero},i}\left(W\right)=\mathrm{2}{\mathit{\zeta }}_{\mathrm{aero},i}\left(W\right)\sqrt{C_{ii}(M_{ii}+A_{ii}(\mathit{\omega }\left)\right)},\end{array}$$

where C_ii, M_ii and A_ii(ω) are taken from the one-DoF oscillator in the corresponding decay test. The table of aerodynamic damping coefficients as a function of wind speed b_aero(W) is stored in a data file, which is loaded into the model. Since the aerodynamic damping coefficients are extracted from simulations with steady wind, but applied in the model in simulations with turbulent wind, an averaging is applied to account for the variability in the wind speed in turbulent conditions. Given the time series of wind speed at hub height V(t), the probability density function (PDF) of a normal distribution given by $\mathcal{N}(\stackrel{\mathrm{‾}}{V},{\mathit{\sigma }}_V)$ is used to estimate the probability of occurrence within V(t) of each discrete value of W. Then the aerodynamic coefficient for the given turbulent wind conditions and the ith DoF is

$$\begin{array}{}\text{(35)}& b_{\mathrm{aero},i}=\sum _{j=\mathrm{1}}^{N_{\mathrm{W}}}\mathrm{PDF}\left(W_j\right)b_{aero,i}\left(W_j\right).\end{array}$$

Figure 6 Surge mooring stiffness as a function of wind speed.

5.8.3 Mooring loads

The equations that provide the loads on a catenary cable depend nonlinearly on the fairlead position. In dynamic mooring models the drag forces on the mooring cables are also included; therefore, the mooring loads also depend on the square of the relative velocity seen by the lines. These nonlinear effects can easily be captured by time-domain models, but cannot be directly accommodated in a linear frequency-domain model. In QuLAF, the mooring system is represented by a linearized stiffness matrix for each wind speed, which is extracted from the SoA model and where hydrodynamic loads on the mooring lines are neglected. The dependence of the mooring matrix on wind speed is necessary because different mean wind speeds generally produce different mean thrust forces, which displace the floating wind turbine to different equilibrium states. The stiffness of the mooring system is different at each equilibrium position because of the nonlinear force-displacement behaviour of the catenary mooring lines.

For each wind speed a first SoA simulation is needed with steady, uniform wind and no waves, where only the tower fore-aft and floater surge, heave and pitch DoFs are enabled. After some time the floating wind turbine settles at its equilibrium position $({\mathit{\xi }}_{\mathrm{eq},\mathrm{1}},{\mathit{\xi }}_{\mathrm{eq},\mathrm{3}},{\mathit{\xi }}_{\mathrm{eq},\mathrm{5}})$, which is stored. These simulations should be just long enough so that the equilibrium state is reached (600 s in this case). Then, a new short SoA simulation with all DoFs disabled is run, where the floater initial position is the equilibrium with a small positive perturbation in surge, $({\mathit{\xi }}_{\mathrm{eq},\mathrm{1}}+\mathrm{\Delta }{\mathit{\xi }}_{\mathrm{1}},{\mathit{\xi }}_{\mathrm{eq},\mathrm{3}},{\mathit{\xi }}_{\mathrm{eq},\mathrm{5}})$. This simulation should be just long enough for the mooring lines to settle at rest (120 s in this case). The global mooring forces in surge and heave and the global mooring moment in pitch are stored, namely $(F_{\mathrm{moor},\mathrm{1}}^{\mathit{\xi }\mathrm{1}+},F_{\mathrm{moor},\mathrm{3}}^{\mathit{\xi }\mathrm{1}+},{\mathit{\tau }}_{\mathrm{moor},\mathrm{5}}^{\mathit{\xi }\mathrm{1}+})$. The process is repeated now with a negative perturbation in surge $({\mathit{\xi }}_{\mathrm{eq},\mathrm{1}}-\mathrm{\Delta }{\mathit{\xi }}_{\mathrm{1}},{\mathit{\xi }}_{\mathrm{eq},\mathrm{3}},{\mathit{\xi }}_{\mathrm{eq},\mathrm{5}})$, giving $(F_{\mathrm{moor},\mathrm{1}}^{\mathit{\xi }\mathrm{1}-},F_{\mathrm{moor},\mathrm{3}}^{\mathit{\xi }\mathrm{1}-},{\mathit{\tau }}_{\mathrm{moor},\mathrm{5}}^{\mathit{\xi }\mathrm{1}-})$. All this information is enough to compute the first column of the mooring matrix C_moor for the wind speed W. Perturbations in heave ±Δξ₃ and pitch ±Δξ₅ provide the necessary information to compute the rest of the columns, and therefore the full matrix:

$$\begin{array}{ll}\text{(36)}& {\displaystyle }& {\displaystyle }{\mathbf{C}}_{\mathrm{moor}}\left(W\right)={\displaystyle }& {\displaystyle }-\left[\begin{array}{cccc}{\displaystyle \frac{F_{\mathrm{moor},\mathrm{1}}^{\mathit{\xi }\mathrm{1}+}-F_{\mathrm{moor},\mathrm{1}}^{\mathit{\xi }\mathrm{1}-}}{\mathrm{2}\mathrm{\Delta }{\mathit{\xi }}_{\mathrm{1}}}}& {\displaystyle \frac{F_{\mathrm{moor},\mathrm{1}}^{\mathit{\xi }\mathrm{3}+}-F_{\mathrm{moor},\mathrm{1}}^{\mathit{\xi }\mathrm{3}-}}{\mathrm{2}\mathrm{\Delta }{\mathit{\xi }}_{\mathrm{3}}}}& {\displaystyle \frac{F_{\mathrm{moor},\mathrm{1}}^{\mathit{\xi }\mathrm{5}+}-F_{\mathrm{moor},\mathrm{1}}^{\mathit{\xi }\mathrm{5}-}}{\mathrm{2}\mathrm{\Delta }{\mathit{\xi }}_{\mathrm{5}}}}& \mathrm{0}\\ {\displaystyle \frac{F_{\mathrm{moor},\mathrm{3}}^{\mathit{\xi }\mathrm{1}+}-F_{\mathrm{moor},\mathrm{3}}^{\mathit{\xi }\mathrm{1}-}}{\mathrm{2}\mathrm{\Delta }{\mathit{\xi }}_{\mathrm{1}}}}& {\displaystyle \frac{F_{\mathrm{moor},\mathrm{3}}^{\mathit{\xi }\mathrm{3}+}-F_{\mathrm{moor},\mathrm{3}}^{\mathit{\xi }\mathrm{3}-}}{\mathrm{2}\mathrm{\Delta }{\mathit{\xi }}_{\mathrm{3}}}}& {\displaystyle \frac{F_{\mathrm{moor},\mathrm{3}}^{\mathit{\xi }\mathrm{5}+}-F_{\mathrm{moor},\mathrm{3}}^{\mathit{\xi }\mathrm{5}-}}{\mathrm{2}\mathrm{\Delta }{\mathit{\xi }}_{\mathrm{5}}}}& \mathrm{0}\\ {\displaystyle \frac{{\mathit{\tau }}_{\mathrm{moor},\mathrm{5}}^{\mathit{\xi }\mathrm{1}+}-{\mathit{\tau }}_{\mathrm{moor},\mathrm{5}}^{\mathit{\xi }\mathrm{1}-}}{\mathrm{2}\mathrm{\Delta }{\mathit{\xi }}_{\mathrm{1}}}}& {\displaystyle \frac{{\mathit{\tau }}_{\mathrm{moor},\mathrm{5}}^{\mathit{\xi }\mathrm{3}+}-{\mathit{\tau }}_{\mathrm{moor},\mathrm{5}}^{\mathit{\xi }\mathrm{3}-}}{\mathrm{2}\mathrm{\Delta }{\mathit{\xi }}_{\mathrm{3}}}}& {\displaystyle \frac{{\mathit{\tau }}_{\mathrm{moor},\mathrm{5}}^{\mathit{\xi }\mathrm{5}+}-{\mathit{\tau }}_{\mathrm{moor},\mathrm{5}}^{\mathit{\xi }\mathrm{5}-}}{\mathrm{2}\mathrm{\Delta }{\mathit{\xi }}_{\mathrm{5}}}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]\end{array}$$

The first element of the mooring matrix C_moor,11 is shown as a function of wind speed in Fig. 6. It is observed that the stiffness in surge reaches its maximum around rated wind speed (11.4 m s⁻¹), where the thrust is also maximum and the floating wind turbine is the furthest from the equilibrium position with no wind.

In the method applied here the linearization of the mooring system has been done with the SoA model. However, in a real design study where the mooring characteristics change, the above procedure can be made significantly faster by direct static analysis of the nonlinear mooring reactions around the floater equilibrium positions.

5.9 Estimation of extreme responses: a spectral approach

Classical Monte Carlo analysis of response to stochastic loads entails running a simulation, extracting the peaks from the response time series, sorting them in ascending order and assigning an exceedance probability to each peak based on their position in the sorted list. Several simulations of the same environmental conditions with different random seeds provide a family of curves in the exceedance probability plot, which can be used to estimate the expected response level for a given exceedance probability. We note that the extracted exceedance probability curves are based on the assumption that the peaks are independent, which may not always be the case. Yet, in this section the linear nature of the simplified model will be further exploited to obtain an estimation of the extreme responses to wave loads by solely using the wave spectrum and the system transfer function, thus eliminating the need for a response time series and the bias introduced by a particular random seed. An extension of the method to wind and wave forcing is further presented and discussed.

In a Gaussian, narrow-banded process, the peaks follow a Rayleigh distribution. In linear stochastic sea states, the free-surface elevation η(t) is a Gaussian random variable R_η with zero mean. Thus, within the narrow-banded assumption, which often applies to good approximation, the crest heights follow a Rayleigh distribution (Longuet-Higgins, 1956) given by

$$\begin{array}{}\text{(37)}& P(R_{\mathit{\eta }}>\mathit{\eta })=e^{-\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{\mathit{\eta }}{{\mathit{\sigma }}_{\mathit{\eta }}}\right)}^{\mathrm{2}}},\end{array}$$

where the variance in η(t) is ${\mathit{\sigma }}_{\mathit{\eta }}^{\mathrm{2}}$, which can be obtained from the integral of the wave spectrum,

$$\begin{array}{}\text{(38)}& {\mathit{\sigma }}_{\mathit{\eta }}^{\mathrm{2}}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{S}_{\mathit{\eta }}\left(\mathit{\omega }\right)d\mathit{\omega }.\end{array}$$

If we only consider linear wave forcing, the response is also Gaussian for the linear system in Eq. (3). If the response is also narrow-banded, its exceedance probability can be found via the standard deviation of the response, which in turn can be obtained by integration of the response spectrum. From Eq. (3) we have

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\widehat{\mathit{\xi }}\left(\mathit{\omega }\right)=\mathbf{H}\left(\mathit{\omega }\right)\widehat{\mathit{X}}\left(\mathit{\omega }\right)\widehat{\mathit{\eta }}\left(\mathit{\omega }\right),\\ \text{(39)}& {\displaystyle }& {\displaystyle }\Rightarrow {\widehat{\mathit{\xi }}}_{\mathrm{glob}}\left(\mathit{\omega }\right)={\mathbf{T}}_{\mathrm{glob}}\mathbf{H}\left(\mathit{\omega }\right)\widehat{\mathit{X}}\left(\mathit{\omega }\right)\widehat{\mathit{\eta }}\left(\mathit{\omega }\right)\equiv {\mathbf{TF}}_{\mathit{\eta }\to \mathit{\xi }}\left(\mathit{\omega }\right)\widehat{\mathit{\eta }}\left(\mathit{\omega }\right),\end{array}$$

where TF_η→ξ(ω) is a direct transfer function from surface elevation to global response. The global response spectra S_ξ,glob(ω) is related to the wave spectrum S_η(ω) in a similar way (Naess and Moan, 2013),

$$\begin{array}{}\text{(40)}& {\mathbf{S}}_{\mathit{\xi },\mathrm{glob}}\left(\mathit{\omega }\right)={\mathbf{TF}}_{\mathit{\eta }\to \mathit{\xi }}\left(\mathit{\omega }\right)S_{\mathit{\eta }}\left(\mathit{\omega }\right){\mathbf{TF}}_{\mathit{\eta }\to \mathit{\xi }}^{*T}\left(\mathit{\omega }\right).\end{array}$$

Here ^*T indicates the transpose and complex conjugate. By virtue of Eq. (37), the exceedance probability of, for example, the surge response ξ₁ is known from the variance in the surge response ${\mathit{\sigma }}_{\mathit{\xi },\mathrm{1}}^{\mathrm{2}}$, which is given by

$$\begin{array}{}\text{(41)}& {\mathit{\sigma }}_{\mathit{\xi },\mathrm{1}}^{\mathrm{2}}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{S}_{\mathit{\xi },\mathrm{glob},\mathrm{11}}\left(\mathit{\omega }\right)d\mathit{\omega }.\end{array}$$

For nacelle acceleration we can write the response as a function of the global nacelle displacement ξ_glob,4; therefore,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\widehat{\ddot{\mathit{\xi }}}}_{\mathrm{glob},\mathrm{4}}\left(\mathit{\omega }\right)=-{\mathit{\omega }}^{\mathrm{2}}{\widehat{\mathit{\xi }}}_{\mathrm{glob},\mathrm{4}}\left(\mathit{\omega }\right),\\ \text{(42)}& {\displaystyle }& {\displaystyle }\Rightarrow {\mathit{\sigma }}_{\ddot{\mathit{\xi }},\mathrm{4}}^{\mathrm{2}}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\mathit{\omega }}^{\mathrm{4}}{S}_{\mathit{\xi },\mathrm{glob},\mathrm{44}}\left(\mathit{\omega }\right)d\mathit{\omega }.\end{array}$$

The turbulent part of the wind speed can also be considered a Gaussian random variable (Longuet-Higgins, 1956). On the other hand, aerodynamic loads are not a linear function of wind speed. Therefore the response to wind loads cannot be assumed to be Gaussian, and the approach shown above is not valid. However, the method above can be applied to cases with wind and wave forcing, bearing in mind that the results may not be accurate since the necessary assumptions are not fulfilled. If wind and wave forcing are considered, Eq. (3) can be written as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\widehat{\mathit{\xi }}\left(\mathit{\omega }\right)=\mathbf{H}\left(\mathit{\omega }\right)\widehat{\mathit{F}}\left(\mathit{\omega }\right),\\ \text{(43)}& {\displaystyle }& {\displaystyle }\Rightarrow {\widehat{\mathit{\xi }}}_{\mathrm{glob}}\left(\mathit{\omega }\right)={\mathbf{T}}_{\mathrm{glob}}\mathbf{H}\left(\mathit{\omega }\right)\widehat{\mathit{F}}\left(\mathit{\omega }\right)\equiv {\mathbf{TF}}_{F\to \mathit{\xi }}\left(\mathit{\omega }\right)\widehat{\mathit{F}}\left(\mathit{\omega }\right),\end{array}$$

where TF_F→ξ(ω) is a direct transfer function from load to global response. The global response spectra S_ξ,glob(ω) is now given by (Naess and Moan, 2013)

$$\begin{array}{}\text{(44)}& {\mathbf{S}}_{\mathit{\xi },glob}\left(\mathit{\omega }\right)={\mathbf{TF}}_{F\to \mathit{\xi }}\left(\mathit{\omega }\right){\mathbf{S}}_F\left(\mathit{\omega }\right){\mathbf{TF}}_{F\to \mathit{\xi }}^{*T}\left(\mathit{\omega }\right).\end{array}$$

Here S_F(ω) is the spectra of the total loads (hydrodynamic and aerodynamic),

$$\begin{array}{}\text{(45)}& {\mathbf{S}}_F\left(\mathit{\omega }\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}d\mathit{\omega }}}\widehat{\mathit{F}}\left(\mathit{\omega }\right){\widehat{\mathit{F}}}^{*T}\left(\mathit{\omega }\right).\end{array}$$

This method provides the exceedance probability of the dynamic part of the response, therefore the static part should be added after applying Eq. (37). Exceedance probability results from this method are compared in the next section to the traditional way of peak extraction from response time series.

5.10 Integration of QuLAF in optimization loops

The main purpose of QuLAF is to provide a quick assessment of loads, response and natural frequencies early in the design phase, where several variations of the baseline design are to be evaluated. The efficiency in the model is achieved by (i) considering only a few DoFs; (ii) solving the linear EoMs in the frequency domain; and (iii) pre-computing the aerodynamic loads and aerodynamic damping coefficients. The application of the model to an optimization loop can be divided into two stages: a preparation stage, which needs to be done only once for a given baseline floating wind turbine, and a calculation stage, which can be repeated for each variation in the baseline design. After the optimal design has been found through optimization, it should be verified by running a complete load basis in a SoA model.

• Preparation stage. Once the wind turbine, the baseline floater design and the design basis are defined, the preparation stage entails the following.

  1. Computation of time series of aerodynamic loads at the shaft for the needed wind speeds and turbulence random seeds, considering rigid blades and fixed nacelle. The wind turbine controller should be active and tuned according to the pitch frequency of the baseline design.

  2. Extraction of aerodynamic damping coefficients for the needed wind speeds, by carrying out decay tests in steady wind of the surge, pitch and clamped tower DoFs.

  3. Storage of the aerodynamic loads and damping coefficients in a database that can be reused for several candidate designs.

• Calculation stage. The calculation stage is done for each candidate design in the pre-design optimization loop by following these steps.

  1. Computation of the radiation-diffraction solution in, for example, WAMIT.

  2. Extraction of structural mass and stiffness properties.

  3. For each wind speed, calculation of the equilibrium position and linearization of the mooring system around it.

  4. Prediction of the natural frequencies, response and loads for several environmental conditions using QuLAF.

When compared to the same number of simulations in the SoA model, the advantage of the simplified model resides in the low computational cost of applying the calculation step 4 to several environmental conditions and different variations in the baseline design. The extra work needed to achieve the speed up comes from the aerodynamic pre-computations in the preparation stage and from the linearization of the mooring system (step 3) in each iteration of the calculation stage. However, the aerodynamic loads need to be extracted from the SoA model only once for a given wind turbine, while the aerodynamic damping coefficients can also be reused for different variations of the baseline floating wind turbine, provided that the system natural frequencies do not change significantly between different design iterations. The linearization of the mooring system and the computation of the radiation-diffraction solution may also be automated. Alternatively, for slender, simpler geometries (such as spars), a Morison approach may be implemented in QuLAF, thus eliminating step 1 in the calculation stage. For the present study, however, the radiation-diffraction solution was chosen due to the shape and size of the floating substructure in consideration, and a comparison to the Morison-based alternative has not been conducted.

Table 3 Natural frequencies and periods obtained in FAST and QuLAF.

6 Validation of the QuLAF model

We now compare and discuss the QuLAF and FAST responses to the same environmental conditions (see Table 4) representative of the Gulf of Maine (Krieger et al., 2015). The cases considered include five irregular sea states with and without turbulent wind, with a single realization for each sea state. In all cases the total simulated time was 5400 s in both models. The first 1800 s were neglected to discard initial transient effects in the time-domain model. The free-surface elevation of irregular sea states was computed in FAST from a Pierson–Moskowitz spectrum, and the turbulent wind fields in TurbSim from an IEC Kaimal spectrum. Since the turbulent wind fields used in the SoA simulations are the same employed for the pre-computation of aerodynamic loads, and the free-surface elevation signal in the cascaded model is also taken from the FAST simulation, a deterministic comparison of time series is possible for all cases. In the plots shown in this section (Figs. 7 and 9), the left-hand side shows a portion of the time series of wind speed at hub height, free-surface elevation, floater surge, heave and pitch, and nacelle acceleration; and the right-hand side shows the PSDs of the same signals. The PSD signals were smoothened with a moving-average filter of 20 points to ease the spectral comparison between models. The short blue vertical lines in the PSD plots indicate the position of the system natural frequencies predicted by the simplified model (see Table 3). In addition, exceedance probability plots of the responses with both models are shown (Figs. 8 and 10), based on peaks extracted from the time series. The peaks were sorted and assigned an exceedance probability based on their position in the sorted list. The exceedance probability of the extracted peaks is compared to the one estimated with the method described in Sect. 5.9, labelled as “Rayleigh”.

Figure 7 Response to irregular waves in the time and frequency domains.

6.1 System identification

The system natural frequencies were calculated in QuLAF by solving the eigenvalue problem in Eq. (20). In FAST, decay simulations were carried out with all DoFs active, where an initial displacement was introduced in each relevant DoF and the system was left to decay. A PSD of the relevant response revealed the natural frequency of each DoF. A comparison of natural frequencies and periods found with the two models is given in Table 3, where it is shown that all floater natural frequencies in the simplified model are within 1.3 % error compared to the SoA model. On the other hand, the tower frequency is 8.6 % below the one estimated in FAST. This difference is due to the absence of flexible blades in the simplified model, which are known to affect the coupled tower natural frequency. With rigid blades, the SoA model predicts a coupled tower natural frequency of 0.684 Hz, only 0.3 % above the tower frequency in QuLAF.

Figure 8 Exceedance probability of the response to irregular waves.

The model presented here may be calibrated against other numerical or physical models if needed, by introducing user-defined additional restoring and damping matrices. For the present study, however, no calibration against the SoA model was applied, in order to keep the model calibration free and assess its suitability for optimization loops.

6.2 Response to irregular waves

The response to irregular waves with H_s=6.14 m and T_p=12.5 s (case “Waves 5” in Table 4) is shown in Fig. 7. On the frequency side, all motions show response mainly at the wave frequency range, and there is a very good agreement between both models for surge and heave. In pitch – and consequently in nacelle acceleration – the QuLAF model shows a lower level of excitation at the wave frequency range when compared to FAST. This deviation was traced to the absence of viscous forcing in the simplified model, since the two pitch responses are almost identical if viscous effects are disabled in both models. As expected, the agreement is better for milder sea states, where viscous forcing is less important. In surge and pitch some energy is visible at the natural frequencies, but only in the FAST model. Since the peaks lie out of the wave spectrum and are not captured by QuLAF, they could originate from nonlinear mooring effects or from the drag loads, which are also nonlinear.

Figure 9 Response to irregular waves and turbulent wind in the time and frequency domains.

Figure 8 shows exceedance probability plots of the response to irregular waves. The Rayleigh curves fit well to the responses given by the simplified model, which is expected, given that the free-surface elevation and the hydrodynamic forcing are linear in the model, and the response can be considered narrow-banded. In the comparison between the two models, the surge and heave peaks are very well estimated by QuLAF. In nacelle acceleration and especially in pitch, however, the model underpredicts the response, with a difference of about 30 % in pitch and about 8 % in nacelle acceleration for the largest peak when compared to FAST. These observations in extreme response are consistent with the spectral results of Fig. 7 discussed above.

Figure 10 Exceedance probability of the response to irregular waves and turbulent wind.

6.3 Response to irregular waves and turbulent wind

The response to irregular waves with H_s=6.14 m and T_p=12.5 s and turbulent wind at W=22 m s⁻¹ (case “Waves + wind 5” in Table 4) is shown in Fig. 9. The surge motion is dominated by the surge natural frequency, which is clearly excited by the wind forcing. The linear model slightly underpredicts this resonance of the wind forcing with the surge natural frequency. Heave is dominated by the wave forcing, and the response of both models agree well. In pitch, resonance with the natural frequency also exists in both models, although QuLAF predicts more energy at that frequency than FAST. Both surge and pitch responses are resonant, thus they are especially sensitive to the amount of damping. The overprediction of pitch motion also leaves a footprint on the PSD of nacelle acceleration, which shows energy at the pitch natural frequency, the wave frequency range and the tower natural frequency. The level of excitation of the tower mode at 0.682 Hz, however, is slightly underpredicted by QuLAF, likely due to an overestimation of the aerodynamic damping on the tower DoF.

The associated exceedance probability plots are shown in Fig. 10. In this case the Rayleigh curves generally do not fit the responses predicted by the linear model, as the extreme peaks are no longer Rayleigh-distributed. This is because the nonlinear nature of the wind loads makes the response non-Gaussian, and in some cases broad-banded with distinct frequency bands excited (e.g. the tower response cannot be considered narrow-banded here). The best fit is seen for heave, which is mainly excited by linear wave loads and is also narrow-banded. When compared to FAST, however, QuLAF shows a good agreement with errors in the largest response peaks of approximately 8 % in surge, 12 % in pitch and 4 % in nacelle acceleration.

6.4 Comparison of fatigue damage-equivalent loads

Table 4 shows a summary of fatigue DELs for a wider range of environmental conditions. Each case is defined by the significant wave height H_s, the wave peak period T_p and the mean wind speed W. The fatigue damage-equivalent bending moment at the tower base estimated with the two models is presented, as well as the error for the simplified model. Finally, the last column shows the ratio between the simulated time and the CPU time in QuLAF, T_rel. The cases labelled as “5” correspond to the results discussed in the previous section. The two DEL columns in Table 4 are also shown in Fig. 11 as a bar plot.

Table 4 Summary of environmental conditions (Krieger et al., 2015) and DEL results obtained in FAST and QuLAF.

For the cases with waves only, the model underpredicts the DEL at the tower base with errors from 0.2 % to 11.3 % that increase with the sea state, as observed in Fig. 11. The significant wave height also increases with the sea state, as do the associated nonlinear effects of position-dependent mooring stiffness and viscous hydrodynamic forcing, which are both included in FAST. QuLAF does not include viscous hydrodynamic forcing, and as a linear model, its accuracy is bound to the assumptions of small displacements around the equilibrium point. Hence, it is expected that the linear model performs worse for the environmental conditions where nonlinear effects are not negligible. This observation is also consistent with the discussion around Fig. 7, which corresponds to the most severe sea state considered here.

Figure 11 Damage-equivalent bending moment at the tower base for different environmental conditions.

For the cases with wind, the errors range from 1.5 % to 6.9 %, but the trend is not as clear. The predictions seem to be worst for the environmental condition corresponding to rated wind speed. Around rated speed the wind turbine operation switches between the partial- and the full-load regions, which correspond to very distinct regimes of the generator torque and blade pitch controller. The complexity of the dynamics involved in this transition zone is not well captured by the simplified model. The vibration of the tower is also more likely to be excited around rated wind speed, where the thrust is maximum. As the coupled tower natural frequency is different for the two models, this will also have an impact on the resulting DEL. This effect has been quantified for rated wind speed (”Waves + wind 3”), where the DEL error becomes −5.6 % when the FAST simulation is carried out with rigid blades, which indicates that the difference in coupled tower frequency has some impact on the DEL error. In addition, the aerodynamic damping – which plays an important role in the resonant response of the tower – is dependent on the frequency at which the rotor moves in and out of the wind. Since the aerodynamic damping on the tower is extracted from a SoA simulation with fixed foundation and rigid blades, it corresponds to a tower natural frequency of 0.51 Hz, different to the coupled tower frequency observed when the floater DoFs are active (0.682 Hz in QuLAF, 0.746 in FAST). This difference in the frequencies at which the aerodynamic damping is extracted and applied is likely to lead to an overprediction of the aerodynamic damping, and an underprediction of the tower vibration and the DEL. This observation is consistent with the level of tower response at the coupled tower frequency shown in Fig. 9. On the other hand, the aerodynamic simplifications in the cascaded model seem to work best for wind speeds above rated, likely due to the thrust curve being flatter in this region. The last column of Table 4 shows that the ratio between simulated time and CPU time is between 1300 and 2700 for a standard laptop with an Intel Core i5-5300U processor at 2.30 GHz and 16 GB of RAM. In other words, all the simulations in Table 4 together, 1.5 h long each, can be done in about half a minute.

7 Conclusions

A model for Quick Load Analysis of Floating wind turbines, QuLAF, has been presented and validated. The model is a linear, frequency-domain tool with four planar degrees of freedom (DoFs): floater surge, heave, pitch and tower modal deflection. The model relies on higher-fidelity tools from which hydrodynamic, aerodynamic and mooring loads are extracted and cascaded. Hydrodynamic and aerodynamic loads are pre-computed in WAMIT and FAST, respectively, while the mooring system is linearized around the equilibrium position for each wind speed using MoorDyn. A simplified approach for viscous hydrodynamic damping was implemented, and the decay-based extraction of aerodynamic damping of Schløer et al. (2018) was extended to multiple DoFs. Without introducing any calibration, a case study with a semi-submersible 10 MW configuration showed that the model is able to predict the motions of the system in stochastic wind and waves with acceptable accuracy. The damage-equivalent bending moment at the tower base is estimated with errors between 0.2 % and 11.3 % for all the five load cases considered in this study, covering the operational wind speed range. The largest errors were observed for the most severe wave climates in wave-only conditions and for turbine operation around rated wind speed for combined wind and wave conditions, due to three main limitations in the model: (i) underprediction of hydrodynamic loads in severe sea states due to the omission of viscous drag forcing; (ii) difficulty to capture the complexity of aerodynamic loads around rated wind speed, where the controller switches between the partial- and full-load regions; and (iii) errors in the estimation of the tower response due to underprediction of the coupled tower natural frequency and overprediction of the aerodynamic damping on the tower. The computational speed in QuLAF is between 1300 and 2700 times faster than real time. Although not done in this study, introducing viscous hydrodynamic forcing and calibration of the damping against the SoA model would likely result in improved accuracy, but at the expense of lower CPU efficiency and less generality in the model formulation.

It has been shown that the model can be used as a tool to explore the design space in the preliminary design stages of a floating substructure for offshore wind. The model can quickly give an estimate of the main natural frequencies, response and loads for a wide range of environmental conditions with aligned wind and waves, which makes it useful for optimization loops. Although a better performance may be achieved through calibration, a calibration-free approach was used here to emulate the reality of an optimization loop, where calibration is not possible. In such a process, once an optimized design has been found, a full aero-hydro-servo-elastic model is still necessary to assess the performance in a wider range of environmental conditions, including nonlinearities, transient effects and real-time control. Since the model is directly extracted from such a SoA model, this step can readily be taken. While the SoA model should thus still be used in the design verification, the present model provides an efficient and relatively accurate complementary tool for rational engineering design of offshore wind turbine floaters. In addition, the QuLAF and FAST models presented in this study have been recently used in the LIFES50+ project for a broader analysis of different design-driving load cases, including normal operation, extreme and transient events (Madsen et al., 2018). Generally, the results of the broader study and the conclusions drawn are aligned with the ones presented here, as well as the limitations observed in the simplified model when compared to its SoA counterpart.

Given the model limitations observed in this study and in Madsen et al. (2018), possible improvements of QuLAF may involve (i) inclusion of viscous drag forcing, (ii) modelling the effect of blade flexibility on the tower natural frequency, (iii) improvement of the extraction of aerodynamic damping from the SoA model, and (iv) extension of the model to out-of-plane DoFs to make it applicable to cases with misaligned wind and waves.