Numerical Study of the Aerodynamic Loads on Offshore Wind Turbines under Typhoon with Full Wind Direction

Abstract

Super typhoon activity is likely to make the electric power network fail, or blow the wind-measuring device off, which all lead to the yaw control system of wind turbine being inactive. Under this condition, blades can be blown by the violent side wind from unfavorable directions, and the aerodynamic loads on the wind turbine will be increased by a large amount, which can lead to low-cycle fatigue damage and other catastrophic collapses. So far, not enough consideration has been given to the above problems in wind turbine design. Using the transient computational fluid dynamics (CFD), this study investigates the wind load characteristics of offshore wind turbines under typhoon condition with 360-degree full wind directions. Two primary influence factors of the aerodynamic characteristics of wind turbines are clarified: variation of the wind direction and different parking positions of the wind rotor. Using 3D-numerical simulation results, this study provides detailed references for the ultimate strength and fatigue life in wind turbine design, and presents the best parking position of the wind turbine with a free yawing strategy.

1. Introduction

Wind supplies competitive renewable energy for mankind’s economic development, especially in the face of the issues of global warming and the shortage of energy supplies. The recent technology trend in the wind energy industry has been towards offshore wind power projects [1]. Offshore wind farms have distinct advantages over onshore wind farms, as well as various new situations and challenges, among which the typhoon has been one of the unavoidable problems [2]. For the past few years, wind farms in China have suffered huge economic losses due to typhoon attacks. In 2003, typhoon Dujuan struck the Honghaiwan wind farm, nine wind turbine blades were irreparably damaged and the yaw control systems of six wind turbines were destroyed, even though the wind speeds in this typhoon activity were much less than the design extreme wind speed [3]. In 2006, five wind turbine towers collapsed and 33 blades were broken when the super typhoon Saomai struck the Hedingshan wind farm, which resulted in the loss of 70 million RMB [4]. In 2010, six wind turbines in the Liuao wind farm were destroyed by typhoon Megi in Fujian Province of China [5]. In 2013, the severe tropical storm Usagi caused 11 blades to break off near their roots or at the roots, excluding the blades of the eight collapsed towers [6]. Therefore, it is necessary to enhance the wind resistance designs and optimize the control strategy for those wind farms that lie in typhoon activity zones [7,8].

In emergency conditions such as a typhoon, wind turbines can automatically start up their brake devices due to the violent wind speed. Under the braking-stop state, blades remain in the feathering position with a 90° pitch angle. During a super typhoon activity, the gale and rainstorm are likely to make the electric power network fail, or blow the wind-measuring device (including wind vanes and anemometers) off, which all lead to the yaw control system being inactive [9]. Under this condition, the wind rotor position can no longer be automatically adjusted to face the favorable wind direction by the yaw system, and violent wind can blow over the wind turbine along unfavorable directions, which will cause rapid increased wind loads to act on the blades, and then, the wind turbine structures are prone to being damaged. However, there are almost no definite stipulations in design standards, nor are there precise descriptions or strength analysis methods of this complex loading condition at present. In this case, protecting blades from being damaged depends heavily on the material, component bond type, and production technology of the manufacturers [10].

In this work, 3D simulations of the FX93-2500 offshore wind turbine were performed based on the wind turbine being in the braking-stop state and its yaw control system being inactive under typhoon conditions, and the wind direction was varied in the interval (−180°, 180°) to analyze the aerodynamic characteristics of wind turbines under typhoon with 360° full wind directions. The aerodynamic characteristics of blades are sensitive to wind direction changes because of the special constructional features of blades. When a typhoon blows in the direction shifting from the feathering position, the vortex-induced vibration caused by boundary layer separation of blade becomes one of the major reasons for wind turbine failures. Vortices caused by boundary layer separation often induce large fluctuating pressures on the rotor that may lead to blades, towers and foundations being subject to accumulated damage by low-cyclic fatigue or resonance. To accurately simulate the formation and development of the separation vortex, the mesh of the computational domain must be fine enough to capture the geometric detail of blades, because the vortices caused by boundary layer separation are directly related to the blade’s geometrical shape. It was finally determined that the total number of the structured mesh elements in the flow field computational domain of the 2.5 MW large-scale offshore wind turbine was almost 14 million. When considering the turbulence intensity of typhoon, it is necessary to utilize large-eddy simulation (LES) with a pulsating velocity inflow boundary condition [11,12,13,14]. Together with the enormous quantity of mesh cells, this LES method is very difficult to realize because of the high computational resources required and the length of time involved. Thus, as an approximation to the actual situation, the average wind speed variation with altitude was applied in the inflow boundary, and this paper is focused on the analysis of aerodynamic loads caused by blade boundary layer separation under typhoon with 360° full wind directions, aiming to provide detailed references for the ultimate strength and fatigue life in wind turbine design, and attempting to find the best parking position of the wind turbine with a free yawing strategy.

2. Numerical Method

2.1. Physical Model

The rotor blade considered in this work is the FX93-2500 wind turbine, which is used in practical engineering in the coastal area of Southeast China. The FX93-2500 wind turbine is a horizontal-axis, pre-bent three-blade, upwind and atmospheric temperature grid-connected wind turbine generator system. The main parameters of the FX93-2500 wind turbine are given in

Table 1. The configuration of the FX93-2500 wind turbine.

Parameter	Value	Parameter	Value
Rated power (MW)	2.5	Length of the blade (m)	45.3
Power modulation	variable pitch and speed	Cut-in wind speed (m·s−1)	3
		Cut-out wind speed (m·s−1)	25
