Analysis of leading edge erosion effects on turbulent flow over airfoils

The surface of wind turbine blades are prone to degradation due to exposure to the elements. Rain, hail, insects are among the many causes of turbine blade degradation or erosion. Surface degradation of the wind turbine blades leads to a reduction in the aerodynamic performance, resulting in power losses. The effect of surface degradation is studied by modeling the turbine blade as a rough surface. Surface roughness can be positive (insects or other foreign objects) or negative (erosion, delamination). The individual roughness elements are however very small and it is not always feasible to study the actual degraded surface. Thus various roughness models have been proposed in the literature which eliminate the need to accurately model the degraded surface by representing erosion with a virtual surface and modeling the effect of erosion on the flow quantities near the eroded surface. In this study, the reduction in performance of airfoils due to leading edge roughness is quantified. Different roughness models are investigated and evaluated against theoretical models. Additionally, the effect of roughness on different integral boundary layer quantities like displacement thickness, momentum thickness and skin friction are presented. © 2021 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Leading edge erosion is an issue of growing concern in the wind turbine industry in recent years. The combination of growth in the size of wind turbines, increased offshore installations, especially in locations with more adverse weather conditions, has made this subject crucial to the industry [1]. Erosion of turbine blades are largely caused by rain, hailstones, accumulation of contaminants and tends to change the shape of the airfoils. This leads to a reduction in aerodynamic performance of the affected sections. Han [2] presented the effects of contamination of the airfoil used at blade tips on a 5 MW NREL turbine blade using CFD simulations. They report a worst case scenario where the Annual Energy Production (AEP) drops by 3:7%. Herring [1] presents a thorough review on the growing importance of leading edge erosion and different coating alternatives to reduce the impact of erosion. A wide range of drop in AEP, from about 25% to about 3:7%, is reported and the authors suggest it is due to different operating conditions and roughness levels used to evaluate the impact of erosion. The authors also note that repair of moderate erosion can recover the AEP by about 2%. Keegan et al. [3] review some of the leading causes of erosion of wind turbine blades including raindrops, hailstones, exposure of composite materials to UV surfaces among others. The authors note that the increasing tip speeds of turbine blades make them increasingly vulnerable to erosion by raindrops and hailstones.
In order to quantify the adverse effects of roughness the flow around the turbine blades should be investigated. Laminar flow tends to transition to turbulent flow prematurely in presence of roughness. A review of experimental approaches to model roughness and its effect on transition can be found in Ehrmann et al. [4] Langel et al. [5] performed experiments on two airfoils by adding cut vinyl decals and focused on 100 < Re k < 400, where Re k is the Reynolds number based on roughness height k. They also present a numerical approach to model the effect of roughness on transition by adding a scalar field variable. The new scalar variable is used to modify the g À Re q transition model [6]. Sareen et al. [7] note that most of the experimental studies on roughness use strips or zigzag tapes to simulate real roughness and not many studies exist on negative roughness like erosion where material is lost from the blade.
Apart from causing early transition, the nature of the turbulent boundary layer also changes due to roughness. Skin friction increases and a shift in the velocity profile in the inner part of the boundary layer is observed. The additional dissipation near the roughness elements leads to thickening of the boundary layer which can make the boundary layer prone to early separation.
The concept of equivalent sand grain roughness is widely used in turbulence models to account for the effect of roughness on turbulent boundary layers. Nikuradse [8] performed experiments to measure pressure losses across pipes due to roughness, which forms the basis of the sand grain roughness concept. Nikuradse provided relations for the loss in pressure head (friction) and the velocity shift as a function of sand grain roughness heights. Real roughness is first converted to equivalent sand grain roughness when using the roughness models for Reynolds averaged Navier Stokes (RANS) turbulence models. Typically the rough surface is replaced by a smooth surface and the effect of roughness is modeled as extra dissipation in the inner boundary layer.
Integral boundary layer based tools like RFOIL [9] are used extensively in the wind energy community for quick and accurate analysis of airfoil performance, especially in combination with other methods like Blade Element Momentum theory, to obtain power output of wind turbines in a relatively inexpensive manner. However, it is restricted mainly to clean airfoils due to lack of research on developing roughness models for integral boundary layer methods. Olsen et al. [10] recently proposed a new closure relation for skin friction in the presence of roughness. The authors also note that further work is necessary to refine their study.
In this study, roughness models for SA and SST kÀ u turbulence models are implemented in the open source tool SU2 [11]. The grid requirements and the accuracy of the two models are examined and validated against experimental data. Two airfoils are considered -NACA 65 2 215 and a popular wind turbine airfoil DU-96-W-180. The NACA 65 2 215 airfoil has been used for validating roughness models earlier [12,13]. Sareen et al. [7] performed experiments on the DU-W-96-180 with 'negative' roughness. Thus different ways to obtain equivalent sand grain roughness for 'negative' roughness are also examined in this paper. The numerical solution of the RANS equations is then used to analyze the behavior of the turbulent boundary layer and the various integral boundary layer quantities in the presence of roughness as well as to analyze the integral boundary layer parameters in order to improve roughness modeling in integral boundary layer methods.
The organization of the paper is as follows: section 2 gives information about SU2, the CFD solver used for the numerical simulations. In Section 3, two different roughness models for RANS are presented with some validation cases in section 4. Based on the results in section 4, the SA roughness model is validated against experiments on airfoils in section 5. In section 6, the effect of roughness on various integral boundary layer properties is analysed. The conclusions are presented in section 7.

