Analysis of the E ﬀ ect of Vortex Generator Spacing on Boundary Layer Flow Separation Control

: During the operation of wind turbines, ﬂow separation appears at the blade roots, which reduces the aerodynamic e ﬃ ciency of the wind turbine. In order to e ﬀ ectively apply vortex generators (VGs) to blade ﬂow control, the e ﬀ ect of the VG spacing ( λ ) on ﬂow control is studied via numerical calculations and wind tunnel experiments. First, the large eddy simulation (LES) method was used to calculate the ﬂow separation in the boundary layer of a ﬂat plate under an adverse pressure gradient. The large-scale coherent structure of the boundary layer separation and its evolution process in the turbulent ﬂow ﬁeld were analyzed, and the e ﬀ ect of di ﬀ erent VG spacings on suppressing the boundary layer separation were compared based on the distance between vortex cores, the ﬂuid kinetic energy in the boundary layer, and the pressure loss coe ﬃ cient. Then, the DU93-W-210 airfoil was taken as the research object, and wind tunnel experiments were performed to study the e ﬀ ect of the VG spacing on the lift–drag characteristics of the airfoil. It was found that when the VG spacing was λ / H = 5 ( H represents the VG’s height), the distance between vortex cores and the vortex core radius were approximately equal, which was more beneﬁcial for ﬂow control. The ﬂuid kinetic energy in the boundary layer was basically inversely proportional to the VG spacing. However, if the spacing was too small, the vortex was further away from the wall, which was not conducive to ﬂow control. The wind tunnel experimental results demonstrated that the stall angle-of-attack (AoA) of the airfoil with the VGs increased by 10 ◦ compared to that of the airfoil without VGs. When the VG spacing was λ / H = 5, the maximum lift coe ﬃ cient of the airfoil with VGs increased by 48.77% compared to that of the airfoil without VGs, the drag coe ﬃ cient decreased by 83.28%, and the lift-to-drag ratio increased by 821.86%.


Introduction
During the operation of wind turbines, flow separation occurs at the blade root, which reduces the wind turbine efficiency [1]. Flow control technologies, such as VGs [2], plasma flow control [3], microflaps [4], microtabs [5], blowing and suction [6], synthetic jets [7], and flexible walls [8], have been increasingly applied to the design or optimization of wind turbine blades, aiming to improve their aerodynamic performance. Barlas [9] and Johnson [10] compared these flow control methods, while Lin [11] and Wang [12] found that VGs are one of the most effective devices for improving the aerodynamic performance of blades.
In 1947, Taylor [13] was the first to apply VGs to wing flow control. The basic principle of flow separation control using VGs is that when the fluid passes through a VG, a concentrated vortex will be generated downstream. Under the action of the concentrated vortex, the high-energy fluid outside the boundary layer will enter the bottom of the boundary layer, increasing the mixing of the fluids inside and outside the boundary layer. This will increase the kinetic energy of the fluid at the bottom of the boundary layer and delay its separation [14,15]. Due to their simple structure and convenient installation, VGs have been widely applied in fields such as flow control [16][17][18] and heat transfer [19,20]. In 1995, Reuss et al. [21] conducted two-dimensional wind tunnel tests on the National Advisory Committee for Aeronautics (NACA) 4415 airfoil in the subsonic wind tunnel of the Aerospace Research Laboratory. The aim of these tests was to record the cross-sectional lift and bending moment characteristics under different airflow conditions. The effect of the airfoil leading-edge roughness and VGs on the aerodynamic performance of the airfoil was investigated. In 2006, Godard and Stanislas [22][23][24] experimentally studied the control effect of triangular VGs on the boundary layer flow separation on a curved plate. Their experimental set-up could change the installation angle of VGs and measure the wake structure of small wings. In 2010, Yang et al. [25] performed numerical simulations to study a blunt trailing edge airfoil with and without VGs, and the performance of the airfoil was analyzed. The interaction between vortices generated by the VGs and wake vortices or separation vortices was investigated, while a mechanism for suppressing the boundary layer separation using VGs was revealed. In 2011, Zhen et al. [26] studied the effect of VGs on the aerodynamic performance of airfoils using numerical simulations, where the Spalart-Allmaras (SA) turbulence model and the standard wall function were employed. Their study results demonstrated that when the VGs were close to the separation point, the airfoil had a high maximum lift coefficient, and the performance of the rectangular and curved edge VGs was better than that of the triangular VGs. Lishu et al. [27] designed two types of VG combinations and conducted wind tunnel experiments in order to investigate airfoil flow control in depth. Their results revealed that the reasonable combination of a flap and VGs can significantly improve the aerodynamic performance of the airfoil, providing an ideal form of combined control. In 2012, Delnero et al. [28] experimentally investigated the effects of triangular VGs on airfoil aerodynamic performance. The VGs were located 10% and 20% away from the airfoil leading edge. The results showed that VGs can improve the airfoil aerodynamic characteristics, and the effect is better when the VGs are arranged 20% away from the leading edge. In 2013, Velte and Hansen [29] performed wind tunnel experiments using three-dimensional particle image velocimetry to observe the effect of VGs on the near-stall flow behavior of the DU91-W2-250 airfoil. Their aim was to study the flow structures induced by the VGs on the wing. It was found that the VGs transfer the high momentum fluid to the near-wall region, and the variation of the velocity component direction in the boundary layer along the flow increases. In 2014, Sørensen et al. [30] conducted numerical simulations on the FFA-W3-301 and FA-W3-360 airfoils. The calculated results were compared with wind tunnel measurements, and the variation of lift and drag with the angle of attack was obtained. Although their method was not able to capture the flow details accurately, it can be used to qualitatively compare the effects of different VG settings on airfoil aerodynamic performance. In 2015, Manolesos and Voutsinas [31] employed passive VGs as a simple and effective method to delay or suppress airfoil flow separation, and carried out an experimental study on the optimized blades. Flow visualization of triangular VGs through three-dimensional particle image velocimetry was carried out. Overall, the results demonstrated that the maximum lift increased by 44%, while the drag increased by 0.2%. Lin et al. [17] experimentally evaluated a boundary layer separation control method using miniature surface VGs. Wind tunnel tests were carried out in a low-turbulence pressure tunnel at the NASA (The National Aeronautics and Space Administration) Langley research center. The measured parameters included lift, drag, surface pressure, wake profile, and surface heat flux fluctuations. The results indicated that small VGs can effectively reduce the flow separation in the boundary layer of the flaps. Gao et al. [32] performed computational fluid dynamics (CFD) simulations to study the effect of the physical characteristics and magnitudes of flow from VGs on the aerodynamic performance of the DU97-W-300 blunt-trailing-edge wing. The effect of the VG size was analyzed in terms of the trailing edge thickness and the VG length. The results demonstrated that the drag loss was more sensitive to the increase in the VG's height than lift, while Appl. Sci. 2019, 9, 5495 3 of 22 an increase in the VG's length had a negative impact on both lift and drag. Analysis of the streamlines and vortices further revealed the flow field characteristics in the wake region. In 2016, Omar et al. [33] experimentally studied the NACA 4415 airfoil equipped with VGs. The results showed that triangular VGs were more suitable for controlling the boundary layer separation, while micro VGs can effectively control the flow. Zhang et al. [34] studied the effect of VGs on the aerodynamic performance of thick blades. In particular, the effects of height and chord position of VGs and double-row VGs on airfoil aerodynamic performance were investigated. The results demonstrated that after the VGs were installed on the airfoil leading edge, the maximum lift coefficient increased. In addition, the lift coefficient of the double-row VGs could be further improved at high angles of attack due to the different installation positions and sizes. In 2017, Wang et al. [12] studied the aerodynamic performance of the S809 airfoil with and without VGs using CFD simulations. The results showed that the VGs can effectively improve the aerodynamic performance of the S809 airfoil, reduce the thickness of the boundary layer, and delay stall. The double-row VG arrangement demonstrated a good performance in controlling flow separation, which further improved the aerodynamic performance of the S809 airfoil. Martínez-Filgueira et al. [35] studied the motion of micro VGs in the boundary layer of a flat plate using numerical simulations. The installation angle of the VGs was 18.5 • . The trajectory and size of the vortices on the plate were analyzed. Brüderlin et al. [36] improved airfoil aerodynamic efficiency by means of passive VGs placed on the control surface. First, the aerodynamic characteristics of the airfoil were simulated using a Reynolds-averaged Navier-Stokes solver (RANS). Then, through the parameterization study of different VG sizes, the effectiveness of the VGs under certain flow conditions were studied. In 2018, Gageik et al. [37] used a numerical visualization method to investigate whether appropriate micro VGs can suppress the pressure wave. The flow around the VGs was characterized using numerical schlieren visualization. Baldacchino et al. [38] performed wind tunnel experiments on the DU97-W-300 airfoil, where the oil flow visualization method was used to verify the effect of VGs on airfoil stall separation suppression. The results demonstrated that VGs can delay stall, while under a stall condition, the presence of VGs increases the load fluctuation. In 2019, Kundu et al. [39] used the numerical simulation method to study the performance of the S1210 airfoil with VGs under improved trailing edge conditions. The results showed that by adding counter-rotating VGs near the airfoil trailing edge, the lift coefficient of the wing can be increased by 17% and the stall angle can be delayed from 10 • to 12 • . Moreover, a rounder and thicker trailing edge can increase the lift coefficient by 13.5%, thus improving the performance. Nowadays, VGs are becoming an important part of wind turbine blade design. However, the challenges involved in calculating the flow field around VGs have not yet been satisfactorily addressed. Mereu et al. [40] used the scale decomposition method to simulate the stall behavior of the DU97-W-300 airfoil at high Reynolds numbers. In addition, a detailed sensitivity analysis was performed to evaluate the effects of time integration, grid width resolution, and region width on the simulation accuracy. Using this method, the flow on the airfoil surface with VGs was numerically simulated. Manolesos et al. [41] proposed a numerical model of VGs for RANS simulations. The performance and application of the VG numerical model were also investigated. A wind turbine impeller was taken as an example to study the reliability of the model. Without considering the mesh-induced errors, the results showed that the vortices generated in the numerical simulation were weaker than those in the fully analytical calculation.
The technical development of VGs as a means of flow control has been analyzed since the 1990s. Although many scholars have investigated the flow control mechanism of VGs, in the practical application of VGs, there are still some problems to be solved. When VGs are used to control the flow of wind turbine blades, an optimal VG spacing needs to be determined since there is flow interference between VGs and the spacing between VGs may affect the development and dissipation of concentrated vortices. The effect of VG spacing on vortex characteristics and wind turbine airfoil aerodynamic characteristics has not yet been reported. Therefore, in this paper, the effect of VG spacing on vortex characteristics and the aerodynamic characteristics of special airfoils for wind turbines is investigated using numerical simulations and wind tunnel experiments.

