1 Introduction

In recent years, wind energy, as a form of environmentally friendly and sustainable energy, has been used worldwide, and its applications are expected to increase in the future, especially in Northwest China. However, sand-wind environments can significantly affect the aerodynamic characteristics of wind turbines. With the passage of time, wind turbine blades will be eroded by the dust in air, which will decrease the power generation capacity and service life of the wind turbines. Therefore, it is necessary to study the effects of erosion on wind turbine airfoils.

Many researchers have investigated this issue over the past 20 years. Khakpour et al.^[1] pointed out that the aerodynamic performance of wind turbine airfoils changed significantly in different dusty environments. Huang et al.^[2] found out that the surface roughness greatly affected the aerodynamic performance of an airfoil. When the surface roughness was high, the lift coefficient (C_l) would decrease by about 40% and the drag coefficient (C_d) would increase by 10 times, compared with a smooth one. They showed that the surface around the leading edge was the most prone to erosion. Sareen et al.^[3-4] and Sayer et al.^[5] studied the effects of leading edge erosion on wind turbine airfoils.

Khalfallah and Koliub^[6] measured the process of erosion in Hurghada wind farm, and pointed out that the blade damage leaded to a significant decrease in the electrical power generated by a wind turbine, and this decrease would emerge over several months due to the effect of a severe sand-wind environment. Alden and Diab^[7] used the discrete phase modeling capability in ANSYS-FLUENT along with the Det Norske Veritas (DNV) erosion model to predict the erosion rates of wind turbine blades, used the erosion rates to calculate the erosion depth, and obtained a result similar to that of Khalfallah and Koliub^[6]. Ren and Ou^[8] arrived at a similar result by using a numerical method. Tabakoff et al.^[9-11], Tabakoff and Metwally^[12], Tabakoff^[13], and Hussein and Tabakoff^[14-15] studied the particle trajectory and the erosion mechanism of the coating used on the surface of wind turbine blades. Grant and Tabakoff^[16] proposed a method to predict the physical process of turbomachinery erosion by means of Monte Carlo simulation. Fiore and Selig^[17-20], Fiore^[21], and Fiore et al.^[22] focused on the erosion mechanism for wind turbine blades with a simulation method. Diab et al.^[23] and Alden and Diab^[24] studied the effects of dust contamination on the annual energy production loss in sand-wind environments.

Li et al.^[25] proposed that there existed a critical Stokes number for wind turbine airfoil erosion, and once the particle Stokes number was larger than the critical value, the airfoil would be eroded. Notably, the angle of attack (AOA), thickness, and airfoil series have little effect on this critical value. It is innovative to connect airfoil erosion with the particle Stokes number, which is widely used in gas-solid flow and can be used to characterize particle flows. Besides, the specific relationship between airfoil erosion and particle Stokes number remains unknown. Therefore, in this paper, we aim to establish this relationship.

2 Mathematical model

The Euler-Lagrange model is adopted in this work. The air is treated as a continuous phase by solving the time-averaged Navier-Stokes equations in the Eulerian reference frame, while inert particles (sand/dust particles, discrete phase) are modeled in a Lagrangian frame of reference to track the particle trajectories.

2.1 Continuous phase model

The flow around an airfoil is modeled with the Navier-Stokes (N-S) equation for three-dimensional (3D) viscous incompressible flow, and the continuity and momentum equations based on the Reynolds averaged N-S equations are solved by using the shear stress transport (SST) k-ω turbulence model. Though this turbulence model has been widely used in numerical simulation, it has shown in experience that this model cannot predict the aerodynamic performance of airfoils very well under stall conditions. Therefore, Wen et al.^[26] and Zhong and Wang^[27] proposed to mitigate these deficiencies by revising the constants a₁ and β^* to improve the capability of the model so as to predict the stall characteristics. According to Wen et al.^[26], a₁ and β^* can be revised by varying the AOA (see Fig. 1).

2.2 Discrete phase model

The discrete phase model (DPM) traces the particles in the discrete phase, and predicts the trajectory of each particle based on the balance of the forces involved. A surface injection method with 960×80 points, where the particles are released, is adopted (see Fig. 2). The initial speed of particles is set to be equal to the wind speed in the present study.

2.3 Erosion model

The erosion rate ε, which is defined as a dimensionless ratio of the mass of the removed eroded material to the erodent particles (sand), can be expressed as follows:

where m_p represents the mass of the impacting particles, c(d) is a function of the particle diameter, f(α) is a function of the impact angle, and α is the impact angle particle path with the target area A_f. Currently, mostly polyurethane derivatives are used for coating the blade of wind turbines^[28], and such materials primarily undergo plastic erosion behavior with the maximum erosion rate at α = 30º^[29]. v is relative particle velocity, and b(v) is a function of the relative particle velocity. Notably, all these functions depend on the eroded material and the material of the erodent particles. In this study, c=1.8×10^-9, b=2.6, and the values of α and f(α) are listed in Table 1.