SU2
SU2, the CFD solver used in this study, is an open-source collection of Cþþ based software tools for performing Partial Differential Equation (PDE) analysis and solving PDE-constrained optimization problems [11]. Originally developed for aerospace applications, the solver has been extended for incompressible flows [14,15]. In this study, we use the low Mach preconditioned incompressible flow solver [14]. The governing equations of SU2 solved on a domain U are written in the general form as where U is the vector of conservative variables, F c i are the convective fluxes, F v i are the viscous fluxes and Q is a source term defined as Here u i are the components of the velocity vector, r is the density, P is the dynamic pressure and the viscous stresses are t ij ¼ The total viscosity coefficient, m tot is the sum of the dynamic viscosity m dyn and turbulent viscosity m tur , which is computed via a turbulence model. The Spalart-Allmaras (SA) [16] and the Mean Shear Stress Transport (SST) [17] turbulence models can be used to compute m tur . More details on the low Mach number preconditioning method can be found in Ref. [14].

Spatial discretization
The spatial discretization is performed on an edge based dual grid using a finite volume approach. The control volumes are constructed using a median-dual (vertex-based) scheme. An upwind Flux Difference Splitting (FDS) scheme is used to compute the convective flux residual. The MUSCL scheme in combination with the van Albada slope limiter is used to obtain second order accuracy. The gradients for flux reconstruction are computed using the Weighted Least Squares method. The gradients required to evaluate the viscous fluxes are computed using either the Least Squares method or the Green Gauss theorem.

Time discretization
Steady state problems are also solved using a pseudo-time stepping approach where the solution is marched in time until convergence. Time integration is carried out using the Implicit Euler method.

Boundary conditions
For the test cases considered in this study, only the no-slip wall boundary condition on the surface of the airfoil and the far field boundary condition on the external domain is necessary. Since SU2 uses a vertex-based dual grid approach, the implementation of the no slip boundary condition is relatively straightforward. For the momentum equations, a no slip condition is enforced strongly by setting the velocity on the wall to zero (since no wall movement is necessary) and a Neumann boundary condition is used for the other equations. For the far field boundaries, the flux across the face is computed in a similar manner to internal faces where the neighboring states are assumed to be the internal solution and the free stream value.

Spalart-allamaras (SA)
The SA model [16] with no trip term can be written in the general form of equation 3 The turbulent viscosity is then computed as Here n ¼ m r is the kinematic viscosity and d is the distance to the nearest wall. The definitions of the other model constants can be found in the literature [16,18].

SST k-u
Following the general form of the equations in equation (1), the corresponding terms for the SST k-u [17] model are Here the production term, P ¼ t ij vu vxj , where t ij is defined earlier in section 2.1, r is density, n t ¼ m t =r is the kinematic turbulent viscosity and m is dynamic viscosity. The turbulent eddy viscosity is computed as where U is the vorticity magnitude and F 2 is a model constant.
More information can be found in the literature [17].

Boundary conditions
At the far field the boundary conditions for the SA and SST kÀ u model are respectively n t;∞ ¼ 0:210438n ∞ to 1:294234n ∞ ; Here n ∞ is the kinematic viscosity in the free stream, V ∞ is the free stream velocity magnitude, r is the density and TI is the turbulent intensity. The ratio m t =m lam and turbulent intensity TI are specified as inputs. On solid walls, the boundary conditions for clean walls are defined below. n t ¼ 0; Dd is the distance to the nearest normal neighbor and b 1 is a model constant. Model constant definitions can be found in the literature [11,19].

