Analysis of the blade aeroelastic effect on the floating offshore wind turbine wake in a focusing wave with a hybrid model

Abstract

With the recent trend toward larger floating offshore wind turbines (FOWT), the influence of blade deformation becomes more obvious. Modelling this fully coupled fluid-structure interaction (FSI) system, especially at extreme sea level with large-scale motion of the entire system, is an essential challenge. A hybrid numerical model integrating the potential-viscous flow model (qaleFOAM), the unstable actuator line method (UALM), the Legendre spectral finite element model (BeamDyn), and the Lumped Mass Mooring Model (MoorDyn) is employed in the current research. This model is capable of effortlessly and accurately reflecting the blade aeroelastic effect on the dynamics of a FOWT system under the circumstances of a focused wave and uniform wind. The result reveals that the aeroelastic impacts have been identified to not only increase the fatigue damage on the turbine blades but also affect the wake field distribution. Based on the temporal and spatial distribution of the velocity fields, it is found that the coupled effects of the large wave elevation and the FOWT high-frequency motions will disturb the velocity in the wake region. The slight variation of the wind speed in the wake region will be exacerbated by about 0.21% to 18.6%, and it transforms the shape of the wake region when the blade aeroelastic effect participates in the coupled effect. This phenomenon may further affect the downstream FOWT system features through the upstream FOWT wake field.

KEYWORDS:

• FOWT
• hybrid model
• aeroelastic response
• fatigue damage
• wake field distribution

1. Introduction

There are numerous hurdles for wind turbines growing in deep water and large-scale deployment. Larger blades will boost wind farm production capacity per unit area while lowering the cost of supporting infrastructure. The 5 MW reference wind turbine set by NREL (National Renewable Energy Laboratory), for example, has more than 60 m of blades (Jonkman, 2007). The DTU (Technical University of Denmark) 10 MW wind turbine has a rotor diameter of about 178 m (Bak et al., 2013). The deflection and vibration of the turbine blade can interfere with the fatigue life and wake flow characteristics (Yu et al., 2020). Especially for the floating offshore wind turbine (FOWT), its dynamic behaviours are particularly challenging in the diverse and variable marine environment. All the coupled physics should be accounted for, including the wind turbine aerodynamic loads, turbine blades’ deformation, the dynamic behaviour of the floating platform and mooring system, and the wake flow field behind the turbine rotor. It is critical in the applied sector to appropriately resolve these interactions (Hao et al., 2023).

The influence of the wind turbine wake and the nonlinear problems related to the supporting platform are disregarded in the Blade Element Momentum (BEM) method or the potential flow theory, which is crucial for determining and assessing FOWT system characteristics (Sebastian & Lackner, 2012). Therefore, the application of computational fluid dynamics (CFD) in FOWT system modelling is a sensible choice because it solves the above problems and leads to more accurate findings. Due to the direct solution of the Navier-Stokes (N-S) equations (Zheng, Ma, et al., 2014; Zheng, Shao, et al., 2017), the comprehensive flow field information, including the wind turbine wakes and the enormous free surface elevation, can be obtained easily. As a result, there has been a lot of analysis based on the CFD tools about FOWT systems. For the wind turbine aerodynamics, Wen et al. (2017) and Tran and Kim (2015) investigated the aerodynamic loads on the wind turbine with given motion. Apsley and Stansby (2020) compared the unsteady thrust of the oscillating wind turbine between BEM method and AL (actuator line) method. Xie (2021) improved the velocity sampling and projection in AL model to better predict the wind turbine performance. For the blade aeroelastic response, Ma et al. (2019) developed the Actuator Line Finite-element Beam method to discuss the influence of the blade deformation on the fixed wind turbine wake characteristics. Liu et al. (2019) combined the CFD and MBD (multi-body dynamics) approach to analyze the dynamic response of the NREL wind turbine. Zheng et al. (2023) coupled STAR CCM+ and ABAQUS to study the influence of turbulence intensity of a wind turbine. For the wind turbine wakes, Miao et al. (2017) studied the wake characteristics interaction between two fixed wind turbines. Ravensbergen et al. (2020) adopted the AL method in the residual-based variational multiscale fluid modelling framework to study the airflow around the NREL wind turbine. Yuan et al. (2023) proposed CFD (ALM)-IDWM to study the wind turbine wake features, which achieved enhanced computational efficacy. It should be noted that such high-fidelity approaches typically require large computing costs. It is inappropriate for further modelling the FOWT arrays.

On the other hand, it is a challenge for the FOWT system modelling because the marine circumstances are so complicated and unpredictable. To overcome such a fully coupled FOWT system, initial research usually simplifies some conditions, i.e. it is focused on the aeroelastic of the fixed wind turbine (Li et al., 2015) or considers the coupling of regular waves and uniform winds with rigid blades (Yu et al., 2022). As technology has advanced, an increasing number of studies have added complex conditions to FOWT simulation. The investigations on FOWT dynamics gradually raise the level of satisfaction to realise the physics situation by taking into account specific conditions, such as the Atmospheric Boundary Layer (ABL), turbulent wind, focused wave, or irregular wave conditions (Zhou, Xiao, Liu. et al., 2019; Zhou, Xiao, Peyrard, et al., 2021; Band et al., 2021; Maduka & Li, 2021). Especially for the FOWT arrays, the accuracy of the downstream FOWT dynamics is dependent on the upstream FOWT’s prediction and its wake distribution.

