Optimization on Airfoil of Vertical Axis Wind Turbine Based on CST Parameterization and NSGA-II Aigorithm

Optimizing the NACA0015 airfoil which is widely applied in small-scale vertical axis wind turbine to make it has a better aerodynamic performance. In the optimization process, using CST parameterization method to perturb the airfoil geometry, the thickness and camber of the airfoil are selected as the constraint, and the value of the maximum tangential force coefficient is chosen as the objective function, the genetic algorithm based on non-dominated sorting (NSGA-II)is selected as an optimization method, calculates the aerodynamic performance of the airfoil by applying the approach of combining XFOIL program and Viterna-Corrigan post-stall mode ,and establishes the optimizing process by the optimization software modefrontier for NACA0015 airfoil’s muti-point optimization, validate the airfoil’s performance with CFD finally. The result illustrates that, by comparing with the NACA0015 airfoil, the optimized airfoil’s lift to drag ratio is improved over a wide range of attack angles, the stall performance is more gentle. The maximum lift coefficient, the maximum lift-drag ratio and the maximum tangential force coefficient are increased by 7.5%,9 and 8.87%, respectively. The optimized airfoil has a wide variable condition performance, more suitable for the operating conditions of a vertical axis wind turbine. Finally, predict the rotor efficiency with optimized airfoil and NACA0015 airfoil for different tip speed ratios and different solidities with multiple streamtube model, the result shows the rotor with optimized airfoil has a higher efficiency. Keyword: vertical axis wind turbine; CST parameterization; NSGA-II; airfoil; optimization; multiple streamtube model


Introduction
Energy is the foundation of economic and social development, is the driving force of economic and social progress, closely related to human life and the living environment.Wind power is a clean and renewable energy, and its application technology is more and more mature.Due to the VAWT (Vertical Axis Wind Turbine), especially the H-VAWT(Straight Bladed Vertical Axis Wind Turbine) has many advantages: such as omni directional，slightly vibration，low sound emission, high safety，simply structure，easy to install and convenient to control and keep in repair , that it is commonly and highly concentrated and became a enthusiastically studied item.Wind turbines rely on the blades to draw wind energy, so airfoil aerodynamic performance directly affects the utilization of wind energy of the wind turbine.Sheldahl et al. [1] reported experimental investigations of the aerodynamic characteristics of NACA-0012, NACA-0015, NACA-0018 and NACA-0021 airfoils for use in VAWT; Sheldahl and Klimas thought Vertical-axis wind turbines have typically used symmetric airfoils, notably the NACA four-series (in particular, the NACA 0012, NACA 0015 and NACA 0018).But M.C.Claessens et al. [2]reported common NACA0015 and NACA0018 symmetric airfoil is not entirely suitable for VAWT, the original airfoil must be improved.The NACA0018 airfoil had been improved by M.C.Claessens, the thickness was increased by 2% and the camber was increased by 0.8%, and got the DU 06-W-200 laminar flow airfoil.Under the negative angle of attack, the aerodynamic performance of DU 06-W-200 and NACA0018 was fairly；under the positive angle of attack，DU 06-W-200 airfoil has a higher CLmax and a wide range of low drag; researchers modified the NACA0012 airfoil and showed that increasing the airfoil's camber can improve the stall characteristics [3]；Liu Xiong et al. [4] found that appropriately thickening the airfoil's trailing edge has a little effect on the aerodynamic performance of the wind turbine, the output power was enhanced slightly with the thicking airfoil.
Md Farhad Ismail et al. [5] investigated the effect of profile-modifications on a NACA-0015 aerofoil used in VAWTs (vertical axis wind turbines).The profile-modifications being investigated consist of a combination of inward semi-circular dimple and Gurney flap at the lower surface of the NACA-0015 aerofoil。Carlos Simão Ferreira et al. [6] demonstrated that the optimization based on lift slope was the correct objective function for power performance of the VAWT.Some previous studies on airfoil optimization were focus on maximizing the ratios of the lift to drag coefficients to improve the performance of wind turbine blades [7,8].As tangential force is responsible for the power produced by VAWT, this study maximizes the tangential force coefficient.This paper aims to provide an optimizing approach to optimize NACA0015 airfoil based on the CST parametric method and NSGA-II.As the angle of attack of the blades in a VAWT changes continuously, this paper aims to maximize the tangential force coefficient at seven angles of attack rather than to maximize the lift coefficient at a single angle of attack.The optimized airfoil's aerodynamic performance has been greatly improved, the maximum lift coefficient increased by 7.5%, the maximum lift-drag ratio increased by 9%, the maximum tangential force coefficient increased by 8.87%.The optimized airfoil has a better aerodynamic performance, the stall performance is more gentle.Compared with NACA0015 airfoil, the optimized airfoil has a wide variable condition performance, the rotor with optimized airfoil has a higher efficiency, so the optimized airfoil is more suitable for the operating conditions of a vertical axis wind turbine.

