A Short Review on the Time-Domain Numerical Simulations for Structural Responses in Horizontal-Axis Offshore Wind Turbines

Abstract

In addition to a carbon-neutral vision being recognized worldwide, the utilization of wind energies via horizontal-axis wind turbines, especially in offshore areas, has been intensively investigated from an academic perspective. Numerical simulations play a significant role in the design and optimization of offshore wind turbines. The current review focuses on studies concerning the numerical simulations of offshore wind turbine dynamics, including the modelling of the aerodynamic and hydrodynamic conditions of the environment and the reduced-order modelling of the wind turbine dynamic responses. In detail, the functions and mechanisms of each module in the numerical simulation of the wind turbine dynamics are articulated, which in turn demonstrates its importance for the design of offshore wind turbines, and hence the development of the offshore wind industry. Based on this review, it is argued that the vertical variations in wind velocities, the blade element momentum theory, the wave dynamic models, and the reduced-order model for structural dynamics are the major concerns for the numerical simulation of wind turbines. Consequently, such directions should be emphasized in future studies.

1. Introduction

Due to the rising anxiety about the environmental consequences and energy crisis related to economic volatility and even geopolitical instability, efforts have been continuously made to shift energy consumption from conventional sources such as fossil fuels to a more eco-friendly energy mix. More specifically, the supply of energy from renewable sources has increased significantly over the past decade due to the rapid growth in solar and wind resource acquisition [1,2]. It is widely acknowledged that renewable energy is key to the sustainable development of society, and hence it is a topic for academic studies concerning sustainability. In fact, wind energy is considered the most promising in the process of replacing conventional fossil fuel sources to sustain the healthy development of society. Given its renewable nature and close-to-zero carbon emissions, the wind energy industry is an important module for a carbon-neutral vision and for the realization of sustainability. As a crucial component of the wind energy industry, the exploitation of offshore wind sources plays a significant role as it provides an ample and steady energy supply close to the energy consumption centre. Therefore, offshore wind has gradually gained momentum recently as it has been actively promoted by many countries [3]. In fact, investment in offshore wind energy worldwide has already exceeded that of offshore oil and gas since 2021, and its total grid-connected installed capacity reached 64.3 GW by the end of 2022.

Although the exploitation of offshore wind sources is a widely acknowledged module for a carbon-neutral vision, the operational strategy and power generation control of offshore wind turbines face various challenges [4]. In the operation of wind turbines installed in a wind farm, a well-designed algorithm is required to stabilize the power generation and mitigate the structural vibration, which improves its general serviceability [5,6]. Currently, the control algorithm for a Horizontal-Axis Wind Turbine (HAWT) adjusts the torque of the generator and the pitch angle of the blade to optimize the turbine performance according to environmental wind velocities [6,7]. While the generator torque is controlled to maximize the power generation when the environmental wind velocity is in the range between the cut-in and rated velocities, the blade pitch angle is adjusted to stabilize the power generation when the environmental wind velocity exceeds the rated velocity of a wind turbine.

Numerical simulation plays a significant role in the development and implementation of novel control algorithms, as the wind turbine state is required as an input for the algorithm to devise an effective control strategy. On the one hand, the numerical simulation provides the overall information on a wind turbine operated under various wind velocities. On the other hand, it delivers the structural responses of a wind turbine, which is helpful for the algorithm to balance the power generation and the vibration mitigation.

The HAWT is a highly complex system consisting of a tower, a series of blades, a hub, a nacelle, a drivetrain, and a platform/foundation. While the tower and the blades are flexible in nature, other parts of a wind turbine are generally modelled as rigid bodies [8]. Compared to the motion simulation of rigid bodies, the flexible simulation of towers and blades is more complex. At present, the structural dynamics of the tower and blades are commonly simulated using the finite element method or a reduced-order model [9]. In fact, high-fidelity simulations, which are widely used in the design of control algorithms, typically employ finite element models of blades and towers to provide more accurate estimates of their structural responses [10]. Efforts have been made to advance the dynamic modelling of the blade because of its importance in the overall dynamics of the wind turbine. For example, finite element tools such as NuMAD [11] and FOCUS [9] are recommended to simulate the vibration and deformation of the blades. Since geometric nonlinearity should be taken into consideration [12,13], Outer Mold Layer shell elements [12,13,14] are suggested to simulate the HAWT blade. For instance, the S4R-type element in the ABAQUS library, which is essentially a four-node quadrilateral shell element, is frequently adopted to capture the blade dynamics [12,15]. Because the use of finite element methods requires relatively high computational costs, the reduced-order model is applied to discern the key dynamic characteristics of a wind turbine in order to develop and implement a novel control algorithm [16]. In addition, a reduced-order model is employed to assess the influence of different control algorithms on the structural vibration of the wind turbine.

The optimization of offshore wind turbines has become critical as the demand for renewable energy rises. Our study focuses on time-domain numerical simulations of the structural responses of HAWTs. In fact, the present study reviews advances in the numerical simulations of the HAWT dynamics, which provide indications for further development. The review points in the direction of increasing the operational lifespan of these turbines and lowering the financial and environmental costs associated with the offshore wind industry.

In the following pages, Section 2 discusses the wind profile modelling techniques, and Section 3 presents a comprehensive review of simplified models for Fluid–Structure Interaction (FSI), aimed at improving the estimation of aerodynamic loads. The hydrodynamic load simulations are reviewed in Section 4. Section 5 reviews the reduced-order model. Concluding remarks are provided in Section 6.

2. Wind Field Modelling