2.4 Computation domain

A C-shaped computational domain is formed by use of a semicircle and a rectangle, where the spanwise length is 0.5 m (see Fig. 2). The chord length (C) of each airfoil used in this study is 1 m. The airfoils used in Subsections 3.1 and 3.2 are S809, and the airfoils used in Subsection 3.3 will be introduced at the beginning of the section. A structural mesh is adopted in this study. The height of the first grid point close to the airfoil surface is set to be 1×10^-5 m to ensure y⁺≤ 1 for Re = 1×10⁶.

2.5 Model validation

The numerical method for clear inflow is validated by comparing the simulation results with the experimental data from the National Renewable Energy Laboratory (NREL) in Ref. [30], and good agreement is obtained (see Fig. 3).

For DPM validation, the experiment about the effect of the gas-solid flow on the aerodynamic performance of the airfoil is basically none. Therefore, an experiment about the circular cylinder flow in a gas-solid environment proposed by Luo et al.^[31] is used to validate the DPM model. As shown in Figs. 4 and 5, the numerical results agree well with the results of the particle image velocitmetry (PIV) test, where d_c in refers to the cylinder diameter.

3 Results and discussion

The particle Stokes number S_tk is the ratio of the particle relaxation time to the flow characteristic time, and it is expressed by

where τ_p is the particle response time, and

d and ρ are the diameter and the density of the particles, respectively. μ is the dynamic viscosity of the air. T_f is the fluid response time, and

L is the characteristic length of the flow. u_∞ is the characteristic speed of the flow. In this study, L and u_∞ denote the chord length of the airfoil and the inflow velocity (V), respectively. The smaller the value of S_tk is, the shorter the time required by the particles to respond to the changes in the flow is.

3.1 Validation of the critical particle Stokes number range

Unless stated otherwise, the flow velocity is 14.61 m/s, the Reynolds number Re = 1×10⁶, the particles density is 920 kg/m³, and the AOA is 6.1º in this study. The particle mass flow rate accounts for 1% of the air mass flow rate.

To ascertain the critical particle Stokes number S_Ctk for the airfoil erosion when erosion is initiated on the airfoil surface, we start by seeking the critical particle diameter for the initiation of erosion. As shown in Fig. 6, when d=10 μm, erosion does not occur, while when d=20 μm, erosion occurs. Therefore, S_Ctk is initially determined to range from 0.004 2 to 0.017 0. Then, dichotomy is used to narrow the range of S_Ctk. As shown in Fig. 6, erosion does not occur when d=15 μm or 18 μm, but occurs when d=19 μm. Then, many values of d are tested, and it is determined with certainty that erosion will occur when d≥19 μm, and it will not occur if d≤ 18 μm. Therefore, the critical particle diameter range should be from 18 μm to 19 μm, and the corresponding S_Ctk range is from 0.013 5 to 0.015 1. It is certain that a value of S_Ctk exists in the range from 18 μm to 19 μm. If the value of S_tk is larger than this S_Ctk, airfoil erosion will occur.

The range of S_Ctk obtained herein is close to the range obtained in Ref. [25]. However, the chord length of the airfoil in Ref. [25] is 0.283 m, while the chord length in this study is 1 m. It is inferred that the range of S_Ctk for S809 is from 0.013 5 to 0.015 1 at an AOA of 6.1º.

The range of S_Ctk can be used as a criterion for the estimation whether an airfoil will undergo erosion (see Ref. [25] for the detail). It is confirmed that the range of S_Ctk will be not affected by the AOA, the relative thickness of the airfoil, etc. However, the particle Stokes number should be determined based on V, d, and ρ. Therefore, whether S_Ctk can be used to determine the erosion state of an airfoil requires further verification.

The applicability of the range of S_Ctk is demonstrated by the variations of V, d, and ρ. The values of the controlling parameters in each case are listed in Table 2. The results are shown in Fig. 7.

As shown in Figs. 6 and 7, erosion appears in Cases 2, 3, and 5, while no erosion appears in Cases 1, 4, and 6. In all cases, erosion occurs only when the value of S_tk of the particles is higher than the maximum value in the range of S_Ctk. Therefore, the range of S_Ctk is applicable to any dusty environment. In summary, the range of S_Ctk can be regarded as a criterion to judge whether a wind turbine airfoil will be eroded by comparing the value of S_tk with the range of S_Ctk.

