Effect of Single-Row and Double-Row Passive Vortex Generators on the Deep Dynamic Stall of a Wind Turbine Airfoil †

Passive vortex generators (VGs) have been widely applied on wind turbines to boost the aerodynamic performance. Although VGs can delay the onset of static stall, the effect of VGs on dynamic stall is still incompletely understood. Therefore, this paper aims at investigating the deep dynamic stall of NREL S809 airfoil controlled by single-row and double-row VGs. The URANS method with VGs fully resolved is used to simulate the unsteady airfoil flow. Firstly, both single-row and double-row VGs effectively suppress the flow separation and reduce the fluctuations in aerodynamic forces when the airfoil pitches up. The maximum lift coefficient is therefore increased beyond 40%, and the onset of deep dynamic stall is also delayed. This suggests that deep dynamic-stall behaviors can be properly controlled by VGs. Secondly, there is a great difference in aerodynamic performance between single-row and double-row VGs when the airfoil pitches down. Single-row VGs severely reduce the aerodynamic pitch damping by 64%, thereby undermining the torsional aeroelastic stability of airfoil. Double-row VGs quickly restore the decreased aerodynamic efficiency near the maximum angle of attack, and also significantly accelerate the flow reattachment. The second-row VGs can help the near-wall flow to withstand the adverse pressure gradient and then suppress the trailing-edge flow separation, particularly during the downstroke process. Generally, double-row VGs are better than single-row VGs concerning controlling deep dynamic stall. This work also gives a performance assessment of VGs in controlling the highly unsteady aerodynamic forces of a wind turbine airfoil.