Since the environmental wind velocity is the key piece of information enabling the algorithm to control the wind turbine and the numerical simulation to assess the dynamic responses of a wind turbine [17], it is usually determined based on the in situ measurements and an empirical model of wind profiles [18,19]. Because the wind sensor is usually installed at a certain height and the wind velocity varies vertically in the atmospheric boundary layer, the wind profile model is necessary to convert the measured wind velocity to various heights for the calculation of wind loads acting on the blades and the tower. Although the computational fluid dynamic simulation could produce the wind profile without empirical models, the relatively high computational cost hinders its application for discerning the wind profile for developing the control strategy of a wind turbine [20]. Consequently, different empirical or semi-empirical models are commonly used to depict the vertical variation in wind velocities.

In fact, the wind flow is usually assumed to be in a steady and horizontally homogeneous state in the wind profile models [21]. In other words, the temporal and horizontal variations in the atmospheric boundary layer are averaged in the empirical wind profile models. From the boundary layer theory, the vertical profile of the averaged wind field is determined by the roughness of the underlying terrain, and hence possesses a logarithmic shape. In addition to the logarithmic model, the vertical variation in wind velocities is also depicted by the power law and presents a similar shape. In fact, the wind profile can be described as

$$U(z)=\frac{u_{\ast }}{k}\ln \left(\frac{z-z_h}{z_0}\right)$$

(1)

$$U(z)=U_{ref}{\left(\frac{z}{z_{ref}}\right)}^{\alpha }$$

(2)

In Equation (1), corresponding to the log-law profile model, $k$ is the Von Karman constant, which commonly takes the value of 0.41, $z$ is the height above the ground, $z_h$ is the zero plane displacement, $z_0$ is the roughness length of the underlying terrain, and $u_{\ast }$ is the shear velocity, which shows the drag from the underlying terrain [21]. In Equation (2), corresponding to the power-law model, $z_{ref}$ is the reference height, which usually takes the hub height of the wind turbine, and $U_{ref}$ is the wind velocity at the reference height. In the power-law model, the shear exponent of $\alpha$ is related to the roughness of the underlying terrain, and hence determines the general shape of the wind profile.

In addition to the mean wind profile, the turbulence also influences the aerodynamic loads acting on the wind turbine. Specifically, the comprehensive modelling of the atmospheric boundary layer wind field adds the turbulent fluctuations described by certain power spectral density models on the mean wind profile. Several models are available to calculate the power spectral density of the wind field inside the atmospheric boundary layer, and the Kaimal model obtained from analysing in situ measurements of wind velocities is widely accepted to present the turbulent characteristics of winds. The Kaimal model of the wind power spectral density shows the following:

$$S(f,z)=\frac{{I_u}^{2}U\left(z\right)l}{{\left(1+1.5\frac{fl}{U\left(z\right)}\right)}^{\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}}$$

(3)

In Equation (3), $S\left(f\right)$ is the power spectral density with the frequency of $f$, $U\left(z\right)$ is the mean wind velocity, $l$ represents the turbulent length scale, and $I_u$ gives the turbulence intensity in the longitudinal direction, which is defined as

$$I_u=\frac{{\sigma }_u}{U\left(z\right)}$$

(4)

In Equation (4), ${\sigma }_u$ is the turbulent wind velocity in the longitudinal direction, which is calculated as the standard deviation for a time series of longitudinal wind velocities.

Given the vertical profile and power spectral density of mean and turbulent winds, a pseudo stochastic wind field $u(z,t)$ can be generated via summing the series cosine functions with random phases, as follows:

$$u(z,t)=U(z)+\sum _{i=1}^{N/2}\sqrt{4\pi \Delta fS\left(f_i\right)}\cos (2\pi f_{i}t-{\phi }_i)$$

(5)

In Equation (5), the power spectra are discretized at the frequencies of $f_i,i=\mathrm{1,2},\cdots ,N$ with an equal frequency step of $\Delta f$, and the randomness is introduced via the phases ${\phi }_i$, varying within the range of $0-2\mathsf{\pi },$ corresponding to various frequencies [22].

Equation (5) presents the generation of pseudo stochastic wind velocities at a specific point, and the time domain simulation of wind turbine dynamics also requires the spatial variations in the wind field as the input. Therefore, spatial correlations of wind velocities should be modelled in addition to the single-point wind time series. More specifically, the coherence is introduced to show the spatial correlation of wind velocities in the frequency domain. Based on measurements accumulated in the field of meteorology, the coherence is suggested to be modelled as follows:

$$\mathrm{c}\mathrm{o}\mathrm{h}(L,f)=\mathrm{e}\mathrm{x}\mathrm{p}\left(-12\frac{fL}{U}\right)$$

(6)

In Equation (6), $\mathrm{c}\mathrm{o}\mathrm{h}$ shows the coherence as functions of the spatial distance $L$ and the frequency of $f$. With the help of the coherence defined in Equation (6), the wind field is not only determined by the power spectral density but also by the cohesive power spectral density, as follows:

$$C(\Delta x,f)=\mathrm{c}\mathrm{o}\mathrm{h}(L,f)\sqrt{S(z_1,f)S(z_2,f)}$$

(7)

In Equation (7), $\Delta x$ shows the vector difference between two points in the space at the heights of $z_1$ and $z_2$, and their spatial distance is $L$.

