**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 gus<sup>t</sup> factor, based on measured data. Ishizaki [3] proposed relationships between turbulence intensity, mean wind speed, ground height, and gus<sup>t</sup> 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

261

**Citation:** Lin, T.-Y.; Yang, C.-Y.; Chau, S.-W.; Kouh, J.-S. Dynamic Amplification of Gust-Induced Aerodynamic Loads Acting on a Wind Turbine during Typhoons in Taiwan. *J. Mar. Sci. Eng.* **2021**, *9*, 352. https://doi.org/10.3390/jmse9040352

Academic Editors: Dong-Sheng Jeng and Raúl Guanche García

Received: 23 February 2021 Accepted: 22 March 2021 Published: 24 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

speed and remains almost constant at high wind speeds, with a reported gus<sup>t</sup> 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 typhoonprone 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 typhooninduced 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 gus<sup>t</sup> 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 gus<sup>t</sup> 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.

**Figure 1.** Framework of this study.

#### **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.

**Figure 2.** Meteorological mast on the Zhangbin coastal area.

In this study, the turbulence, wind shear, and gus<sup>t</sup> 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*ref is the reference height, and *V*ref is the velocity at the reference height. The power law exponent of the meteorological mast *α* is deduced as 0.3, according to the regression fit to the anemometers at different heights [13].

> *zz*ref *α*

*<sup>V</sup>*(*z*) =

$$V(z) = V\_{\rm ref}(z)$$
 
$$V(z) = \frac{\int\_0^{z\_{\rm ref}} dz}{z}$$
 
$$\approx \frac{\int\_0^{z\_{\rm ref}} dz}{z}$$
 
$$\approx \frac{\int\_0^{z\_{\rm ref}} dz}{z}$$
 
$$\approx \frac{\int\_0^{z\_{\rm ref}} dz}{z}$$

**Figure 3.** Wind shear profile.

A gus<sup>t</sup> 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 gus<sup>t</sup> magnitude in terms of the gus<sup>t</sup> factor *GF*, which is defined as the ratio of the peak *τ*-second average wind speed *Vτ* to the corresponding *T*0-second average wind speed *VT*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, *τ* was set as 3 s and *T*0 was taken as 60 s to compute *GF*.

$$GF(\tau, T\_0) \equiv \frac{V\_\tau}{V\_{T\_0}}\tag{2}$$

(1)

Following from the local wind measurements, the gus<sup>t</sup> 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 gus<sup>t</sup> 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 argumen<sup>t</sup> as *GF*.

**Figure 4.** Wind measurements. (**a**) Gust factor; (**b**) turbulence intensity; (**c**) wind speed correlation.

#### *2.2. Extreme Wind Condition*

In order to find the extreme condition of the gus<sup>t</sup> factor and turbulence intensity, the gus<sup>t</sup> 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:


**Figure 5.** Data processing procedure for wind measurements.

The extreme values of the gus<sup>t</sup> factor and the turbulence intensity during typhoons are shown in Figure 6. The fitting curve for the gus<sup>t</sup> 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 \cdot V\_{60} \, ^{-0.9702}\tag{3}$$

$$I = 1.0231 \cdot V\_{600}{}^{-0.5715} \tag{4}$$

**Figure 6.** Extreme value analysis during typhoons. (**a**) Gust factor; (**b**) turbulence intensity.

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 gus<sup>t</sup> 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 gus<sup>t</sup> model predicts about 20 m/s higher maximum wind speeds. The stronger gus<sup>t</sup> 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.

**Figure 7.** Wind condition comparison between local statistics and design standard. (**a**) Gust factor; (**b**) turbulence intensity.

#### *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 gus<sup>t</sup> 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 gus<sup>t</sup> factor *GF* was subsequently determined as 1.26.

There is a difficulty in reconstructing this extreme *GF* to a transient gus<sup>t</sup> 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 gus<sup>t</sup> model (EOG) in IEC 61400-1 [2] is only applied to describe the gus<sup>t</sup> 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*(*t*), 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(t) = \begin{cases} \begin{array}{c} V\_0[1 - K \cdot \mathcal{g}(t)], \; 0 \le t < T\\ V\_0, \; t \ge T \end{array}, \text{ where } \mathcal{g}(t) = \sin \frac{3\pi}{T} t \left(1 - \cos \frac{2\pi}{T} t\right) \end{array} \tag{5}$$

Then, in order to express *K* in terms of *GF*, the peak average wind speed *V<sup>τ</sup>*, and the average wind speed over gus<sup>t</sup> duration, *VT*, are subsequently expressed by Equations (6) and (7). Additionally, *VT* has to be transformed into the base reference observation period time *VT*0 , such that it is the same as in the definition of *GF*, as per Equation (8). Substituting *Vτ* and *VT*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}{\pi} \int\_{\frac{T-\tau}{2}}^{\frac{T+\tau}{2}} V(t)dt = V\_{0} \Big(1 - K \frac{G(\tau)}{\overline{\tau}}\Big), \text{ where} \\ G(\tau) = \int\_{\frac{T-\tau}{2}}^{\frac{T+\tau}{2}} g(t)dt = -\frac{T}{\pi} \Big(\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}\Big) \end{array} \tag{6}$$

