Unsteady Aerodynamic Characteristics Simulations of Rotor Airfoil under Oscillating Freestream Velocity

: The dynamic stall characteristics of rotor airfoil are researched by employing unsteady Reynolds-Averaged Navier-Stokes (RANS) method under oscillating freestream velocity conditions. In order to simulate the oscillating freestream velocity of airfoil under dynamic stall conditions, the moving-embedded grid method is employed to simulate the oscillating velocity. By comparing the simulated dynamic stall characteristics of two-dimensional airfoil and three-dimensional rotor, it is indicated that the dynamic stall characteristics of airfoil under oscillating freestream velocity reﬂect the actual dynamic stall characteristics of rotor airfoil in forward ﬂight more accurately. By comparing the simulated results of OA209 airfoil under coupled freestream velocity / pitching oscillation conditions, it is indicated that the dynamic stall characteristics of airfoil associate with the critical value of Cp peaks (i.e., the dynamic stall characteristics of OA209 airfoil would be enhanced when the maximum negative pressure is larger than − 1.08, and suppressed when this value is smaller than − 1.08). By comparing the characteristics of vortices under di ﬀ erent oscillating velocities, it indicates that the dissipation rate of leading edge vortex presents as exponent characteristics, and it is not sensitive to di ﬀ erent oscillating velocities.


Introduction
In forward flight conditions, the freestream velocity of a helicopter rotor blade in the advancing side is larger than that of the retreating side due to the effect of forward flight. As a result, the inflow velocity of rotor airfoil varies with the rotor azimuth. Therefore, investigations on unsteady aerodynamic characteristics of rotor airfoil under oscillating freestream velocity are important, especially for the dynamic stall phenomenon [1,2]. To discover the physical essence of a rotor airfoil dynamic stall, a lot of research under steady freestream velocity (SFV) conditions has been accomplished by using the experimental method [3][4][5][6], theoretical method [7][8][9][10][11][12], and numerical method [13][14][15]. Hence, the dynamic stall of airfoil under SFV conditions is well-studied at present. However, the dynamic stall of rotor airfoil under oscillating freestream velocity (OFV) conditions was seldom taken into account in early researches.
In order to discover the physical phenomenon of dynamic stall under OFV condition, Favier [16] preliminarily researched dynamic stall characteristics of the NACA0012 under the condition of oscillating velocity by using the experimental method. After that, the aerodynamic loads of the NACA0012 airfoil were measured by Gursul [17,18] at a fixed angle of attack (AoA) under unsteady freestream conditions in a vertical unsteady water tunnel with a cross-sectional area of 45.7 × 45.7 cm. Based on this experiment, it was found that the lift coefficient (C l ) can be one order of magnitude higher under OFV conditions than that under SFV conditions. Gharali [19] investigated the dynamic stall of NACA0012 airfoil under variational velocity by employing the ANSYS Fluent software, and the results were indicated that the leading edge vortex (LEV) formation could be obviously affected by the phase angle between the oscillating pitch and the oscillating freestream velocity. However, the maximum AoA and maximum velocity of airfoil were superimposed in this research, which is much different from the actual working environment of a helicopter rotor in forward flight. As a result, the conclusions of the dynamic stall characteristics of two-dimensional (2D) airfoil may not coincide with that of a helicopter rotor airfoil. At the same time, a new experiment [20,21] was performed by the researchers of Ohio State University (OSU) to test the dynamic stall characteristics of SSC-A09 airfoil under oscillating velocity conditions. By comparing the measured C l and C m (pitching moment coefficient) under OFV conditions and SFV conditions, it was found that the lift curve slope and stall angle increased under the OFV conditions. However, due to the restriction of the experimental equipment, the maximum oscillating velocity (0.08 Ma) was much smaller than the normal forward flight velocity of the helicopter (Mach number of 0.1-0.3). Moreover, the aspect ratio of the blade model used in these experiments was only 1.0, which would generate obvious three-dimensional (3D) effects, and the characteristics of dynamic stall vortices would be influenced seriously. What is more, the added-mass effect due to the velocity variation would also influence the measured results in this research as mentioned in another study [22].
At the same time, Jones et al. [23][24][25] researched the dynamic stall characteristics of an airfoil under OFV conditions in horizontal free-surface water tunnel of the U.S. Air Force Research Laboratory. However, these studies were focused on aerodynamic characteristics of airfoil at fixed AoAs, which was much different from the working environment of a helicopter rotor. Meanwhile, the freestream velocities were also smaller than the rotational velocity of the helicopter rotor due to the restriction of water tunnel. Therefore, the dynamic stall characteristics of the rotor airfoil under the real environment of helicopter are still not included in these studies, and this issue is valuable to research more deeply.
As mentioned above, the dynamic stall characteristics of rotor airfoil under an actual working environment of a helicopter are not well studied. Previous research seldom focused on the characteristics of a dynamic stall vortex under coupled freestream velocity/pitching oscillation conditions directly. Therefore, the purpose of this study is to explore the effects of OFV on dynamic stall characteristics of rotor airfoil. In this research, the moving-embedded grid method [26,27] is adopted to simulate the oscillating velocity and pitching of an airfoil. The unsteady RANS equations coupling with third-order Roe-MUSCL spatial discretization scheme are chosen as the governing equations to predict the unsteady flowfield of an airfoil, and the highly-efficient implicit scheme of lower-upper symmetric Gauss-Seidel (LU-SGS) is adopted for temporal discretization. To capture the separated airflow and vortex more accurately, the SST k − ω turbulence model is employed in the Computational Fluid Dynamics (CFD) code. Based on these methods, the dynamic stall characteristics of airfoil are researched under coupled conditions. It is indicated from the simulated results that the dynamic stall characteristics of airfoil associate with a critical value of Cp peak. By studying the vortices characteristics, it is illustrated that the dissipation rate of the LEV presents as exponent characteristics. However, this dissipation rate is not sensitive to oscillating velocity.