CST Parameterization Method
During the airfoil optimization process, the airfoil geometry shape must be perturbed before the optimization.In order to ensure the perturbed point smoothly and continuously connected, the perturbation should be made on a certain mathematical carrier, that is the parameterization method of the perturbation。The merits of the parameterization method has a very important impact on the final optimization result, it is the key factor to decide the efficiency and result of the optimization.The parameterization method used in this paper is a common method based on class function / shape function transformation(CST)，it was proposed by Kulfan in US Boeing Airplane Company.The CST parameterization method has clear geometric meaning， less control parameters, good adaptable and better accuracy.Kulfan noted that the method can be used to the parameterization of the airfoil coordinate and other geometric features parameterization such as the camber and arc.The CST parameterization method has a good prospect in Aerodynamic optimization of two-dimensional airfoils and three-dimensional wings, so it has also been widely applied to the airfoil optimization of wind turbines [9].
The Bernstein polynomials is used as the shape function of CST parameterization method, which is typically expressed as: For airfoil, the expressions of geometric perturbation parameterization of the upper and lower surfaces are given as: ( ) ( ) Where, the subscripts u, l are denoted as the upper and lower surfaces respectively, yTEU, yTEL are denoted as the y coordinates of the upper and lower surfaces trailing edge.C(x) is the class function, it is used to limit the airfoil's type, different types of airfoil can be gotten by adjusting the value of N1、N2，the class function is defined as： for the round leading and sharp trailing airfoil , the value of N1 and N2 are adopted as 0.5 and 1.0 respectively.S(x)is the Shape function used to modify the basic shape formed by the class function, it is given as: ( ) ( ) Where, Aui, Ali are the coefficients to be determined，Si(x) is the Bernstein polynomials.The airfoil shape can be determined through determining the coefficients with the least squares fitting.

The Genetic Algorithm Based on Non-Dominated Sorting(NSGA-II)
Genetic algorithm is a semi-random search optimization algorithms based on the "survival of the fittest" of Darwinian, its specific mechanism constituted with the selection, crossover and mutation operations make the optimization has a robustness.The genetic algorithm based on non-dominated sorting was proposed by professor Srinivas and Deb, it has a better distribution than PAES and SPEA algorithm(another two multi-objective genetic algorithm with elitist strategy, these two algorithms have the advantage of creating a variety of Pareto optimal levels), and its convergence is closer to the actual Pareto optimal level [10].The comparison was made between the NSGA-II and PAES and SPEA algorithm under the same test functions, the test functions were given as: MOP1： ( ) ( ) MOP2: ( ) ( ) ( ) ( ) ( ) ( ) Where，A1=0.5sin1-2cos1＋sin2-1.5cos2 A2=1.5sin1-cos1＋2sin2-0.5cos2 B1=0.5sinx-2cosx＋siny-1.5cosy B2=1.5sinx-cosx＋2siny-0.5cosy MOP3: ( ) ( ) ( ) ( ) ( ) ( ) Comparing the real distance to the Pareto optimal boundary and its standard deviation of three algorithms, the result was calculated as the following table: It can be seen from the above table that the convergence of NSGA-II was closer to the actual Pareto optimal level and the deviation was smaller.

Process
The aerodynamic calculating of airfoil is the important step to do the airfoil's optimum designing.
For the small VAWTs， they usually work under low Reynolds number (Re=0.5×10 5 ～3×10 5 ) and the related airfoil adopted is very sensitive to the changing of the turbulence strength, the roughness of the airfoil surface and it's self vibration, etc. XFOIL is a viscous and non-viscous iteration programme which is widely used for analyzing the aerodynamic performance of airfoil with low Reynolds number.It could not only solve the nonlinear couple of the single bubble which occurred among the viscous, transition and non-viscous formulas, but also could illustrate the complex physical phenomena such as the transiting bubble, etc. NACA63-615 、 NACA63-415 、 NACA63-421 and other five airfoils were calculated by XFOIL and experimented in the wind tunnel, the result showed that the datas obtained by the two methods had a good coherency before stall.
When airfoil in stall condition, the aerodynamic performance is similar to the flat disturbed flow in a high angle of attack and no longer depend on airfoil's geometry.Therefore, the Viterna-Corrigan post-stall mode is used to calculate the aerodynamic performance after the airfoil stall.Literature [11] showed that the datas calculated with Viterna-Corrigan post-stall mode had a good fit with the experimental data.The Viterna-Corrigan post-stall mode is specifically described as: ( ) sin cos 1.11 0.018 sin 1.11 0.018 sin cos 2 c o ss i n 1.11 0.018 sin 1.11 0.018 sin cos cos Where, Cl , Cd are the lift and drag coefficients，Cls , Cds are the lift and drag coefficients of stall，α ， αs are the angle of attack and stall angle of attack，μ is the aspect ratio.

