Accurate RANS Simulation of Wind Turbine Stall by Turbulence Coefficient Calibration

Abstract

Stall, a complex phenomenon related to flow separation, is difficult to be predicted accurately. The motivation of the present study is to propose an approach to improve the simulation accuracy of Reynolds Averaged Navier–Stokes equations (RANS) for wind turbines in stall. The approach is implemented in three steps in simulations of the S809 airfoil and the NREL (National Renewable Energy Laboratory) Phase VI rotor. The similarity between airfoil and rotor simulations is firstly investigated. It is found that the primary reason for the inaccuracy of rotor simulation is not the rotational effect or the 3-D effect, but the turbulence-related problem that already exists in airfoil simulation. Secondly, a coefficient of the SST turbulence model is calibrated in airfoil simulation, ensuring the onset and development of the light stall are predicted accurately. The lift of the airfoil in the light stall, which was overestimated about 30%, is reduced to a level consistent with experimental data. Thirdly, the calibrated coefficient is applied to rotor simulation. That makes the flow patterns on the blade properly simulated and the pressure distribution of the blade, as well as the torque of the rotor, are predicted more accurately. The relative error of the predicted maximum torque is reduced from 34.4% to 3.2%. Furthermore, the procedure of calibration is applied to the MEXICO (Model Experiments in Controlled Conditions) rotor, and the predicted pressure distributions over blade sections are better than the CFD (Computational Fluid Dynamics) results from the Mexnext project. In essence, the present study provides an approach for calibrating rotor simulation using airfoil experimental data, which enhances the potential of RANS in accurate simulation of the wind turbine aerodynamic performance.

1. Introduction

Aerodynamics is one of the most important topics in wind turbine technology. Accurate calculation of the aerodynamic loads is essential for the power prediction and the structural design. A large number of computational tools has been developed based on the blade element momentum (BEM) theory [1,2,3,4], the prescribed/free vortex wake models (VWM) [5,6], or the computational fluid dynamics (CFD) method [7,8]. BEM and VWM are specially developed for rotary machinery, such as propellers, helicopter rotors, and wind turbine rotors. CFD is a general method that can be used to solve almost all common flow problems based on detailed numerical solutions of the Navier–Stokes equations and gives physical quantities of the whole computational domain. This brings many benefits to the computational research of wind turbines. For example, the rotational effect and the three-dimensional (3-D) effect are naturally included in CFD, and therefore the additional models are not required as in BEM. On the other hand, CFD still faces some problems, among which turbulence is the most intractable one, since it is too computationally expensive to be directly simulated.

Turbulence plays an important role in the aerodynamic performance of wind turbines. It exists in the wind, in the blade boundary-layer, and in the wake behind the rotor. There are mainly two ways to deal with turbulence in modern CFD applications. One is the Reynolds Averaged Navier–Stokes equations (RANS) simulation, the other is the Large Eddy Simulation (LES). Although LES provides a much better description of turbulence compared to RANS, it is much more computationally expensive. In order to make the computational cost acceptable, LES is usually used together with the actuator disk (AD) or the actuator line (AL) for wind turbine CFD. The methods of AD-LES and AL-LES have made great progress in the study of rotor wake [9,10]. However, they do not describe the geometry of rotors, and therefore the detailed flow around blades is not simulated. For simulations with blade geometry, RANS simulation is the most practical option.

The simplification of turbulence in RANS increases the inaccuracy of CFD, especially for phenomena directly related to turbulence. The prediction of stall, which is determined by the simulation of turbulence in the boundary-layer, becomes one of the most challenging problems for RANS. Stall is related to flow separation that occurs on the suction side of an airfoil (or blade) when the angle of attack (AoA) exceeds a certain value. It can be divided into two stages of light stall and deep stall, corresponding to trailing edge flow separation and complete flow separation, respectively. Light stall decreases the slope of lift curve, while deep stall causes the lift to drop sharply from the peak value.

A typical example of the inaccuracy of stall prediction is the RANS simulations of the S809 airfoil [11]. A large number of simulations using various turbulence models have been carried out since the end of the twentieth century. As early as 1994, Yang [12] made RANS simulation of the S809 airfoil using a 2-D incompressible solver and the Baldwin–Lomax algebraic turbulence model. The predicted lift coefficient significantly exceeds the experimental data when light stall occurs. Researchers later attempted to improve the simulation by using various turbulence models such as k-ω used by Yang [13], S-A, SST and SA-DES used by Benjanirat [14], and SST and RNG k-ε used by Guerri [15]. Their results repeatedly showed the same inaccuracy: predicting a late onset of stall and therefore overestimating the maximum lift coefficient. Instead of fully turbulent simulations, laminar/turbulent transitional simulations were carried out by other researchers. Wolfe [16] specified a transition point, while Bertagnolio [17] and Langtry [18] employed the Michel model and the γ-Reθ model [19], respectively. Their results showed that stall prediction depends on the simulation of transition, but it is still not accurate enough.