Number of blades	3	Rated wind speed (m·s−1)	12
Rotor diameter (m)	93.4	Rated turning speed (rpm)	17
Hub height (m)	86	Useful life (year)	20

.

With the rapid development of the computer performance, the computational fluid dynamics (CFD) method, which, based on the Navier-Stokes equation, is more useful to utilize as an efficient tool for the wind turbine design and analysis [15]. Comparing with the traditional BEM method [16,17], the CFD simulation possesses high accuracy and can get the details of the blade surface pressure distribution and the flow field information, resulting in a better understanding of the aerodynamic phenomena on the rotor flow field [18,19,20,21,22,23,24].

The wind rotor is the key component of the wind turbine, of which the aerodynamic performance directly determines the wind energy extracted from the wind by the wind turbine and the main aerodynamic loads of the whole wind turbine system. The aerodynamic performance of the wind rotor is the core content of wind power technology [25]. Therefore, in order to reduce the complexity of the problem and the time spent on computational analysis, this study ignores the influence of the nacelle and tower, given that they are not the focus of the research [19,26].

During a typhoon activity, wind turbines are parked in the feathered position. To facilitate the analysis and calculation, the rotor coordinate system {X, Y, Z} and blade coordinate system {X_B, Y_B, Z_B} were established, which can be seen from Figure 1. The origin of the rotor coordinate system is located in the center of the wind rotor, the X-axis is the rotation axis of the wind rotor, the Y-axis lies on the horizontal plane and is perpendicular to the X-axis, the Z-axis is sticking straight up, and YOZ is the rotating plane of the wind rotor. As shown in Figure 1c, the origin of the blade coordinate system is located in the root center of the corresponding blade, the X_B-axis is parallel to the X-axis, the Z_B-axis is along the blade span, and the Y_B-axis is approximately perpendicular to the airfoil chord-wise direction as the blade is in the feathering position [27]. Figure 1 also shows the parking position of the wind rotor during a typhoon. The three blades of the wind rotor are indicated as B1, B2 and B3, the angle between blade B3 and the Z-axis is θ, and the wind direction is measured by the inflow angle β. In this study, we analyzed three parking positions of wind rotor, about which the values of θ were 0°, 30° and 60°. Figure 1d shows that the wind direction β varies in the internal (−180°, 180°), among which β = 0° represents the wind blowing along the positive direction of the X-axis, β = 90° represents the wind blowing along the positive direction of the Y-axis, and β = −90° represents the wind blowing along the negative direction of the Y-axis.

2.2. Computational Domains and Mesh Generation

The scale of the computational domain for this study is shown in Figure 2a based on the wind rotor diameter D. The distances from the wind turbine to the boundaries all around were 5D, the distances from the wind turbine to the upper and lower boundary were 3D and 86 m (hub height) respectively.

To obtain better convergence and precise aerodynamic results, full-hexahedral mesh grids of reasonable quality subdivided the entire computational domain by applying the block topology, which can be seen in Figure 2. By adopting the O-grid generation technique [28], mesh refinement was generated near the blade surface to describe with sufficient precision the boundary layer flow, and the size of the first element in the wall-normal direction was fine enough to ensure y⁺ < 5. This y⁺ value made it suitable for the use of the enhanced wall function [29]. As a result of the thin twist blades with their highly complex shape of the FX93-2500 wind turbine, the difficulty and quantity of the structured mesh generation are greatly increased. The fluid computational domain contained, on average, 13,953,703 hexahedral mesh elements, which was several times larger than the previous wind farm CFD studies [1,18,19,22,30,31].

2.3. Numerical Simulation Method Based on CFD

There is no definite stipulation regarding the typhoon wind profile in existing design standards. Powell has indicated that the average wind speed in the boundary layer of a typhoon within the near-surface range of 200 m changes with the altitude and the change agrees with the power-law distribution [32]. Based on the results of typhoon wind-field full-scale measurements in China coastal areas, the Li Shengqiu research group suggested a set of reference designators regarding the typhoon wind field, in that the wind profile approximately obeyed the power-law distribution at altitudes less than 200 m. However, the fitting power-law exponent was different from the normal wind condition, which was suggested to be α = 0.142 [33,34].

In this study, the wind turbine flow fields were simulated based on the three-dimensional Navier-Stocks equation, applying a high order accurate solver method. The RNG k-ε turbulence model was adopted, which is widely used and has more data accumulation and considerable precision [29,35,36]. Additionally, the enhanced wall function was utilized to obtain the flow field characteristics in the wall region of the turbulent boundary layer. Therefore, the time-dependent aerodynamic loads on the surface of wind rotor were computed by solving the unsteady Reynolds-averaged Navier-Stocks equation [37]. According to suggestions given by the Li Shengqiu research group, the typhoon wind profile applied to the inflow boundary can be expressed as:

$$V=V_{\text{ref}}{\left(\frac{z}{z_{\text{ref}}}\right)}^{\alpha }$$

(1)

where V is the wind speed according to the height; V_ref is the wind speed at the reference height, the maximum average wind speed at the height of 10 m is 30.3 m/s based on the measured results of typhoon Chanthu [34]; z is the height, z_ref is the reference height, α is the fitting power-law exponent, which is 0.142 in this work.

Based on the above, the time step was chosen to be Δt = 5 × 10⁻⁴ s, and the aerodynamic characteristics of the wind turbine in the typhoon were simulated with the wind direction β varying in the range (−180°, 180°). The computations were carried out in a parallel computing environment, using two high-performance servers, taking approximately 20 months to complete simulations for all of the cases.

3. Analysis of Flow Field Characteristics