Design variables and constraints
During the Airfoil Optimization Design, the selection of the design variables has a significant impact on the optimization process.The thickness and camber of the airfoil are very important geometrical parameters.Under low Reynolds number, the lift coefficient, drag coefficient and stall angle can be increased by increasing the camber of the airfoil, the sensitivity of the airfoil to the leading edge roughness can be reduced at the same time.Appropriately increasing the airfoil thickness can increase the starting torque of VAWT, enhance the structural strength of the blade, improve the aerodynamic performance of the airfoil.For the VAWT airfoil , M.C.Claessens showed that the airfoil with the 15%-18% thickness which related to the chord length and the camber which did not exceed the value of 6% had better aerodynamic performance.Therefore, the design variables were identified as the airfoil thickness t and camber c, the constraints were 15% ≤d≤18% and 0% ≤c≤6%, the Reynolds number was 300,000, and made the multi-point optimization for NACA0015 airfoil at seven angles of attack.The orientation of the blades relative to the true wind direction is known as the azimuthal angle (θ).The rotation of the VAWT results in a change in the tangential and normal component of the wind (relative to the blade).As a result, the chordal (tangential) velocity (Vc ) and normal velocity (V n ) can be written in terms of the induced velocity (U) and tip speed ratio as Figure 1.

Determination of the objective function
( ) Where, λ=ωR/U∞ is the tip speed ratio.
In the absence of flow restrictors, the induced velocity can be assumed equal to the free-stream velocity (U∞), i.e.U= U∞.The effective wind velocity (W) and the angle of attack (α) are given as: During the airfoil optimization design, lift coefficient, drag coefficient, lift-drag ratio and other parameters can be selected as the optimization goal, the optimization process is to make the selected parameters achieve optimal solutions.
For VAWT, the tangential force coefficient Ct is always selected to evaluate the aerodynamic performance of VAWT airfoil at rated operation。Therefore, the objective function is determined as max Ct .The tangential force coefficient Ct is defined to represent the non-dimensional tangential force in terms of the lift coefficient (Cl ) and drag coefficient (Cd ) as:  The optimized airfoil with the 16.65% thickness and the 2% camber which related to the chord length.The shape of the optimized airfoil and NACA0015 airfoil had been shown in Figure 3.As can be seen from Figure 3, the geometric shapes of the optimized airfoil had obvious changes, the thickness, camber and leading edge radius increased in varying degrees.
The aerofoil having a chord length c=0.4 m is considered, a sufficiently long domain (50c) is chosen to avoid the effects of the outlet boundary condition.The 2D simulated NACA0015 airfoil was meshed with a fully structured C-type grid with 800 nodes over the surface as shown in Figure 4: (a) Airfoil magnified view grid (b) Computational grid  It can be seen from Figure 5 that the two-equation k-ω SST turbulence model gave closer prediction of lift coefficient both in pre-stall and post-stall region than the SA turbulence models hence it was considered the best model.