To overcome the above-mentioned issues, a partially simplified technique, the hybrid model utilised in the present investigation, has been established in our previous work (Yu et al., 2023). In this model, the potential-viscous flow model (qaleFOAM) (Ma & Yan, 2006; Yan et al., 2020), which handles the hydrodynamics and wake flow of FOWT, can efficiently reduce the size of the viscous flow domain. To avoid the boundary layer grid surrounding the blade surface and the sliding or overset mesh methodology for blade rotation, the unsteady actuator line method (UALM) is built to replace the influence of the physical turbine blade on the flow field (Yu et al., 2018). To calculate the structural dynamics of the blades and mooring system, the BeamDyn (Wang et al., 2016) and MoorDyn (Hall & Goupee, 2015) are introduced. Several cases have been used to validate and verify the capability of the present model with some simplification. The predictions from this model agree reasonably well with the experimental data and other numerical results. It has been demonstrated that this model is capable of handling coupled aero-hydro-mooring-elastic FOWT systems and is more efficient than the CFD approach with real blades (Yu et al., 2023).

The current research will expand on the conditions mentioned above by involving the focused wave and elastic blade. To guide future studies on the FOWT arrays modelling, the influences on the FOWT wake field are investigated. Furthermore, the OC4 DeepCwind semi-submersible NREL 5 MW floating offshore wind turbine (Coulling et al., 2013), a popular FOWT type, is utilised. The focused wave and incoming wind are set up according to Beaufort Sea State 6. All dynamic behaviours are considered by the hybrid model (Yu et al., 2023), including the wind turbine aerodynamics, the blade structural dynamics, the floating platform motion responses, and mooring dynamics. As these dynamic behaviours have been validated well with the experimental and numerical results from other publications (Yu et al., 2023), the grid convergence test and uncertainty analysis will be further demonstrated. Then, the influence of blade deformation on the wind turbine wake features under such extreme sea state will be addressed and discussed.

2. Numerical methodology

2.1. Governing equations

2.1.1. Hydrodynamic model

In the hydrodynamic model, the potential-viscous flow hybrid approach is employed. The majority of the wave field is handled in the FNPT domain, where the fluid is assumed to be incompressible, inviscid, and irrotational. QALE-FEM (Ma & Yan, 2006; Yan et al., 2020) is employed for describing wave generation and propagation. The governing equation is defined by velocity potential Φ as follows: (1) $\begin{array}{r}{\mathrm{\nabla }}^2\mathrm{\Phi }=0\end{array}$(1) where the corresponding boundary conditions satisfy $\frac{\mathrm{\partial }\mathrm{\Phi }}{\mathrm{\partial }n}=\mathit{n}\cdot \mathit{U}$on the wavemaker, $\frac{\mathrm{\partial }\mathrm{\Phi }}{\mathrm{\partial }n}=0$ on fixed boundaries, $\frac{D\mathit{r}}{Dt}=\mathrm{\nabla }\mathrm{\Phi }$ and $D$ on the free surface, where n and U are the normal vectors of the wavemaker boundary and the fluid velocity, respectively; r is the position vector (x, y, z) of the free surface; and g is the gravity acceleration.

