Dynamic Amplification of Gust-Induced Aerodynamic Loads Acting on a Wind Turbine during Typhoons in Taiwan

Abstract

Typhoons, such as Soudelor, which caused the collapse of several onshore wind turbines in 2015, pose a considerable challenge to Taiwan’s wind energy industry. In this study the characteristics of the aerodynamic loads acting on a wind turbine due to wind gusts in a typhoon are studied with a view to providing a proper definition of the S-Class wind turbine proposed in International Electrotechnical Commission (IEC) 61400-1. Furthermore, based on analysis of wind data during typhoons, as obtained from the meteorological mast in the Zhangbin coastal area, an extreme wind speed and gust model corresponding to the typhoon wind conditions in Taiwan are herein proposed. Finally, the flow fields around a parked wind turbine experiencing both an unsteady gust and a steady extreme wind were simulated by a numerical approach. Numerical results show that the aerodynamic shear force and overturning moment acting on the target wind turbine in a steady wind are significantly lower than those under an unsteady gust. The gust-induced amplification factors for aerodynamic loadings are then deduced from numerical simulations of extreme wind conditions.

1. Introduction

Taiwan, benefitted by its unique geographical location, has rich wind resources particularly in the region of the Taiwan Strait. This geographical advantage clearly favors large-scale wind farm development, but the frequent incidence of typhoons in Taiwan is a real threat to wind farm safety, having already caused significant damage to installed wind turbines in recent years. For example, Typhoon Soudelor in 2015 collapsed six wind turbines and seriously damaged the blades of a seventh wind turbine in the Taichung wind farm. The maximum wind speed at a nearby meteorological mast was reported as 62.6 m/s by Liu and Chen [1]. This event indicates that a wind turbine may not intactly survive the extreme winds of a typhoon if this scenario is not well-defined. Although both steady and turbulent wind conditions are considered in the Extreme Wind Speed Model (EWM) of the International Electrotechnical Commission (IEC) 61400-1 [2] standard, the corresponding wind characteristics and induced aerodynamic loads on wind turbines under such extreme conditions still require further investigations and inspired this study to address the dynamic effect of gusts during extreme typhoons.

In order to define the proper design requirements for wind turbines operating in typhoon-prone areas, several studies have, in recent decades, investigated the characteristic wind conditions of typhoons, such as wind profile, turbulence intensity, and gust factor, based on measured data. Ishizaki [3] proposed relationships between turbulence intensity, mean wind speed, ground height, and gust factor from a statistical analysis of typhoon measurements. Cao et al. [4] analyzed the wind conditions of Typhoon Maemi through the wind speed samples measured by nine vane-type and seven sonic-type anemometers at a height of 15 m. They found that the turbulence intensity decreases with increasing wind speed and remains almost constant at high wind speeds, with a reported gust factor of 1.6. Clausen et al. [5] proposed a method to characterize tropical cyclones and to derive the structural design wind speed at a given site based on the existing and publicly available cyclone data. Garciano and Koike [6] employed an extreme wind speed delivered by generalized extreme value distribution to estimate the buckling strength by a buckling capacity model recommended by the ISO. They further assessed the probability of buckling failure of a wind turbine considering extreme wind speed distributions from both typhoon-prone and non-typhoon-prone areas, as well as the buckling resistance of a tower as a function of wind speed. They proposed a reference wind speed for wind turbines in typhoon-prone areas based on the annual extreme wind speeds measured at fifty weather stations around the Philippines. Their reference wind speed is about 16% higher than that suggested for wind turbines of Class I.

Several studies have focused on the damage inflicted to wind turbines by typhoon-induced wind loads, with a specific focus on the fluid-structure interactions. Ishihara et al. [7] analyzed the damage to the wind turbines on Miyakojima Island after the onslaught of typhoon Maemi, where the estimated maximum gusts were over 70 m/s. They employed a finite element method (FEM) to forecast the displacement at the tower top together with the bending moment at the turbine foundation. Uchida et al. [8] investigated the cause of blade damage to wind turbines in southern Honshu when Typhoon Melor struck Japan in 2009. The WRF-ARW meteorological model was employed to predict the mesoscale flow behavior, followed by a large eddy simulation (LES) model, i.e., RIAM-COMPACT, to describe the near-field flow features, and, finally, a Reynolds-Averaged Navier–Stokes (RANS) model was coupled with an FEM model to determine the flow field around the blades and the resulting stress on the blades. In a recent study [9], the aerodynamic loads on a 5-MW wind turbine during an extreme gust, modelled in accordance with the IEC 61400-1 standard, were assessed via an unsteady RANS method, and the increase of blade loading was quantified and verified. The dynamic response of this 5-MW wind turbine on a floating platform in irregular seas was investigated via an aero-servo-elastic modeling approach [10]. From the perspective of dynamic analysis, Amaechi et al. [11] analyzed submarine hoses attached to a mooring buoy and suggested a dynamic amplification factor (DAF) of 2 on the environmental loads. Haddadin et al. [12] defined the DAF of a lattice structure as the ratio between the peak total response and the peak quasi-static response, which is the definition adopted in the present study for the aerodynamic response of a wind turbine.