The problem also appears in RANS simulations of the well-known NREL Phase VI rotor [20] and MEXICO rotor [21]. A large number of simulations, including those in the blind comparison [22], were carried out for the NREL Phase VI rotor using various solvers such as OVERFLOW-D [23], elsA [24], EllipSys3D [25], CFX [26], FLUENT [27], THETA [28], etc. The MEXICO rotor was simulated by participants of the Mexnext project [29] and researchers using a variety of RANS codes [30,31,32]. Many of the above simulations failed to accurately predict the stall. The inaccuracy of stall prediction remains an unsolved problem [33,34].

Stall frequently occurs on stall-regulated wind turbines. Also, it may occur on variable-pitch wind turbines during operation in yaw or in extreme wind conditions, and particularly on the blade inboard sections where blade incidence angle is usually high. Therefore, the accurate stall prediction is of great significance to wind turbine load computation, as well as the optimal design of a wind turbine. The present study is devoted to achieving accurate stall prediction in wind turbine RANS simulations. Firstly, the similarity between airfoil simulation and rotor simulation is investigated, which makes it feasible to obtain accurate rotor simulation by making a calibration in airfoil simulation. Then, a coefficient ${\beta }_{\infty }^{*}$ of the SST turbulence model is calibrated according to airfoil experimental data. Finally, the calibrated ${\beta }_{\infty }^{*}$ is applied to rotor simulation, and more accurate prediction of stall is realized. The above approach enhances the potential of RANS in accurate simulation of the wind turbine aerodynamic performance.

2. Simulation Methods

In the present study, 2-D airfoil simulation and 3-D rotor simulation are carried out using ANSYS FLUENT on the S809 airfoil and the NREL (National Renewable Energy Laboratory) Phase VI rotor, respectively. The involved governing equations are the incompressible RANS equations. A comparison between an incompressible solver and a compressible solver was made by Länger-Möller [28] for the NREL Phase VI rotor. The results of the two solvers were almost identical. Therefore, it is reasonable to ignore the compressibility of air in the present study.

Two kinds of 2-D computational domains are generated for the airfoil simulation. One contains the walls of the wind tunnel’s experimental section, the other does not. The outer boundaries of the domain without tunnel walls are set to more than 20 times the chord length away from the airfoil. In each domain, 430 grid nodes are properly distributed on the airfoil surface, and about 100,000 grid cells are employed. The initial height of the first grid layer measured from the airfoil surface is about 5 × 10⁻⁶ of the chord length. The values of y⁺ (see Figure 1) are found less than 1.0 on most areas of the airfoil surface, indicating that the initial height of the first grid layer is small enough to capture the boundary layer. The computational settings of the airfoil simulation are summarized in

Table 1. Computational settings of the airfoil simulation.

Items	Options	Settings
	Space	2D
Models	Time	Steady
	Viscous	SST k-omega
Pressure-Velocity Coupling	Type	Coupled
	Courant Number	50
Discretization Scheme	Pressure	Second Order
	Momentum	Second Order Upwind
	Turbulent Kinetic Energy	Second Order Upwind
	Specific Dissipation Rate	Second Order Upwind

.

The computational domain for the rotor simulation is shown in Figure 2. A pair of rotational periodic surfaces is defined to reduce the computational domain to a half, taking advantage of the axial symmetry of the rotor geometry. The grid nodes on the two periodic surfaces would overlap one to one if one surface were rotated 180° around the axis. The far boundaries are set to more than 15 times the rotor diameter away from the rotor, with no consideration of wind tunnel walls. The computational settings of the rotor simulation are summarized in

Table 2. Computational settings of the rotor simulation.

Items	Options	Settings
Models	Space	3D
	Time	Steady
	Viscous	SST k-omega
	Motion type	Moving reference frame
	Periodic type	Rotationally Periodic
Pressure-Velocity Coupling	Type	Coupled
	Courant Number	50
Discretization Scheme	Pressure	Second Order
	Momentum	Second Order Upwind
	Turbulent Kinetic Energy	Second Order Upwind
	Specific Dissipation Rate	Second Order Upwind

.

A high-quality mesh with enough cells and properly distributed nodes is essential for a reliable numerical simulation. In the present study, the mesh independency has been tested for both the airfoil and the rotor simulations. The five meshes involved in the rotor test are shown in

Table 3. The meshes involved in the test of mesh independency for the NREL Phase VI rotor.

