Prediction of Roughness Effects on Wind Turbine Aerodynamics

Wind turbines may suffer power loss after a long operation period due to degradation of the surface structure around the leading edge, by accumulation of contaminants and/or by erosion. For understanding the underlying physics, and validating the available prediction methods, the aerodynamics of smooth and artificially roughened airfoil profiles are investigated by different modelling procedures. A spectrum of turbulence models are applied encompassing RANS, URANS and DES. Two approaches are used to model the roughness. In one approach, the roughness is not geometrically resolved but its effect is modelled via the wall-functions of the turbulence models. In the alternative approach, the roughness structures are geometrically resolved, which necessitates a three-dimensional formulation, whereas the wallfunctions based approach may also be applied in two-dimensions. Computational results are compared with the experimental results of other authors, and the predictive capability of different modelling procedures are assessed.


Introduction
Wind energy is being increasingly used an in power generation from the renewable energy resources. Accumulation of contaminants and erosion can occur in a zone near the leading edge of the turbine, after a rather long operation period, which, then, causes a degradation of the performance. This can be caused by different environmental parameters, such as dust, rain and hail. Beyond the existence of contaminant and/or erosive material in the atmosphere, the degree of contamination and/or erosion is dependent on temperature, humidity, and, of course, on the wind speed.
Contamination and erosion effects can cause substantial amounts of performance loss. It was reported by Corten and Veldkamp [1] that the accumulated insect debris around the leading edge of wind turbine blades was the main cause for approximately 25% power loss observed on wind farms in California. Khalfallah and Koliub [2] analyzed the effect of dust accumulation at the leading edge of a 300kW stall regulated Horizontal Axis Wind Turbine in a dusty region in Egypt and found that a power loss about 50% in nine months can occur. Performance degradation of wind turbines can also be caused by ice accretion [3], due to similar reasons, which is also a serious problem for aircraft flight.
Surface roughness has a rather complex influence on airfoil aerodynamics, since it can also effect the transition behavior. For a smooth surface, the boundary layer transition usually occurs with involvement of the Tollmein-Schlichting waves as the primary instabilities [4]. It is known that surface roughness alters the transition behavior due to the changes in the underlying mechanism and can also trigger an earlier transition, the modes of which being strongly dependent on the type of roughness. Due to its significance, the effect of surface roughness has been investigated experimentally and computationally over several decades.
Experimental investigations on the performance degradation of airfoils with leading edge roughness was provided by Zhang et al. [5]. Further experimental studies were presented by of Ehrmann et al. [6]. Recent experimental investigations were performed by Sareen et al. [7]. It was estimated [7] that an 80% increase in drag, which was caused by a rather small degree of leading edge erosion, can result in approx. 5% loss in annual energy production, and for an increase in drag of 400-500% coupled with the loss in lift, as observed for many of the moderate-to-heavy erosion cases, the loss in annual energy production could be as high as 25%.
A computational analysis of the effect of surface roughness on dynamic stall of a wind turbine blade was presented by Noura et al. [8], where the Shear Stress Transport (SST) model of Menter [9,10] was used as turbulence model, and the roughness was modelled via corresponding wall-functions, assuming an equivalent sand-grain roughness, assumed to hold for the whole blade surface. Effect of roughness near the trailing edge of airfoils was computationally and experimentally investigated by Dhiliban et al. [11]. The standard k-ε turbulence model [9,12] was used. The quite large, regular, 2D, sawtooth like roughness elements were resolved, without applying a roughness model. CFD calculations for the effect of surface roughness on the aerodynamic performance of a turbine blade cascade was presented by Bai et al. [13], where the 4-equation transition model, i.e. the SST-γ-Reθ model (γ: intermittency, Reθ: transition momentum thickness) of Menter [9,14,15] was used. Roughness was modelled via wall-functions assuming equivalent sand-grain roughness, and applied uniformly over the whole blade surface. The SST turbulence model was applied by Zhang et al. [16] to study the effect of roughness on a blunt trailing-edge airfoil in two-dimensions, assuming sawtooth like roughness shapes, resolved by the grid. Schramm et al. [17] computationally studied erosion effects in 2D using the 3-equation SST-γ model of Menter [9]. The roughness was modelled via wall functions for low erosion while for high erosion, the blade shape was modified. Han et al. [18] applied the SST-γ-Reθ turbulence model to calculate the roughness effects in 2D, modelling the roughness by wall functions.
In the present work, we computationally investigate the leading edge roughness effects on airfoil aerodynamics applying different modelling approaches and compare the results with the measurements of Zhang et al. [5]. The main purpose is a comprehensive and coherent validation of turbulence and roughness modelling approaches. An important distinguishing feature of the present work from the previous work is the application of three-dimensional surface roughness resolving approach, which was not investigated within this context before, along with roughness modelling via wall functions and assessing the performance of the latter by comparisons. Furthermore, a wide range of turbulence models are applied with and without transition modelling, and their performance is assessed.