In the N-S domain, a small region near FOWT is considered. The two-phase flow model is used through the VOF (volume of fluid) approach (Hirt & Nichols, 1981) to capture the free surface. The governing continuity, momentum, and volume fraction (α) equations are summarised below. (2) $\begin{array}{rl}& \frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \rho \mathit{U}=0\end{array}$(2) (3) $\begin{array}{rl}& \frac{\mathrm{\partial }\rho \mathit{U}}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left[\rho \right(\mathit{U}\mathbf{-}{\mathit{U}}_g\text{)}\mathit{U}]\\ & =-\mathrm{\nabla }p-\rho \mathit{g}+\mathrm{\nabla }\cdot \left(\mu \mathrm{\nabla }\mathit{U}\right)+\mathrm{\nabla }\mathit{U}\cdot \mathrm{\nabla }\mu +\rho {\mathit{f}}_{UALM}\end{array}$(3) (4) $\begin{array}{rl}& \frac{\mathrm{\partial }\alpha }{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \alpha (\mathit{U}-{\mathit{U}}_g)+\mathrm{\nabla }\cdot {\mathit{U}}_c\alpha (1-\alpha )=0\end{array}$(4) in which U_g is the velocity of the computational grid; p is the dynamic pressure; f_UALM denotes the source item solved by the aerodynamic model and will be discussed in the following subsection; U_c is the artificial velocity field; α = 0 and α = 1 represent the liquid phase and gas phase, respectively; 0 < α < 1 means the free surface between air and water; ρ = ρ_air (1-α) + ρ_water α is the fluid density; μ = μ_air (1-α) + μ_water α is the dynamic viscosity; and the subscripts ‘air’ and ‘water’ refer to the quantities for air and water, respectively.

2.1.2. Aerodynamic model

In the aerodynamic model, the unsteady actuator line method (UALM) substitutes the real blade with an actuator line with the body forces distributed on it. The actuator line is then divided into a number of distinct actuator elements. To describe the influence of the aerodynamic loads on the flow field, they are projected onto the background mesh and act as the body forces on the right side of the momentum equations, as shown in Equation (3). The range of the projection is determined by the Gauss function (Yu et al., 2023), which can smooth the distribution of the forces and avoid the numerical singular problem. (5) $\begin{array}{r}{\mathit{f}}_{UALM}=-F_{tip}{\mathit{f}}_{aero}\otimes {\eta }_{ϵ}\left(d\right)\end{array}$(5) where d is the distance from the measured point to the aerodynamic centre; F_tip is the tip loss coefficient and is determined by the Prandtl-Glauert tip loss model (Shen et al., 2005); f_aero is the aerodynamic load of each actuator element, which can be described as ${\mathit{f}}_{aero}=(\mathit{L},\mathit{D})=\frac{1}{2}{\rho }_{air}{A}_{elem}|{\mathit{U}}_{rel}{|}^2(C_l{\mathit{e}}_L+C_d{\mathit{e}}_D)$, in which L and D denote the aerodynamic lift and drag, respectively; ρ_air is the air density and is 1.225 kg/m³ in this research; A_elem is the blade element planform area; C_l and C_d are the lift and drag coefficients, respectively; e_L and e_D are the unit vectors in the lift and drag directions; U_rel is the relative velocity of each actuator element and is defined as ${\mathit{U}}_{rel}={\mathit{U}}_{in}+\mathit{\Omega }\times \mathit{r}-{\mathit{U}}_{rot}+{\mathit{U}}_M+{\mathit{U}}_{def}$, where U_in represents the inflow velocity component that is vertical to the rotor plane; U_rot is the flow velocity induced by the rotating blade; Ω × r is the speed of this airfoil caused by the blades’ rotation; U_M is the resultant velocity induced by the motions of the FOWT, which can be further divided into translational velocity and linear velocity caused by rotation around the FOWT centre; and U_def is the velocity because of the blade deflection.

Through above equations, the UALM takes into consideration both FOWT motion and the flexible blade effect. All actuator line elements must be rotated around the hub to their positions in the following time step after the flow field at the current time step has been solved. The hub and tower of the wind turbine in this research are also replaced by the actuator line except the blades. The airfoils for these two parts are treated as equivalent cylinders. However, they are treated as rigid bodies in the present study. Moreover, it is significant for the modelling of the FOWT system in the extreme condition that the motion amplitude of FOWT is large enough to cause some airfoils’ attack angles to enter the stall range. The dynamic stall model (Sheng et al., 2008) must be considered here to correct the lift and drag coefficients.

2.1.3. Structural model

The turbine blade structural dynamics are solved by the open source code BeamDyn (Wang et al., 2016), and it uses the geometrically exact beam theory (GEBT) with the Legendre-spectral-finite-element (LSFE) implementation. The LSFE is a high-order finite element and more computationally efficient than low-order finite elements. It can deal with geometric nonlinear beam problems for dynamic analyses and take full elastic coupling effects into account. Because of the state-space form, the code can be easily modified and coupled with other code.

The motion of each beam element is governed by (6) $\begin{array}{rl}& \frac{\mathrm{\partial }{\mathit{h}}_e}{\mathrm{\partial }t}-\frac{d{\mathit{F}}_s}{d\mathit{s}}={\mathit{f}}_s\end{array}$(6) (7) $\begin{array}{rl}& \frac{\mathrm{\partial }{\mathit{g}}_e}{\mathrm{\partial }t}+\frac{\mathrm{\partial }\tilde{\mathit{u}}}{\mathrm{\partial }t}{h}_e-\frac{d{\mathit{M}}_s}{d\mathit{s}}-({\tilde{\mathit{x}}}_0+\frac{\mathrm{\partial }\tilde{\mathit{u}}}{\mathrm{\partial }\mathit{s}}){\mathit{F}}_s={\mathit{m}}_s\end{array}$(7) where, h_e is the linear momenta of the beam element; g_e is the corresponding angular momenta; F_s is the force for the beam section; M_s is the corresponding moment; u is the displacement at initial undeformed time; x₀ is the vector of the initial position; the superscript ‘∼’ is regarded as the skew-symmetric tensor; f_s is the external distributed force on the beam element; and m_s is the corresponding external moment. They are converted into distribution loads by dividing the distance between two adjacent actuator elements from the aerodynamic loads f_aero in Equation (5). Linear interpolation is then used to transform the forces acting from the actuator element to structural element. In addition, the motion velocities of the FOWT from the hydrodynamic solver are also needed to update the boundary condition in the structural solver.

The Newton–Raphson method is used to solve the governing equations of the entire beam. It is discretized into multiple Gauss-lobatto-Legendre points by the LSFE method. the generalized-α scheme with second-order accuracy implemented for the time integration. Furthermore, a linearisation process is needed because the original governing equations are nonlinear. More detailed deductions can be found in the research of Wang et al (Wang et al., 2016).

2.1.4. Mooring line model

The mooring dynamics are solved individually and act as the restraints for the FOWT system. The dynamic model, MoorDyn (Hall & Goupee, 2015), considers the incorporates axial elasticity, hydrodynamic loading by Morison's equation, and bottom contact. It adopts a lumped-mass approach to discretize the mooring cable. The mooring line is broken up into several evenly sized sections, and the weight of each part is lumped at one node. According to the force equilibrium equations, the forces and motion of each node should satisfy as the following, (8) $\begin{array}{r}[{\mathit{m}}_i+{\mathit{a}}_{pi}+{\mathit{a}}_{qi}]\frac{{\mathrm{\partial }}^2{\mathit{r}}_i}{\mathrm{\partial }{t}^2}={\mathit{T}}_i+{\mathit{C}}_i+{\mathit{W}}_i+{\mathit{B}}_i+{\mathit{D}}_{pi}+{\mathit{D}}_{qi}\end{array}$(8) where, subscript i represents node i; a_p is the added mass in the transverse direction; a_q is the added mass in the tangential direction; m_i means the mass of node i; T_i denotes the tension in cable; C_i represents the internal damping force,; W_i means the net buoyancy of node i; B_i denotes the force when the node contacts the seabed; D_pi is the transverse drag force; and D_qi is the and tangential drag force.

2.2. Coupling strategies

The variable of data exchange between each model is marked red in Figure 1. A two-way coupling approach is used in the following processes, which are essential for capturing the coupling effect of the FOWT modelling.

Figure 1. Flow chart of the floating offshore wind turbine hybrid model.

Between the aerodynamic and hydrodynamic solvers, the motion data of the FOWT is essential for the coupling strategy. If the motion velocities of the FOWT are not considered, the impact of the platform’s hydrodynamic response on aerodynamic performance cannot be captured. The aerodynamic solver transfers to the hydrodynamic solver with aerodynamic forces on turbine blades and tower. The hydrodynamic solver feedbacks the 6-DOF displacement and motion velocities of the FOWT to the aerodynamic solver. The displacements are used for updating the actuator line element position to keep the wind turbine following the FOWT motion. And the motion velocities are considered when calculating the relative velocities for each airfoil in the present improved UALM (see Section 2.1).

Between aerodynamics and structural models, BeamDyn library shares the blade deformation for updating the actuator line element position and corresponding velocities to the angle of attack process calculation for each airfoil in the present ALM model. The BeamDyn code is translated and edited as a dynamic link library. And it is called here via a proprietary namespace inheriting the OpenFOAM in-house Foam namespace. It should be noted that the structural solver typically has a smaller time step than the fluid solver. If the larger time step is adopted, the calculation accuracy of the structural solution will be low accurate or even lead to divergent results. Or if a smaller time step is used, the computational cost will increase. Therefore, an inner iteration is added to solve this problem. In addition, the motion velocities of the FOWT from the hydrodynamic solver are also needed to update the boundary conditions in the structural solver.

The FOWT motion uses the body reference system in the hydrodynamic model, but the mooring line solver adopts its inertial frame of reference for its governing equation. Hence, the displacement and motion velocities of the FOWT need to be transferred to the MoorDyn code at each time step. The hydrodynamic solver shares the 6-DOF displacement and motion velocities of the FOWT. According to the motion data, the location of each fairlead of the floating platform is calculated. Therefore, the present state of the mooring line can be updated through the new fairleads and seabed position. The mooring dynamic solver (MoorDyn) returns the mooring line tension as the restraint forces acting on the fairleads of the platform.

2.3. Wave generation and propagation

The focused wave used in this work is generated by the second order wavemaker (Schaffer, 1996; Yan et al., 2015) far away from the FOWT with a group of wave components. Each component is determined by the JONSWAP spectrum (Katsidoniotaki, 2019). The spectrum S(f_i) utilised in this paper is expressed as, (9) $\begin{array}{rl}S\left(f_i\right)& =a{H}_s^2{f}_p^4{f}_i^{-5}exp\left[-\frac{5}{4}{\left(\frac{f_p}{f_i}\right)}^4\right]{\gamma }^{\text{β}}\end{array}$(9) (10) $\begin{array}{rl}a& =\frac{0.06238}{0.230+0.0336\gamma -\frac{0.185}{1.9+\gamma }}\left(1.094-0.01915\ln \gamma \right)\end{array}$(10) (11) $\begin{array}{rl}\text{β}& =\exp \left[\frac{-{(f-f_p)}^2}{2{\sigma }^2}{f}_p^2\right]\end{array}$(11) (12) $\begin{array}{rl}\sigma & =\{\begin{array}{ll}0.07& f\le f_p\\ 0.09& f>f_p\end{array}\end{array}$(12) where H_s is the significant wave height; f_p represents the peak frequency; f_i denotes the wave frequency; the subscripts i represent the component number; and γ is the peak enhanced factor and is 3.3 in this paper.

The generated wave propagates into the N-S domain through the wave inlet boundary and the attached relaxation zone technique. At the wave inlet boundary of the N-S domain, the fluid velocity u is corrected by. (13) $\begin{array}{r}\mathit{u}(x,y,z)=\{\begin{array}{cc}\mathrm{\nabla }\mathrm{\Phi }(x,y,z)& z\le \eta \\ (1-R_z)\mathrm{\nabla }\mathrm{\Phi }(x,y,\eta )\\ +R_z{\mathit{u}}_w(x,y,z)& z>\eta \end{array}\end{array}$(13) in which, η is the wave elevation; R_z is the weighting function, ranging from 0 to 1, and it secures a smooth transition of the fluid velocity from the water phase to the air phase; and u_w is the wind speed.

In the relaxation zone, the N-S solution f (velocity u and pressure p) is corrected by f (x,y,z) = f_NS (x,y,z) w + f_FNPT (x,y,z) (1-w), where w is the weighting function that is 0 on the inlet boundary of N-S domain and 1 on the other boundary of the relaxation zone, and the subscripts NS and FNPT stand for N-S solution and QALE-FEM solution, respectively.

3. Numerical model

The NREL 5 MW OC4 semi-submersible floating offshore wind turbine is adopted. The details of the geometry parameters for the platform and the mooring model, air foil data for the actuator line model, structural properties for the elastic model can be found in the NREL reports (Jonkman, 2007) or our previous studies (Yu et al., 2018; Yu, Hu, et al. 2020; Yu, Yan, et al. 2022; Yu, Ma, et al. 2023). Figure 2 depicts the whole simulation domain, which includes the detailed size of the FNPT and N-S domains. The simulation FNPT domain has a total length of 3000 m, and the distance from the potential wavemaker to the inlet boundary of N-S domain is 950 m. The N-S box domain is defined as 650 m × 200 m × 460 m. The inlet and outlet boundaries are about 1D and 4D of the turbine rotor away from the centre of the platform. The 6-DOF movement of the FOWT is dealt with by the dynamic mesh technique introduced above, and the size of the mesh morphing zone is 100 m × 150 m × 200 m. The thickness of the mesh deformation zone is 80 m. The water depth is 200 m.

Figure 2. Sketch of the whole computational domain.

The incident focused wave has a significant wave height of 6 m and a wave peak period of 10 s, which is determined by the Beaufort Sea State 6. Its corresponding wind speed range is from 10.8 m/s to 13.8 m/s, and the rated wind speed of 11.4 m/s is just located in this range. Due to the transient effect, it takes some time for a FOWT system to reach a steady state (Yu et al., 2022). As a result, the focused position is set at x = 1100 m, where the FOWT system centre is located, and the focused time is 150 s in the N-S domain. The frequency band is from 0.06 Hz to 0.2 Hz, with 63 wave components. For this case, the time history of the wave elevation at the focusing position without considering FOWT is shown in Figure 3, and the results are close to those of the qaleFEM. From Figure 3, the maximum wave height H_max is about 8.96 m. The time step is 0.025 s, which corresponds to the blade rotation of approximately 2 degrees (12.1 rpm) in the uniform wind speed of 11.4 m/s. On this basis, the time step in the structural solver is 0.005 s. As is shown in Figure 2, the turbine blade, hub, and tower are replaced by several lines because they are solved by UALM and do not exist in the simulation domain.

Figure 3. Comparisons of the wave elevation time history at focusing location.

Figure 4 presents the computational mesh of the N-S domain, which is generated by the built-in mesh generation tool in OpenFOAM 5.x. The size of the mesh near the wind turbine is 1 m, and the coarse grid size of 2 m is extended to the outlet of the N-S domain. Therefore, the background grid size is 4 m. The grid size near the free surface will be discussed in the following section. The total mesh for the present computation is almost 5,749,412. Furthermore, the combination algorithm PIMPLE (PISO + SIMPLE) is utilised to solve the pressure-velocity coupling. The LES with WALE (Wall-Adapting Local Eddy-viscosity) model is used to deal with the flow condition at a high Reynolds number around the wind turbine.

Figure 4. Grid structure for the N-S domain.

4. Validation and verification

The accuracy, reliability, and efficiency of the hybrid model have been demonstrated in several cases compared with the numerical and experimental data, including the aerodynamic loads of the wind turbine, the blade deformation, the hydrodynamic responses of the supporting platform, and the mooring line tensions (Yu et al., 2023). On this basis, the grid convergence test and uncertainty analysis will be further verified in this section.

The physical conditions of the present verification are the same as in the previous study (Yu et al., 2022) and are not described again. Three different meshes are generated for the grid sensitivity test of the free surface mesh, as listed in Table 1. The refinement ratio between neighbouring levels is set to$\sqrt{2}$. The time history of the wave elevation for different grid levels is shown in Figure 5. It seems that the predicted results are not sensitive to the present three mesh densities. Meanwhile, considering the computational efficiency, the medium-sized of grid will be utilised in the following simulation.

Figure 5. Free surface mesh sensitivity tests for focused wave generation.

Table 1. Mesh strategy for sensitivity tests near the free surface under the focused wave.

Grid level	Refinement	λHmax/Δx	Hmax/Δz
1	Fine	406	26
2	Medium	287	18
3	Coarse	204	13

The uncertainty is analyzed through a progress recommended by ITTC (2017) and listed as the following steps:

1. Make sure the convergence ratio R_ζ through Equation (14). (14) $\begin{array}{r}{R}_{\zeta }=\frac{{ϵ}_{\zeta 21}}{{ϵ}_{\zeta 32}}\end{array}$(14) where ϵ_ζ21 and ϵ_ζ32 are the difference of the physical quantity f between neighbour levels, which are equal to f_ζ2 – f_ζ1 and f_ζ3 – f_ζ2, respectively.

2. Calculate the accuracy order p_ζ and the first order Richardson extrapolation δ*_REζ (15) $\begin{array}{rl}{p}_{\zeta }& =\frac{\ln \left({ϵ}_{\zeta 32}/{ϵ}_{\zeta 21}\right)}{\ln r_{\zeta }}\end{array}$(15) (16) $\begin{array}{rl}{\delta }_{RE\zeta }^{\ast }& =\frac{{ϵ}_{\zeta 21}}{r_{\zeta }^{p_{\zeta }}-1}\end{array}$(16) where r_ζ is refinement ratio between neighbour levels.

3. The uncertainty coefficient U_ζ can be computed by Equation (17)

(17) $\begin{array}{r}{U}_{\zeta }=(|C_{\zeta }|+|1-C_{\zeta }|)|{\delta }_{RE\zeta }^{\ast }|\end{array}$(17) where C_ζ is the correction factor and defined as Equation (18)

(18) $\begin{array}{r}{C}_{\zeta }=\frac{r_{\zeta }^{p_{\zeta }}-1}{r_{\zeta }^2-1}\end{array}$(18) Three sets of cell sizes are studied and summarised in Table 2. The refinement ratio of the grid, not only in the wake region but also near the turbine blades, is kept constant at$\sqrt{}2$. The corresponding cell numbers of levels 1, 2, and 3 are about 8.73, 3.33, and 1.23 million, respectively. In these cases, the FOWT system is subjected to the focused wave, tested above, and the uniform wind of 11 m/s. It should be noted that the system needs some time to balance in the form of free decay under uniform wind loading, but the focused time is too early. However, it saves computational cost and increases the maximum amplitude of the motion response, which examines the capability of the dynamic mesh motion.

Table 2. Mesh strategy of FOWT system for uncertainty analysis and convergence test.

Grid level	Refinement	Horizontal mesh size near free surface (m)	Mesh size in the wake region (m)	Total grid numbers
1	Fine	2.839	0.71	8.73 M
2	Medium	4.014	1.00	3.33 M
3	Coarse	5.676	1.41	1.23 M

The test results are presented in the maximum amplitudes of pitch (f₁), surge (f₂), heave (f₃), thrust (f₄), and power (f₅) from time history near the focused time. Table 3 illustrates the errors and uncertainty of grids. It indicates that all the uncertainties in properties are less than 2.3%, which indicates the results are acceptable. From the table, the accuracy except surge and pitch motion is more than third order. The pitch and surge motion have nearly first order. In addition, these five parameters are nondimensionalized and plotted in Figure 6 with different grid numbers. It can be concluded that the difference in results between the medium and fine refinement grids is almost negligible. Taking the computational cost into account, the medium size (level 2) should be considered suitable, and this form of mesh generation will be adopted in the following study.

Figure 6. Convergence tests of grid: (a) pitch; (b) surge; (c) heave; (d) thrust; (e) power.

Table 3. Results of error and mesh uncertainty analysis.

Properties	Rζ	pζ	δ*REζ	Cζ	Uζ	Uζ / fζ1
Pitch (f1)	0.636	1.305	0.091	0.749	0.091	1.563%
Surge (f2)	0.956	0.730	−0.468	0.068	0.468	2.253%
Heave (f3)	0.220	4.366	−0.007	1.050	0.008	0.062%
Thrust (f4)	0.015	12.186	0.006	1.000	0.006	0.002%
Power (f5)	0.290	3.573	10.223	1.079	11.831	1.022%

5. Results and discussions

5.1. Wind turbine aeroelastic responses

There are mainly two aspects of the blade’s aeroelastic influence on the aerodynamics. The deformations of the turbine blades transform their location, which leads the flow field near the blade to be different from that of the rigid blade. On the other hand, the deflection velocity affects the relative velocity in each actuator element (referred to Equation (5)). In this section, the FOWT with the rigid blade under the same focused wave and uniform wind is used for comparison. Figure 7 shows the time history and power spectrum of the aero thrust and power near the focused time (90 s–210 s). The power spectra are computed by the Periodogram, which estimates the power from the amplitude of the Fourier transformed data. Besides, the method of time-integral squared amplitude (TISA) is adopted to normalise the power spectrum in each domain. Table 4 summarises the comparison of the maximum and minimum values in Figure 7. We can find that the overall impact of blade aeroelasticity on thrust and power is not significant. Especially for the extreme values, their differences are less than 3%. On this basis, the blade elasticity effect is further weakened after being transmitted to the platform through the tower, and such an insensitive influence (see Figure 8) on the platform hydrodynamics will not be further discussed in the present study. However, the spectral power of the elastic blade is significantly higher than that of the rigid blade by about 10² orders, referring to Figure 7(b and d). It indicates that the aeroelasticity of the blade still plays a role in the wind turbine's aerodynamics.

Figure 7. Comparison of the floating offshore wind turbine system aerodynamics between elastic and rigid blade: (a) time history of thrust; (b) spectra of thrust; (c) time history of power; (d) spectra of power.

Figure 8. Comparisons of the floating offshore wind turbine hydrodynamic responses: (a) surge; (b) pitch.

Table 4. Comparison of the thrust and power with and without blade elastic.

Properties	Maximum thrust (kN)	Minimum thrust (kN)	Maximum power (MW)	Minimum power (MW)
Elastic blade	763.54	598.93	5.83	3.56
Rigid blade	775.39	603.77	5.79	3.46
Difference (error)	11.85 (1.56%)	4.84 (0.81%)	0.04 (0.69%)	0.10 (2.8%)

To quantify the aeroelastic effect on aerodynamics, the equivalent fatigue load on turbine blades is discussed. The concept of equivalent fatigue loads D_el (Freebury & Musial, 2000; Meng et al., 2018) is defined as below, (19) $\begin{array}{r}{D}_{el}=\frac{{\sum }_{i=1}^{N_t}{n}_i{\left(F_i\right)}^m}{N_t}\end{array}$(19) in which i is the load case number; n_i is the current number of cycles under the loads F_i acting on the turbine blades, and it can be calculated by the rain-flow counting algorithm (Yeter et al., 2015); N_t is the total number of load cases; m is the fatigue damage coefficient, which is 10 in this research.

Analyzing Table 5, the equivalent fatigue load of the elastic turbine blade is obviously increased by nearly 33% compared with the rigid turbine blade, which is significant for the fatigue analysis and cannot be neglected.

Table 5. Comparisons of the turbine blade equivalent fatigue load with and without blade elastic effects.

Properties	equivalent fatigue load (kN)
Elastic blade	128.24
Rigid blade	86.29
Difference (error)	41.95 (32.71%)

The blade tip flapwise vibration is then investigated. To compare with the present results, an additional case is modelled in which the FOWT is fixed under the same wind and wave conditions. As shown in Figure 9, the motion of the FOWT system disturbs the blade deformation, which is irregular compared with the fixed situation. Figure 9(b) illustrates the corresponding spectra. Several frequencies of the FOWT coupling effect can be captured. Peak A (0.01 Hz) and peak B (0.0366 Hz) are the surge and pitch natural frequencies of the FOWT system. Peak C (approximately 0.1 Hz) corresponds to the incident wave peak frequency. Peak D (0.2 Hz) represents the wind turbine rotation frequency. In addition, it is found that the average of the blade tip flapwise deflection varies insignificantly after considering the system motion. Although the wave elevation achieves its peak at t = 150 s, the blade structural response reaches its peak delayed. In their corresponding power spectra, the amplitude of blade rotation frequency is nearly 10¹ orders higher than that of other frequency points. This phenomenon indicates that the aerodynamic performance of the wind turbine plays a primary role and the system’s hydrodynamic response plays a second primary role. Moreover, both influences cannot be ignored.

Figure 9. Comparisons of the blade tip flapwise deflection between fully coupled and fixed FOWT: (a) time history; (b) the corresponding spectra.

5.2. The raised local velocity in the wake field

Since the influence between the upstream and downstream wind turbines is mainly reflected through the velocity deficit, the velocity distribution of the wake field will be discussed in this section. Figures 10 and 11 are the contours of the velocity for the fixed wind turbine and the FOWT with the elastic blade, respectively. To observe the influence of the large wave elevation and the FOWT motion, the sampled time doesn’t use one cycle time of the rotor rotation but according to the time history shown in Figures 3 and 8. There is an obvious velocity loss in the wake region due to the utilisation of wind energy by the wind turbine. Compared with the fixed wind turbine, the local velocity in certain regions behind FOWT has a significant increase around t = 144 s, t = 150 s, and t = 160 s, as shown in Figures 10 and 11(b–e). From Figure 8, the high frequency motion of FOWT is found during the occurrence of the wave focus. It leads the vibration of the wake region outline behind the turbine rotor (Figure 11(c and e)). Moreover, the wave elevation at the FOWT position reaches extremum at this moment, which further influences the wake flow field by reducing the distance between the turbine rotor and the free surface. The flow velocity is inversely proportional to the cross-sectional area. Therefore, the local velocity is raised and the vortex shedding from the blade tip is also extruded. Consequently, the phenomenon of the increased local velocity behind the turbine rotor can be concluded by the coupled effect between the high frequency motion of FOWT and the large wave elevation. This effect propagates in the downstream direction and increases the flow field velocity in a wide range of areas greater than those not disturbed by the coupled effect. As summarised in Table 6, one may observe that only large wave elevation (fixed FOWT with the same wave condition in Figure 10) or only FOWT motion (fully coupled motion of FOWT before and after the focused time in Figure 11(a and f)) cannot cause the raised local velocity.

Figure 10. Velocity fields of the fixed wind turbine with elastic blade at midplane: (a) t = 100 s; (b) t = 144 s; (c) t = 150 s; (d) t = 160 s; (e) t = 170 s; (f) t = 300 s.

Figure 11. Velocity fields of the fully coupled wind turbine with elastic blade at midplane: (a) t = 100 s; (b) t = 144 s; (c) t = 150 s; (d) t = 160 s; (e) t = 170 s; (f) t = 300 s.

Table 6. The summary of the raised local velocity reason in the wake field.

Large wave elevation	High frequency motion of FOWT	Raised local velocity	Evidence
√	×	×	Figure 10
×	√	×	Figure 11(a and f)
√	√	√	Figure 11(b–e)

5.3. The blade aeroelastic effect on the wake velocity distribution

As discussed in Figure 8, the motions of the platform between the rigid and elastic blades are almost the same. To clearly observe the effect of the blade aeroelasticity in the wake region, the velocity difference along the wind direction between the deformable and the rigid blades of the fully coupled FOWT is shown in Figure 12. Figure 12(a and f) demonstrate that the influence of the elastic blade on the velocity in the wake field is minor when the wave is insignificant. However, when the above-mentioned coupled effect of the FOWT motion and the high wave elevation happened near the focused time (Figure 12(b–e)), there is a significant difference near the free surface in the wake flow field after considering the deformation of the blade. It is caused by the blade flap motion, which participates in the coupled effect, and the local velocity is accelerated or delayed in the flow field. And such a slight variation near the turbine blade tip will be exacerbated by the wake propagation along the incoming wind direction.

Figure 12. The velocity field difference of the fully coupled wind turbine at midplane between the rigid and elastic blade: (a) t = 100 s; (b) t = 144 s; (c) t = 150 s; (d) t = 160 s; (e) t = 170 s; (f) t = 300 s.

To study the deficit of velocity in the wake region, the axial velocity distribution at 100, 150, and 200 m behind the rotor is investigated in Figure 13, in which the longitudinal coordinate origin (x = 0) represents the initial location of the wind turbine. Each curve is the local velocity distribution profile and represents the percentage of the unit value, which is labelled in the figure. The shapes of the velocity distribution at t = 100 s and t = 300 s behind the rotor are almost the same and partly cover each other. However, the velocity pattern is much different near the focus time. The velocity acceleration region about 150 m to 200 m behind the wind turbine and closed to the free surface is significantly affected by the blade aeroelasticity. There is an underestimation of the velocity recovery in this region when the deformation of the wind turbine blade is neglected. It disturbs the axial velocity distribution of the wake field and greatly changes the shape of the wake contour.

Figure 13. Comparisons of the velocity distribution in the fully coupled wind turbine wake between the rigid and elastic blade: (a) t = 100 s; (b) t = 144 s; (c) t = 150 s; (d) t = 160 s; (e) t = 170 s; (f) t = 300 s.

The shape of the wake region can also be observed in the horizontal section, as described in Figure 14. The wake fields at x = 200 m behind the turbine rotor and platform are compared between the elastic and rigid blades. The original circle of the fixed wind turbine wake is changed into an oval shape when the FOWT motion is considered (referred to Figure 14(a and f)). However, the wake patterns are similar between elastic and rigid blades without the influence of the FOWT motion and the high wave elevation. Around the focused time, such a coupled effect disturbs the wake field and generates a raised velocity region. In this region, the elastic blade has a more serious influence and will be further analyzed next. Consequently, the velocity profiles and distribution spanwise analyzed above demonstrated that the accuracy of the rigid model is not enough for engineering applications when focused wave events are encountered.

Figure 14. Comparisons of the velocity field at x = 200 m in the fully coupled wind turbine wake between the elastic and rigid blade: (a) t = 100 s; (b) t = 144 s; (c) t = 150 s; (d) t = 160 s; (e) t = 170 s; (f) t = 300 s.

To quantify the phenomena discussed above, the flow velocities at several positions behind the rotor are sampled, as shown in Figure 15. On this basis, the time histories of these wake field velocities are plotted in Figure 16. From Figure 16(a, c, and e), the wind speed is steady due to the absence of the free surface in the coupled effect. The fluctuation (point A’) is mainly caused by the motion of FOWT. The distance between adjacent points (50 m) divided by the time corresponding to the velocity peak is exactly equal to the local wind speeds (approximately 8.44, 7.96, and 7.70 m/s, respectively). It indicates that the propagation velocity in the wake fields is directly proportional to the local wind speed.

Figure 15. The velocity sampling location of the wind turbine wake field.

Figure 16. Comparisons of the wind turbine wake field velocity time history: (a) point A; (b) point B; (c) point C; (d) point D; (e) point E; (f) point F.

The slight variety in the wake will be amplified since the large free surface and FOWT high-frequency motions (see Figure 16(b, d, and f), and Table 7). The maximum impact ratios of the blade aeroelasticity on the velocity peak values are approximately −8.63%, 18.60%, and 5.17%, respectively. It proves that the influence of the blade deformation on the wake velocity is related to the locations in the wake field. It will reduce the local wind speeds in the area within 1D (D is the diameter of the rotor) behind the wind turbine. Nevertheless, this effect accelerates the local velocities significantly beyond the incoming wind speeds near 1D behind the rotor. As the distance increases, this influence will also propagate backward. Although the influence gradually decreases, this non-uniformity in the flow field has a great impact on the aerodynamic loads of the downstream wind turbine blades and even increases their fatigue damage.

Table 7. Comparisons of wake field velocity ratio between elastic and rigid blades.

Point	Time	Elastic blade	Rigid blade	Error
B	B’	0.973	1.057	−8.63%
	C’	1.137	1.099	3.34%
	D’	0.934	1.012	−8.35%
D	B’	1.403	1.142	18.60%
	C’	1.027	1.041	−1.36%
	D’	0.931	0.961	−3.22%
F	B’	1.316	1.248	5.17%
	C’	0.942	0.944	−0.21%
	D’	0.902	0.927	−2.77%

5.4. Wake characteristics of wind turbine

The blades’ elasticity will change their position and the corresponding velocity in the flow field compared with the rigid blade. Especially for the blade tip, its deformation and vibration are obviously larger than those of other parts of the turbine blade. Due to this, the vortex shedding from the blade tip will also be affected, as discussed below. Figure 17 shows the wake structure of FOWT with the elastic blade, which is visualised by using the Q criterion. The large wave has not formed in Figure 17(a), and the FOWT has been travelling steadily after the focused time in Figure 17(f). These two conditions can be regarded as only FOWT motion. According to Figure 17(a and f), only vortex shedding from the turbine hub can propagate up to the outlet boundary of the N-S domain, whereas the vortex shedding from the blade tip cannot. However, when the focused wave happens, the vortex shedding from the blade tip can also propagate up to the outlet boundary of the N-S domain. As a result, because of the initial location and velocity of the vortex shedding from the blade tip vary, the extruded vortex by the free surface spreads farther away than that generated by only FOWT motion. It will affect the production of the downstream wind turbine and increase its’ blade fatigue damage.

Figure 17. Wake structures of the fully coupled wind turbine with elastic blade: (a) t = 100 s; (b) t = 144 s; (c) t = 150 s; (d) t = 160 s; (e) t = 170 s; (f) t = 300 s.

6. Conclusion

In this study, the investigation of the OC4 Semi-submersible NREL 5 MW FOWT subjected to focused wave and uniform wind is carried out. The modelling was achieved by using the hybrid model, which couples the aerodynamics of a wind turbine, the blade structural dynamics, the hydrodynamics of a floating platform, and mooring dynamics. As this hybrid model has been evidenced, the fully aero-elastic-hydro-mooring-wake dynamics interaction of FOWT can be well captured. The uncertainty and the convergence test of grids are analyzed in this paper, which can demonstrate that the results are credible.

The influence of blade aeroelasticity on the FOWT’s dynamic behaviour is discussed. For the turbine aerodynamics, it increases fatigue damage but does not affect the trend of the loading on the turbine blades. In the present physical condition, the deformation of flexible blades influences the thrust and power by less than 5%. This conclusion is consistent with the previous research with either the CFD-MBD method (Liu et al., 2019) or the CFD-CSD method (Zheng et al., 2023). As a result, this effect is further weakened by the turbine tower and cause the platform hydrodynamics are insensitive to the deformable blades.

For the wake flow fields behind the turbine rotor, the large wave elevation, the motion of FOWT, and the blade elasticity are coupled together. They transform the flow field velocity and the wake structure. Especially during the focused time, the coupled effect will be obvious and further disturb the wake flow. The slight variety in the wake will be amplified by about 0.21% to 18.6%. Based on the temporal and spatial distribution of the velocity fields, the influence of the wind turbine blade deformations on the wake flow field is mainly concentrated in the area at a height of 10 m to 40 m from the horizontal plane, approximately 1D behind the wind turbine. Although this difference in quite small extent of time and space will have a quite small effect on the operation of the wind turbine, it may be amplified as the turbine wake propagates backwards and has a great impact on the downstream wind turbine. In practice, the FOWT system is primarily subjected to random wind and waves. The occurrence of the transient focusing wave occurs frequently during the FOWT’s operating life. It leads to the above impact in the wake flow region repeating multiple times and influencing the downstream FOWT characteristic, which will be explored in the future.

Acknowledgement

This research work was funded by the National Key Research and Development Program of China (No.2020YFB1506701), the National Natural Science Foundation of China (Nos. 51739001; 51879051); Natural Science Foundation of Heilongjiang Province in China (LH2020E071), Open Fund of Zhejiang Provincial Key Laboratory of Wind Power Technology (ZOE2020007). Ziying Yu thanks Qingwei Ma and Shiqiang Yan at City, University of London for the support in software qaleFoam during his study.

Data availability

The data that supports the findings of this study are available within the article.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This research work was funded by the National Key Research and Development Program of China (No. 2020YFB1506701), the National Natural Science Foundation of China (Nos. 51739001; 51879051); Natural Science Foundation of Heilongjiang Province in China (LH2020E071), Open Fund of Zhejiang Provincial Key Laboratory of Wind Power Technology (ZOE2020007).