Numerical Investigations of the Effects of the Rotating Shaft and Optimization of Urban Vertical Axis Wind Turbines

The central shaft is an important and indispensable part of a small scale urban vertical axis wind turbines (VAWTs). Normally, it is often operated at the same angular velocity as the wind turbine. The shedding vortices released by the rotating shaft have a negative effect on the blades passing the wake of the wind shaft. The objective of this study is to explore the influence of the wake of rotating shaft on the performance of the VAWT under different operational and physical parameters. The results show that when the ratio of the shaft diameter to the wind turbine diameter (α) is 9%, the power loss of the wind turbine in one revolution increases from 0% to 25% relative to that of no-shaft wind turbine (this is a numerical experiment for which the shaft of the VAWT is removed in order to study the interactions between the shaft and blade). When the downstream blades pass through the wake of the shaft, the pressure gradient of the suction side and pressure side is changed, and an adverse effect is also exerted on the lift generation in the blades. In addition, α = 5% is a critical value for the rotating shaft wind turbine (the lift-drag ratio trend of the shaft changes differently). In order to figure out the impacts of four factors; namely, tip speed ratios (TSRs), α, turbulence intensity (TI), and the relative surface roughness value (ks/ds) on the performance of a VAWT system, the Taguchi method is employed in this study. The influence strength order of these factors is featured by TSRs > ks/ds > α > TI. Furthermore, within the range we have analyzed in this study, the optimal power coefficient (Cp) occurred under the condition of TSR = 4, α = 5%, ks/ds = 1 × 10−2, and TI = 8%.


Introduction
Wind turbines can be categorized into horizontal axis wind turbines (HAWTs) and vertical axis wind turbines (VAWTs) [1,2].The HAWTs have been extensively employed in the large scale wind farm far away from urban areas due to the loading noise that is created by the large rotors [3,4].In order to increase the power of the HAWTs, the size of the turbine is getting larger than before and this brings difficulties to the maintenance of the wind farm [5].In contrast to HAWTs, VAWTs have less commercial applications; the power generation capacity of VAWTs is less than that of HAWTs [6][7][8].However, VAWTs also have some noteworthy advantages, for instance, the VAWTs can capture unstable and turbulent wind energy from any direction, and they also produce less noise.In addition, the manufacturing and maintenance of VAWTs are much easier [9][10][11][12].Therefore, the VAWTs are more suitable for urban environment relative to traditional HAWTs.In recent years, the amount of relevant studies on small scale urban VAWTs are increasing.Overall, a review of recent studies on VAWTs reveals that the research directions with regard to VAWTs can be divided into four directions: the studies of a single airfoil profile, a single turbine, a pair of wind turbines, and multiple turbines.For the study of airfoil profile, most researches are devoted to propose a new design method to reduce the unsteady aerodynamic phenomenon of dynamic stall in order to improve the lift-drag ratio of the airfoil.For example, Zhong et al. [13] added a cylinder at the front of leading edge of the airfoil in order to delay the dynamic stall angle, and the lift-drag ratio can be also increased by 30.7%.For the research of a single wind turbine, although there is only a single VAWT in the flow field and there are no interactions between the wind turbine and other wind turbines in the process of research, the flow phenomenon in these researches is also complex, for instance, the phenomenon of dynamic stall during the operation of the wind turbine.Most of these researches apply the improved blade profile to the wind turbine, and they investigate power output and the variation of lift and drag coefficient of each blade in the wind turbine.Zamani et al. [14] investigated the variations of power coefficient of a VAWT with new J-shaped blade by three-dimensional (3D) CFD (Computational fluid dynamics) method, and compared the results to the transitional wind turbine with normal symmetrical airfoil.For a pair of wind turbines, in order to find the condition of optimal power output of the dual turbine system, the Taguchi method was employed in the research of Chen et al. [15].They found that the power coefficient of the dual turbine system will be enhanced by 9.97% when compared to the single wind turbine when the wind turbine operated under the optimal condition, as summarized by the Taguchi method, which is an analysis method of combination of experiment and data statistic developed by Genichi Taguchi [16].For multiple turbines, in order to increase the production of energy in existing wind farms, Xie et al. [17] conducted a 3D numerical simulation, where 20 small VAWTs were added around each large HAWT in a traditional wind farm.They found that the modified vertically staggered wind farm produced up to 32% increment of power as compared to a traditional wind farm.
The attention of the effect of the tower shadow on the performance of the large wind turbine was paid in the previous studies [18][19][20].Besides, lots of previous studies on large scale VAWTs, the ratio of the shaft diameter to turbine diameter is very small, so that the effect of the shaft can be negligible.However, for the small scale urban VAWTs, the ratio of the shaft diameter to the turbine diameter is relative large.There will be significant effect acting on the downstream blades, due to the blades passing through the wake of the shaft.Rezaeiha et al. [21] find that the air flow passing around the central shaft of the VAWT produces shedding vortices, which exert significant impact on the downstream blades.In addition, the influences of the wake of the shaft on the downstream blades are also studied by Chigliaro [22].According to research results of Rezaeiha [21] and Chigliaro [22], many factors, including the tip speed ratios (TSRs) of the operational wind turbine, the turbulence intensity (TI) of the flow, and the ratio of the shaft diameter to the wind turbine diameter (α) all can affect the wake of the central rotating shaft of the VAWT.Besides the relative surface roughness value (k s /d s ) of the shaft is also an important influencing factor to the power output of a wind turbine.
In this paper, the research objectives of this current study mainly include the following points: i. Analyze the reasons for power coefficient of wind turbines decrease as the ratio of the shaft diameter to the wind turbine diameter (α) increase from 0% to 9% quantitatively.ii.
Investigate the vortex variations behind the rotating shaft of the VAWT and the variations of the boundary separation point as α increase from 0% to 9%.iii.
Explore the wind velocity distribution ruler of the flow field when the vertical axis wind turbine operates under different operating conditions (including TI and TSRs) and physical parameters (including α and k s /d s ).iv.
Many factors can have impacts on the wake effect of the rotating central shaft on the performance of urban VAWTs, for instance, α, TSRs, k s /d s , and TI.One of the purposes of this research is to figure out the influence strength order of these factors, and to find out an optimal operational condition for the VAWT.
The outline of this present study obeys the following order: detail description of wind turbine model, grid generation method, boundary conditions, and verifications for both central cylinder and wind turbine are presented in Section 2. In Section 3.1, the variations of power coefficient are investigated in detail when α increases from 0% to 9%.The flow characteristics around the rotating shaft and the variations of the boundary separation point are discussed in Sections 3.2-3.4present a quantitative analysis to investigate the wind velocity distribution at the different downstream positions behind the rotating shaft and the pressure variations on the both sides of the blade surface.In Section 4, the performance analysis of the VAWT is presented with the Taguchi method.Besides, the influence strength order of α, TSRs, k s /d s , and TI is also found out in Section 4.1, and an optimal operational condition for the VAWT is suggested in Section 4.2.

Model Geometry
To investigate the impact of the rotating shaft on the aerodynamic performance, a three-bladed H-type VAWT equipped with symmetric NACA0022 airfoils is used; the chord length (C) of the blade is 0.04 m.The VAWT has a diameter (D) of 0.7 m and the diameter of the wind turbine's shaft (d s ) is 27 mm, meanwhile, d s is the characteristic value of the diameter-based Reynolds number (Re s ).The blade length (H) of the VAWT is 0.6 m.The detail features of the VAWT are summarized in Table 1.The experimental power coefficient of this VAWT can be found in the previous study of Danao [23].In order to investigate the influence of the ratio of the shaft diameter to the wind turbine diameter (α) on the aerodynamic performance of VAWT, five different values are chosen for α, as follows: α = 3.9%, α = 5%, α = 6%, α = 7%, and α = 9%, respectively.As α increases from 3.9% to 9%, the Re s increases from 1.31 × 10 4 to 2.9 × 10 4 .The initial positions and serial numbers of the three blades are shown in Figure 1.θ represents the azimuth angle of the blade.The position of blade 1 is defined as an initial position of one cycle.In this study, this initial position is corresponding to θ = 0 • .Figure 1.Schematic of different mesh zones in two-dimensional (2D) computational domain.