Large Eddy Simulation
A large eddy simulation (LES) is a spatial average of turbulent fluctuations, that is, large-scale vortices and small-scale vortices are separated using a filtering function [42]. More specifically, large-scale vortices are simulated directly, while small-scale vortices are closed using models. The theoretical basis of an LES is the existence of inertial sub-scale vorticity in turbulent flows with high Reynolds numbers. The vorticity of such scale has statistical isotropy. It transfers the energy of vortices with an energy scale to vortices with a dissipation scale. After filtering the Navier-Stokes equation with a filter function, the governing equation of the LES method for an incompressible flow can be derived: where σ ij is the stress tensor and τ ij is the subgrid stress. In this paper, the Smagorinsky-Lilly model was used for the subgrid model. This model was first proposed by Smagorinsky [43]. The eddy viscosity can be defined as follows: where L s is the mixed length of the grid and s 2s ij s ij , L s = min kd, C s V 1 3 , where k is a constant, d is the nearest distance from the wall, C s is the Smagorinsky constant (in this study, C s = 0.1), and V is the volume of the calculation unit.

Mesh Development
The ICEM software (ANSYS 16.0, Canonsburg, PA, USA) was employed to mesh the computational domain. The number of mesh nodes corresponding to the length, width, and height of the computational domain was 220 × 120 × 100, respectively, and the total number of elements was 2,640,000. The height of the first layer mesh was 0.001 mm and Y + < 1, which met the calculation requirements. The mesh around the VGs consisted of 50 nodes along the length direction and 40 nodes along the height direction. The mesh distribution can be seen in Figure 1.
The applied boundary conditions were the velocity inlet and pressure outlet. The upper boundary of the computational domain was set as an Euler wall, the bottom boundary and the VGs were set as adiabatic non-slip walls, and the two sides of the computational domain were set as symmetrical boundaries. The dimensions of the computational domain are shown in Figure 2. A turbulent velocity boundary condition was set at the inlet, where the inlet turbulence velocity type is given as follows: where u is the local velocity, δ is the thickness of the boundary layer, U is the main velocity, and y is the normal height of the wall.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 24 The applied boundary conditions were the velocity inlet and pressure outlet. The upper boundary of the computational domain was set as an Euler wall, the bottom boundary and the VGs were set as adiabatic non-slip walls, and the two sides of the computational domain were set as symmetrical boundaries. The dimensions of the computational domain are shown in Figure 2. A turbulent velocity boundary condition was set at the inlet, where the inlet turbulence velocity type is given as follows: where u is the local velocity, δ is the thickness of the boundary layer, U is the main velocity, and y is the normal height of the wall.

Computational Setup
The FLUENT software (ANSYS 16.0, Canonsburg, PA, USA) was used for the simulation. The finite volume method was used to discretize the equations and the finite centered difference scheme was used to discretize the space. Because for LES computing, the SIMPLE (Semi-Implicit Method for Pressure Linked Equations) algorithm saves more computing resources than the PISO (Pressure Implicit with Splitting of Operators) algorithm, the pressure-velocity coupling was based on the SIMPLE algorithm in this paper. The time step was set as Δt = 3 × 10 −5 s. After statistical stability in the flow field was achieved, time-averaged statistics were performed. In particular, the six cycles after the flow field was stabilized were time-averaged with the fluid sweeping across the lower wall  The applied boundary conditions were the velocity inlet and pressure outlet. The upper boundary of the computational domain was set as an Euler wall, the bottom boundary and the VGs were set as adiabatic non-slip walls, and the two sides of the computational domain were set as symmetrical boundaries. The dimensions of the computational domain are shown in Figure 2. A turbulent velocity boundary condition was set at the inlet, where the inlet turbulence velocity type is given as follows: where u is the local velocity, δ is the thickness of the boundary layer, U is the main velocity, and y is the normal height of the wall.