$$V\_T = \frac{1}{T} \int\_0^T V(t)dt = V\_0 \left(1 + \frac{8}{15\pi}K\right) \tag{7}$$

$$V\_{T\_0} = \frac{V\_0(T\_0 - T) + \overline{V}\_T T}{T\_0} = V\_0 \left(1 + K \frac{8}{15\pi} \frac{T}{T\_0}\right) \tag{8}$$

$$=-\frac{GF(\tau, T\_0) - 1}{\frac{G(\tau)}{\tau} + \frac{8}{15\pi} \frac{T}{T\_0} GF(\tau, T\_0)}\tag{9}$$

After all the required parameters for Equation (5) were determined, the extreme gus<sup>t</sup> transient profile during typhoons in the Zhangbin area was illustrated, as in Figure 8. The related wind speeds following this extreme gus<sup>t</sup> profile are summarized in Table 1, 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 gus<sup>t</sup> 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.

*K*

**Figure 8.** Time history of wind speed in a gus<sup>t</sup> period.

**Table 1.** Wind speed parameters.


#### **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. 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 × 10<sup>6</sup> and 1.6 × 107, 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*x,*F*y,*F*z) denote the resultant forces acting on the wind turbine along three coordinate directions, (*Q*x,*Q*y,*Q*z) denote the resultant moments with respect to three coordinate axes, *Q*A denotes the rotor moment with respect to the rotor axis, and (*Q*p1,*Q*p2,*Q*p3) 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 *θ<sup>a</sup>*. 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*A and pitch moment of the blade (*Q*p1,*Q*p2,*Q*p3) are expected to be similar to a sine function. The peak loads at certain azimuthal angles are then obtained in the steady simulations.

$$
\Delta \mathfrak{a} = \tan^{-1} \left( \frac{\sin \theta\_l \sin \theta\_d}{\cos \theta\_l} \right) \tag{10}
$$

**Table 2.** Main particulars of the target wind turbine.

**Figure 9.** National Renewable Energy Laboratory (NREL) 5-MW wind turbine. (**a**) Tilt angle; (**b**) blade pitch angle; (**c**) aerodynamic loads.

**Figure 10.** Change of angle of attack due to tilt in a blade rotation. (**a**) Rotor view; (**b**) section view.

#### *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 generalpurpose 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 *ui* is the mean velocity component in the *xi* direction, *p* is the pressure, *ρ* and *μ* are the density and viscosity of air, respectively, *<sup>ρ</sup>uiuj* is the Reynolds stress, and *ρbi* is the gravitational force component in the *xi* direction. The SST *k*-*<sup>ω</sup>* 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, *ω* is the dissipation rate of *k*, *τij* is the shear stress, *μt* is the turbulent viscosity, and (*β*<sup>∗</sup>, *σk*, *γ*, *β*, *σω*, *δ*, *σω*2) are the equation constants.

$$\frac{\partial u\_i}{\partial x\_i} = 0 \tag{11}$$

$$\frac{\partial u\_i}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_j} (u\_i u\_j) = -\frac{1}{\rho} \frac{\partial p}{\partial \mathbf{x}\_i} + \frac{1}{\rho} \frac{\partial}{\partial \mathbf{x}\_j} \left[ \mu \left( \frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{\partial u\_j}{\partial \mathbf{x}\_i} \right) \right] - \frac{\partial \left( \overline{u\_i' u\_j'} \right)}{\partial \mathbf{x}\_j} + b\_i \tag{12}$$

$$\frac{\partial k}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_{\dot{j}}} (\mu\_{\dot{j}} k) = \frac{1}{\rho} \tau\_{\dot{l}\dot{j}} \frac{\partial \mu\_{\dot{i}}}{\partial \mathbf{x}\_{\dot{j}}} - \beta^\* \omega k + \frac{1}{\rho} \frac{\partial}{\partial \mathbf{x}\_{\dot{j}}} \left[ (\mu + \sigma\_{\dot{k}} \mu\_{t}) \frac{\partial k}{\partial \mathbf{x}\_{\dot{j}}} \right] \tag{13}$$

