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

2.1. Expression method of blunt trailing-edge airfoil profile

The coordinates of blunt trailing-edge airfoil profiles before 0.4$c$ (with $c$ the chord length) on the upper surface and 0.5$c$ on the lower surface are given by the profile integration theory:

where $x$ and $y$ are the abscissa and longitudinal coordinates of the airfoil, $a$ is a quarter of the chord length, and $\rho \left(\theta \right)$ is the shape function. $\rho \left(\theta \right)=C_0+C_1\theta +C_2{\theta }^2+\cdots +C_r{\theta }^r$, $r=$1, 2, 3,…, in which $\theta$ is the argument, $C_0$, $C_1$, $C_2$,…, $C_r$ are the coefficients, and $C_0=\text{1}$.

The coordinates after 0.4$c$ on the upper surface and 0.5$c$ on the lower surface are given by the B-spline curve:

where $s=$ 1 and 2 represents the upper and lower profiles, and $t$ is relative position of the curve.

The last points formed by the profile integration theory, and the blunt trailing-edge points $(1,h\times k)$ and $(1,-h\times \left(1-k\right))$, go through $P_{\mathrm{0,3}}^s\left(0\right)$ and $P_{\mathrm{0,3}}^s\left(1\right)$, respectively, where $h$ and $k$ are the trailing-edge thickness and the ratio of trailing-edge thickness of upper surface to that of airfoil. So $P_0^s$ and $P_3^s$ are calculated by the Eqs. (2), and Eq. (3) is obtained:

2.2. 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 $(x_a,y_a)$ and $(x_b,y_b)$ are the coordinates of rime ice feature points before the interpolation, $(x_d,y_d)$ are the coordinates of points obtained after the interpolation, $d$ is the $d$-th key point of rime ice shape, and $\beta$ is the angle between the $x$ axis and connection line of the interpolation point and coordinate origin.

2.3. 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.

The second to twelfth coefficients of shape function [16], the blunt trailing-edge thickness and its distribution ratio, and the B-spline control parameters are selected as design variables:

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 $C_L^{\text{'}}$ should be higher than that of sharp trailing-edge airfoil $C_L$:

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:

2.4. Optimization example

The optimal design is performed for S809 airfoil with the relative thickness of 21 % at 0.395$c$ and relative camber of 0.99 % at 0.823$c$. According to Ref. [17], the icing parameters are given in Table 1 and $c$ 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$c$ and a rectangle of size 30$c$×20$c$, 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 $y^{+}=$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.

Fig. 1Computational domain of airfoil

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 $k$-$\omega$ 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 used to solve the pressure and velocity coupling equation. The convergence criteria of continuity component, velocity component, $k$ and $\omega$ are all 10^-4.

Fig. 2Grids of icy airfoil

a) Global grids

b) Local grids

c) Grids near the leading-edge

The initial optimization conditions are as follows: the Reynolds number Re is 10⁶, 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 $C_1$ and $C_2$ are both 2, the inertia weight $w$ is given by $w=w_{\mathrm{m}\mathrm{a}\mathrm{x}}-\frac{tn}{t{n}_{\mathrm{m}\mathrm{a}\mathrm{x}}}(w_{\mathrm{m}\mathrm{a}\mathrm{x}}-w_{\mathrm{m}\mathrm{i}\mathrm{n}})$, $tn$ is the current iteration number, the maximum iteration number $t{n}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is 300, the maximum inertia weight $w_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is 0.9, and the minimum inertia weight $w_{\mathrm{m}\mathrm{i}\mathrm{n}}$ is 0.4. The S809-BT airfoil with the trailing-edge thickness of 0.0445$c$ and thickness distribution ratio of 1:13.35 is obtained, as shown in Fig. 3.

Fig. 3S809 and S809-BT airfoils

3.1. 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 $R$ 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 $x_v$ and $y_v$ is two-dimensional coordinates of the $v$-th point on each cross-section before the transformation, $z_u$ is span-wise coordinate of the $u$-th cross-section,$x_{uv}^{\mathrm{\text{'}}}$, $y_{uv}^{\mathrm{\text{'}}}$, $z_{uv}^{\mathrm{\text{'}}}$ is three-dimensional coordinates of the $v$-th point on the $u$-th cross-section after the transformation, and $c_u$ and ${\theta }_u$ are chord length and torsion angle of the $u$-th cross-section, respectively.

Fig. 4Three-dimensional model of sharp trailing-edge blade

3.2. 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 $l$ is divided into $m$ segments along the $z$-axis, and $n$ key points are gained after fitting the airfoil’s rime ice. So the $n$ ice shape space curves along the span-wise direction form, which have discrete points of $m+1$ group, and the polynomial fitting formulas are as follows:

where $x_{ij}$, $y_{ij}$, $z_{ij}$ are coordinates of the $j$-th point on the $i$-th space curve along the chord-wise, perpendicular to chord-wise and span-wise directions, respectively, and $p_1$, $p_2$, $p_3$,…, $p_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.

Fig. 5Sharp and blunt trailing-edge blades with rime ice (Blue represents rime ice)

a) Icy sharp trailing-edge blade

b) Part of icy sharp trailing-edge blade

c) Icy blunt trailing-edge blade