Mesh ID	M01	M02	M03	M04	M05
Initial height of cells adjacent to blades (m)	5 × 10−6	5 × 10−6	5 × 10−6	5 × 10−6	5 × 10−6
number of Nodes along the blade span	40	45	50	55	60
number of Nodes along the blade chord	140	190	240	290	320
Growth ratio of grid size away from blade	1.2	1.15	1.1	1.08	1.065
Total number of cells (million)	1.36	2.29	4.24	6.97	9.83

, where they are denoted by M01, M02, M03, M04, and M05, respectively. The initial height of the grid cells adjacent to the blade is set to 5 × 10⁻⁶ meter equally for all of the five meshes, according to the criterion of y⁺ < 1.0. The simulation results of M03, M04, and M05 are found very close to each other, and the M03 mesh is preferred, since it has the smallest number of cells among the three. For a further verification of sufficient mesh resolution, the Richardson extrapolation [35,36] is applied to the predicted rotor torque from the simulations using M01, M02, and M03.

Table 4. Results of the mesh independency study based on Richardson extrapolation ¹.

Wind Speed	Mesh ID			Richardson Extrapolation
	M01	M02	M03	p	R*	TE
7 m/s	744.2	756.5	757.4	4.18	0.073	757.5
10 m/s	1597.4	1778.6	1802.6	3.19	0.132	1806.5
15 m/s	850	1309.5	1408.8	2.36	0.216	1439.1
20 m/s	1047.6	1006.4	1001.4	3.33	0.121	1000.7

¹ TE represents the extrapolated torque, R* the ratio of errors, and p the order of accuracy.

shows that the order of accuracy of the Richardson extrapolation is greater than 2 in conditions of various wind speeds. That means the M03 mesh (demonstrated in Figure 3) is fine enough for the present simulation based on the second order upwind scheme, and therefore it is finally adopted.

3. Investigation of the RANS Simulation Error

3.1. Simulation of the S809 Airfoil

The wind tunnel experiments of the S809 airfoil, which were conducted in Delft University of Technology (DUT) [11], Ohio State University (OSU) [37], Colorado State University (CSU) [38], and University of Glasgow [39], respectively, cover a range of Reynolds number (Re) between 0.3 × 10⁶ and 2 × 10⁶. In the present study, full turbulence RANS simulation of the S809 airfoil is firstly made at Re = 1 × 10⁶. It is close to the Reynolds number experienced by the blades of the NREL Phase VI rotor. (The Reynolds number varies between 0.7 × 10⁶ and 1.4 × 10⁶ at the blade root and between 1.0 × 10⁶ and 1.1 × 10⁶ at the blade tip, at the wind speeds between 7 m/s and 25 m/s [40].)

The simulation results of the airfoil lift coefficient (C_L) for a wide range of AoA between 0° and 40° are shown in Figure 4. The experimental data from two different sources are referenced in the figure. One is from the OSU [37] experiment, the other is the measured data of the r/R = 0.63 section of a parked NREL Phase VI blade [41]. The r/R = 0.63 section is close to the mid-span of the blade, and thus its performance can be approximately recognized equal to the 2-D S809 airfoil. The two sets of experimental data are consistent when AoA < 20° but deviate from each other when AoA > 20° where deep stall occurs. It is worth asking why such deviation occurs and which set of data is the appropriate reference for the simulation results in deep stall. Rooij [41] once made an analysis of the parked NREL Phase VI blade. The data collected in his article showed that the deviation does not become more notable at the r/R = 0.8 and r/R = 0.95 sections, although the 3-D effect is more significant there. Therefore, the primary reason for the deviation might not be attributed to the 3-D effect of the blade.

The present study pays attention to the tunnel blockage effect. The size of the test section of the NASA Ames wind tunnel, in which the NREL Phase VI experiment was conducted, is 24.4 m × 36.6 m. The chord length of the NREL Phase VI blade (between 0.355 m and 0.737 m) is very small compared to the wind tunnel, and therefore the tunnel blockage effect can be ignored in standstill cases. The wind tunnel experiment of the S809 airfoil was faced with a completely different situation. Figure 5 shows the present simulation result of the flow field around the airfoil at AoA = 40° in the wind tunnel of OSU. It is found that the actual blockage in deep stall is determined by the size of the vortices rather than the airfoil geometry. Therefore, the blockage effect correction applied to the experimental data was inaccurate, since only the size of the airfoil was considered [37].

Figure 6 shows the lift coefficient from the airfoil simulation with tunnel walls compared with the measured data of the OSU experiment. The blockage effect correction has been applied to the simulation data, according to the method used in the experiment. That leads to an appropriate comparison between the simulation result and the experimental data. The shown curves can be divided into three stages according to the simulation error. At stage “A”, the computational lift coefficient agrees with the measured data and increases linearly with the angle of attack. The lift coefficient is significantly overestimated by the simulation at stage “B” where light stall occurs, while it is underestimated at stage “C” where deep stall occurs. It is worth noting that the underestimation in deep stall is not as serious as that of the conventional comparison shown in Figure 4.

