Effects of blunt trailing-edge optimization on aerodynamic characteristics of NREL phase VI wind turbine blade under rime ice conditions

Abstract. To reduce the adverse effects of the ice on aerodynamic characteristics, a new NREL Phase VI wind turbine blade which is suitable to rime ice environments is developed through the blunt trailing-edge optimization. The parametric control equations of blunt trailing-edge airfoil are established by adopting the airfoil profile integration theory and B-spline curve, and the curve fitting of the airfoil’s rime ice from LEWICE software is carried out using the linear interpolation algorithm with equidistant and equiangular step lengths. The S809 airfoil under rime ice conditions is optimized to maximize the lift coefficient by the particle swarm optimization (PSO) coupled with GAMBIT and FLUENT, and a NREL Phase VI blade is formed with the optimized airfoil S809-BT (with BT the blunt trailing-edge). The blade’s rime ice is obtained through using the polynomial fitting to deal with projection point coordinates of airfoils’ ice shapes in lagging and flapping surfaces, and the pressure coefficient, flow characteristics, torque and output power of icy sharp and blunt trailing-edge blades are investigated. The results indicate that in rime ice conditions, compared with those of sharp trailing-edge blade, the pressure difference and vortex size of blunt trailing-edge blade become larger, and the torque and output power increase by 4.36 %, 1.55 % and 2.88 % at v 7 m/s, 15 m/s and 20 m/s, respectively. The research provides significant guidance for improving the aerodynamic performance of wind turbine blade considering the icing effects.


Introduction
The wind energy resources are mainly concentrated in the contraction zones of coast and open continent. Owing to the cold climate in these regions, the ice often occurs on the wind turbine blade. The icing decreases the aerodynamic performance, wind energy utilization, and annual power generation [1][2][3][4][5][6], and increases the static and dynamic loads, and vibrational frequency [7,8]. It is necessary to study a reasonable and effective method to improve the aerodynamic and structural characteristics of icy blades.
The airfoil, that is, the cross-section of the blade, directly affects the performance. The blunt trailing-edge modification and direct optimization can improve the aerodynamic characteristics, structural strength and stiffness of rough airfoils. Baker et al. [9] investigated the airfoils with symmetrical trailing-edge thickness through experiments, and found that the lift-drag ratio improved and the sensitivity to leading-edge roughness reduced by increasing appropriate trailing-edge thickness. Yang et al. [10] used the Computational fluid dynamics (CFD) method to study the aerodynamic performance of blunt trailing-edge airfoil, and the results showed that the maximum lift coefficient increased and the effects of the leading-edge roughness on lift characteristics decreased. Zhang et al. [11,12] numerically studied the effects of the blunt trailing-edge modification on the aerodynamic performance enhancement of the airfoil with sensitive roughness height and position.
Miller et al. [13] performed the multi-objective aero-structural optimization of thick flatback airfoils using a genetic algorithm, and found that the obtained airfoils displayed a relative insensitivity to roughness. Chen et al. [14] presented an optimization method of wind turbine airfoil to improve the aerodynamic performance under typical rime ice conditions. Zhang et al. [15] carried out the optimization design of blunt trailing-edge airfoil by PSO coupled with XFOIL, then added a lug boss to obtain the rough airfoil, and numerically analyzed the lift-drag ratios, pressure coefficients and flow characteristics.
The above studies separately analyzed the effects of the blunt trailing-edge modification and direct optimization on aerodynamic characteristics of rough airfoils. In view of the advantages of two methods, this paper carries out the optimization design of blunt trailing-edge airfoil under rime ice conditions and investigates its impact on the aerodynamic performance enhancement of icy blades. Section 2 establishes the control equations of blunt trailing-edge airfoil profile, fits the airfoil's rime ice obtained by LEWICE software, and optimizes the blunt trailing-edge airfoil with rime ice. Section 3 uses the optimized airfoil as the cross-sections from 46.6 % (with the radius of wind wheel) to 75 % to construct the NREL Phase VI blade, and applies the polynomial fitting to deal with projection point coordinates of airfoil's ice shape in lagging and flapping surfaces. Section 4 studies the aerodynamic characteristics of sharp and blunt trailing-edge blades before and after icing.