The test case
Experiments of Zhang [5] are considered, where wind tunnel measurements were performed for the NACA GA(W)-1 airfoil. Different Reynolds numbers (Re) and angles of attack (AoA) were considered. In the present study, we investigate the case with Re=169,000 and AoA=12°. Leading edge roughness was simulated by thin plastic strips with hemispherical roughness elements, which were placed in a region 5% chord length (c) from the leading edge in in-line or staggered arrangements, where the staggered configuration is considered in the present work. For the roughness height (k), two values were considered: k/c=0.0025, k/c=0.005.

Modelling
Incompressible flow of air described by Navier-Stokes equations [19] is computationally modelled within the framework of the general-purpose CFD software ANSYS Fluent 18.0 [9]. No heat transfer effects [20] and constant material properties are considered (ρ: density). Turbulence is modelled by different approaches. Although different versions of the k-ε model has been successfully applied in various turbomachinery applications [21], the SST model [10] is used in the present work, due to its favorable properties in modelling the near-wall flow, which has also been observed also in different applications [22]. In applying SST, a production limiter is applied along with the Kato-Launder model [9] to prevent excessive turbulence energy generation near the stagnation point, which is typical for general two-equation models. Curvature correction [9] is always applied, although it observed, in several comparisons, not to have a significant effect. The effect of the low Reynolds number (LowRe) corrections within the framework of the SST model was also investigated. Additionally, the 4-Equation SST-γ-Reθ transitional turbulence model [9], as well as the 3-Equation SST-γ transitional model [9] are applied. Besides the surface resolving approach, capturing the real roughness geometry, wall-functions based roughness modelling approach [9] is also used, using different methods [23] for estimating the equivalent sand grain roughness.
For the velocity-pressure coupling, the SIMPLEC [9] scheme is used. In unsteady calculations, a second-order accurate bounded backward differencing scheme is used for time discretization [9], choosing time step-sizes to assure cell Courant numbers smaller than unity. For the discretization of convective terms, the second order accurate upwind scheme [9], for all variables with the exception of momentum equations within the framework of DES, where a bounded central differencing is used.

Solution domain and boundary conditions
A 2D view of the solution domain is presented in Figure  1, with indication of boundary types. At the inlet, uniform profiles are applied for the velocity components and turbulence quantities of the uni-directional flow. (U: inlet velocity) Assuming a low turbulence level at the inlet, a turbulence intensity of 0.5% and turbulent-tolaminar viscosity ratio of 1 is assumed to derive the inlet boundary conditions of the turbulence quantities. For the intermittency, γ=1.0 is applied [9]. At walls (the airfoil), the no-slip conditions apply, whereas at the outlet the static pressure is prescribed along with the zero-gradient conditions for all diffusive-convectively transported variables. At the upper and lower boundaries (Fig. 1), symmetry planes are assumed. In 3D, the 2D domain is extruded for a distance of 0.04c. Boundaries occurring in the third dimension, are prescribed to be periodic.

Grid generation
In problems with roughness, grid generation is very challenging due to the huge scale disparity, if the geometries of roughness structures are intended to be resolved by the computational grid, as it is presently the case. In the present work, a grid generation strategy is applied, consisting of three steps. As much as possible, the grids are generated applying structured/block-structured, arranged as C-type around the airfoil (Fig. 1). In all grids, it is ensured that the nondimensional wall distance y + is always smaller than unity (roughened regions while using a correlation based roughness modelling is an exception, where a "y + -shift" is applied by the model [9]) The first step has been the generation of an adequate grid for two-dimensional flow.

