Separated Flow Control of Small Horizontal-Axis Wind Turbine Blades Using Dielectric Barrier Discharge Plasma Actuators

: The ﬂow control over the blades of a small horizontal-axis wind turbine (HAWT) model using a dielectric barrier discharge plasma actuator (DBD-PA) was studied based on large-eddy simulations. The numerical simulations were performed with a high-resolution computational method, and the e ﬀ ects of the DBD-PA on the ﬂow ﬁelds around the blades were modeled as a spatial body force distribution. The DBD-PA was installed at the leading edge of the blades, and its impacts on the ﬂow ﬁelds and axial torque generation were discussed. The increase in the ratios of the computed, cycle-averaged axial torque reasonably agreed with that of the available experimental data. In addition, the computed results presented a maximum of 19% increase in the cycle-averaged axial torque generation by modulating the operating parameters of the DBD-PA because of the suppression of the leading edge separation when the blade’s e ﬀ ective angles of attack were relatively high. Thus, the suppression of the leading edge separation by ﬂow control can lead to a delay in the breakdown of the tip vortex as a secondary e ﬀ ect. The above-mentioned studies have demonstrated the applicability of the DBD-PA-based ﬂow control to wind turbine blades of di ﬀ erent sizes and ﬂow conditions. However, the e ﬀ ects of DBD-PAs on the three-dimensional ﬂow structures around the blades have not been adequately understood yet. The objectives of the current study are to investigate the control mechanisms of the separated ﬂow around the wind turbine blades of a HAWT and analyze the controlled ﬂow structures around the blades using a DBD-PA. To this end, LES of the unsteady separated ﬂow around the blades


