Investigations on the Unsteady Aerodynamic Characteristics of a Horizontal-Axis Wind Turbine during Dynamic Yaw Processes

Wind turbines inevitably experience yawed flows, resulting in fluctuations of the angle of attack (AOA) of airfoils, which can considerably impact the aerodynamic characteristics of the turbine blades. In this paper, a horizontal-axis wind turbine (HAWT) was modeled using a structured grid with multiple blocks. Then, the aerodynamic characteristics of the wind turbine were investigated under static and dynamic yawed conditions using the Unsteady Reynolds Averaged Navier-Stokes (URANS) method. In addition, start-stop yawing rotations at two different velocities were studied. The results suggest that AOA fluctuation under yawing conditions is caused by two separate effects: blade advancing & retreating and upwind & downwind yawing. At a positive yaw angle, the blade advancing & retreating effect causes a maximum AOA at an azimuth angle of 0◦. Moreover, the effect is more dominant in inboard airfoils compared to outboard airfoils. The upwind & downwind yawing effect occurs when the wind turbine experiences dynamic yawing motion. The effect increases the AOA when the blade is yawing upwind and vice versa. The phenomena become more dominant with the increase of yawing rate. The torque of the blade in the forward yawing condition is much higher than in backward yawing, owing to the reversal of the yaw velocity.


Introduction
The use of wind power has grown rapidly over the past few decades. The aerodynamics of the wind turbine are a core subject of wind power generation. Owing to the intermittency of wind, the orientation of a wind turbine must be frequently adjusted, typically using the yaw mechanism, to ensure the rotor disk is perpendicular to the direction of the wind. Yaw control is one of the most important strategies for improving the efficiency of wind energy conversion. Under the yawed condition, velocity component of the upstream wind is parallel to the plane of rotation. Thus, the angle of attack of airfoil along blade spanwise exhibits periodic variation in one revolution [1]. The aerodynamic performance of a wind turbine under yaw, as well as the surrounding flow, show highly unsteady characteristics, resulting in high fatigue loads that can dramatically reduce the operational life of a wind turbine. Therefore, understanding the unsteady aerodynamic characteristics of wind turbines under yaw is important for improving their design [2]. In particular, accurate predictions of unsteady aerodynamic loads under yawed conditions are important for blade design and operational control.
Knowledge of complex three-dimensional (3D) flow phenomena over wind turbine blades and near wake under the yawing condition are still lacking. Thus, improved numerical models that Yaw is a continuous rotation process with a constant or variable velocity of rotation. Qiu et al. [21] investigated the shaft torque of the blade and instantaneous wake flow under the dynamic yawing process with a yaw rate between 5 and 20 • /s for the NREL Phase VI wind turbine. For larger wind turbines, such as the NREL 5-MW et al., these yaw rate seem to be impossible owing to the larger moment of inertia of the rotor-nacelle subsystem due to enlargement of the blade and tower length. Leble et al. [35] carried out a series of CFD simulations to investigate the performance of the DTU 10-MW wind turbine and modeled the turbine as a structured grid with moving boundaries. Simulations were performed with a yaw rate of 0.3 • /s and 3 • yaw angle and the results showed that overall power under the dynamic yawing condition was much larger than for the yawed case. Then, Wen et al. [36] investigated the NREL 5-MW yawing dynamic sinuous motion with an averaged yaw rate of 1.2 • /s and 2.4 • /s for the case of f = 0.1 Hz and f = 0.2 Hz, respectively. The dynamic yaw motion induced the upwind & downwind yawing effect, which considerably influenced the AOA of blade sections. As the blade is yawing upwind, the AOA increases and vice versa. In summary, the unsteady aerodynamic characteristics of wind turbine blades under yaw considering full 3D rotational effects are still unclear. In the present paper, unsteady dynamic CFD simulations were performed to investigate the aerodynamic characteristics of a wind turbine blade during the rotational revolutions using a full 3D wind turbine model. With the assumption of the rigid body, the main contribution of this paper is the proposal of a new grid methodology for analyzing the overall performance using URANS simulation, aerodynamic loads, and flow field of wind turbines in the yawed and yawing processes. The method can be used to analyze interactions between transient aerodynamic phenomena associated with the wind turbine control system. The other main objective of this study was to investigate the effects of yaw on the dynamic output power, rotor thrust, and the blade sectional aerodynamic characteristics caused by continuous changes in the yaw angle. The main novelty of the current work is the inclusion of the dynamic yawing simulation of wind turbine with the multiple structured domain sliding mesh. The rest of this paper is organized as follows: Section 2 describes the geometry model and computational method. In Section 3, simulation results under several yaw angles are compared to experimental data to validate the model. The discussion mainly focuses on the sectional azimuthal variations of aerodynamic force coefficients, as well as the pressure distribution on the blade. Finally, the conclusions are summarized in Section 4.

Numerical Model
The 5-MW reference wind turbine (RWT) designed by the NREL was used in this study [37]. The blade uses the Delft University (DU) airfoil family with a relative thickness ranging from 21 to 40%. The blade chord, relative thickness, and twist angle all have non-linear distributions. The rotor diameter is 126 m and the wind turbine operates at a wind speed of 11.4 m/s with a rotational velocity of 12.1 rpm, resulting in the tip speed ratio of λ = 7. The definitions of yaw angle and azimuth angle are shown in Figure 1. Yaw is a continuous rotation process with a constant or variable velocity of rotation. Qiu et al. [21] investigated the shaft torque of the blade and instantaneous wake flow under the dynamic yawing process with a yaw rate between 5 and 20°/s for the NREL Phase VI wind turbine. For larger wind turbines, such as the NREL 5-MW et al., these yaw rate seem to be impossible owing to the larger moment of inertia of the rotor-nacelle subsystem due to enlargement of the blade and tower length. Leble et al. [35] carried out a series of CFD simulations to investigate the performance of the DTU 10-MW wind turbine and modeled the turbine as a structured grid with moving boundaries. Simulations were performed with a yaw rate of 0.3°/s and 3° yaw angle and the results showed that overall power under the dynamic yawing condition was much larger than for the yawed case. Then, Wen et al. [36] investigated the NREL 5-MW yawing dynamic sinuous motion with an averaged yaw rate of 1.2°/s and 2.4°/s for the case of f = 0.1 Hz and f = 0.2 Hz, respectively. The dynamic yaw motion induced the upwind & downwind yawing effect, which considerably influenced the AOA of blade sections. As the blade is yawing upwind, the AOA increases and vice versa. In summary, the unsteady aerodynamic characteristics of wind turbine blades under yaw considering full 3D rotational effects are still unclear. In the present paper, unsteady dynamic CFD simulations were performed to investigate the aerodynamic characteristics of a wind turbine blade during the rotational revolutions using a full 3D wind turbine model. With the assumption of the rigid body, the main contribution of this paper is the proposal of a new grid methodology for analyzing the overall performance using URANS simulation, aerodynamic loads, and flow field of wind turbines in the yawed and yawing processes. The method can be used to analyze interactions between transient aerodynamic phenomena associated with the wind turbine control system. The other main objective of this study was to investigate the effects of yaw on the dynamic output power, rotor thrust, and the blade sectional aerodynamic characteristics caused by continuous changes in the yaw angle. The main novelty of the current work is the inclusion of the dynamic yawing simulation of wind turbine with the multiple structured domain sliding mesh. The rest of this paper is organized as follows: Section 2 describes the geometry model and computational method. In Section 3, simulation results under several yaw angles are compared to experimental data to validate the model. The discussion mainly focuses on the sectional azimuthal variations of aerodynamic force coefficients, as well as the pressure distribution on the blade. Finally, the conclusions are summarized in Section 4.