3.2. Simulation of the NREL Phase VI Rotor

Full turbulence simulation of the NREL Phase VI rotor in non-yawed condition is carried out at the wind speeds between 7 m/s and 25 m/s. Figure 7 shows the predicted rotor torque versus the wind speed, with a comparison to the measured data. The curves can also be divided into three stages. At stage “A”, the predicted torque agrees with the measured data and increases linearly with the wind speed. The torque is overestimated at stage “B”, while it is underestimated at stage “C”. That is similar to the results of the airfoil simulation shown in Figure 6, which raises a question: Why is there such a similarity of simulation error between the airfoil lift coefficient (see Figure 6) and the rotor torque (see Figure 7)?

Figure 8 shows the simulated limiting streamlines on the suction side of the NREL Phase VI blade at different wind speeds. At 7 m/s, all streamlines are directed from the leading edge to the trailing edge, indicating the flow on the entire blade is attached and no stall occurs. At 10 m/s, the flow pattern on the blade is a combination of trailing edge flow separation (light stall) and attached flow (no stall). At 15 m/s, it becomes a combination of complete flow separation (deep stall) and trailing edge flow separation (light stall). When the wind speed increases to 20 m/s, complete flow separation (deep stall) covers the entire blade. The flow patterns at five sections of the blade are listed in

Table 5. The flow patterns of five sections of the NREL Phase VI blade at various wind speeds.

Wind Speed	r/R = 0.3	r/R = 0.47	r/R = 0.63	r/R = 0.8	r/R = 0.95
V = 7 m/s	No stall	No stall	No stall	No stall	No stall
V = 10 m/s	Light stall	Light stall	Light stall	Light stall	No stall
V = 15 m/s	Deep stall	Deep stall	Deep stall	Light stall	Light stall
V = 20 m/s	Deep stall	Deep stall	Deep stall	Deep stall	Deep stall

. It is found that the simulation error of the rotor torque in Figure 7 is closely related to the flow patterns. At 7 m/s, no stall occurs, and the torque is well predicted. At 10 m/s, most sections are in light stall, and the torque is remarkably overestimated. At 15 m/s, three sections are under deep stall, while two sections are under light stall, and the overestimation of the torque is reduced. At 20 m/s, deep stall occurs on the entire blade, and the torque is underestimated. The above results indicate that the simulation accuracy of rotor torque is much affected by the flow patterns on the blade. The similarity of the simulation error between the airfoil lift coefficient (see Figure 6) and the rotor torque (see Figure 7) is raised by the fact that light stall leads to overestimation of the aerodynamic force, while deep stall leads to underestimation for both the airfoil and the blade sections.

Figure 9 shows the C_t′ versus the wind speed, in which C_t′ is the dimensionless coefficient of the aerodynamic force contributing positive torque of the blade. The shown three sections of r/R = 0.3, r/R = 0.47, and r/R = 0.95 represent the root, the mid-span, and the tip of the blade, respectively. The predicted curves exceed the experimental ones at the wind speeds between about 9 m/s and 17 m/s. The simulation error at the root section (r/R = 0.3) and the tip section (r/R = 0.95) is found not greater than that at the mid-span section (r/R = 0.47), although the flow at the root and the tip is more complex due to the rotational effect and the 3-D effect. It suggests that the rotational effect and the 3-D effect are not the primary cause of the inaccuracy of the rotor simulation.

3.3. The Similarity between Rotor Simulation and Airfoil Simulation

Figure 10 shows the C_t′ versus AoA at the three sections of r/R = 0.3, r/R = 0.47, and r/R = 0.95. For comparison, it also shows the corresponding curves of the S809 airfoil, which are obtained by replacing the aerodynamic force of each blade section with that of the airfoil. The shadow regions with sloping lines are determined by the positive simulation error (computational data subtracts experimental data). The points S (S′), E (E′), and P (P′) represent the starting, the ending, and the peak of the positive error, respectively. The shadow regions for both the airfoil and the blade sections are Λ-shaped, which is evidence of the similarity of simulation error between them. The shadow region for the root section (r/R = 0.3) is located at the upper right of the shadow for the airfoil, due to the rotational effect, which makes the stall delayed and thus makes the simulation error delayed as well. The two shadow regions at the mid-span section (r/R = 0.47) have large overlapping area, which is further evidence of the similarity of error between the airfoil and rotor simulations. The shadow region for the tip section (r/R = 0.95) becomes smaller than that for the airfoil, suggesting that the 2-D airfoil simulation may not be more accurate than the 3-D rotor simulation. Figure 11 shows the magnitude of the simulation error. The peak error of the blade sections is found lower than that of the airfoil. One reasonable explanation may be that the rotational effect and the 3-D effect produce a gentler slope in the process from light stall to deep stall.