Roughness modeling
To motivate the roughness model used in this study, a brief introduction of turbulent boundary layers and the impact of roughness is presented below.
The turbulent boundary layer can be broadly divided into two regions [20,21]: the inner region where viscous dissipation is comparable to the turbulent dissipation and the outer region where turbulence dissipation dominates completely. The inner region can be further subdivided into three regions -the viscous sub-layer where viscous effects dominate and turbulent effects are absent, a buffer region where the turbulent stresses start to grow and finally an overlap region or a logarithmic region where the turbulent and viscous dissipation match. The overlap region leads into the outer layer of the boundary layer where viscous effects are minimal. The velocity profile in the viscous sub-layer and logarithmic region can be written respectively as The region of the boundary layer between 5 y þ 30 is the buffer region. In the above relations, y þ is the non dimensional wall normal coordinate and u þ is the normalized velocity defined as Here u t is known as the wall friction velocity and is used as the velocity scale close to the wall, t w is the wall shear stress, r is the density, u is the local velocity and n is the kinematic viscosity. The constant in equation (7) for a smooth wall is known to be C ¼ 5:0. The presence of surface roughness on the wall alters the nature of the velocity distribution near the wall. The roughness elements will introduce new turbulent fluctuations in the flow increasing the skin friction. Typically, a standardized notion of roughness known as the "equivalent sand grain roughness height (k s )" is used to denote roughness of a wall [8,21,22]. A given physical roughness distribution is converted into the "equivalent sand grain roughness height" using empirical correlations like [23e25]. A more detailed review is presented in section 5.2.3. Based on the non dimensional roughness height, k þ s ¼ k s u t =n, three regimes of roughness can be identified [21]. If the roughness elements are within the viscous sub-layer (k þ s 5, hydraulically smooth), the effect of roughness is not relevant and there is no difference with the smooth velocity profile. As the height of the roughness element increases (5 k þ s 70, transitionally rough), a shift in the velocity profile is observed. Once the roughness elements are fully within the overlap region (k þ s > 70, fully rough), the viscous sub-layer plays no part and the flow is in the fully rough regime. It must be noted here that the equivalent sand grain roughness concept typically applies only to the commonly observed distributed roughness (kÀ type roughness [26]) and not to isolated roughness elements. To reproduce the proper shift Du þ in the boundary layer velocity profiles, turbulence models typically increase the eddy viscosity dissipation within the inner part of the boundary layer [22]. Aupoix et al. [22], identify two methods to accomplish this with eddy viscosity based turbulence models (e.g SA and SST): 1. Finite eddy viscosity at the wall which can be interpreted as using a virtual wall to represent roughness and 2. Zero eddy viscosity at the wall where the origin of the wall is at the bottom of roughness but turbulence damping in the wall region is reduced.
With this background on roughness modeling in turbulent boundary layers, roughness models for the SA and SST turbulence models are presented.

Roughness modification for SA model
The roughness modification proposed by Boeing [18,22] is considered in this section. An alternate modification was also proposed by ONERA in Aupoix et al. [22], but is not considered since it requires the additional input of friction velocity. The effect of roughness is accounted for by shifting the virtual wall to the top of the roughness element. This can be achieved by offsetting the distance to the wall everywhere. The changes to the turbulence model are with c R1 ¼ 0:5. The eddy viscosity at the wall is now changed from n ¼ 0 to a non-zero value by using a mixed (Robin) boundary condition at the wall, vñ vn j wall ¼ñ wall 0:03k s ; where vñ vn is the gradient ofñ in the direction normal to the wall.

Roughness modification for SST model
The effect of roughness can be accounted for in the kÀ u SST turbulence model by modifying the boundary conditions at the wall as [19].
From equation (5), the eddy viscosity remains zero at the wall, but there is an increase in turbulence dissipation compared to the clean boundary conditions. Here k rough is the turbulent kinetic energy and k þ s is the non dimensional equivalent sand grain roughness height.
The two roughness models are implemented in SU2 and are validated below.

. Grid refinement study
Turbulent flow over a flat plate with different roughness heights is simulated with the SA and the SST turbulence models and their respective roughness corrections presented above. The flat plate domain is 2m long and 1m high and Re ¼ 6:0 Â 10 6 . A grid refinement study is carried out for the geometry under clean and three roughness levels. There are 57, 113 and 225 points on the surface of the 2-D flat plate for the three grids. The minimum grid spacing is Dy 1 z2 Â 10 À6 m. A second set of grids are made with same geometry and same number of points on the surface but with a minimum grid spacing of Dy 2 z3 Â 10 À8 m for the SST roughness model. A growth ratio of 1.09 is used in the normal direction. The skin friction values computed at x ¼ 0:93m are tabulated in Table 1.
Three different roughness heights, k s ¼ 1:23 Â 10 À4 m, k s ¼ 2:46 Â 10 À4 m and k s ¼ 9:48 Â 10 À4 m are tested. The k þ s values are around 24, 50 and 200 respectively. With a grid spacing of Dy 1 the y þ z0:3 at x ¼ 0:93 under clean conditions. As seen in Table 1, a grid independent solution is obtained for the clean case for both the turbulence models with this minimum grid spacing. The SA roughness model gives largely grid independent result for all the roughness heights. However, a grid independent solution is not possible under rough conditions with the SST model. The variation is marginal at low roughness heights and increases as the roughness height increases. With the grid spacing of Dy 2 , grid independent solutions at different roughness heights are obtained with the SST model as well. The y þ under clean conditions for this grid spacing is 0.006. The first two roughness heights are in the transitional roughness regime while the third roughness height is in the fully rough regime. The SA roughness modification gives a grid independent solution with a minimum y þ z0:3 whereas the SST roughness model fails to do so in the fully rough regime. This is likely due to how the roughness modification is introduced in the two models. The eddy viscosity at the wall is directly modified in SA but in SST it remains zero. In the fully rough regime, there is a non-zero eddy viscosity in the inner region of the boundary layer where the viscous sub layer previously existed. Since the eddy viscosity is still zero at the wall for the SST roughness modification, to capture the steep increase in eddy viscosity a finer mesh is likely required compared to the SA roughness modification.