When the wind direction β = 0°, the violent wind of typhoon was blowing along the favorable direction with the blade angle of attack at zero degree; the air flowed around the blade keeping the streamlines running smoothly, and there was no boundary layer separation. With the increase of the wind direction β, the blade angle of attack also grew, and then a vortex started to appear on the leeward side of the blade, which was the boundary layer separation. As the wind direction β continued to increase, the eddy zone expanded rapidly, and the flow field situation continued to grow worse. Figure 3 shows the blade positions when the wind rotor is in the braking-stop state. It is obvious that blades at different positions have different angles of attack.

The inflow velocity is V, and it can be decomposed in the X-axis direction and the Y-axis direction:

$$V_X=V\cos \mathsf{\beta }, V_Y=V\sin \mathsf{\beta }$$

(2)

Supposing that there is no interaction among the three blades of the wind rotor, and ignoring the twist of the blade, Figure 3b shows the wind speed triangle and aerodynamic forces on the airfoil section when the airflow gets around the blade in the feathering position, among which, airfoil chord parallels with the X_B-axis, F_L and F_D are the aerodynamic lift and drag forces. In the plane X_BO_BY_B, V_X and V_Y can also be decomposed in the X_B-axis direction and Y_B-axis direction:

$$V_{r{X}_B}=V_X,V_{r{Y}_{\mathrm{B}}}=V_Y\sin \varphi$$

(3)

where $\varphi$ is the angle between the blade and the positive direction of the Y-axis, which increases in the clockwise direction as shown in Figure 3a; V_rXB is the projection of V_X and V_Y onto the X_B-axis direction, V_rYB is the projection of V_X and V_Y onto the Y_B-axis direction.

As shown in Figure 3b, V_r is the resultant velocity of V_rXB and V_rYB, and the angle between V_r and the X_B-axis is actually the blade angle of attack $\phi$, which can be expressed as:

$$\phi =\arctan \frac{V_{r{Y}_{\mathrm{B}}}}{V_{r{X}_{\mathrm{B}}}}=\arctan (\frac{V_Y}{V_X}\sin \varphi )=\arctan (\tan \mathsf{\beta }\sin \varphi )$$

(4)

It shows that the angle of attack reaches maximum on the blade in the vertical position ($\varphi$ = 90° or 270°), and the attack angle of the blade in the horizontal position ($\varphi$ = 0° or 180°) is approximately equal to zero. Moreover, under the favorable wind direction β = 0°, the attack angles of all blades are approximately equal to zero; for the side wind direction β = ±90°, the attack angles of all blades are also ±90°.

Figure 4 shows the streamline distribution on the leeward side of blade, from which we can see that the boundary layer separation conditions agree with Equation (4). The blade B3 of the θ = 0° parking position and blade B1 of the θ = 60° parking position were vertical, and they had greater attack angles than other blades, therefore, the boundary layer separation appeared first with the increase of the wind direction β. Meanwhile, blade B2 of the θ = 30° parking position was located in the favorable horizontal position, and there was almost no boundary layer separation in any wind direction. Moreover, when the attack angle $\phi$ < 0°, the upper surface of the blade was facing the wind, and then, the flow state tended to get worse than that the lower surface against the wind ($\phi$ > 0°).

There exist obvious interactions among the three blades of the wind rotor, which can be seen in Figure 5. A vortex is a circular or spiral set of streamlines, and a vortex core is a special type of isosurface that displays a vortex. Here we apply the λ₂-criterion to determine the vortex core structure in the flow field. As shown in Figure 5, for the same wind rotor, it is clear that the downwind blades have a contraction effect on the trailing vortex of the upwind blades. For the wind rotor in the θ = 0° parking position, the interaction between blades B1 and B2 was conspicuous: the downwind blade narrowed the range of the vortex region of the upwind blade, ameliorating the condition of the flow field around the upwind blade, while the flow state around the downwind blade got worse. Blades B2 and B3 of the θ = 60° parking position is similar to that described above. In particular, for the wind rotor in the θ = 30° parking position, blade B2, which was in the horizontal position, significantly affected both blades B1 and B3: for the wind direction β < 0°, blade B2 posed a strong contraction effect on the trailing vortex of the other two blades, and then the flow state around blades B1 and B3 was significantly improved.

4. Analysis of Aerodynamic Characteristics

The variable wind directions under typhoon activity may lead to large areas of blades being affected by the violent wind, which can cause large increments of wind loads to act on the blades. The strong wind loads on blades can bring huge impacts on the regulating mechanism and the yaw control system, resulting in failures of wind turbine structures. Additionally, in the case of large angle of attack, the wind loads on blades are evidently unsteady, which will have an important influence on the fatigue of the wind turbine.

4.1. Total Aerodynamic Forces Acting on the Wind Rotor

The aerodynamic force is one of the major loads for the wind turbine. The total aerodynamic force acting on the wind rotor can be passed to the tower and foundation through the transmission shaft, which will significantly affect the loading on the tower and foundation. Figure 6 shows the relationship between the total aerodynamic forces of wind rotor and the wind directions, among which F_X, F_Y and F_Z are the three-dimensional total forces in the wind rotor coordinate system, and the plus and minus signs here just indicate the loading direction.