Unfortunately, no local wind characteristics during typhoons in Taiwan are explicitly disclosed in the aforementioned references. For this reason, the present study first applied a generalized extreme value analysis to local wind measurements during typhoons to propose an extreme wind speed and a gust model for the Taiwan region. Figure 1 illustrates the framework of this paper. A measurement-based approach to describe the extreme wind conditions of typhoons is given in Section 2. Then, Section 3 uses this extreme condition to conduct steady and transient RANS simulations, and examines the aerodynamic loads on a parked wind turbine. Because of a lack of in-situ load measurements of a wind turbine in such an extreme typhoon in Taiwan, a benchmark study of the target wind turbine under its rated condition was performed to ensure the reliability of the RANS results. Finally, Section 4 discusses the load amplification factors due to gust effects based on the comparison of the steady and unsteady simulation results. Note that the amplification induced by the turbulent wind in EWM is not performed in the present study, so only the partial dynamic behaviour of this wind turbine is resolved by the simulations.

2. Typhoon Wind Condition

2.1. Wind Speed Measurement

Wind speed data measured by a meteorological mast located in the Zhangbin coastal area (24°06′ N 120°23′ E), from 2007 to 2015, were analyzed to determine the wind characteristics on the Taiwanese west coast. Vane-type anemometers, with sampling rates of 10 Hz and data logging intervals of one minute, were installed at heights of 10 m, 30 m, 50 m, and 70 m above ground, as shown in Figure 2. The accessible wind speed data consisted of ten-minute average $V_{600}$, one-minute average $V_{60}$, the maximum of 3-s averages taken over a minute $V_{\max }$, and the 10-min standard deviation.

In this study, the turbulence, wind shear, and gust during typhoons were considered. Turbulence is the phenomenon of random fluctuations of flow. The turbulence intensity describes the strength of the turbulence and is defined as the ratio of the wind speed standard deviation to the corresponding average wind speed. In IEC 61400-1, the wind shear is defined as the variation of wind speed with height above ground and is modelled by the power law. Equation (1) [2] and Figure 3 give the mathematical and graphical description of wind shear, respectively, where $V$ is the wind speed, $z$ is the height above the ground, $z_{\mathrm{ref}}$ is the reference height, and $V_{\mathrm{ref}}$ is the velocity at the reference height. The power law exponent of the meteorological mast $\alpha$ is deduced as 0.3, according to the regression fit to the anemometers at different heights [13].

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

(1)

A gust is meteorologically defined as the wind condition where the instantaneous wind speed is at least 5 m/s faster than the corresponding average wind speed. The World Meteorological Organization (WMO) [14] describes gust magnitude in terms of the gust factor $GF$, which is defined as the ratio of the peak $\tau$-second average wind speed $V_{\tau }$ to the corresponding $T_0$-second average wind speed $V_{T_0}$, as shown in Equation (2), where $T_0$ is a base reference observation period. To comply with the specifications of the data logger in the device, $\tau$ was set as 3 s and $T_0$ was taken as 60 s to compute $GF$.

$$GF\left(\tau ,T_0\right)\equiv \frac{V_{\tau }}{V_{T_0}}$$

(2)

Following from the local wind measurements, the gust factor $GF$ and the turbulence intensity $I$ of the recorded data from the meteorological mast in the Zhangbin area during typhoons from 2007 to 2014 are shown in Figure 4a,b. According to Equation (2), the gust factors were calculated from the maximum 3-s average wind speeds and 1-min average wind speeds, and the turbulence intensities were calculated from the wind speed 10-min standard deviations and the 10-min average wind speeds. Since $GF$ and $I$ refer to different averaging durations of wind speed, a correlation between $V_{60}$ and $V_{600}$, Figure 4c, is further required so that we can convert $I$ into terms of $V_{60}$, so as to have the same argument as $GF$.

2.2. Extreme Wind Condition

In order to find the extreme condition of the gust factor and turbulence intensity, the gust factor and the turbulence intensity of the recorded data during typhoons were processed by the following steps, as shown in the flowchart of Figure 5:

• Filter 1: The raw data were firstly filtered by removing zeros and empty fields. High-frequency noise was trimmed and low gusts that did not meet this definition were removed.
• Grouping: The recorded data were divided into 1 m/s bins, employing the 1-min average wind speeds for gust factors and 10-min average wind speeds for turbulence intensities.
• Filter 2: The data in each interval were assumed to fit a normal distribution. Unreasonable values and outliers of each interval were filtered by box-whisker plot. This filter performed the calculation of the first (Q1) and third (Q3) quartiles of the gust factor and turbulence intensity, and then the interquartile range (IQR) was computed. Data values that exceeded Q3 by three times IQR or were less than Q1 by three times IQR were identified as outliers and removed.
• Chi-squared test: The maximum cumulative probability of all the intervals filtered by chi-squared test was applied to estimate the extreme value of each interval by the assumption of normal distribution with a confidence interval of 95%.
• Fitting: The fitting curve of the estimated extreme values of each wind speed group was calculated by an exponential function with a least-squares fitting algorithm.

The extreme values of the gust factor and the turbulence intensity during typhoons are shown in Figure 6. The fitting curve for the gust factor at 99.999% extreme values is given by Equation (3), while Equation (4) expresses the fitting curve for 99.99999% extreme values of turbulence intensity.

$$GF=1+8.4177·V_{60}{}^{-0.9702}$$

(3)

$$I=1.0231·V_{600}{}^{-0.5715}$$

(4)

Equations (3) and (4) are the calibrated curves based on local measurement in the Zhangbin area, and the corresponding comparison to the IEC 61400-1 design standards for gust factor and turbulence intensity are shown in Figure 7. One can easily see that the maximum wind speeds during typhoons in the Zhangbin area clearly exceed those suggested in the IEC standards. The present gust model predicts about 20 m/s higher maximum wind speeds. The stronger gust fluctuations might impose higher wind loads on the turbine, and this is the main reason we proceeded further with numerical simulations based on Equation (3), instead of the value defined in IEC standards. Additionally, the turbulence intensities at hub height also differ from values adopted in IEC standards. During typhoons, the turbulence intensity is found higher than any of the design turbulence types at low wind speeds below 30 m/s, but it decreases quickly with increasing wind speed. Taking a wind speed of 36 m/s, for example, the corresponding turbulence intensity is 0.13, which is close to the type B defined in the IEC standards.

2.3. Reconstruction of Transient Gust

In the IEC standards, the extreme wind speed is defined as the maximum value of 3-s averages taken over a minute, which was 45.36 m/s as obtained from the typhoon wind measurements in the Zhangbin area. The local gust duration time T, suggested by a related study [15], was 6 s. By utilizing the extreme wind conditions, i.e., Equations (2) and (3), $V_{60}$ was iteratively solved and the corresponding gust factor $GF$ was subsequently determined as 1.26.

There is a difficulty in reconstructing this extreme $GF$ to a transient gust time series, which is explicitly indicated by the extreme wind speed model (EWM) in the IEC 61400-1 standard that unsteadiness should be considered. The anemometer data logger did not record the short transient state of a gust. Furthermore, the extreme operating gust model (EOG) in IEC 61400-1 [2] is only applied to describe the gust for a wind turbine in operating condition, but not in the parked condition, such as during typhoons. To overcome the mismatch, the transient profile of the EOG model, $g\left(t\right)$, was adopted, and the model constant K was calibrated by the local measurements of $GF$ obtained in the Zhangbin area, as per Equation (5), where $t$ is time, $T$ is the duration of gust, and $V_0$ is the initial wind speed.