Introduction
Although passive vortex generators (VGs) are very simple, they have been proven to suppress the flow separation effectively and then boost the aerodynamic performance of horizontal axis wind turbines (HAWTs) [1]. Conventional VGs are composed of some pairs of vanes sticking out from the surface, angled to the incoming flow [2]. The height of VGs is close to the boundary-layer thickness. The fundamental principle of VGs is to produce streamwise vortices. These vortices can reenergize the boundary layer to resist the adverse pressure gradient.
The effectiveness of VG designs is primarily determined by the evolution of streamwise vortices. This vortex evolution is further impacted by various VG parameters. Godard and Stanislas [3] measured the boundary layer flow of a two-dimensional bump with VGs, using stereo particle image velocimetry (PIV) and hot-film probes. They found the triangular VGs better than rectangular VGs in decreasing the drag penalty at low angles of attack (AOAs). The counter-rotating configuration was also found better than the co-rotating one in generating upwash and downwash wake regions. Mueller-Vahl et al. [4] carried out wind-tunnel measurements of the NACA 63 (3) -618 airfoil equipped with triangular VGs. Their research implied that decreasing the spanwise spacing of VGs could not only delay the onset of static stall, but also cause a high drag penalty. Baldacchino et al. [5] systematically studied the effect of VG parameters on the aerodynamic performance of DU97-W-300 airfoil using wind-tunnel experiments. They found that the vane height and chordwise location of VGs are the main factors in the airfoil performance. The chordwise location plays a significant role in the post-stall behavior of airfoil. Positioning VGs too downstream can result in an early abrupt stall, because VGs become very prone to be submerged in the separation zones [4,5]. Wang et al. [6] studied the effect of rectangular VGs on the performance of NREL S809 airfoil by URANS simulations. Compared to single-row VGs, they found double-row VGs to further suppress the flow separation and then further delay the static stall.
RANS-based simulations of airfoil flow with VGs are by far the most common, although some researches were performed by highly expensive DNS/LES-type simulations [7,8]. Nevertheless, RANS methods also need a large quantity of computational cost, because the boundary layer with VGs requires adequately fine resolution. To reduce the cost of fully resolved RANS method, the VG modelling is often simplified [9]. The idea is to add a flow-dependent forcing term to the momentum equations based on the thin airfoil theory. This method can successfully predict the aerodynamic performances of both the airfoil and blade with VGs [10,11]. However, this simplified modelling has to calibrate the key coefficient first.
VGs have succeeded in aerospace engineering and have been practically applied in wind turbine engineering ( Figure 1). However, the HAWT blade flow controlled by VGs remains unclear. Most studies have focused on the two-dimensional steady airfoil flow controlled by VGs. In contrast, the blade flow is three-dimensional, rotational, and often becomes unsteady. The unsteady operating conditions are attributed to complicated environmental effects such as wind gust, turbulent inflow, and yaw misalignment [12,13]. The blade sections therefore undergo a time-varying AOA. If the AOA variation is dramatic enough, dynamic stall of the rotating blade will occur [14].
Energies 2020, 13, x FOR PEER REVIEW 2 of 12 velocimetry (PIV) and hot-film probes. They found the triangular VGs better than rectangular VGs in decreasing the drag penalty at low angles of attack (AOAs). The counter-rotating configuration was also found better than the co-rotating one in generating upwash and downwash wake regions. Mueller-Vahl et al. [4] carried out wind-tunnel measurements of the NACA 63(3)-618 airfoil equipped with triangular VGs. Their research implied that decreasing the spanwise spacing of VGs could not only delay the onset of static stall, but also cause a high drag penalty. Baldacchino et al. [5] systematically studied the effect of VG parameters on the aerodynamic performance of DU97-W-300 airfoil using wind-tunnel experiments. They found that the vane height and chordwise location of VGs are the main factors in the airfoil performance. The chordwise location plays a significant role in the post-stall behavior of airfoil. Positioning VGs too downstream can result in an early abrupt stall, because VGs become very prone to be submerged in the separation zones [4,5]. Wang et al. [6] studied the effect of rectangular VGs on the performance of NREL S809 airfoil by URANS simulations. Compared to single-row VGs, they found double-row VGs to further suppress the flow separation and then further delay the static stall. RANS-based simulations of airfoil flow with VGs are by far the most common, although some researches were performed by highly expensive DNS/LES-type simulations [7,8]. Nevertheless, RANS methods also need a large quantity of computational cost, because the boundary layer with VGs requires adequately fine resolution. To reduce the cost of fully resolved RANS method, the VG modelling is often simplified [9]. The idea is to add a flow-dependent forcing term to the momentum equations based on the thin airfoil theory. This method can successfully predict the aerodynamic performances of both the airfoil and blade with VGs [10,11]. However, this simplified modelling has to calibrate the key coefficient first.
VGs have succeeded in aerospace engineering and have been practically applied in wind turbine engineering ( Figure 1). However, the HAWT blade flow controlled by VGs remains unclear. Most studies have focused on the two-dimensional steady airfoil flow controlled by VGs. In contrast, the blade flow is three-dimensional, rotational, and often becomes unsteady. The unsteady operating conditions are attributed to complicated environmental effects such as wind gust, turbulent inflow, and yaw misalignment [12,13]. The blade sections therefore undergo a time-varying AOA. If the AOA variation is dramatic enough, dynamic stall of the rotating blade will occur [14]. Dynamic stall is characterized by the shedding and passage of a strong vortical disturbance over the suction surface, thereby causing a highly nonlinear fluctuating pressure field [17]. Dynamic stall often means the unsteady blade loads and there is noticeable aerodynamic hysteresis. These unsteady aerodynamic forces are directly linked to the structure failures, reduced turbine life, and increased operating maintenance.
Therefore, some methods were proposed to control dynamic stall, including aerodynamic blowing [18], trailing-edge flap [19], co-flow jet [20], and plasma actuator [21]. These existing ways can be classified as active control techniques, which will introduce auxiliary power equipment. This leads to a more complicated design of blades. Consequently, active control techniques are often limited to wind turbine blades [22]. In contrast, passive control techniques can also improve the wind turbine performance without external energy expenditure, among which VGs are very cost-effective. Dynamic stall is characterized by the shedding and passage of a strong vortical disturbance over the suction surface, thereby causing a highly nonlinear fluctuating pressure field [17]. Dynamic stall often means the unsteady blade loads and there is noticeable aerodynamic hysteresis. These unsteady aerodynamic forces are directly linked to the structure failures, reduced turbine life, and increased operating maintenance.
Therefore, some methods were proposed to control dynamic stall, including aerodynamic blowing [18], trailing-edge flap [19], co-flow jet [20], and plasma actuator [21]. These existing ways can be classified as active control techniques, which will introduce auxiliary power equipment. This leads to a more complicated design of blades. Consequently, active control techniques are often limited to wind turbine blades [22]. In contrast, passive control techniques can also improve the wind turbine performance without external energy expenditure, among which VGs are very cost-effective. Nevertheless, the effect of VGs on dynamic stall has been rarely investigated and hence is still poorly understood. Our previous works [23,24] demonstrated that VGs could effectively suppress the flow separation of oscillating wind turbine airfoil, thereby attenuating the aerodynamic hysteresis. In this regard, double-row VGs are more effective than single-row VGs. Our previous works only focused on the light dynamic stall controlled by VGs. However, the deep dynamic stall is often accompanied by stronger vortex motions and severer flow separation.
This work aims to investigate the effect of single-row and double-row VGs on the deep dynamic stall. The URANS method is used to identify the unsteady airfoil flow characteristics with and without VGs. The aerodynamic hysteresis loops, flow structures, and boundary-layer velocity profiles are analyzed in detail to reveal the effect of VGs on deep dynamic stall. Figure 2 illustrates the geometry of VGs used on the NREL S809 airfoil. The main VG parameters are h = 5 mm, d/h = 3.5, D/h = 7, L/h = 3, and β = 18 • , based on the VG design methodology [4,5]. Two chordwise locations are considered: x VG /c = 15% (single-row); x VG /c = 15% and 40% (double-row).

