A Numerical Study of Wind Turbine Airfoil Combined Oscillations Considering Phase Difference under Yaw Loads

 Abstract: In order to further improve the operating efficiency of wind turbines and explore the aerodynamic performance of the complex motion of wind turbine blades under yaw loads. In this study, the change in the angle of attack of the blade section airfoil under yaw load can be modeled as an oscillating airfoil and combined with the blade's flapwise motion. The NREL S809 airfoil are chosen for the research, based on the SST k-ω turbulence model with transition correction, under the condition of Reynolds number of 10 6 . The effect of phase difference on its aerodynamic performance under combined flapwise and pitching motion in various flapwise amplitudes and working conditions were analyzed. For the combined oscillations, the effects of the flapwise amplitude ( h ) in the range of 0.2≤ h ≤0.5 are investigated with the phase differences of Φ=±3π/4, ±π/2, ±π/4, 0. The results show that the phase difference between the pitching motion and the flapping motion and the different flapping amplitudes can have a large impact on the aerodynamic performance of the airfoil during dynamic stall, but the degree of influence is greatly different in different situations.


Introduction
With the rapid development of wind power technology in recent years, the size of wind turbines has gradually increased. Wind turbines are in complex wind fields all year round. When the wind turbine is subjected to yaw  [1]. These changes may cause dynamic stall (DS) of the airfoil [2][3].
Meanwhile, the blades also vibrate when they are rotating. These vibrations are mainly manifested in the bending vibration (edgewise) of the blade in the plane of rotation and the bending vibration (flapwise) perpendicular to the plane of rotation [4].
In the past few decades, the DS of the airfoil and the vibration of the blades has been studied by many researchers. The most studied is the DS phenomenon of the airfoil pitching motion. McCroskey [5] experimentally studied the DS phenomenon of airfoil pitching, and they found that the stall AOA of the airfoil in pitch motion is larger than that of the static airfoil, resulting in higher maximum values of lift and drag coefficients (CL and CD). However, compared with the static airfoil, the reattachment of the dynamic airfoil boundary layer will occur at a lower angle of incidence, resulting in a hysteresis of the aerodynamic coefficient. On the other hand, the periodic load due to DS and the unstable wind turbine wake have a negative impact on the rotor performance of the wind turbine [6][7], leading to fatigue failure [8][9]. In the 1990s, Ramsay et al. [10] systematically studied the aerodynamic performance of the airfoil under pitching oscillation motions through wind tunnel tests of the airfoil dedicated to wind turbines. Choudhry et al. [11] experimentally investigated the effects of reduced frequency (k), Mach number (Ma), Reynolds number (Re), and airfoil geometry on the DS of the airfoil. In addition, they analyzed the behavior of the lift curve slope of an airfoil subjected to continuous pitching DS to explore the mechanism of instability. Gharali et al. [12][13] studied the aerodynamic load of an airfoil when the flow incident velocity and pitch ·3· incident angle oscillated at the same frequency but with a certain phase difference. Subsequently, the effects of different reduced frequencies on the above problems were studied, and it was found that the reduced frequency and phase difference can significantly change the aerodynamic load during DS [14].
However, as the power of modern wind turbines gradually increases, the blades become thinner and lighter, so the blades become more flexible. Inevitably, the blades of large wind turbines will experience more severe vibrations during operation [15]. Considering the actual operating conditions of the wind turbine, Xiong et al. [16] used the computational fluid dynamics (CFD) method to study the DS characteristics of the S809 airfoil under the coupling of flapwise, edgewise and pitching. It is found that the airfoil has a non-negligible DS intensity under the flapwise motion, although it is weaker than that under the pitch motion. Kun et al. [17] studied in more detail the aerodynamic performance of wind turbine airfoil when flapwise, edgewise and coupling the two under different amplitudes.
When a wind turbine is subjected to yaw loads, the AOA of the blade section airfoil changes with the position of the blade. Under idealized yaw loads, the motion of the blade section airfoil can be approximated as pitch oscillation [1], but the blade itself also has vibrations in the flapping and edgewise direction during operation. Due to the different types of wind disturbances that cause yaw loads, the phase of the airfoil pitch oscillations and other vibrations may be different, resulting in changes in their aerodynamic performance. Existing literature has generally studied the pitching motion, flapwise and edgewise motion of wind turbine airfoils, or several coupled motions [18][19][20][21], but did not consider the phase of these motions during the superposition process. In this study, the aerodynamic performance of a wind turbine airfoil when the pitching motion and the flapwise motion at the same frequency but with different phase differences under different degrees of stall conditions are studied in detail. Figure 1 shows that the three different wind disturbances that cause the yaw load of a horizontal axis wind turbine. In the Figure 1 (a-c), the down direction of the x-axis is Angle of attack π π Angle of attack π π 2π