Figure 6a provides the average values of the total aerodynamic forces as the wind direction β varies in the interval (−180°, 180°). It is clear that the horizontal force F_Y is the greatest, followed by the axial force F_X and the vertical force F_Z. Overall, the axial force F_X and the horizontal force F_Y approximately experienced a process of increasing first and then decreasing, with sharp variance in the range of |β| = 0°–30°. The maximum of F_X occurred in the range |β| = 75°–105°, and the maximum of F_Y occurred in the range |β| = 30°–45°. However, the vertical force F_Z fluctuated with wind direction, and F_Z on the wind rotors in the θ = 0° and θ = 60° parking positions had a similar magnitude but acted in the opposite direction. It is worth mentioning that the average values of F_Y and F_Z are almost equal to zero at the favorable wind direction β = 0°; in contrast, the axial force F_X is greater at this time, and acts in the positive direction only when β = −5°–5° because of the special aerodynamic form of the blade.

Figure 6b gives the root-mean-square (RMS) amplitude of the total aerodynamic forces, in which the RMS amplitude is the root mean square of the fluctuating load that equals the aerodynamic load minus its average value. It is shown that the RMS amplitudes of the total aerodynamic forces increase quickly when the wind direction |β| rises above 30°, which agrees with the flow field characteristics shown in Figure 4. A large vortex started to appear on the leeward side of the blade when the wind direction rose to |β| ≥ 30°, as a result, the fluctuating load rapidly increased. Moreover, Figure 6b also shows that the peak RMS amplitudes of the total aerodynamic forces occurred in the range of |β| = 105°–150°, in these cases, the fluctuating load can be much greater than that produced by the turbulence of typhoon [34,38,39,40], and it may have a large effect on the fatigue properties of the tower and foundation, potentially causing devastating damage to them.

Furthermore, for wind rotors in different parking positions, there was little difference in the average values of the total aerodynamic forces except F_Z, but a great difference in the fluctuation characteristic of the total aerodynamic force. The RMS amplitudes of F_X and F_Y acting on the wind rotor in the θ = 30° parking position were significantly smaller than that of F_X and F_Y acting on the wind rotors in the other two parking positions, but the RMS amplitude of the vertical force F_Z on the wind rotor in the θ = 30° parking position was larger. The RMS amplitude of F_X reached its maximum value of 9.89 kN at the θ = 0° parking position, which was 2.04 and 1.20 times that of the θ = 30° and θ = 60° parking positions; the RMS amplitude of F_Y reached its maximum value of 9.80 kN also at the θ = 0° parking position, which was 2.27 and 1.07 times that of the θ = 30° and θ = 60° parking positions; the RMS amplitude of F_Z reached its maximum value of 4.28 kN at the θ = 30° parking position, which was 2.77 and 2.94 times that of the θ = 0° and θ = 60° parking positions. Figure 7, Figure 8 and Figure 9 provide the time-domain waveforms and the frequency spectrograms of the total aerodynamic forces acting on the wind rotor at several typical wind directions. Based on the figures, under the same wind flow direction, for the three parking positions of the wind rotor, there is obviously difference in the magnitude and dominant frequency of the fluctuating load in the X-axis direction, from which we can observe that the parking position significantly influence the fluctuation characteristic of the total aerodynamic force.

Table 2. Range of the dominant frequency of the total aerodynamic force acting on the wind rotor.

Dominant Frequency (Hz)		FX	FY	FZ
θ = 0°	β > 0°	0.080–14.960	0.080–14.960	0.080–14.960
	β < 0°	0.165–13.740	0.178–13.740	0.165–13.740
θ = 30°	β > 0°	0.127–6.250	0.127–6.250	0.085–14.480
	β < 0°	0.289–5.202	0.289–6.650	0.180–6.650
θ = 60°	β > 0°	0.055–6.480	0.055–6.480	0.055–6.480
	β < 0°	0.108–13.480	0.108–13.480	0.125–13.480

gives the range of the dominant frequency of the total aerodynamic force; the result data indicates that the dominant frequency of the total aerodynamic force is not only related to wind direction but also rotor parking position. Those feathering wind directions of −5°–5° or 180° are ignored here to observe clearly the effects of wind direction and rotor parking position on the wind load dominant frequency. It shows that the ranges of the dominant frequency of the three-dimensional total forces acting on the wind rotor are nearly the same size, and the dominant frequency of the total aerodynamic force changes within a large scale from low frequency to high frequency. Most forces with high frequency appeared in the range of |β| = 15°–30° or 150°–165°, and forces with low frequency appeared in the range of |β| = 60°–120°. For the wind rotor in the θ = 0° parking position, the range of the dominant frequency was larger than that of the other two parking positions; for the wind rotor in the θ = 60° parking position, the dominant frequency of the total force was more significantly affected by wind direction, and there existed no dominant frequency greater than 7 Hz when the wind direction β > 0°; for the wind rotor in the θ = 30° parking position, the range of the dominant frequency was the minimum, except that of the dominant frequency of F_Z.

The coordinate system {X_G, Y_G, Z_G} was built in the center of the foundation, whose coordinate axes were the same as the rotor coordinate system {X, Y, Z}; based on that, we could calculate the aerodynamic loads acting on the foundation passed from the wind rotor. As shown in Figure 10, M_ZG is the torsional moment passed from the wind rotor, and M_XYG is the overturning moment passed from the wind rotor, which can be expressed as:

$$M_{XY\text{G}}=\sqrt{M_{X_G}^2+M_{Y_G}^2}$$

(5)

Similar to the total aerodynamic force, the average values of M_XYG and M_ZG also increased first and then decreased gradually, reached their maximum in the range of |β| = 30°–45°. For wind rotors in different parking positions, the averages of M_XYG were similar, the maxima of which were 28.52, 28.11 and 27.01 MN·m in the θ = 0°, θ = 30° and θ = 60° parking positions, respectively; however, the averages of the torsional moment M_ZG had great changes with the parking positions, and the curve of M_ZG varied with wind direction in the θ = 30° parking position was similar to the curve which was obtained by translating the curve of M_ZG in the θ = 0° or θ = 60° parking position 1 MN·m upwards. Then, the torsional moment M_ZG of the θ = 30° parking position was greater than that of the other two parking positions when the wind direction β < 0°, and it became smaller when the wind direction β > 0°.