Numerical Model
The 5-MW reference wind turbine (RWT) designed by the NREL was used in this study [37]. The blade uses the Delft University (DU) airfoil family with a relative thickness ranging from 21 to 40%. The blade chord, relative thickness, and twist angle all have non-linear distributions. The rotor diameter is 126 m and the wind turbine operates at a wind speed of 11.4 m/s with a rotational velocity of 12.1 rpm, resulting in the tip speed ratio of λ = 7. The definitions of yaw angle and azimuth angle are shown in Figure 1.

Computational Domain and Boundary Conditions
The computational domain, as shown in Figure 2a, is a rectangular zone that is 330D in length and 240D in width, where D is the rotor diameter. The computational domain is made up of two main zones, the far-field zone and internal zone, as depicted in Figure 2b. The far-field is used to set the inflow conditions. The internal zone can be further subdivided into two parts: (1) cylinder yaw zone with a radius of 6D and a height of 7.5D; (2) rotor roll zone with 1.5R radius and 1.6-m height. The cylinder zone was used to configure the dynamic yaw motion. The distance between the inlet and rotational cylinder was 120D. The rotor is located at the center of the rotor roll zone. The computational domain, as shown in Figure 2a, is a rectangular zone that is 330D in length and 240D in width, where D is the rotor diameter. The computational domain is made up of two main zones, the far-field zone and internal zone, as depicted in Figure 2b. The far-field is used to set the inflow conditions. The internal zone can be further subdivided into two parts: (1) cylinder yaw zone with a radius of 6D and a height of 7.5D; (2) rotor roll zone with 1.5R radius and 1.6-m height. The cylinder zone was used to configure the dynamic yaw motion. The distance between the inlet and rotational cylinder was 120D. The rotor is located at the center of the rotor roll zone.  Figure 3 shows the boundary conditions used in the simulations. A uniform wind speed of 11.4 m/s was set at the inlet of the domain with a turbulence intensity of 5%. Boundary conditions for the up and bottom planes in the domain are free-slip and no-slip walls. The blade surface was set as a no-slip wall. The pressure outlet condition was assigned to the outlet of the domain. To simplify the static yaw simulation, the rotational axis was fixed, while the inlet velocity direction was varied to generate different angles mimicking different yaw condition. In Fluent, it is possible to select a mesh that can process non-conformal interfaces, namely boundaries interfaces between cell zones, for which the mesh node locations are not identical. In the current simulation case, the matching mesh interface option was applied to three interfaces, for example, the interface between the far-field zone and yawing zone, which is more accurate than the mapped interface option.   Figure 3 shows the boundary conditions used in the simulations. A uniform wind speed of 11.4 m/s was set at the inlet of the domain with a turbulence intensity of 5%. Boundary conditions for the up and bottom planes in the domain are free-slip and no-slip walls. The blade surface was set as a no-slip wall. The pressure outlet condition was assigned to the outlet of the domain. To simplify the static yaw simulation, the rotational axis was fixed, while the inlet velocity direction was varied to generate different angles mimicking different yaw condition. In Fluent, it is possible to select a mesh that can process non-conformal interfaces, namely boundaries interfaces between cell zones, for which the mesh node locations are not identical. In the current simulation case, the matching mesh interface option was applied to three interfaces, for example, the interface between the far-field zone and yawing zone, which is more accurate than the mapped interface option. The computational domain, as shown in Figure 2a, is a rectangular zone that is 330D in length and 240D in width, where D is the rotor diameter. The computational domain is made up of two main zones, the far-field zone and internal zone, as depicted in Figure 2b. The far-field is used to set the inflow conditions. The internal zone can be further subdivided into two parts: (1) cylinder yaw zone with a radius of 6D and a height of 7.5D; (2) rotor roll zone with 1.5R radius and 1.6-m height. The cylinder zone was used to configure the dynamic yaw motion. The distance between the inlet and rotational cylinder was 120D. The rotor is located at the center of the rotor roll zone.  Figure 3 shows the boundary conditions used in the simulations. A uniform wind speed of 11.4 m/s was set at the inlet of the domain with a turbulence intensity of 5%. Boundary conditions for the up and bottom planes in the domain are free-slip and no-slip walls. The blade surface was set as a no-slip wall. The pressure outlet condition was assigned to the outlet of the domain. To simplify the static yaw simulation, the rotational axis was fixed, while the inlet velocity direction was varied to generate different angles mimicking different yaw condition. In Fluent, it is possible to select a mesh that can process non-conformal interfaces, namely boundaries interfaces between cell zones, for which the mesh node locations are not identical. In the current simulation case, the matching mesh interface option was applied to three interfaces, for example, the interface between the far-field zone and yawing zone, which is more accurate than the mapped interface option.

Grid Setup
Building multiple blocks is recommended for studying the yawing dynamics of coupled effects within two axes (yawing axis and rotational axis) since the relative rotation between zones can be investigated. Figure 4a,b show grids in the three different zones of the turbine (far-field, yaw cylinder, and rotation) generated using the ICEM software. The NREL 5-MW rotor was modeled including the hub (without the tower) [31]. The grid number for the far-field was 3.35 million, which is enough to transition from the inflow condition to the internal zone. Dynamic yaw motion was set, along with the rotational velocity in the cylinder yaw zone, with 3.25 million grids. Figure 4a illustrates the grids of the rotor zone, generated using AutoGrid tool in the NUMECA software. The grid number for this zone was 3.76 million, the first layer wall normal distance is about 10 −4 m. The total grid number was about 10.27 million. The work of Tran and Kim [38] employed a 6 million-cell grid for a 5-MW wind turbine. The blade surface was resolved with 93 cells along the span. In order to resolve the boundary layer, the grids around the blades were refined to keep y+ less than 5 on the blade surface, which satisfies the needs of the T-SST turbulent model. Building multiple blocks is recommended for studying the yawing dynamics of coupled effects within two axes (yawing axis and rotational axis) since the relative rotation between zones can be investigated. Figure 4a-b show grids in the three different zones of the turbine (far-field, yaw cylinder, and rotation) generated using the ICEM software. The NREL 5-MW rotor was modeled including the hub (without the tower) [31]. The grid number for the far-field was 3.35 million, which is enough to transition from the inflow condition to the internal zone. Dynamic yaw motion was set, along with the rotational velocity in the cylinder yaw zone, with 3.25 million grids. Figure 4a illustrates the grids of the rotor zone, generated using AutoGrid tool in the NUMECA software. The grid number for this zone was 3.76 million, the first layer wall normal distance is about 10 -4 m. The total grid number was about 10.27 million. The work of Tran and Kim [38] employed a 6 million-cell grid for a 5-MW wind turbine. The blade surface was resolved with 93 cells along the span. In order to resolve the boundary layer, the grids around the blades were refined to keep y+ less than 5 on the blade surface, which satisfies the needs of the T-SST turbulent model.