Computational Domain and Mesh Generation
The two dimensional computational domain is employed for the simulation, as shown in Figure 1, the width of the computational domain is 20D (D is the rotating diameter of the VAWT).In order to minimize the effects of blockage on the simulation results, the rotor axis distance from both the inlet and outlet side are 10D [24].The whole computational domain consists of a rotating region where the wind turbine is located and a stationary outside region surrounding the rotating region.To simulate the rotation of the VAWT, the two regions are connected by an interface boundary condition and the inside rotating region enables the rotation of the VAWT.In order to ensure the continuity of air flow during the simulation process, the diameter of the rotating region is set to 1.5D [24,25].
Grid generation is a critical part of the numerical simulation.Good grid quality can improve the accuracy of wind turbine simulation.Accordingly, a structured O-grid of quadrilateral elements is generated around the three blades of the VAWT to ensure the quality of the mesh around the blade and a structured grid is chosen for the inner rotating region and the outer stationary domain in all cases.As shown in Figure 2, in order to clearly observe the influence of existence of the shaft on the performance of the VAWT obviously, we generate another grid without a shaft for comparison calculation.In addition, the accuracy of the simulation with Transition SST turbulence model in this study are related to the dimensionless wall distance y + of blade surface and shaft surface [26,27], and the y + can be defined as y + = ρuτy/μ, where ρ is density, uτ is the friction velocity, y is the boundary layer length, and μ is dynamic viscosity.In order to adapt to the requirements of the turbulence model and to ensure the accuracy of simulation, the boundary layer grid is employed on the shaft surface and the surface of the airfoil, while the heights of the first layer grid of the shaft and airfoil are 10 −4 D and 10 −5 C, respectively, in order to ensure that y + < 1 [13,26].

Computational Domain and Mesh Generation
The two dimensional computational domain is employed for the simulation, as shown in Figure 1, the width of the computational domain is 20D (D is the rotating diameter of the VAWT).In order to minimize the effects of blockage on the simulation results, the rotor axis distance from both the inlet and outlet side are 10D [24].The whole computational domain consists of a rotating region where the wind turbine is located and a stationary outside region surrounding the rotating region.To simulate the rotation of the VAWT, the two regions are connected by an interface boundary condition and the inside rotating region enables the rotation of the VAWT.In order to ensure the continuity of air flow during the simulation process, the diameter of the rotating region is set to 1.5D [24,25].
Grid generation is a critical part of the numerical simulation.Good grid quality can improve the accuracy of wind turbine simulation.Accordingly, a structured O-grid of quadrilateral elements is generated around the three blades of the VAWT to ensure the quality of the mesh around the blade and a structured grid is chosen for the inner rotating region and the outer stationary domain in all cases.As shown in Figure 2, in order to clearly observe the influence of existence of the shaft on the performance of the VAWT obviously, we generate another grid without a shaft for comparison calculation.In addition, the accuracy of the simulation with Transition SST turbulence model in this study are related to the dimensionless wall distance y + of blade surface and shaft surface [26,27], and the y + can be defined as y + = ρu τ y/µ, where ρ is density, u τ is the friction velocity, y is the boundary layer length, and µ is dynamic viscosity.In order to adapt to the requirements of the turbulence model and to ensure the accuracy of simulation, the boundary layer grid is employed on the shaft surface and the surface of the airfoil, while the heights of the first layer grid of the shaft and airfoil are 10 −4 D and 10 −5 C, respectively, in order to ensure that y + < 1 [13,26].

Figure 2.
Meshing details: entire mesh in computational domain, mesh in the inner region and mesh near the airfoil.

Boundary Conditions and Numerical Settings
To ensure the accuracy of the simulation results, the setting of boundary conditions and the methods used in the solution are based on a previous experiment that was reported by Danao [23].In this study, the left side of the computational domain is the velocity-inlet boundary (V = 7 m/s) with turbulent intensity of 8% and the turbulence viscosity ratio is set to 14, the turbulence parameters in the velocity-inlet boundary are set according to the experiment measure reported by Danao [23].The right side of the domain is the pressure-outlet and the gauge pressure is zero.The upper and lower edges of the stationary domain are defined as the symmetry boundary conditions.No-slip boundary condition is specified at surfaces of the shaft and airfoil.
In ANSYS Fluent, The SIMPLE algorithm is employed for the coupling of pressure and velocity.Second-order implicit transient formulation is used and second-order upwind spatial discretization is employed in the pressure, momentum, and turbulence equations.The under relaxation factors for the turbulent kinetic energy, specific dissipation rate, intermittency, and turbulent viscosity are set to 0.8, 0.8, 0.8, 1.During the solution process, the convergence criteria are set to 10 −5 .In the simulation conducted as part of this study, we considered that the wind turbine performed a rotation of 30 cycles, in all cases.When the monitoring parameters changed periodically, we investigated the last revolution.
According to previous research, for the wind turbine model in this study, the results of the simulation using the four-equation transition SST turbulence model are more accurate in comparison to using the k-ω SST turbulence model [28,29].The transition SST turbulence model is a very effective tool for simulating the various transition processes that were proposed by Menter [30].The transition SST model consists of two additional transport equations and two k-ω SST transport equations.Therefore, we also use the four-equation transition SST turbulence model in order to simulate all the cases in this study [31,32].

Boundary Conditions and Numerical Settings
To ensure the accuracy of the simulation results, the setting of boundary conditions and the methods used in the solution are based on a previous experiment that was reported by Danao [23].In this study, the left side of the computational domain is the velocity-inlet boundary (V = 7 m/s) with turbulent intensity of 8% and the turbulence viscosity ratio is set to 14, the turbulence parameters in the velocity-inlet boundary are set according to the experiment measure reported by Danao [23].The right side of the domain is the pressure-outlet and the gauge pressure is zero.The upper and lower edges of the stationary domain are defined as the symmetry boundary conditions.No-slip boundary condition is specified at surfaces of the shaft and airfoil.
In ANSYS Fluent, The SIMPLE algorithm is employed for the coupling of pressure and velocity.Second-order implicit transient formulation is used and second-order upwind spatial discretization is employed in the pressure, momentum, and turbulence equations.The under relaxation factors for the turbulent kinetic energy, specific dissipation rate, intermittency, and turbulent viscosity are set to 0.8, 0.8, 0.8, 1.During the solution process, the convergence criteria are set to 10 −5 .In the simulation conducted as part of this study, we considered that the wind turbine performed a rotation of 30 cycles, in all cases.When the monitoring parameters changed periodically, we investigated the last revolution.
According to previous research, for the wind turbine model in this study, the results of the simulation using the four-equation transition SST turbulence model are more accurate in comparison to using the k-ω SST turbulence model [28,29].The transition SST turbulence model is a very effective tool for simulating the various transition processes that were proposed by Menter [30].The transition SST model consists of two additional transport equations and two k-ω SST transport equations.Therefore, we also use the four-equation transition SST turbulence model in order to simulate all the cases in this study [31,32].

Validation Study of Central Cylinder
A two-dimensional (2D) computational domain is employed to the validation study of central cylinder, as shown in Figure 3a, the width of the computational domain is 20d (d is the cylinder diameter, d = 0.027 m) and the cylinder distances from both the inlet and outlet side are 10d and 15d, respectively.A structured O-grid of quadrilateral elements is generated around the cylinder.Moreover, the first layer grid of the cylinder is 10 −5 d in order to ensure that y + < 1.The overall grid diagram and the cylinder grid enlargement diagram are shown in Figure 3b.The parameter settings in the CFD simulation are the same as that in Section 2.3.