Computational Setup
The FLUENT software (ANSYS 16.0, Canonsburg, PA, USA) was used for the simulation. The finite volume method was used to discretize the equations and the finite centered difference scheme was used to discretize the space. Because for LES computing, the SIMPLE (Semi-Implicit Method for Pressure Linked Equations) algorithm saves more computing resources than the PISO (Pressure Implicit with Splitting of Operators) algorithm, the pressure-velocity coupling was based on the SIMPLE algorithm in this paper. The time step was set as Δt = 3 × 10 −5 s. After statistical stability in the flow field was achieved, time-averaged statistics were performed. In particular, the six cycles after the flow field was stabilized were time-averaged with the fluid sweeping across the lower wall

Computational Setup
The FLUENT software (ANSYS 16.0, Canonsburg, PA, USA) was used for the simulation. The finite volume method was used to discretize the equations and the finite centered difference scheme was used to discretize the space. Because for LES computing, the SIMPLE (Semi-Implicit Method for Pressure Linked Equations) algorithm saves more computing resources than the PISO (Pressure Implicit with Splitting of Operators) algorithm, the pressure-velocity coupling was based on the SIMPLE algorithm in this paper. The time step was set as ∆t = 3 × 10 −5 s. After statistical stability in the flow field was achieved, time-averaged statistics were performed. In particular, the six cycles after the flow field was stabilized were time-averaged with the fluid sweeping across the lower wall of the flat plate as a cycle. Figure 3 illustrates the structural parameters of VGs, where l is the VG length, H is the VG height, S is the distance between VGs, and λ is the spacing of the VGs (see Figure 3 for the difference between S and λ). The present study focused on the effect of the VG spacing on flow separation control. Table 1 lists the geometric parameters of the VGs, while the inflow velocity was set as U = 82 m/s, and the Reynolds number based on the VG height was about 3 × 10 4 .
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 24 of the flat plate as a cycle. Figure 3 illustrates the structural parameters of VGs, where l is the VG length, H is the VG height, S is the distance between VGs, and λ is the spacing of the VGs (see Figure  3 for the difference between S and λ). The present study focused on the effect of the VG spacing on flow separation control. Table 1 lists the geometric parameters of the VGs, while the inflow velocity was set as U = 82 m/s, and the Reynolds number based on the VG height was about 3 × 10 4 .    The pressure loss coefficient at the inlet and outlet of the calculation domain was defined as C ∆p , which can be calculated according to Equation (5), where P inlet is the total pressure at the inlet and P outlet the total pressure at the outlet. The pressure difference coefficient reflects the magnitude of pressure difference loss in the computational domain. For flow control, the smaller the pressure loss difference, the smaller the corresponding energy loss.
In order to verify the mesh independence, three meshes with different densities were developed with 1,320,000 (coarse), 1,980,000 (medium), and 2,640,000 (fine) elements. The Richardson extrapolation method was used to verify the convergence of the mesh, and the results are presented in Table 2, where RE is the value obtained using a Richardson extrapolation, p is the order of accuracy, and R is the ratio of errors. The convergence conditions are 0 < R < 1 (monotonic convergence), R < 0 (oscillation convergence), and 1 < R (divergence). As it can be seen, the R = 0.22 result indicates monotonic convergence. Therefore, in this paper, the finest mesh was used for the numerical simulations.