In practice, the space to generate the pseudo stochastic wind field is usually discretized into a mesh, and the power spectral densities and cohesive spectral densities are usually organized into a matrix corresponding to the grids of the mesh. By decomposing such a matrix to specify the magnitudes of cosine variations as shown in Equation (5), the full-set pseudo stochastic wind field can be generated for the time domain simulation of the wind turbine dynamics. It is noted that such a pseudo stochastic wind field only shows the turbulent fluctuations, and the mean wind profile should be added to drive the numerical time-domain simulation of the HAWT dynamics.

Although the spectral characteristics summarized based on the in situ measurements provide the most common base for generating the pseudo-stochastic wind field as the input for the time-domain simulation of wind turbine dynamics, other approaches have also been suggested in previous studies. For example, the spectral tensor model has been suggested within the framework of the rapid distortion theory to show the development of the natural wind field from a hypothetical, isotropic state [23]. The discretization of the spectral tensor therefore shows spatial variations in the natural wind field, which are then used to generate the pseudo-stochastic wind field given a proper spectral tensor.

In terms of numerical tools for the generation of a pseudo-stochastic wind field for the purpose of numerically simulating the wind turbine dynamics, TurbSim, developed by the National Renewable Energy Laboratory (NREL), is frequently used [24,25,26,27,28,29]. TurbSim describes the natural wind field based on the Kaimal power spectral density model and provides the user with various options to calculate the cohesive power spectral densities, including the common von Karman model and the Riso-Smooth-Train model [30].

3. Aerodynamic Modelling for Fluid–Structure Interaction (FSI)

Once the turbulent wind field has been generated, the aerodynamic loads acting on the blades and the tower can be estimated via the Blade Element Momentum (BEM) theory [31,32,33]. More specifically, the temporally varying aerodynamic loads can be estimated according to the wind velocity time series generated from the pseudo-stochastic wind field model, which then leads to the estimation of the ultimate torque acting on the generator and bending moment at the tower base [34]. In the estimation of the aerodynamic loads acting on the blades, the blade is segmented spanwise into a series of sections. The drag and lift forces acting on the sections along an individual blade are accumulated as the linear and angular momentum acting on the blade root, which consequently presents the aerodynamic torque for the estimation of the power generation and also internal forces in the blades for the structural vibration assessment [6,35].

As for the floating HAWT, the aerodynamic loads acting on the blades and the tower are not solely determined by the inflow wind velocities [6,36]. In fact, the motion of the floating foundation should be considered, and the relative velocity between the blade motion should also be taken into consideration for the floating foundation motions and the inflow wind velocity to determine the aerodynamic loads acting on the blades of the floating HAWT [37]. Figure 1 presents the motion responses of a floating HAWT, and hence illustrates the necessity of using the relative velocity to calculate the aerodynamic forces on the floating HAWT [38].

Based on the rigid assumption of the wind turbine, the relative velocities of a blade element with a distance of $r$ from the rotational center of the blades can be calculated as follows [38]:

$$u=U_{\infty }(\cos {\theta }_y-a)+\dot{z} \cos {\theta }_y+(H+r\cos \varphi )\dot{{\theta }_y}$$

(8)

$$w=r{\omega }_r(1+a’)-(\dot{z}+U_{\infty })\sin {\theta }_y$$

(9)

In Equations (8) and (9), $U_{\infty }$ is the inflow velocity perpendicular to the rotation plane of the HAWT and ${\theta }_y$ and $\dot{z}$ show the pitch and heave responses of the floating HAWT. While $H$ is the height of the rotational center and $r$ is the radial distance of the blade section from the rotational center, $\varphi$ presents the azimuth angle of a certain blade. $a$ and $a’$, on the other hand, are the axial and tangential induction factors, which are calculated in Equations (16) and (17). It is noted that the aerodynamic loads acting on the blade rely primarily on the relative velocity component perpendicular to the blade ($u$) and in the vertical direction ($w$), and hence the component along the blade ($v$) is neglected. The inflow angle, which is key to calculating the angle of attack (AOA) and then for the induction factors, is therefore determined as follows:

$$\tan \phi =\frac{w}{u}$$

(10)

The aerodynamic loads acting on the blades, specifically the lift and drag forces, are mainly determined by the AOA, the Mach number, and the Reynolds number. While the Mach number is typically small enough in the case of HAWT (usually less than 0.3) to be ignored, the Reynolds number interacts with the airfoil design of the blade to impact the aerodynamic load calculation. Nonetheless, it is commonly acknowledged that the AOA (${\alpha }_{attack}$), as shown in Figure 2, is the most influential factor in the estimation of the drag and lift forces acting on the blade [39,40], defined as the difference between the inflow angle and the $\theta$ (sum of the pitch angle ($\beta$) and the twisting angle ($\vartheta$) of the specific section), as follows:

$${\alpha }_{attack}=\phi -(\beta +\vartheta )$$

(11)

With the known AOA, the drag ($C_d$) and lift ($C_l$) coefficients can be derived from a table containing the data accumulated from the experiment or full-set numerical simulation results [41]. In fact, it is found from the literature that Reynolds numbers [42], the general geometry of the blade [43], the use of end-plates [44], and other factors all influence the aerodynamic characteristics of the blade. Hence, it is necessary to obtain the drag and lift coefficient database specifically for the blade design under investigation. Given the inflow angle of $\phi$, the normal ($C_n$) and tangential ($C_t$) load coefficients are calculated as follows:

$$C_n=C_l\cos \phi +C_d\sin \phi$$

(12)

$$C_t=C_l\cos \phi -C_d\sin \phi$$

(13)

In addition, the thrust and torque acting on the rotating centre induced by a single blade are calculated as follows:

$$F=\frac{1}{2}\rho {\int }_{R_h}^{R_b}{C}_n(u^2+w^2)c\left(r\right)dr$$