Validation Study of Central Cylinder
A two-dimensional (2D) computational domain is employed to the validation study of central cylinder, as shown in Figure 3a, the width of the computational domain is 20d (d is the cylinder diameter, d = 0.027 m) and the cylinder distances from both the inlet and outlet side are 10d and 15d, respectively.A structured O-grid of quadrilateral elements is generated around the cylinder.Moreover, the first layer grid of the cylinder is 10 −5 d in order to ensure that y + < 1.The overall grid diagram and the cylinder grid enlargement diagram are shown in Figure 3b.The parameter settings in the CFD simulation are the same as that in Section 2.3.To ensure the simulation results, the effectiveness of the mesh generation method and the rationality of parameter setting in solving process are reliable, the flow field around a twodimensional smooth cylinder with a Reynolds number (Re) of 1.32 × 10 4 is simulated.The simulation results of cylinder are compared against the previous published experimental results by Okamoto and Yagita [33] and Norberg [34] and the simulation results by Prsic [35].Besides, the Strouhal number (St) is chosen as the measuring standard to test the accuracy of simulation.The Strouhal number is expressed as the following Equation, where d is the diameter of the cylinder, f is the frequency of shedding vortex, and U∞ is free wind velocity.In this study, the frequency of the shedding vortex can be obtained from the spectral analysis of the lift force fluctuation; the f in this study is 55.79 Hz.We can calculate that St is equal to 0.2152.
The comparison results are shown in Table 2.By observing Table 2, we can find that the present simulation result shows a good agreement with the previous experiment and simulation results.The St calculated from the present simulation is only 2.5% and 7.6% higher than that of the experiment result of Okamoto [33] and Norberg [34].It can be found by comparing with the previous research results that the simulation result in the present study is accurate and reliable.To ensure the simulation results, the effectiveness of the mesh generation method and the rationality of parameter setting in solving process are reliable, the flow field around a two-dimensional smooth cylinder with a Reynolds number (Re) of 1.32 × 10 4 is simulated.The simulation results of cylinder are compared against the previous published experimental results by Okamoto and Yagita [33] and Norberg [34] and the simulation results by Prsic [35].Besides, the Strouhal number (St) is chosen as the measuring standard to test the accuracy of simulation.The Strouhal number is expressed as the following Equation, where d is the diameter of the cylinder, f is the frequency of shedding vortex, and U ∞ is free wind velocity.In this study, the frequency of the shedding vortex can be obtained from the spectral analysis of the lift force fluctuation; the f in this study is 55.79 Hz.We can calculate that St is equal to 0.2152.The comparison results are shown in Table 2.By observing Table 2, we can find that the present simulation result shows a good agreement with the previous experiment and simulation results.The St calculated from the present simulation is only 2.5% and 7.6% higher than that of the experiment result of Okamoto [33] and Norberg [34].It can be found by comparing with the previous research results that the simulation result in the present study is accurate and reliable.In this study, the numerical model of equivalent sand is used to simulate the different central rotating shaft's surface roughness.This model simulates that there are some sand particles with different diameters on the cylinder surface.As shown in Figure 4, k s represents the diameter of the sand model and k s /d s is the relative surface roughness value, where d s is the smooth shaft diameter of the VAWT [36,37].In this study, the numerical model of equivalent sand is used to simulate the different central rotating shaft's surface roughness.This model simulates that there are some sand particles with different diameters on the cylinder surface.As shown in Figure 4, ks represents the diameter of the sand model and ks/ds is the relative surface roughness value, where ds is the smooth shaft diameter of the VAWT [36,37].

Grid Independence and Numerical Validation
In order to ensure the accuracy of simulation in this study, the boundary conditions and the solution parameters are all setting according to the previous experiment reported by Danao [23] with a wind velocity of 7 m/s.Appropriate CFD numerical simulation validation investigation has been carried out in this section.
First, in order to complete the verification of grid independence, three different mesh systems, including coarse, medium, and fine meshes have been considered in this section.The detailed descriptions of the characteristic sizes of the three different quality grids are presented in Table 3.We mainly compare the power coefficient to the experimental power coefficient of the wind turbine for the three grid resolutions at TSR = 4.The tip speed ratio (TSR) is defined as where ω is the rotor angular velocity of the VAWT, R is the rotating radius of VAWT, and U∞ is the free wind velocity.According to the results, when the grid quality is coarse, the power coefficient of the simulation prediction is not accurate.The simulation results are similar for the medium and fine grids; therefore,