Geometry and Mesh Generation
x VG is measured between the leading edges of airfoil and VGs, and c is the airfoil chord length.
Energies 2020, 13, x FOR PEER REVIEW 3 of 12 Nevertheless, the effect of VGs on dynamic stall has been rarely investigated and hence is still poorly understood. Our previous works [23,24] demonstrated that VGs could effectively suppress the flow separation of oscillating wind turbine airfoil, thereby attenuating the aerodynamic hysteresis. In this regard, double-row VGs are more effective than single-row VGs. Our previous works only focused on the light dynamic stall controlled by VGs. However, the deep dynamic stall is often accompanied by stronger vortex motions and severer flow separation.
This work aims to investigate the effect of single-row and double-row VGs on the deep dynamic stall. The URANS method is used to identify the unsteady airfoil flow characteristics with and without VGs. The aerodynamic hysteresis loops, flow structures, and boundary-layer velocity profiles are analyzed in detail to reveal the effect of VGs on deep dynamic stall. Figure 2 illustrates the geometry of VGs used on the NREL S809 airfoil. The main VG parameters are h = 5 mm, d/h = 3.5, D/h = 7, L/h = 3, and β = 18°, based on the VG design methodology [4,5]. Two chordwise locations are considered: xVG/c = 15% (single-row); xVG/c = 15% and 40% (double-row). xVG is measured between the leading edges of airfoil and VGs, and c is the airfoil chord length. The present numerical modelling is essentially the same as that in our previous works [23,25]. Table 1 gives the main features of the mesh. The computational mesh includes only one pair of VGs (Figure 3), and the translational periodic boundary condition is used on spanwise boundaries. The mesh dependency study has been done by the General Richardson Extrapolation method [23]. The selected mesh with single-row VGs has about 1000 cells on each VG vane. There are 200 × 190 × 80 points in the wrap-around, normal, and spanwise directions, respectively. The Reynolds number is 1 × 10 6 (i.e., c = 0.457 m and U0 = 33.68 m/s), where U0 is the freestream velocity.   The present numerical modelling is essentially the same as that in our previous works [23,25]. Table 1 gives the main features of the mesh. The computational mesh includes only one pair of VGs (Figure 3), and the translational periodic boundary condition is used on spanwise boundaries. The mesh dependency study has been done by the General Richardson Extrapolation method [23]. The selected mesh with single-row VGs has about 1000 cells on each VG vane. There are 200 × 190 × 80 points in the wrap-around, normal, and spanwise directions, respectively. The Reynolds number is 1 × 10 6 (i.e., c = 0.457 m and U 0 = 33.68 m/s), where U 0 is the freestream velocity.

URANS settings
Dynamic stall of the airfoil is obtained by sinusoidal pitch oscillation about the quarter-chord axis. The sliding mesh method [26] is used to simulate the dynamic motion of airfoil. The instantaneous AOA follows the sinusoidal variation: The reduced frequency is defined as k = πfc/U0. In this work, the deep dynamic stall condition of αm = 18.75°, A = 10.3°, and k = 0.078 is simulated following the wind-tunnel experiments [27].
In the mesh motion simulated by sliding mesh method, nodes rigidly move in a given dynamic zone, but the cells defined by these nodes will not deform. A sliding interface is also introduced to connect multiple cell zones. The sliding interface is updated and synchronized with the mesh motion to reflect the new positions. Therefore, the computational domain is divided into two subdomains. The inner is a rotating region, and the outer a stationary region. The interaction between these two regions is made through a cylindrical sliding interface at radius of 4.4c in this work.
The commercial software ANSYS/FLUENT 16.0 [26] is used to numerically solve the URANS equations. Table 2 provides the main URANS settings. Ekaterinaris and Platzer [28] found that the proper consideration of transitional flow effect can improve the predictive accuracy of aerodynamic hysteresis. Therefore, the turbulence is simulated by the SST k-ω eddy viscosity model [29] incorporated with the γ-Reθ transition model [30]. This turbulence modelling has been proven to be reliable in simulating the dynamic stall of wind turbine airfoils [31]. The time step is set to assure 540 steps over each cycle and 20 inner iterations per time step, based on our previous works [13]. Iterative convergence criterion is met by assuring the cycle-to-cycle force variations negligible.

Spatial Discretization
Third-order MUSCL convection scheme

Temporal Discretization
Bounded second-order implicit scheme

Pressure-Velocity Coupling
Coupled algorithm Time Steps Per Cycle 540 Inner Iterations 20