Multiple Streamtube Model
In order to analyze the aerodynamic performance of the vertical axis wind turbine, the free wake model, rigid wake model and streamtube model had beenbuilt up.The free wake model has been successfully applied to the propeller and airfoil design, while due to its complexity and the calculation requires a lot of computing time, it is unsuitable to predict the performance of a vertical axis wind turbine; As the single streamtube passes through the rotor, the wind velocity is assumed to be everywhere constant.The forces on the airfoil blades are then computed,using this uniform velocity.The wind velocity in the streamtube at the rotor is then related to the undisturbed freestream velocity by equating the drag force on the rotor to the change in fluid momentum through the rotor.While this approach is somewhat elegant in its simplicity and predicts overall performance rather well for lightly loaded blades, it is incapable of adequately predicting information which requires a more precise knowledge of wind velocty variations across the rotor.
Multiple streamtubes model was firstly proposed by Strickland in 1975.The calculation process of the multiple streamtubes model is relatively simple and has a higher accuracy than the single streamtube [13].rotor phase angle is θ.The free stream velocity is disturbed when passes through the streamtube, the velocity after passing through the streamtube is denoted by V.
Using the time averaged streamwise momentum equation and Glauert blade elements, the average streawise force ⎯Fx exerted by the blade elements as they pass through the streamtube can be written as: ( ) Where ρ is the fluid density and AS is the streamtube cross-sectional area.The average force⎯Fx in the streamtube can be related to the streamwise force Fx exerted by an individual blade element as it passes through the streamtube by noting that each of N blade elements spend Δθ/π percent of their time in the streamtube.Therefore the average force can be denoted as: Eliminating⎯Fx from equations 22 and 23 yields: As the Figure 7 shows, the streamwise force Fx can be decomposed into the force FT which is tangent to the airfoil chord line and FN which is normal to the chord line.So the resultant strearmwise force Fx can be given by: ( ) Where the meridional angle δ indicates the angle between the blade and the horizontal plane.
The forces FT and FN can be expressed in terms of the fluid density ρ, the airfoil chord length c and the relative velocity W: In non-dimensional form, these forces can be written as: where VT is the maximum tip speed at the rotor equator.
The coefficients CT and CN are related to the more common airfoil lift and drag coefficients CL The angle of attack and associated relative velocity in the plane of the airfoil cross section can be obtained from figure 8.The angle of attack α is given by: sin sin arctan cos As the figure 11 shows, the relative velocity W in the plane of the airfoil cross section can be expressed as: For convenience, the left hand side of equation 23 can be denoted by Fx * : Defining an interference factor by: The streamwise momentum equataion can be denoted as: The formula can calculate the rotor induced velocity based on the iterative solution of the streamtube momentum equation.
Once the streamtube momentum equation has been solved, the torque produced by a rotor blade element as it passes through the streamtube can be obtained by: The torque on a complete blade is thus given by: Where NS is the number of blade segments, each blade segment is assumed to be of a length Δh/sinδ.
To obtain the average torque produced on the rotor by all of the N blades, the value of TB must be time averaged and multiplied by N .If values of TS are obtained at Nt values of θ in increments of π/Nt , then the average rotor torque can be written as: The calculations were made at every 10° intervals in θ and at intervals in Z equal to one-tenth of the rotor height, where Z is measured from the rotor base along the vertical axis.Therefore, NS = 10 and Nt = 19.
The rotor power coefficient in terms of the average rotor torque is given by: Combining equations 37,38,39 and 40, one obtains: Where X=Rω/V∞ is the tip speed ratio of the local radius r.
The computer program code associated with Strickland multiple streamtube model known as DART, DART is used to calculate the performance of small wind turbine .Wind tunnel tests of two 2-meter diameter rotors were conducted in the LTV wind tunnel in May, 1975.A three bladed rotor with a value of NC/R = 0.27 were tested.The aluminum rotor blades were NACA 0012 airfoils.The tests were conducted with freestream velocities of 7, 9, and 11 meters per second.For the 9 meter per second windspeed, blade Reynold numbers on the rotor tip range from about 0.10×10 6 to 0.36×10 6 for tip to wind speed ratios of 2 and 7 respectively.Data to be used in the DART model were selected from reference [14].The comparison of DART with test results for a blade Reynold number of 0.3×10 6 is shown in Figure 9. Figure 9 shows the relatively good agreement between wind tunnel measurement of the rotor power coefficient and the DART model predictions for a blade Reynold number of 0.3×10 6 .The failure to agree exactly on the left hand portion of the curves is at least partially due to the difference in blade Reynold numbers between test and analysis.On the right hand side of the curve, the test Reynold numbers at the rotor tip and the Reynold numbers used in the DART analysis are nearly the same.The DART prediction is again somewhat high which may be in part due to blade Reynold numbers toward the rotor hub which are again less than that used in the analysis.