Numerical Methods
The unsteady Reynolds-averaged Navier-Stokes equations (URANS) were solved in ANSYS Fluent. The k-ω transitional SST turbulent model was used for turbulence modeling considering both laminar and transitional effects at the blade surface. The sliding mesh technique was used for data exchange between the rotational zone and the stationary zone. To ensure time-accurate calculations, a dual-time implicit time integration algorithm based on linearized second-order Euler backward differencing was used. The time-step size was equivalent to a rotor azimuthal increment of 5° coupled with 20 pseudo-time sub-iterations. The residual of the continuity equation was reduced by at least three orders throughout the calculations. For comparison, BEM computations using the Fast software were also performed. Work of Tran and Kim [38] employed 2° increments of the azimuth angle. The yawed case was achieved by changing the inlet velocity component to implement yaw inflow.
The yaw velocity can be constant, shown as the black dashed line in Figure 5. In fact, yawing is a process that changes gradually, similar to the definition of 2-s (or 4-s) start or stop duration. The selected variation of angular velocity for studying the influence of start-stop duration on the dynamic aerodynamic characteristics under yawing is shown a solid black line in Figure 6, which can be separated into three part: 1. Start duration with sinusoidal variation law of yaw velocity, 2. the stage with constant yaw velocity of 0.3°/s, 3. Stop duration with sinusoidal law of yaw velocity.

Numerical Methods
The unsteady Reynolds-averaged Navier-Stokes equations (URANS) were solved in ANSYS Fluent. The k-ω transitional SST turbulent model was used for turbulence modeling considering both laminar and transitional effects at the blade surface. The sliding mesh technique was used for data exchange between the rotational zone and the stationary zone. To ensure time-accurate calculations, a dual-time implicit time integration algorithm based on linearized second-order Euler backward differencing was used. The time-step size was equivalent to a rotor azimuthal increment of 5 • coupled with 20 pseudo-time sub-iterations. The residual of the continuity equation was reduced by at least three orders throughout the calculations. For comparison, BEM computations using the Fast software were also performed. Work of Tran and Kim [38] employed 2 • increments of the azimuth angle. The yawed case was achieved by changing the inlet velocity component to implement yaw inflow.
The yaw velocity can be constant, shown as the black dashed line in Figure 5. In fact, yawing is a process that changes gradually, similar to the definition of 2-s (or 4-s) start or stop duration. The selected variation of angular velocity for studying the influence of start-stop duration on the dynamic aerodynamic characteristics under yawing is shown a solid black line in Figure 6, which can be separated into three part: 1. Start duration with sinusoidal variation law of yaw velocity, 2. the stage with constant yaw velocity of 0.3 • /s, 3. Stop duration with sinusoidal law of yaw velocity. In the present investigation, the yawing process of the NREL 5-MW turbine was calculated under two different start-stop durations, 2-s and 4-s. Variation of the yaw velocity is shown in Figure  6, and changes more gradually with the 4-s duration of the start-stop process compared to the 2-s duration. During the start-stop yawing process, the sinusoidal variation of yaw velocity was applied to achieve transient yawing of the wind turbine, as follows: where rpm = 0.3°/s, according to the description of the NREL5 WM design manual [37], and duration is either 2-s or 4-s. Up to 0.3°/s, the wind turbine changes from transient to a yaw angle of approximately 20° with the maximum rpm.

Analysis of Results
The NREL 5-MW reference wind turbine (RWT) was studied [28]. The yawed and yawing wind turbines were considered. The parameters used in the simulation are listed in Table 1 for each case. Two of the cases assumed uniform inflow and the blades were assumed to be rigid in the model. Physically, the flow-field around a rotating HAWT is significantly influenced by the existence of wind shear, turbulence, gust, and yaw motion of the nacelle. For a yawing wind turbine, flow characteristics are more complex than those of a static yawed wind turbine, and the additional wind contribution effects transmitted to the rotor due to the nacelle motion must be considered. Therefore, accurate prediction of the unsteady aerodynamics calculated by many conventional numerical approaches is still questionable for a yawing wind turbine. In this study, unsteady CFD simulations In the present investigation, the yawing process of the NREL 5-MW turbine was calculated under two different start-stop durations, 2-s and 4-s. Variation of the yaw velocity is shown in Figure 6, and changes more gradually with the 4-s duration of the start-stop process compared to the 2-s duration. During the start-stop yawing process, the sinusoidal variation of yaw velocity was applied to achieve transient yawing of the wind turbine, as follows: where rpm = 0.3 • /s, according to the description of the NREL5 WM design manual [37], and duration is either 2-s or 4-s. Up to 0.3 • /s, the wind turbine changes from transient to a yaw angle of approximately 20 • with the maximum rpm. In the present investigation, the yawing process of the NREL 5-MW turbine was calculated under two different start-stop durations, 2-s and 4-s. Variation of the yaw velocity is shown in Figure  6, and changes more gradually with the 4-s duration of the start-stop process compared to the 2-s duration. During the start-stop yawing process, the sinusoidal variation of yaw velocity was applied to achieve transient yawing of the wind turbine, as follows: where rpm = 0.3°/s, according to the description of the NREL5 WM design manual [37], and duration is either 2-s or 4-s. Up to 0.3°/s, the wind turbine changes from transient to a yaw angle of approximately 20° with the maximum rpm.

Analysis of Results
The NREL 5-MW reference wind turbine (RWT) was studied [28]. The yawed and yawing wind turbines were considered. The parameters used in the simulation are listed in Table 1 for each case. Two of the cases assumed uniform inflow and the blades were assumed to be rigid in the model. Physically, the flow-field around a rotating HAWT is significantly influenced by the existence of wind shear, turbulence, gust, and yaw motion of the nacelle. For a yawing wind turbine, flow characteristics are more complex than those of a static yawed wind turbine, and the additional wind contribution effects transmitted to the rotor due to the nacelle motion must be considered. Therefore, accurate prediction of the unsteady aerodynamics calculated by many conventional numerical approaches is still questionable for a yawing wind turbine. In this study, unsteady CFD simulations