Validation of numerical modelling
Due to the lack of experimental data of dynamic stall with VGs (unsteady-controlled), present numerical modelling has been validated against two sets of available experimental data: steadycontrolled and unsteady-uncontrolled.
For steady-controlled data, the numerical modelling can reliably predict the pressure distributions of the DU97-W-300 airfoil with and without triangular VGs ( Figure 4). Figure 4 also demonstrates that VGs are effective in suppressing the trailing-edge separated flow at α = 15°, thereby leading to a high leading-edge suction and greatly increasing the lift coefficient Cl. VGs, however, have a marginal effect on the pressure distribution when the flow is fully attached at α = 10°.

URANS Settings
Dynamic stall of the airfoil is obtained by sinusoidal pitch oscillation about the quarter-chord axis. The sliding mesh method [26] is used to simulate the dynamic motion of airfoil. The instantaneous AOA follows the sinusoidal variation: The reduced frequency is defined as k = πfc/U 0 . In this work, the deep dynamic stall condition of α m = 18.75 • , A = 10.3 • , and k = 0.078 is simulated following the wind-tunnel experiments [27].
In the mesh motion simulated by sliding mesh method, nodes rigidly move in a given dynamic zone, but the cells defined by these nodes will not deform. A sliding interface is also introduced to connect multiple cell zones. The sliding interface is updated and synchronized with the mesh motion to reflect the new positions. Therefore, the computational domain is divided into two subdomains. The inner is a rotating region, and the outer a stationary region. The interaction between these two regions is made through a cylindrical sliding interface at radius of 4.4c in this work.
The commercial software ANSYS/FLUENT 16.0 [26] is used to numerically solve the URANS equations. Table 2 provides the main URANS settings. Ekaterinaris and Platzer [28] found that the proper consideration of transitional flow effect can improve the predictive accuracy of aerodynamic hysteresis. Therefore, the turbulence is simulated by the SST k-ω eddy viscosity model [29] incorporated with the γ-Re θ transition model [30]. This turbulence modelling has been proven to be reliable in simulating the dynamic stall of wind turbine airfoils [31]. The time step is set to assure 540 steps over each cycle and 20 inner iterations per time step, based on our previous works [13]. Iterative convergence criterion is met by assuring the cycle-to-cycle force variations negligible.

Validation of Numerical Modelling
Due to the lack of experimental data of dynamic stall with VGs (unsteady-controlled), present numerical modelling has been validated against two sets of available experimental data: steady-controlled and unsteady-uncontrolled.
For steady-controlled data, the numerical modelling can reliably predict the pressure distributions of the DU97-W-300 airfoil with and without triangular VGs (Figure 4). For unsteady-uncontrolled data, the obtained results also show a good agreement with the experimental data of the NREL S809 airfoil in light dynamic stall [23]. Moreover, Figure 5 suggests that the calculated aerodynamic hysteresis loops generally agree with the experimental data [27] and Johansen's CFD results [31] in deep dynamic stall. The hysteresis loops also show noticeable fluctuations at the high AOAs. This is due to the severe vortex shedding and passage over the suction surface. The unsteady aerodynamic forces are accurately predicted during the flow separation and flow reattachment processes. Consequently, present numerical modelling of the deep dynamic stall of the NREL S809 airfoil with VGs should be adequately correct.   Figure 6 shows the calculated aerodynamic hysteresis loops with and without VGs. During the upstroke process, the aerodynamic coefficients of the airfoil with single-row and double-row VGs were relatively close. VGs significantly delayed the onset of dynamic stall. The Cl of clean airfoil began to diverge from the linear regime at α = 16°, implying the start of flow separation (Figure 6a). However, the Cl with VGs well followed the linear theory until α = 22°. At the high AOAs, the strong dynamic stall vortex motion caused large fluctuations in the aerodynamic coefficients of the clean airfoil. Figure 6 also indicates that the degree of these fluctuations can be decreased by VGs.