4.2. Total Aerodynamic Moments of the Wind Rotor

In the rotor coordinate system {X, Y, Z}, the total aerodynamic moments of the wind rotor can be divided into the wind torque M_X, pitching moment M_Y and yawing moment M_Z. Figure 11, Figure 12, Figure 13, Figure 14, Figure 15 and Figure 16 show the curves of variation of the average total aerodynamic moments with wind direction, and Figure 17 shows the fluctuation characteristic of the total aerodynamic moment. Unlike the total aerodynamic force, the total aerodynamic moments are significantly affected by the parking positions, in both average values and fluctuation characteristics.

Figure 11 shows the curves of the average wind torque M_X varied with wind direction. To show the variation clearly, we divided the wind direction interval into two sections. In the wind direction interval (−90°, 90°), the average wind torque M_X increased first and then decreased approximately, reaching its maximum when the wind direction |β| = 30°–45°; the wind torque M_X of wind rotors in the θ = 0° and θ = 60° parking positions were opposite for most wind directions, and they were all greater than that of the rotor in the θ = 30° parking position. In the wind direction intervals (−180°, −90°) and (90°, 180°), there were little differences in the average wind torque M_X among the three parking positions, and the peak values appeared at |β| = 165°–180°. In full wind direction, the average wind torque M_X of the rotor in the θ = 0° parking position reached the maximum value of 1341.66 kN·m at the wind direction of β = 30°; the average wind torque M_X of the rotor in the θ = 30° parking position reached the maximum value of 1240.91 kN·m at the wind direction of β = ±180°; and the average wind torque M_X of the rotor in the θ = 60° parking position reached the maximum value of 2279.74 kN·m at the wind direction of β = −30°, which was 1.37 times that of the rated operating state.

To better understand the variation of the wind torque M_X with wind direction, Figure 12 shows the contribution of each blade to the wind torque. We use $M_X^b$, $M_Y^b$ and $M_Z^b$ here to represent the moments of a single blade in the rotor coordinate system, in order to distinguish them from the total aerodynamic moments of the wind rotor. It is clear that the moment $M_X^b$ of each blade is subject to a sinusoidal distribution, which directly relates to the attack angle of blade as shown in Equation (4). The plus-minus of the attack angle determines the direction of the blade moment $M_X^b$, and the size of the attack angle also approximately determines the size of the blade moment $M_X^b$. With little effect of the hub, the wind torque M_X is approximately equal to the sum of the moments $M_X^b$ of the three blades of the wind rotor. Therefore, the wind torques of the rotors in the θ = 0° and θ = 60° parking positions were opposite for most wind directions because of the different windward sides (opposite algebraic signs of the attack angle); for the wind rotor in the θ = 30° parking position, the moments$M_X^b$ of blades B1 and B3 were nearly the same in size and opposite in direction, thus making the wind torque smallest. Moreover, as the moments $M_X^b$ of the three blades of the wind rotor had the same direction when the wind direction |β| = 150°–180°, a peak of the wind torque appeared at the wind direction angle of ±180° nearby.

Figure 13 shows the curves of the average pitching moment M_Y varied with wind direction. The average pitching moment M_Y shows approximately a linear relationship with the wind direction β in the two sections of the wind direction interval. In the range of small wind direction angles (−15° < β < 15° and 150° < |β| < 180°), the pitching moments M_Y of rotors in the θ = 0° and θ = 60° parking positions are very sensitive to the changes in wind direction, which dropped to zero quietly at the favorable wind direction (β = 0° or ±180°). The influence of rotor parking positions on the pitching moment M_Y was also obvious: the pitching moment M_Y of the rotor in the θ = 0° parking position was almost in the negative direction, the pitching moment M_Y of the rotor in the θ = 60° parking position was almost in the positive direction, and the pitching moment M_Y of the rotor in the θ = 30° parking position fell in between the above two. In full wind direction, the average pitching moment M_Y of the rotor in the θ = 30° parking position reached the maximum value of 795.95 kN·m at the wind direction of β = 90°; the average pitching moment M_Y of the rotor in the θ = 60° parking position reached the maximum value of 1388.66 kN·m at the wind direction of β = 105°; the average pitching moment M_Y of the rotor in the θ = 0° parking position reached the maximum value of 1587.76 kN·m at the wind direction of β = −105°, which was 7.63 times that of the rated operating state, even close to the value of the wind torque in the rated operating state, that would pose a great threat to the regulating mechanism.

Figure 14 shows the contribution of each blade to the pitching moment. Unlike operating states, the pitching moment has almost nothing to do with the wind shear, but directly relates to the blade position. The pitching moment of the rotor in the θ = 0° parking position was dominated by blade B3, the pitching moment of the rotor in the θ = 60° parking position was dominated by blade B1, and the moments $M_Y^b$ of the above two blades were opposite because of the different windward sides, so the pitching moments of the rotors in these two parking positions were also opposite; moreover, the moment $M_Y^b$ of each blade changed sharply in the range of small attack angles due to the aerofoil lift-drag force characteristic, and then the pitching moment became very sensitive to the wind direction changes in the range of small wind direction angles. Alternately, the pitching moment of the rotor in the θ = 30° parking position was dominated by blades B1 and B3 jointly, among which the moment $M_Y^b$ of blade B1 was similar to that of blade B1 of the rotor in the θ = 60° parking position, and the moment $M_Y^b$ of blade B3 was similar to that of blade B3 of the rotor in the θ = 0° parking position, it was the reason why the pitching moment of the rotor in the θ = 30° parking position fall in between that of the other two parking positions, with a smooth variation and no large gradient.