Experimental Setup
An experiment of an airfoil with VGs was carried out in the wind tunnel of North China Electric Power University (NCEPU), which contains a 2-D airfoil measurement section. The size of the test section is 1.5 × 3 × 4.5 m 3 , and the maximum wind speed of the wind tunnel is 62 m/s. The measurement section is equipped with 256-and 64-channel electronic scanning valves and a PXI (Platform PCI eXtensionsfor Instrumentation) data acquisition system. A pitot tube was installed on the wall of the upstream wind tunnel of the airfoil and provided a measurement of the total pressure of the inflow. The sidewall of the pitot tube was punctured to provide air holes, which were used to measure the static pressure of the incoming flow. There were 96 pressure taps in the middle section of the airfoil to measure the static pressure distribution on the surface of the airfoil. The static pressure coefficient (C p ) was the time-averaged result found with a sampling time interval of 0.1 seconds. The experimental sampling time was 10 seconds for each working condition. There were 111 total pressure piezometric holes and 5 static pressure piezometric holes (called the wake rake) downstream of the airfoil, which were used to calculate the drag coefficient of the airfoil. The geometric dimensions of the wind tunnel and test airfoils are shown in Figure 4.
The test airfoil used in the experiment is a special wind turbine airfoil: the DU93-W-210 airfoil. The airfoil was developed and designed by Delft University. The DU93-W-210 airfoil is very representative of the aerodynamic characteristics of wind turbine airfoils in general. The thickness of the airfoil was 21% of C, the chord length of the tested airfoil was 0.8 m, and the span of the airfoil was 1.5 m. The VGs were arranged at 0.2 times the chord length from the leading edge of the airfoil. Figure 5 shows the test airfoil with the VGs installed. The geometric dimensions of the VGs are shown in Table 3. The tunnel wind speed was set to 18.2 m/s, and the Reynolds number based on the chord of the airfoil was 1 × 10 6 . The AoA of the airfoil ranged from −10 • to 25 • .
The static pressure coefficient (Cp) was the time-averaged result found with a sampling time interval of 0.1 seconds. The experimental sampling time was 10 seconds for each working condition. There were 111 total pressure piezometric holes and 5 static pressure piezometric holes (called the wake rake) downstream of the airfoil, which were used to calculate the drag coefficient of the airfoil. The geometric dimensions of the wind tunnel and test airfoils are shown in Figure 4 The test airfoil used in the experiment is a special wind turbine airfoil: the DU93-W-210 airfoil. The airfoil was developed and designed by Delft University. The DU93-W-210 airfoil is very representative of the aerodynamic characteristics of wind turbine airfoils in general. The thickness of the airfoil was 21% of C, the chord length of the tested airfoil was 0.8 m, and the span of the airfoil was 1.5 m. The VGs were arranged at 0.2 times the chord length from the leading edge of the airfoil. Figure 5 shows the test airfoil with the VGs installed. The geometric dimensions of the VGs are shown in Table 3. The tunnel wind speed was set to 18.2 m/s, and the Reynolds number based on the chord of the airfoil was 1 × 10 6 . The AoA of the airfoil ranged from −10° to 25°.    The test airfoil used in the experiment is a special wind turbine airfoil: the DU93-W-210 airfoil. The airfoil was developed and designed by Delft University. The DU93-W-210 airfoil is very representative of the aerodynamic characteristics of wind turbine airfoils in general. The thickness of the airfoil was 21% of C, the chord length of the tested airfoil was 0.8 m, and the span of the airfoil was 1.5 m. The VGs were arranged at 0.2 times the chord length from the leading edge of the airfoil. Figure 5 shows the test airfoil with the VGs installed. The geometric dimensions of the VGs are shown in Table 3. The tunnel wind speed was set to 18.2 m/s, and the Reynolds number based on the chord of the airfoil was 1 × 10 6 . The AoA of the airfoil ranged from −10° to 25°.  Table 3. VG geometric parameters and arrangement. 20  17  5  5 3, 5, 7, 9   Figure 6 shows the wind tunnel test results of the DU93-W-210 airfoil without the VGs. NCEPU provided the test results in this paper. NPU is the wind tunnel test results from the Northwest Polytechnic University, and Delft is the wind tunnel test results from Delft University. By comparing the results of the three wind tunnels, it can be seen that the wind tunnel test results in this paper were close to those obtained by NPU and Delft before a stall AoA (8 • ) for the airfoil lift coefficient. After the stall AoA, the lift coefficient measured in the wind tunnel was slightly higher than that measured in the other two wind tunnels. For the drag coefficient, before the stall AoA, the drag coefficient matched the results of the other wind tunnels. However, after the stall AoA, the wind tunnel test results in this paper were closer to the Delft University wind tunnel test results. Since the aerodynamic characteristics of airfoils fluctuate after a stall, the lift and drag coefficients must be revised, and the mechanism of each correction method is different. However, the present study is focused on the aerodynamic characteristics of the airfoil before the stall. measured in the other two wind tunnels. For the drag coefficient, before the stall AoA, the drag coefficient matched the results of the other wind tunnels. However, after the stall AoA, the wind tunnel test results in this paper were closer to the Delft University wind tunnel test results. Since the aerodynamic characteristics of airfoils fluctuate after a stall, the lift and drag coefficients must be revised, and the mechanism of each correction method is different. However, the present study is focused on the aerodynamic characteristics of the airfoil before the stall.  Figure 7 shows time-averaged axial velocity contours on the symmetrical plane of the computational domain. As it can be observed, the channel expanded between 0-35H, which led to an adverse pressure gradient (APG), and a continuous decrease in fluid velocity. Under the action of the APG and wall viscosity, flow separation occurred at the lower wall, and a low-velocity zone appeared. The flow separation position was approximately at x1 ≈ 20H. Since the free-slip wall boundary condition was adopted on the upper wall of the computational domain, no flow separation occurred. In the 35-135H section, there was a zero pressure gradient (ZPG), and flow reattachment occurred at x2 ≈ 60H. The axial velocity contours with VGs revealed that the VGs inhibited flow separation in all four cases. When the VG spacing was λ/H = 7, a small low-speed region appeared at x ≈ 30H. When the VGs spacing was λ/H = 9, the area of the low-speed region at the same location increased. This showed that when λ/H = 9, the effectiveness of the concentrated vortices at this location became weak, while the hydrodynamic energy was already low under the inverse pressure gradient due to the large spacing between VGs. No low-velocity zone appeared at the other two installations with a smaller spacing.  Figure 7 shows time-averaged axial velocity contours on the symmetrical plane of the computational domain. As it can be observed, the channel expanded between 0-35H, which led to an adverse pressure gradient (APG), and a continuous decrease in fluid velocity. Under the action of the APG and wall viscosity, flow separation occurred at the lower wall, and a low-velocity zone appeared. The flow separation position was approximately at x 1 ≈ 20H. Since the free-slip wall boundary condition was adopted on the upper wall of the computational domain, no flow separation occurred. In the 35-135H section, there was a zero pressure gradient (ZPG), and flow reattachment occurred at x 2 ≈ 60H. The axial velocity contours with VGs revealed that the VGs inhibited flow separation in all four cases. When the VG spacing was λ/H = 7, a small low-speed region appeared at x ≈ 30H. When the VGs spacing was λ/H = 9, the area of the low-speed region at the same location increased. This showed that when λ/H = 9, the effectiveness of the concentrated vortices at this location became weak, while the hydrodynamic energy was already low under the inverse pressure gradient due to the large spacing between VGs. No low-velocity zone appeared at the other two installations with a smaller spacing. Figure 8 shows the distribution of the pressure coefficient (C p ) on the bottom wall of the computational domain. In the no VG case, the inflection point of C p appeared at x ≈ 20H, and the slope of the C p change decreased. In the cases where the boundary layer flow was controlled by the VGs, the C p on the wall increased in the APG section. Since no flow separation occurred, no alteration in the slope of C p change in the APG section was observed. The C p in most positions of the ZPG segment was basically 0. The difference between different VG spacings was primarily reflected in the turning position of the C p curve. In particular, the larger the VG spacing, the smaller the slope of the wall pressure drop and the worse the flow control effect. After the fluid passed through the expansion section, the C p increased. According to Figure 7, the flow separation occurred at x ≈ 20H. After the flow separation, the slope of the C p drop decreased, while after the flow reattachment, no significant changes in C p were observed in the ZPG section. Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 24      Although in some simple flows, the existence of vortices can be determined by intuition and visualization, in three-dimensional viscous flows, especially in complex flows, a large number of experimental or numerical simulation data are required to show the structure, evolution, and interaction of vortices. Therefore, it is necessary to give an objective discrimination of the vortices' criteria. At present, the Q criterion is more commonly used, where the region of Q > 0 is defined as a vortex, that is Ω 2 > E 2 . Its physical meaning is that the fluid rotation in the region of the vortex plays a leading role in comparing the strain rate. The specific formula is as follows:

Analysis of the CFD Calculation Results
where E 2 = e ij e ji and Ω 2 = Ω ij Ω ji = 1 2 |ω| 2 , e ij , Ω ij are the strain rate tensor and the vorticity tensor, respectively.
In this study, the Q criterion was used to identify the three-dimensional vortices in the flow field. Figure 9 shows the Q isosurfaces colored according to the velocity in the flow field (Q = 1 × 10 5 ). Many disordered vortices were observed in the turbulent flow field without VGs. However, typical turbulent coherent structures could be distinguished. In the initial flow stage, under the fluid viscosity and shear stress, the vortices were rolled up near the wall (as indicated by A). Due to self-rotation of the vortices, the radius of the vortex cores increased gradually downstream. At the same time, disturbed by the three-dimensional flow field, the vortices fluctuated extensively, twisted, and protruded upward along the normal direction. The protrusion was induced by the self-rotation of the spanning vortices, which intensified the protrusion height (as indicated by B). Under the action of a shear flow, the vortices continued to stretch and lift, and gradually developed into closed upper end hairpin vortices with an opening at the lower end (as indicated by C). Moving backwards, the legs of the vortices extended and the head of the hairpin vortices became larger, forming a local high-shear layer. Under the action of the mainstream fluid, the vortex was stretched along the flow direction into a tube (as indicated by D). This indicated that in turbulent flows, the average length of the vortex tube was generally increased. The energy transport process from larger to smaller vortices and the viscous dissipation into heat is called the cascade principle.
plays a leading role in comparing the strain rate. The specific formula is as follows: Ω are the strain rate tensor and the vorticity tensor, respectively.
In this study, the Q criterion was used to identify the three-dimensional vortices in the flow field. Figure 9 shows the Q isosurfaces colored according to the velocity in the flow field (Q = 1 × 10 5 ). Many disordered vortices were observed in the turbulent flow field without VGs. However, typical turbulent coherent structures could be distinguished. In the initial flow stage, under the fluid viscosity and shear stress, the vortices were rolled up near the wall (as indicated by A). Due to self-rotation of the vortices, the radius of the vortex cores increased gradually downstream. At the same time, disturbed by the three-dimensional flow field, the vortices fluctuated extensively, twisted, and protruded upward along the normal direction. The protrusion was induced by the self-rotation of the spanning vortices, which intensified the protrusion height (as indicated by B). Under the action of a shear flow, the vortices continued to stretch and lift, and gradually developed into closed upper end hairpin vortices with an opening at the lower end (as indicated by C). Moving backwards, the legs of the vortices extended and the head of the hairpin vortices became larger, forming a local high-shear layer. Under the action of the mainstream fluid, the vortex was stretched along the flow direction into a tube (as indicated by D). This indicated that in turbulent flows, the average length of the vortex tube was generally increased. The energy transport process from larger to smaller vortices and the viscous dissipation into heat is called the cascade principle.  Figure 10 shows Q isosurfaces colored according to the velocity in the flow field (Q = 1 × 10 5 ) as a result of the VG control. First, by comparing Figures 9 and 10, it can be found that without the VG control, many irregular vortices were generated downstream of the computational domain, while this number was reduced when the flow was controlled by VGs. By comparing the vortices under different VG spacings, it can be found that with the increase of the VG spacing, the height of the vortices increased. As it can be seen in Figure 10b, the presence of VGs created longitudinal vortices,  Figure 10 shows Q isosurfaces colored according to the velocity in the flow field (Q = 1 × 10 5 ) as a result of the VG control. First, by comparing Figures 9 and 10, it can be found that without the VG control, many irregular vortices were generated downstream of the computational domain, while this number was reduced when the flow was controlled by VGs. By comparing the vortices under different VG spacings, it can be found that with the increase of the VG spacing, the height of the vortices increased. As it can be seen in Figure 10b, the presence of VGs created longitudinal vortices, and this fact organized the local vortex structure and further hindered the development and movement of the spanwise vortices. The vortices generated between two VGs rotated in opposite directions, thus the outward vortices of the VGs were "suppressed." On the other hand, the upward forces at the inner side of the VGs were toward the center. Under the action of the two concentrated vortices, the inner extended vortices were "raised" to form hairpin-like vortices, and as the flow developed downstream, the concentrated vortices gradually merged with the leg of the hairpin-like vortices. The legs of the concentrated and hairpin-like vortices played an important role in the exchange of kinetic energy and energy dissipation near the wall. When λ/H = 3, the centralized vortices generated by the two VGs were so close to each other that the spanning vortices inside the VGs were not "raised," but instead were fused together at a downstream location under a reverse pressure gradient, forming hairpin-like vortices. When λ/H = 5, the outer spreading vortices of the VGs were "suppressed," while the inner spreading vortices were "raised," forming regular hairpin-like vortices that moved downstream. When λ/H = 7 and λ/H = 9, no regular hairpin-like vortices were developed on the inner side due to the large VG spacing. Such disordered turbulent structures were not conducive to the energy exchange in the near-wall region, which became prone to flow separation.
vortices. The legs of the concentrated and hairpin-like vortices played an important role in the exchange of kinetic energy and energy dissipation near the wall. When λ/H = 3, the centralized vortices generated by the two VGs were so close to each other that the spanning vortices inside the VGs were not "raised," but instead were fused together at a downstream location under a reverse pressure gradient, forming hairpin-like vortices. When λ/H = 5, the outer spreading vortices of the VGs were "suppressed," while the inner spreading vortices were "raised," forming regular hairpin-like vortices that moved downstream. When λ/H = 7 and λ/H = 9, no regular hairpin-like vortices were developed on the inner side due to the large VG spacing. Such disordered turbulent structures were not conducive to the energy exchange in the near-wall region, which became prone to flow separation.  Figure 11 shows vorticity isolines and spatial streamlines at five cross-sectional locations downstream of the VGs. First, as it can be seen in Figure 11a, the concentrated vortices generated by the VGs scrolled around the vortex core center and rotated in spirals. When λ/H = 7, at x > 20H, the streamlines began to bend due to the counter-pressure gradient effect, and the flow became unsmooth. When λ/H = 9, at x > 20H, the streamlines were twisted, which indicated that the VGs were not able to adequately suppress the flow separation. In Figure 11b, it can be observed that with the increase of the VG spacing, the distance between the two vortices gradually increased. Downstream of x/H = 15, when λ/H = 7 and λ/H = 9, the shape of the vortices was irregular with no obvious vortex center, while the area surrounding the vortices increased. When λ/H = 3 and λ/H = 5, the two vortices were closer to each other, and the shape of the concentrated vortices was still intact. Since the concentrated vortices rotated upward at the inner side of the VGs and were subjected to upward forces, the highest vortex core height from the wall was observed at λ/H = 3. When the distance between the vortices was small enough, the upward rotation force made the vortices lift off the wall. This implies that in order to restrain the flow separation using VGs, the vortices need to be closer to the wall to give a stronger energy exchange effect in the near-wall area, and the more beneficial the flow control. When λ/H = 3, the spacing was somewhat small, which was not conducive to flow control.
unsmooth. When λ/H = 9, at x > 20H, the streamlines were twisted, which indicated that the VGs were not able to adequately suppress the flow separation. In Figure 11b, it can be observed that with the increase of the VG spacing, the distance between the two vortices gradually increased. Downstream of x/H = 15, when λ/H = 7 and λ/H = 9, the shape of the vortices was irregular with no obvious vortex center, while the area surrounding the vortices increased. When λ/H = 3 and λ/H = 5, the two vortices were closer to each other, and the shape of the concentrated vortices was still intact. Since the concentrated vortices rotated upward at the inner side of the VGs and were subjected to upward forces, the highest vortex core height from the wall was observed at λ/H = 3. When the distance between the vortices was small enough, the upward rotation force made the vortices lift off the wall. This implies that in order to restrain the flow separation using VGs, the vortices need to be closer to the wall to give a stronger energy exchange effect in the near-wall area, and the more beneficial the flow control. When λ/H = 3, the spacing was somewhat small, which was not conducive to flow control.  Figure 12 shows the relationship between the vortex core radius and vortex core spacing with flow distance, where r represents the vortex radius and Δz represents the distance between the two  Figure 12 shows the relationship between the vortex core radius and vortex core spacing with flow distance, where r represents the vortex radius and ∆z represents the distance between the two vortex cores. With the increase of the VG spacing, the ∆z increased gradually. The vortex core radius represents the radial range of action. Theoretically, when VGs are used to suppress flow separation, the optimal distance is when the vortex radius r is equal to the distance between the two vortex cores ∆z/2 (half the distance). However, the radius of the vortices increased along the flow direction, and the vortices slightly deviated toward the flow direction. Therefore, the radius of the vortices and spacing between vortices cannot be equal everywhere along the flow direction. When the distance between vortices was ∆z/2 < r, VGs hindered the development of vortices. When ∆z/2 > r, a region where the vortices could not act was developed, which led to the decrease of the added value of fluid kinetic energy in the boundary layer. By comparing the four different VG spacing results, it can be seen that, except at x/H = 1, the distance between the vortices ∆z was smaller than the vortex radius r at the other locations when λ/H = 3. This shows that this VG spacing hindered the development of vortices at most positions along the flow direction, and that the VG spacing was too small to facilitate flow control. When λ/H = 7 and 9, ∆z > r for each direction, which reduced the effective range of the vortices and was not conducive to flow control. When λ/H = 5, it was found that ∆z > r before x/H = 10, while after x/H = 10, ∆z < r. Consequently, a VG spacing of λ/H = 5 was more suitable for flow control. the optimal distance is when the vortex radius r is equal to the distance between the two vortex cores Δz/2 (half the distance). However, the radius of the vortices increased along the flow direction, and the vortices slightly deviated toward the flow direction. Therefore, the radius of the vortices and spacing between vortices cannot be equal everywhere along the flow direction. When the distance between vortices was Δz/2 < r, VGs hindered the development of vortices. When Δz/2 > r, a region where the vortices could not act was developed, which led to the decrease of the added value of fluid kinetic energy in the boundary layer. By comparing the four different VG spacing results, it can be seen that, except at x/H = 1, the distance between the vortices Δz was smaller than the vortex radius r at the other locations when λ/H = 3. This shows that this VG spacing hindered the development of vortices at most positions along the flow direction, and that the VG spacing was too small to facilitate flow control. When λ/H = 7 and 9, Δz > r for each direction, which reduced the effective range of the vortices and was not conducive to flow control. When λ/H = 5, it was found that Δz > r before x/H = 10, while after x/H = 10, Δz < r. Consequently, a VG spacing of λ/H = 5 was more suitable for flow control. Second-order spatial correlation is a statistical method commonly used to study turbulence. Two-point spatial correlation can reflect the spatial relationship of one or more pulsations in the turbulent flow field. This method can be used to analyze the correlation of the structural characteristics of the downstream VG vortices. The correlation function is defined as follows: where ui and uj are components of the pulsating velocity, and x and r are vectors. When i = j, this expresses the autocorrelation of a certain pulsation, and when i ≠ j, it means that the two pulsations are correlated. The autocorrelation diagram is shown in Figure 13. The above formula can be transformed into dimensionless correlation coefficients as follows: The correlation coefficient magnitude represents the statistical similarity between two random variables or the probability that the fluctuation takes the same sign at different positions in the flow field. Second-order spatial correlation is a statistical method commonly used to study turbulence. Two-point spatial correlation can reflect the spatial relationship of one or more pulsations in the turbulent flow field. This method can be used to analyze the correlation of the structural characteristics of the downstream VG vortices. The correlation function is defined as follows: where u i and u j are components of the pulsating velocity, and x and r are vectors. When i = j, this expresses the autocorrelation of a certain pulsation, and when i j, it means that the two pulsations are correlated. The autocorrelation diagram is shown in Figure 13.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 15 of 24 that they are completely unrelated. For uniform turbulence, the correlation coefficient is independent of the coordinate position, and related only to the correlation distance. In order to investigate the correlation of the turbulent momentum downstream of the VGs, the correlation characteristics of the turbulent momentum v′ along the flow direction at certain points at different heights on a pair of centrosymmetric planes of the VGs were compared.   (Figure 14a). The minimum correlation distance was 0.5 mm, the maximum correlation distance was 675 mm, and the VGs were placed at the 0 position. According to Figure  14b, the coherence coefficient at the reference point is 1, indicating that the reference point was correlated with itself. In addition, the larger the correlation coefficient near the reference point, the larger the correlation coefficient amplitude, and the higher the correlation between the fluctuation and the reference point. A comparison between the following four cases demonstrated that when The above formula can be transformed into dimensionless correlation coefficients as follows: The correlation coefficient magnitude represents the statistical similarity between two random variables or the probability that the fluctuation takes the same sign at different positions in the flow field. R ij = ±1 denotes that two random variables are completely correlated, while R ij = 0 denotes that they are completely unrelated. For uniform turbulence, the correlation coefficient is independent of the coordinate position, and related only to the correlation distance. In order to investigate the correlation of the turbulent momentum downstream of the VGs, the correlation characteristics of the turbulent momentum v along the flow direction at certain points at different heights on a pair of centrosymmetric planes of the VGs were compared. Figure 14 illustrates the autocorrelation coefficient curve R v v (r 1 ) for the calculation of the fluctuating velocity v at different positions on the central symmetry plane of the domain. Five reference points at different heights and different points along the flow direction were selected for coherence analysis (Figure 14a). The minimum correlation distance was 0.5 mm, the maximum correlation distance was 675 mm, and the VGs were placed at the 0 position. According to Figure 14b, the coherence coefficient at the reference point is 1, indicating that the reference point was correlated with itself. In addition, the larger the correlation coefficient near the reference point, the larger the correlation coefficient amplitude, and the higher the correlation between the fluctuation and the reference point. A comparison between the following four cases demonstrated that when λ/H = 3, the magnitude of the correlation coefficient was the largest at most locations, which indicated that at this VG spacing, the structure of turbulence in the downstream domain preserved a high organization. When λ/H = 7 and 9, the magnitude of the correlation coefficient was larger at most locations, while that of λ/H = 5 was closer to R v v (r 1 ) = 0 at most locations. The correlation between the turbulent vortices downstream of the VGs and the reference point was poor, indicating that this VG spacing (λ/H = 5) was more suitable for suppressing flow separation. The width between the correlation coefficient curve and the two intersections reflected the scale of the coherent structure in the downstream turbulent flow field. It can be observed that the scale of the vortices increased gradually as the flow developed downstream. Figure 15 shows the average velocity distribution at different positions along the downstream direction. <u>/U ∞ represents the ratio of the time-averaged velocity to the mainstream velocity and y/H is the ratio of the normal wall height to the VG height. It can be observed that the increase in the kinetic energy of the fluid at the bottom of the boundary layer corresponded to a decrease in the kinetic energy of the fluid above the boundary layer. The S-type velocity profile took the center of the vortex core as the symmetrical point, which is the main principle of using VGs to control flow separation. In particular, due to the rotation of the concentrated vortices, the high-energy fluid outside of the boundary layer was transferred into the bottom of the boundary layer, inhibiting the boundary layer separation. By comparing the velocity profiles developed under the four VG spacings, it can be seen that the kinetic energy of the fluid at the bottom of the boundary layer increased with the decrease of the VG spacing. However, when λ/H = 3, downstream of x/H = 15, the center of the vortex core was higher than the wall, and the high-energy fluid in the boundary layer was also higher than the wall height. When λ/H = 5, the kinetic energy at the bottom of the boundary layer was the largest.
between the turbulent vortices downstream of the VGs and the reference point was poor, indicating that this VG spacing (λ/H = 5) was more suitable for suppressing flow separation. The width between the correlation coefficient curve and the two intersections reflected the scale of the coherent structure in the downstream turbulent flow field. It can be observed that the scale of the vortices increased gradually as the flow developed downstream.   Figure 15 shows the average velocity distribution at different positions along the downstream direction. <u>/U∞ represents the ratio of the time-averaged velocity to the mainstream velocity and y/H is the ratio of the normal wall height to the VG height. It can be observed that the increase in the kinetic energy of the fluid at the bottom of the boundary layer corresponded to a decrease in the kinetic energy of the fluid above the boundary layer. The S-type velocity profile took the center of the vortex core as the symmetrical point, which is the main principle of using VGs to control flow separation. In particular, due to the rotation of the concentrated vortices, the high-energy fluid outside of the boundary layer was transferred into the bottom of the boundary layer, inhibiting the boundary layer separation. By comparing the velocity profiles developed under the four VG spacings, it can be seen that the kinetic energy of the fluid at the bottom of the boundary layer increased with the decrease of the VG spacing. However, when λ/H = 3, downstream of x/H = 15, the center of the vortex core was higher than the wall, and the high-energy fluid in the boundary layer was also higher than the wall height. When λ/H = 5, the kinetic energy at the bottom of the boundary  Figure 16 shows the relationship between the pressure loss coefficient CΔP and the VG spacing. The percentage decrease of CΔP with VGs compared to that without VGs (clean) can be observed. When VGs were installed, the CΔP in the calculation domain decreased compared to that without VGs. Among them, the CΔP was the smallest when λ/H = 5, and the maximum CΔP was reduced by 30.95%. When λ/H = 9, the CΔP was the largest, which was reduced by 2.37% compared to that without VGs. Therefore, λ/H = 5 was more suitable for flow control.   Figure 17 shows the lift-drag coefficient and lift-drag ratio curves of the airfoil with and without VGs. Figure 17a demonstrates the lift coefficient curve. It can be seen that the lift coefficients of the airfoil with and without VGs were almost the same for an angle-of-attack (AoA) below 8°. Above 8°, the VGs could effectively increase the lift of the airfoil. In Table 4, it can be seen that when λ/H = 9, the stall AoA of the airfoil was 16°, while under the other VG spacings, the stall AoA was 18°. The maximum lift coefficient was obtained when λ/H = 5, which was increased by 48.77% compared to that of the clean airfoil. When λ/H = 3, λ/H = 7, and λ/H = 9, the maximum lift coefficient increased by 45.67%, 48.32%, and 30.20%, respectively. Therefore, the maximum lift coefficient of the airfoil was the largest when λ/H = 5. On the other hand, the drag coefficient was the worst when λ/H = 9. In Table 3, it can be seen that when λ/H = 3, the drag coefficient of the airfoil decreased by 82.06% compared to that of the clean airfoil at an AoA of 18°. In addition, it decreased by 83.28% when λ/H = 5, by 78.43% when λ/H = 7, and by 13.06% when λ/H = 9. The lift-drag ratio of the airfoil increased by 741.03% compared to that of the clean airfoil at an AoA of 18° when λ/H = 3, by 821.86% when λ/H = 5, by 612.06% when λ/H = 7, and by 55.06% when λ/H = 9.  Figure 16 shows the relationship between the pressure loss coefficient C ∆P and the VG spacing. The percentage decrease of C ∆P with VGs compared to that without VGs (clean) can be observed. When VGs were installed, the C ∆P in the calculation domain decreased compared to that without VGs. Among them, the C ∆P was the smallest when λ/H = 5, and the maximum C ∆P was reduced by 30.95%. When λ/H = 9, the C ∆P was the largest, which was reduced by 2.37% compared to that without VGs. Therefore, λ/H = 5 was more suitable for flow control.  Figure 16 shows the relationship between the pressure loss coefficient CΔP and the VG spacing. The percentage decrease of CΔP with VGs compared to that without VGs (clean) can be observed. When VGs were installed, the CΔP in the calculation domain decreased compared to that without VGs. Among them, the CΔP was the smallest when λ/H = 5, and the maximum CΔP was reduced by 30.95%. When λ/H = 9, the CΔP was the largest, which was reduced by 2.37% compared to that without VGs. Therefore, λ/H = 5 was more suitable for flow control.  Figure 17 shows the lift-drag coefficient and lift-drag ratio curves of the airfoil with and without VGs. Figure 17a demonstrates the lift coefficient curve. It can be seen that the lift coefficients of the airfoil with and without VGs were almost the same for an angle-of-attack (AoA) below 8°. Above 8°, the VGs could effectively increase the lift of the airfoil. In Table 4, it can be seen that when λ/H = 9, the stall AoA of the airfoil was 16°, while under the other VG spacings, the stall AoA was 18°. The maximum lift coefficient was obtained when λ/H = 5, which was increased by 48.77% compared to that of the clean airfoil. When λ/H = 3, λ/H = 7, and λ/H = 9, the maximum lift coefficient increased by 45.67%, 48.32%, and 30.20%, respectively. Therefore, the maximum lift coefficient of the airfoil was the largest when λ/H = 5. On the other hand, the drag coefficient was the worst when λ/H = 9. In Table 3, it can be seen that when λ/H = 3, the drag coefficient of the airfoil decreased by 82.06% compared to that of the clean airfoil at an AoA of 18°. In addition, it decreased by 83.28% when λ/H = 5, by 78.43% when λ/H = 7, and by 13.06% when λ/H = 9. The lift-drag ratio of the airfoil increased by 741.03% compared to that of the clean airfoil at an AoA of 18° when λ/H = 3, by 821.86% when λ/H = 5, by 612.06% when λ/H = 7, and by 55.06% when λ/H = 9.  Figure 17 shows the lift-drag coefficient and lift-drag ratio curves of the airfoil with and without VGs. Figure 17a demonstrates the lift coefficient curve. It can be seen that the lift coefficients of the airfoil with and without VGs were almost the same for an angle-of-attack (AoA) below 8 • . Above 8 • , the VGs could effectively increase the lift of the airfoil. In Table 4, it can be seen that when λ/H = 9, the stall AoA of the airfoil was 16 • , while under the other VG spacings, the stall AoA was 18 • . The maximum lift coefficient was obtained when λ/H = 5, which was increased by 48.77% compared to that of the clean airfoil. When λ/H = 3, λ/H = 7, and λ/H = 9, the maximum lift coefficient increased by 45.67%, 48.32%, and 30.20%, respectively. Therefore, the maximum lift coefficient of the airfoil was the largest when λ/H = 5. On the other hand, the drag coefficient was the worst when λ/H = 9. In Table 3, it can be seen that when λ/H = 3, the drag coefficient of the airfoil decreased by 82.06% compared to that of the clean airfoil at an AoA of 18 •   In Figure 18, the surface pressure coefficient (Cp) distribution of the airfoil is presented. It can be seen that the curves of Cp distribution on the airfoil surface with and without VGs coincided at 0° and 5° AoA. At AoAs of 10°, 15°, and 18°, the VGs postponed the pressure plateau to some extent. However, when λ/H = 9, the effect of the VGs flow control was the worst, and the pressure plateau appeared earlier than in the other three cases. It can be observed that the Cp curves of the other three VG spacings basically coincided. At an AoA of 20°, the Cp curves and the position of the pressure platform of λ/H = 9 were similar to those of the clean airfoil, while the Cp curves of the other three VGs spacings were similar. At 22°, the Cp curves of λ/H = 3 and λ/H = 5 basically coincided with the clean airfoil, while λ/H = 7 exhibited the best effect. When the AoA was equal to 25°, the Cp curves of the airfoil with VGs coincided with those of the clean airfoil.  In Figure 18, the surface pressure coefficient (C p ) distribution of the airfoil is presented. It can be seen that the curves of C p distribution on the airfoil surface with and without VGs coincided at 0 • and 5 • AoA. At AoAs of 10 • , 15 • , and 18 • , the VGs postponed the pressure plateau to some extent. However, when λ/H = 9, the effect of the VGs flow control was the worst, and the pressure plateau appeared earlier than in the other three cases. It can be observed that the C p curves of the other three VG spacings basically coincided. At an AoA of 20 • , the C p curves and the position of the pressure platform of λ/H = 9 were similar to those of the clean airfoil, while the C p curves of the other three VGs spacings were similar. At 22 • , the C p curves of λ/H = 3 and λ/H = 5 basically coincided with the clean airfoil, while λ/H = 7 exhibited the best effect. When the AoA was equal to 25 • , the C p curves of the airfoil with VGs coincided with those of the clean airfoil. Figure 19 shows the curve of the total pressure coefficient distribution at the wake rake. It can be seen that at AoAs of α = 0 • , 5 • , and 10 • , for an airfoil with VGs, the total pressure loss at the wake rake was larger because the generators increased the thickness of the boundary layer and thus increased the viscous losses. When the AoA was 15 • and 18 • , the total pressure loss of the airfoil with VGs was smaller than that of the clean airfoil. Among the four VG spacings, the total pressure loss of the airfoil with VGs was the largest when λ/H = 9, and under the other three spacings, the total pressure loss was similar. At α = 20 • , the total pressure loss was the largest when λ/H = 9 and the smallest when λ/H = 7. At α = 22 • , the total pressure loss was the largest when λ/H = 7, which was larger than that of the clean airfoil, while the total pressure loss of λ/H = 7 was the smallest. When α = 25 • , the total pressure loss was the same, with or without VGs. appeared earlier than in the other three cases. It can be observed that the Cp curves of the other three VG spacings basically coincided. At an AoA of 20°, the Cp curves and the position of the pressure platform of λ/H = 9 were similar to those of the clean airfoil, while the Cp curves of the other three VGs spacings were similar. At 22°, the Cp curves of λ/H = 3 and λ/H = 5 basically coincided with the clean airfoil, while λ/H = 7 exhibited the best effect. When the AoA was equal to 25°, the Cp curves of the airfoil with VGs coincided with those of the clean airfoil.  Figure 19 shows the curve of the total pressure coefficient distribution at the wake rake. It can be seen that at AoAs of α = 0°, 5°, and 10°, for an airfoil with VGs, the total pressure loss at the wake rake was larger because the generators increased the thickness of the boundary layer and thus increased the viscous losses. When the AoA was 15° and 18°, the total pressure loss of the airfoil with VGs was smaller than that of the clean airfoil. Among the four VG spacings, the total pressure loss of the airfoil with VGs was the largest when λ/H = 9, and under the other three spacings, the total pressure loss was similar. At α = 20°, the total pressure loss was the largest when λ/H = 9 and the smallest when λ/H = 7. At α = 22°, the total pressure loss was the largest when λ/H = 7, which was larger than that of the clean airfoil, while the total pressure loss of λ/H = 7 was the smallest. When α = 25°, the total pressure loss was the same, with or without VGs.  Figure 19. Distribution of total pressure coefficient at the wake rake.