Results of the grid independence study (RANS-SST)
showing the variation of the predicted lift coefficient (CL) by different grids relative to the finest grid are presented in Table 1. One can see that the variations are small, and even the coarsest grid (114,000 nodes) can be seen to provide sufficient grid independence.
For cases, where the surface roughness is to be resolved by the grid, a three-dimensional formulation is necessary. A further reason for a three-dimensional formulation is given through the unsteady formulation, since three-dimensional flow structures can temporarily emerge, although the time-averaged flow may be twodimensional.
Keeping the structured grid topology while retaining the resolution on the blade surface including the third dimension would result in too large grid sizes. For this purpose, as the second step of the grid generation strategy, a locally refined block-structured grid with nonconformal block interfaces is generated, which provides a comparable accuracy to the structured grid from the first step. The finest structured grid from the first step (G5, Table 1) and the resulting non-conformal blockstructured grid (G-blocks) are displayed in Fig. 2.
The grid sizes of the grids and the deviation in the predicted lift coefficient are presented in Table 2. The results of the G-blocks can be considered to be sufficiently close to those of G5, indicating the adequacy of the former.
In generating the three-dimensional grid, where the 2D geometry is extruded by a distance of 0.04c in the third direction, the G-blocks grid is taken as basis and extruded. For the grid resolution in the third dimension, An exception to the previously described structure is the treatment of a tiny region at the leading edge, where the roughness elements are placed, for the threedimensional roughness-resolving calculations. This tiny region is discretized by an imbedded unstructured block which corresponds to approx. 300 subdivisions in the third direction on the blade surface. Based on this "base grid", local grid refinements are additionally applied during the calculations, to ensure the local near-wall y + values to remain smaller than unity. A section of the 3D grid in is presented in Fig. 3, where this region can be recognized (in comparison to Fig. 2). The number of nodes of the resulting three-dimensional grids are between 8·10 6 and 9·10 6 , depending in the roughness type considered. A detail view of the surface grid generated (before local grid refinements) for the threedimensional surface roughness-resolving approach for the case k/c=0.0025 with staggered arrangement is provided in Fig. 4. For the resolution of the third direction, no grid independence study is performed. It is assumed that the surface resolution near the leading edge is sufficient, and the need for an eventual finer resolution in downstream is assessed by inspecting the resulting distributions.

Overview of applied modelling approaches
A large number of modelling approaches are investigated. Table 3

Field distributions
A general view of the flow pattern is provided in Figure  5, showing velocity magnitude (u) predicted by SST-S.   Near the leading edge, the stagnation region on the pressure side, as well as the acceleration over the nose, towards the suction side can be observed (Fig. 5).
For the smooth blade, the predicted streamlines by different models are compared with the measured ones in Figure 6. One can see that the measurements indicate a separation bubble on the suction side near the leading edge (Fig. 6a). As the SST model (SST-S, Fig. 6b) cannot predict this, the transitional SST models (SST-S-3TR, Fig. 6d, SST-S-4TR, Fig. 6e), but also the SST model with LowRe correction (SST-S-LR, Fig. 6c) can qualitatively predict the phenomenon.
For the rough blade, distribution of the wall shear stress (τW) near the leading edge predicted by DES is  Figure 7, for an instant of time, for two roughnesses. It can be seen that shear stress is distributed quite unevenly over the roughness elements.

Lift, Drag, Pressure Coefficients
For the smooth blade, for Re=169,000 and AoA=12°, the variations of the pressure coefficient (Cp) along the blade surface are predicted by the applied twodimensional modelling approaches are compared with the measured values (EXP) in Figure 8. One can see that all presented predictions show a generally quite fair agreement with the measurements. The variations between the models are observed rather on the suction side of the blade, where SST-3TR and SST-4TR show a better agreement with the experiments.
Predicted coefficients of lift (CL) and drag (CD) for the smooth blade are compared with experimental values [5] (EXP) for Re=169,000 and AoA=12°, in Table 4. Percentage deviations of predictions from the measured value are also indicated in parentheses. One can see that CL is over-predicted by about 18% by the SST. A very good prediction of CL is provided by SST-3TR, which is even better than that of SST-4TR. The 3D, unsteady IDDES modelling provides a quite good CL prediction, which is, but, not as good as that of SST-3TR.

Conclusions
Aerodynamics of smooth and at the leading edge roughened wind turbine blades are calculated by different turbulence modelling approaches. For the smooth blade, SST based transitional turbulence models applied in 2D, or in 3D within a IDDES formulation are observed to provide quite good predictions of the lift coefficient. For the rough blade, the 3D, IDDES based surface resolving approach is observed to perform remarkably better than the empirical roughness modelling approaches in 2D.