Expression method of blunt trailing-edge airfoil profile
The coordinates of blunt trailing-edge airfoil profiles before 0.4 (with the chord length) on the upper surface and 0.5 on the lower surface are given by the profile integration theory: where and are the abscissa and longitudinal coordinates of the airfoil, is a quarter of the chord length, and ( ) is the shape function. ( ) = + + + ⋯ + , = 1, 2, 3,…, in which is the argument, , , ,…, are the coefficients, and = 1. The coordinates after 0.4 on the upper surface and 0.5 on the lower surface are given by the B-spline curve: where = 1 and 2 represents the upper and lower profiles, and is relative position of the curve. The last points formed by the profile integration theory, and the blunt trailing-edge points (1, ℎ × ) and (1, −ℎ × (1 − )), go through , (0) and , (1), respectively, where ℎ and are the trailing-edge thickness and the ratio of trailing-edge thickness of upper surface to that of airfoil.

Fitting of airfoil's rime ice
The quantity and location of airfoil's ice shape feature points vary with the airfoil type, chord length, and icing conditions. The linear interpolation algorithm in the first and fourth quadrants with equidistant step length and the second and third quadrants with equiangular step length is used to fit the rime ice so as to remain the number of feature points unchanged and realize the automatical mesh partition in the optimization. The coordinate formulas of interpolation points using the equal angle as the step length are as follows: where ( , ) and ( , ) are the coordinates of rime ice feature points before the interpolation, ( , ) are the coordinates of points obtained after the interpolation, is the -th key point of rime ice shape, and is the angle between the axis and connection line of the interpolation point and coordinate origin.

Optimization design
PSO coupled with GAMBIT and FLUENT is applied to optimize the blunt trailing-edge airfoil under rime ice conditions. The .jou file realizes the parameterized geometric modeling, computational domain establishment and mesh generation in GAMBIT, and the automatic computation in FLUENT.
In order to ensure that the profile has the airfoil's shape characteristics, the boundary constraints of optimization variables are: Moreover, the lift coefficient of optimized blunt trailing-edge airfoil should be higher than that of sharp trailing-edge airfoil : The lift coefficient is a significant index to measure the aerodynamic performance of wind turbine airfoil [17]. So the design objective is to maximize the lift coefficient:

Optimization example
The optimal design is performed for S809 airfoil with the relative thickness of 21 % at 0.395 and relative camber of 0.99 % at 0.823 . According to Ref. [17], the icing parameters are given in Table 1 and is 0.267 m. The .jou file in GAMBIT is compiled to automatically establish the computational domain and generate the grid, and the results are stored in the .msh file. The computational domain consists of a semicircle with a diameter of 30 and a rectangle of size 30 ×20 , and is divided into quadrilateral surfaces ANIH, NMJI, MLKJ, LFGK, FEDG, DCHG and CBAH, as shown in Fig. 1. The chord line of the airfoil is horizontal, and its end coincides with the center of the semicircle. The division method of each edge is given in Table 2, and the seven surfaces are discretized by the C-type structured grid, as shown in Fig. 2. The first boundary layer height is 0.000010512, and = 1. The 200 nodes are arranged on the airfoil with rime ice, and the overall grids consist of 13520 quadrilateral and 13832 nodes. The edges BA, ANMLF and FE are defined as the velocity inlet, and the inlet velocity direction is from left to right. The pressure outlet is applied for BC, CD and DE, and the atmospheric pressure is zero. The edges HIJKG and GH are set to be no-slip and static wall boundary conditions.  The .jou file calls the .msh file to initialize and calculate the flow field while saving the results of lift and drag coefficients. The SST -turbulence model suitable to the wake flow of blunt trailing-edge airfoil is chosen to close the governing equations. The second-order upwind difference scheme is adopted for the discretization of each equation. The SIMPLE algorithm is The initial optimization conditions are as follows: the Reynolds number Re is 10 6 , the lift coefficient of S809 airfoil with rime ice is 0.3133 which is in close agreement with the result of Ref. [18], the sizes of population and dimension of variables are both 20, the learning factors and are both 2, the inertia weight is given by = − ( − ), is the current iteration number, the maximum iteration number is 300, the maximum inertia weight max is 0.9, and the minimum inertia weight min is 0.4. The S809-BT airfoil with the trailing-edge thickness of 0.0445 and thickness distribution ratio of 1:13.35 is obtained, as shown in Fig. 3.