Grid Independence and Numerical Validation
In order to ensure the accuracy of simulation in this study, the boundary conditions and the solution parameters are all setting according to the previous experiment reported by Danao [23] with a wind velocity of 7 m/s.Appropriate CFD numerical simulation validation investigation has been carried out in this section.
First, in order to complete the verification of grid independence, three different mesh systems, including coarse, medium, and fine meshes have been considered in this section.The detailed descriptions of the characteristic sizes of the three different quality grids are presented in Table 3.We mainly compare the power coefficient to the experimental power coefficient of the wind turbine for the three grid resolutions at TSR = 4.The tip speed ratio (TSR) is defined as where ω is the rotor angular velocity of the VAWT, R is the rotating radius of VAWT, and U ∞ is the free wind velocity.According to the results, when the grid quality is coarse, the power coefficient of the simulation prediction is not accurate.The simulation results are similar for the medium and fine grids; therefore, for the sake of independency of the solution to the meshes, the medium mesh can be chosen as the most suitable grid in the further simulations.
To investigate the effect of the time step on the simulation results of the VAWT, three different time steps are considered, as listed in Table 4.Then, simulation validation is performed using these three different time steps.Thirty revolutions are simulated for the VAWT at TSR = 4 and the last revolution is chosen to analyze the instantaneous torque coefficient (C m ) of blade 1 of the wind turbine.The instantaneous torque coefficient C m can be defined as [38] where T is the torque of blade, A is the project area of VAWT (A = RH, H is the blade length), ρ is the air density, R is the rotor radius of VAWT, and U ∞ is the free wind velocity.As can be seen in Figure 5, when the time step size is 0.109018 milliseconds (corresponds to the time of wind turbine rotates 0.5 • at a time), the maximum value of C m is 0.0022 higher than that of the simulation results that were obtained from the time step size of 0.21817 milliseconds (corresponds to the time of wind turbine rotates 1 • at a time), and it is only 0.0011 lower than that of the simulation results that were obtained from the time step size of 0.054541 milliseconds (corresponds to the time of wind turbine rotates 0.25 • at a time).This also results the power coefficient C p of the calculated result of time step of 0.109018 milliseconds is 0.7% and 0.21% higher than that of the calculated result of time step of 0.21817 milliseconds and 0.054541 milliseconds.The power coefficient C p can be defined as Energies 2018, 11, x 8 of 24 for the sake of independency of the solution to the meshes, the medium mesh can be chosen as the most suitable grid in the further simulations.
To investigate the effect of the time step on the simulation results of the VAWT, three different time steps are considered, as listed in Table 4.Then, simulation validation is performed using these three different time steps.Thirty revolutions are simulated for the VAWT at TSR = 4 and the last revolution is chosen to analyze the instantaneous torque coefficient (Cm) of blade 1 of the wind turbine.The instantaneous torque coefficient Cm can be defined as [38] 2 0.5 where T is the torque of blade, A is the project area of VAWT (A = RH, H is the blade length), ρ is the air density, R is the rotor radius of VAWT, and U∞ is the free wind velocity.As can be seen in Figure 5, when the time step size is 0.109018 milliseconds (corresponds to the time of wind turbine rotates 0.5° at a time), the maximum value of Cm is 0.0022 higher than that of the simulation results that were obtained from the time step size of 0.21817 milliseconds (corresponds to the time of wind turbine rotates 1° at a time), and it is only 0.0011 lower than that of the simulation results that were obtained from the time step size of 0.054541 milliseconds (corresponds to the time of wind turbine rotates 0.25° at a time).This also results the power coefficient Cp of the calculated result of time step of 0.109018 milliseconds is 0.7% and 0.21% higher than that of the calculated result of time step of 0.21817 milliseconds and 0.054541 milliseconds.The power coefficient Cp can be defined as  The calculation methods of TSR and C m can be found in Equations ( 2) and (3).Due to the C p calculated by the time step size of 0.109018 milliseconds is only 0.7% and 0.21% higher than that of time step size of 0.21817 and 0.05451 milliseconds, the differences between the three different time step sizes are small.This simulation results are in good agreement with the experimental and simulated results of Danao [23], as shown in Figure 6 (in order to change the TSRs, the simulations and the experiments all keep the wind speed constant and change the angular velocity).
To ensure that the simulation data were adequate, and to reduce the time that is required for the calculation, we chose the time step of 0.109018 milliseconds with respect to an azimuthal angle of 0.5 • .The results of investigating the time step are in good agreement with the results of Rezaeiha [24,39].It can be observed from Figure 6 that when the TSR < 4, the simulation results are in good agreement with the results of experiment both in overall trend and the value of C p .However, when TSR = 4, the deviation of simulation prediction results show an increasing trend, the prediction C p of simulation is 33.2% higher than that of the experiment.There is an obvious difference between the results of experiment and the simulation at the range of 4 < TSR < 5.For the aspects of optimal value prediction, the maximum C p for the experimental results is 0.21 occurred at TSR = 4, however, the maximum C p of the present simulation is 0.321 when TSR = 4.5, and this brings a deviation of 52.8% for the prediction of maximum C p .
The calculation methods of TSR and Cm can be found in Equations ( 2) and (3).Due to the Cp calculated by the time step size of 0.109018 milliseconds is only 0.7% and 0.21% higher than that of time step size of 0.21817 and 0.05451 milliseconds, the differences between the three different time step sizes are small.This simulation results are in good agreement with the experimental and simulated results of Danao [23], as shown in Figure 6 (in order to change the TSRs, the simulations and the experiments all keep the wind speed constant and change the angular velocity).To ensure that the simulation data were adequate, and to reduce the time that is required for the calculation, we chose the time step of 0.109018 milliseconds with respect to an azimuthal angle of 0.5°.The results of investigating the time step are in good agreement with the results of Rezaeiha [24,39].It can be observed from Figure 6 that when the TSR < 4, the simulation results are in good agreement with the results of experiment both in overall trend and the value of Cp.However, when TSR = 4, the deviation of simulation prediction results show an increasing trend, the prediction Cp of simulation is 33.2% higher than that of the experiment.There is an obvious difference between the results of experiment and the simulation at the range of 4 < TSR < 5.For the aspects of optimal value prediction, the maximum Cp for the experimental results is 0.21 occurred at TSR = 4, however, the maximum Cp of the present simulation is 0.321 when TSR = 4.5, and this brings a deviation of 52.8% for the prediction of maximum Cp.The possible reasons for the simulation error when the TSRs are in the range of 1 to 5 are concluded, as follows: (1) In the process of measuring the experimental data, a method of spin down tests is used, it leads the actual wind speed acted on the wind turbine is lower than the set wind speed as the presence of blockage effects.Therefore, the measured results will still be small, although some corrections are made during the subsequent data processing.(2) Since a simplified 2D geometrical model is employed in the simulation and the influence of the blade span and junctions connected the blades and the shaft also not considered in the process of simulating, a deviation might be existed between the experimental results and the simulation results.It is noteworthy that for high TSRs (4 < TSRs < 5), there is a mass of drag can be produced by the junctions connected the blades and the shaft in the actual operation of the VAWT.In addition, the shedding vortices that are generated by the internal junctions during the operation of the VAWT have not been accurately simulated with a 2D simplified model.These reasons have great influence on the prediction results of numerical simulation under high TSRs.(3) In the course of numerical calculation, the accuracy of using simplified 2D unsteady Reynoldsaverage Navier-Stokes (URANS) to solve the 3D problem with complex flow characteristics is limited [40].Especially for the study of blade-wake interactions in this paper.This is also one of The possible reasons for the simulation error when the TSRs are in the range of 1 to 5 are concluded, as follows: (1) In the process of measuring the experimental data, a method of spin down tests is used, it leads the actual wind speed acted on the wind turbine is lower than the set wind speed as the presence of blockage effects.Therefore, the measured results will still be small, although some corrections are made during the subsequent data processing.(2) Since a simplified 2D geometrical model is employed in the simulation and the influence of the blade span and junctions connected the blades and the shaft also not considered in the process of simulating, a deviation might be existed between the experimental results and the simulation results.It is noteworthy that for high TSRs (4 < TSRs < 5), there is a mass of drag can be produced by the junctions connected the blades and the shaft in the actual operation of the VAWT.In addition, the shedding vortices that are generated by the internal junctions during the operation of the VAWT have not been accurately simulated with a 2D simplified model.These reasons have great influence on the prediction results of numerical simulation under high TSRs.
(3) In the course of numerical calculation, the accuracy of using simplified 2D unsteady Reynolds-average Navier-Stokes (URANS) to solve the 3D problem with complex flow characteristics is limited [40].Especially for the study of blade-wake interactions in this paper.This is also one of the possible reasons for the differences between the simulated predicted results and experimental measurements.

Taguchi Method
The Taguchi method is widely used in engineering practice, and it is a combination of a statistical approach and experimental design [16,41].The Taguchi method can not only analyze the optimization scheme for a product or a production process, but it can also find out the influence strength order of the influent factors in the process of the analysis [42].
When compared to the traditional statistical method, the Taguchi method has the following advantages: the Taguchi method emphasizes the employment of loss functions, it can also help researchers to understand the relationships between each response and the factors preferably.In addition, due to the Taguchi method assuming that there are no interactions among the control factors, this leads researcher can investigate more influences of control factors during a small number of experiments.Simultaneously, since a data statistical method of orthogonal array is employed to the factors, the experiment cost can be decreased.The disadvantages of Taguchi method are as follows: only the optimal scheme can be found by Taguchi method from the specified parameter level combinations, so the feasible solution space will be constrained once the parameter levels are determined.Besides, the Taguchi method cannot find the optimal solution when the variable of process parameters is continuous [43].
The following steps should be followed when using the Taguchi method for researches [44]: (g) Select the optimal levels of each control factor.
To apply the Taguchi method in this study, a simulation scheme of the L 9 (3 4 ) orthogonal array is created firstly (including four factors, four levels are considered in each factor).Calculate the C p and the signal-noise ratio (S/N ratio) of the VAWT system of each combination in the orthogonal array scheme.Typically, the employment of the S/N ratio is to evaluate the quality characteristics deviating from the desired value.The mean S/N ratios of the four factors are further calculated to investigate the influence strength order for the four factors.Furthermore, the profiles of the mean S/N ratio will be further analyzed to find out the optimal combination within the range of values that we have studied.
There are many possible factors that can have impact on the wake effect of the rotating shaft.According to previous studies, four alternative factors are investigated and analyzed, namely, tip speed ratios (TSRs), shaft diameter to wind turbine diameter ratios (α), the turbulent intensity of incoming flow (TI), and the relative surface roughness of value of the shaft.[21,45].
The Taguchi method is widely used for optimizing industrial processes, it is also a reasonable tool to investigate and optimize the power output of the VAWT's system [15,46].The method is based on the orthogonal array of the control factors, and measure the loss of the product quality quantify through the concept of the loss function, in addition, the signal-noise ratio (S/N ratio) is employed to analyze the quality characteristics of the product or process parameters [47,48].The loss function can be defined as the following quadratic form: where L is the loss value, K is a constant which depends on the magnitude of the characteristic, and T is the target response.There are three different kinds of S/N ratio, nominal-the better (NB), the larger-the-better (LB), and the small-the-better (SB), are defined in the Taguchi method.For the present study, the larger S/N ratio will correspond to a better wind turbine performance, therefore the LB S/N ratio is calculated based on the following function [15,49], In the present study, we use power coefficient to measure the performance of VAWT, the higher the power coefficient means the better the single VAWT system.Accordingly, in Equation ( 6), T is the highest value of a single VAWT and its value is 0.593 [50], while y is the predicted value of the power coefficient of simulation.