Yaw Loads
Angle of attack (a) Positive horizontal cross (b) Positive horizontal shear (c) Positive vertical wind Figure 1 Three types of wind disturbances causing yaw loads ·4· zero azimuth (Ψ=0), and the direction of rotation of the wind wheel is positive. Figure 1(a) shows that under positive horizontal cross wind Vy, when the blades are in π/2<Ψ<3π/2, the direction of the air from the blade rotation (rΩ) is the same as the direction of horizontal cross flow (Vy). While when the blades are in 3π/2<Ψ<π/2, the direction of the air from the blade rotation (rΩ) is the same as the direction of horizontal cross flow (Vy). The last figure in Figure 1(a) shows that the blade section airfoil has the smallest AOA at the azimuth π. At Ψ=0, the blade section has the largest angle of incidence. Shear wind (only horizontal shear wind is shown here, Figure 1(b)) is generated along the rotor's horizontal inflow velocity with spatial changes. The main reasons for this situation are the changes in the atmosphere itself, environmental and geographical factors, or a combination of factors. In this case, the maximum incident angle to which the airfoil of the blade cross section occurs at Ψ=π/2, and the minimum incident angle occurs at the lowest value of horizontal flow velocity [12]. The positive vertical wind is shown in Figure 1(c). The minimum incident angle of the blade section airfoil appears at 3π/2, and the maximum incident angle appears at π/2 [1].
In practical conditions, the yaw load on wind turbines is more complicated. Under different yaw loads or their combinations, the initial phase of the pitch oscillation can be changed over a wide range. In the following parts, the effects of different initial phases (Φ) of pitch oscillations when coupled with other vibrations will be discussed.

Simulated Cases
The S809 airfoil designed for horizontal axis wind turbine blades, with a relative thickness of 12%, has been studied for more than two decades, and the experimental data are sufficient. Hence this kind of airfoil was chosen in this study. The studied airfoil motion type is shown in Figure 2 [16], which includes three types of pitching motion, flapwise motion and edgewise motion.

Figure 2 Airfoil motion types
The pitching motion usually makes the airfoil sine oscillate around the axis 1/4 chord from the leading edge. As previously discussed, idealized yaw loads can be classified as periodic loads. Meanwhile the AOA of the wind turbine blade airfoil changes periodically under the yaw load, so it can be modeled as an pitch oscillating airfoil.
Pitching airfoil oscillation law: Where α(t) is the instantaneous AOA. αmean is the mean AOA. αamp is the AOA amplitude. f is the oscillation frequency. t is the time. Φ is the phase difference. The value of the Φ are ±3π/4, ±π/2, ±π/4, 0. Since the main research is DS, the mean AOA is 14°, the AOA amplitudes are 5.5° and 10°, which are consistent with the parameters set by Ramsay [10] in the Ohio State University (OSU) wind tunnel test.
The displacement of the blade in the flapwise direction is [22][23]: The displacement of the blade in the edgewise direction is [24]: Where A and B are the airfoil flapwise amplitude and edgewise amplitude respectively, f is the oscillating frequency.
For flapwise motion, the influence of airfoil on its aerodynamic performance is relatively small [25][26][27]. Therefore, for the sake of calculation, only the compound movements under the flapwise motion was investigated in this study. For the selection of flapwise motion parameters, Tuncer et al. [28] used the concept of relative flapwise amplitude (h) in the study, calculated as follows: Where A is the airfoil flapwise amplitude. c is the airfoil chord length.
They studied the aerodynamics of NACA0012 airfoil at h=0.4~0.7, Reynolds number Re=1×10 6 , and reduced frequency k=0.8. Reference [16] studied the DS characteristics of S809 airfoil at relative flapwise and edgewise amplitudes of 0.5 and 0.227. Referring to [16][17][28][29], combined with the actual running status of the fan, it is also convenient for comparative analysis to determine the relative flapwise amplitude h=0.2~0.5. The movement laws and main parameters of each example are shown in Table 1.