Analysis of Results
The NREL 5-MW reference wind turbine (RWT) was studied [28]. The yawed and yawing wind turbines were considered. The parameters used in the simulation are listed in Table 1 for each case. Two of the cases assumed uniform inflow and the blades were assumed to be rigid in the model. Physically, the flow-field around a rotating HAWT is significantly influenced by the existence of wind shear, turbulence, gust, and yaw motion of the nacelle. For a yawing wind turbine, flow characteristics are more complex than those of a static yawed wind turbine, and the additional wind contribution effects transmitted to the rotor due to the nacelle motion must be considered. Therefore, accurate prediction of the unsteady aerodynamics calculated by many conventional numerical approaches is still questionable for a yawing wind turbine. In this study, unsteady CFD simulations based on the sliding structure mesh technique were performed to analyze the yawing motion of the wind turbine due to the nacelle motion. Thus, to investigate the effects of vortex-wake-blade interactions on the aerodynamic performance of the wind turbine, the yawing motion of the rotating turbine blades due to the nacelle motion was considered.
The CFD method was first verification with grid independence, turbulent model studies and time step studies. Then the computational results are validated using both fixed yaw and yawing of the NREL 5-MW wind turbine. The unsteady aerodynamic loads of the yawing wind turbine were more sensitive to change as the yaw angle was varied. Nine simulation conditions were performed based on a fixed yaw angle of 0 • , 5 • , 10 • , 15 • , 20 • , 25 • , and 30 • .
Next, the total aerodynamic performance, which was used to validate the results of the simulation, will be discussed. Thereafter, aerodynamic characteristics and three-dimensional stall characteristics at different span sections will be investigated. Finally, the results of two separate simulations of the yawing start-stop stage will be used to examine the mechanism of the yawing effect. Table 2 shows the different torques calculated for five different numbers of cells. Every mesh was refined along the airfoil circumferential direction. The grid numbers were computed in the T-SST turbulent mode and simulations were performed with the following grid numbers: 1.6 million, 3.09 million, 3.67 million, 5.31 million, and 6.05 million. In addition, two grid numbers (1.6 million and 3.67 million) were investigated under different yaw angles (10 • and 20 • ), as shown in Table 2. For non-yaw condition, the design torque at 11.4 m/s is about 4.08 × 10 6 N·m, which is taken as a reference value. The relative errors for the results using other four mesh are 1.9~0.49%. The relative error for mesh of 3.67 million cells is smaller than 1% and the computational cost using this mesh is moderate. Thus this mesh was selected for further investigation.

Turbulent Model Studies
In this section, we present results of the unsteady independence study based on four turbulent models: SST, k-kl-ω, Reynold stress, and T-SST. The rotor torque under axial flow obtained for each turbulent model is presented in Table 3. The computed torque was similar for each of the four turbulent models. The simulation cases using the SST model and Reynold stress model were used to interpret the fully developed flow mechanism without transient processes, resulting less torque compared to using Transient SST turbulent model. When the wind turbine operates under low wind speeds, transient phenomena occur on the suction side of the blade; therefore, the T-SST turbulent model was selected for static and yawing process studies.

Time Step Studies
In the unsteady numerical simulation, the time step could influence the overall performance of the model, therefore, it is necessary to investigate time-step independence. Table 4 lists the averaged rotor torques during one revolution for non-yawed case obtained using the T-SST turbulent model with three different time steps (72, 144, and 360). The comparisons show slight differences among the results (less than 3%). The result using 72 time steps in one revolution gave a better agreement comparing to the designed value. Considering the computational cost, a time step of 72 during one revolution was selected for subsequent simulations.

The Validation of Numerical Simulation Results for Yawed Wind Turbine
After the verification of the simulation, the subsections below will show some aerodynamic analysis of wind turbine under yawed and yawing simulation. Figure 7 shows the overall performance of the rotor under yaw. As shown in Figure 7a, only a small amount of variation in rotor power can be observed under small yaw angles (≤5 • ). When the yaw angle exceeds 5 • , the rotor torque significantly decreases as the yaw angle increases. The torque value computed using the BEM method, with or without the Beddoes stall model, is higher than the value computed by CFD since the BEM computation does not consider flow separation and flow transient phenomena. Zhu [39] extensively investigated the combined effects of rotational augmentation and dynamic stall, and found that the hysteresis loop of aerodynamic load is much larger compared to 2D simulations. In fact, for the NREL 5-MW wind turbine and a wind speed of 11.4 m/s, flow separation and 3D radial flow mainly occur on the inner board region, resulting in a lower torque and lower thrust. The deviation in the power and thrust between CFD and BEM was 4.1% and 5.91%, respectively. Compared with the axial free inflow, the three functions on cos(γ), cos 2 (γ), cos 3 (γ) of the power and thrust are shown with dashed, dotted and dash-dotted line, respectively. The variation of averaged power in yaw conditions will decrease by cos 2 (γ); the averaged thrust agrees well with cos(γ).

Overall Performance Analysis for Different Yaw Angles
Energies 2019, 12, x FOR PEER REVIEW 9 of 23

Time Step Studies
In the unsteady numerical simulation, the time step could influence the overall performance of the model, therefore, it is necessary to investigate time-step independence. Table 4 lists the averaged rotor torques during one revolution for non-yawed case obtained using the T-SST turbulent model with three different time steps (72, 144, and 360). The comparisons show slight differences among the results (less than 3%). The result using 72 time steps in one revolution gave a better agreement comparing to the designed value. Considering the computational cost, a time step of 72 during one revolution was selected for subsequent simulations.

The Validation of Numerical Simulation Results for Yawed Wind Turbine
After the verification of the simulation, the subsections below will show some aerodynamic analysis of wind turbine under yawed and yawing simulation.
Overall Performance Analysis for Different Yaw Angles Figure 7 shows the overall performance of the rotor under yaw. As shown in Figure 7a, only a small amount of variation in rotor power can be observed under small yaw angles (≤5°). When the yaw angle exceeds 5°, the rotor torque significantly decreases as the yaw angle increases. The torque value computed using the BEM method, with or without the Beddoes stall model, is higher than the value computed by CFD since the BEM computation does not consider flow separation and flow transient phenomena. Zhu [39] extensively investigated the combined effects of rotational augmentation and dynamic stall, and found that the hysteresis loop of aerodynamic load is much larger compared to 2D simulations. In fact, for the NREL 5-MW wind turbine and a wind speed of 11.4 m/s, flow separation and 3D radial flow mainly occur on the inner board region, resulting in a lower torque and lower thrust. The deviation in the power and thrust between CFD and BEM was 4.1% and 5.91%, respectively. Compared with the axial free inflow, the three functions on cos(γ), cos 2 (γ), cos 3 (γ) of the power and thrust are shown with dashed, dotted and dash-dotted line, respectively. The variation of averaged power in yaw conditions will decrease by cos 2 (γ); the averaged thrust agrees well with cos(γ).