Influence of Ratio of the Shaft Diameter to the Wind Turbine Diameter on Wind Turbine Performance
To analyze the influence of the wake effect that is caused by the central rotating shaft of the VAWT on the power coefficient of the blade, we compare the instantaneous power coefficient of blade 1 of different α from 3.9% to 9% in one rotation, to a wind turbine under a hypothetical no-shaft condition, as shown in Figure 7.It can be seen that the power coefficient of the blade fluctuates fiercely when the azimuthal angle is between 250 • and 290 • .The instantaneous power coefficient of the blades obviously declines in the vicinity of an azimuthal angle of 270 • , due to the effect of the central rotating shaft wake vortex.Then, the instantaneous power coefficient of the blades returns to normal after they pass through the region that influenced by the shaft wake.However, the azimuthal angle of lowest point of instantaneous power coefficient of blade 1 goes slightly backwards when α > 7%.In order to further study the reason why the instantaneous power of blade 1 declines obviously in the vicinity of an azimuthal angle of 270 • , the instantaneous contours of velocity is observed in Figure 8 to investigate the state of the wake at this moment.However, it is important to emphasize that the state of the wake is dependent on the frequency of the shedding vortices.When the VAWT operates to the next revolution, the frequency of the shedding vortices is substantially different from that of the current revolution.Therefore, the wake of the rotating shaft in this revolution is not same as the next revolution.As can be seen in Figure 8, the existence of the central rotating shaft increases the low speed area near the blade, while the large velocity change on the pressure side of the blade results in the change of the blade's power coefficient.As the diameter of the wind turbine's central rotating shaft increases, it can be seen that a large-scale vortex struck the surface of the blade.After α ≥ 7%, a large scale vortex appears due to the violent wake effect that is exerted by the central rotating shaft.Consequently, the large-scale vortex does not hit the downstream blades at an azimuthal angle of 270 • .Therefore, at that moment, the lowest point of the instantaneous power coefficient curve of blade 1 is pushed backwards when α ≥ 7%.When the diameter of the central rotating shaft increases again, the scale of the vortex shedding is observed to be larger, and it results in the significant backward shift of the power coefficient curve's minimum point when α = 9%.The change of power coefficients of these VAWTs with the increase of α is shown in Figure 9.It can be seen that the wind turbine's power coefficient exhibits a downward trend as α increases.When α increases to 3.9%, the power coefficient of the wind turbine decreases intensively.As α increase from 0% to 9%, the power loss of the wind turbine increased from 0% to 25%.

Analysis of the Flow Characteristics around the Rotating Shaft
As shown in Figure 7, for α < 7%, the minimum power coefficient occurs when the blade operates at the vicinity of azimuthal angle of 270 • .However, in the case of α ≥ 7%, the minimum value of the power coefficient moves backwards, it appears at the vicinity of azimuthal angle of 280 • .A partially enlarged view of the time-averaged streamlines near shaft (for the last 10 operating revolutions of VAWT) and the instantaneous vortices shedding process of the shaft are shown in Figure 10.As shown in Figure 10, for α = 3.9%, the top and bottom surface flows of the central rotating shaft are symmetrical and without any clear boundary layer separation.Furthermore, the top and bottom surfaces of the central rotating shaft alternate the periodic release of shedding vortices.However, the spacing between the front and back vortices that are produced by the shaft is relatively small, and the vortex intensity is weak, otherwise, the average lift coefficient is approximately zero.When α starts to increase to 5%, the flow become asymmetrical, the separation of the boundary layer on the top surface of the central rotating shaft is apparent, and the wake vortex intensity also become stronger.The range of the effect on the downstream wind turbine blades increases, and the average lift coefficient becomes negative.As α increase, the separation of the boundary layer becomes more obvious, while the separation point keeps moving forward.For α = 9%, the angle between the two separation points is almost close to 180 • and this led to the gradual reduction of the average lift coefficient.The shedding vortex intensity of the central rotating shaft increases gradually, and it exerts a wider effect on the downstream blades of the wind.

Analysis of the Flow Characteristics around the Rotating Shaft
As shown in Figure 7, for α < 7%, the minimum power coefficient occurs when the blade operates at the vicinity of azimuthal angle of 270°.However, in the case of α ≥ 7%, the minimum value of the power coefficient moves backwards, it appears at the vicinity of azimuthal angle of 280°.A partially enlarged view of the time-averaged streamlines near shaft (for the last 10 operating revolutions of VAWT) and the instantaneous vortices shedding process of the shaft are shown in Figure 10.As shown in Figure 10, for α = 3.9%, the top and bottom surface flows of the central rotating shaft are symmetrical and without any clear boundary layer separation.Furthermore, the top and bottom surfaces of the central rotating shaft alternate the periodic release of shedding vortices.However, the spacing between the front and back vortices that are produced by the shaft is relatively small, and the vortex intensity is weak, otherwise, the average lift coefficient is approximately zero.When α starts to increase to 5%, the flow become asymmetrical, the separation of the boundary layer on the top surface of the central rotating shaft is apparent, and the wake vortex intensity also become stronger.The range of the effect on the downstream wind turbine blades increases, and the average lift coefficient becomes negative.As α increase, the separation of the boundary layer becomes more obvious, while the separation point keeps moving forward.For α = 9%, the angle between the two separation points is almost close to 180° and this led to the gradual reduction of the average lift coefficient.The shedding vortex intensity of the central rotating shaft increases gradually, and it exerts a wider effect on the downstream blades of the wind.The lift and drag coefficients of the shaft of VAWT are defined as [51], where the Fl and Fd are the lift and drag forces, respectively, and ds is the diameter of the shaft.The amplitudes in the fluctuation of the shaft are defined, respectively, as [52] ,max ,min ' By examining the vortices' contours of the central rotating shaft shown in Figure 10, it can be observed that the sequential shedding vortices that are released alternately from the top and bottom surfaces of the central rotating shaft at α = 3.9%, and the vortex intensity, are both weak.The intensity of the shedding vortices releases from the wind turbine's central rotating shaft became stronger from the point of α = 5%.Several main vortices are formed and shed backwards one by one.The distance between the two main anteroposterior vortices is large.Moreover, as shown in Figure 11, during the increase of α from 3.9% to 5%, the amplitude of the lift coefficient changes significantly, and it then remains almost unchanged.The amplitude of the drag coefficient is also obvious during the process of α increasing from 3.9% to 5%.The amplitude curve of the drag coefficient increases slowly.Therefore, α = 5% is assessed as a critical state.The changes in the average lift and drag coefficient ratio are shown in Figure 12.Additionally, it can be confirmed that the ratio of the average lift coefficient to the average drag coefficient reaches the maximum value at the critical value.When α > 5%, the large scale shedding vortices exerts a wider effect on the downstream blades.Moreover, the shedding vortex intensity is large, which exerted a significant influence on the instantaneous power coefficient of the wind turbine's downstream blades.The lift and drag coefficients of the shaft of VAWT are defined as [51], where the F l and F d are the lift and drag forces, respectively, and d s is the diameter of the shaft.
The amplitudes in the fluctuation of the shaft are defined, respectively, as By examining the vortices' contours of the central rotating shaft shown in Figure 10, it can be observed that the sequential shedding vortices that are released alternately from the top and bottom surfaces of the central rotating shaft at α = 3.9%, and the vortex intensity, are both weak.The intensity of the shedding vortices releases from the wind turbine's central rotating shaft became stronger from the point of α = 5%.Several main vortices are formed and shed backwards one by one.The distance between the two main anteroposterior vortices is large.Moreover, as shown in Figure 11, during the increase of α from 3.9% to 5%, the amplitude of the lift coefficient changes significantly, and it then remains almost unchanged.The amplitude of the drag coefficient is also obvious during the process of α increasing from 3.9% to 5%.The amplitude curve of the drag coefficient increases slowly.Therefore, α = 5% is assessed as a critical state.The changes in the average lift and drag coefficient ratio are shown in Figure 12.Additionally, it can be confirmed that the ratio of the average lift coefficient to the average drag coefficient reaches the maximum value at the critical value.When α > 5%, the large scale shedding vortices exerts a wider effect on the downstream blades.Moreover, the shedding vortex intensity is large, which exerted a significant influence on the instantaneous power coefficient of the wind turbine's downstream blades.