$$\frac{\partial \overline{\boldsymbol{\omega}}}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_j} (\boldsymbol{\mu}\_l \boldsymbol{\omega}) = \frac{\gamma}{\mu\_t} \boldsymbol{\tau}\_{lj} \frac{\partial \boldsymbol{u}\_i}{\partial \mathbf{x}\_j} - \beta \boldsymbol{\omega}^2 + \frac{1}{\rho} \frac{\partial}{\partial \mathbf{x}\_j} \left[ (\mu + \sigma\_\omega \mu\_t) \frac{\partial \boldsymbol{\omega}}{\partial \mathbf{x}\_j} \right] + 2(1 - \delta) \sigma\_\omega \frac{1}{\omega} \frac{\partial \boldsymbol{k}}{\partial \mathbf{x}\_j} \frac{\partial \boldsymbol{\omega}}{\partial \mathbf{x}\_j} \tag{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. 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 gus<sup>t</sup> 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.

**Figure 11.** (**a**) Computational domain; (**b**) boundary conditions.

**Table 3.** Geometrical parameters of computational domain.


**Table 4.** Boundary conditions.


#### *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 *α* = 0.3 was used. The secondorder upwind scheme was used for spatial discretization. Since the real geometry of the wind turbine is considered, a multi-block unstructured body-fitted grid arrangemen<sup>t</sup> 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*y about the tower base was selected as the target function in terms of grid number. Figure 13a displays the dependence of *Q*y on grid number, where *Q*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%.

**Figure 12.** Body-fitted grid for steady wind case. (**a**) Refined cylindrical region; (**b**) surface grid.

**Figure 13.** (**a**) Dependence of aerodynamic load on grid number; (**b**) discretization error.

The proposed numerical model is then validated for the rotor torque *Q*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*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*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.

**Figure 14.** (**a**) Rotor torque at rated condition; (**b**) axial velocity contour and in-plane vector on the midplane between the rotor plane and the tower, and pressure contour on the blade surface for the blade configuration of 0◦–120◦–240◦.

#### **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 gus<sup>t</sup> in transient mode. To enhance the numerical stability of the transient simulations, a precursor calculation of constant wind speed was first conducted. Then, one gus<sup>t</sup> 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 gus<sup>t</sup> simulation.

#### *4.1. Steady Aerodynamic Loadings*

The simulation results of steady extreme wind speed are summarized in Table 5. The drag force *F*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* 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*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◦ ≤ *θa* ≤ 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◦ ≤ *θa* ≤ 360◦. As expected, the blade pitch moment at *θ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*y is much higher than its counterpart *Q*x that constantly changes its direction in a rotor rotation period. Additionally, *Q*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*y. The twisting moment *Q*z is mostly governed by the rotor but its magnitude is about one-order of magnitude smaller than *Q*x and approximately two-order of magnitude less than *Q*y.

## *4.2. Dynamic Amplification*

The gus<sup>t</sup> simulation results for the blade configuration of 0◦–120◦–240◦ are shown in Table 6, which lists the aerodynamic loading at the time of peak wind speed, i.e., *t* =3s, 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.


**Table 5.** Wind loads under steady extreme wind.

**Table 6.** Wind loads under unsteady extreme gus<sup>t</sup> for the blade configuration of 0◦–120◦–240◦.


In order to investigate the correlation between the aerodynamic loads and wind speed in the gus<sup>t</sup> cases, the amplification factor, *DAF*, is defined as follows [12]:

$$DAF = \frac{\phi}{|\phi|\_{\text{max}}} \tag{15}$$

where *φ* denotes any physical quantity, and |*φ*|max is the maximum absolute value of the physical quantity *φ*. 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 gus<sup>t</sup> 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*A behaves similar to the rotor drag force and *Q*A precedes the gus<sup>t</sup> 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*y and *Q*x show a leading phase of the gus<sup>t</sup> peak while the phase of *Q*z lags the gus<sup>t</sup> peak by 0.16 s. The phase inconsistency between the maximum aerodynamic loadings and the gus<sup>t</sup> 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 gus<sup>t</sup> period where the superscript \* represents the physical quantity normalized by its maximum value in a gus<sup>t</sup> period. Figure 17 indicates that the wind velocity just experiences a sharp deceleration process as the flow velocity slowly reaches the gus<sup>t</sup> 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 gus<sup>t</sup> 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 summarizes the amplification factor of aerodynamic loads given in Tables 5 and 6. The ratio of the peak tower drag in gus<sup>t</sup> to the steady tower drag is approximately proportional to the square of the gus<sup>t</sup> factor, *GF*, but the corresponding ratio for rotor drag varies quite linearly with the gus<sup>t</sup> 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*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 *GF*<sup>2</sup> (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 *λ*-shaped). The gus<sup>t</sup> exhibits a positive impact on *Q*x with a peak reduction by 64% in a *λ*-shaped blade configuration. Contrary to *Q*x, the extreme moments along other two axes, i.e., *Q*y and *Q*z, are governed by the shear effect where the corresponding amplification factor is close to the gus<sup>t</sup> factor.

**Figure 16.** Time history of wind load for the blade configuration of 0◦–120◦–240◦: (**a**) *F*x; (**b**) *Q*p; (**c**) *Q*y.

**Figure 17.** Time history of amplification factor for *F*x.

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