3.2 Erosion extent of airfoils at the same particle Stokes number

Several operating conditions at a given S_tk are set up by changing the parameters, e.g., the inflow velocity, the particle diameter, and the particle density. The erosion extent under different operating conditions is compared so as to establish a relationship between the particle Stokes number and the airfoil erosion.

Five sets of controlling parameters at the same particle Stokes number of 0.15 for the AOA of 6.1º are listed in Table 3.

Erosion occurs mainly near the leading edge, and it usually spreads from the leading edge to the trailing edge along the pressure surface. Figure 8 shows that at the AOA of 6.1º, the erosion areas of the airfoil are essentially the same, and the maximum erosion rates are almost identical when the values of S_tk are the same even if the controlling parameters are different. To compare the size of the erosion area more intuitively, the ratio (A) of the erosion area to the area of the pressure surface and the maximum erosion rate (ε_max) in each of the cases are listed in Table 4. From the table, we can see that, the erosion extent of the airfoil is almost the same when the particle Stokes number is 0.15 in each of the cases.

To further verify the effect of the Stokes number on the airfoil erosion, Table 5 lists five sets of controlling parameters, where the Stokes number is 0.267, and the AOA is 11.21º. From the table, we can see that no erosion occurs on the suction surface in these five cases. The erosion areas are enlarged compared with those in Table 4. The erosion areas and the maximum erosion rates of the five cases are almost the same. This confirms that the above conclusion is applicable when the AOA is 11.21º (see Fig. 9 and Table 6).

When the AOA changes to 16.2º, the S809 airfoil will stall. The controlling parameters are the same as those in Table 5. The results are shown in Fig. 10 and listed in Table 7. From the results, we can see that the above conclusion is applicable at this stall AOA as well. A comparison among the cases with the AOA of 11.21º shows that the erosion area is enlarged slightly, and the maximum erosion rate increases slightly. This is because that when the AOA increases, the impact area increases considerably, and more particles crash onto the airfoil surface.

The results show that the erosion area on the airfoil and the maximum erosion rate are almost equal under the same particle Stokes number and AOA. This shows that it is feasible to predict the erosion extent on an airfoil surface based on the value of S_tk. According to the measurements and calculations, the particle Stokes number can be obtained and used to estimate in advance whether a wind turbine blade will undergo erosion, and the erosion extent can be predicted.

3.3 Effects of the airfoil shape on erosion

If the AOA changes, the flow around the airfoil will definitely change, which will greatly affect the flow of the particles, and the erosion extent of the airfoil will also change. The range of S_Ctk for the S809 airfoil is ascertained previously in dusty environments at the AOA of 6.1º. In this subsection, the same simulation will be repeated at the AOA of 18.2º. Figure 11 shows the particle tracks when erosion is initiated at the AOAs of 6.1º and 18.2º. As shown in the figure, at 18.2º, there is no collision between the particles and the airfoil surface, even though the particle diameter is 24 μm. Therefore, erosion is not initiated. When the particle diameter is changed to 25 μm, a few collisions could be observed between the particles and the airfoil surface, which will lead to erosion. The corresponding range of S_Ctk is from 0.024 to 0.026. This indicates that the airfoil is prone to erosion at small AOAs, and it is not eroded easily at large AOAs. Figure 11 shows that the particles can bypass the airfoil surface increasingly easily when the AOA increases, which could prevent collisions. Thus, the value of S_Ctk increases when the AOA increases.

The above finding seems to be inconsistent with the conclusions obtained in Ref. [25], where the range of S_Ctk did not change with the AOA. It is speculated that the geometric shape greatly affects the erosion characteristics of airfoils, because the simulation object in Ref. [25] is NACA 4418, the geometric shape of which is considerably different from that of the S809 airfoil.

Figure 12 shows the profiles of the airfoils used in this section. To investigate the effects of the geometric parameters of the airfoil on the erosion characteristics, the radius of the leading edge (R), the maximum thickness (T), and its position with respect to the chord (x(T)/C) are considered.

Figure 13 shows the particles tracks when erosion is initiated on the surface of the NACA 4418 airfoil. The range of the critical particle diameter for the initial erosion on the airfoil is from 15 μm to 16 μm for the AOA of 6.1º and from 18 μm to 19 μm for the AOA of 18.2º, and the corresponding ranges of S_Ctk are from 0.009 4 to 0.011 0 and from 0.013 5 to 0.015 1, respectively. These two ranges of S_Ctk essentially coincide with the range of S_Ctk, i.e., from 0.007 8 to 0.014 0 in Ref. [25]. This shows that the NACA 4418 and the S809 airfoils exhibit different erosion characteristics. It is found that, for S809, T is larger and R is smaller, i.e., T = 21%, and R≈1%, and thus the shape of the pressure surface near the leading edge is flat. For NACA 4418, the opposite is true, where