Figure 15 provides the curves of the average yawing moment M_Z varied with wind direction. Different from the wind torque or pitching moment, the yawing moment M_Z of the rotor in the θ = 30° parking position is much greater than that of the other two parking positions and surpasses 600 kN·m for almost all wind directions, reaching its maximum value of 1189.68 kN·m at the wind direction of β = 90°. This maximum was 2.38 and 2.24 times that of the other two parking positions, and reached approximately 53 times that of the rated operating state, which would make a large impact on the yaw control system, possibly resulting in damage. Moreover, in the range of small wind direction angles (−30° < β < 30° and 135°< |β| < 180°), the yawing moment M_Z of the rotor in the θ = 30° parking position was very sensitive to the wind direction changes, which dropped to zero quietly at the favorable wind direction (β = 0° or ±180°). The variation of the yawing moment can be further described in Figure 16, which gives the contribution of each blade to the yawing moment. For wind rotors in the θ = 0° and θ = 60° parking positions, the moment $M_Z^b$ of the blade in the vertical position was approximately zero, and the moments $M_Z^b$ of the other two blades of the same rotor were similarly equal and opposite, therefore, the yawing moments of the rotors in these two parking positions were very small. In contrast, for the wind rotor in the θ = 30° parking position, the moments $M_Z^b$ of the three blades of the same rotor were almost in the same direction, resulting in the enormous yawing moment.

Based on the above characteristics of the yawing moment, it can be worthwhile to adopt the typhoon resistance strategy without electrical power, in order to avoid the yaw control system becoming inactive when the electric power network fails or the wind-measuring device breaks down under typhoon activity, i.e., starting the typhoon resistance strategy after the forecast of a typhoon, setting the blades in the feathering position, and loosening the yaw brake device to make the wind rotor act as a tail vane to achieve the free yawing. Comparing the different parking positions, the yawing moment of the rotor in the θ = 30° parking position was the greatest and very sensitive to wind direction changes, the intensity of which was enough to conquer the friction moment at the top of the tower, assuring the free yawing of the wind turbine with wind direction changes. Accordingly, these results suggest turning the wind rotor to the θ = 30° parking position, loosening the yaw brake, and the wind rotor will finally be adjusted to the favorable wind direction conditions (β = 0° or ±180°) under the function of the great yawing moment M_Z. At this moment, it can significantly reduce the typhoon wind loadings on the wind rotor, and protect the safety of wind turbine structures.

Figure 17 provides the RMS amplitude of the total aerodynamic moment. Similar to the total aerodynamic force, the peak RMS amplitudes of the wind torque and pitching moment also occurred in the range of |β| = 105°–150°, but the peak RMS amplitudes of the yawing moment occurred in the ranges of β = 75°–135° and −90°–−120°. Of those, the maximum of the RMS amplitude of the wind torque appeared in the θ = 0° parking position, which was 111.16 kN·m at the wind direction of β = −150°, 30.89% and 6.01% greater than that of the θ = 30° and θ = 60° parking positions, respectively; the maximum of the RMS amplitude of the pitching moment also appeared in the θ = 0° parking position, which was 148.29 kN·m at the wind direction of β = 105°, 82.17% and 18.57% greater than that of the θ = 30° and θ = 60° parking positions, respectively; the maximum of the RMS amplitude of the yawing moment appeared in the θ = 30° parking position, which was 35.74 kN·m at the wind direction of β = 105°, 3.11 and 2.95 times that of the θ = 0° and θ = 60° parking positions respectively. According to these results, under typhoon side wind, not only do the averages of the total aerodynamic moments increase by a large amount, the fluctuation intensity is also greatly augmented, which may cause accumulated damage of wind turbine structures under low-cycle fatigue.

For the wind rotor in the θ = 30° parking position, the average value and RMS amplitude of the wind torque and pitching moment were all smaller, but the average value and RMS amplitude of the yawing moment were far greater than that of the other two parking positions. Moreover, as shown in Figure 17, it is clear that the effect of the interactions among the three blades of the wind rotor on the total loads is more significant in the θ = 30° parking position, which causes the RMS amplitudes of the total aerodynamic moments in the wind direction interval (−180°, 0°) to be much less than that in the wind direction interval (0°, 180°). In agreement with Figure 5, for the wind rotor in the θ = 30° parking position, when the wind direction β < 0°, the downwind blade B2 has an obviously contraction effect on the trailing vortex of the upwind blades B1 and B3, which greatly improves the flow state around blades B1 and B3, obviously decreasing the fluctuating wind load. What is more, the RMS amplitude of the total aerodynamic forces of the rotor in the θ = 30° parking position had a similar variation characteristic, but it was not as significant as that of the total aerodynamic moments. Therefore, the θ = 30° parking position is the optimum for wind turbines that adopt free yawing strategy, especially when the wind direction β < 0°. However, the fluctuation amplitude of the yawing moment of the rotor in the θ = 30° parking position increases rapidly in the wind direction interval (0°, 180°), which will add great difficulties to the implementation of the free yawing strategy; meanwhile, the fast changing wind directions also increase the frequency of yawing, which will also produce fatigue loads on the wind turbine. Thus, it is very necessary to elaborately match the yaw dampers in the nacelle.