Aerodynamic Load Analysis along the Span of the Blade
Theoretical analysis is necessary to determine the mechanism behind the yawed effect. Figure 8 shows the velocity diagram under the yawed condition. A velocity component exists in the blade revolution plane and can be projected into radial and chordwise components denoted V r and V c , as shown in Figure 8a. The retreating and advancing effects occur because of V c , which causes the sectional load and AOA to fluctuate periodically. The maximum AOA and maximum load occur when the blade is in the 12 o'clock direction in the current simulation setup. Additionally, loads on the retreating side are much larger than on the advancing side. Figure 9b,c quantify the influence of V c , which is calculated by: Energies 2019, 12, x FOR PEER REVIEW 10 of 23

Aerodynamic Load Analysis along the Span of the Blade
Theoretical analysis is necessary to determine the mechanism behind the yawed effect. Figure 8 shows the velocity diagram under the yawed condition. A velocity component exists in the blade revolution plane and can be projected into radial and chordwise components denoted Vr and Vc, as shown in Figure 8a. The retreating and advancing effects occur because of Vc, which causes the sectional load and AOA to fluctuate periodically. The maximum AOA and maximum load occur when the blade is in the 12 o'clock direction in the current simulation setup. Additionally, loads on the retreating side are much larger than on the advancing side. Figure 9b,c quantify the influence of Vc, which is calculated by:   (1) Time-averaged load analysis on spanwise section under yawed condition Figures 9-11 show the analysis of the variation of averaged AOA, C n , and C t under one rotor revolution. The value of C n and C t is calculated by: where C p-up and C p-down denote the pressure coefficients in the upper and down side of the airfoil, respectively; y0 and y1 mean the leading edge position and trailing edge position of the airfoil. The AOA distribution presented in Figure 9 takes into account the advancing and retreating effects, as previously shown by Castellani [20] in the yawed simulation of a HAWT using the BEM model. Comparing the simulation results for CFD and FAST, the same trends can be observed along the blade span, whereas the AOA value is much larger than the FAST result. Figures 10 and 11 show the variation of aerodynamic load spanwise along the blade. For the axial flow, some large abnormal fluctuations in the AOA and aerodynamic loads can be observe d. Differences in the results of the CFD and FAST methods mainly occur along the inner board, owing to flow separation, which lead to smaller aerodynamic loads compared to those computed by the BEM model that does not consider flow separation. The AOA distribution presented in Figure 9 takes into account the advancing and retreating effects, as previously shown by Castellani [20] in the yawed simulation of a HAWT using the BEM model. Comparing the simulation results for CFD and FAST, the same trends can be observed along the blade span, whereas the AOA value is much larger than the FAST result. Figures 10 and 11 show the variation of aerodynamic load spanwise along the blade. For the axial flow, some large abnormal fluctuations in the AOA and aerodynamic loads can be observe d. Differences in the results of the CFD and FAST methods mainly occur along the inner board, owing to flow separation, which lead to smaller aerodynamic loads compared to those computed by the BEM model that does not consider flow separation.  The AOA distribution presented in Figure 9 takes into account the advancing and retreating effects, as previously shown by Castellani [20] in the yawed simulation of a HAWT using the BEM model. Comparing the simulation results for CFD and FAST, the same trends can be observed along the blade span, whereas the AOA value is much larger than the FAST result. Figures 10 and 11 show the variation of aerodynamic load spanwise along the blade. For the axial flow, some large abnormal fluctuations in the AOA and aerodynamic loads can be observe d. Differences in the results of the CFD and FAST methods mainly occur along the inner board, owing to flow separation, which lead to smaller aerodynamic loads compared to those computed by the BEM model that does not consider flow separation.  The AOA distribution presented in Figure 9 takes into account the advancing and retreating effects, as previously shown by Castellani [20] in the yawed simulation of a HAWT using the BEM model. Comparing the simulation results for CFD and FAST, the same trends can be observed along the blade span, whereas the AOA value is much larger than the FAST result. Figures 10 and 11 show the variation of aerodynamic load spanwise along the blade. For the axial flow, some large abnormal fluctuations in the AOA and aerodynamic loads can be observe d. Differences in the results of the CFD and FAST methods mainly occur along the inner board, owing to flow separation, which lead to smaller aerodynamic loads compared to those computed by the BEM model that does not consider flow separation. advancing and retreating effect, which is partly caused by the inflow velocity component under yaw. After post-processing, the AOA using the combination of the CFD sectional airfoil aerodynamic loads and the BEM method, the maximum and minimum AOA during one revolution occur at an azimuth angle of 0 • and 180 • , respectively. Wen [36] found that under the impact of non-uniform effects (due to variations of the induced factor caused by the radial position and azimuth angle in the rotational plane), the maximum AOA tends to occur at 90 • for inboard airfoils. The maximum and minimum aerodynamic loads occur at an azimuth angle of 90 • and 270 • , respectively. Thus, future computations of the AOA should consider non-uniform effects.
Energies 2019, 12, x FOR PEER REVIEW 12 of 23 processing, the AOA using the combination of the CFD sectional airfoil aerodynamic loads and the BEM method, the maximum and minimum AOA during one revolution occur at an azimuth angle of 0° and 180°, respectively. Wen [36] found that under the impact of non-uniform effects (due to variations of the induced factor caused by the radial position and azimuth angle in the rotational plane), the maximum AOA tends to occur at 90° for inboard airfoils. The maximum and minimum aerodynamic loads occur at an azimuth angle of 90° and 270°, respectively. Thus, future computations of the AOA should consider non-uniform effects.    processing, the AOA using the combination of the CFD sectional airfoil aerodynamic loads and the BEM method, the maximum and minimum AOA during one revolution occur at an azimuth angle of 0° and 180°, respectively. Wen [36] found that under the impact of non-uniform effects (due to variations of the induced factor caused by the radial position and azimuth angle in the rotational plane), the maximum AOA tends to occur at 90° for inboard airfoils. The maximum and minimum aerodynamic loads occur at an azimuth angle of 90° and 270°, respectively. Thus, future computations of the AOA should consider non-uniform effects.  Figure 14 shows the variation of the averaged AOA with aerodynamic load at five typical spanwise sections with respect to yaw angle. Fluctuations of the aerodynamic performance are clearly observed and become larger as the yaw angle increases. The mean aerodynamic performance remains relatively constant with less decrease due to the yaw effect. Variation of the AOA and aerodynamic loads exhibit much higher fluctuations in the  Figure 14 shows the variation of the averaged AOA with aerodynamic load at five typical spanwise sections with respect to yaw angle. Fluctuations of the aerodynamic performance are clearly observed and become larger as the yaw angle increases. processing, the AOA using the combination of the CFD sectional airfoil aerodynamic loads and the BEM method, the maximum and minimum AOA during one revolution occur at an azimuth angle of 0° and 180°, respectively. Wen [36] found that under the impact of non-uniform effects (due to variations of the induced factor caused by the radial position and azimuth angle in the rotational plane), the maximum AOA tends to occur at 90° for inboard airfoils. The maximum and minimum aerodynamic loads occur at an azimuth angle of 90° and 270°, respectively. Thus, future computations of the AOA should consider non-uniform effects.   compared to the middle and outer board. With the increasing of the yaw angle, the circumferential loads and AOA also cause higher maximum loads.