Aerodynamic hysteresis loops
During the downstroke process, the aerodynamic coefficients showed a clear difference between single-row and double-row VGs. The downstroke process from αmax to αmin can be further divided into three parts: =15.216  For unsteady-uncontrolled data, the obtained results also show a good agreement with the experimental data of the NREL S809 airfoil in light dynamic stall [23]. Moreover, Figure 5 suggests that the calculated aerodynamic hysteresis loops generally agree with the experimental data [27] and Johansen's CFD results [31] in deep dynamic stall. The hysteresis loops also show noticeable fluctuations at the high AOAs. This is due to the severe vortex shedding and passage over the suction surface. The unsteady aerodynamic forces are accurately predicted during the flow separation and flow reattachment processes. Consequently, present numerical modelling of the deep dynamic stall of the NREL S809 airfoil with VGs should be adequately correct.
Energies 2020, 13, x FOR PEER REVIEW 5 of 12 For unsteady-uncontrolled data, the obtained results also show a good agreement with the experimental data of the NREL S809 airfoil in light dynamic stall [23]. Moreover, Figure 5 suggests that the calculated aerodynamic hysteresis loops generally agree with the experimental data [27] and Johansen's CFD results [31] in deep dynamic stall. The hysteresis loops also show noticeable fluctuations at the high AOAs. This is due to the severe vortex shedding and passage over the suction surface. The unsteady aerodynamic forces are accurately predicted during the flow separation and flow reattachment processes. Consequently, present numerical modelling of the deep dynamic stall of the NREL S809 airfoil with VGs should be adequately correct.   Figure 6 shows the calculated aerodynamic hysteresis loops with and without VGs. During the upstroke process, the aerodynamic coefficients of the airfoil with single-row and double-row VGs were relatively close. VGs significantly delayed the onset of dynamic stall. The Cl of clean airfoil began to diverge from the linear regime at α = 16°, implying the start of flow separation (Figure 6a). However, the Cl with VGs well followed the linear theory until α = 22°. At the high AOAs, the strong dynamic stall vortex motion caused large fluctuations in the aerodynamic coefficients of the clean airfoil. Figure 6 also indicates that the degree of these fluctuations can be decreased by VGs.

Aerodynamic hysteresis loops
During the downstroke process, the aerodynamic coefficients showed a clear difference between single-row and double-row VGs. The downstroke process from αmax to αmin can be further divided into three parts:  Figure 6 shows the calculated aerodynamic hysteresis loops with and without VGs. During the upstroke process, the aerodynamic coefficients of the airfoil with single-row and double-row VGs were relatively close. VGs significantly delayed the onset of dynamic stall. The C l of clean airfoil began to diverge from the linear regime at α = 16 • , implying the start of flow separation (Figure 6a). However, the C l with VGs well followed the linear theory until α = 22 • . At the high AOAs, the strong Energies 2020, 13, 2535 6 of 13 dynamic stall vortex motion caused large fluctuations in the aerodynamic coefficients of the clean airfoil. Figure 6 also indicates that the degree of these fluctuations can be decreased by VGs.


From αmax to α = 25°, double-row VGs quickly restored the decreases in Cl and Cd in comparison with single-row VGs. This suggests that the second-row VGs impacted greatly on the massive flow separation when the airfoil began to pitch down.  From α = 25° to α = 13°, the Cl with single-row VGs kept high but was accompanied by high hysteresis intensities of the Cd and Cm. In contrast, the Cl with double-row VGs decreased gradually at first and then increased slowly.  From α = 13° to αmin, single-row VGs produced considerable increases in hysteresis intensities. The second-row VGs significantly helped the Cl readjust to the linear regime, so that double-row VGs led to low hysteresis intensities.  Table 3 provides the dynamic-stall parameters extracted from the hysteresis loops in Figure 6. The definition of aerodynamic pitch damping ζCm is given by:  During the downstroke process, the aerodynamic coefficients showed a clear difference between single-row and double-row VGs. The downstroke process from α max to α min can be further divided into three parts:

•
From α max to α = 25 • , double-row VGs quickly restored the decreases in C l and C d in comparison with single-row VGs. This suggests that the second-row VGs impacted greatly on the massive flow separation when the airfoil began to pitch down. • From α = 25 • to α = 13 • , the C l with single-row VGs kept high but was accompanied by high hysteresis intensities of the C d and C m . In contrast, the C l with double-row VGs decreased gradually at first and then increased slowly. • From α = 13 • to α min , single-row VGs produced considerable increases in hysteresis intensities. The second-row VGs significantly helped the C l readjust to the linear regime, so that double-row VGs led to low hysteresis intensities. Table 3 provides the dynamic-stall parameters extracted from the hysteresis loops in Figure 6. The definition of aerodynamic pitch damping ζ Cm is given by: A high ζ Cm also implies a high torsional aeroelastic stability. If the ζ Cm decreases to a negative value, the amplitude of airfoil pitch will increase rapidly, and then the flutter occurs unless restrained [17]. Single-row and double-row VGs increased the C l,max of NREL S809 airfoil by 40% and 49%, respectively. This is because double-row VGs could further delay the onset of dynamic stall. Both single-row and double-row VGs dramatically reduced the C m,min by almost 70%. The reason is that VGs hindered the forward motion of the center of pressure with the trailing-edge flow separation effectively suppressed. Table 3 also indicates a large decrease of 64% in the ζ Cm of the airfoil with single-row VGs. Therefore, single-row VGs can reduce the torsional aeroelastic stability, thereby likely causing the airfoil flutter. In this regard, double-row VGs are better to only reduce the ζ Cm from 0.153 to 0.124. The clean airfoil flow also showed wider separation zones during the downstroke process than during the upstroke process. This manifests that the aerodynamic hysteresis was attributed to the retarded flow reattachment when the airfoil pitched down.