The above results suggest that there is a high degree of similarity of error between the airfoil and the rotor simulations, and the rotor simulation error is approximately at the same level of the airfoil simulation error. That leads to a conclusion that the main source of the rotor simulation error already exists in airfoil simulation. Therefore, it is feasible to achieve better rotor simulation by calibrating the RANS method in airfoil simulation.

4. Simulations with Turbulence Coefficient Calibration

4.1. Closure Coefficients of the SST Turbulence Model

Turbulence is a key determinant of flow separation and stall. The SST turbulence model used in the present study is recognized as one of the best linear eddy viscosity turbulence models providing a relatively good prediction of turbulence for both attached and separated flows. There are two transport equations for the turbulence kinetic energy k and the specific dissipation rate ω in the SST model. More than a dozen closure coefficients are involved in the equations. Their default values are

$$\begin{array}{c}{\sigma }_{k,1}=1.176, {\sigma }_{\omega ,1}=2.0, {\sigma }_{k,2}=1.0, {\sigma }_{\omega ,2}=1.168, {\alpha }_1=0.31, {\beta }_{i,1}=0.075, {\beta }_{i,2}=0.0828, {\alpha }_{\infty }^{*}=1.0,\\ {\alpha }_0=1/9, {\beta }_{\infty }^{*}=0.09, R_{\beta }=8.0, R_k=6, R_{\omega }=2.95, {\xi }^{*}=1.5\end{array}$$

Most of these values were determined by the observed turbulence properties in certain conditions. That limits the universality of the turbulence model. The turbulence simulation may be more suitable in some cases but less applicable in other cases. Calibrating the values of these coefficients is a way to achieve a better result for a given case. After a lot of computational tests, the coefficient ${\beta }_{\infty }^{*}$ is adopted to be calibrated in the present study. It is a key coefficient responsible for the dissipation of k and the production of ω in the transport equations. The variation of the value of ${\beta }_{\infty }^{*}$ eventually leads to a change of the turbulent viscosity in RANS simulations. Therefore, the calibration of ${\beta }_{\infty }^{*}$ can be considered as an appropriate adjustment of turbulent viscosity.

4.2. Airfoil Simulation with Different Values of ${\beta }_{\infty }^{*}$

The airfoil simulation of the S809 airfoil has been made using various values of ${\beta }_{\infty }^{*}$. The optimal value can easily be obtained by comparison with the airfoil experimental data. The simulation results show that ${\beta }_{\infty }^{*}$ = 0.11 is the optimal value for the S809 airfoil at Re = 1 × 10⁶. The predicted lift coefficient is shown in Figure 12 for ${\beta }_{\infty }^{*}$ = 0.09 (the default value) and ${\beta }_{\infty }^{*}$ = 0.11. The results for other values of ${\beta }_{\infty }^{*}$ are not displayed here for simplicity. In general, a higher ${\beta }_{\infty }^{*}$ leads to a lower lift coefficient in the range of light stall.

Table 6. The relative simulation error of lift coefficient for different values of ${\beta }_{\infty }^{*}$.

AoA (°)	8	10	12	14	16	18	20	$\overline{\mathit{\delta }}$
${\beta }_{\infty }^{*}$ = 0.09	5.8%	18.1%	23.0%	18.9%	22.7%	39.4%	82.2%	30.0%
${\beta }_{\infty }^{*}$ = 0.11	2.2%	1.3%	−0.3%	0.5%	6.1%	17.5%	−3.9%	4.5%

lists the relative simulation error defined as

$$\delta =\frac{C_{L,sim}-C_{L,\exp }}{C_{L,\exp }}\times 100\%$$

(1)

in which $C_{L,sim}$ and $C_{L,\exp }$ are the lift coefficients obtained from the simulation and the experiment, respectively. An averaged relative error is defined as following to measure the overall error level in a range of angles of attack,

$$\overline{\delta }=\frac{1}{n}{\displaystyle \sum _{n=1 }^i\left|{\delta }_i\right|}$$

(2)

in which i is the index and n is the count.

It can be seen from the figure and the table that ${\beta }_{\infty }^{*}$ = 0.11 leads to significantly better results in the range of 8° < AoA < 20° where light stall occurs. The averaged relative error is reduced to 4.5% from 30.0%. On the other hand, the calibration results in little change of lift coefficient in the range of AoA > 20° where deep stall occurs.