$$V\left(t\right)=\{\begin{array}{c}{V}_0\left[1-K·g\left(t\right)\right], 0\le t<T\\ V_0 , t\ge T\end{array}, \mathrm{where} g\left(t\right)=\sin \frac{3\pi }{T}t\left(1-\cos \frac{2\pi }{T}t\right)$$

(5)

Then, in order to express K in terms of $GF$, the peak average wind speed $V_{\tau }$, and the average wind speed over gust duration, $V_T$, are subsequently expressed by Equations (6) and (7). Additionally, $V_T$ has to be transformed into the base reference observation period time $V_{T_0}$, such that it is the same as in the definition of $GF$, as per Equation (8). Substituting $V_{\tau }$ and $V_{T_0}$ (Equations (6) and (8)) into Equation (2), and expressing the model constant K explicitly yields Equation (9). Finally, using the extreme $GF$ in the Zhangbin area as calculated in the previous section, 1.26, $K$ was calculated by Equation (9) as 0.4.

$$\begin{array}{c}{V}_{\tau }=\frac{1}{\tau }{{\displaystyle \int }}_{\frac{T-\tau }{2}}^{\frac{T+\tau }{2}}V\left(t\right)dt=V_0\left(1-K\frac{G\left(\tau \right)}{\tau }\right), \mathrm{where}\\ G\left(\tau \right)={{\displaystyle \int }}_{\frac{T-\tau }{2}}^{\frac{T+\tau }{2}}g\left(t\right)dt=-\frac{T}{\pi }\left(\sin \frac{\pi }{2}\frac{\tau }{T}+\frac{2}{3}\sin \frac{3\pi }{2}\frac{\tau }{T}+\frac{1}{5}\sin \frac{5\pi }{2}\frac{\tau }{T}\right)\end{array}$$

(6)

$$V_T=\frac{1}{T}{{\displaystyle \int }}_0^{T}V\left(t\right)dt=V_0\left(1+\frac{8}{15\pi }K\right)$$

(7)

$$V_{T_0}=\frac{V_0\left(T_0-T\right)+{\overline{V}}_{T}T}{T_0}=V_0\left(1+K\frac{8}{15\pi }\frac{T}{T_0}\right)$$

(8)

$$K=-\frac{GF\left(\tau ,T_0\right)-1}{\frac{G\left(\tau \right)}{\tau }+\frac{8}{15\pi }\frac{T}{T_0}GF\left(\tau ,T_0\right)}$$

(9)

After all the required parameters for Equation (5) were determined, the extreme gust transient profile during typhoons in the Zhangbin area was illustrated, as in Figure 8. The related wind speeds following this extreme gust profile are summarized in

Table 1. Wind speed parameters.

Wind Speed (m/s)	z = 70 m (Met. Height)	z = 90 m (Hub Height)
$V_0$	35.75	38.55
$V_{60}$	36.00	38.82
$V_3$	45.36	48.91
$V_1$	61.41	66.22
$V_{\max }$	64.85	69.93

, where $V_{\max }$ was used in the steady simulation, i.e., the extreme wind condition. The measured maximum 1-s average wind speed of the Soudelor typhoon in 2015 was reported as 59.4 m/s [1], which agrees well with the present gust model with an error of 3%. With the help of Equation (1), the wind speeds from the meteorological mast height of 70 m were able to be extrapolated to the hub height of 90 m. So, the spatial and temporal distributions of wind speed were known for computing aerodynamic loads on a wind turbine. The turbulence intensity during a typhoon in the Zhangbin area was estimated as 0.13 from Equation (4) as input to the following simulations.

3. Wind Load Assessment

3.1. Wind Turbine Model

The target wind turbine is a 5-MW horizontal-axis wind turbine proposed by the National Renewable Energy Laboratory (NREL) [16], and its geometric and operational parameters are given in

Table 2. Main particulars of the target wind turbine.

Rotor type	Upwind, 3 blades
Diameter ($D$)	126 m
Hub height ($H_{\mathrm{hub}}$)	90 m
Hub diameter ($D_{\mathrm{hub}}$)	3 m
Tilt angle (${\theta }_t$)	5°
Pitch angle (${\theta }_p$)	Operating: 0°~23.5°; Parked: 90°
Rated wind speed	11.4 m/s
Rated rotor speed	12.1 rpm

. The tilt angle is defined as the angle between the rotor’s rotation axis and the ground level (Figure 9a) and the pitch angle is defined as the angle between the rotor’s rotation plane and the chord line of the blade tip (Figure 9b). The detailed section profiles for each radial position of the blade are given in [16]. According to the blade geometry and operating conditions of the target wind turbine, the Mach number is smaller than 0.3 and the Reynolds number ranges between $5\times {10}^6$ and $1.6\times {10}^7$, which justifies the turbulent flow simulation employed in this study.

During typhoons, the wind turbine is designated as “parked”. The rotor is then stationary with a blade pitch angle of $90°$ and the yaw system keeps the wind turbine facing into the wind to reduce the induced aerodynamic loads. The investigated aerodynamic loads acting on the wind turbine are given in Figure 9c, where ($F_{\mathrm{x}}$,$F_{\mathrm{y}}$,$F_{\mathrm{z}}$) denote the resultant forces acting on the wind turbine along three coordinate directions, ($Q_{\mathrm{x}}$,$Q_{\mathrm{y}}$,$Q_{\mathrm{z}}$) denote the resultant moments with respect to three coordinate axes, $Q_{\mathrm{A}}$ denotes the rotor moment with respect to the rotor axis, and ($Q_{\mathrm{p}1}$,$Q_{\mathrm{p}2}$,$Q_{\mathrm{p}3}$) denote three blade pitch moments with respect to the blade axis.

Owing to a nontrivial tilt angle, the rotor axis is principally not parallel to the incoming wind, so that the angle of attack at each blade section is not constant about the rotor axis. Figure 10 illustrates the section profile in the A-A plane, together with the relative velocity vector triangle, of a blade with azimuthal position ${\theta }_a$. The change of the angle of attack due to the tilt angle is described by Equation (10) that indicates the change of the angle of attack of the blade section decreasing with an azimuthal angle between 0° and 180° (on the right-hand side of sweeping disk), and increasing with an azimuthal angle between 180° and 360°. Therefore, the distribution of the rotor moment $Q_{\mathrm{A}}$ and pitch moment of the blade ($Q_{\mathrm{p}1}$,$Q_{\mathrm{p}2}$,$Q_{\mathrm{p}3}$) are expected to be similar to a sine function. The peak loads at certain azimuthal angles are then obtained in the steady simulations.

$$\Delta \alpha ={\tan }^{-1}\left(\frac{\sin {\theta }_t\sin {\theta }_a}{\cos {\theta }_t}\right)$$

(10)

3.2. Numerical Method

The characteristics of the aerodynamic loads acting on a wind turbine during typhoons were studied by numerical simulations. The flow solver ANSYS Fluent [17], a general-purpose finite volume and RANS code, was applied to simulate the flow field around the target wind turbine. The wind conditions during a typhoon were described by the results of the statistical analysis of the recorded data from the meteorological mast in the Zhangbin area. The flow around wind turbines is considered an incompressible and turbulent flow. The corresponding governing equations for conservation of mass and momentum are, respectively, as per Equations (11) and (12), where $u_i$ is the mean velocity component in the $x_i$ direction, $p$ is the pressure, $\rho$ and $\mu$ are the density and viscosity of air, respectively, $\rho \overline{u_i^{\prime }{u}_j^{\prime }}$ is the Reynolds stress, and $\rho b_i$ is the gravitational force component in the $x_i$ direction. The SST k-ω model, as per Equations (13) and (14), has been widely accepted for simulating flow past airfoils [9] and was applied in this study to compute the Reynolds stress term in Equation (12), where $k$ is the turbulent kinetic energy, $\omega$ is the dissipation rate of $k$, ${\tau }_{ij}$ is the shear stress, ${\mu }_t$ is the turbulent viscosity, and (${\beta }^{*}$, ${\sigma }_k$, $\gamma$, $\beta$, ${\sigma }_{\omega }$, $\delta$, ${\sigma }_{\omega 2}$) are the equation constants.

$$\frac{\partial u_i}{\partial x_i}=0$$

(11)

$$\frac{\partial u_i}{\partial t}+\frac{\partial }{\partial x_j}\left(u_i{u}_j\right)=-\frac{1}{\rho }\frac{\partial p}{\partial x_i}+\frac{1}{\rho }\frac{\partial }{\partial x_j}\left[\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\right]-\frac{\partial \left(\overline{u_i^{\prime }{u}_j^{\prime }}\right)}{\partial x_j}+b_i$$

(12)

$$\frac{\partial k}{\partial t}+\frac{\partial }{\partial x_j}\left(u_{j}k\right)=\frac{1}{\rho }{\tau }_{ij}\frac{\partial u_i}{\partial x_j}-{\beta }^{*}\omega k+\frac{1}{\rho }\frac{\partial }{\partial x_j}\left[\left(\mu +{\sigma }_k{\mu }_t\right)\frac{\partial k}{\partial x_j}\right]$$

(13)

$$\frac{\partial \omega }{\partial t}+\frac{\partial }{\partial x_j}\left(u_j\omega \right)=\frac{\gamma }{{\mu }_t}{\tau }_{ij}\frac{\partial u_i}{\partial x_j}-\beta {\omega }^2+\frac{1}{\rho }\frac{\partial }{\partial x_j}\left[\left(\mu +{\sigma }_{\omega }{\mu }_t\right)\frac{\partial \omega }{\partial x_j}\right]+2\left(1-\delta \right){\sigma }_{\omega 2}\frac{1}{\omega }\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j}$$