Velocity profiles
Despite the finer mesh the skin friction values predicted by the SST roughness model does not match those from the SA model Table 1 Skin friction (C f ) at x ¼ 0:93m for different grid resolutions and roughness levels with SA and SST. N denotes the number of points on the surface of the flat plate.  we can see that the clean case matches the viscous sub-layer and log law in the overlap region closely for both the SA and SST models. Further, increasing the equivalent roughness height has the predicted effect of a shift of the velocity profile away from the clean case and once k þ s > 70, the viscous sub-layer disappears. To further verify the two results, a comparison is made with the empirical shift in velocity profile as proposed by Nikuradse [27] shown below.
where k ¼ 0:40 and the shift B is given by Comparing the empirical predictions of the velocity shift (Fig. 1), a slight overprediction is observed in the transitionally rough region by the SA roughness model. This was also reported in Knopp et al. [13]. The SST roughness model does not perform as well as the SA model especially in the fully rough regime (Fig. 2) despite using a much finer grid. The limitations in the k À u SST roughness model are also reported elsewhere [13,27,28]. It must be noted that various corrections for the SST roughness model have been proposed (for example [13,27]) but are not investigated in the current study.

Blanchard experiments
In this section, the two roughness models are compared to the experimental data from Blanchard obtained from Aupoix et al. [22]. The sand grain roughness height was 4:25 Â 10 À4 m. With an incoming velocity of 45ms À1 , the simulation is carried out on a 2m long flat plate. The resulting Reynolds number is Re ¼ 6:46 Â 10 6 . The y þ of the mesh used is less than 0.4 throughout the domain for the SA roughness model and less than 0.007 for the SST roughness model. The comparison is shown in Fig. 3. Both the SA and SST models predict a higher skin friction compared to the clean flat plate but the results from the SA roughness model are significantly closer to the experimental data. The resulting k þ s z150 makes the flow fully rough. As seen in Fig. 2, the SST roughness model performs poorly in this regime which results in an underprediction of the skin friction.

Roughness on airfoil sections
Determining the damage caused by erosion on turbine blades is an ongoing field of study. Analytical, numerical and experimental studies [29e31] have been carried out to determine the extent of erosion caused by raindrops. The most widely used approach is the droplet impact model. Eisenberg et al. [30] use analytical methods to determine the extent of erosion damage over time due to raindrops (Fig. 4). The effect of continuous exposure of turbine blades to rain is represented by the cumulative number of rain drop impacts and the material removed because of the impacts are modeled based on experimental research. A review on different approaches to study the erosion caused by weather can be found in Keegan et al. [3].
In this section, the focus will be on validating the roughness models against experimental data where the roughness heights are already determined. To that end, two cases are considered: NACA 64 2 215 airfoil at a Reynolds number of Re ¼ 2:6 Â 10 6 and the DU-96-W-180 airfoil at a Reynolds number of Re ¼ 1:85 Â 10 6 . As seen in section 4.1.1, a very fine grid in the wall normal direction is required for the SST roughness model compared to the SA roughness model, which gives grid independent results with meshes comparable to the clean cases. Additionally, despite the fine grid

NACA 65 2 215
In this section the SA model is further validated against the NACA 65 2 215 airfoil. The Reynolds number is Re ¼ 2:6 Â 10 6 and the roughness covers the entire upper surface and on the lower surface from the leading edge up to x=c ¼ 0:15. Three roughness heights k s =c ¼ 1:54 Â 10 À4 , k s =c ¼ 3:08 Â 10 À4 and k s = c ¼ 1:23Â 10 À3 are considered here. Clean experiments were performed by Abbot and von Doenhoff [32]. Ljungstrom performed experiments with different roughness heights on the NACA 65 2 A215 airfoil, a closely related airfoil. These experiments have been used to validate roughness models by Knopp [13] and Hellsten [12] previously. The experimental data are also extracted from Knopp and Hellsten.

Grid details
A two dimensional C-grid topology (Fig. 5) is used for all the simulations. A grid refinement study is carried out at an angle of attack of 8 + on meshes with 150, 250 and 450 nodes on the airfoil surface. A y þ z0:3 is maintained for the three grids. A growth ratio of 1.09 is used within the boundary layer. The computational domain extends to 150 chord lengths in all directions. The grid is shown in Fig. 6. The resulting lift and drag coefficients are listed in Table 2. Since no appreciable difference is observed between the results on the grids with 250 and 450 points (see Table 2), the grid with 250 points on the airfoil was used for further computations. The far field and wall boundary conditions are applied at the edge of the domain and on the airfoil respectively. Fig. 7 shows the comparison of the numerical results from the SA model under clean conditions. The results from SU2 compare very well against results from RFOIL [9] and the experiments from Abbot [32] at lower angles of attack, but SU2 overpredicts the maximum lift. This could be due to a later prediction of the flow separation by the SA turbulence model compared to the experiments. Since no experimental pressure data is available, this cannot be confirmed. However, the lift values reported by Ljungstrom are significantly lower. Since the two airfoils under consideration are supposed to be very similar, Hellsten [12] concludes that lift values reported by Ljungstrom are too low likely due to imperfections from a retracted flap in the airfoil geometry setup. The absolute values of the lift coefficients do not compare well against experimental data from Ljungstrom, but considering the comments of Hellsten the absolute lift coefficient values are not comparable either under clean or rough conditions. The maximum lift is    observed around an angle of attack of 16 + for the clean case in both numerical and experimental data.