Results and discussion
Applying  of pressure between the upper and lower surface.When the attack angle is zero, the pressure of NACA0015 airfoil upper and lower surface is symmetry, so there is no lift; while the pressure of the optimized airfoil lower surface is larger than upper surface, so the optimized airfoil has a lift.The pressure distribution on the airfoil lower surface showed positive values which produce positive lift force when the angle of attack is larger than zero.Pressure coefficient plots shows that there is a high pressure at the leading edge and low at trailing edge.The pressure distribution on the airfoil lower surface shows positive values which produce positive lift force while the upper surface shows positive values which produce negative lift force.Due to the optimized airfoil's curvature of the leading edge and the camber of trailing edge are increased, the pressure difference on the optimized airfoil upper and lower surface is larger than NACA0015, so the lift coefficient of the optimized airfoil is greater than NACA0015.In the postmedian of the airfoil, the pressure difference on the optimized airfoil upper and lower surface has a smooth transition, this indicates that the load pressure gradient is reduced uniformly from the airfoil's center to the trailing edge, which would make the airfoil have a good mechanical performance.In addition, the pressure gradient on For different tip speed ratios and different solidities, the rotor efficiency with optimized airfoil is higher than NACA0015 airfoil.When the tip speed ratio is less than 1, the effect is not obvious; when the tip speed ratio is greater than 1, the rotor efficiency with optimized airfoil is significantly higher than NACA0015 airfoil; especially when the tip speed ratio is 2-4, the efficiency increased considerably; when the tip speed ratio is greater than 4, the increase has been reduced.When the value of NC/R is 0.18, the maximum efficiency of the rotor is gotten at the tip speed ratio of 4.5; while the value of NC/R is 0.27, the maximum efficiency of the rotor is gotten at the tip speed ratio of 4; the increase of the highest efficiency point is 4.88% and 9.5% respectively.Thus, the rotor with optimized airfoil has a higher efficiency, the optimized airfoil is more suitable for vertical axis wind turbine operation.

Conclusions
This paper has proposed a process to optimize the NACA0015 airfoil to get the maximum Ct This paper uses the joint optimization approach with airfoil camber and thickness, the thickness and camber of the optimized airfoil has a increase compared with NACA0015 airfoil, the camber and thickness of the leading edge also increased.After the optimization, the differential pressure of Upper and lower surfaces increased, the lift performance is improved and the maximum lift coefficient increased by 7.5%.Due to the multi-point optimization, therefore the optimized airfoil having a wide variable condition performance, the lift to drag ratio has also been improved over a wide range of the angle of attack, the maximum lift-drag ratio increased 9%.The stall angle of attack is both 15 degrees of the two airfoils, while the separation point of the optimized airfoil is backward and the separation area decreases, the optimized airfoil stalled more gently, the stall performance is greatly improved.After the airfoil stalled, the lift coefficient, lift to drag ratio and tangential force coefficient of the optimized airfoil are much higher than NACA0015 airfoil, the optimized airfoil is more suitable for VAWT.
Finally, this paper verifies the accuracy of the multiple streamtube model and predict the rotor efficiency with optimized airfoil and NACA0015 airfoil for different tip speed ratios and different solidities.The result shows the rotor with optimized airfoil has a higher efficiency, the optimized airfoil is more suitable for vertical axis wind turbine operation.

Conflict of Interest
The authors declare no conflict of interest.

Figure 1 .
Figure 1.Force analysis of a vertical axis wind turbine

20) 1 . 4 . 3
Establishment of optimization processModefrontier is a general multi-objective and multi-disciplinary optimization software developed by the Italian ESTECO company, it established the corresponding optimization process based on the differences in the actual optimization problems.The entire workflow including data flow and logic flow: the data flow is the progress to get the output variables, input variables are imported through the data interface of the integrated module, after the calculations are complete, the data are exported to the next integrated modules to calculate; the logic flow is the progress to generate the input variables, control the optimization process and assess the output variables.The decision variables, constraints and objective functions of the optimal mathematical model are correspond to input variables, constraints and output variables of the workflow[12].The airfoil optimization process has been shown below.

Figure 2 .
Figure 2. The optimizing process of modefrontier

Figure 4 .
Figure 4. Computational and Airfoil magnified view grids

Figure 5 .
Figure 5.Comparison of CFD and experimental lift coefficients of NACA0015 airfoil

Figure 9 .
Figure 9. comparison of DART with Sandia test data k-ω SST model, simulations for the pressure distribution and flow field structure on NACA0015 airfoil at wind speed of 11 m/s were carried out at various angles of attack.The simulation results were analyzed in various stages by the CFD-post features in Tecplot 360.

Figure 13 Figure 14 .
Figure 13 comparison of lift coefficients, lift-drag ratios, tangential force coefficient with CST and NSGA-II.In this paper, CFD simulations are used to obtain the data needed after the optimization (at Reynolds number，Rec was 300000.The turbulence model used in the CFD study was k-ω SST model because it gave closer prediction of lift coefficient both in pre-stall and post-stall region.

Table 1
the real distance to the Pareto optimal boundary and its standard deviation