(14)

$$M=\frac{1}{2}\rho {\int }_{R_h}^{R_b}{C}_t(u^2+w^2)rc\left(r\right)dr$$

(15)

In Equations (14) and (15), $F$ is the thrust, $M$ is the torque, $c\left(r\right)$ shows the chord length as a function of the radial distance ($r$), and $R_h$ and $R_b$ are the radius of the hub and the radial distance of the blade tip, respectively. The axial and tangential induction factors are functions of the relative velocity, solidity, normal, and tangential load coefficients, as follows:

$$a=\frac{1}{1+\frac{4{\sin }^2\phi }{\sigma C_n}}$$

(16)

$$a’=\frac{1}{1+\frac{4\sin \phi \cos \phi }{\sigma C_t}-1}$$

(17)

In Equations (16) and (17), $\sigma$ is the solidity of the blade section and is defined as the ratio of the swept area and the corresponding control volume.

It is understandable that the aerodynamic loads estimated according to classic BEM theory are inaccurate in a number of aspects [45]. Consequently, there are several corrections reported in the literature to reduce the bias of the assumptions employed in the classic BEM theory (relevant details are available in [46]). For example, the Prandtl correction factor has been suggested to correct the unrealistic infinite blade section assumption [47], and the Glauert correction is suggested for cases where the axial induction factor is beyond 0.4 [6,34].

The calculation of the aerodynamic loads acting on the HAWT according to the BEM theory is in fact an iterative process starting with an assuming axial and tangential factor [48]. The AOA can then be calculated according to Equation (11), which leads to the estimates of lift and drag coefficients. Afterwards, the normal and tangential load coefficients are delivered to give the estimates of the induction factors in the current iteration, as shown in Equations (16) and (17). Such an iterative process continues until the residuals corresponding to both induction factors reduce below the specific thresholds. The ultimate estimates of the thrust and torque acting on the blades are then calculated according to Equations (14) and (15). Alongside this conventional approach, different algorithms have also been suggested, such as the simplification reported in [49].

4. Hydrodynamic Modelling for Offshore Wind Turbines

Along with the exploitation of wind resources moving from onshore sites to offshore sites, and further from shallow waters to deep-sea areas, the numerical simulation of floating HAWT has attracted attention from both academia and industry [6]. For a floating HAWT, the hydrodynamic loads should be realistically and accurately modelled for the numerical simulation of its overall dynamic responses [5]. Therefore, many methods have been suggested to estimate the hydrodynamic loads acting on various components of the floating HAWT [50].

Among the various factors inducing the hydrodynamic loads, waves play the most significant role as they are associated with the most destructive effect [51]. For estimating the wave loads acting on the floating HAWT, the velocity potential theory, which is used to describe a flow without viscosity and rotation, is widely adopted to show the kinetics, and hence dynamics, of waves and their loading effect [7]. Specifically, a potential function ($\phi$) is introduced which summarizes the flow velocity components ($V$) in three orthogonal directions via spatial gradients as $V=\nabla \phi$ [52,53]. Given the velocity potential function, the continuity equation of the flow is reduced to the Laplace equation of ${\nabla }^2\phi =0$. Solving the Laplace equation with specific boundary conditions yields the field of flow velocities, which can then be used to estimate hydrodynamic loads following the Bernoulli equation. Although the introduction of the velocity potential simplifies the calculation of hydrodynamic loads, the track of free surfaces ($\xi (x,y,t)$), which is key in the specification of the upper boundary condition of the Laplace equation, is challenging [51,54,55]. The linearization of the boundary condition at the free surface is therefore necessary to derive the analytical solutions describing the kinematics of waves. It is noted that the combination of the Laplace equation and the linearized boundary condition leads to the description of regular waves, but the realistic waves in the sea are irregular in nature [56]. In order to depict the natural wave field realistically, wave spectra are suggested based on long-term observations obtained worldwide. Among various wave spectra models, the JONSWAP and P-M models are most commonly used in the field of ocean engineering, and hence are frequently adopted for calculating the hydrodynamic loads acting on the floating HAWT. The JONSWAP spectra model shows the following:

$$S(f)=0.3125H_s{T}_p{\left(\frac{f}{f_p}\right)}^5\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{5}{4}{\left(\frac{f}{f_p}\right)}^{-4}\right](1-0.287\mathrm{l}\mathrm{n}\gamma ){\gamma }^{\mathrm{e}\mathrm{x}\mathrm{p}\left[\frac{{(\omega -{\omega }_p)}^2}{2{\sigma }^2{{\omega }_p}^2}\right]}$$

(18)

In Equation (18), $T_p$ and $f_p$ are the peak wave period and frequency, $H_s$ is the significant wave height, and both are estimated based on long-term in situ observations. $\sigma =0.07$ in the cases where $f<f_p$, and $\sigma =0.07$ corresponding to $f\ge f_p$. $\gamma$ is the peak parameter of the JONSWAP spectra model, and can be calculated as follows:

$$\gamma =\left\{\begin{array}{cc}5& \frac{T_p}{\sqrt{H_s}}\le 3.6\\ \mathrm{e}\mathrm{x}\mathrm{p}\left(5.75-1.15\frac{T_p}{\sqrt{H_s}}\right)& 3.6<\frac{T_p}{\sqrt{H_s}}\le 5\\ 1& \frac{T_p}{\sqrt{H_s}}>5\end{array}\right.$$

(19)

Given the wave spectra model, the elevations and then velocities of the wave field can be numerically simulated via superimposing the cosine functions with different frequencies and random phases, as follows:

$$\eta (t)={\sum }_{j=1}^N{A}_j\sin \left({\omega }_{j}t-k_{j}x+{\varphi }_j\right)$$

(20)

$$A_j=\sqrt{2S\left({\omega }_j\right)\Delta \omega }$$

(21)

$$u={\sum }_{j=1}^N{\omega }_j{A}_j\frac{\cosh \left[k(z_w+d_w)\right]}{T_p\sinh \left(k{d}_w\right)}\sin \left({\omega }_{j}t-k_{j}x+{\varphi }_j\right)$$

(22)

In Equations (20)–(22), $\eta \left(t\right)$ is the wave elevation time series, ${\omega }_j$ is the angular frequency related to the frequency ($f$) as ${\omega }_j=2\pi f_j$, $A_j$ is the amplitude of a regular wave component at a certain frequency, and ${\varphi }_j$ is the phase angle, taking the value in the range of $[0-2\mathsf{\pi }]$, which is usually modelled as a random variable with uniform distribution. $d_w$ is the water depth and $z_w$ shows the vertical distance from the still water surface. It is noted that the wave number $k$ is related to the angular wave frequency $\omega$ through the dissipation equation as follows:

$$k\tanh k{d}_w=\frac{{\omega }^2}{g}$$

(23)

With the kinematics either solved analytically from the velocity potential theory or simulated numerically from the empirical wave spectra model, the hydrodynamic loads are commonly calculated from the Morrison equation shown in Equation (24).

$$F=C_M\rho V_0\frac{du}{dt}+\frac{1}{2}{C}_d\rho A_0|u|u$$

(24)

In Equation (24), $u$ is the velocity of water particles in a wave field, and is obtained from the wave kinematics; $C_M$ and $C_d$ are empirical coefficients generally extracted from in situ observations or scaled-down laboratory experiments; and $\rho$, $V_0$ and $A_0$ show the density, the characteristic volume, and the area of the floater under investigation, respectively. It is noted that the Morrison equation is only applicable to estimate the wave forces experienced by slender structures whose horizontal dimension is generally less than 0.15 times the wavelength corresponding to the dominant wave frequencies [54]. When the dimension of the offshore structure, or its key component, exceeds the above-mentioned threshold, the Morrison equation is no longer applicable [57]. In such cases, solving the Laplace equation with complex solid boundary conditions corresponding to the offshore structure is unavoidable. In practice, the boundary element method is commonly used to deliver the velocity, and hence the dynamic pressure, around the object of interest by numerically solving the Laplace equation. Within the framework of the boundary element method, a branch of software (such as WAMIT Ver. 7.0 [58]) has been developed for estimating the dynamic loads induced by waves acting on the offshore structures with large dimensions.

Although the wave force is the main reason for the platform to move in the system of the floating HAWT, the mooring lines subjected to both wave and current loads also influence the dynamics of the floating platform [59,60,61]. Understandably, the Morrison equation is utilized to estimate the wave load acting on the mooring lines due to the small cross-section in the comparison with the wavelength [59]. From the hydrodynamic forces calculated from the environmental conditions, the dynamic responses of the mooring system are simulated using the Kane or Euler–Lagrange methods [62]. Commonly, the mooring lines are modelled as continuous beams vibrated under the excitation of both wave and current loads, in which the bending and torsional stiffness are taken into consideration. As shown in Figure 3, the mooring system in a floating HAWT connects the floating platform at the fairlead and the seabed at the anchor point [62].

5. The Reduced Order Model

After the aerodynamic and hydrodynamic loads acting on the HAWT have been calculated, and the motion of the floating platform has been restrained by the mooring system, the structural and motion responses of the floating HAWT can be numerically simulated following the procedures sketched in Figure 4, once the structural responses of the wind turbine have been simulated via the reduced order model.

Conventionally, in the reduced order models, the blade and the tower are modelled as elastic beams in the numerical simulation of the dynamic responses of the HAWT, according to either the Euler–Bernoulli or Timoshenko theory [63,64]. While the Euler–Bernoulli beam theory is inadequate for the thick beam since it neglects the shear deformation and rotatory inertia [65], the Timoshenko beam theory describes the deformation via a set of partial differential equations with infinite degrees of freedom (DOFs), which are mathematically difficult to solve [5,54]. Hence, the Galerkin method is introduced to reduce the complex partial differential equation into a series of ordinary differential equations with a limited number of DOFs by separating the variable, describing the dynamic characteristics of the beam linearly.

Based on the Galerkin simplification of the Timoshenko beam theory [16], the blade and the tower of the HAWT are modelled as thin beams whose in-plane and out-plane vibrations can be decoupled (shown in Figure 5) as follows [16,66]:

$$u_{in,j}={\sum }_{i=1}^N{\phi }_{in,i}{q}_{in,i}, j=\mathrm{1,2},3$$

(25)

$$u_{out,j}={\sum }_{i=1}^N{\phi }_{out,i}{q}_{out,j}, j=\mathrm{1,2},3$$

(26)

In Equations (25) and (26), ${\phi }_{in,i}$ and ${\phi }_{out,i}$ are the $i^{th}$ mode shape for the vibration of the blade/tower corresponding to the in-plane and out-plane directions, respectively; $q_{in,i}$ and $q_{out,j}$ are the participation coefficients of the corresponding modes; and $u_{in,j}$, $u_{out,j}$ show the resulting vibration of the specific blade or the tower. It is noted that the mode shape ${\phi }_{in,i}$ and ${\phi }_{out,i}$ are estimated based on the geometric and structural characteristics of the blade or the tower, and usually normalized either at the blade tip or the tower top [67]. Therefore, the reliable estimates of the vibration modes depend on the precise data concerning the twist, stiffness, and density of the blade or the tower.

Taking the NREL 5 MW baseline wind turbine as an example, the reliable estimates of the vibration mode shape of both blades and the tower, which in turn lead to the correct simulation of its dynamic responses under wind and wave loads, require the general and aerodynamic properties as summarized in

Table 1. The general properties of the NREL 5 MW baseline wind turbine [68].

Rating	5 MW
Rotor orientation, Configuration	Upwind, 3 blades
Control	Variable speed, Collective pitch
Drivetrain	High speed, Multiple-stage gearbox
Rotor, Hub diameter	126 m, 3 m
Hub height	90 m
Cut-In, Rated, Cut-out wind speed	3 m, 11.4 m, 25 m
Rated tip speed	80 m/s
Overhang, Shaft tilt, Precone	$5 \mathrm{m}, 5°$$, 2.5°$
Rotor Mass	110,000 Kg
Nacelle Mass	240,000 Kg
Tower Mass	347,460 Kg
Coordinate location of overall CM	(−0.2 m, 0.0 m, 64.0 m)

and

Table 2. The blade aerodynamic properties of the NREL 5 MW baseline wind turbine.

Node	Rnode	Aero-Twist	Chord	Airfoil Type
1	2.8667	13.308	3.542	Cylinder 1
2	5.6000	13.308	3.854	Cylinder 2
3	8.3333	13.308	4.167	Cylinder 2
4	11.7500	13.308	4.557	DU40_A17
5	15.8500	11.480	4.652	DU35_A17
6	19.9500	10.162	4.458	DU35_A17
7	24.0500	9.011	4.249	DU30_A17
8	28.1500	7.795	4.007	DU25_A17
9	32.2500	6.544	3.748	DU25_A17
10	36.3500	5.361	3.502	DU21_A17
11	40.4500	4.188	3.256	DU21_A17
12	44.5500	3.125	3.010	NACA64_A17
13	48.650	2.319	2.764	NACA64_A17
14	52.7500	1.526	2.518	NACA64_A17
15	56.1667	0.863	2.313	NACA64_A17
16	58.9000	0.370	2.086	NACA64_A17
17	61.6333	0.106	1.419	NACA64_A17

, respectively [68].

The HAWT is a complex multi-body system consisting of the blades, the tower, the drivetrain, and the hub [69]. Even when the dynamic vibration of the blades and the tower is analysed with the reduced order method, the interactions within this multi-body system still put forward a heavy computational burden for the dynamic assessment of the HAWT. Consequently, the numerical simulation of the overall dynamics of the HAWT is still in need of reduced order models. Generally, the Euler–Lagrange or the Kane approaches have been employed in the simulation of the HAWT to date [16,67,70,71,72].

5.1. The Euler–Lagrange Approach

In a multi-body system, the interactions among different bodies, and their responses under external excitation, could be modelled via the generalized coordinate system following the Euler–Lagrange concept. More specifically, the dynamic system without constrain is governed by the following:

$$\frac{d}{dt}\left(\frac{\partial L}{\partial \left\{q\right\}}\right)-\frac{\partial L}{\partial \left\{q\right\}}=\{Q\}$$

(27)

In Equation (27), $L$ is the Lagrangian variable, $\left\{q\right\}$ shows the vector of the generalized coordinates, and $\left\{Q\right\}$ is the vector of the generalized force. It is apparent that the complexity of the dynamic system relies on the length of the generalized coordinate vector. In other words, the computational burden for solving Equation (27) is related to the number of DOFs considered in the dynamic system. For a bottom-fixed offshore HAWT, there are usually 11 DOFs constituting the generalized coordinate vector [73]. Specifically, there are six DOFs for the blade, two for the tower, two for the drivetrain, and one for the pitch system to rotate the blade, and these are employed to describe the dynamics of the entire wind turbine system. For a floating HAWT, an additional six DOFs are introduced to describe the motion responses of the floating platform [74]. In order to model the wind turbine dynamics at a more realistic and reliable level, more DOFs could be included in the generalized coordinate vector shown in Equation (27). For example, two modes for the in-plane (edgewise) vibration of the blade, one mode for the flap-wise vibration of the blade, two for the side–side vibration of the tower, and two for the fore–aft vibration of the tower are added [75]. The Lagrangian variable $L$ is usually defined as the difference between the kinematic and potential energy of the dynamic system as functions of the generalized coordinates, as follows:

$$L=T(\dot{\left\{q\right\}}, \{q\left\}\right)-V\left(\right\{q\left\}\right)$$

(28)

In Equation (28), $T(\dot{\left\{q\right\}}, \{q\left\}\right)$ is the kinematic energy as the general coordinates $\left\{q\right\}$ and generalized velocities $\dot{\left\{q\right\}}$ and $V\left(\right\{q\left\}\right)$ gives the potential energy depending on the generalized coordinate.

Generally, the dynamic system governed by Equation (27) is constrained by either holonomic or non-holonomic equations of the generalized coordinates. While the holonomic constraints depend only on the generalized coordinates, the non-holonomic constraints are the linear combination of the generalized coordinates and generalized velocities.

Within the framework of the Euler–Lagrange approach, the dynamics of the HAWT are described via a series of independent generalized coordinates. The total kinetic energy of the HAWT is the summation of the kinetic energies of the blades, the tower, the nacelle, the drivetrain, and the floating platform, in the case of floating HAWT. More specifically, the total kinetic energy of the HAWT could be expressed as

$$T=\frac{1}{2}\sum _{i=1}^3{\int }_0^R{\mu }_b{v_{b,i}}^{2}dr+\frac{1}{2}{\int }_0^H{\mu }_t{v_t}^{2}dh+\frac{1}{2}{M}_n{\dot{q_{n,in}}}^2+\frac{1}{2}{M}_n{\dot{q_{n,out}}}^2$$

(29)

In Equation (29), ${\mu }_b$ and ${\mu }_t$ are the densities of the blade and the tower, $v_{b,i}$ and $v_t$ are the velocities of vibration along the blade and the tower, and $M_n$ gives the mass of the nacelle. $R$ is the radius of the blade and $H$ is the height of the tower. It is noted that the vibration velocities of various components in a HAWT are determined in accordance with the coordinate system. Taking a floating HAWT of the spar-type as an example [76] (as illustrated in Figure 6), the coordinate system contains at least eight different sets: (1) an inertial coordinate system $z_1{z}_2{z}_3$ that is established and fixed at a certain point at the sea surface; (2) a coordinate system $a_1{a}_2{a}_3$ that is attached to the tower base; (3) a coordinate system $b_1{b}_2{b}_3$ that is attached to the nacelle; (4) a coordinate system $c_1{c}_2{c}_3$ that is fixed to the low-speed shaft and rotated with it; (5) a coordinate system $d_1{d}_2{d}_3$ that is attached to the low-speed shaft and rotated with the nacelle in the yaw direction; (6) a coordinate system $e_1{e}_2{e}_3$ that is attached to the low-speed shaft and rotated in the tilt direction; (7) a coordinate system $g_1{g}_2{g}_3$ that is aligned with the rotor and rotated with it; and (8) a coordinate system $f_1{f}_2{f}_3$ that is attached to the blade and takes the cone angle of the blade into consideration [55,76,77,78]. Although the deformation and vibration of important components in a HAWT could be described by the local coordinate system, the kinetic energy calculation shown in Equation (29) requires the local motions to be transferred into the inertial coordinate system of $z_1{z}_2{z}_3$ via rotational matrices [55]. In fact, the velocities of the key components shown in Equation (29) can be obtained by transferring the corresponding position vectors from the local coordinate system into the inertial coordinate system and taking time derivatives [79,80]. The potential energies in association with the blades and the tower are the summation of the strain energy of the material, the centrifugal energy of the blade rotation, and the gravitational energy of the axial members. The potential energy could be calculated based on the linear Euler–Bernoulli beam theory and the small deformation assumption [16,66,70,79,80,81]. Substituting the kinetic and potential energies calculated into Equation (28), the Lagrange variable $L$ is ready for further use.

In addition to the Lagrange variable, the generalized force vector $\left\{Q\right\}$ shown in Equation (30) is determined by the virtual displacement method. Specifically, the virtual work ($\delta W$) for the dynamic system to deform in the generalized coordinate system with a series of virtual displacements ($\delta \left\{q\right\}$) is calculated, which yields the results of the generalized forces as

$$\{Q\}=\frac{\partial \left(\delta W\right)}{\partial \left(\delta \right\{q\left\}\right)}$$

(30)

In a dynamic system of the floating HAWT, the generalized forces are essentially linear combinations of the known external loads, such as aerodynamic loads, hydrodynamic loads, gravity, and hydrostatic restoring and damping forces.

5.2. The Kane Approach

In order to model the complex dynamics of a multi-body system, Kane equations have been proposed based on the simplification of the Gibs–Appel approach based on the Gibs function [82]. In fact, the Kane approach shows the dynamics via acceleration vectors instead of the scalar function of energy, as follows [83,84,85]:

$${\sum }_{i=1}^N\left(m_i{a}_{Gi}\frac{\partial V_{Gi}}{\partial \dot{q_k}}+\dot{H_{Gi}}\frac{\partial {\omega }_{Gi}}{\partial \dot{q_k}}\right)=Q_k, k=\mathrm{1,2},\cdots ,n$$

(31)

In Equation (31), $N$ is the amount of flexible/rigid bodies in the multi-body system, and $m_i$, $a_{Gi}$, $H_{Gi}$, and ${\omega }_{Gi}$ are the mass, acceleration of the centroid, angular momentum, and angular velocity of the $i^{th}$ body. $k$ shows the counts of the DOFs in the system and $Q_k$ is the generalized force corresponding to the specific DOF. The generalized force is calculated as

$$Q_k={\sum }_{i=1}^N\left(F_i\frac{\partial V_{Gi}}{\partial \dot{q_k}}+M_{Gi}\frac{\partial {\omega }_{Gi}}{\partial \dot{q_k}}\right)$$

(32)

In Equation (32), $F_i$ and $M_{Gi}$ are the external forces and moments acting on the $i^{th}$ body.

Compared to the Euler–Lagrange approach, the Kane approach is highly effective for modelling the complex multi-body dynamic system with non-holonomic constraints [84], and thus is suitable for large dynamic systems containing an excessive number of DOFs. Since the Kane approach deals with the individual components separately, the dynamic responses in association with the blades and tower can be estimated using the reduced-order Galerkin method, in which only a selective set of modes are considered to participate in the vibration. The Kane method is widely adopted by the numerical simulations of wind turbines. For example, the well-known numerical tool of FAST, an open-access code developed by NREL for the global dynamic simulation of onshore and offshore HAWT, models the aero-elastic loads and the dynamic responses of the wind turbine based on the Kane method. Table 3 presents a complete review of the existing approaches for HAWT modelling.

Table 3. Summary of the suggested HAWT modelling methods.

Ref	Method	Turbine	DOFs	Description
[86,87,88,89,90]	Euler–Lagrange	NREL 5 MW/Offshore monopole	14	six DOFs for the blades, six DOFs for the tower, and two DOFs for the drivetrain.
[74,91]	Euler–Lagrange	NREL 5 MW/floating	17	six DOFs for the blade, three DOFs for the tower, two DOFs for the drivetrain, and six DOFs for the spar type platform.
[77,78]	Kane	NREL 5 MW/Floating	22	nine DOFs for the blades (six for flapwise and three for edgewise), six DOFs for the platform, four DOFs for the tower (two for fore–aft and two for side–side), two DOFs for the drivetrain, one DOF for the nacelle yaw
[79,80,81]	Euler–Lagrange	NREL 5 MW/offshore monopile	12	six DOFs for the blades, two DOFs for the tower, and four DOFs for the monopile (translation and rotation)
[92]	Euler–Lagrange	NREL 5 MW/onshore	8	six DOFs for the blade (three for the edgewise and three for the flapwise) and two DOFs for the tower.
[93]	Euler–Lagrange	NREL 5 MW/onshore	3	three DOFs for the modelling of edgewise vibration of the blade.

Table 3. Summary of the suggested HAWT modelling methods.

Ref	Method	Turbine	DOFs	Description
[86,87,88,89,90]	Euler–Lagrange	NREL 5 MW/Offshore monopole	14	six DOFs for the blades, six DOFs for the tower, and two DOFs for the drivetrain.
[74,91]	Euler–Lagrange	NREL 5 MW/floating	17	six DOFs for the blade, three DOFs for the tower, two DOFs for the drivetrain, and six DOFs for the spar type platform.
[77,78]	Kane	NREL 5 MW/Floating	22	nine DOFs for the blades (six for flapwise and three for edgewise), six DOFs for the platform, four DOFs for the tower (two for fore–aft and two for side–side), two DOFs for the drivetrain, one DOF for the nacelle yaw
[79,80,81]	Euler–Lagrange	NREL 5 MW/offshore monopile	12	six DOFs for the blades, two DOFs for the tower, and four DOFs for the monopile (translation and rotation)
[92]	Euler–Lagrange	NREL 5 MW/onshore	8	six DOFs for the blade (three for the edgewise and three for the flapwise) and two DOFs for the tower.
[93]	Euler–Lagrange	NREL 5 MW/onshore	3	three DOFs for the modelling of edgewise vibration of the blade.

Figure 6. Coordinate sets of a floating HAWT model [94].

Figure 6. Coordinate sets of a floating HAWT model [94].

6. Conclusions

Along with the development of offshore wind energy, the numerical simulation of offshore wind turbines, especially the floating HAWT, plays an important role in the design, manufacture, installation, operation, and maintenance of the wind turbines. Given its importance, a considerable number of previous studies have been devoted to either developing computer-aided design tools or discerning dynamic characteristics of offshore wind turbines. Therefore, the present review summarizes relevant studies concerning the numerical simulation of the offshore wind turbines in the following aspects:

• The wind profile model serves as the fundamental basis for accurately simulating the dynamics of HAWTs. Empirical wind profile models, such as the ones employed in the widely acknowledged code of TurbSim, show that the vertical variations in wind speeds could be described using the logarithm law or the power law.
• The estimates of aerodynamic loads acting on HAWTs concern the interactions between turbulent winds and the blades and the tower. The blade element momentum theory and its corrections are primarily used with consideration of floating platform motion on induced velocity and inflow angle.
• The classic wave models and hydrodynamic load calculations are necessary for the numerical simulation of floating HAWTs. Specifically, the wave kinematics shown from the potential flow theory and the empirical stochastic models are key for determining the hydrodynamic loads.
• Reduced-order models are applied to discern the structural responses of HAWTs. Specifically, the Euler–Lagrange approach and the Kane approach exhibit different computational efficiencies and result reliability, and hence should be employed under different situations.

The transition to wind energy is in line with the urgent necessity for a worldwide move towards sustainable energy sources in response to the growing problems brought on by climate change and the finite nature of conventional fossil fuels [95,96]. The market has already acknowledged the long-term advantages of renewable energy, but the use of renewable energy also faces considerable obstacles, such as the enormous amount of initial investment required, the uncertainty of the power generation, and the intricate dynamics of funding distribution [97]. The need to counteract climate change, lower carbon emissions, and promote energy independence makes the switch to wind energy inevitable. As recently reported by the International Energy Agency (IEA), switching to wind energy is not only a need but also a choice that must be made in order to meet global climate targets [98].

Challenges arise from the complexity of international markets and the variety of energy dependences. The shift especially entails a number of financial challenges, from the large capital costs associated with developing infrastructure to the unpredictability of investment returns and financing distribution. A number of studies have already tried to address these issues, which is where numerical simulations have come in quite handy. As shown by Li et al. [99] and Miller et al. [100], numerical simulations are a powerful tool for reducing financial uncertainty related to wind energy transitions. Firstly, computer simulations offer a virtual testing ground in which wind turbine designs can be optimized. Secondly, numerical simulations make it possible to thoroughly examine the dynamic responses of wind turbines, enhancing the estimation of their operation and maintenance costs. Considering the importance of numerical simulation in financially assessing the deployment of offshore wind turbines, the relevant studies deserve a dedicated review, which is planned as a subsequent work after the present study.