Rough results
In Fig. 8 the predicted lift coefficients with different roughness heights are shown. With increasing roughness, the maximum lift value and the angle at which this occurs decrease. Based on the computed skin friction values at an angle of attack of 8 + , k þ s varies from 70 to about 850. These values suggest the flow is likely to be fully rough but it will vary depending on the flow conditions. As noted earlier, the absolute values of the lift coefficients do not match but the relative drop of lift from SU2 matches closely with the experiments (Table 3). However, SU2 predicts a higher value for the angle at which the maximum lift occurs compared to experiments. This is again likely due to the later prediction of the separation location by the SA model.

DU 96-W-180
In this section, the SA roughness model is applied to the DU96-W-180 airfoil. DU96-W-180 is an 18% thick airfoil developed by Delft University [33] and is widely used in the wind energy community. Sareen et al. [7] performed experiments on this airfoil at different Reynolds numbers under different stages and types of erosion. Sareen et al. [7] determine the levels of erosion based on photographs of eroded blades. In this study, the leading edge erosion due to pits and gouges (see Fig. 9) under the two most severe stages are considered at Re ¼ 1:85 Â 10 6 . These cases correspond to Type B stage 3 and stage 4 as reported in Ref. [7].
The depths of pits and gouges are respectively 0:51mm and 2:54mm. The pits and gouges have equal depths and diameters. The chord-wise extent of the pits and gouges are 10% on the upper surface and 13% on the lower surface. The number of pits and gouges on the lower surface is also 1.3 times that on the upper surface. In stage 3 there are 400 pits and 200 gouges on the upper surface and in stage 4 there are 800 pits and 400 gouges on the upper surface.

Grid details
As seen in section 4.1.1, the SA roughness model requires a wall normal grid spacing that corresponds to y þ z0:3 under clean conditions to obtain grid converged results in rough conditions. Thus, this minimum grid spacing is maintained. A grid refinement study is carried out at an angle of attack of 8 + with N ¼ 125; 250; 500 and 750 points along the airfoil. A growth ratio of 1.09 is used in the normal direction. The resulting lift and drag coefficients are listed in Table 4 along with the fully turbulent results obtained from RFOIL [9]. Based on these results the grid with N ¼ 500 points on the airfoil is chosen for further analysis.

Clean results
A baseline case of fully turbulent flow is considered before rough simulations. A transition model is not considered since the effect of roughness on transition is not implemented. Fig. 10 shows the lift coefficient at different angles of attack from SU2 and RFOIL under fully turbulent conditions compared to experimental data. Since no mention of tripping the flow is made in Ref. [7], it is likely that the flow is not fully turbulent but transitional, especially at lower angles of attack. Consequently, a consistent underprediction of lift is observed in both numerical tools. The results from SU2 and RFOIL match closely in the linear region and deviate at higher angles of attack.     Fig. 11 shows the comparison of lift and drag coefficients of the two numerical results from SU2 and RFOIL with the experimental data. Once again, since the experimental flow conditions were not fully turbulent there is a consistent overprediction of the drag coefficient by both SU2 and RFOIL. As seen in Fig. 10, there is good agreement between the numerical results at lower angles of attack. However, RFOIL predicts increasing flow separation to occur from an AoA ¼ 9 + , which is close to what is observed in the experiment but is not predicted by SU2. This is likely due to the poor prediction of separation by the SA turbulence model, which was also observed earlier.