The Influence of Rotating Central Shaft on the Flow Field inside the Wind Turbine
In order to analyze the variations of velocity inside the wind turbine under different α, the velocity profile at four different positions has been analyzed, as shown in Figure 13. Figure 13a provides a diagram of the four positions located at x/R = 0.2, 0.4, 0.6, and 0.8.In Figure 13b-e, the variations of the velocity inside the wind turbine are investigated in detail.Besides, the velocities and lateral position are non-dimensionalised by the inlet velocity, U∞, and the rotor radius R, respectively.The influence of different α is investigated in Figure 13b when TSR = 4, TI = 8%, and the shaft is smooth at this time, α is increased from 3.9% to 9%.The influence of surface roughness is discussed in Figure 13c when α = 5%, TSR = 4 and TI = 8%, the ks/ds changes from 10 −4 to 10 −2 .Under the condition of α = 5%, TSR = 4, and the shaft is smooth, we study the effects of TI by changing values of TI from 4% to 12%, the comparison results are shown in Figure 13d.In Figure 13e, the influence of TSRs is shown when α = 5%, TI = 8%, the shaft of wind turbine is smooth; the values of TSRs are 2, 3, and 4, respectively.For a single VAWT system, the factors affecting the power coefficient of the VAWT include the TSRs of the VAWT, shaft diameter to wind turbine diameter ratios (α), the surface roughness of the central shaft, and the turbulent intensity of incoming flow (TI).Therefore, four factors and three levels are considered in the simulations to account for their impacts on the performance of the urban VAWT.Furthermore, the optimum operating conditions of the VAWT for maximizing the performance are obtained by Taguchi method in Section 4. It can be observed from Figure 13b that when α = 3.9% and α = 5%, the impact of rotating central shaft on flow field inside the wind turbine is almost the same.From Figure 13c-e, we can see that the surface roughness, turbulence intensity of incoming flow, and TSRs all have an effect on the downstream flow field of the VAWT.The TSRs have a great impact, the TI has a weak influence, and the surface roughness of the central shaft only influences the wake of the shaft and affects the downstream velocity field.The specific effects of these factors on power output of the VAWT will also be discussed in Section 4.

The Influence of Rotating Central Shaft on the Flow Field inside the Wind Turbine
In order to analyze the variations of velocity inside the wind turbine under different α, the velocity profile at four different positions has been analyzed, as shown in Figure 13. Figure 13a provides a diagram of the four positions located at x/R = 0.2, 0.4, 0.6, and 0.8.In Figure 13b-e, the variations of the velocity inside the wind turbine are investigated in detail.Besides, the velocities and lateral position are non-dimensionalised by the inlet velocity, U ∞ , and the rotor radius R, respectively.The influence of different α is investigated in Figure 13b when TSR = 4, TI = 8%, and the shaft is smooth at this time, α is increased from 3.9% to 9%.The influence of surface roughness is discussed in Figure 13c when α = 5%, TSR = 4 and TI = 8%, the k s /d s changes from 10 −4 to 10 −2 .Under the condition of α = 5%, TSR = 4, and the shaft is smooth, we study the effects of TI by changing values of TI from 4% to 12%, the comparison results are shown in Figure 13d.In Figure 13e, the influence of TSRs is shown when α = 5%, TI = 8%, the shaft of wind turbine is smooth; the values of TSRs are 2, 3, and 4, respectively.For a single VAWT system, the factors affecting the power coefficient of the VAWT include the TSRs of the VAWT, shaft diameter to wind turbine diameter ratios (α), the surface roughness of the central shaft, and the turbulent intensity of incoming flow (TI).Therefore, four factors and three levels are considered in the simulations to account for their impacts on the performance of the urban VAWT.Furthermore, the optimum operating conditions of the VAWT for maximizing the performance are obtained by Taguchi method in Section 4. It can be observed from Figure 13b that when α = 3.9% and α = 5%, the impact of rotating central shaft on flow field inside the wind turbine is almost the same.From Figure 13c-e, we can see that the surface roughness, turbulence intensity of incoming flow, and TSRs all have an effect on the downstream flow field of the VAWT.The TSRs have a great impact, the TI has a weak influence, and the surface roughness of the central shaft only influences the wake of the shaft and affects the downstream velocity field.The specific effects of these factors on power output of the VAWT will also be discussed in Section 4.

The Influence of Rotating Central Shaft on the Flow Field inside the Wind Turbine
In order to analyze the variations of velocity inside the wind turbine under different α, the velocity profile at four different positions has been analyzed, as shown in Figure 13. Figure 13a provides a diagram of the four positions located at x/R = 0.2, 0.4, 0.6, and 0.8.In Figure 13b-e, the variations of the velocity inside the wind turbine are investigated in detail.Besides, the velocities and lateral position are non-dimensionalised by the inlet velocity, U∞, and the rotor radius R, respectively.The influence of different α is investigated in Figure 13b when TSR = 4, TI = 8%, and the shaft is smooth at this time, α is increased from 3.9% to 9%.The influence of surface roughness is discussed in Figure 13c when α = 5%, TSR = 4 and TI = 8%, the ks/ds changes from 10 −4 to 10 −2 .Under the condition of α = 5%, TSR = 4, and the shaft is smooth, we study the effects of TI by changing values of TI from 4% to 12%, the comparison results are shown in Figure 13d.In Figure 13e, the influence of TSRs is shown when α = 5%, TI = 8%, the shaft of wind turbine is smooth; the values of TSRs are 2, 3, and 4, respectively.For a single VAWT system, the factors affecting the power coefficient of the VAWT include the TSRs of the VAWT, shaft diameter to wind turbine diameter ratios (α), the surface roughness of the central shaft, and the turbulent intensity of incoming flow (TI).Therefore, four factors and three levels are considered in the simulations to account for their impacts on the performance of the urban VAWT.Furthermore, the optimum operating conditions of the VAWT for maximizing the performance are obtained by Taguchi method in Section 4. It can be observed from Figure 13b that when α = 3.9% and α = 5%, the impact of rotating central shaft on flow field inside the wind turbine is almost the same.From Figure 13c-e, we can see that the surface roughness, turbulence intensity of incoming flow, and TSRs all have an effect on the downstream flow field of the VAWT.The TSRs have a great impact, the TI has a weak influence, and the surface roughness of the central shaft only influences the wake of the shaft and affects the downstream velocity field.The specific effects of these factors on power output of the VAWT will also be discussed in Section 4.

Influence of Central Rotating Shaft on Flow Field near Blade
In order to illustrate the wake effect of the wind turbine's central rotating shaft on the blades, we observe the speed variations of in the downstream position of x/R = 1, as shown in Figure 14a.The comparison results are shown in Figure 14b,c, where it can be seen that the velocity and distance are normalized by the inlet velocity U∞ and wind turbine rotation radius R, respectively.As shown in Figure 14b, we assume that there is no-shaft in the wind turbine, the speed variations near the pressure surface of the blade change slightly.As α increase, the wake of the central rotating shaft exerts increasing influence on the downstream blades.When α > 5%, two peaks appear in the speed curve.As shown in Figure 14c, when α increase to 7%, the blade rotates at an azimuthal angle of 280°, and the influence that is caused by the intersection of the blade trajectory and shaft wake is also very drastic.When α = 9%, the same phenomenon occur.The velocity variations near pressure surface of the blade at an azimuthal angle of 270° and 280° also explains the phenomenon where the instantaneous power of the blade occur a minimum (phenomenon in Figure 7).