Conclusions
When VGs are used for the flow control of wind turbine blades, they are always installed under a certain arrangement. In this paper, the LES method was used to simulate the boundary layer separation flow with and without VGs in an adverse pressure gradient environment. The large-scale coherent structure of the boundary layer separation caused by the adverse pressure gradient and its evolution process in the turbulent flow field were analyzed. The effect of different VG spacings on the suppression of the boundary layer separation were compared through the distance between vortex cores, the kinetic energy of the fluid in the boundary layer, and the pressure loss coefficient. Subsequently, the effect of VG spacings on the lift-drag characteristics, airfoil surface pressure distribution, and airfoil pressure loss were investigated in a wind tunnel using the DU93-W-210 airfoil as the research object.
According to the relationship between the vortex core radius and the distance between the vortex cores, when the VG spacing was λ/H = 3, the distance between the vortex cores Δz/2 was smaller than the vortex core radius r in most positions along the flow direction. Therefore, this spacing inhibited the development of vortices. When λ/H = 7 and 9, Δz/2 was always larger than r. When λ/H = 5, at x/H = 10, the Δz/2 and the r were approximately equal. Thus, the λ/H = 5 spacing was more suitable for flow control. The increase in fluid kinetic energy in the boundary layer was inversely proportional to the VG spacing. However, when the VG spacing was λ/H = 3, the distance between vortices was small and the vortex core was far from the wall, resulting in the high-energy fluid in the downstream boundary layer being higher from the wall. In the near-wall region, the kinetic energy of the fluid was the largest when λ/H = 5, especially at the x/H = 35 cross-section. When λ/H = 5, the CΔP between the inlet and outlet of the computational domain was the smallest, which was by 30.95% lower than that in the computational domain without VGs.
In the wind tunnel experimental results, it was found that the stall AoA of the airfoil with VGs was increased by 10° compared to that of the clean airfoil. When λ/H = 5, the maximum lift coefficient of the airfoil with VGs increased by 48.77%, which was the largest increase. Furthermore, at an AoA of 18°, when λ/H = 5, the drag coefficient decreased by 83.28%, which was the largest reduction. At an AoA of 18°, the lift-drag ratio of the airfoil increased by 741.03% when λ/H = 3, Figure 19. Distribution of total pressure coefficient at the wake rake.