Equivalent sand grain roughness
Roughness height, k, is usually defined as the height or depth of roughness elements on the surface, for example, the depth of pits and gouges in Fig. 9. Determination of the equivalent sand grain roughness height, k s , from the roughness height k is an active area of research. The roughness density parameter, L s , is widely used in literature as a means of relating geometric surface roughness with equivalent sand grain roughness where S is the total wall area where roughness is present, S f is the roughness frontal area, A f is the frontal area of a single roughness element, and A s is the surface area of a single roughness element in the direction of wind flow. Based on data from Schlichtling's experiments, Danberg and Sigal [34] proposed the following relations for 2-D Van Rij et al. [35] generalized the roughness shape factor A f =A s for irregular 3-D roughness elements as S f =S w where S f is the total frontal area of all roughness elements and S w is the total wetted area of all roughness elements and proposed the following relation McClain [36] used the discrete element method approach and proposed a single relation as However, these correlations are mainly derived by adding roughness elements like spheres, cones and hemispheres and their validity for negative roughness like pits and gouges is not clear. Various researchers have used statistical representations of rough surfaces in combination with experiments and numerical simulations using LES and DNS simulations to obtain more general correlations based on rms height (k rms ), skewness (Sk) and higher order moments of the rough surface height probability density functions. The k rms and skewness Sk can be computed as where k i are the heights or depths of individual roughness elements, for example a pit, and N is the total number of such roughness elements. Flack and Schultz [37] proposed k s ¼ 4:43k rms ð1 þ SkÞ 1:37 (21) but note that it is not very general since it does not include information about roughness density. In a more recent study, Flack et al. [38] proposed different relations for different types of skewness as k s ¼ 2:48k rms ð1 þ SkÞ 2:24 (22) for positive skewness, k s ¼ 2:73k rms ð2 þ SkÞ À0:45 (23) for negative skewness and k s ¼ 2:11k rms (24) for zero skewness. They also note that negatively skewed surfaces like those with pits and gouges have a smaller impact on drag than positive skewness due to roughness elements. Flack and Schultz [39]. and Forooghi et al. [40] also note that another parameter that accounts for sparse roughness is necessary and propose a relation of the form where ES is the effective slope which is related to the solidity of roughness (l) as ES ¼ 2l and k z is related to k rms as k z ¼ 4:4k rms . Note that solidity is defined as the ratio of total roughness frontal area (S f ) to total wall area (S). They recommend FðSkÞ ¼ 0:67Sk 2 þ 0:93Sk þ 1:3 (26) and GðESÞ ¼ 1:07 1 À e À3:5ES : (27) In this study, equations (25)e (27) suggested by Ref. [40] are used.

Roughness definition
Sareen et al. [7] create different amounts of roughness on the upper and lower surface with the lower surface being 1.  Table 5. Fig. 12 shows the comparison of the lift coefficient as a function of the angle of attack between SU2 and experiments under stage 3 erosion. There is a small underprediction of lift at lower angles of attack, similar to what was observed in the clean case. This is likely due to the flow still being mildly transitional at lower angles of attack. With increasing angle of attack, the prediction from SU2 matches the experimental data quite closely likely due to the flow becoming fully turbulent in the experiment. Fig. 13 shows the drag and lift coefficients. Once again, the numerical results from SU2 overpredicts the drag compared to the experimental data. Flow separation starts to occur before AoAz 8 + in the experiments whereas SU2 does not predict separation till after AoAz9 + . Fig. 14 shows the comparison of the lift coefficient as a function of the angle of attack between SU2 and experiments under stage 4 erosion. The numerical results agree with the experiments more closely compared to stage 3 likely due to the flow being fully turbulent due to the higher roughness level. Fig. 15 shows the drag and lift coefficients. Once again, numerical results from SU2 overpredict the drag compared to the experimental data. Flow separation is also predicted better with this extent of erosion compared to stage 3.

Rough results
Figs. 13 and 15 also show the lift and drag values in clean conditions. The increase in drag even at lower angles of attack can be seen clearly. The maximum lift also decreases in rough conditions for both roughness levels considered. However, since Sareen et al. [7] do not report lift and drag values at higher angles of attack, the magnitude of reduction cannot be compared. It is very likely that the airfoil will stall earlier for both the roughness cases compared to the clean conditions.
Discussion. In this section the SA roughness model was first validated against experiments on the NACA 65 2 215 airfoil with a given equivalent sand grain roughness. The SA model predicted the drop in lift very closely compared to the experiments. Subsequently, the SA model was used on the DU-96-W-180 airfoil with 'negative' roughness. It was seen that a statistical description of the surface is required to accurately calculate the equivalent sand grain roughness. Results under clean conditions differed from the experiments likely due to the absence of a transition scheme, but the Table 5 Roughness definition for DU-96-W-180 based on Sareen et al. [7].  Table 5). Fig. 13. Comparison of lift coefficient (C l ) against drag coefficient (C d ) for fully turbulent flow against experimental data (stage 3 see Table 5).  Table 5). numerical results, especially lift coefficient, matched closely with the experimental data under roughness when the flow is likely fully turbulent. It was seen that roughness causes a considerable reduction in lift and increase in drag and can lead to premature stalling of the airfoils.

Boundary layer analysis
Since the NACA 65 2 215 airfoil has a larger rough surface than the DU-96-W-180, it is chosen for the boundary layer analysis. The boundary layer parameters for SU2 are computed by extracting the velocity vector along surface normals at various points along the airfoil. The edge of the boundary layer is assumed to be at the location where the ratio of the magnitude of the vorticity at that location to the value at the wall is less than 10 À4 . In this section the effect of roughness on the boundary layer properties of airfoils will be investigated.