Introduction
Wind turbines convert the kinematic energy of wind into electrical energy. A horizontal-axis wind turbine (HAWT) is a common type of wind turbine, and usually consists of three main components, namely, a rotor, generator, and surrounding structure. The rotor includes the blades, where the aerodynamics of each blade have been sufficiently studied for more efficient and stable power generation [1]. However, the flow around wind turbine blades is sometimes separated due to the environmental flow around them, and this separation often causes unsteadiness and vibrations. Hence, flow separation is one of the main issues that disrupts the stable and robust power generation. Thus, in the current study, we investigate the control of the separated flow over wind turbine blades using a dielectric barrier discharge plasma actuator (DBD-PA). Generally, a DBD-PA has a greater advantage in structure, energy consumption, size, and weight in comparison with conventional devices, such as jet blowing devices [2]. In addition, it can be easily installed on the surfaces of existing objects and mechanical devices. Figure 1a shows a schematic of a typical-type DBD-PA, which consists of two thin electrodes and a dielectric material. When an alternating-current (AC) voltage with high volts (an order of 10 3 ) and driving frequencies (an order of 10 3 -10 4 ) is applied between the two electrodes, atmospheric-pressure non-equilibrium plasma is generated, which produces a jet flow above the dielectric material [3,4]. The relationship between the velocity of the induced wall jet and the parameters of the DBD-PA was experimentally investigated. Moreover, the parameters of the DBD-PA include operational parameters, such as the input voltage, frequency, and waveform, and geometric parameters, such as the electrode height, length, and gap [5][6][7][8][9][10][11][12]. The results [5][6][7][8][9][10][11][12] demonstrate that if the electrodes of the DBD-PA are geometrically optimized, then the time-averaged maximum velocity of the wall jet can reach up to 7 m/s, and if the multi-frequency of modulation is adopted as the input waveform (we called it a burst modulation, see Figure 1b), then the induced flow shows a rich spectrum in the frequency characteristics. Figure 1b shows the waveform of input AC voltage in the case of burst modulation. In Figure 1b, V pp , f base , f burst , T on , and T burst represent peak-to-peak magnitude and base frequency of the input AC voltage, the burst frequency, the duration in which the actuation of the DBD-PA is on, and the time period of burst modulation (T burst = 1/f burst ), respectively.
With knowledge of the fundamental characteristics of DBD-PAs and the results of their usage [5][6][7][8][9][10][11][12], a DBD-PA-based flow control was applied to stalled stationary airfoils where massive flow separation occurs. Improvements in the aerodynamic performance of the stalled airfoils have been reported in previous investigations [7,8,[11][12][13][14][15][16][17]. For example, Sato et al. [13][14][15] presented the effective operational parameters of DBD-PAs based on the results of more than 200 large-eddy simulations (LES) and the mechanisms of suppression of separated flows from the leading edge of an airfoil under a low Reynolds number (O(10 4 )) condition. Moreover, the suppression of the separated flow from the leading edge of the airfoil using a DBD-PA with burst modulations at a chord-based Reynolds number of 10 5 [15,16] and 10 6 [15,17] was demonstrated based on the results of LES. Recently, Fujii [18] provided a comprehensive review based on previous efforts [13][14][15][16][17] and also presented a further possibility of airfoil-flow control using a DBD-PA-based flow control based on the results of LES and other related experiments.
Energies 2020, 13, x FOR PEER REVIEW 2 of 16 objects and mechanical devices. Figure 1a shows a schematic of a typical-type DBD-PA, which consists of two thin electrodes and a dielectric material. When an alternating-current (AC) voltage with high volts (an order of 10 3 ) and driving frequencies (an order of 10 3 -10 4 ) is applied between the two electrodes, atmospheric-pressure non-equilibrium plasma is generated, which produces a jet flow above the dielectric material [3,4]. The relationship between the velocity of the induced wall jet and the parameters of the DBD-PA was experimentally investigated. Moreover, the parameters of the DBD-PA include operational parameters, such as the input voltage, frequency, and waveform, and geometric parameters, such as the electrode height, length, and gap [5][6][7][8][9][10][11][12]. The results [5][6][7][8][9][10][11][12] demonstrate that if the electrodes of the DBD-PA are geometrically optimized, then the timeaveraged maximum velocity of the wall jet can reach up to 7 m/s, and if the multi-frequency of modulation is adopted as the input waveform (we called it a burst modulation, see Figure 1b), then the induced flow shows a rich spectrum in the frequency characteristics. Figure 1b shows the waveform of input AC voltage in the case of burst modulation. In Figure 1b, Vpp, fbase, fburst, Ton, and Tburst represent peak-to-peak magnitude and base frequency of the input AC voltage, the burst frequency, the duration in which the actuation of the DBD-PA is on, and the time period of burst modulation (Tburst= 1/fburst), respectively. With knowledge of the fundamental characteristics of DBD-PAs and the results of their usage [5][6][7][8][9][10][11][12], a DBD-PA-based flow control was applied to stalled stationary airfoils where massive flow separation occurs. Improvements in the aerodynamic performance of the stalled airfoils have been reported in previous investigations [7,8,[11][12][13][14][15][16][17]. For example, Sato et al. [13][14][15] presented the effective operational parameters of DBD-PAs based on the results of more than 200 large-eddy simulations (LES) and the mechanisms of suppression of separated flows from the leading edge of an airfoil under a low Reynolds number (O(10 4 )) condition. Moreover, the suppression of the separated flow from the leading edge of the airfoil using a DBD-PA with burst modulations at a chord-based Reynolds number of 10 5 [15,16] and 10 6 [15,17] was demonstrated based on the results of LES. Recently, Fujii [18] provided a comprehensive review based on previous efforts [13][14][15][16][17] and also presented a further possibility of airfoil-flow control using a DBD-PA-based flow control based on the results of LES and other related experiments. Recently, a DBD-PA-based flow control was applied to the blades of HAWTs [19][20][21][22][23][24][25][26][27]. Tables 1  and 2 show that each study evaluated the performance of the DBD-PA-based flow control technique while taking into consideration different configurations of HAWTs during different flow conditions. In Tables 1 and 2, there are several nondimensional parameters (i.e., the Reynolds number Rec, tipspeed ratio , reduced frequency k, Mach number M, and duty ratio (burst ratio, BR). They are defined as follows: (1)  Recently, a DBD-PA-based flow control was applied to the blades of HAWTs [19][20][21][22][23][24][25][26][27]. Tables 1  and 2 show that each study evaluated the performance of the DBD-PA-based flow control technique while taking into consideration different configurations of HAWTs during different flow conditions. In Tables 1 and 2, there are several nondimensional parameters (i.e., the Reynolds number Re c , tip-speed ratio λ, reduced frequency k, Mach number M, and duty ratio (burst ratio, BR). They are defined as follows: where ρ ∞ , µ ∞ , U ∞ , U tip , a ∞ , c, b, r root , and f rot , represent the density of air, the viscosity of air, the velocity of air, the velocity of blade tip, speed of sound of air, the chord length of blade, the radial length at the root of blade from the center of rotation, and the rotational frequency of blade, respectively, and subscript ∞ indicates the quantities in the freestream condition. Moreover, certain improvements in torque generation due to the DBD-PA-based flow control technique have been reported [19][20][21][22][23][24][25][26][27].
Of note, all the cases attached the DBD-PAs near the leading edges of the turbine blades and utilized the burst modulation with the duty cycle. Moreover, Matsuda et al. [19,20] and Mitsuo et al. [25] measured the change in the output power due to the variation of the yaw angle (Φ yaw ) of the freestream coming to the rotational plane of a HAWT. Their results showed that decreases in the average power output took place when the yaw angle was high. In addition, a rotating stall around a small-scale HAWT [26] was experimentally studied recently, and the smart rotor control [27] of DBD-PAs was characterized and assessed without considering the yaw angle of the freestream. The flow field around the blade was also visualized using a phase-locked particle image velocimetry. Their results showed that a DBD-PA placed along the span produces a chord-wise body force and has very little effect on flow separation, even with the burst modulation. By contrast, a DBD-PA installed along the blade chord generates a body force into the radial directions, and this case successfully mitigates the rotating stall [26]. The torque due to the aerodynamic drag was reduced by up to 22% at a tip-speed ratio of 3.7, suppressing the stall over the outboard on 50% of the blade [26]. Moreover, Greenblatt et al. [28] studied pulsed, leading edge DBD-PAs on a vertical-axis wind turbine (VAWT). The control of the dynamic stall of VAWT blades was successfully demonstrated, leading to up to a 38% increase in the turbine power output. The above-mentioned studies have demonstrated the applicability of the DBD-PA-based flow control to wind turbine blades of different sizes and flow conditions. However, the effects of DBD-PAs on the three-dimensional flow structures around the blades have not been adequately understood yet. The objectives of the current study are to investigate the control mechanisms of the separated flow around the wind turbine blades of a HAWT and analyze the controlled flow structures around the blades using a DBD-PA. To this end, LES of the unsteady separated flow around the blades Energies 2020, 13, 1218 4 of 16 and the flow induced by the installed DBD-PA at the leading edge of the blade were conducted. The outline of this paper is as follows. First, the model information, flow conditions, and operational conditions of a DBD-PA are explained. Second, the numerical methods adopted in this study are briefly described. Then, the results computed in the cases with/without flow control are presented, and the key observations of the flow control effects using a DBD-PA are discussed. Finally, the conclusion is drawn.

Flow Conditions and Wind Turbine Blade Model
The current study considers a wind turbine blade configuration and a flow condition that are based on a model reported by Mitsuo et al. [25]. The experimental study [25] was mainly designed to analyze the effects of flow control around a simple HAWT blade model using a DBD-PA. Figure 2a shows the configuration of the HAWT blade considered in the current study. The cross-sectional shape of the blade is a NACA0012 airfoil. The chord length of the blade c was 0.2 m, and its span length b was 1.0 m. The blade was installed with a geometrical angle of attack α g of 10 degrees with regard to the rotational plane, and it revolved around the x-axis. Its rotational frequency f rot was 174 rpm. The position of the blade was represented using the azimuth angle φ(t). For example, φ of 0 degree and 90 degrees corresponds to a blade with a 12 and 9 o'clock direction from the upstream viewpoint, respectively. The freestream velocity U ∞ was set as 10 m/s, and the tilt angle of the freestream was 40 degrees (i.e., yaw angle, Φ yaw ) with respect to the rotational axis. Of note, although this flow condition was consistent with that of a previous experiment [25], it may or may not often occur in practice. Thus, for simplification purposes, the current study did not include a nacelle and a pole in all the computations (see Figure 2a). All the nondimensional parameters were computed based on Equation (1), and Re c , λ, k, and M ∞ were equal to 1.3 × 10 5 , 2.369, 0.182, and 0.07, respectively.
Considering the aerodynamic analysis of the wind turbine blade, information on the instantaneous effective angle of attack α e (t) at an arbitrary blade span location is useful. Thus, α e (t) at the radial position of the blade from the center of the rotational point (r) was estimated using the following formula based on the blade element theory: where r, φ(t), u a , u v , Φ yaw , and ω represent the radial location of the blade from the center of the rotation, azimuth angle of the blade, flow velocity normal to the rotational plane, effective incoming flow at the cross-section of the blade, yaw angle, and angular frequency, respectively. Figure 2b plots the time histories of α e for the selected span positions. The effective angle of attack α e reduces when the azimuth angle φ varies from 0 degree to 180 degrees. The cycle-averaged effective angles of attack at 50% and 75% spans from the root were 28.7 and 15.8 degrees, respectively. In the stationary airfoil cases, the stall angle of the NACA0012 at Re c = 1.0 × 10 5 was approximately 11 degrees [29].

of 16
Energies 2020, 13, x FOR PEER REVIEW 5 of 16 the azimuth angle  varies from 0 degree to 180 degrees. The cycle-averaged effective angles of attack at 50% and 75% spans from the root were 28.7 and 15.8 degrees, respectively. In the stationary airfoil cases, the stall angle of the NACA0012 at Rec= 1.0 × 10 5 was approximately 11 degrees [29].

Governing Equations
Three-dimensional Navier-Stokes (NS) equations, nondimensionalized by the freestream density, freestream velocity, and chord length of the blade, were employed as the governing equations. In the nondimensional form, the governing equations are represented as follows: , where xi, ui, qi, , p, e, ij, Si, ij, , t, and Dc denote the nondimensional forms of the positional vector, velocity vector, heat flux vector, density, pressure, total energy per unit volume, stress tensor, body force vector, Kronecker delta, ratio of specific heats, time, and nondimensional parameter for the DBD-PA body force, respectively. The Prandtl number Pr is defined as follows:

Governing Equations
Three-dimensional Navier-Stokes (NS) equations, nondimensionalized by the freestream density, freestream velocity, and chord length of the blade, were employed as the governing equations. In the nondimensional form, the governing equations are represented as follows: and D c denote the nondimensional forms of the positional vector, velocity vector, heat flux vector, density, pressure, total energy per unit volume, stress tensor, body force vector, Kronecker delta, ratio of specific heats, time, and nondimensional parameter for the DBD-PA body force, respectively. The Prandtl number Pr is defined as follows: where κ ∞ and C p represent the thermal conductivity and constant pressure specific heat, respectively. The last terms on the right-hand side of Equations (6) and (7) represent the momentum and energy added to the unit volume by the DBD-PA body force model. The description of these terms is presented in the next subsection. The system is closed and the equation for the ideal gas of state is written as follows: where p and γ represent the static pressure and the specific-heat ratio.

DBD-PA Body Force Model
The effects of the DBD-PA body force were represented as D c S i and D c S i u k in the NS equations. There were several DBD-PA body force models [30][31][32][33][34], and this study adopted one of them [33]. Detailed information regarding the adopted DBD-PA body force model can be found in Aono et al. [33], and a brief description of the model is given below. Figure 3 shows the sketch of the modeled DBD-PA, contour of the magnitude of the nondimensional body force, and body force vectors near the exposed electrode. The spatial distribution of the body force was similar to that of the Suzen-Huang's model [32]. The body force vector in the component along the blade span was set to zero because the uniform geometry of the DBD-PA along the blade span was considered. The body force model rotated 90 degrees counterclockwise about the blade span axis and then attached to the blade surface at the leading edge of the blade.
where p and  represent the static pressure and the specific-heat ratio.

DBD-PA Body Force Model
The effects of the DBD-PA body force were represented as DcSi and DcSiuk in the NS equations. There were several DBD-PA body force models [30][31][32][33][34], and this study adopted one of them [33]. Detailed information regarding the adopted DBD-PA body force model can be found in Aono et al. [33], and a brief description of the model is given below. Figure 3 shows the sketch of the modeled DBD-PA, contour of the magnitude of the nondimensional body force, and body force vectors near the exposed electrode. The spatial distribution of the body force was similar to that of the Suzen-Huang's model [32]. The body force vector in the component along the blade span was set to zero because the uniform geometry of the DBD-PA along the blade span was considered. The body force model rotated 90 degrees counterclockwise about the blade span axis and then attached to the blade surface at the leading edge of the blade. The magnitude of the DBD-PA body force was represented by multiplying the nondimensional parameter (Dc) by the body force vectors (Si) in the NS equations (Equations (6) and (7)). The definition of Dc is as follows: , (2) zoom-up view near the exposed electrode.
The magnitude of the DBD-PA body force was represented by multiplying the nondimensional parameter (D c ) by the body force vectors (S i ) in the NS equations (Equations (6) and (7)). The definition of D c is as follows: where Q c and E i denote the electric charge and electric-field vector, respectively, and subscript ref denotes the reference value. Moreover, we set Q c,ref and E ref as the maximum values of Q c and E in the pre-computation of the model, respectively [32]. In addition, the temporal variation in the body force distribution was formulated by the multiplication of the body force vectors (S i ) and sin 2 (2πF base t).
Here, F base is the nondimensional base frequency of the input AC voltage (F base = f base c/U ∞ ), and was set at 170, which is higher than that of the fluctuation of flow. This model assumed that the temporal variation of the body force can be characterized as push/push [35]. The input waveform of the burst modulation was generated based on the nondimensional burst frequency (F + = f burst c/U ∞ ) and BR as shown in Figure 1b.

Note of Reliability of Approach and Adopted DBD-PA Body Force Model
In this subsection, the reliability of the approach adopted in this study and adopted DBD-PA body force model are noted. As described in Sections 3.1 and 3.2, in our approach, before conducting flow-control simulations, the partial differential equations of the Suzen-Huang's model [32] with a static boundary condition were numerically solved in order to calculate the body force vectors (S i ) that represent the maximum forcing. After that, we incorporated this electro-hydrodynamic body force vector into the right-hand side of the NS equations as the body force terms, with an assumption that the magnitude of the distribution body force was in proportion to the boundary condition. The reliability the approach and the body force model were investigated in our previous studies [13][14][15]33,34,36]. In particular, the flow induced by the DBD-PA under the quiescent air condition [33,34] and the flow structures around the airfoil at the chord-based Reynolds number of 63,000 [13][14][15] and the circular disk at the diameter-based Reynolds number of 5000 [36] were reported. Setting the appropriate range of nondimensional parameter D c was important in order to obtain the realistic range of time-averaged maximum induced velocity by the Suzen-Huang's model with the burst modulations under the quiescent air condition and the reliable results of controlled effects using the DBD-PA. Therefore, we carried out several pre-computations of the flow field induced by the DBD-PA in a quiescent field in order to find the appropriate range of D c . As a result, the maximum time-averaged induced wall-parallel velocities were close to 0.25U ∞ and 0.4U ∞ in the cases of burst modulation with D c of 0.5, a BR of 10% and F + of 1 and continuous actuation with D c of 0.125, respectively. From this analysis and considering previous efforts [13][14][15]33,34,36], we think that induced flow and effects of control using the DBD-PA on the blades discussed in this study are sufficiently reliable.

Numerical Methods
In this study, an in-house fluid dynamic solver (LANS3D [37]) was used, and the NS equations were solved in the generalized curvilinear coordinates (ξ, η, ζ). The spatial derivatives of the convective, viscous terms, symmetric-conservative metrics [38], and Jacobian [38] were evaluated with a sixth-order compact difference scheme [39]. In addition, near the computational boundaries, including the inner boundaries resulting from the domain decomposition, a second-order, one-sided, and central formulas at the first and second grid points off the boundaries were used to evaluate the spatial derivatives because of the unavailability of a complete stencil of the six-order compact difference scheme. A tenth-order filtering [39,40] was used with a filtering coefficient of 0.4 to stabilize the computation. As for the time integration, an alternating direction implicit symmetric Gauss-Seidel method [41,42] was used. Moreover, a backward second-order difference formula was utilized for time integration, and three sub-iterations [43] were adopted to ensure the accuracy of time integration. The computational time step was 3.5 × 10 −5 c/U ∞ , which corresponds to the maximum local Courant number that is approximately equal to 30. All the grids were rigidly rotated around the origin of the coordinate, and the rotations of the blades were simulated. The coherent structure model [44] was used as a subgrid scale model of the LES, and a parallel computation was conducted using message passing interfaces. All the simulations were carried out using 18,800 cores of the K computer, which is based on distributed memory architecture of more than 640,000 cores [45]. The subroutines related to the compact finite differencing scheme in LANS3D were tuned for the K computer [46]. Moreover, efficient evaluation techniques of the symmetric conservative metrics [47] were adopted to reduce the computational time.

Computational Grids and Boundary Conditions
This study used the overset grid technique [48] to handle the region of the body force and the three wind turbine blades. Here, a single grid is called a zone. Table 3 shows the number of grid points of each grid (zone), where the total number of grid points is approximately 700 million. Figure 4 shows the computational grids around the model. The grid system requires four single grids per wind turbine blade to represent the blade surface and DBD-PA, where these grids are composed of the blade-main, blade-tip, blade-root, and body-force subgrids. The body force distribution was computed based on the grid corresponding to the body force model region, and the body force was mapped to zones 6, 10, and 14. Two background grids surrounding the blade volume grids were employed, leading to the completion of the overset grid system. The length between the outer boundary and center of rotation was 50c. At the outflow boundaries, zero-gradient conditions were adopted, and at the inflow boundaries, the freestream velocity was assigned. As for the surface of the blade, no-slip and adiabatic-wall conditions were adopted. The minimum grid size in the direction normal to the airfoil surface was set as 1.369 × 10 −4 c. This grid spacing is so small that at least 20 grid points can exist within the laminar boundary layer around the leading edge of the blade. It was impossible for us to perform a Energies 2020, 13, 1218 8 of 16 grid sensitivity analysis in this study due to the required highly expensive computational resources. Therefore, the current study qualitatively discusses the effects of DBD-PAs on the aerodynamics of wind turbine blades and the applicability of DBD-PA-based flow control techniques to wind turbine blades. Although the grid convergence study is absent, the qualitative trend of an increase in the ratio of the cycle-averaged axial torque shows a reasonable agreement with the experimental data (see Figure 5) [25]. Furthermore, the used fluid dynamics solver in the current study was validated for flow control simulations around the airfoil in the Re c range of O(10 4 ) to O(10 6 ) [13][14][15][16][17].

Computational Cases
All the cases are tabulated in Table 4. In the current study, one computation without flow control using a DBD-PA and three computations with flow control using a DBD-PA with different operational conditions were carried out. The operational parameters of some of the controlled cases were set based on the experimental conditions of the previous study [25], and the duty cycle/BR was fixed to 0.1. The DBD-PA was switched on at 2 degrees or 0.33 degrees per 21 degrees of the azimuth angle () with F + of 1 and 6, respectively.

Computational Cases
All the cases are tabulated in Table 4. In the current study, one computation without flow control using a DBD-PA and three computations with flow control using a DBD-PA with different operational conditions were carried out. The operational parameters of some of the controlled cases were set based on the experimental conditions of the previous study [25], and the duty cycle/BR was fixed to 0.1. The DBD-PA was switched on at 2 degrees or 0.33 degrees per 21 degrees of the azimuth angle (φ) with F + of 1 and 6, respectively.

Data Processing
The computed data were averaged based on 72 phases (i.e., 5 degrees of the azimuth angle φ per phase), and a phase number of one meant that the instantaneous flow data around the blades with an azimuth angle between 0 degrees and 5 degrees was averaged. The results at the two phases, where a remarkable difference between the cases with and without the DBD-PA was found, and the maximum torque generation were highlighted. Moreover, α e was the same at both phases. The cycle-averaged axial torque T a was calculated by adding the cycle-averaged axial torque values of the blades. Two revolutions of the blades were also used to obtain the cycle-averaged value. Figure 5 shows an increase in the ratios of the cycle-averaged axial torque (∆T a ) generated by one blade due to the flow control of the DBD-PA, and the computed results indicate that the flow control using DBD-PAs installed at the leading edges of wind turbine blades can increase the axial torque generation by more than 10% in comparison with blades without flow control. The maximum ∆T a was approximately 19% when F + = 1 and D c = 0.5 (D c 0.5F + 1). When cases with the same D c were considered, a higher ∆T a was obtained with a higher F + (D c 0.125F + 6). This qualitative trend of the improvement in ∆T a with respect to the change of F + is similar to that reported in a previous experiment [25].

Effect of DBD-PA on the Cycle-Averaged Axial Torque Output
Energies 2020, 13, x FOR PEER REVIEW 10 of 16 improvement in Ta with respect to the change of F + is similar to that reported in a previous experiment [25].    Figure 6 presents the phase-averaged torque coefficient of the blade (C Ta ) as a function of the azimuth angle φ in the cases with and without control. The variation of C Ta with respect to the change of the azimuth angle φ is similar in all cases. Unlike the stationary airfoil-flow control cases [13][14][15][16][17], the effects of F + on the aerodynamic performance were not clearly identified. However, when the azimuth angle φ approached 90 degrees, the difference in the value of C Ta with and without flow control was clear, and it became minimum when C Ta reached its maximum value. The trend of the difference in the value of C Ta is most likely due to differences in the flow structures around the blades because of the variation in the effective angle of attack α e . Next, the flow structures around the blades at the selected two phases that correspond to the azimuth angle φ of 112.5 degrees (see Figure 6 (i) refers to near minimum axial torque generation) and 247.5 degrees (see Figure 6 (ii) refers to near maximum axial torque generation) were analyzed. as a function of F . Figure 6 presents the phase-averaged torque coefficient of the blade (CTa) as a function of the azimuth angle  in the cases with and without control. The variation of CTa with respect to the change of the azimuth angle  is similar in all cases. Unlike the stationary airfoil-flow control cases [13][14][15][16][17], the effects of F + on the aerodynamic performance were not clearly identified. However, when the azimuth angleapproached 90 degrees, the difference in the value of CTa with and without flow control was clear, and it became minimum when CTa reached its maximum value. The trend of the difference in the value of CTa is most likely due to differences in the flow structures around the blades because of the variation in the effective angle of attack e. Next, the flow structures around the blades at the selected two phases that correspond to the azimuth angle  of 112.5 degrees (see Figure 6 (i) refers to near minimum axial torque generation) and 247.5 degrees (see Figure 6 (ii) refers to near maximum axial torque generation) were analyzed.  Figure 7 shows the three-dimensional flow structures around the blades at two selected phases (= 112.5 degrees and 247.5 degrees, respectively). The isosurfaces of the second invariant of the velocity gradient tensor (Q-criterion [49], Q= 200) were visualized and colored by the flow velocity magnitude in the x-direction. In the first phase (Figure 7 (i)), the flow around blade 1 was separated in both cases. The suppression of the flow separation around blade 1 due to the flow control can be  Figure 7 shows the three-dimensional flow structures around the blades at two selected phases (φ = 112.5 degrees and 247.5 degrees, respectively). The isosurfaces of the second invariant of the velocity gradient tensor (Q-criterion [49], Q = 200) were visualized and colored by the flow velocity magnitude in the x-direction. In the first phase (Figure 7 (i)), the flow around blade 1 was separated in both cases. The suppression of the flow separation around blade 1 due to the flow control can be found in the region between 50% of the blade's span and the blade tip. By contrast, in the second phase (Figure 7 (ii)), the flow around blade 1 near the tip was weakly separated. Thus, the suppression of the flow separation due to the flow control was lesser compared with that in the first phase. Thus, in the first phase ( Figure 6 (i) and Figure 7 (i)), the increase in the axial torque generation from the flow control using the DBD-PA was more present in comparison with the second phase ( Figure 6 (ii) and Figure 7 (ii)).

Phase-Averaged Flow and Aerodynamic Characteristics
The flow visualization illustrates that the effects of the flow control using DBD-PAs on the flow structure around turbine blades vary with the span direction. Here, we further analyzed the relationship between the flow structures and the contribution of the axial torque based on the span position of the blade. Figure 8 plots the contribution of the axial torque along the span of blade 1. According to the flow visualization shown in Figure 7, the three regions of the blade span (i.e., from the root to 50% span, from 50% span to 75% span, and from 75% span to the tip) were considered. In the first phase (Figure 8 (i)), the axial torque generated by the span position from 50% to the blade tip improved because of the flow control using the DBD-PA, and it was also equal to the axial torque generated by the blade span 50% to the blade root. In the second phase (Figure 8 (ii)), the axial torque generated by the span position from the blade root to 50% span improved due to the control of flow by the DBD-PA.
found in the region between 50% of the blade's span and the blade tip. By contrast, in the second phase (Figure 7 (ii)), the flow around blade 1 near the tip was weakly separated. Thus, the suppression of the flow separation due to the flow control was lesser compared with that in the first phase. Thus, in the first phase (Figures 6 (i) and 7 (i)), the increase in the axial torque generation from the flow control using the DBD-PA was more present in comparison with the second phase (Figures 6 (ii) and 7 (ii)). The flow visualization illustrates that the effects of the flow control using DBD-PAs on the flow structure around turbine blades vary with the span direction. Here, we further analyzed the relationship between the flow structures and the contribution of the axial torque based on the span position of the blade. Figure 8 plots the contribution of the axial torque along the span of blade 1. According to the flow visualization shown in Figure 7, the three regions of the blade span (i.e., from the root to 50% span, from 50% span to 75% span, and from 75% span to the tip) were considered. In the first phase (Figure 8 (i)), the axial torque generated by the span position from 50% to the blade tip improved because of the flow control using the DBD-PA, and it was also equal to the axial torque generated by the blade span 50% to the blade root. In the second phase (Figure 8 (ii)), the axial torque generated by the span position from the blade root to 50% span improved due to the control of flow by the DBD-PA. Furthermore, Figure 9 shows a comparison of the pressure distributions over the blade section at span positions of 50% and 75%. Here, the surface pressure coefficient was normalized based on the dynamic pressure of the effective flow velocity, and l in Figure 9 indicates the local chord length from the leading edge of the blade. Figure 9 presents the difference between the cases with and without flow control using the DBD-PA except for the results at the 75% span section in the second phase (Figure 9b (ii) 75% span). The differences in the surface pressure distribution were related to the axial torque generation. The suction peaks were enhanced with the DBD-PA, and the surface pressure distribution at the 50% span section in the first phase (Figure 9a (i) 50% span) indicated a vortex convection near the blade surface. In the second phase, the surface pressure distributions at the 50% and 75% span sections with the DBD-PA indicated the presence of a laminar separation bubble on the suction side of the blade.  Furthermore, Figure 9 shows a comparison of the pressure distributions over the blade section at span positions of 50% and 75%. Here, the surface pressure coefficient was normalized based on the dynamic pressure of the effective flow velocity, and l in Figure 9 indicates the local chord length from the leading edge of the blade. Figure 9 presents the difference between the cases with and without flow control using the DBD-PA except for the results at the 75% span section in the second phase (Figure 9b (ii) 75% span). The differences in the surface pressure distribution were related to the axial torque generation. The suction peaks were enhanced with the DBD-PA, and the surface pressure distribution at the 50% span section in the first phase (Figure 9a (i) 50% span) indicated a vortex convection near the blade surface. In the second phase, the surface pressure distributions at the 50% and 75% span sections with the DBD-PA indicated the presence of a laminar separation bubble on the suction side of the blade. dynamic pressure of the effective flow velocity, and l in Figure 9 indicates the local chord length from the leading edge of the blade. Figure 9 presents the difference between the cases with and without flow control using the DBD-PA except for the results at the 75% span section in the second phase (Figure 9b (ii) 75% span). The differences in the surface pressure distribution were related to the axial torque generation. The suction peaks were enhanced with the DBD-PA, and the surface pressure distribution at the 50% span section in the first phase (Figure 9a (i) 50% span) indicated a vortex convection near the blade surface. In the second phase, the surface pressure distributions at the 50% and 75% span sections with the DBD-PA indicated the presence of a laminar separation bubble on the suction side of the blade. In addition, we examined the turbulent kinetic energy (TKE) distribution by collecting the maximum value of TKE in the direction normal to the surface (see Figure 10). The rapid increase in the TKE and higher maximum TKE (TKEmax) was presented in the cases with flow control using the DBD-PA at the 50% span section. The results implicate that DBD-PAs promote the laminar-toturbulent tension of separated leading edge shear layers, which has a similar influence on the separated leading edge shear layer of the airfoil that was reported in stationary airfoil cases [13][14][15][16][17]. In addition, we examined the turbulent kinetic energy (TKE) distribution by collecting the maximum value of TKE in the direction normal to the surface (see Figure 10). The rapid increase in the TKE and higher maximum TKE (TKE max ) was presented in the cases with flow control using the DBD-PA at the 50% span section. The results implicate that DBD-PAs promote the laminar-to-turbulent tension of separated leading edge shear layers, which has a similar influence on the separated leading edge shear layer of the airfoil that was reported in stationary airfoil cases [13][14][15][16][17]. The results so far are the observations and findings related to the leading edge separation of the wind turbine blade. Here, the wake structures can be seen, especially the tip vortices. Recalling Figure  7, the tip vortex around blade 1 remained longer when using a DBD-PA for flow control. The differences in the tip-vortex length around blade 1 in cases with control and without control were computed based on iso-Q surfaces. The differences were approximately 4.0c at = 112.5 degrees and almost zero at = 247.5 degrees, respectively. In the case without flow control, the separated flow and vortices were convected to the downstream and interacted with the tip vortex. However, such interactions can cause a breakdown of the tip vortex. The results in the cases with flow control show a suppression of the leading edge separation near the tip region, which leads to very weak vortexflow interactions between the tip vortex and the separated flow from the leading edge. This is the secondary effect of the flow control of DBD-PAs, and it may play a role in the reduction of sound generation. Of note, the considered primary effect of DBD-PAs in this study is the suppression or delay of the leading edge separation of wind turbine blades.

Discussion about Flow Control around the Sophisticated Wind Turbine Airfoil Cross-Section
In the previous subsection, the effects of the DBD-PA-based flow control on wind turbine blades were presented and discussed. In comparison with those reported in a previous study [26], the The results so far are the observations and findings related to the leading edge separation of the wind turbine blade. Here, the wake structures can be seen, especially the tip vortices. Recalling Figure 7, the tip vortex around blade 1 remained longer when using a DBD-PA for flow control. The differences in the tip-vortex length around blade 1 in cases with control and without control were computed based on iso-Q surfaces. The differences were approximately 4.0c at φ = 112.5 degrees and almost zero at φ = 247.5 degrees, respectively. In the case without flow control, the separated flow and vortices were convected to the downstream and interacted with the tip vortex. However, such interactions can cause a breakdown of the tip vortex. The results in the cases with flow control show a suppression of the leading edge separation near the tip region, which leads to very weak vortex-flow interactions between the tip vortex and the separated flow from the leading edge. This is the secondary effect of the flow control of DBD-PAs, and it may play a role in the reduction of sound generation. Of note, the considered primary effect of DBD-PAs in this study is the suppression or delay of the leading edge separation of wind turbine blades.

Discussion about Flow Control around the Sophisticated Wind Turbine Airfoil Cross-Section
In the previous subsection, the effects of the DBD-PA-based flow control on wind turbine blades were presented and discussed. In comparison with those reported in a previous study [26], the discrepancy observed can be due to three possible reasons. First, the difference in the airfoil shape affected the separation characteristics and hysteresis of the flow around the blade. In particular, for the dynamic stalled case, the previous study presented a considerable difference between the NACA0012 and NACA63 6 -618 airfoils in terms of the separation point in time [50]. In addition, in the case of the NACA0012 airfoil, the separation point was fixed at the leading edge while in the case of the NACA63 6 -618 airfoil the separation point changed in time [50]. Moreover, the difference in the airfoil geometry can result in changes in the behavior of flow separation, reattachment, and secondary separation. In a well-designed airfoil, such as the NRELS822, the trailing-edge separation was more prominent compared with that in the NACA0012 airfoil. In such flow fields, a DBD-PA installed at the leading edge cannot directly affect the flow near the trailing edge, so the aerodynamic performances of the airfoil cannot be effectively enhanced. Second, the yaw angle of the freestream causes a time variation of the effective angle of attack of the blade; hence, the flow becomes more unsteady in time and three-dimensional in space. As the yaw angles were not considered in the previous study [26], the flow dynamics around the blades were most likely different, and less leading edge separation was expected. Third, the effectiveness of the operational condition of the DBD-PA depends on the flow structures around the blades. Therefore, it is necessary to further study and prove the above-mentioned three points in the near future.

Conclusions
The effects of DBD-PAs on axial toque generation and flow structures around simple experimental horizontal-axis wind turbine blades under a low tip-ratio condition were numerically investigated and reported in this paper. A DBD-PA was modeled as a spatial body force distribution with time variations, and the flow around the rotating blades was simulated based on LES with the compact finite differencing scheme. Moreover, large-scale parallel computations were conducted using 18,800 cores of the K computer. As the low tip-ratio condition and an inclined incoming flow were considered, the computed flow around the blades was dynamically changed with the time variation of the azimuth angles. Moreover, the effects of the DBD-PA on axial torque generation were evaluated based on three revolutions after the DBD-PA was switched on. The results show that revolution-averaged axial torque approximately increased by 11%-19% compared to the case without a DBD-PA. Here, the DBD-PA seemed to mainly affect the separated leading edge shear layer, leading to a delayed leading edge separation of the wind turbine blades. Thus, DBD-PAs can improve the axial torque generation. Further experimental studies to validate the model for the present test case are required for future research.