Aerodynamic Environment of Helicopter Rotor
The Figure 1a illustrates the velocity distribution of a helicopter rotor with blade tip, Mach number of 0.6, and advance ratio of 0.3. It is indicated that the rotor blade velocity varies with the rotor azimuth in forward flight. As a result, the relative freestream velocity of airfoil at different spanwise sections of rotor blade could be expressed as where, M B represents the rotational velocity of the blade, M F denotes the forward flight velocity of the helicopter, ω and t represent the angular velocity and time, respectively.
Aiming at illustrating the relationship between the basic freestream velocity and the oscillating velocity, a normalized oscillating velocity λ is defined as The AoA variation induced by cyclic pitching of the rotor blade can be written as where, α m and ∆α represent the mean AoA and amplitude of pitch oscillation, respectively. In this environment, the variations of AoA (α m = 8.0 Aiming at illustrating the relationship between the basic freestream velocity and the oscillating velocity, a normalized oscillating velocity λ is defined as

Numerical Simulation Method
In order to generate the computational grids around the OA209 airfoil, the Poisson equations [28] are chosen as the governing equations under 2D conditions. By solving these equations, the Ctopology grids with 459 × 80 points are generated. To satisfy the requirement of OFV simulation, the far-field boundary of airfoil grids is fixed at 25 times the airfoil chord, and the y + of the grids near the wall surface is smaller than 1.2. In this research, the background grid is composed of 300 × 300 points with a far-field boundary of 50 c (airfoil chord), and the inverse map method [29] is employed in this code to search donor element. As shown in Figure 2, the oscillating velocity is simulated by moving the airfoil back and forth periodically [30,31].

Numerical Simulation Method
In order to generate the computational grids around the OA209 airfoil, the Poisson equations [28] are chosen as the governing equations under 2D conditions. By solving these equations, the C-topology grids with 459 × 80 points are generated. To satisfy the requirement of OFV simulation, the far-field boundary of airfoil grids is fixed at 25 times the airfoil chord, and the y + of the grids near the wall surface is smaller than 1.2. In this research, the background grid is composed of 300 × 300 points with a far-field boundary of 50 c (airfoil chord), and the inverse map method [29] is employed in this code to search donor element. As shown in Figure 2, the oscillating velocity is simulated by moving the airfoil back and forth periodically [30,31]. In order to simulate the unsteady compressible flowfield around an airfoil, the integral form of the Navier-Stokes equations is employed in this work, i.e., where W represents the conserved variables, c F and v F denote convective fluxes and viscous fluxes, respectively, i.e., V is contravariant velocity. ij τ denotes the viscous stresses, and i Θ is the term describing the heat condition in the fluid.
Meanwhile, the unsteady condition simulation is accomplished by employing the dual timestepping approach. Then Equation (4) would be changed as where * W is the approximation for 1 n + W , τ represents a pseudo time variable, and the unsteady residual is defined as To improve the computational efficiency, the implicit LU-SGS [32,33] scheme is employed in this work. The two equation SST k ω − turbulence model is employed to simulate the viscous stresses of the flowfield.
In order to verify the accuracy of the present grid system and flowfield solver, a case of the NACA4412 airfoil is presented in this section. In this case, the pressure coefficient (Cp) simulated by the present CFD method compared with the test data [34] is shown in Figure 3a. In this case, three different grids (699 × 100, 459 × 80, 319 × 60) are used to analyze the grid sensitivity. By comparing with the test data, it illustrated that the Cp simulated by the CFD method correlates well with test data for a different grid. However, the convergence of a larger grid is lower than coarse grid, but it still converges at least five orders of magnitude. It is indicated that the grid used in this research could satisfy the calculation requirement. Meanwhile, it is illustrated in Figure 3b that the velocity profiles are also close to the test data. As a result, it is indicated that the grid used in this research is suitable. In order to simulate the unsteady compressible flowfield around an airfoil, the integral form of the Navier-Stokes equations is employed in this work, i.e.,

∂ ∂t
where W represents the conserved variables, F c and F v denote convective fluxes and viscous fluxes, respectively, i.e., where V r = V − V t , V is absolute velocity, V t is contravariant velocity. τ ij denotes the viscous stresses, and Θ i is the term describing the heat condition in the fluid.
Meanwhile, the unsteady condition simulation is accomplished by employing the dual time-stepping approach. Then Equation (4) would be changed as where W * is the approximation for W n+1 , τ represents a pseudo time variable, and the unsteady residual is defined as To improve the computational efficiency, the implicit LU-SGS [32,33] scheme is employed in this work. The two equation SST k − ω turbulence model is employed to simulate the viscous stresses of the flowfield.
In order to verify the accuracy of the present grid system and flowfield solver, a case of the NACA4412 airfoil is presented in this section. In this case, the pressure coefficient (Cp) simulated by the present CFD method compared with the test data [34] is shown in Figure 3a. In this case, three different grids (699 × 100, 459 × 80, 319 × 60) are used to analyze the grid sensitivity. By comparing with the test data, it illustrated that the Cp simulated by the CFD method correlates well with test data for a different grid. However, the convergence of a larger grid is lower than coarse grid, but it still converges at least five orders of magnitude. It is indicated that the grid used in this research could satisfy the calculation requirement. Meanwhile, it is illustrated in Figure 3b that the velocity profiles are also close to the test data. As a result, it is indicated that the grid used in this research is suitable.  , and the reduced frequency (k) is 0.151. As shown in Figure 4, it is indicated that the simulated Cl by the present CFD method is correlated well with the test data of National Aeronautics and Space Administration (NASA) [35]. Meanwhile, the calculated results of the present CFD method are better than that calculated by Fluent software with SST k ω − turbulence model.
Since we lack appropriate test data of rotor airfoil under dynamic stall conditions coupled with variational velocity, the theory of Isaacs [36] is employed to verify the accuracy of the present CFD method under OFV condition, as shown in Figure 5. It is illustrated that the numerical results correlate well with the theoretical data at could be attributed to the fact that the airflow compressibility and viscidity are neglected in the Isaacs theory, where Cl0 represents the lift force at steady freestream velocity. Consequently, it is indicated that the present CFD method is suitable for simulating unsteady characteristics of a rotor airfoil under dynamic stall conditions coupled with oscillating freestream velocity.   In order to verify the accuracy of the present CFD method used in the dynamic stall characteristic simulations of rotor airfoil, a typical case of dynamic stall about the NACA0012 airfoil is presented under the SFV condition in this section. The freestream Mach number is 0.283, the variation of AoA is α = 14.91 • + 9.88 • sin(ωt), and the reduced frequency (k) is 0.151. As shown in Figure 4, it is indicated that the simulated C l by the present CFD method is correlated well with the test data of National Aeronautics and Space Administration (NASA) [35]. Meanwhile, the calculated results of the present CFD method are better than that calculated by Fluent software with SST k − ω turbulence model. Since we lack appropriate test data of rotor airfoil under dynamic stall conditions coupled with variational velocity, the theory of Isaacs [36] is employed to verify the accuracy of the present CFD method under OFV condition, as shown in Figure 5. It is illustrated that the numerical results correlate well with the theoretical data at λ = 0.4, and the numerical results with λ = 0.8 are basically close to the theoretical value. The deviation of the C l /C l0 at λ = 0.8 could be attributed to the fact that the airflow compressibility and viscidity are neglected in the Isaacs theory, where C l0 represents the lift force at steady freestream velocity. Consequently, it is indicated that the present CFD method is suitable for simulating unsteady characteristics of a rotor airfoil under dynamic stall conditions coupled with oscillating freestream velocity.  , and the reduced frequency (k) is 0.151. As shown in Figure 4, it is indicated that the simulated Cl by the present CFD method is correlated well with the test data of National Aeronautics and Space Administration (NASA) [35]. Meanwhile, the calculated results of the present CFD method are better than that calculated by Fluent software with SST k ω − turbulence model.
Since we lack appropriate test data of rotor airfoil under dynamic stall conditions coupled with variational velocity, the theory of Isaacs [36] is employed to verify the accuracy of the present CFD method under OFV condition, as shown in Figure 5. It is illustrated that the numerical results correlate well with the theoretical data at could be attributed to the fact that the airflow compressibility and viscidity are neglected in the Isaacs theory, where Cl0 represents the lift force at steady freestream velocity. Consequently, it is indicated that the present CFD method is suitable for simulating unsteady characteristics of a rotor airfoil under dynamic stall conditions coupled with oscillating freestream velocity.

The Dynamic Stall Characteristics of Airfoil Under 3D Conditions
Aiming at researching the dynamic stall characteristics of rotor airfoil with OFV under actual rotor working conditions, a two-blade untwisted rotor model with a rectangular planform is designed based on the OA209 airfoil, and the aspect ratio of the rotor blade is 15 Table 1. In order to verify the accuracy of the present CFD method used in the 3D rotor condition, a case of SA349/2 rotor with an advance ratio of 0.26 is present in this section, just as shown in Figures 6 and 7. The comparisons of normal force (Cn) and Cp distributions (0.75 r/R) are correlated well with the test data [37]. It is illustrated that the present CFD method can effectively simulate the unsteady aerodynamic characteristics of a rotor.

The Dynamic Stall Characteristics of Airfoil Under 3D Conditions
Aiming at researching the dynamic stall characteristics of rotor airfoil with OFV under actual rotor working conditions, a two-blade untwisted rotor model with a rectangular planform is designed based on the OA209 airfoil, and the aspect ratio of the rotor blade is 15 Table 1. In order to verify the accuracy of the present CFD method used in the 3D rotor condition, a case of SA349/2 rotor with an advance ratio of 0.26 is present in this section, just as shown in Figures 6 and 7. The comparisons of normal force (Cn) and Cp distributions (0.75 r/R) are correlated well with the test data [37]. It is illustrated that the present CFD method can effectively simulate the unsteady aerodynamic characteristics of a rotor.
The comparisons of dynamic stall characteristics of the OA209 airfoil between the 2D and 3D conditions are shown in Figure 8. It is illustrated that the simulated Cl under the OFV condition is much closer to that simulated under the 3D rotor condition compared with the simulated results under SFV conditions at sections of 0.6 R and 0.7 R. Meanwhile, it can be noticed that the simulated results under the 2D OFV conditions are smaller than that of the 3D simulation between the azimuth The comparisons of dynamic stall characteristics of the OA209 airfoil between the 2D and 3D conditions are shown in Figure 8. It is illustrated that the simulated Cl under the OFV condition is much closer to that simulated under the 3D rotor condition compared with the simulated results under SFV conditions at sections of 0.6 R and 0.7 R. Meanwhile, it can be noticed that the simulated results under the 2D OFV conditions are smaller than that of the 3D simulation between the azimuth The comparisons of dynamic stall characteristics of the OA209 airfoil between the 2D and 3D conditions are shown in Figure 8. It is illustrated that the simulated C l under the OFV condition is much closer to that simulated under the 3D rotor condition compared with the simulated results under SFV conditions at sections of 0.6 R and 0.7 R. Meanwhile, it can be noticed that the simulated results under the 2D OFV conditions are smaller than that of the 3D simulation between the azimuth of 240.0 • to 360.0 • , because the airflow separation is restricted by the spanwise flow under the 3D condition [38]. Therefore, the lift stall would be postponed. From this case, it is illustrated that the dynamic stall characteristics of an airfoil under the OFV condition could reflect the actual unsteady aerodynamic characteristics of rotor airfoil in forward flight more accurately. Therefore, the dynamic stall characteristics of rotor airfoil need to be researched under the OFV conditions more deeply. of 240.0° to 360.0°, because the airflow separation is restricted by the spanwise flow under the 3D condition [38]. Therefore, the lift stall would be postponed. From this case, it is illustrated that the dynamic stall characteristics of an airfoil under the OFV condition could reflect the actual unsteady aerodynamic characteristics of rotor airfoil in forward flight more accurately. Therefore, the dynamic stall characteristics of rotor airfoil need to be researched under the OFV conditions more deeply.

Aerodynamic Loads of Airfoil under Coupled Condition
Generally, the adverse pressure gradient near the leading edge of airfoil would be enlarged by increasing the AoA at fixed freestream velocity. On the contrary, it would be reduced by decreasing the freestream velocity at fixed AoA. Therefore, the dynamic stall characteristics of airfoil would be affected by both AoA and freestream velocity. In order to explore these characters more deeply, two cases of dynamic stall coupled with OFV are presented in this research.

Aerodynamic Loads of Airfoil under Coupled Condition
Generally, the adverse pressure gradient near the leading edge of airfoil would be enlarged by increasing the AoA at fixed freestream velocity. On the contrary, it would be reduced by decreasing the freestream velocity at fixed AoA. Therefore, the dynamic stall characteristics of airfoil would be affected by both AoA and freestream velocity. In order to explore these characters more deeply, two cases of dynamic stall coupled with OFV are presented in this research.  of 240.0° to 360.0°, because the airflow separation is restricted by the spanwise flow under the 3D condition [38]. Therefore, the lift stall would be postponed. From this case, it is illustrated that the dynamic stall characteristics of an airfoil under the OFV condition could reflect the actual unsteady aerodynamic characteristics of rotor airfoil in forward flight more accurately. Therefore, the dynamic stall characteristics of rotor airfoil need to be researched under the OFV conditions more deeply.

Aerodynamic Loads of Airfoil under Coupled Condition
Generally, the adverse pressure gradient near the leading edge of airfoil would be enlarged by increasing the AoA at fixed freestream velocity. On the contrary, it would be reduced by decreasing the freestream velocity at fixed AoA. Therefore, the dynamic stall characteristics of airfoil would be affected by both AoA and freestream velocity. In order to explore these characters more deeply, two cases of dynamic stall coupled with OFV are presented in this research.  The streamlines and pressure contour around the OA209 airfoil under the coupled conditions are shown in Figure 10. It can be seen that the LEV at the AoA of 12.84 • during the upstroke process is not obvious under the SFV condition. However, the LEV is observed near the leading edge of airfoil for the state of λ = 0.5. It is indicated that the LEV would form earlier under the condition of OFV. Due to the induction of LEV and the influence of trailing edge of airfoil, the trailing edge vortex (TEV) forms near the trailing edge of airfoil. The AoAs of LEV and TEV formation are shown in Table 2. It is illustrated that the AoAs of LEV and TEV formation are both decreased with the increase of normalized oscillating velocity. The streamlines and pressure contour around the OA209 airfoil under the coupled conditions are shown in Figure 10. It can be seen that the LEV at the AoA of 12.84° during the upstroke process is not obvious under the SFV condition. However, the LEV is observed near the leading edge of airfoil for the state of 0.5 λ = . It is indicated that the LEV would form earlier under the condition of OFV. Due to the induction of LEV and the influence of trailing edge of airfoil, the trailing edge vortex (TEV) forms near the trailing edge of airfoil. The AoAs of LEV and TEV formation are shown in Table 2. It is illustrated that the AoAs of LEV and TEV formation are both decreased with the increase of normalized oscillating velocity.   Figure 11 shows the pressure coefficient (Cp2, normalized by 2 / 2 a ρ ∞ ∞ ) of the OA209 airfoil on the upper surface through a pitching cycle. It is illustrated that the Cp2 near the leading edge of airfoil under the OFV condition is smaller than that under the SFV condition at the same AoA between the azimuth of 60° to 180°, because the increase of freestream velocity under the OFV condition would enlarge the pressure. What is more, another feature can be noticed in that the negative pressure peak is restricted at Cp2 = −1.08 due to the effect of the LEV.   Figure 11 shows the pressure coefficient (Cp 2 , normalized by ρ ∞ a 2 ∞ /2) of the OA209 airfoil on the upper surface through a pitching cycle. It is illustrated that the Cp 2 near the leading edge of airfoil under the OFV condition is smaller than that under the SFV condition at the same AoA between the azimuth of 60 • to 180 • , because the increase of freestream velocity under the OFV condition would enlarge the pressure. What is more, another feature can be noticed in that the negative pressure peak is restricted at Cp 2 = −1.08 due to the effect of the LEV.  Table 3, it is illustrated that the peak of Cl with 0.5 λ = is 14.5% higher than that of 0.0 λ = . Meanwhile, the divergence angles of the Cm and Cd are decreased when the oscillating velocity increases. Additionally, the peaks of Cm and Cd loops under the OFV condition are also larger than that of the SFV condition, and the peak values are shown in Table 3. Therefore, it is illustrated that the unsteady aerodynamic characteristics of airfoil under the OFV condition can be aggravated due to the effects of dynamic stall vortex.   The aerodynamic loads of the OA209 airfoil are shown in Figure 12. It can be noticed that the C l is increased by enlarging the normalized oscillating velocity during the upstroke process, and the lift curve slopes at a large AoA, for the different oscillating velocities are 0.371 (line a for λ = 0.5), 0.273 (line b for λ = 0.25) and 0.228 (line c for λ = 0.0), respectively. As a result, the peaks of C l are also enlarged under the OFV condition. This characteristic indicates that lift induced by LEV would be enlarged as the oscillating velocity increases. From Table 3, it is illustrated that the peak of C l with λ = 0.5 is 14.5% higher than that of λ = 0.0. Meanwhile, the divergence angles of the C m and C d are decreased when the oscillating velocity increases. Additionally, the peaks of C m and C d loops under the OFV condition are also larger than that of the SFV condition, and the peak values are shown in Table 3. Therefore, it is illustrated that the unsteady aerodynamic characteristics of airfoil under the OFV condition can be aggravated due to the effects of dynamic stall vortex.  Table 3, it is illustrated that the peak of Cl with 0.5 λ = is 14.5% higher than that of 0.0 λ = . Meanwhile, the divergence angles of the Cm and Cd are decreased when the oscillating velocity increases. Additionally, the peaks of Cm and Cd loops under the OFV condition are also larger than that of the SFV condition, and the peak values are shown in Table 3. Therefore, it is illustrated that the unsteady aerodynamic characteristics of airfoil under the OFV condition can be aggravated due to the effects of dynamic stall vortex.    The dynamic stall under the coupled condition with k = 0.05 and α = 10.0 presented in this case. The comparisons of aerodynamic loads are shown in Figure 13. It is illustrated that the dynamic stall characteristics of the OA209 airfoil is obviously alleviated under the OFV condition (i.e., the hysteresis loop area of the dynamic stall is decreased). Meanwhile, this tendency is enhanced by increasing the oscillating velocity (i.e., the stall AoA of C l is postponed from 14.   By comparing the simulated results in these two cases, it is indicated that the dynamic stall characteristics of airfoil are associated with the pressure near the leading edge of airfoil, which would induce airflow separation and form a LEV. If the maximum pressure exceeds the critical value of leading edge separation under dynamic stall conditions, the unsteady characteristics would be enhanced with increased oscillating velocity. Otherwise, the dynamic stall characteristics would be inhibited, because the increment of adverse pressure gradient due to the increase of AoA would be restricted by decreasing freestream velocity.

Circulation of Flowfield
Aimed at investigating the characteristics of a separated vortex, the circulation of flowfield is estimated from the local vorticity field via numerical integration based on Stokes' theorem [39]:   By comparing the simulated results in these two cases, it is indicated that the dynamic stall characteristics of airfoil are associated with the pressure near the leading edge of airfoil, which would induce airflow separation and form a LEV. If the maximum pressure exceeds the critical value of leading edge separation under dynamic stall conditions, the unsteady characteristics would be enhanced with increased oscillating velocity. Otherwise, the dynamic stall characteristics would be inhibited, because the increment of adverse pressure gradient due to the increase of AoA would be restricted by decreasing freestream velocity.

Circulation of Flowfield
Aimed at investigating the characteristics of a separated vortex, the circulation of flowfield is estimated from the local vorticity field via numerical integration based on Stokes' theorem [39]: where A represents the integral window encompassing the vortex, ω z represents the vorticity. In order to research the vorticity variation of LEV and TEV, the circulations are integrated in the simulated flowfield (ω z < 0.0 for the LEV, and ω z > 0.0 for the TEV). Meanwhile, the vorticity in the boundary layer is neglected in the integration. Since the second vortex induced by the LEV usually has the same rotation as the TEV, the TEV integral window is started behind 0.6c position of airfoil to eliminate the interference of the second vortex, just as shown in Figure 15. By comparing the simulated results in these two cases, it is indicated that the dynamic stall characteristics of airfoil are associated with the pressure near the leading edge of airfoil, which would induce airflow separation and form a LEV. If the maximum pressure exceeds the critical value of leading edge separation under dynamic stall conditions, the unsteady characteristics would be enhanced with increased oscillating velocity. Otherwise, the dynamic stall characteristics would be inhibited, because the increment of adverse pressure gradient due to the increase of AoA would be restricted by decreasing freestream velocity.

Circulation of Flowfield
Aimed at investigating the characteristics of a separated vortex, the circulation of flowfield is estimated from the local vorticity field via numerical integration based on Stokes' theorem [39]: for the TEV). Meanwhile, the vorticity in the boundary layer is neglected in the integration. Since the second vortex induced by the LEV usually has the same rotation as the TEV, the TEV integral window is started behind 0.6c position of airfoil to eliminate the interference of the second vortex, just as shown in Figure 15.  The variations of circulation at different oscillating velocities on upstroke process are shown in Figure 16, and the simulated state is the same as case 1 (i.e., k = 0.075, α = 12.0 • + 8.0 • sin(ωt + π)).
It can be seen from Figure 16a that the circulation of LEV at λ = 0.0 is less than that under the OFV condition when the AoA is smaller than 15.5 • , because the flowfield under the SFV condition is more stable. Meanwhile, it is indicated that the AoAs of LEV formation are decreased by increasing the oscillating velocity, since the variation of freestream velocity could increase the instability of the boundary layer of airfoil, which induces the airflow separation easily. Meanwhile, it can be seen from Figure 16b that the LEV generates earlier under OFV conditions compared with SFV conditions. However, the accumulating rate of LEV circulation at λ = 0.25 and λ = 0.5 is smaller than the circulation at λ = 0.0, because the freestream velocity is decreased when the AoA increases under the OFV condition, which could reduce the energy supplement from the airflow into the LEV. As a result, the circulation would be decreased with the increase of oscillating freestream velocity. In addition, it also can be seen that the circulation of TEV has the same variational characters due to the circulation conservation of flowfield. Therefore, it is indicated that the circulation variations under the OFV condition are milder compared with the SFV condition.
Appl. Sci. 2018, 8, x FOR PEER REVIEW 13 of 18 oscillating velocity, since the variation of freestream velocity could increase the instability of the boundary layer of airfoil, which induces the airflow separation easily. Meanwhile, it can be seen from Figure 16b that the LEV generates earlier under OFV conditions compared with SFV conditions. However, the accumulating rate of LEV circulation at 0.25 λ = and 0.5 λ = is smaller than the circulation at 0.0 λ = , because the freestream velocity is decreased when the AoA increases under the OFV condition, which could reduce the energy supplement from the airflow into the LEV. As a result, the circulation would be decreased with the increase of oscillating freestream velocity. In addition, it also can be seen that the circulation of TEV has the same variational characters due to the circulation conservation of flowfield. Therefore, it is indicated that the circulation variations under the OFV condition are milder compared with the SFV condition. In order to analyze the dynamic stall vortex characteristics deeply, Figure 17 shows the circulation variations of case 2 (i.e., k = 0.05, In addition, it also can be seen that the circulation of TEV is also stable for this case due to the circulation conservation of the flowfield. As a result, the dynamic stall characteristics of airfoil would be alleviated under the OFV condition with light dynamic stall, and this variational characteristic is similar to the aerodynamic loads.  In order to analyze the dynamic stall vortex characteristics deeply, Figure 17 shows the circulation variations of case 2 (i.e., k = 0.05, α = 10.0 • + 6.0 • sin(ωt + π)), at different oscillating velocities. It is illustrated in this figure that the circulation of LEV at λ = 0.0 obviously increases when the AoA is larger than 12.0 • . However, the circulation of LEV at λ = 0.5 almost keeps stable at different AoAs ( Figure 17a) or times (Figure 17b). It is indicated that the LEV is basically eliminated in the state of λ = 0.5. In addition, it also can be seen that the circulation of TEV is also stable for this case due to the circulation conservation of the flowfield. As a result, the dynamic stall characteristics of airfoil would be alleviated under the OFV condition with light dynamic stall, and this variational characteristic is similar to the aerodynamic loads. oscillating velocity, since the variation of freestream velocity could increase the instability of the boundary layer of airfoil, which induces the airflow separation easily. Meanwhile, it can be seen from Figure 16b that the LEV generates earlier under OFV conditions compared with SFV conditions. However, the accumulating rate of LEV circulation at 0.25 λ = and 0.5 λ = is smaller than the circulation at 0.0 λ = , because the freestream velocity is decreased when the AoA increases under the OFV condition, which could reduce the energy supplement from the airflow into the LEV. As a result, the circulation would be decreased with the increase of oscillating freestream velocity. In addition, it also can be seen that the circulation of TEV has the same variational characters due to the circulation conservation of flowfield. Therefore, it is indicated that the circulation variations under the OFV condition are milder compared with the SFV condition. In order to analyze the dynamic stall vortex characteristics deeply, Figure 17 shows the circulation variations of case 2 (i.e., k = 0.05, 10 In addition, it also can be seen that the circulation of TEV is also stable for this case due to the circulation conservation of the flowfield. As a result, the dynamic stall characteristics of airfoil would be alleviated under the OFV condition with light dynamic stall, and this variational characteristic is similar to the aerodynamic loads.

Dissipation Process of Separated Vortex
Since the separated vortex has significant influence on airfoil unsteady aerodynamic force characteristics, the vortex circulation near the upper surface should be researched deeply. As a result, in order to explore the dissipation process of the separated vortex, the integration with fixed windows is accomplished in this research, and the region of integral windows is 0.1 c × 0.08 c for LEV, and 0.08 c × 0.08 c for TEV at different AoAs, just as shown in Figure 18. The dissipation of vortex circulation can be expressed as where Γ 0 denotes original vortex circulation, κ denotes dissipation coefficient of LEV.
In order to analyze the influence of the integral windows on the dissipation process of the vortex, the circulation of LEV is estimated in different integral windows (i.e., 0.1 c × 0.08 c (large window) and 0.05 c × 0.04 c (small window)), then the variations of the LEV circulations of case 1 with λ = 0.5 are shown in Figure 19. By fitting the samples of LEV circulation, the dissipation coefficients (κ) of LEV for different integral windows are 0.447 for 0.1 c × 0.08 c and 0.451 for 0.05 c × 0.04 c, respectively. It is indicated that the dissipation rate of LEV with different integral windows is similar.

Dissipation Process of Separated Vortex
Since the separated vortex has significant influence on airfoil unsteady aerodynamic force characteristics, the vortex circulation near the upper surface should be researched deeply. As a result, in order to explore the dissipation process of the separated vortex, the integration with fixed windows is accomplished in this research, and the region of integral windows is 0.1 c × 0.08 c for LEV, and 0.08 c × 0.08 c for TEV at different AoAs, just as shown in Figure 18. The dissipation of vortex circulation can be expressed as where 0 Γ denotes original vortex circulation, κ denotes dissipation coefficient of LEV.
In order to analyze the influence of the integral windows on the dissipation process of the vortex, the circulation of LEV is estimated in different integral windows (i.e., 0.1 c × 0.08 c (large window) and 0.05 c × 0.04 c (small window)), then the variations of the LEV circulations of case 1 with 0.5 λ = are shown in Figure 19. By fitting the samples of LEV circulation, the dissipation coefficients ( κ ) of LEV for different integral windows are 0.447 for 0.1 c × 0.08 c and 0.451 for 0.05 c × 0.04 c, respectively. It is indicated that the dissipation rate of LEV with different integral windows is similar.  The dissipation process of LEV circulation in Figure 20a illustrates that the circulation of LEV under the OFV condition is smaller than that under the SFV condition at the same AoA. Therefore, it is indicated that the LEV forms and sheds at smaller AoAs under the OFV condition, and this tendency accelerates with the increase in oscillating velocity. In order to explore the dissipation rate of LEV, the variations of LEV circulation with time are shown in Figure 20b. It is illustrated that the dissipation rates of different LEVs are almost the same. As a result, a new dissipation coefficient for different freestream velocities is calculated by fitting the three group vortexes (i.e., 0.426 ). It is indicated that the oscillating velocity has no obvious influence on the dissipation

Dissipation Process of Separated Vortex
Since the separated vortex has significant influence on airfoil unsteady aerodynamic force characteristics, the vortex circulation near the upper surface should be researched deeply. As a result, in order to explore the dissipation process of the separated vortex, the integration with fixed windows is accomplished in this research, and the region of integral windows is 0.1 c × 0.08 c for LEV, and 0.08 c × 0.08 c for TEV at different AoAs, just as shown in Figure 18. The dissipation of vortex circulation can be expressed as where 0 Γ denotes original vortex circulation, κ denotes dissipation coefficient of LEV.
In order to analyze the influence of the integral windows on the dissipation process of the vortex, the circulation of LEV is estimated in different integral windows (i.e., 0.1 c × 0.08 c (large window) and 0.05 c × 0.04 c (small window)), then the variations of the LEV circulations of case 1 with 0.5 λ = are shown in Figure 19. By fitting the samples of LEV circulation, the dissipation coefficients ( κ ) of LEV for different integral windows are 0.447 for 0.1 c × 0.08 c and 0.451 for 0.05 c × 0.04 c, respectively. It is indicated that the dissipation rate of LEV with different integral windows is similar.  The dissipation process of LEV circulation in Figure 20a illustrates that the circulation of LEV under the OFV condition is smaller than that under the SFV condition at the same AoA. Therefore, it is indicated that the LEV forms and sheds at smaller AoAs under the OFV condition, and this tendency accelerates with the increase in oscillating velocity. In order to explore the dissipation rate of LEV, the variations of LEV circulation with time are shown in Figure 20b. It is illustrated that the dissipation rates of different LEVs are almost the same. As a result, a new dissipation coefficient for different freestream velocities is calculated by fitting the three group vortexes (i.e., 0.426 ). It is indicated that the oscillating velocity has no obvious influence on the dissipation  The dissipation process of LEV circulation in Figure 20a illustrates that the circulation of LEV under the OFV condition is smaller than that under the SFV condition at the same AoA. Therefore, it is indicated that the LEV forms and sheds at smaller AoAs under the OFV condition, and this tendency accelerates with the increase in oscillating velocity. In order to explore the dissipation rate of LEV, the variations of LEV circulation with time are shown in Figure 20b. It is illustrated that the dissipation rates of different LEVs are almost the same. As a result, a new dissipation coefficient for different freestream velocities is calculated by fitting the three group vortexes (i.e., κ = 0.426(Γ 0 = 0.228)). It is indicated that the oscillating velocity has no obvious influence on the dissipation rate of LEV. The variations of TEV circulation are shown in Figure 21. It is illustrated that the dissipation rate of TEV under the condition of OFV is slower than that under the SFV condition. This tendency could be attributed to the reason that the freestream velocity under the OFV condition is smaller than that under the SFV condition when the TEV is forming. Therefore, more energy would be supplied to the TEV from the lower surface of the airfoil since the TEV has a longer attachment time. Consequently, the dissipation rate of TEV decreased under the OFV condition.  Figure 21. It is illustrated that the dissipation rate of TEV under the condition of OFV is slower than that under the SFV condition. This tendency could be attributed to the reason that the freestream velocity under the OFV condition is smaller than that under the SFV condition when the TEV is forming. Therefore, more energy would be supplied to the TEV from the lower surface of the airfoil since the TEV has a longer attachment time. Consequently, the dissipation rate of TEV decreased under the OFV condition.

Conclusions
An unsteady numerical method incorporated with the moving-embedded grid method is developed to simulate dynamic stall characteristics of a rotor airfoil under OFV conditions, and two conclusions can be summarized as: (1). By comparing the simulated result, it is indicated that the dynamic stall characteristics of the OA209 airfoil would be enhanced with increased oscillating velocity when the maximum negative pressure exceeds the critical value (−1.08). On the contrary, it would be inhibited when the maximum negative pressure is less than this critical value. (2). By comparing the dissipation processes of vortex circulation at a fixed integral window, it indicates that the dissipation of LEV presents as exponent characteristics under different  Figure 21. It is illustrated that the dissipation rate of TEV under the condition of OFV is slower than that under the SFV condition. This tendency could be attributed to the reason that the freestream velocity under the OFV condition is smaller than that under the SFV condition when the TEV is forming. Therefore, more energy would be supplied to the TEV from the lower surface of the airfoil since the TEV has a longer attachment time. Consequently, the dissipation rate of TEV decreased under the OFV condition.

Conclusions
An unsteady numerical method incorporated with the moving-embedded grid method is developed to simulate dynamic stall characteristics of a rotor airfoil under OFV conditions, and two conclusions can be summarized as: (1). By comparing the simulated result, it is indicated that the dynamic stall characteristics of the OA209 airfoil would be enhanced with increased oscillating velocity when the maximum negative pressure exceeds the critical value (−1.08). On the contrary, it would be inhibited when the maximum negative pressure is less than this critical value. (2). By comparing the dissipation processes of vortex circulation at a fixed integral window, it indicates that the dissipation of LEV presents as exponent characteristics under different

Conclusions
An unsteady numerical method incorporated with the moving-embedded grid method is developed to simulate dynamic stall characteristics of a rotor airfoil under OFV conditions, and two conclusions can be summarized as: (1). By comparing the simulated result, it is indicated that the dynamic stall characteristics of the OA209 airfoil would be enhanced with increased oscillating velocity when the maximum negative pressure exceeds the critical value (−1.08). On the contrary, it would be inhibited when the maximum negative pressure is less than this critical value. (2). By comparing the dissipation processes of vortex circulation at a fixed integral window, it indicates that the dissipation of LEV presents as exponent characteristics under different oscillating velocities. Meanwhile, the dissipation rate of LEV is not sensitive to different oscillating velocities. However, the dissipation rate of TEV would be reduced by increasing the oscillating velocity.