Sharp trailing-edge blade
The NREL Phase VI wind turbine has two blades which use the S809 airfoil from the 25 % radial section to blade tip. The radius of wind wheel is 5.029 m, and the pitch angle of blade tip is 3°. The cut-in wind speed is 5 m/s, the rated power is 19.8 kW, and the rotational speed is 71.63RPM. The airfoil's aerodynamic center is selected as the torsion center. According to the linear distributions of chord length and torsion angle along the span-wise direction in Ref. [19], the conversion formulas of space coordinates of each cross-section are as follows, and the three-dimensional model is shown in Fig. 4:  where and is two-dimensional coordinates of the -th point on each cross-section before the transformation, is span-wise coordinate of the -th cross-section, , , is three-dimensional coordinates of the -th point on the -th cross-section after the transformation, and and are chord length and torsion angle of the -th cross-section, respectively.

Sharp and blunt trailing-edge blades with rime ice
The rime ices of multiple cross-sections are obtained from LEWICE and processed by the linear interpolation algorithm. The feature point coordinates of ice shapes of adjacent sections are mapped into lagging and flapping surfaces, and the space curves of rime ices on the blade surface are formed by using the polynomial fitting deal with projection points.
The blade with the length of is divided into segments along the -axis, and key points are gained after fitting the airfoil's rime ice. So the ice shape space curves along the span-wise direction form, which have discrete points of + 1 group, and the polynomial fitting formulas are as follows: = 1 × 2 + 2 × + 3 , = 4 × 2 + 5 × + 6 , = 1,2, ⋯ , , = 1,2, ⋯ , + 1, where , , are coordinates of the -th point on the -th space curve along the chord-wise, perpendicular to chord-wise and span-wise directions, respectively, and 1 , 2 , 3 ,…, 6 are polynomial fitting parameters. Fig. 5 gives two NREL Phase VI blades with S809 airfoil between 25 % radial section to blade tip and optimized S809-BT airfoil between 46.6 % and 75 % radial sections after icing, respectively.

Computational domain of blade
Only half of the flow field is analyzed because the wind wheel rotates periodically around the -axis. The blade's computational domain is a semi-cylindrical body and divided into three static outer domains B, C, D, and one rotating inner domain A, as shown in Fig. 6. The blade is located in the center of domain A, and the blade root is 0.508 m away from the wheel center.

Grid generation of blade
The division method of edges and arcs of blade's computational domain is given in Table 3. The structured and unstructured quadrilateral grids are used to divide the surface, top and bottom of the blade. The inner and outer domains are divided by unstructured and hexahedral structured grids, respectively. For C domain, the division in the order of the edge, surface and volume is carried out to obtain hexahedral structured grids. The common surfaces between B domain and A domain and between B domain and C domain are defined as mapping source surfaces to establish the volume mesh of B domain, and the same mesh method is applied in D domain. About 4 million grids are generated, as shown in Fig. 7, and is less than 2.