Influence of Central Rotating Shaft on Flow Field near Blade
In order to illustrate the wake effect of the wind turbine's central rotating shaft on the blades, we observe the speed variations of in the downstream position of x/R = 1, as shown in Figure 14a.The comparison results are shown in Figure 14b,c, where it can be seen that the velocity and distance are normalized by the inlet velocity U ∞ and wind turbine rotation radius R, respectively.As shown in Figure 14b, we assume that there is no-shaft in the wind turbine, the speed variations near the pressure surface of the blade change slightly.As α increase, the wake of the central rotating shaft exerts increasing influence on the downstream blades.When α > 5%, two peaks appear in the speed curve.As shown in Figure 14c, when α increase to 7%, the blade rotates at an azimuthal angle of 280 • , and the influence that is caused by the intersection of the blade trajectory and shaft wake is also very drastic.When α = 9%, the same phenomenon occur.The velocity variations near pressure surface of the blade at an azimuthal angle of 270 • and 280 • also explains the phenomenon where the instantaneous power of the blade occur a minimum (phenomenon in Figure 7).The pressure coefficient (CoP) distribution of the wind turbine's blade surface is observed in order to further analyze the effect of the power coefficient of wind turbine blades.The pressure coefficient CoP is defined by, The pressure coefficient (CoP) distribution of the wind turbine's blade surface is observed in order to further analyze the effect of the power coefficient of wind turbine blades.The pressure coefficient CoP is defined by, where P is the pressure of the blade surface, ρ is the air density, and U ∞ is the free wind velocity.
Figure 15a shows that when α < 7% and the blade rotates at an azimuthal angle of 270 • , a pressure gradient reversal occurs at the trailing edge of the blade.This has negative impact on the blade's lift coefficient.As α increase, the overall variations of the blade's pressure gradient exhibit a decreasing tendency, which leads to the instantaneous power coefficient exhibiting a decreasing trend when the blade rotates to the vicinity of an azimuthal angle of 270 • , as shown in Figure 7. Similarly, as shown in Figure 15b, when α = 9%, the instantaneous power coefficient is much less than that of the blade when α = 7%, due to the blade surface pressure reversing in advance, and because the overall pressure gradient of the blade only changes slightly.
Energies 2018, 11, x 18 of 24 2 0.5 where P is the pressure of the blade surface, ρ is the air density, and U∞ is the free wind velocity.Figure 15a shows that when α < 7% and the blade rotates at an azimuthal angle of 270°, a pressure gradient reversal occurs at the trailing edge of the blade.This has negative impact on the blade's lift coefficient.As α increase, the overall variations of the blade's pressure gradient exhibit a decreasing tendency, which leads to the instantaneous power coefficient exhibiting a decreasing trend when the blade rotates to the vicinity of an azimuthal angle of 270°, as shown in Figure 7. Similarly, as shown in Figure 15b, when α = 9%, the instantaneous power coefficient is much less than that of the blade when α = 7%, due to the blade surface pressure reversing in advance, and because the overall pressure gradient of the blade only changes slightly.

Power Coefficients of VAWT in Taguchi Approach
Details of factors, control parameters, and levels are shown in Table 5.If all cases are taken into account, a total of 3 4 (=81) runs are required.The Taguchi method uses an L9 (3 4 ) orthogonal array to

Power Coefficients of VAWT in Taguchi Approach
Details of factors, control parameters, and levels are shown in Table 5.If all cases are taken into account, a total of 3 4 (=81) runs are required.The Taguchi method uses an L 9 (3 4 ) orthogonal array to analyze the performance of a single VAWT system where only nine runs are needed, as shown in Table 6.The predicted values of C p and S/N ratio are shown in Table 7.The higher the value of C p , the higher the S/N ratio corresponds to a better performance of VAWT.The maximum value of C p is 0.26224, which occurs at Run 7, whereas the minimum value of −0.1359 is found at Run 1, the S/N ratios of Run 7 and Run 1 are 9.610 and 2.747, respectively.These results clearly demonstrate that the combination of factors and levels in Run 7 is beneficial to the performance of VAWT, and the combination of factors and levels also play an important role in terms of improving the power output of the VAWT.According to the values of the S/N ratios that were obtained in Table 7, the profiles of the mean S/N ratios of the four factors are shown in Figure 16a.The value of the mean S/N ratio at Factor A and Level 1 can be calculated by the values of the S/N ratio from the three Level 1 values of Factor A in Table 6 (Runs 1 to 3), namely, the mean S/N ratios of Factor A and Level 1 can be calculated by the average of 2.747, 2.879, and 2.844.For another instance, in order to calculate the value of the mean S/N ratio at Factor B and Level 2, as shown in Table 6, the mean S/N ratios of Factor B and Level 2 can be calculated by the average of 2.879, 5.792, and 8.810.The mean S/N ratios of other Factors and Levels are also calculated as the same method.In order to investigate the degree of influences on the performance of VAWT of different factors, the effect value of the each factor is calculated by the maximum mean S/N ratio minus the minimum mean S/N ratio.The effect values of the four factors are shown in Figure 16b.We can infer from Figure 16b that the influence strength of the effect of rotating shaft on the performance of VAWT is ranked as: TSRs > k s /d s > α > TI.As shown in Figure 16b, the effect value of Factor A is much higher than that of other factors.This means that as the influence of TSRs act on the whole process of operation of VAWT, the TSRs play an important role in the power output of VAWT, the similar conclusion can be also found in the literature [53], for VAWT, there is a strong dependence on for TSRs.In addition, the effect value of Factor C is also higher than that of Factor B and D. This is mainly because the boundary layer separation point can be delayed by increasing the ks/ds, and this results in the downstream blades being less affected by the shedding vortices.However, the influence of ks/ds and α only act on the rotating shaft primarily, so the changes of these two factors mainly affect the flow field behind the shaft.Therefore, the influence strength of these two factors is relatively lower when compared to the influence of TSRs.The effect value of Factor D is the lowest among the four factors.It has also been reported [54] that the effect of TI on the performance of VAWT is weak at low TSRs, but the influence is more evident at higher TSRs.This may be the reason why the TSRs play the most important role on Cp of VAWT, and it is insensitive to Factor D. According to the calculations of factor effect value and the previous research results, it can be seen that the rank of influence strength of these factors in As shown in Figure 16b, the effect value of Factor A is much higher than that of other factors.This means that as the influence of TSRs act on the whole process of operation of VAWT, the TSRs play an important role in the power output of VAWT, the similar conclusion can be also found in the literature [53], for VAWT, there is a strong dependence on for TSRs.In addition, the effect value of Factor C is also higher than that of Factor B and D. This is mainly because the boundary layer separation point can be delayed by increasing the k s /d s , and this results in the downstream blades being less affected by the shedding vortices.However, the influence of k s /d s and α only act on the rotating shaft primarily, so the changes of these two factors mainly affect the flow field behind the shaft.Therefore, the influence strength of these two factors is relatively lower when compared to the influence of TSRs.The effect value of Factor D is the lowest among the four factors.It has also been reported [54] that the effect of TI on the performance of VAWT is weak at low TSRs, but the influence is more evident at higher TSRs.This may be the reason why the TSRs play the most important role on C p of VAWT, and it is insensitive to Factor D. According to the calculations of factor effect value and the previous research results, it can be seen that the rank of influence strength of these factors in this section is precise and reasonable.

Optimum Operation
In order to find the optimum operation of the VAWT in the range of the data we studied in this study, and minimize the negative impact of the rotating shaft on the downstream blades, the profiles of the mean S/N ratio will be further analyzed.Figure 16a reflects that the combination of A3, B1, C3, and D2 is able to maximize the power coefficient of the VAWT within the range that we have studied in Section 4.1, and this combination corresponds to the optimal combination is TSR = 4, α = 5%, k s /d s = 1 × 10 −2 , and TI = 8%.Based on the optimal combination, the value of the power coefficient is 0.262.