Flow Field Developments
During the upstroke process, both single-row and double-row VGs eliminated the TE separation vortex at α = 18.75 • . The C l of airfoil with VGs was therefore dramatically increased (Figure 6a). At α = 27.67 • (↑), there were three separation vortices on the upper surface of the clean airfoil. Two small separation vortices were located on the first half chord, and one large separation vortex was near the trailing edge. Furthermore, single-row and double-row VGs produced a fourth small TE separation vortex. This small vortex crowded out the large separation vortex, thereby leading to a high TE suction peak (Figure 8a). Surprisingly, the second-row VGs seemed to bring about an undesirable effect to reduce the TE suction peak.
During the downstroke process, the leading-edge (LE) and TE separation vortices shed into the wake alternately. At α = 27.67 • (↓), the airfoil flow field with single-row VGs was highly distorted, because the LE separation vortex was hardly attached to the surface. Consequently, the suction value with single-row VGs was greatly decreased, even lower than that of the clean airfoil (Figure 8b). Double-row VGs, however, suppressed the LE flow separation effectively, and hence kept a high suction airfoil flutter. In this regard, double-row VGs are better to only reduce the ζCm from 0.153 to 0.124. Figure 7 illustrates the flow field developments around the airfoil with and without VGs. Three AOAs of 9.83°, 18.75°, and 27.67° were chosen to represent the three typical degrees of flow separation: fully attached flow, trailing-edge (TE) separated flow, and massively separated flow, respectively. The clean airfoil flow also showed wider separation zones during the downstroke process than during the upstroke process. This manifests that the aerodynamic hysteresis was attributed to the retarded flow reattachment when the airfoil pitched down.   During the upstroke process, both single-row and double-row VGs eliminated the TE separation vortex at α = 18.75°. The Cl of airfoil with VGs was therefore dramatically increased (Figure 6a). At α = 27.67° (↑ ), there were three separation vortices on the upper surface of the clean airfoil. Two small separation vortices were located on the first half chord, and one large separation vortex was near the trailing edge. Furthermore, single-row and double-row VGs produced a fourth small TE separation vortex. This small vortex crowded out the large separation vortex, thereby leading to a high TE suction peak (Figure 8a). Surprisingly, the second-row VGs seemed to bring about an undesirable effect to reduce the TE suction peak.

Flow field developments
During the downstroke process, the leading-edge (LE) and TE separation vortices shed into the wake alternately. At α = 27.67° (↓ ), the airfoil flow field with single-row VGs was highly distorted, because the LE separation vortex was hardly attached to the surface. Consequently, the suction value with single-row VGs was greatly decreased, even lower than that of the clean airfoil (Figure 8b). Double-row VGs, however, suppressed the LE flow separation effectively, and hence kept a high suction on the first half chord. Additionally, both single-row and double-row VGs avoided the secondary separation vortex near the trailing edge.
At α = 18.75° (↓ ), a large LE separation vortex appeared in the clean airfoil flow. The decreased height of separation vortex and the downstream movement of the vortex core suggest that the clean airfoil flow began to reattach gradually. However, both single-row and double-row VGs caused a tertiary vortex. This separation vortex was even detached from the upper surface due to double-row VGs, so that the suction on the upper surface and the Cl of airfoil with double-row VGs were vastly reduced (Figure 8c). Interestingly, although double-row VGs were positioned on the upper side, they significantly affected the Cp distribution on the lower side. This could decrease the pressure difference between the upper and lower sides. At α = 9.83° (↓ ), Figure 7 also implies that the second-row VGs further accelerated the flow reattachment and hence resulted in a high LE suction peak.