Integral boundary layer methods
The combination of integral boundary layer (IBL) methods with panel methods known as viscous-inviscid interaction (VII) are very commonly used in aerodynamic analysis of airfoils. XFOIL [41] and RFOIL [9] are some of the widely used tools based on this approach. Viscous-inviscid interaction methods give very accurate results at a very low computational cost compared to standard RANS simulations. They are commonly used during the design process of wind turbine blades in combination with blade element momentum (BEM) theory and other rotor design methods. A roughness model for the integral boundary layer methods will allow for a more accurate and quick assessment of the effect of roughness on turbine performance during the design phase.
Integral boundary layer equations are obtained by integrating the boundary layer equations in the direction normal to the wall. More details on deriving the governing equations can be found in Refs. [42,43]. The new integral quantities introduced are displacement thickness d * , momentum thickness q and kinetic energy thickness d k .
Here u is the local velocity, d is the boundary layer thickness, u e is the velocity magnitude at the edge of the boundary layer and y is the wall normal direction. Further, the following shape parameters are defined The governing equations resulting from integration of the continuity and momentum equations used in RFOIL are Note that other formulations of the integral boundary layer equations are used in other tools [43]. In order to close the equations, closure relations [42,43] are defined for the kinetic energy shape factor H k , the skin friction coefficient C f and the dissipation coefficient C D . The closure relations are different for laminar and turbulent flows. For turbulent flows, an additional equation for lag in Reynolds shear stress (C t ) is also solved. H ** is a shape factor based on the variation of density within the boundary layer and M e is the Mach number of the external flow. Both can be ignored for incompressible flows. These closure relations are defined in terms of the shape factors introduced earlier and the Reynolds number based on momentum thickness Re q , where Re q ¼ u e q=n with n being the kinematic viscosity. In the following sections, the effect of roughness on the different thicknesses, shape factors and closure relations are examined.

Clean results
First the calculated integral boundary layer quantities from SU2 under clean conditions are compared against the RFOIL results. It must be noted that the X À axis of all the plots in this section range from x=c ¼ 0:025 to x=c ¼ 1 to avoid the stagnation region. Fig. 16 shows the displacement thickness on both the upper and lower sides at angles of attack of 0 + and 4 + . The calculated displacement thickness matches closely with the values from RFOIL with some deviation near the trailing edge in both cases.
The momentum thickness is slightly overpredicted by SU2 after x=c ¼ 0:4 at an angle of attack of 0 + but matches closely for an angle of attack of 4 + as seen in Fig. 17.
The comparisons of the shape factors are shown in Fig. 18. The shape factor is larger for AoA ¼ 4 + compared to AoA ¼ 0 + indicating a thicker boundary layer as the angle of attack increases. While the computed shape factors from RFOIL and SU2 do not match exactly, both display similar behavior initially decreasing towards the middle of the airfoil and increasing near the trailing edge.
Here Re q is the local Reynolds number based on momentum thickness q. Fig. 19 shows the comparison of the skin friction coefficient between RFOIL, the values reported by SU2 originally by the RANS computation (denoted as 'SU2 original') and the skin friction calculated based on the computed integral boundary layer  Table 5).

Rough results
Since the entire upper surface is rough the results for the upper surface only are presented in this section. Figs. 20 and 21 shows the displacement and momentum thickness for different roughness levels compared to the clean case. As expected, these thicknesses increase with increasing roughness. A very steep increase is observed in the momentum thickness near the trailing edge for the largest roughness. Fig. 22 shows the shape factors for different roughness levels compared to the clean case at AoA ¼ 0 + and AoA ¼ 4 + . The shape factor increases for all roughness levels with the largest increase for k s ¼ 1:23 Â 10 À3 . The maximum k þ s values varies with angle of attack. At an angle of attack of 0 + , the k þ s are 25, 75 and 286 indicating that the flow is in the transitional rough region for the two lower roughness levels and is fully rough for the highest roughness level. However, at an angle of attack of 4 + , the maximum k þ s values are 75, 180 and 750 indicating that the flow is fully rough for the k s = c ¼ 3:08 Â 10 À4 case also. From Fig. 22 it is seen that the behavior of the shape factor in the k s =c ¼ 3:08 Â 10 À4 case is similar for both angles of attack despite one being transitionally rough and the other fully rough.

Skin friction coefficient
Equation (32) will not be valid here as the properties of the boundary layer change due to roughness. Olsen et al. [10] suggested a new closure relation for skin friction for rough surfaces including the Reynolds number based on roughness height, Re k ¼ u e k= n as ðjlog 10 Re q À log 10 Re k þ 1:11jÞ 2:45À0:15H : Fig. 23 shows the skin friction from equations (32) and (33), clean and rough SU2 results at angles of attack of 0 + (top) and 4 + (bottom). Clearly equation (32) is not valid for rough surfaces. The new closure relation proposed by Olsen et al. appears to overpredict the skin friction. However, since the computed Re k for the first two roughness levels are approximately 400 and 800, it is outside the range of the data used by the authors in their study. The third roughness level has an average Re k z3000 and is within the valid range of data used to derive the model. The authors report convergence difficulties when roughness was applied to regions before x=c ¼ 0:6 and from Fig. 23 it can be seen that C f is overpredicted by a significant amount in that region and is closer to the values reported by SU2 after x=c ¼ 0:6.