Conclusions
In the present study, the effects of the wake of the rotating shaft on the wind turbine performance are investigated numerically.Dynamic analysis of the VAWT is carried out by the method of Computational Fluid Dynamics (CFD).The computational results are obtained by solving the 2D unsteady Navier-Stokes equations using the Transition SST turbulence model.The results are also validated with the experimental results.Moreover, the validation study is also carried out to verify the accuracy of central cylinder simulation.The influence factors, including the tip speed ratios (TSRs) of the operational wind turbine, the turbulence intensity (TI) of the flow, the ratio of the shaft diameter to the wind turbine diameter (α), and the relative surface roughness value k s /d s , are comprehensive considered and studied in detail in this current study.The main conclusions of this current study can be summarized, as follows:

•
A minimum instantaneous power coefficient value appears when the wind turbine blade operates at the vicinity of an azimuthal angle of 270 • , which indicates that the occurrence of the wake effect that is generated by the central rotating shaft.As α increases to 9%, the wind turbine has a minimum power coefficient value, and the power loss of the wind turbine in one revolution increased from 0% to 25%, relative to that of the no-shaft wind turbine.

•
The influence of the wake effect on the downstream velocity field is very strong during the VAWT operation of VAWT.The effect of the shedding vortex, which is exerted by the wind turbine's rotating shaft, becomes stronger and decreases the pressure gradient between the suction side and pressure side.In addition, the lift that is produced by the blades can be also reduced, and this results in the reduction of the C p of the wind turbine in one revolution.

•
As the value of α increase, the phenomenon of boundary layer separation will become increasingly evident.For the wind turbine model in the present study, α = 5% is a critical point, the shedding vortex released by the rotating shaft alternates periodically when α < 5%.The phenomenon is that the distance between the two vortices at the front and back is small.Subsequently, the front shedding vortex exerts a small effect on the back vortex, and the intensity of the vortices is weak.The shedding vortices released by the central rotating shaft are mostly a single main vortex when α > 5%.The distance between the two vortices in the front and back is large, while the influence exertion on each other is weak.

•
In order to investigate the influence strength order of four factors (TSRs, α, TI, and k s /d s ) to the wake effect of central rotating shaft on the performance of VAWT, the Taguchi method has been employed in this study, the analyses of signal-to-noise ratios (S/N ratios) suggest that the influence strength order of the factors is: TSRs > k s /d s > α > TI.In addition, the analysis of the S/N

Figure 2 .
Figure 2. Meshing details: entire mesh in computational domain, mesh in the inner region and mesh near the airfoil.

Figure 3 .
Figure 3. (a) Schematic of the 2D computational domain; and, (b) the mesh of overall computational domain and enlarge grid around the cylinder.

Figure 3 .
Figure 3. (a) Schematic of the 2D computational domain; and, (b) the mesh of overall computational domain and enlarge grid around the cylinder.

Figure 4 .
Figure 4. Schematic of the equivalent sand model.

Figure 4 .
Figure 4. Schematic of the equivalent sand model.

Figure 5 .
Figure 5.Time step effects on torque coefficient for blade one in one revolution.

Figure 5 .
Figure 5.Time step effects on torque coefficient for blade one in one revolution.

Figure 6 .
Figure 6.Comparison of numerical simulation results to experimental data and other Computational Fluid Dynamics (CFD) results.

Figure 6 .
Figure 6.Comparison of numerical simulation results to experimental data and other Computational Fluid Dynamics (CFD) results.

( a )
The quality characteristics and control factors need explicit.(b) Identify the number of levels for the control factors and the interactions between the factors.(c) List the table of orthogonal array according to the numbers of factors and the levels of each factor.(d) Implement the experiments based on the arrangement of orthogonal array.(e) Analyze the experimental results with the signal-to-noise ratio.(f) Obtain the effect strength order of each influent factor.

Figure 7 .
Figure 7.Comparison of instantaneous power coefficient in one revolution with respect to different values of wind turbine diameter (α).

Figure 9 .
Figure 9. Power coefficient and its relative change with respect to α = 0 in one turbine revolution.

Figure 7 .
Figure 7.Comparison of instantaneous power coefficient in one revolution with respect to different values of wind turbine diameter (α).

Figure 7 .
Figure 7.Comparison of instantaneous power coefficient in one revolution with respect to different values of wind turbine diameter (α).

Figure 9 .
Figure 9. Power coefficient and its relative change with respect to α = 0 in one turbine revolution.

Figure 9 .
Figure 9. Power coefficient and its relative change with respect to α = 0 in one turbine revolution.

Figure 9 .
Figure 9. Power coefficient and its relative change with respect to α = 0 in one turbine revolution.

Figure 10 .
Figure 10.Instantaneous vortices contours of central rotating shaft and the time-averaged streamlines for the last 10 operating revolutions of VAWT with respect to different values of α.

Figure 10 .
Figure 10.Instantaneous vortices contours of central rotating shaft and the time-averaged streamlines for the last 10 operating revolutions of VAWT with respect to different values of α.

Figure 11 .
Figure 11.Comparison of central rotating shaft fluctuation amplitude with respect to different values of α.

Figure 11 .
Figure 11.Comparison of central rotating shaft fluctuation amplitude with respect to different values of α.

Figure 12 .
Figure 12.Variation of central rotating shaft lift-drag coefficient ratio with respect to different α.

Figure 12 .
Figure 12.Variation of central rotating shaft lift-drag coefficient ratio with respect to different α.

Energies 2018, 11, x 15 of 24 Figure 12 .
Figure 12.Variation of central rotating shaft lift-drag coefficient ratio with respect to different α.

Figure 13 .
Figure 13.(a) The four measured positions of wind velocity inside the wind turbine; (b) fluctuation of wind velocity in downstream region of x/R = 0.2, 0.4, 0.6, 0.8 at different α; (c) fluctuation of wind velocity in downstream region of x/R = 0.2, 0.4, 0.6, 0.8 at different ks/ds; (d) fluctuation of wind velocity in downstream region of x/R = 0.2, 0.4, 0.6, 0.8 at different turbulence intensity (TI); and, (e) fluctuation of wind velocity in downstream region of x/R = 0.2, 0.4, 0.6, 0.8 at different tip speed ratios (TSRs).

Figure 13 .
Figure 13.(a) The four measured positions of wind velocity inside the wind turbine; (b) fluctuation of wind velocity in downstream region of x/R = 0.2, 0.4, 0.6, 0.8 at different α; (c) fluctuation of wind velocity in downstream region of x/R = 0.2, 0.4, 0.6, 0.8 at different k s /d s ; (d) fluctuation of wind velocity in downstream region of x/R = 0.2, 0.4, 0.6, 0.8 at different turbulence intensity (TI); and, (e) fluctuation of wind velocity in downstream region of x/R = 0.2, 0.4, 0.6, 0.8 at different tip speed ratios (TSRs).

Figure 15 .
Figure 15.(a) Pressure coefficient distribution for blade 1 surface at azimuth angle of 270°; and, (b) pressure coefficient distribution for blade 1 surface at azimuth angle 280°.

Figure 15 .
Figure 15.(a) Pressure coefficient distribution for blade 1 surface at azimuth angle of 270 • ; and, (b) pressure coefficient distribution for blade 1 surface at azimuth angle 280 • .

Figure 16 .
Figure 16.(a) Profiles of mean S/N ratio; and, (b) factor effect value versus factors.

Figure 16 .
Figure 16.(a) Profiles of mean S/N ratio; and, (b) factor effect value versus factors.

Table 1 .
Geometrical characteristics of vertical axis wind turbines (VAWT) and computational domain.

Table 2 .
Comparisons of Strouhal number between present simulation and experimental research.

Table 2 .
Comparisons of Strouhal number between present simulation and experimental research.

Table 3 .
Main detail sizes of three different grid resolutions.

Table 3 .
Main detail sizes of three different grid resolutions.

Table 4 .
Different time steps for verification.

Table 4 .
Different time steps for verification.

Table 5 .
Factors, control parameter, and levels.

Table 7 .
Predicted results of C p and signal-to-noise (S/N) ratio in L 9(3 4) orthogonal array.