Table 3. Range of the dominant frequency of the total aerodynamic moment of the wind rotor.

Dominant Frequency (Hz)		MX	MY	MZ
θ = 0°	β > 0°	0.080–14.960	0.080–14.960	0.080–14.960
	β < 0°	0.178–13.740	0.165–13.740	0.165–13.740
θ = 30°	β > 0°	0.085–14.480	0.174–6.087	0.127–6.087
	β < 0°	0.111–13.590	0.111–13.590	0.289–13.590
θ = 60°	β > 0°	0.055–11.930	0.055–11.930	0.055–6.480
	β < 0°	0.108–13.480	0.108–13.480	0.125–13.480

gives the range of the dominant frequency of the total aerodynamic moment, which is basically the same as that of the total aerodynamic force. It requires that the natural frequency of the wind rotor be outside these ranges in design to avoid resonance failures.

4.3. Aerodynamic Loads of the Blade

The blade vibration characteristics are directly related to the blade root bending moments [41], which are measured in the blade coordinate system {X_B, Y_B, Z_B}, as shown in Figure 18. The blade root bending moments can be divided into the flap-wise moment M_XB, edge-wise moment M_YB and twisting moment M_ZB. The flap-wise moment M_XB causes the flexural vibrations of blades in the direction perpendicular to the airfoil chord; the edge-wise moment M_YB causes the flexural vibrations of blades in the direction of the airfoil chord; and the twisting moment M_ZB causes the twisting vibrations of blades around the longitudinal central shaft, which can be seen as the pitch moment.

The variation of the average value of the flap-wise moment M_XB with wind direction is similar to a sinusoidal distribution. In the range of small wind direction angles, there is no stall on the blade, so the flap-wise moment M_XB rapidly increases with the increasing of the wind direction |β|; when |β| increases to a certain degree (30°–45°), a stall occurs, and then the flap-wise moment M_XB begins to decrease with the increasing of the wind direction |β|. The RMS amplitude of the flap-wise moment M_XB reached its maximum at the wind directions β = 60°–150° and −90°–−150°. For the wind rotor in the θ = 0° parking position, the average and RMS amplitude of the flap-wise moment M_XB on blade B3 were greater, the maximum average of 3545.43 kN·m appeared at the wind direction of β = −45°, and the maximum RMS amplitude of 99.70 kN·m appeared at the wind direction of β = −150°; for the wind rotor in the θ = 30° parking position, the average and RMS amplitude of the flap-wise moment M_XB on blade B3 were also greater, the maximum average of 3179.68 kN·m appeared at the wind direction of β = −30°, and the maximum RMS amplitude of 68.69 kN·m appeared at the wind direction of β = 120°; for the wind rotor in the θ = 60° parking position, the average and RMS amplitude of the flap-wise moment M_XB on blade B1 were greater, the maximum average of 2960.10 kN·m appeared at the wind direction of β = 45°, and the maximum RMS amplitude of 98.32 kN·m appeared at the wind direction of β = −105°.

The average of the edge-wise moment M_YB is significantly affected by the windward side of the blade. When the lower surface of the blade was facing the wind, the average edge-wise moment M_YB was greater than that when the upper surface of the blade facing the wind. Moreover, similar to the total aerodynamic force F_X, the edge-wise moment M_YB acted in the negative direction only when β = −5°–5° because of the special aerodynamic form of the blade. For the wind rotor in the θ = 0° parking position, the average and RMS amplitude of the edge-wise moment M_YB on blade B3 were greater, the maximum average of 1740.43 kN·m appeared at the wind direction of β = −105°, and the maximum RMS amplitude of 137.27 kN·m appeared at the wind direction of β = 105°; for the wind rotor in the θ = 30° parking position, the average and RMS amplitude of the edge-wise moment M_YB on blade B1 were greater, the maximum average of 1393.16 kN·m appeared at the wind direction of $\mathsf{\beta }=90°$, and the maximum RMS amplitude of 63.30 kN·m appeared at the wind direction of β = 105°; for the wind rotor in the θ = 60° parking position, the average and RMS amplitude of the edge-wise moment M_YB on blade B1 were greater, the maximum average of 1576.00 kN·m appeared at the wind direction of β = 105°, and the maximum RMS amplitude of 116.07 kN·m appeared at the wind direction of β = −105°.

The plus or minus sign of the twisting moment M_ZB that is used to indicate the direction is approximately identical to that of the flap-wise moment, as shown in Figure 18. With the increasing of |β|, the twisting moment M_ZB also increased first and then decreased, and the average value of the twisting moment M_ZB at the wind direction of β = ±180° was more than double that at the wind direction of β = 0°. For the wind rotor in the θ = 0° parking position, the average and RMS amplitude of the twisting moment M_ZB on blade B3 were greater, the maximum average of 135.40 kN·m appeared at the wind direction of β = 120°, which was 19.48 times that of the rated operating state, and the maximum RMS amplitude of 7.94 kN·m appeared at the wind direction of β = 105°; for the wind rotor in the θ = 30° parking position, the average and RMS amplitude of the twisting moment M_ZB on blade B3 were greater, the maximum average of 116.17 kN·m appeared at the wind direction of β = 120°, which was 16.71 times that of the rated operating state, and the maximum RMS amplitude of 4.28 kN·m appeared at the wind direction of β = 120° as well; for the wind rotor in the θ = 60° parking position, the average and RMS amplitude of the twisting moment M_ZB on blade B1 were greater, the maximum average of 117.32 kN·m appeared at the wind direction of β = −120°, which was 16.88 times that of the rated operating state, and the maximum RMS amplitude of 6.69 kN·m appeared at the wind direction of β = −105°. Thus, the design process of the pitch control system must consider not only the requirement of normal pitch control but also resistance to the twisting moment in side wind of typhoon, in order to prevent blades from being rotated into an operation position with a small pitch angle under the huge twisting moment, in that case, wind turbines may be running at high speed even if they are in braking status, which will result in the brake disks shattering, even causing the disintegration of the wind turbine. Moreover, it also should try to improve the torsional rigidity of the blades to avoid the torsional divergence in large attack angle conditions.