Kinetic energy shape factor
As seen above the closure sets for skin friction are not valid for rough airfoils. The other closure relation defined in terms of H and The computed H k based on equation (29) (denoted by symbols) and those based on the closure relations in equations (35) and (36) (denoted by solid lines) are shown in Fig. 24. The computed values agree with the closure relations closely for the clean case and also for the two lowest roughness cases. However, as the level of roughness increases the closure relation does not predict H k accurately. The Re ks of the first two roughness cases are approximately 400 and 800, indicating that the closure sets are likely valid for small roughness levels but deviate for higher roughness levels. The deviation observed in the third roughness level is also much less than the deviation observed for the skin friction coefficient. Fig. 25 shows the variation of H k for a higher angle of attack of 12 + . From the top figure it is seen that the behavior of H k is similar to that observed for lower angles of attack when the flow is attached. However, as the bottom figure shows, the deviation increases for all roughness levels when the flow separates. The wiggles observed are likely an artefact of how the edge of the boundary layer is detected during the post processing.
Since the closure relations for the dissipation coefficient (C D ) and for turbulent flows the Reynolds shear stress coefficient (C t ) are based on H, Re q , C f and H k , all of which change with roughness, new closure relations need to be defined. Thus, in order to model roughness in integral boundary layer method based tools like RFOIL, new closure relations need to be derived for all of the above quantities.      (29)) shown as symbols and result from the closure relations (equations (35) and (36)) as solid lines.

Discussion
In this section, the different integral boundary layer quantities and closure relations used in RFOIL are examined under clean and rough conditions. Three different roughness levels were considered corresponding to Re ks of approximately 400, 800 and 3000. The boundary layer thicknesses increase due to roughness and the shape factor is also higher. A larger shape factor typically implies a thicker boundary layer that is prone to separation. The shape factor remained less than 2 for the lower two roughness levels but for the highest roughness level the shape factor neared the value for separation even at low angles of attack. Additionally, it was seen that the variation of the shape factor followed the variation of Re ks more closely than the variation of k þ s . Closure relations are crucial for an accurate solution in integral boundary layer methods. The performance of the closures for skin friction and kinetic energy shape factor was examined under rough conditions. The closure relation for skin friction underpreformed significantly for all cases. The new closure relation proposed by Olsen et al. [10] was observed to overpredict the skin friction. The kinetic energy shape factor closure relation was less sensitive to roughness and showed significant deviation only for the largest roughness case and separated flow.

Conclusions and future outlook
Leading edge erosion causes a reduction in aerodynamic performance of wind turbine blades. Depending upon the extent of roughness, a drop in maximum lift of up to 30% can be observed. Skin friction increase was observed in all cases leading to an increase in drag. Effect of roughness is modeled using the equivalent sand grain roughness height. Roughness models for the two RANS turbulence models SA and SST were implemented in SU2 and the accuracy was examined via grid refinement. The models were validated against empirical models for the shift in velocity profiles in the boundary layer and experimental skin friction data on flat plates. It was seen that the SST roughness model required a much finer grid compared to the SA roughness model to give a grid independent solution. However, despite the finer grid the results from the SST roughness model did not match the experimental data or the empirical models under fully rough conditions, unlike the SA roughness model. Based on these results the SA roughness model was further validated against experimental data on two different airfoils. The SA model predicted the reduction in lift for different roughness levels accurately for the NACA 65 2 215 airfoil. The SA model was also validated against an experiment with negative roughness (pits and gouges) on the DU-96-W-180 airfoil.
Encouraging results were observed for both roughness levels tested. The statistical method to determine the equivalent sand grain roughness proved to be accurate. Some differences were observed in the clean simulation, most likely due to the fact the simulations were run under fully turbulent conditions, unlike the experiments.
Further, the behavior of different integral boundary layer properties like displacement thickness, momentum thickness, shape factors and closures were investigated for the NACA 65 2 215 airfoil. The existing skin friction closure relations for clean surfaces greatly underpredict the skin friction (C f ) and are not valid for rough surfaces. However, the closure relation for the kinetic energy shape factor (H k ) performed well for low roughness levels (Re ks < 1000Þ) but deviated at higher roughness levels and under separation. The deviation was only marginal compared to the skin friction closure relation. However, since the closure relations for other quantities like the dissipation coefficient and the Reynolds shear stress coefficient depend on C f and H k new closure relations will be needed in order to simulate rough surfaces in integral boundary layer tools like RFOIL.
The main focus of this study was on the effect of roughness on turbulent boundary layers. For laminar boundary layers, roughness leads to premature transition and for turbulent boundary layers. In order to fully capture the effect of roughness, the effect on transition will also be considered in the future. Further, more boundary layer data at different roughness levels are needed to derive new closure relations for integral boundary layer methods.