Numerical Simulation of Dynamic Yawing Wind Turbine
This section presents results of the numerical analysis of the NREL 5-MW RWT under the dynamic yawing process. Figure 15 shows the distribution of the wind rotor torque under different start-stop yaw velocities. During the yawing process, fluctuation of the wind rotor torque is small. At the beginning of dynamic yaw, larger torques occur due to changes in the rotor position. As the yaw angle increases, the torque gradually decreases and is 5 × 10 5 N·m larger in the case of the 2-s duration compared to the 4-s duration. The 2-s duration may induce a larger rotor torque since the yaw angle changes more quickly than with 4-s duration. The reasons are given below. At the start and stop stage, the yaw velocity is changed with the sinusoidal variation law related to the frequency, the higher frequency of yaw velocity caused much larger of the power and thrust, which is similar to the variation of power under platform pitching [38,40]. When the wind turbine is yawing with the constant yaw velocity, the additional velocity inducing by yawing is the same, and the power is only decreased with the square cosine of yaw angle. When wind turbine yawed to the stage of stop, the power under the 2-s case decrease much faster than under 4-s cases, which is similar to the yaw start period. The mean aerodynamic performance remains relatively constant with less decrease due to the yaw effect. Variation of the AOA and aerodynamic loads exhibit much higher fluctuations in the inner board compared to the middle and outer board. With the increasing of the yaw angle, the circumferential loads and AOA also cause higher maximum loads.

Numerical Simulation of Dynamic Yawing Wind Turbine
This section presents results of the numerical analysis of the NREL 5-MW RWT under the dynamic yawing process. Figure 15 shows the distribution of the wind rotor torque under different start-stop yaw velocities. During the yawing process, fluctuation of the wind rotor torque is small. At the beginning of dynamic yaw, larger torques occur due to changes in the rotor position. As the yaw angle increases, the torque gradually decreases and is 5 × 10 5 N·m larger in the case of the 2-s duration compared to the 4-s duration. The 2-s duration may induce a larger rotor torque since the yaw angle changes more quickly than with 4-s duration. The reasons are given below. At the start and stop stage, the yaw velocity is changed with the sinusoidal variation law related to the frequency, the higher frequency of yaw velocity caused much larger of the power and thrust, which is similar to the variation of power under platform pitching [38,40]. When the wind turbine is yawing with the constant yaw velocity, the additional velocity inducing by yawing is the same, and the power is only decreased with the square cosine of yaw angle. When wind turbine yawed to the stage of stop, the power under the 2-s case decrease much faster than under 4-s cases, which is similar to the yaw start period. The results of 58 rotor revolutions were analyzed. During each revolution, 12 torque measurements were collected. Then, the fast Fourier transform was used to obtain the rotor aerodynamic frequency. Figure 16 shows the torque power spectra of the rotor under 2-s and 4-s duration yawing process. The blade passing frequency of the NREL 5-MW turbine is about 0.2017 Hz (1P fluctuation), and is clearly the main frequency of the two yawing processes. Both are 0.6 Hz (which is approximately a 3P rotor fluctuation), which is similar to results presented by Castellani [20]. Figure16 a shows a secondary frequency of 0.2 Hz (1P fluctuation). Due to the less sampling data in the sinusoidal stage, the yawing start-stop frequency (2-s and 4-s duration with a corresponding main frequency of 1/8 Hz and 1/16 Hz, respectively) cannot be captured in the computation of the torque power structure. The results of 58 rotor revolutions were analyzed. During each revolution, 12 torque measurements were collected. Then, the fast Fourier transform was used to obtain the rotor aerodynamic frequency. Figure 16 shows the torque power spectra of the rotor under 2-s and 4-s duration yawing process. The blade passing frequency of the NREL 5-MW turbine is about 0.2017 Hz (1P fluctuation), and is clearly the main frequency of the two yawing processes. Both are 0.6 Hz (which is approximately a 3P rotor fluctuation), which is similar to results presented by Castellani [20]. Figure 16a shows a secondary frequency of 0.2 Hz (1P fluctuation). Due to the less sampling data in the sinusoidal stage, the yawing start-stop frequency (2-s and 4-s duration with a corresponding main frequency of 1/8 Hz and 1/16 Hz, respectively) cannot be captured in the computation of the torque power structure.  Figure 17 shows the variation of torque under different yaw rates in two periodic processes of start-stop yawing. Similarly, the torque of the blade under a 2-s duration is larger than the torque of the 4-s duration. Interestingly, the torque of the blade in the forward yaw stage (shown in Figure 18 within 0~66.8 min, yaw angle from 0° to 20°) is larger than that of the backward yaw stage (shown in Figure 18 within 66.8~133.6 min, yaw angle from 20° back to 0°) since the dynamic yawing effect, which generates the dynamic velocity, reduces the relative velocity.

Wake Flow Characteristics
Wake effects are important in analyzing wind turbine aerodynamics. For convenience, the dynamic yaw was classified into six cases, as illustrated in Figure 18. Cases 1,3,4,6 include the forward-start yaw, forward-stop yaw, backward-start yaw, and backward-stop yaw, respectively. For Cases 2 and 5, the yaw angle is 10° and they represent forward yaw and backward yaw, respectively. The grey line in the figure illustration the variation of yaw angular velocity during the simulation, and shows that except the Cases 2 and 5, other cases are all in the simulation about the variation of yaw angular velocity.  Figure 17 shows the variation of torque under different yaw rates in two periodic processes of start-stop yawing. Similarly, the torque of the blade under a 2-s duration is larger than the torque of the 4-s duration. Interestingly, the torque of the blade in the forward yaw stage (shown in Figure 18 within 0~66.8 min, yaw angle from 0 • to 20 • ) is larger than that of the backward yaw stage (shown in Figure 18 within 66.8~133.6 min, yaw angle from 20 • back to 0 • ) since the dynamic yawing effect, which generates the dynamic velocity, reduces the relative velocity.  Figure 17 shows the variation of torque under different yaw rates in two periodic processes of start-stop yawing. Similarly, the torque of the blade under a 2-s duration is larger than the torque of the 4-s duration. Interestingly, the torque of the blade in the forward yaw stage (shown in Figure 18 within 0~66.8 min, yaw angle from 0° to 20°) is larger than that of the backward yaw stage (shown in Figure 18 within 66.8~133.6 min, yaw angle from 20° back to 0°) since the dynamic yawing effect, which generates the dynamic velocity, reduces the relative velocity.