(14)

The computational domain of the numerical model is shown in Figure 11a. There is a distance of two times the rotor diameter between the upstream boundary and the wind turbine, six times the rotor diameter between the wind turbine and the downstream boundary, 1.5 times the rotor diameter between the wind turbine and each of the side boundaries, and 3 times the rotor diameter between the bottom and top boundaries, as listed in

Table 3. Geometrical parameters of computational domain.

Lh	Lw	L1	L2	DII	HII
3D	3D	2D	6D	1.27D	1.31D

. Figure 11b illustrates the locations for the adopted boundary conditions. The bottom boundary and the wind turbine surface are defined as no-slip walls. The upstream, top, and side faces are defined as inlet, where the inflow condition is specified according to the calibrated gust condition. Assuming that the distance between the wind turbine and the downstream boundary is far enough, a vanishing velocity gradient is employed at the outlet. The boundary conditions for velocity are summarized in

Table 4. Boundary conditions.

Boundary	Type	Velocity Condition
abcd (upstream)	Inlet	Velocity specified
abfe, aehd, bcgf (far field)	Inlet	Velocity specified
ghef (downstream)	Outlet	Zero velocity gradient
cdhg (ground)	No-slip wall	Zero velocity
Wind Turbine Surface	No-slip wall	Zero velocity

.

3.3. Verification and Validation