d) Part of icy blunt trailing-edge blade

3.3. Computational domain of blade

Only half of the flow field is analyzed because the wind wheel rotates periodically around the $x$-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.

Fig. 6Computational domain of blade

3.4. 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 $y^{+}$ is less than 2.

Fig. 7Grid of icy blade

3.5. Boundary conditions and solving control parameters

The faces A₁B₁C₁D₁E₁F₁, A₁A₂F₂F₁, A₂A₃F₃F₂ and A₃A₄F₄F₃ are velocity inlet boundaries, and the face A₄B₄C₄D₄E₄F₄ is a pressure outlet boundary. The face B₂B₃E₃E₂ 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 $k$-$\omega$ 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.

4.1. Pressure coefficient

Fig. 8(a) and 8(b) show that after the sharp trailing-edge blade is covered with rime ice, the pressure coefficient $C_P$ and the area enclosed by the pressure curve are almost unchanged, except for the pressure coefficients on the suction surfaces of 95 % radial section at $v=$7 m/s and 57 % radial section at $v=$10 m/s. The pressure coefficient on the suction surface increases for $x=$–0.12~0 m and 0.075~0.26 m of 95 % radial section at $v=$7 m/s and $x=$–0.12~0.05 m of 57 % radial section at $v=$10 m/s, and is basically unchanged elsewhere. The pressure curve area reduces a little, which shows that the lift coefficient decreases slightly. In Fig. 8(c), the pressure coefficients and curve areas at 25 % and 57 % radial sections are also unchanged at $v=$15 m/s. For 30 %, 80 % and 95 % radial sections, the pressure coefficient on the pressure surface decreases and that on the suction surface increases at $x=$–0.2~–0.05 m, –0.14~0.24 m and –0.12~0.24 m, respectively, and so the curve area evidently decreases. The reason of these phenomena is that the rime ice makes the blade surface be rough, which increases the flow disturbance.

Under rime ice conditions, compared with the sharp trailing-edge blade, the pressure coefficient on the pressure surface of blunt trailing-edge blade increases at $v=$7 m/s, decreases at $v=$10 m/s, and increases for $x=$–0.02~0.34 m at $v=$15 m/s. The pressure coefficient on the suction surface decreases for $x=$–0.18~–0.08 m and 0.09~0.4 m at $v=$7 m/s, $x=$ –0.04~0.31 m and 0.37~0.4 m at $v=$10 m/s, and is almost unchanged at $v=$15 m/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.

4.2. Flow distribution

Fig. 9 shows when $v=$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 $v=$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 $v=$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.

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 $v=$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 $v=$10 m/s and 15 m/s, and is larger than that on the sharp trailing-edge blade.

Fig. 8Pressure coefficients at different cross-sections before and after icing

a)$v=$7 m/s

b)$v=$10 m/s

c)$v=$ 15 m/s

Fig. 9Flow distribution at different cross-sections before and after icing

a) Sharp trailing-edge blade ($v=$7 m/s)

b) Icy sharp trailing-edge blade ($v=$7 m/s)

c) Icy blunt trailing-edge blade ($v=$ 7 m/s)

d) Sharp trailing-edge blade ($v=$10 m/s)

e) Icy sharp trailing-edge blade ($v=$10 m/s)

f) Icy blunt trailing-edge blade ($v=$10 m/s)

g) Sharp trailing-edge blade ($v=$15 m/s)

h) Icy sharp trailing-edge blade ($v=$15 m/s)

i) Icy blunt trailing-edge blade ($v=$15 m/s)

4.3. Torque and output power

Fig. 10 shows for the torque $T$ and output power $P$ 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. $T$ and $P$ 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. $T$ and $P$ reduce obviously after icing, and the largest decline is 22.78 % at $v=$15 m/s for two parameters which have the same rate of change according to $P=T\omega$. Compared with those of sharp trailing-edge blade under rime ice conditions, $T$ and $P$ of blunt trailing-edge blade increase by 4.36 %, 0.60 %, 1.55 %, 2.88 % and 0.40 % at $v=$7 m/s, 10 m/s, 15 m/s, 20 m/s and 25 m/s, respectively. Overall, the icing makes $T$ and $P$ decrease, but the effect of icing can be reduced through the blunt trailing-edge optimization.

Fig. 10T and P of sharp and blunt trailing-edge blades with rime ice

a)$T$

b)$P$