Wake Flow Characteristics
Wake effects are important in analyzing wind turbine aerodynamics. For convenience, the dynamic yaw was classified into six cases, as illustrated in Figure 18. Cases 1,3,4,6 include the forward-start yaw, forward-stop yaw, backward-start yaw, and backward-stop yaw, respectively. For Cases 2 and 5, the yaw angle is 10° and they represent forward yaw and backward yaw, respectively. The grey line in the figure illustration the variation of yaw angular velocity during the simulation, and shows that except the Cases 2 and 5, other cases are all in the simulation about the variation of yaw angular velocity.

Wake Flow Characteristics
Wake effects are important in analyzing wind turbine aerodynamics. For convenience, the dynamic yaw was classified into six cases, as illustrated in Figure 18. Cases 1,3,4,6 include the forward-start yaw, forward-stop yaw, backward-start yaw, and backward-stop yaw, respectively. For Cases 2 and 5, the yaw angle is 10 • and they represent forward yaw and backward yaw, respectively. The grey line in the figure illustration the variation of yaw angular velocity during the simulation, and shows that except the Cases 2 and 5, other cases are all in the simulation about the variation of yaw angular velocity.
(1) Forward yaw-start stage (Case 1). The velocity contours at t = 1/8T (T = 8 s or 16 s) at the beginning of dynamic yaw are shown in Figure 19. The near wake velocity flow structure shows the same status at the beginning of the yawing process for both the 2-s and 4-s start-stop duration.
The different yaw velocities have very little influence on the velocity distribution in the 2D range downstream. (1) Forward yaw-start stage (Case 1). The velocity contours at t = 1/8T (T = 8 s or 16 s) at the beginning of dynamic yaw are shown in Figure 19. The near wake velocity flow structure shows the same status at the beginning of the yawing process for both the 2-s and 4-s start-stop duration. The different yaw velocities have very little influence on the velocity distribution in the 2D range downstream.
(a) 2-s ( b) 4-s (2) Yaw angle of 10° under forward yawing stage (Case 2). Sketches of the velocity streamlines of the two dynamic yawing processes (yaw angle of 10°) are shown in Figure 20a and 20b are similar to those obtained for the yawed case, as shown in Figure 20c. The dynamic yaw rotates with a fixed yaw velocity of 0.3°/s, and similar results are observed for the dynamic process. More energy intermediate effects can be observed between the wake zone and the main flow zone than in the yawed case. This may be due to the effects of dynamic stall.  (3) Forward-yaw-stop stage (Case 3). Figure 21 shows the instantaneous velocity contours for the dynamic yawing and yawed cases. Both dynamic yawing processes result in a much larger wake     (3) Forward-yaw-stop stage (Case 3). Figure 21 shows the instantaneous velocity contours for the dynamic yawing and yawed cases. Both dynamic yawing processes result in a much larger wake zones than under the yawed case. Meanwhile, some deflection occurs in the velocity wake, which is similar to the wake deflection effect reported in the work of Qian [18], which took into account the velocity deficit and turbulent intensity using a Gaussian-based wake model.
zones than under the yawed case. Meanwhile, some deflection occurs in the velocity wake, which is similar to the wake deflection effect reported in the work of Qian [18], which took into account the velocity deficit and turbulent intensity using a Gaussian-based wake model. (4) Backward-yaw start stage (Case 4). Figure 22 illustrates the velocity wake in the backward-yaw start stage. The wake zone velocity field retains almost the same flow structure. To refine the simulation results, large eddy simulations can be used to improve the interpretation of the flow process.
(a) 2-s ( b) 4-s  (6) Backward-yaw stop stage (Case 6). As shown in Figure 24, the wind turbine returns to its initial state, the wind direction is normal to the rotor rotational plane, and the three-dimensional flow structure is symmetrical in streamwise. (4) Backward-yaw start stage (Case 4). Figure 22 illustrates the velocity wake in the backward-yaw start stage. The wake zone velocity field retains almost the same flow structure. To refine the simulation results, large eddy simulations can be used to improve the interpretation of the flow process. zones than under the yawed case. Meanwhile, some deflection occurs in the velocity wake, which is similar to the wake deflection effect reported in the work of Qian [18], which took into account the velocity deficit and turbulent intensity using a Gaussian-based wake model. (4) Backward-yaw start stage (Case 4). Figure 22 illustrates the velocity wake in the backward-yaw start stage. The wake zone velocity field retains almost the same flow structure. To refine the simulation results, large eddy simulations can be used to improve the interpretation of the flow process.
(a) 2-s ( b) 4-s  (6) Backward-yaw stop stage (Case 6). As shown in Figure 24, the wind turbine returns to its initial state, the wind direction is normal to the rotor rotational plane, and the three-dimensional flow structure is symmetrical in streamwise.  zones than under the yawed case. Meanwhile, some deflection occurs in the velocity wake, which is similar to the wake deflection effect reported in the work of Qian [18], which took into account the velocity deficit and turbulent intensity using a Gaussian-based wake model. (4) Backward-yaw start stage (Case 4). Figure 22 illustrates the velocity wake in the backward-yaw start stage. The wake zone velocity field retains almost the same flow structure. To refine the simulation results, large eddy simulations can be used to improve the interpretation of the flow process.
(a) 2-s ( b) 4-s  (6) Backward-yaw stop stage (Case 6). As shown in Figure 24, the wind turbine returns to its initial state, the wind direction is normal to the rotor rotational plane, and the three-dimensional flow structure is symmetrical in streamwise. (6) Backward-yaw stop stage (Case 6). As shown in Figure 24, the wind turbine returns to its initial state, the wind direction is normal to the rotor rotational plane, and the three-dimensional flow structure is symmetrical in streamwise. In summary, dynamic yawing based on two start-stop durations shows the approximate flow structure, and the upwind and downwind effect induced by yawing process expands the wake zone earlier and much larger than the yawed condition. In summary, dynamic yawing based on two start-stop durations shows the approximate flow structure, and the upwind and downwind effect induced by yawing process expands the wake zone earlier and much larger than the yawed condition.