As the study utilizes a computational model, the verification process is essential to mitigate modelling errors within tolerance. For the purpose of verification, a steady wind speed of 60 m/s at hub height with the exponential profile $\alpha$ = 0.3 was used. The second-order upwind scheme was used for spatial discretization. Since the real geometry of the wind turbine is considered, a multi-block unstructured body-fitted grid arrangement was employed. The grid distribution is shown in Figure 12, where a refined cylindrical region was generated around the turbine. Several prismatic layers were fit around the body surfaces, and the thickness of the first layer was tuned to achieve y⁺ less than 10 for the requirement of the selected SST turbulence model. A grid dependence test, consisting of five levels of systematic grid refinement, was conducted to verify that the solution is numerically convergent. The aerodynamic moment $Q_{\mathrm{y}}$ about the tower base was selected as the target function in terms of grid number. Figure 13a displays the dependence of $Q_{\mathrm{y}}$ on grid number, where $Q_{\mathrm{y}}$ approached a constant value with an increasing grid number. The discretization error $E$ for a method of second-order accuracy is estimated through Richardson’s extrapolation, where $dx$ denotes the dimensionless grid size. Figure 13b indicates that the model with grid number over 10 million is capable of reaching a discretization error less than 2%.

The proposed numerical model is then validated for the rotor torque $Q_{\mathrm{A}}$ at rated condition, namely wind speed 11.4 m/s, rotor speed 12.1 rpm, and zero blade pitch angle where an unsteady simulation was performed. The time step was 5 ms and it corresponded to 1000 steps per revolution or a blade rotation of 0.36° per time step. The time step was small enough to accurately resolve the temporal variation of wind turbine aerodynamics. A sliding mesh approach was applied to model the rotating rotor disk to cope with the relative motion between the rotor and tower in the unsteady flow simulation. Figure 14a depicts the evolution of $Q_{\mathrm{A}}$ with respect to time. After a simulation time of 40 s, i.e., 10 turns of the rotor, the rotor torque converges. The mean value of $Q_{\mathrm{A}}$ over the 11th revolution of the rotor was 4.09 MN-m, which was about 2% less than the benchmark value of 4.18 MN-m. Figure 14b shows the velocity contours and vectors on the midplane between the rotor plane and the tower as well as the pressure contours on the blades for the blade configuration of 0°–120°–240°. The axial wake region induced by the rotor is clearly seen and a radial wake expansion is implied by outward vectors across the rotor disk boundary.

4. Results and Discussion

To investigate the characteristics of the aerodynamic loads acting on the wind turbine, several cases with different wind conditions were simulated. The predicted cases were separated into two categories: in the first category, the flow field around the parked wind turbine was simulated with constant extreme wind speed. These simulations were conducted in steady mode since the rotor was not moving and a mean wake was of main interest. Flow separation may occur behind the tower, but it is expected to be insignificant, especially in high wind conditions. In the other category, the flow field around the parked target wind turbine was simulated with a gust in transient mode. To enhance the numerical stability of the transient simulations, a precursor calculation of constant wind speed was first conducted. Then, one gust period with a varying wind speed was simulated using the precursor flow field as an initial condition. Four blade configurations were studied in the steady simulation, while only one blade configuration was investigated in the gust simulation.

4.1. Steady Aerodynamic Loadings

The simulation results of steady extreme wind speed are summarized in

Table 5. Wind loads under steady extreme wind.

Load (kN/MN-m)	Part	Blade Configuration
		0°–120°–240°	30°–150°–270°	60°–180°–300°	90°–210°–330°
$F_{\mathrm{x}}$	Tower	276.0	280.3	306.0	284.3
	Rotor	172.5	175.7	146.6	177.3
	Total	451.6	458.9	454.5	465.5
$Q_{\mathrm{A}}$	Rotor	3.80	5.00	1.20	2.50
$Q_{\mathrm{p}}$	Blade 1	−0.20	−0.23	−0.23	−0.21
	Blade 2	−0.17	−0.17	−0.10	−0.12
	Blade 3	−0.13	−0.12	−0.13	−0.17
$Q_{\mathrm{x}}$	Tower	−3.7	−1.6	−10.8	−1.2
	Rotor	−1.0	5.2	1.8	4.0
	Total	−5.0	3.7	−9.1	3.1
$Q_{\mathrm{y}}$	Rotor	18.1	18.2	15.8	17.8
	Tower	13.8	13.9	15.9	14.2
	Total	32.3	32.5	31.9	32.5
$Q_z$	Rotor	0.60	1.30	−0.09	0.50
	Tower	<0.001
	Total	0.60	1.30	−0.06	0.60