The pressure distributions over the airfoil at four angles of attack of 8°, 12°, 16°, and 20° are shown in Figure 13. According to the NREL Phase VI experiment, light stall is about to occur at 8°, is in progress at 12° and 16°, and converts to deep stall at 20°. The pressure distributions for ${\beta }_{\infty }^{*}$ = 0.11 are consistent with the experimental data at all the four angles of attack. In contrast, ${\beta }_{\infty }^{*}$ = 0.09 only performs well at 8° where no stall occurs. Figure 14 shows the simulated streamlines around the airfoil at angles of attack of 16° and 20°. At 16°, the flow patterns for ${\beta }_{\infty }^{*}$ = 0.09 and ${\beta }_{\infty }^{*}$ = 0.11 are both trailing edge flow separation (light stall). The difference is that ${\beta }_{\infty }^{*}$ = 0.11 leads to a larger separation vortex. At 20°, ${\beta }_{\infty }^{*}$ = 0.11 and ${\beta }_{\infty }^{*}$ = 0.09 gives complete flow separation (deep stall) and trailing edge flow separation (light stall), respectively.

The airfoil simulation has also been carried out at the Reynolds numbers of 0.75 × 10⁶ and 1.5 × 10⁶, which are the lowest and the highest Reynolds number of the OSU experiment [37], respectively. The results are shown in Figure 15. ${\beta }_{\infty }^{*}$ = 0.11 makes the simulation accuracy significantly improved at the both Reynolds numbers, suggesting that the optimal ${\beta }_{\infty }^{*}$ calibrated at a certain Reynolds number stays in a relatively wide range of Reynolds numbers.

4.3. Rotor Simulation Using the Calibrated Coefficient

${\beta }_{\infty }^{*}$ = 0.11, the optimal value for the S809 airfoil, is applied to the simulation of the NREL Phase VI rotor. The predicted rotor torque versus the wind speed is shown in Figure 16. Both the results for ${\beta }_{\infty }^{*}$ = 0.09 and ${\beta }_{\infty }^{*}$ = 0.11 are consistent with the experimental data when the wind speed is lower than 9 m/s. They diverge from each other at the wind speeds between 9 m/s and 20 m/s. ${\beta }_{\infty }^{*}$ = 0.09 leads to a notable overestimation of the maximum torque, while ${\beta }_{\infty }^{*}$ = 0.11 makes it accurately predicted.

Table 7. The relative simulation error of the rotor torque for different values of ${\beta }_{\infty }^{*}$.

Wind Speed	9	10	13	15	17	20
${\beta }_{\infty }^{*}$ = 0.09	4.7%	34.4%	29.0%	20.2%	0.6%	−9.8%
${\beta }_{\infty }^{*}$ = 0.11	−2.3%	3.2%	−4.5%	−15.7%	−18.8%	−9.6%

shows the relative simulation error of the rotor torque. At 10 m/s, where the maximum torque is achieved for both the predicted curves, the relative error is reduced to 3.2% from 34.4% by the application of ${\beta }_{\infty }^{*}$ = 0.11. At 17 m/s, the result for ${\beta }_{\infty }^{*}$ = 0.09 looks better than that of ${\beta }_{\infty }^{*}$ = 0.11. However, it is a mere coincidence due to the change from positive error to negative error of the curve for ${\beta }_{\infty }^{*}$ = 0.09 at this wind speed. When the wind speed increases to 20 m/s or higher, the two curves overlap again and are lower than the experimental data. The maximum difference of the two predicted curves is located at the wind speed of 10 m/s. The pressure distributions of five typical blade sections at this wind speed are given in Figure 17. It is shown that the application of ${\beta }_{\infty }^{*}$ = 0.11 leads to accurate prediction of the pressure distributions. Taking the section of r/R = 0.47 as an example, ${\beta }_{\infty }^{*}$ = 0.11 makes the onset of deep stall accurately captured and properly gives the flat pressure distribution on the suction side of the blade.

The limiting streamlines on the suction side of the blade at the wind speeds of 10 m/s, 15 m/s, and 20 m/s are shown in Figure 18. As compared to ${\beta }_{\infty }^{*}$ = 0.09, the results for ${\beta }_{\infty }^{*}$ = 0.11 have two differences at 10 m/s. Firstly, the area of trailing edge flow separation (light stall) is larger and the area of attached flow (no stall) is correspondingly smaller. Secondly, a small piece of complete flow separation (deep stall) occurs at about r/R = 0.47, which has been corroborated by the experimental pressure distribution at this section shown in Figure 17. At 15 m/s, a large part of the blade is under complete flow separation (deep stall) for both ${\beta }_{\infty }^{*}$ = 0.09 and ${\beta }_{\infty }^{*}$ = 0.11. The difference is that the area of deep stall is larger for ${\beta }_{\infty }^{*}$ = 0.11. At 20 m/s, complete flow separation (deep stall) covers the whole blade, and there is little difference between the two results for ${\beta }_{\infty }^{*}$ = 0.09 and ${\beta }_{\infty }^{*}$ = 0.11. The above results suggest that the simulation using ${\beta }_{\infty }^{*}$ = 0.11 predicts an earlier onset of trailing edge flow separation (light stall) and also an earlier conversion to complete flow separation (deep stall) compared to the simulation using ${\beta }_{\infty }^{*}$ = 0.09.