Numerical setup
The S809 wind turbine dedicated airfoil was considered. The Mach number is low in the calculation conditions, it can be regarded as incompressible flow. Two-dimensional continuity equations and two-dimensional incompressible N-S equations are used in this study. The SST k-ω model that based on Menter's [30] turbulence model with transition correction was considered. This study uses the moving mesh method to achieve the airfoil movement. The mesh and boundary conditions are shown in Figure 3. The entire calculation domain has a radius of 20 times the airfoil chord length and is divided into two areas. The outer static area uses a triangular mesh, while the inner moving area uses a triangular and quadrilateral mixed mesh. A total of 560 nodes are set on the airfoil surface, and a C-shaped quadrilateral boundary layer mesh is used around it. In order to ensure accurate simulation of the boundary layer flow [31][32], the first layer mesh layout ensures that y+ is on the order of 1.0, and the forward growth factor is set to 1.08. The calculation time step after time independence verification is set to 0.0005s, and the residual is set to 1×10 -6 . A user-defined function (UDF) is used to control the airfoil movement in Fluent software. At least 5 cycles of calculation are performed for each operating condition.

Validation
In order to prove the validity of the numerical method used, the results of CFD calculations under two operating conditions were compared with the unsteady aerodynamic performance of the S809 airfoil in the case of pitch oscillation in the OSU wind tunnel test [10]. The comparison of the results is shown in Figure 4. Since the average AOA is closer to the static stall angle at 14°and the oscillation amplitude is 5.5° (Figure 4(a)), the calculated value of the CL is higher than the experimental value [33], and the agreement of the drag coefficient is relatively good.
Ekaterinaris et al. [34] also pointed out that when the airfoil AOA is close to the static stall angle and the oscillation amplitude is small, the accuracy of the calculation results is reduced. Although the CL is slightly higher than the experimental value at the dynamic stall angle (αDS), it can be seen that the numerical calculation results are in good agreement with the experimental value as a whole, and the details of the DS are also well predicted, which can indicate that the numerical method can effectively calculate this complex flow.  significantly different from the simple pitch motion. The different phases were divided into three parts for discussion.