Boundary-layer velocity profiles
Figures 9 and 10 show the non-dimensionalized streamwise velocity profiles when the AOA reached 27.67° during the upstroke and downstroke processes, respectively. Note that α = 27.67° means the airfoil flow fell into the deep dynamic-stall process. Boundary-layer velocity profiles at x/c = 10% and x/c = 75% can represent the LE and TE separation vortices, respectively.
During the upstroke process, the streamwise velocity profiles with single-row and double-row VGs were quite close (Figure 9). This is attributed to the similar flow fields around the airfoil with single-row and double-row VGs at α = 27.67° (↑ ) (Figure 7). Figure 9 also indicates that VGs increased the height of LE separation vortex and the severity of TE reverse flow. Nevertheless, the external flow was effectively accelerated by VGs, thereby producing a high suction value (Figure 8a).
During the downstroke process, the streamwise velocity profiles also showed a clear difference between single-row and double-row VGs (Figure 10). At x/c = 10%, although the boundary-layer thicknesses with double-row VGs and without VGs were close, the external flow was effectively accelerated due to double-row VGs (Figure 10a). This suggests that double-row VGs made the LE flow withstand a higher adverse pressure gradient. However, the boundary-layer thickness with single-row VGs was decreased, because the LE separation vortex seemed to be detached from the wall surface (Figure 7). At x/c = 75%, the boundary-layer thicknesses from high to low was in the sequence of single-row VGs, clean, and double-row VGs. This sequence also determined the severity of TE flow separation. Interestingly, double-row VGs could effectively counteract the adverse pressure gradient and then suppress the TE flow separation, but single-row VGs could not ( Figure  8b). This finding highlights the great impact of the second-row VGs during the downstroke process. At α = 18.75 • (↓), a large LE separation vortex appeared in the clean airfoil flow. The decreased height of separation vortex and the downstream movement of the vortex core suggest that the clean airfoil flow began to reattach gradually. However, both single-row and double-row VGs caused a tertiary vortex. This separation vortex was even detached from the upper surface due to double-row VGs, so that the suction on the upper surface and the C l of airfoil with double-row VGs were vastly reduced ( Figure 8c). Interestingly, although double-row VGs were positioned on the upper side, they significantly affected the C p distribution on the lower side. This could decrease the pressure difference between the upper and lower sides. At α = 9.83 • (↓), Figure 7 also implies that the second-row VGs further accelerated the flow reattachment and hence resulted in a high LE suction peak.   Figures 9 and 10 show the non-dimensionalized streamwise velocity profiles when the AOA reached 27.67° during the upstroke and downstroke processes, respectively. Note that α = 27.67° means the airfoil flow fell into the deep dynamic-stall process. Boundary-layer velocity profiles at x/c = 10% and x/c = 75% can represent the LE and TE separation vortices, respectively.

Boundary-layer velocity profiles
During the upstroke process, the streamwise velocity profiles with single-row and double-row VGs were quite close ( Figure 9). This is attributed to the similar flow fields around the airfoil with single-row and double-row VGs at α = 27.67° (↑ ) (Figure 7). Figure 9 also indicates that VGs increased the height of LE separation vortex and the severity of TE reverse flow. Nevertheless, the external flow was effectively accelerated by VGs, thereby producing a high suction value (Figure 8a).
During the downstroke process, the streamwise velocity profiles also showed a clear difference between single-row and double-row VGs (Figure 10). At x/c = 10%, although the boundary-layer thicknesses with double-row VGs and without VGs were close, the external flow was effectively accelerated due to double-row VGs (Figure 10a). This suggests that double-row VGs made the LE flow withstand a higher adverse pressure gradient. However, the boundary-layer thickness with single-row VGs was decreased, because the LE separation vortex seemed to be detached from the wall surface (Figure 7). At x/c = 75%, the boundary-layer thicknesses from high to low was in the sequence of single-row VGs, clean, and double-row VGs. This sequence also determined the severity of TE flow separation. Interestingly, double-row VGs could effectively counteract the adverse pressure gradient and then suppress the TE flow separation, but single-row VGs could not ( Figure  8b). This finding highlights the great impact of the second-row VGs during the downstroke process. During the upstroke process, the streamwise velocity profiles with single-row and double-row VGs were quite close (Figure 9). This is attributed to the similar flow fields around the airfoil with single-row and double-row VGs at α = 27.67 • (↑) (Figure 7). Figure 9 also indicates that VGs increased the height of LE separation vortex and the severity of TE reverse flow. Nevertheless, the external flow was effectively accelerated by VGs, thereby producing a high suction value (Figure 8a).

Conclusions
This paper gives a flow analysis of deep dynamic stall of the NREL S809 airfoil controlled by single-row and double-row VGs. VGs were fully resolved, and URANS simulations were conducted with the transitional SST k-ω eddy viscosity model. Based on this study, several conclusions were reached as follows.