As stated above, the blade root bending moment is different in different positions. The blades in the vertical position are the most risky under typhoon activity, and the average values and RMS amplitudes of the flap-wise moment, edge-wise moment and twisting moment on them are all greater than the other blades. In addition, Figure 18 also shows the interactions among the three blades of the wind rotor, and it is mainly embodied in the RMS amplitude of the blade root bending moment, among which the most significant are blades of the rotor in the θ = 30° parking position. When the wind direction was in the range of β < 0°, the downwind blade B2 had an contraction effect on the trailing vortex of the upwind blades B1 and B3, making the fluctuating wind loads of the two blades obviously decrease, while the fluctuating wind load of blade B2 was always a small value because of its favorable horizontal position, thereby lowering the fluctuating wind loads acting on the whole wind rotor. Figure 19 provides the time-domain waveforms and the frequency spectrograms of the edge-wise moment of the wind rotor in the θ = 30° parking position, from which we can observe clearly the improvement of the load state on blades B1 and B3 by blade B2.

Table 4. Range of the dominant frequency of the blade root bending moment.

Dominant Frequency (Hz)		MXB		MYB		MZB
		β > 0°	β < 0°	β > 0°	β < 0°	β > 0°	β < 0°
θ = 0°	B1	0.292~5.100	0.140~15.860	0.101~5.100	0.140~19.450	0.444~14.960	0.140~19.450
	B2	0.080~25.310	0.140~15.690	0.080~5.100	0.140~6.417	0.080~20.330	0.140~5.268
	B3	0.080~18.400	0.165~13.740	0.080~14.960	0.165~13.740	0.08~18.400	0.165~13.740
θ = 30°	B1	0.085~7.692	0.089~13.590	0.137~6.250	0.557~13.590	0.085~6.250	0.089~13.590
	B2	0.213~6.250	0.180~5.395	0.170~6.250	0.315~5.202	0.236~6.250	0.180~5.202
	B3	0.127~14.480	0.180~26.980	0.127~14.480	0.111~6.852	0.127~14.480	0.180~7.818
θ = 60°	B1	0.055~11.930	0.108~13.480	0.055~11.930	0.108~13.480	0.055~11.930	0.108~13.480
	B2	0.055~22.820	0.077~5.128	0.055~5.061	0.108~5.128	0.055~8.417	0.566~5.128
	B3	0.128~8.003	0.108~6.727	0.323~6.480	0.108~5.128	0.055~6.480	0.108~6.727

gives the range of the dominant frequency of the blade root bending moment, and it requires that the natural frequency of the blade be outside these ranges in design to avoid resonance failures. It also shows that the interactions among the three blades of the wind rotor have an effect on the dominant frequency of the blade wind loads as well, which can cause high dominant frequency to occur on the upwind blade.

5. Conclusions

This study has used the CFD numerical simulation method to analyze the aerodynamic loads on offshore wind turbines under typhoon with full wind direction, based on the situation that the wind turbine was in braking-stop state and its yaw control was inactive due to the typhoon. As the wind direction of the typhoon changed from the feathering position, the total aerodynamic forces acting on the wind rotor approximately increased first and then decreased because of the special aerodynamic form of the blade. The parking positions of the wind rotor had a small effect on the average values of the total aerodynamic forces, but had a great effect on the fluctuation characteristic of the aerodynamic force. However, unlike the aerodynamic force, the average values and fluctuation characteristics of the total aerodynamic moments of the wind rotor were all significantly affected by the rotor parking positions. The wind rotor in the θ = 30° parking position had a better load state than the other parking positions, but the yawing moment was far greater than that of the other parking positions and very sensitive to wind direction changes. Therefore, it can be worthwhile to adopt a typhoon resistance strategy without electrical power, in order to avoid the yaw control system being inactive when the electric power network fails or the wind-measuring device breaks down under typhoon activity. This study suggests turning the wind rotor to the θ = 30° parking position, loosening the yaw brake, and utilizing the wind rotor as a tail vane to achieve free yawing, and then the wind rotor will finally be adjusted to the favorable wind direction under the function of a great yawing moment.

With direction changes of typhoon, it is obvious that blades at different positions have different attack angles, among which the attack angles of the blades in the vertical position are greater, therefore these blades are the most risky under typhoon activity. In addition, there exist obvious interactions among the three blades of the wind rotor. The downwind blades have a contraction effect on the trailing vortex of the upwind blades, greatly decreasing the fluctuating loads of the upwind blades. Especially, the interactions among the three blades of the wind rotor in the θ = 30° parking position is most significant: when the blade in the horizontal position is in the downwind, the flow state around the other two upwind blades will be greatly improved, showing obvious improvement in the load state of the wind turbine.

Limited by the computer capacity and huge computational grid number, there are some shortcomings and regrets in this study. CFD analysis adopting large eddy simulation (LES) with high-performance supercomputers or wind tunnel tests should be carried out in the future, in order to obtain more complete results of wind loads on offshore wind turbines under typhoon.