Figures 11 and 13 show that the particles crash more easily onto the airfoil surface and cause erosion in the case of NACA 4418, owing to the larger curvature of its pressure surface near the leading edge. On the contrary, owing to its smaller curvature, particles easily bypass the airfoil surface in the case for the S809 airfoil.

Compared with NACA 4418, the GOE798 airfoil is characterized by larger R and T, i.e,

but the shapes of the pressure surfaces of the two airfoils essentially coincide for x/C≥0.09.

Figure 14 shows the particle tracks when erosion is initiated on the GOE798 airfoil surface. From the figure, we can see that the values of S_Ctk for the GOE798 airfoil range from 0.008 2 to 0.009 4 at the AOA of 6.1º and from 0.012 0 to 0.013 5 at the AOA of 18.2º. Comparing NACA 4418 and S809, it can be concluded that airfoil is more prone to erosion when R is larger.

It is indicated that airfoil erosion is more likely to occur when the AOA is small. By contrast, the effect of the AOA on airfoil erosion is related to the airfoil shape. The smaller the curvature (close to the leading edge, the pressure surface) is, the more sensitive the value of S_Ctk is to the changes in the AOA.

The S802 (T = 12%, R≈1%), S822 (T = 16%, R ≈1%), and S809 airfoils (T = 21%, R≈1%) are selected to investigate the effects of T on airfoil erosion. These three airfoils have similar x(T)/C, i.e., 40%. The values of S_Ctk of the S802 and S822 airfoils at the AOA of 6.1º range from 0.011 0 to 0.012 1 and from 0.012 1 to 0.013 5, respectively, and they range from 0.024 0 to 0.026 0 at the AOA of 18.2º. A comparison of the particle trajectories shown in Figs. 11 and 15 shows that, although the S_Ctk ranges of the three airfoils coincide at the AOA of 18.2º, the erosion of the S802 airfoil is the most severe, and the erosion of the S809 airfoil is the least severe. This is because of the differences in T. With a similar R, when T increases, the pressure surface (near the leading edge) will be flatter and more slanted, which ensures that the particles could bypass the airfoil surface more easily.

When x(T)/C changes, the impact angle will change. This will greatly affect the erosion characteristics of the airfoil. The S811 (T = 26% and x(T)/C = 30%) and S814 (T = 24% and x (T)/C = 30%) airfoils are selected for the investigation of this effect of x(T)/C. The S_Ctk ranges of the S811 and S814 airfoils at the AOA of 6.1º are from 0.024 0 to 0.026 1 and from 0.022 1 to 0.024 0, respectively, and are from 0.037 6 to 0.040 0 and from 0.035 1 to 0.037 6 at the AOA of 18.2º. In comparison with the corresponding values of the S809 airfoil, a smaller x(T)/C causes the pressure surface to be more slanted (near the leading edge). An airfoil characterized by such a shape might show a strong decrease in the particle velocity, which means that, airfoils, such as S811 and S814 used at the blade root, are difficult to erode.

From the above analysis, we can conclude that, the erosion characteristics of a wind turbine airfoil depends on the airfoil geometry. The curvature of the pressure surface near the leading edge is the key to improving the anti-erosion performance, and it is decided by many geometric parameters, e.g., R, T, and x(T)/C. The goal is to find a balance point between the aerodynamic characteristics and the anti-erosion characteristics of the airfoils. The anti-erosion characteristics can be improved to the extent possibly by modifying the airfoil shape while guaranteeing the aerodynamic characteristics. Then, the service life of wind turbines can be increased, and their power generation capacity can be increased and stabilized. This work can serve as a reference for wind turbine blade design.

4 Conclusions

At the AOA of 6.1º, S809 airfoil is found to undergo erosion when the value of S_tk is greater than 0.015 1, but no erosion occurrs when the value of S_tk is smaller than 0.013 5. This illustrates that there is a critical range for the particle Stokes number for airfoil erosion, and the erosion of the airfoil surface occurs only when the value of S_tk exceeds the maximum value in this range.

The erosion area on the airfoil and the maximum erosion rate are almost the same under the same particle Stokes number and the AOA. This indicates that S_tk can be used to predict the erosion extent of of an airfoil. Meanwhile, erosion will be more severe at larger values of S_tk.

The geometric shape of the pressure surface greatly affects the airfoil erosion, especially in the part near the leading edge. The larger the curvature is, the more likely the airfoil is to undergo erosion.

Airfoil erosion is difficult when the AOA is large. However, the effect of the AOA on airfoil erosion is related to the airfoil shape.