Present numerical modelling can accurately predict the aerodynamic loads of both an airfoil with VGs and an airfoil undergoing deep dynamic stall.  Both single-row and double-row VGs postpone the flow separation from α = 16° to α = 22°, when the airfoil pitches up. Then, the Cl,max of airfoil with VGs is considerably increased beyond 40%.  Although single-row and double-row VGs produce an additional TE separation vortex, they can reduce the fluctuations in aerodynamic coefficients near the αmax.  Single-row VGs bring about a vast decrease in the Cl from 2.2 to 0.3 near the αmax when the airfoil begins to pitch down, implying severe dynamic-stall behaviors. Single-row VGs also seriously retard the flow reattachment near the αmin. Therefore, single-row VGs considerably reduce the ζCm by 67% and hence undermine the torsional aeroelastic stability of airfoil.  Double-row VGs can quickly restore the decrease in Cl and Cd near the αmax in comparison with single-row VGs. Double-row VGs also significantly help the Cl readjust to the linear regime with the flow reattachment effectively accelerated.  Double-row VGs can effectively counteract the adverse pressure gradient and hence suppress the TE flow separation during the downstroke process, but single-row VGs cannot. This explains the clear difference in aerodynamic responses between single-row and double-row VGs. This paper also provides a performance assessment of VGs in controlling highly unsteady aerodynamic loads of a wind turbine airfoil. This study may contribute to understanding the deep dynamic stall controlled by single-row and double-row VGs. Future work should concentrate on the effect of passive VGs on a rotating blade undergoing dynamic stall.
Author Contributions: C.Z. conceived of the research and wrote the manuscript. C.Z. and J.C. conducted the data collection. T.W. and W.Z. contributed technical guidance and revised the manuscript. All authors have read and agreed to the published version of the manuscript. During the downstroke process, the streamwise velocity profiles also showed a clear difference between single-row and double-row VGs (Figure 10). At x/c = 10%, although the boundary-layer thicknesses with double-row VGs and without VGs were close, the external flow was effectively accelerated due to double-row VGs (Figure 10a). This suggests that double-row VGs made the LE flow withstand a higher adverse pressure gradient. However, the boundary-layer thickness with single-row VGs was decreased, because the LE separation vortex seemed to be detached from the wall surface (Figure 7). At x/c = 75%, the boundary-layer thicknesses from high to low was in the sequence of single-row VGs, clean, and double-row VGs. This sequence also determined the severity of TE flow separation. Interestingly, double-row VGs could effectively counteract the adverse pressure gradient and then suppress the TE flow separation, but single-row VGs could not (Figure 8b). This finding highlights the great impact of the second-row VGs during the downstroke process.

Conclusions
This paper gives a flow analysis of deep dynamic stall of the NREL S809 airfoil controlled by single-row and double-row VGs. VGs were fully resolved, and URANS simulations were conducted with the transitional SST k-ω eddy viscosity model. Based on this study, several conclusions were reached as follows.
• Present numerical modelling can accurately predict the aerodynamic loads of both an airfoil with VGs and an airfoil undergoing deep dynamic stall. • Both single-row and double-row VGs postpone the flow separation from α = 16 • to α = 22 • , when the airfoil pitches up. Then, the C l,max of airfoil with VGs is considerably increased beyond 40%.

•
Although single-row and double-row VGs produce an additional TE separation vortex, they can reduce the fluctuations in aerodynamic coefficients near the α max .

•
Single-row VGs bring about a vast decrease in the C l from 2.2 to 0.3 near the α max when the airfoil begins to pitch down, implying severe dynamic-stall behaviors. Single-row VGs also seriously retard the flow reattachment near the α min . Therefore, single-row VGs considerably reduce the ζ Cm by 67% and hence undermine the torsional aeroelastic stability of airfoil. • Double-row VGs can quickly restore the decrease in C l and C d near the α max in comparison with single-row VGs. Double-row VGs also significantly help the C l readjust to the linear regime with the flow reattachment effectively accelerated. • Double-row VGs can effectively counteract the adverse pressure gradient and hence suppress the TE flow separation during the downstroke process, but single-row VGs cannot. This explains the clear difference in aerodynamic responses between single-row and double-row VGs.
This paper also provides a performance assessment of VGs in controlling highly unsteady aerodynamic loads of a wind turbine airfoil. This study may contribute to understanding the deep dynamic stall controlled by single-row and double-row VGs. Future work should concentrate on the effect of passive VGs on a rotating blade undergoing dynamic stall.
Author Contributions: C.Z. conceived of the research and wrote the manuscript. C.Z. and J.C. conducted the data collection. T.W. and W.Z. contributed technical guidance and revised the manuscript. All authors have read and agreed to the published version of the manuscript. Acknowledgments: The authors would like to express their gratitude to the conference chairs of SEGT 2019 for recommendation of this publication in Energies.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
↑ upstroke process ↓ downstroke process α angle of attack (AOA) α Clmax AOA of C l,max α m mean AOA α max maximum AOA α min minimum AOA β geometric vane inflow angle ζ Cm aerodynamic pitch damping A AOA amplitude c chord length C d drag coefficient C l lift coefficient C l,dec C l at α Clmax (↓) C l,max maximum C l C m pitching moment coefficient C m,dec C m at α Clmax (↓) C m,inc C m at α Clmax (↑) C m,min minimum (maximum nose-down) C m C p pressure coefficient D inter-vane spacing d intra-vane spacing f frequency of oscillation h vane height k reduced frequency L vane length S n normal distance away from the wall surface u streamwise velocity U 0 freestream velocity x chordwise location x VG chordwise location measured between the airfoil and VG leading edges