5. Verification of the Calibration to the MEXICO Rotor

5.1. Introduction to the MEXICO Rotor

In order to study how robust the present approach of calibration is for different airfoils and rotors, the calibration is applied to RANS simulation of the turbine used in the MEXICO (Model Experiments in Controlled Conditions) project [21]. The MEXICO turbine consists of a three-bladed rotor with a diameter of 4.5 m. Three different airfoils of DU91-W2-250, RISØ-A1-21, and NACA 64-418 were used in the blade design. The DU91-W2-250 airfoil was applied from 20% to 45.6% span, the RISØ-A1-21 airfoil from 54.4% to 65.6% span, and the NACA 64-418 airfoil outboard to 74.4% span. The turbine was tested in the 9.5 × 9.5 m² open section of the largest European wind tunnel, the Large-Scale Low Speed Facility (LLF) of the German Dutch Wind Tunnels (DNW), in December 2006. The tests were performed at both non-yawed and yawed flows at different wind speeds and resulted in a database of combined blade pressure distributions, loads, and flow field measurements. As a follow up of the first MEXICO campaign, the New MEXICO measurement on the same turbine was carried out in the same wind tunnel between 20th June and 4th July 2014. Several open questions from the first campaign have been resolved, and good agreement has been found between the first and the new measurements [42].

An international research collaborative project (Mexnext) within the framework of IEA Task 29 was created in June 2008. It was carried out in three phases. The first phase from September 2008 until December 2011 included an assessment of the measurement uncertainties and a validation of different categories of aerodynamic models. In the second phase from January 2012 until the end of 2014, unexplored aerodynamic measurements on wind turbines (both in the wind tunnel as well as in the field) were analyzed from a wide variety of sources. In the third phase from January 2015 until December 2017, several rounds of comparisons were made between the New MEXICO measurements and the results from various codes. Some comparisons between CFD simulations and measurements are reported in [43,44,45].

5.2. Introduction to the Present Simulation

The present study investigates the performance of the SST turbulence model using different values of ${\beta }_{\infty }^{*}$ in a non-yawed case with a rotational speed of 424.5 rpm and a wind speed of 24 m/s. Flow separation occurs on the rotor blades in this case. The separation point location was generally predicted closer to the trailing edge than the experimental value, and the suction side pressure was over-predicted by most of the CFD codes used in the Mexnext project [29].

Pressure distributions were measured in the first and the New MEXICO experiments over five sections located at 25%, 35%, 60%, 82%, and 92% span of the blade. Information about the five sections and their Reynolds numbers in the present case (424.5 rpm, 24 m/s) are shown in

Table 8. Information about the five sections and their Reynolds numbers in the present case.

r/R	r (m)	Chord (m)	Profile	Re
0.25	0.563	0.224	DU91-W2-250	0.38 × 106
0.35	0.788	0.193	DU91-W2-250	0.46 × 106
0.6	1.350	0.142	RISØ-A1-21	0.58 × 106
0.82	1.845	0.113	NACA 64-418	0.64 × 106
0.92	2.070	0.099	NACA 64-418	0.62 × 106

.

The detailed experimental aerodynamic coefficients of the airfoils DU91-W2-250, RISØ-A1-21, and NACA 64-418 were given in the project of Mexnext. The available experimental data of DU91-W2-250 and NACA 64-418 were measured at Reynolds numbers of 0.5 × 10⁶ and 0.7 × 10⁶, respectively. They are close to the operating Reynolds numbers of the corresponding sections shown in Table 8. However, the data of RISØ-A1-21 are available only at Reynolds number of 1.6 × 10⁶, which is much higher than the operating Reynolds number of the corresponding section of r/R = 0.6. That may reduce the reliability of the present calibration for the blade segment based on the airfoil RISØ-A1-21.

The grid layout and density for the airfoils are identical to those for the S809 airfoil. Instead of another mesh independency study, the mesh for the MEXICO rotor is further refined compared to the mesh for the NREL Phase VI rotor. The number of grid cells is increased to about 13 million, which is more than three times that of the NREL Phase VI rotor. All other computational settings for the airfoils and the MEXICO rotor are identical to those for the S809 airfoil and the NREL Phase VI rotor, respectively.

5.3. Results and Disscussion