Conclusions
When VGs are used for the flow control of wind turbine blades, they are always installed under a certain arrangement. In this paper, the LES method was used to simulate the boundary layer separation flow with and without VGs in an adverse pressure gradient environment. The large-scale coherent structure of the boundary layer separation caused by the adverse pressure gradient and its evolution process in the turbulent flow field were analyzed. The effect of different VG spacings on the suppression of the boundary layer separation were compared through the distance between vortex cores, the kinetic energy of the fluid in the boundary layer, and the pressure loss coefficient. Subsequently, the effect of VG spacings on the lift-drag characteristics, airfoil surface pressure distribution, and airfoil pressure loss were investigated in a wind tunnel using the DU93-W-210 airfoil as the research object.
According to the relationship between the vortex core radius and the distance between the vortex cores, when the VG spacing was λ/H = 3, the distance between the vortex cores ∆z/2 was smaller than the vortex core radius r in most positions along the flow direction. Therefore, this spacing inhibited the development of vortices. When λ/H = 7 and 9, ∆z/2 was always larger than r. When λ/H = 5, at x/H = 10, the ∆z/2 and the r were approximately equal. Thus, the λ/H = 5 spacing was more suitable for flow control. The increase in fluid kinetic energy in the boundary layer was inversely proportional to the VG spacing. However, when the VG spacing was λ/H = 3, the distance between vortices was small and the vortex core was far from the wall, resulting in the high-energy fluid in the downstream boundary layer being higher from the wall. In the near-wall region, the kinetic energy of the fluid was the largest when λ/H = 5, especially at the x/H = 35 cross-section. When λ/H = 5, the C ∆P between the inlet and outlet of the computational domain was the smallest, which was by 30.95% lower than that in the computational domain without VGs.
In the wind tunnel experimental results, it was found that the stall AoA of the airfoil with VGs was increased by 10 • compared to that of the clean airfoil. When λ/H = 5, the maximum lift coefficient of the airfoil with VGs increased by 48.77%, which was the largest increase. Furthermore, at an AoA of 18 • , when λ/H = 5, the drag coefficient decreased by 83.28%, which was the largest reduction. At an AoA of 18 • , the lift-drag ratio of the airfoil increased by 741.03% when λ/H = 3, 821.86% when λ/H = 5, 612.06% when λ/H = 7, and 55.06% when λ/H = 9. Consequently, it can be deduced that, when the VG spacing was λ/H = 5, the comprehensive effect of airfoil flow control was the best. The optimization of the VG spacing was essentially the optimization of vortex scale and vortex core spacing. The vortex scale and vortex core spacing were related to the VGs installation angle, geometric height, and arrangement spacing. Distance between the trailing edge of the vortex generators in the same direction C ∆p Pressure loss coefficient at the inlet and outlet of the calculation domain P inlet Total pressure at the inlet P outlet Total pressure at the outlet RE The value obtained using the Richardson extrapolation R The ratio of errors p The order of accuracy C p Airfoil surface static pressure coefficient C l