Boundary conditions and solving control parameters
The faces A 1 B 1 C 1 D 1 E 1 F 1 , A 1 A 2 F 2 F 1 , A 2 A 3 F 3 F 2 and A 3 A 4 F 4 F 3 are velocity inlet boundaries, and the face A 4 B 4 C 4 D 4 E 4 F 4 is a pressure outlet boundary. The face B 2 B 3 E 3 E 2 is the interface between stationary and rotating domains, and is set as the interior. The bottom of semi-cylindrical body is set as periodic boundary, and the blade surface is set as static and non-slip wall. The MRF method and SST -model are selected. Based on the pressure solver, the governing equation is solved by the SIMPLE algorithm. The pressure-velocity coupling treatment uses the PRESTO! format. The second-order upwind style is adopted for other discrete schemes. The convergence criteria of all variables are less than 10 -5 . /s. The pressure curve area evidently increases, and the aerodynamic performance improves. This is because that the blunt trailing-edge design reduces the sensitivity of the airfoil to leading-edge roughness, which decreases the flow disturbance induced by the rime ice. Fig. 9 shows when = 7 m/s, there is no flow separation on each cross-section of sharp trailing-edge blade, and 25 % and 30 % radial sections of sharp trailing-edge blade with rime ice. Small separation vortexes occur on the trailing-edge of suction surface at 57 %, 80 % and 95 % radial sections for icy blade, and are similar in size. When = 10 m/s, a vortex starts to occur on the trailing-edge of suction surface at 57 %, 80 % and 95 % radial sections, and its size increases first and then decreases as the cross-section approaches the blade tip before and after icing. The size of vortices at 57 % and 95 % radial sections increases after icing, but basically remains unchanged at the 80 % radial section. When = 15 m/s, there exists a vortex on the trailing-edge of suction surface at the 25 % radial section, and a pair of vortices with opposite direction of rotation and large differences in size exists on the suction surface at 30 %, 57 % and 95 % radial sections before and after icing. At the 80 % radial section, the vortex on the suction surface has shed and begun to reform before icing, but after icing, a pair of vortices with opposite rotation direction and large differences in size has newly formed. That is to say, the vortex on the suction surface increases until it sheds and then reforms and gradually increases as the section approaches the blade tip, and the vortex starts shedding in advance after icing.

Flow distribution
For the blunt trailing-edge blade with rime ice, the flow separation has occurred on the trailing-edge of suction surface, but the vortex has not yet formed at = 7 m/s. A pair of vortices which rotate in opposite directions and have large differences in size occurs on the trailing-edge of suction surface at = 10 m/s and 15 m/s, and is larger than that on the sharp trailing-edge blade.  Fig. 10 shows for the torque and output power of sharp trailing-edge blade, the numerical results are slightly larger than experimental data form Ref. [20] due to neglect the effects of support structure and aero-elasticity, and they have similar variance tendency. and quickly increase first, then decrease gradually and increase slowly with the increase of wind speed for sharp and blunt trailing-edge blades before and after icing. and reduce obviously after icing, and the largest decline is 22.78 % at = 15 m/s for two parameters which have the same rate of change according to =

Torque and output power
. Compared with those of sharp trailing-edge blade under rime ice conditions, and of blunt trailing-edge blade increase by 4.36 %, 0.60 %, 1.55 %, 2.88 % and 0.40 % at = 7 m/s, 10 m/s, 15 m/s, 20 m/s and 25 m/s, respectively. Overall, the icing makes and decrease, but the effect of icing can be reduced through the blunt trailing-edge optimization. 6

Conclusions
The present work firstly performs the blunt trailing-edge optimization of S809 airfoil under rime ice conditions, and then the effects of the rime ice on the pressure coefficient, flow distribution, torque and output power of sharp and blunt trailing-edge NREL Phase VI blades are investigated. The main conclusions are as follows: 1) After the sharp trailing-edge blade is covered with rime ice, the pressure difference between upper and lower surfaces decreases at 95 % ( = 7 m/s), 57 % ( = 10 m/s) and 30 %, 80 %, 95 % ( = 15 m/s) radial sections, respectively, and is unchanged at other cross-sections. The flow separation is more obvious, the size of vortices on the trailing-edge increases, and the vortices move toward the leading-edge. and reduce largely, and the biggest decrease is 22.78 % at = 15 m/s.
2) The blunt trailing-edge optimization makes the pressure difference on the blade surface and the size of vortices larger than those of sharp trailing-edge blade under rime ice conditions. Therefore, and both increase by 4.36 %, 1.55 % and 2.88 % at = 7 m/s, 15 m/s and 20 m/s, respectively, and change little for other wind speeds.
Gege Wang received Master degree in School of Mechanical Engineering from Tiangong University, Tianjin, China, in 2019. Now she works at LM Wind Power Company. Her current research interests include optimal design, aerodynamic performance and structural damage analysis.
Wei Li received Ph.D. degree in School of Mechanical Engineering from Tianjin University, Tianjin, China, in 2010. Now he works at Tianjin Chengjian University. His current research interests include optimal design, aerodynamic performance and structural damage analysis.