Φ=π
Take αamp=10°, Φ=π as an example to discuss in detail the aerodynamic situation of the airfoil in compound motion. Figure 6 (o) and (p) shows the hysteresis curve of the CL and CD of the compound motion in the phase. The pressure coefficient (Cp) diagrams of pitch motion and compound motion are given in figure 7. In this figure, the streamlines are superimposed to present more details. The Cp is defined as： where p is the static pressure on the airfoil surface. p  is the static pressure from the far forward flow. and ρ is the air density. U  is the flow velocity from the far forward. 10 15 Compared with only pitching, the trend of the aerodynamic hysteresis curve of the compound motion under different flapwise amplitudes is consistent. When the airfoil starts to move up, the lift increases linearly until the AOA reaches a certain angle and then becomes flat. During this period, due to the superposition of the flapwise motion, the CL and CD of the compound motion are higher than those of the only pitching motion, and the higher amplitude corresponds to a high aerodynamic coefficient. In the uplift phase of the airfoil, the adhering flow on the airfoil surface of the compound motion and the pitch motion shows a large difference, as shown in Figure 7. The pitching motion continues to a high AOA. However, due to the downward movement of the flapwise during compound movement, the flow pressure near the pressure side of the airfoil is greater than that during the simultaneous pitching movement. The pressure side will create a pressure gradient with the suction side. At this time, due to the effect of the pressure gradient, the airflow on the pressure side of the trailing edge bypasses the trailing edge to the suction side earlier. Subsequently, the downstream fluid is forced to return, causing the fluid in the boundary layer to move away from the wall and face the mainstream. Therefore, the boundary layer is separated earlier (the flow separation of the compound motion occurs near 13°, which is about 3.5° earlier than the pitch movement, Figure 7  With reference to Figure 9 (here, the vorticity is dimensionless) and Figure 7 (c) and (d), When h=0.5 and the AOA reaches around 17°. A leading edge vortex (LEV) begins to form, which is about 3° earlier than the pitch motion. In addition, the generated LEV have higher intensity, resulting in higher aerodynamic loads before stalling. However, when LEV are generated, the vortices near the trailing edge are weaker than the pitching motion. As the AOA continues to increase, the LEV grows rapidly on the upper surface of the airfoil. Near 21.8°, the LEV covers the entire upper surface and DS occurs under the maximum aerodynamic load, as shown in Figure 7 (e) and (f). The angle at which the DS of the compound motion occurred was approximately 1.5° earlier than the pitch motion, and the CL during DS increased by 0.38 (Table 3). From the comparison of the Cp diagram, the upper surface of the airfoil has a lower pressure when the compound motion is stalled, and the vorticity of the leading edge covering the suction side of the airfoil is stronger. After that, the LEV falls off while

Pitching motion
Compound motion αamp=10°, Φ=π, h=0.5 forming the trailing edge vortex (TEV), causing a sudden drop in lift and drag, Figure 7 (g) and (h). Then, when the downstroke moves, a LEV is generated again, and the lift force also increases, as shown in Figure 7 (i) and (j). Because the second LEV is not as strong as the first LEV, it cannot significantly increase lift, and it takes much less time to grow and fall off than the first LEV.
The above situation has similar effects on different h values, and the specific situation will be discussed in detail in the next section. In general, when a Φ=π, during the uplift phase of the airfoil, the vibration velocity in the flapwise direction increases the relative flow velocity and AOA of the airfoil. The LEV generates faster and breaks away from the upper surface of the airfoil earlier, causing the airfoil to enter the stall earlier. The main vortex detached from the airfoil surface earlier than the pitch motion. In the downward pitching phase, the equivalent AOA is reduced, which makes the separation flow reattach earlier than the simple pitching motion. More visually, the flapwise motion at Φ=π is equivalent to deepening the stall. This will allow the airfoil to spend more time in the stall area, and the aerodynamic fluctuations during DS will be more intense. Figures 5 and 6 (o), (p), (a), (b), (c), (d), (e), (f), (g) and (h) respectively show the CL and CD curves of the airfoil at Φ =π, -3π/4, -π/2, -π/4 and 0 under conditions αamp=5.5° and αamp=10°.

-π≤Φ≤0
We can find that when αamp=5.5°, almost all stalls occurred in this interval, and the αDS under two operating conditions is delayed as the phase difference increases, while the CL and CD first increase during DS, then decrease. From Figures 8 and 9, we can find that due to the increase in phase difference, the upstroke of the pitching motion gradually overlaps with the upstroke of the flapwise motion. The upstroke of the flapwise motion can reduce the relative flow velocity and AOA of the airfoil, thereby delaying the airfoil. LEV are generated during the ascending and ascending stages and slow the development of LEV. However, the second half of the pitch upstroke is still in the downstroke flapwise state, and the flapwise velocity of the airfoil in the upstroke near the αDS gradually increases until -π/2 reaches the maximum, the CL and CD s reach the highest, and then the velocity gradually decreases. The CL and CD during DS also gradually decrease, which is also the main reason for αamp=5.5° stall. Figure 8 Airfoil combined motion vorticity contours at different flapwise amplitude at αamp=5.5°, Φ=-π/2, 0, π/2 and π. Figure 9 Airfoil combined motion vorticity contours at different flapwise amplitude at αamp=10°, Φ=-π/2, 0, π/2 and π.

0<Φ<π
In Figures 5 and 6 (g), (h), (i), (j), (k), (l), (m), (n), (o) and (p),with the increase of the value Φ, the flapwise downstroke and the pitch upstroke begin to overlap gradually, and the CL curve of the combined motion is gradually higher than that of the pitch only. Meanwhile the overlapping area of the upstroke of the flapwise motion and the upstroke of the pitching movement is reduced and is still in the period of the upstroke of the flapwise motion when the downstroke of the pitching movement starts. It will delay the DS of αamp=10° from 0 to π/2 compared to pitching only. Until Φ=π/2, the velocity of the upstroke of the flapwise motion near the αDS of the pitching motion reaches the maximum. So that the DS of the combined movement of αamp=10° at Φ=π/2 is severely suppressed, and no DS of the combined movement occurs at h=0.5. With αamp=5.5°, the hysteresis curve becomes narrower than that when only pitching is performed. As the amplitude h increases, the aerodynamic load of the airfoil becomes more stable throughout the cycle. As can be seen from Figures 8 and 9, the vorticity changes of the airfoil during the whole movement. After π/2, the flapwise motion tends to be smooth near the αDS of the pitching motion.
When αamp=10°, the αDS of the combined motion gradually approaches only the pitching motion and gradually advances as the value of Φ continues to increase.

Effect of h on aerodynamic performance of airfoil
In order to compare the influence of various factors on the aerodynamic performance of airfoil composite motion, the maximum lift and drag coefficient (CLmax and CDmax) curves at different flapwise amplitudes ( Figure 10) and the αDS curves at different amplitudes ( Figure 11) and Aerodynamic coefficients and dynamic stall angles of airfoils in different motion modes (Table 2 and 3) are given. It can be seen from Figure 10 that the aerodynamic load of the airfoil compound motion under the two conditions shows the same trend. When Φ=π, -3π/4, -π/2, -π/4 and 0, the CL and CD in all cases show an upward trend with increasing amplitude. This is mainly because the stall area of the airfoil is under the downstroke of the flapwise motion. Increasing the amplitude will increase the relative incoming flow velocity of the airfoil, resulting in changes in the aerodynamic load; At Φ=π/4, π/2, and 3π/4, the airfoil stall area is on the flapwise motion, increasing the flapwise amplitude shows the opposite trend.  In terms of stall angle (Figure 11), for the stall situation that has occurred, the two conditions behave identically. Increasing the flapwise amplitude at αamp=10° and Φ=0, π/4, π/2 can delay the stall angle, and the maximum delay of DS observed at Φ =π/4 is 1.408°. Increasing the flapwise amplitude at αamp=10°, Φ =3π/4 and two conditions will advance the stall angle. Under both conditions, Φ=0, π/4, π/2 reach the maximum advance at Φ =π, which are 1.43° and 1.93°, respectively. It is particularly noted that the increase in flapwise amplitude at Φ=-π/4 has little effect on the αDS.
In the previous discussion, the upstroke of the flapwise motion can reduce the relative inflow velocity and AOA of the airfoil, thereby delaying the generation of the LEV during the pitch upstroke phase and slowing the development of the LEV. On this basis, increasing the amplitude has a more severe inhibition on the generation and development of the LEV, and the flapping velocity of the airfoil stall area plays a decisive role in the stall degree under different amplitudes.

Aerodynamic load sensitivity
As mentioned above, the aerodynamic load of the airfoil during compound motion is significantly dependent on the h and Φ parameters, and the aerodynamic load is one of the key parameters in the design of the wind turbine. Therefore, the analysis of the degree of influence of the load is performed here in order to understand h and Φ the degree of influence on the aerodynamic load more intuitively. As can be seen from Figure 10 and CDmax is more pronounced when excluding αamp=10°, Φ=π/2 and αamp=5.5°, Φ=0 and π. In summary, it can be explained that when the flapwise amplitude is in a low range and in some special cases, such as αamp=10°, Φ=π/2 and αamp=5.5°, Φ=0 and π, both h and Φ can greatly affect the dynamic load. The effect is that the dynamic load is more sensitive to Φ than h when the flapwise amplitude is in a high range.

Conclusion
Under yaw loads, the aerodynamic loads on the blades of a horizontal axis wind turbine show periodic changes. In each cycle, the blade section airfoil undergoes complex motion. The SST k-ω turbulence model with transition correction was used to simulate the aerodynamic performance of the DS of the airfoil. The aerodynamic results obtained were in good agreement with the experimental data. By analyzing the aerodynamic data of the complex movement of the airfoil under yaw loads, the following conclusions were drawn.
(1) When the impeller is under yaw load, the phase of the pitching motion and the flapwise motion can oscillate in a wide range. At -π≤Φ≤0, αamp=5.5° will cause a significant stall phenomenon, and the αDS will be gradually delayed as the phase difference increases, and the aerodynamic load will fluctuate sharply. However, at 0<Φ<π, αamp=5.5° could not be stalled. In the αamp=10° state, the αDS was gradually advanced as the phase difference increased, and the aerodynamic load fluctuation was relatively smooth under the two conditions. (2) The effect of h on the aerodynamic performance of the airfoil is mainly reflected in: in the case of -π≤Φ≤0, as the amplitude increases, the CL and CD both show an increasing trend. In other cases, it shows an opposite trend. In the case of DS, when Φ=0, π/4 and π/2, increasing the flapwise amplitude significantly delays the αDS. At Φ=3π/4, π, -3π/4 and -π/2, increasing the flapwise amplitude will advance the stall angle, while at Φ=-π/4, the flapwise amplitude has little effect on the αDS. (3) In the amplitude range of the flapwise motion studied in this study, the phase difference between the pitching motion and the flapwise motion and different flapwise amplitudes can have a large impact on the aerodynamic performance of the airfoil during DS. However, the impact of the two is very different in different situations. When the flapwise amplitude is in the lower range and in some special cases, such as αamp=10°, Φ=π/2 and αamp=5.5°, Φ=0 and π, both h and Φ can have a great impact on dynamic load. When the flapwise amplitude is in a high range, the dynamic load is more sensitive to Φ than h.