. The drag force $F_{\mathrm{x}}$ acting on the tower as well as rotor is insignificantly impacted by the blade configuration, except the case of 60°–180°–300° when the tower is directly shadowed by a blade. In this case the rotor shows a lower aerodynamic loading because of a small mean blade height implying a low incoming wind speed along with zero angle of attack, whereas the tower loading slightly grows due to a local flow acceleration due to the blockage effect of the blade. However, the total drag of the tower and rotor delivers a similar magnitude. The rotor torque $Q_{\mathrm{A}}$ is relatively sensitive to the blade configuration and fluctuates in the range between 1.2 MN-m and 5.0 MN-m where the blade configuration of 60°–180°–300° clearly shows the smallest value benefitted from its advantage in mean blade height and effective angle of attack. The blade pitch moment $Q_{\mathrm{p}}$ is approximately −0.167 MN-m in average, where the minus sign indicates aerodynamic loading tending to decrease the blade pitch angle, i.e., to lead to an unfavorable blade position to acquire higher blade pitch moment. In the range of $0°\le {\theta }_a\le 180°$ the blade pitch moment generally declines with the blade azimuthal angle, while the blade pitch moment grows with the blade azimuthal angle for $180°\le {\theta }_a\le 360°$. As expected, the blade pitch moment at ${\theta }_a=180°$ is the smallest among all studied blade positions due to a small local wind speed as well as angle of attack. The wind turbine suffers from a higher overturning risk along the y axis than along the x axis because $Q_{\mathrm{y}}$ is much higher than its counterpart $Q_{\mathrm{x}}$ that constantly changes its direction in a rotor rotation period. Additionally, $Q_{\mathrm{y}}$ varies little among different blade configurations and has a mean value of 32.3 MN-m. The overturning moment along the y axis is mainly contributed by the rotor and tower where they deliver a comparative importance in $Q_{\mathrm{y}}$. The twisting moment $Q_{\mathrm{z}}$ is mostly governed by the rotor but its magnitude is about one-order of magnitude smaller than $Q_{\mathrm{x}}$ and approximately two-order of magnitude less than $Q_{\mathrm{y}}$.

4.2. Dynamic Amplification

The gust simulation results for the blade configuration of 0°–120°–240° are shown in

Table 6. Wind loads under unsteady extreme gust for the blade configuration of 0°–120°–240°.

Load (kN/MN-m)	Part	At t = 3 s	Maximum	tmax (s)
$F_{\mathrm{x}}$	Tower	400.3	458.5	2.70
	Rotor	209.8	223.0	2.86
	Total	621.3	693.8	2.74
$Q_{\mathrm{A}}$	Rotor	5.11	5.31	2.87
$Q_{\mathrm{p}}$	Blade 1	−0.19	−0.19	3.01
	Blade 2	−0.19	−0.19	3.00
	Blade 3	−0.16	−0.16	2.93
$Q_{\mathrm{x}}$	Total	0.76	−1.78	2.28
$Q_{\mathrm{y}}$	Total	39.16	42.04	2.80
$Q_{\mathrm{z}}$	Total	0.80	0.84	3.06

, which lists the aerodynamic loading at the time of peak wind speed, i.e., $t$ = 3 s, as well as the time instance of the maximum loading during the unsteady gust. Wind profile is depicted in Figure 15, which also shows the low speed wake behind the blades and tower. The pressure contour at $t$ = 2 s is plotted on the right-hand side of Figure 15 to illustrate a positive pressure gradient along the wind direction during the acceleration phase of the gust.

In order to investigate the correlation between the aerodynamic loads and wind speed in the gust cases, the amplification factor, $DAF$, is defined as follows [12]:

$$DAF=\frac{\varphi }{{\left|\varphi \right|}_{\max }}$$

(15)

where $\varphi$ denotes any physical quantity, and ${\left|\varphi \right|}_{\max }$ is the maximum absolute value of the physical quantity $\varphi$. Both the rotor and tower drags reach their peak values somewhat in advance of the 3-s peak. Figure 16a compares the time history of different drag components in a gust period where top represents the contribution of rotor-nacelle assembly. The aerodynamic loading is found to principally mimic the wind speed variation. The tower drag force plays a major role in the total drag while the contribution of nacelle is relatively trivial. The rotor torque $Q_{\mathrm{A}}$ behaves similar to the rotor drag force and $Q_{\mathrm{A}}$ precedes the gust peak by 0.13 s. Different from the rotor torque, the peak blade pitch moments echo the arrival of the top wind speed. Figure 16b displays the blade pitch moment of the blade configuration of 0°–120°–240°. The wind unsteadiness interestingly leads to a more uniform peak blade pitch moment among the three blades accompanied by a slight increase in the mean value. The overturning moment $Q_{\mathrm{y}}$ and $Q_{\mathrm{x}}$ show a leading phase of the gust peak while the phase of $Q_{\mathrm{z}}$ lags the gust peak by 0.16 s. The phase inconsistency between the maximum aerodynamic loadings and the gust peak is a contribution of time terms in the momentum equations. In the unsteady flow the flow velocity as well as its time derivative contributes to the stagnation pressure whereas the contribution of time derivative vanishes in the steady case. Figure 17 depicts the normalized profiles of wind velocity, the time derivative of wind velocity and the drag force experienced by the wind turbine in a gust period where the superscript * represents the physical quantity normalized by its maximum value in a gust period. Figure 17 indicates that the wind velocity just experiences a sharp deceleration process as the flow velocity slowly reaches the gust peak. The combining effect of a gradually rising wind velocity and a steep velocity gradient in time clearly leads to a peak-load offset away from the time instance of gust peak, especially the aerodynamic load is governed by the stagnation pressure. For the aerodynamic load mainly stemming from viscous shear, such as the blade pitch moment and the twisting moment, this offset behavior is less unnoticeable.