RANS simulations of the three airfoils of the MEXICO blade are carried out, and the resulted lift coefficients are shown in Figure 19. For DU91-W2-250 at Re = 0.5 × 10⁶, the default ${\beta }_{\infty }^{*}$ of 0.09 gives results consistent with the measured data at a wide range of angles of attack. For RISØ-A1-21 at Re = 1.6 × 10⁶, ${\beta }_{\infty }^{*}$ = 0.10 looks better. For NACA 64-418 at Re = 0.7 × 10⁶, ${\beta }_{\infty }^{*}$ = 0.105 gives perfect prediction before deep stall. As a result, in the simulations of the MEXICO rotor, different values of ${\beta }_{\infty }^{*}$ should be applied at different span locations depending on the airfoil used there. That would be an inconvenience for application. A function might have to be established to give a continuous change of ${\beta }_{\infty }^{*}$ in the computational domain, especially along the blade span. The function could be applied to the SST model if an open source solver or a house code were used. However, it cannot easily have been realized in ANSYS FLUENT. As a compromise, three independent rotor simulations are carried out, in which the ${\beta }_{\infty }^{*}$ values of 0.09, 0.10, and 0.105 are applied, respectively.

The rotor simulation results of pressure distributions over the five blade sections of 25%, 35%, 60%, 82%, and 92% span are shown in Figure 20. Measured data from both the MEXICO and the New MEXICO experiments are shown for comparison. For the default ${\beta }_{\infty }^{*}$ value of 0.09, notable deviations between the computational and experimental data appear on the suction side at sections of r/R = 0.35, 0.82, and 0.92. The present simulation error of the three sections is at an average level of the CFD results from the Mexnext project [29]. The Figure 20c–e give the results of r/R = 0.6, 0.82, and 0.92 for the calibrated values of ${\beta }_{\infty }^{*}$, in which the results for ${\beta }_{\infty }^{*}$ = 0.09 are also displayed. At r/R = 0.6, the prediction of = 0.09 is consistent with the New MEXICO data, while the result of ${\beta }_{\infty }^{*}$ = 0.10 is closer to the MEXICO data. At r/R = 0.82 and 0.92, the predicted pressure distributions on the suction side become much closer to the experimental data when ${\beta }_{\infty }^{*}$ = 0.105 is applied. Also, the present results are better than a collection of the CFD results from the Mexnext project [29]. The comparisons in the Figure 20d,e give strong evidence for the contribution of the present calibration to the accurate simulation of wind turbine stall.

6. Conclusions

An approach of turbulence coefficient calibration has been proposed for accurate simulation of wind turbine stall. The similarity of error between the rotor and airfoil simulations is investigated on the S809 airfoil and the NREL Phase VI rotor. It is an important foundation of the present work, which provides a way to improve rotor simulation using airfoil experimental data. The simulation of the NREL Phase VI rotor is significantly improved by making the calibration on the S809 airfoil. Furthermore, the procedure of calibration is applied to the MEXICO rotor, and it also significantly improves the simulation. The following conclusions are drawn from the present study:

• RANS simulations using the SST turbulence model tend to over-predict the aerodynamic force of airfoils and rotors in light stall. The present study involves four airfoils of S809, DU91-W2-250, RISØ-A1-21, and NACA 64–418, three of which face the problem of over-prediction at given Reynolds numbers. The over-prediction of the NREL Phase VI rotor and the MEXICO rotor is dependent on the airfoils of their blades.
• There is a high degree of similarity of error between the airfoil and the rotor simulations. In conditions of stall, the rotor simulation error is approximately at the same level of the airfoil simulation error. That means the rotor simulation error is not mainly caused by the rotational effect or the 3-D effect. For example, the range and amplitude of the simulation error at the mid-span section (r/R = 0.47) of the NREL Phase VI blade overlaps much with those of the S809 airfoil, and the peak error at the root section (r/R = 0.3) and the tip section (r/R = 95) is not higher than that of the airfoil.
• The closure coefficient ${\beta }_{\infty }^{*}$ of the SST turbulence model has a great impact on the prediction of light stall for both airfoils and rotors. The predicted airfoil lift in light stall reduces with the increase of ${\beta }_{\infty }^{*}$. For example, the lift coefficient of the S809 airfoil, which was over-predicted about 30% in light stall conditions, is reduced to a level consistent with experimental data as ${\beta }_{\infty }^{*}$ is adjusted to 0.11. Also, the error of predicted torque of the NREL Phase VI rotor is reduced when the same adjustment is made.
• The optimal ${\beta }_{\infty }^{*}$ obtained from the calibration in airfoil simulation can significantly improve the stall prediction in rotor simulation, and thus the rotor aerodynamic performance is predicted more accurately. For the NREL Phase VI rotor, the relative error of the predicted maximum torque is reduced from 34.4% to 3.2%. For the MEXICO rotor, the predicted pressure distributions over typical sections are better than the CFD results from the Mexnext project.

In essence, the present study provides an approach for calibrating the rotor simulation by using airfoil experimental data and enhances the potential of RANS in accurate simulations of the wind turbine aerodynamic performance. In future work, more simulations of various airfoils and rotors will be carried out to obtain a more comprehensive evaluation of the approach.