Aerodynamic Characteristics along Blade Spanwise Section
To investigate the yawing effect, additional velocity induced by the yawing wind rotor should also be examined. Figure 25 illustrates the velocity triangle along the span of the blade under dynamic yawing. The dashed line of the yawing zone indicates positive yaw (negative yaw was not investigated in the present simulation). Four process variations can be extracted according to the dynamic yawing stage. Note that R = 0 means that initial position of the rotational axis of the wind turbine is the rotor hub position, which is different from the platform yawing process used for the offshore wind turbines. The direction of V dyn is different on both sides of the yaw axis. In the process of yawing counterwise, V dyn and the inflow wind speed create an acute angle in the right side of yaw axis, while V dyn and the inflow wind speed create an obtuse angle in the left one of yaw axis, as shown in Figure 26a. The process of yawing clockwise is just the contrary to the case of yawing counterwise. The absolute formulation of V dyn can be written as: where, V dyn become zero at the azimuth angle of 0 • and 180 • in the current setup of wind turbine. The relative velocity of the section under dynamic yawing can be defined as: where η = 1 if the direction of Vdyn and inflow wind velocity creates an acute angle and η = −1 when the direction of Vdyn and inflow wind velocity creates an blunt angle. In the present yawing The relative velocity of the section under dynamic yawing can be defined as: where η = 1 if the direction of V dyn and inflow wind velocity creates an acute angle and η = −1 when the direction of V dyn and inflow wind velocity creates an blunt angle. In the present yawing simulation, if the wind turbine is rotating in the positive yaw direction, β = 1. Figure 26 shows the variation of the AOA in three typical section (r/R = 0.21, 0.4, or 0.67), with respect to yaw angles 0-20-0-20-0 • . The upper abscissa is the wind rotor yaw angle, while the bottom abscissa is the wind turbine rotational cycle. The vertical coordinate indicates the AOA distribution along the radial direction of the blade. The three sections are indicated by solid lines of different thickness of airfoil. Inner board airfoils are the thickest and outer board one are thinnest.
The AOA gradually decreases along the span of the blade. Due to the effect of dynamic yawing, the AOA oscillates during the rotational period of the rotor, with the combined effects of both retreating & advancing and upwind & downwind. Figure 26a shows the calculated results of 2-s start-stop yaw rate. In this case, the start-stop velocity is faster, and the overall fluctuation are larger than those in the 4-s scenario. The start-stop effect mainly affects the AOA near the outer board but has less influence on other areas. In addition, the start-stop process leads to changes in the AOA.  Figure 27 show the distribution of the normal force coefficient at three typical sections (r/R = 0.21, 0.4, 0.67) within the dynamic process, including two forward yawing and backward yawing states with a yaw angle of 20°. Variation of Cn under dynamic yawing with a 2-s duration causes larger differences at the radial position (r/R = 0.4) than the 4-s duration. Under the yawed case, the average normal force coefficient of the r/R = 0.21 section with a yaw angle of 20° is 1.1, and the maximum and minimum load coefficients are 1.7 and 0.5, respectively (see Figure 14). In the current yawing case, the maximum and minimum load factors are 1.9 and 1.0, which are due to the downwind and upwind effects of the yaw dynamics, similar to the horizontal wind shear effect, but more pronounced than typical horizontal wind shear effects. In the start and stop duration of yawing, the yawing velocity of 2-s case has much higher frequency than 4-s case, resulting to much higher additional velocity and aerodynamic loads. Thus the overall performance of 2-s yawing presents much larger fluctuation in the process of yawing.  Figure 27 show the distribution of the normal force coefficient at three typical sections (r/R = 0.21, 0.4, 0.67) within the dynamic process, including two forward yawing and backward yawing states with a yaw angle of 20 • . Variation of C n under dynamic yawing with a 2-s duration causes larger differences at the radial position (r/R = 0.4) than the 4-s duration. Under the yawed case, the average normal force coefficient of the r/R = 0.21 section with a yaw angle of 20 • is 1.1, and the maximum and minimum load coefficients are 1.7 and 0.5, respectively (see Figure 14). In the current yawing case, the maximum and minimum load factors are 1.9 and 1.0, which are due to the downwind and upwind effects of the yaw dynamics, similar to the horizontal wind shear effect, but more pronounced than typical horizontal wind shear effects. In the start and stop duration of yawing, the yawing velocity of 2-s case has much higher frequency than 4-s case, resulting to much higher additional velocity and aerodynamic loads. Thus the overall performance of 2-s yawing presents much larger fluctuation in the process of yawing.

Conclusions
This paper has presented the results of unsteady numerical simulations investigating the static and dynamic aerodynamic performances of a NREL 5-MW HAWT under yawed and yawing conditions. The blades were considered rigid for cases with a prescribed motion, and the tower was not included in the computational domain. The numerical simulation results based on the T-SST turbulent model were consistent with the BEM data. Fluctuating and alternative loads were observed for the yawed case, resulting in the retreating and advancing effect. As the yawing angle decreased, overall performance also decreased according to the cosine law. The results suggest larger variations in power for the dynamic yaw case, compared to yawed cases. The wind turbine suffers from the coupled effects of the rotational axes (yawing axis and rotating axis). For maximum loading of rotor, the combined effects of the yawed condition and yawing resulted in an azimuth around the first periodic, which became larger as the yaw angle increased.
In the simulations of two different yawing start-stop durations (2-s and 4-s), the faster shift of yaw angle (2-s duration) resulted in a larger torque than the 4-s duration, owing to the additional velocity induced by the yawing effect of different sinusoidal frequency and yaw velocity magnitude. Interestingly, the torque of the rotor and blade was much higher in forward yawing compared to backward yawing, caused by reversing the direction of the yaw velocity. The wake deflection occurs in the near wake flow structure and expands more than in the yawed case. The different start-stop durations under dynamic yawing have little influence on the near wake velocity contour. The next step will be to include deformation of the blade in the simulation of multiple axial angular motions.
Author Contributions: S.K. and H.H. contributed to the conception of the study. X.W. and Z.Y. contributed significantly to analysis and manuscript preparation; Z.Y. performed the data analyses and wrote the manuscript; X.W. helped perform the analysis with constructive discussions.

Conflicts of Interest:
The authors declare no conflict of interest. Normal forces coefficient on the local airfoil C t Tangential forces coefficient on the local airfoil C p-up Pressure coefficient at the upper side of airfoil C p-down Pressure coefficient at the down side of airfoil C l Life forces coefficient C d

Abbreviations
Drag forces coefficient C X Tangent forces coefficient on the rotational plane C Z Axial forces coefficient on the rotational plane C p Pressure coefficient C po Power coefficent C t Thrust coefficent T Wind rotor torque F t Wind rotor thrust