Table 7. Dynamic amplification factor for the wind turbine with a blade configuration of 0°–120°–240°.

Load (kN/MN-m)	Gust	Steady Wind	$DAF$
$F_{\mathrm{x}}$	693.8	451.6	1.54
$Q_{\mathrm{A}}$	5.31	3.80	1.40
$Q_{\mathrm{p}}$	0.18	0.167	1.10
$Q_x$	1.78	5.00	0.36
$Q_{\mathrm{y}}$	42.04	32.27	1.30
$Q_{\mathrm{z}}$	0.84	0.62	1.35

summarizes the amplification factor of aerodynamic loads given in Table 5 and Table 6. The ratio of the peak tower drag in gust to the steady tower drag is approximately proportional to the square of the gust factor, $GF$, but the corresponding ratio for rotor drag varies quite linearly with the gust factor. This is explained by the blunt shape of the tower and the streamlined geometry of the rotor blade where the former is governed by stagnation pressure characterized by the square of flow velocity and the latter is dominated by viscous shear governed by the velocity gradient. The rotor torque $Q_{\mathrm{A}}$ is mainly contributed by the lateral force exerted on the blades where pressure and shear components are both important in the head wind condition. This accounts for an amplification factor of $1.4$ falling between $GF$ (shear dominant) and $G{F}^2$ (pressure dominant). The unsteady amplification of blade pitch moment is not very significant where only 10% increase is found for this blade configuration (0°–120°–240° or $\lambda$-shaped). The gust exhibits a positive impact on $Q_{\mathrm{x}}$ with a peak reduction by 64% in a $\lambda$-shaped blade configuration. Contrary to $Q_{\mathrm{x}}$, the extreme moments along other two axes, i.e., $Q_{\mathrm{y}}$ and $Q_{\mathrm{z}}$, are governed by the shear effect where the corresponding amplification factor is close to the gust factor.

5. Conclusions

Taiwan, due to its geographical location, has an excellent wind conditions and high potential for wind power production as well as strong typhoons which have damaged several wind turbines in recent years. The aerodynamic loads during typhoons were investigated by using the extreme wind conditions prescribed in IEC 61400-1 together with a measurement-calibrated gust factor. The extreme wind conditions proposed in the present study were based on meteorological measurements from the Zhangbin area during typhoon invasions from 2007 to 2014. A statistical approach was employed to convert the raw data into a fitting formula to quantitatively describe the extreme wind conditions for typhoons. Our study shows that the 1-min average wind speed was 36 m/s, the maximum 3-s average wind speed was 45.36 m/s, and the maximum instantaneous wind speed was 64.85 m/s, where the maximum instantaneous wind speed in Taiwan was about 20 m/s higher than the values suggested in the IEC standards. Additionally, the gust factor for typhoons was estimated at 1.26 with a gust period of 6 s. This local gust factor was used to calibrate the EOG model constant, such that a transient gust time series was reconstructed for the unsteady simulation. The turbulence intensity during typhoons was found to be higher than those of turbulence types defined in the IEC standards at low wind speeds below 30 m/s, but it declines quickly with the growth of wind speed.

The aerodynamic loads acting on the NREL 5-MW wind turbine in a steady extreme wind condition as well as an unsteady extreme gust were studied by RANS simulations. The simulation results show that the extreme aerodynamic loads acting on the wind turbine are obviously underestimated when the extreme wind condition only adopts a constant wind speed. The amplification factor of aerodynamic loads is predicted as follows: 1.54 for the drag force, 1.4 for the rotor torque, 1.3 for the overturning moment along the y direction, 1.1 for the blade pitch moment, and the yawing moment by 1.35 times. The amplification factor is governed by the nature of aerodynamic loading. It is close to the square of the gust factor in a stagnation-pressure dominant case while it approximates to the gust factor when viscous shear plays a major role. The quantitative amplification factors may not be universal for different wind turbines in different areas because the present study is limited to the Taiwan local wind conditions and this specific wind turbine. With a concern of lacking validation under extreme typhoons in this study, more load measurements of in-situ wind turbines are expected in the future. Lastly, it is suggested to conduct future studies on the dynamic amplification of this reference turbine by adopting the proposed extreme turbulence intensity in EWM, such that the full aerodynamics can be resolved and compared.