Numerical simulations of ice accretion on wind turbine blades: are performance losses due to ice shape or surface roughness?

Ice accretion on wind turbine blades causes both a change in the shape of its sections and an increase in surface roughness. These lead to degraded aerodynamic performances and lower power output. Here, a high-fidelity multi-step method is presented and applied to simulate a 3 h rime icing event on the National Renewable Energy Laboratory 5 MW wind turbine blade. Five sections belonging to the outer half of the blade were considered. Independent time steps were applied to each blade section to obtain detailed ice shapes. The roughness effect on airfoil performance was included in computational fluid dynamics simulations using an equivalent sand-grain approach. The aerodynamic coefficients of the iced sections were computed considering two different roughness heights and extensions along the blade surface. The power curve before and after the icing event was computed according to the Design Load Case 1.1 of the International Electrotechnical Commission. In the icing event under analysis, the decrease in power output strongly depended on wind speed and, in fact, tip speed ratio. Regarding the different roughness heights and extensions along the blade, power losses were qualitatively similar but significantly different in magnitude despite the well-developed ice shapes. It was found that extended roughness regions in the chordwise direction of the blade can become as detrimental as the ice shape itself.

1 Introduction

Arguably, wind energy will lead the energy transition in Europe. In 2022, it produced 15 % of the total generated electricity (Jones et al., 2023). Its share is due to increase to approximately 30 % by 2030 and 45 % by 2050 to reach carbon neutrality (European Commission, 2018). Future installations of wind turbines will mainly occur in cold regions. Higher wind speeds and air density guarantee a higher wind power density, leading to a competitive cost of energy and viable investments (European Commission, 2012). However, two phenomena may reduce the power output in cold climates when considering a single wind turbine or a wind farm. First, cold climates often lead to a stable atmospheric boundary layer, causing less mixing of the wake of a wind turbine. Thus, wake effects are stronger, and the wind turbines located downwind in wind farms produce less power. Second, ice can form on wind turbines when clouds or super-cooled fog arise at low elevations and temperatures drop below 0 ^∘C. These conditions may persist for days or even weeks, and, due to climate change, they may occur in previously unexpected locations. In February 2021, three consecutive winter storms swept over Texas and caused a more than 80 % reduction in energy output compared to the previous week. It is shown in Fig. 1.

Figure 1 Texas region electricity generation by wind energy in February 2021. Winter storms occurred between 10 and 20 February. Source: US Energy Information Administration, Hourly Electric Grid Monitor, available at https://www.eia.gov/electricity/gridmonitor (last access: 10 January 2022).

Indeed, ice may affect wind turbines in several ways. The first visible effect is the blade aerodynamics degradation, reducing the wind turbine power output. Instrumentation and controller errors may follow. As more ice is accreted, the structural behaviour changes as well, and the fatigue life of the structure can be affected. Ice shedding may also be a severe threat, endangering equipment and people nearby and causing significant load imbalances on the rotor. In many cases, a shutdown may become unavoidable. For these reasons, wind turbines operating in cold regions must be equipped with an ice protection system (IPS). Electro-thermal IPSs provide a possible solution. These devices are energy-consuming, especially if run in anti-icing mode (Caccia et al., 2022). Recent design solutions combine the effects of centrifugal force and IPS heat to remove ice from the blades (Getz and Palacios, 2021). Better predictions of power losses due to icing may lead to improved designs and further energy savings during the operation of such devices.

The problem of ice accretion on wind turbines has been long studied but is still of great interest. In 2002, the International Energy Agency (IEA) established a cooperation group of international experts, Task 19, to study wind energy in cold climates. The working group has been running ever since. A recent report by the Task 19 group reviewed the technologies available for wind turbines operating in cold climates (Lehtomäki et al., 2018). Many fundamental aspects were covered, including ice detection systems, ice protection systems, ice accretion models, ice shedding, operation and maintenance, standards, and testing. Ice accretion models were categorised into two classes, i.e. simplified and advanced. Simplified models are typically used within weather models to estimate icing rate and mass on objects. They are based on empirical relationships to simulate ice growth on cylinders. On the other hand, advanced models are used for aerodynamic analyses by capturing detailed physics of the ice accretion process to produce accurate 2D or 3D ice shapes. State-of-the-art advanced ice accretion models were benchmarked in 2022 at the first AIAA Ice Prediction Workshop for code-to-code and code-to-experiment comparisons, including FENSAP-ICE (Ozcer et al., 2022), GlennICE (Wright et al., 2022), IGLOO3D (Radenac and Duchayne, 2022), and PoliMIce (Morelli et al., 2022). A review of the results was written by Laurendeau et al. (2022). In this work, we analysed the aerodynamics of iced blade sections in detail. So, an advanced ice accretion model was used with the in-flight ice accretion engine PoliMIce (Gori et al., 2015).

Depending on atmospheric conditions, different types of ice can form. Their standard classification is rime, glaze, and mixed ice. The current study focused on rime ice, formed when super-cooled water droplets freeze instantly upon impact. The type of ice and the accretion rate depend on various parameters. According to Etemaddar et al. (2014), such parameters can be divided into two categories, i.e. atmospheric and system parameters. The atmospheric parameters are the freestream temperature (T_∞); the liquid water content (LWC), defining the amount of water dispersed in a reference volume of air; and the median volumetric diameter (MVD), defining the droplet diameter above and below which half the volume of water is contained. The system parameters are the geometric parameters of the object, i.e. shape, orientation, and dimension (shape, angle of attack, chord c, and thickness t of an airfoil), and the relative wind speed V_rel. The formation of rime ice is favoured by low T_∞, low LWC, small MVD, and low V_rel. Moreover, the rate of ice accretion increases with increasing LWC, MVD, and V_rel and decreasing $t/c$. As a result, the ice mass accreted on a wind turbine blade increases from the root to the tip. For a more in-depth analysis, the reader can refer to the works by Etemaddar et al. (2014), Homola et al. (2010a, b), and Virk et al. (2010).

When ice is formed on an airfoil, two main factors alter its performance: the ice shape and the increase in surface roughness. While the former can be assessed numerically, the latter needs to be estimated either experimentally or using empirical correlations. During ice accretion, the roughness height evolves in a complex way with both space and time (Steiner and Bansmer, 2016; McClain et al., 2020). In computational fluid dynamics (CFD) simulations, it is common practice to map the actual roughness distribution into the so-called equivalent sand-grain roughness (Nikuradse, 1950). For this type of roughness, the shift in the velocity profile in the logarithmic region of the boundary layer is known as a function of the roughness height, k_s, in wall units ($k_{\mathrm{s}}^{+}$). This topic will be addressed later.

Several numerical studies on power losses due to ice accretion are available in the literature. An exhaustive review was carried out by Contreras Montoya et al. (2022). The procedure for estimating power losses usually relies on the blade element momentum (BEM) theory. It can be summarised as follows: (1) computation of the aerodynamic coefficients of the clean airfoils, (2) simulation of steady-state ice accretion on relevant 2D sections, (3) computation of the aerodynamic coefficients of the iced airfoils, and (4) computation of the power curves. Indeed, the numerical study of 2D sections can provide an affordable and reliable approximation of 3D ice accretion. In terms of ice shape, when ice accretion is coupled with the BEM solution, the two approaches are almost equivalent (Switchenko et al., 2014). In terms of airfoil performance, an experimental and numerical study on the 3D scan of an iced wind turbine airfoil was carried out by Knobbe-Eschen et al. (2019). The authors found that accurate predictions of the numerical 3D solution can be obtained by studying 2D slices of the actual ice shape, i.e. including the localised macroscopic roughness. On the other hand, the airfoil performance was over-predicted when the ice shape was span-wise averaged, i.e. when macroscopic and microscopic roughness were removed. In this context of 2D BEM simulations of ice accretion, three representative works on the National Renewable Energy Laboratory (NREL) 5 MW reference wind turbine (Jonkman et al., 2009) are reported here.

One of the first efforts to study ice accretion numerically on a wind turbine was carried out by Homola et al. (2012). Five sections belonging to the whole blade span were studied. A BEM code was used to evaluate both boundary conditions for each section and wind turbine performances. The wind turbine was operating with V_∞=10 m s⁻¹, $T=-\mathrm{10}$ ^∘C, LWC=0.22 g m⁻³, and MVD=20 µm for a total time of 1 h. Roughness was applied on the entire blade surface. The equivalent sand-grain roughness height $k_{\mathrm{s}}/c$ was not specified. The power loss was about 25 % between 7 and 11 m s⁻¹ in a steady wind.

In 2013, Turkia et al. (2013) studied power losses on a down-scaled version of the NREL 5 MW turbine to achieve a rated power of 3 MW. Two sections belonging to the outer third of the blade were studied. The considered atmospheric conditions were V_∞=7 m s⁻¹, $T_{\mathrm{\infty }}=-\mathrm{7}$ ^∘C, LWC=0.2 g m⁻³, and MVD=25 µm for a total time of 10 h. Roughness was included in two different ways. A large-scale roughness was applied to the predicted ice shapes, while the effect of small-scale roughness was included by correcting the drag coefficients of the iced airfoils with a method proposed by Bragg (1982) for rime ice. The relation depends on equivalent sand-grain roughness height $k_{\mathrm{s}}/c$, estimated between $\mathrm{0.9}\times {\mathrm{10}}^{-\mathrm{3}}$ and $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}$ using Shin's relation (Shin et al., 1991). However, Bragg's relation was developed to estimate the drag coefficient of an iced airfoil in the presence of leading-edge roughness using the clean C_D as input and not the iced one. As a result, the effect of drag on each section was considered twice. At the end of the icing event, power losses were approximately 25 % in the 6–12 m s⁻¹ range and occurred up to 17 m s⁻¹. No information was provided on the type of wind used as input.

Finally, Etemaddar et al. (2014) simulated a 24 h icing event on the NREL 5 MW turbine, varying the atmospheric conditions in the considered time window and sampling them every 15 min. Six sections in the outer half of the blade were considered. Ice accretion was simulated with LEWICE (Wright, 2008). Roughness was applied to the first 25 % of the airfoil chord, and $k_{\mathrm{s}}/c$ was set to $\mathrm{0.5}\times {\mathrm{10}}^{-\mathrm{3}}$. Wind turbine operation was simulated with HAWC2 (Larsen and Hansen, 2007) with the Mann spectral tensor model for atmospheric turbulence (Mann, 1994). In this case, power losses were about 45 % at cut-in wind speed, reducing to 34 %, 23 %, and 1.8 % at 7, 11, and 16 m s⁻¹, respectively. Icing loads were studied as well.

These results do not show a clear trend in power losses. At 11 m s⁻¹, the reduction in extracted power was comparable regardless of the duration of the icing event. Moreover, the trend found for increasing wind speeds was different. Homola et al. (2012) and Turkia et al. (2013) predicted a slight increase in power loss from 7 to 11 m s⁻¹, while Etemaddar et al. (2014) found a drastic decrease. Indeed, different atmospheric conditions led to different ice shapes, resulting in different aerodynamic performances of the airfoils. However, in each work, the roughness was taken into account profoundly differently, which may have led to this pattern of results.

A study by Switchenko et al. (2014) supports this hypothesis. In their work, the authors numerically simulated a real-world icing event on a 1.5 MW wind turbine. Two important conclusions were drawn from this study. The first one, as already mentioned, was that simulating ice accretion on 2D sections with the BEM methodology provides very similar results to the 3D solution. The second one was that “roughness of the ice can at times be more significant than the actual size, shape and placement of the accreted ice”, and so “more research is needed to evaluate the effect of atmospheric icing on wind turbine blades and their surface roughness characteristics” (Switchenko et al., 2014). However, the authors never computed the airfoil polars since the BEM method was integrated within iterative CFD simulations. Thus, the nature of power losses remained unknown. On the other hand, the detrimental effect of roughness on the aerodynamic coefficients of a wind turbine airfoil was shown by Blasco et al. (2017) through ice accretion experiments and wind tunnel measurements. Ice accretion time was set to 45 min to obtain streamlined ice shapes. The maximum lift decrease and drag increase at an operational angle of attack were about 25 % and 220 %, respectively. However, these results are different from those by Switchenko et al. (2014), in which the effect of roughness was highlighted for long-lasting icing events.

In view of the above, the aim of this work is to (1) perform a high-fidelity ice accretion simulation on the NREL 5 MW wind turbine blade, (2) compute the aerodynamic performances of the iced blade sections as a function of roughness, and (3) assess the effect of icing in operating conditions. The icing event was long enough for streamlined, protruded ice shapes to form and the effects of ice shapes and roughness to be combined. The work was carried out using both open-source software and in-house codes. We presented a preliminary work in this context in Caccia et al. (2021). Original contributions to the state of the art include the introduction of span-dependent time stepping in the ice accretion simulation and an analysis of the sensitivity of the solution to roughness height and extension. It will be shown that roughness can significantly affect airfoil aerodynamics and power production also when complex ice shapes are present.

The paper is structured as follows. The methodology is presented in Sect. 2, together with the set-up of the numerical simulations. In Sect. 3, the numerical set-up is compared with experimental data. The aerodynamic coefficients of the clean airfoils and an icing experiment on a rotating blade section are reproduced. In Sect. 4, the results of ice accretion are presented. Then, the aerodynamic coefficients of the iced sections are computed and used to simulate the power curve of the wind turbine. Finally, conclusions are drawn in Sect. 5.

2 Methodology

The NREL 5 MW reference wind turbine was analysed in this work. The blade structural and aerodynamic designs are based on the Dutch Offshore Wind Energy Converter (DOWEC) project (Kooijman et al., 2003). According to a more detailed design of the blade (Resor, 2013), it is an International Electrotechnical Commission (IEC) Class I and Category B wind turbine.

The aero-servo-elastic response of the wind turbine was modelled with OpenFAST (https://github.com/OpenFAST/openfast, last access: 20 September 2020). Wind turbine aerodynamics was modelled through blade element momentum theory. Thus, 2D, independent sections were analysed throughout the whole work. Within the BEM framework, the blade was discretised with 19 nodes along the blade span. It was made of two cylindrical sections, five Delft University (DU) airfoils, and one National Advisory Committee for Aeronautics (NACA) airfoil. The aerodynamic coefficients of DU airfoils were measured by Ruud van Rooij of the Delft University of Technology at a Reynolds number (Re) of 7×10⁶. NACA 64₃-618 coefficients were taken from Abbott et al. (1945) at a Reynolds number of 6×10⁶. The 2D data of the airfoils are provided in the DOWEC report. A review of the airfoils distribution along the blade, their thickness, and the tested Reynolds number is presented in Table 1.

Table 1 Airfoils composing the wind turbine blade.

The first step involved reproducing these experimental data with a CFD solver, SU2 (Economon et al., 2015). This permitted the validation of the solver set-up. Reynolds-averaged Navier–Stokes (RANS) equations were solved. Then, the numerical data were extrapolated to the entire 360^∘ range of angle of attack (AoA) values and corrected for 3D effects using NREL's tool “AirfoilPrep” to prepare them for the BEM model of the wind turbine. In particular, the Viterna method (Viterna and Janetzke, 1982) was used for extrapolation, considering an aspect ratio of 17. The 3D corrections by Du and Selig (1998) with Eggers C_D adjustment (Eggers et al., 2003) were applied only for positive AoA values and considering the rated rotational velocity as input, being consistent with NREL's aerodynamic data of the wind turbine model (Jonkman et al., 2009).

The second step was the simulation of the icing event. The atmospheric conditions of the icing event on the wind turbine are reported in Table 2. The same conditions were studied by Homola et al. (2012) for 1 h and by Zanon et al. (2018) for approximately 8 h. In the current contribution, the icing event lasted for 3 h. A wind shear exponent of 0.15 was also considered so that, according to the normal wind profile model of the DNV GL guidelines (Germanischer-Lloyd, 2010), the average horizontal wind speed V as a function of the height above the ground z is

$$\begin{array}{}\text{(1)}& V\left(z\right)=V_{\mathrm{hub}}{\left(z/z_{\mathrm{hub}}\right)}^{\mathrm{0.15}},\end{array}$$

where V_hub is the wind speed V at the hub height z_hub. These atmospheric conditions led to the formation of rime ice on the blade surface. During the icing event, wind turbine operation was computed using OpenFAST to find the equilibrium condition of the whole system, considering wind shear and blade deformability. A steady wind was assumed at this stage. A retroaction on the wind turbine operating state due to a change in the blade aerodynamics was included. However, such retroaction would have been influential on a longer timescale or considering a more penalising roughness during ice accretion. Local boundary conditions were evaluated on significant blade sections and used as input for ice accretion. A multi-step approach was adopted, dividing the total accretion time into sub-intervals. For each section, time steps were set independently according to the required time discretisation, taking full advantage of the BEM approximation. The ice accretion engine PoliMIce (Gori et al., 2015) was used for these simulations. The software uses SU2 for computing the aerodynamic field and PoliDrop (Bellosta et al., 2019) for Lagrangian particle tracking. The numerical set-up chosen for ice accretion was validated against nine experimental test cases on a rotating section by Han et al. (2012).

Table 2 Atmospheric conditions during the icing event studied.

The third step was the computation of the aerodynamic coefficients of the iced sections with SU2. Roughness was modelled using an equivalent sand-grain approach. Two values of k_s were compared. Moreover, the effect of the extension of the rough region along the surface of each iced section was assessed. Once more, the numerical data were extrapolated to the 360^∘ AoA range and corrected for 3D effects.

The C_P–TSR (tip speed ratio) curves were then computed using “AeroDyn”, i.e. the aerodynamic module of OpenFAST, using the blade sectional aerodynamic coefficients and the blade geometry. The power coefficient C_P is the ratio between the power extracted by the wind turbine and a conventional freestream available power:

$$\begin{array}{}\text{(2)}& C_{\mathrm{P}}={\displaystyle \frac{P}{\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }\left(\mathit{\pi }{R}^{\mathrm{2}}\right)V_{\mathrm{\infty }}^{\mathrm{3}}}},\end{array}$$

where R is the rotor radius and V_∞ the freestream velocity. By neglecting blade elasticity, it can be shown that the quantity depends only on two system parameters, i.e. the blade pitch angle and the tip speed ratio $\mathrm{TSR}=V_{\mathrm{tip}}/V_{\mathrm{\infty }}$, where V_tip is the blade tip velocity. At a certain V_∞, the controller makes the wind turbine work at a specific TSR and with a specified blade pitch. Thus, the C_P–TSR curve can provide valuable insight into power production and power losses of wind turbines.

The final step was the computation of the power curve of the wind turbine before and after the icing event with OpenFAST. Both steady and turbulent winds were considered. Atmospheric turbulence was modelled as defined by the IEC for the Design Load Case (DLC) 1.1 of a Category B wind turbine (IEC, 2005) using the IEC Kaimal spectral model. The realisations of the turbulent wind were generated with TurbSim (Jonkman, 2009), considering a reference turbulence intensity of 0.16. A Weibull-averaged power was computed to provide a single figure of the severity of the icing event under the two sets of roughness under analysis. This was simply computed as

$$\begin{array}{}\text{(3)}& P_{\mathrm{W}}=\underset{V_{\mathrm{in}}}{\overset{V_{\mathrm{out}}}{\int }}P\left(V\right)f_{\mathrm{w}}\left(V\right)\mathrm{d}V,\end{array}$$

where V_in is the cut-in wind speed, V_out the cut-out wind speed, P(V) the power curve, and f_w(V) the Weibull probability density function. The latter is defined as

$$\begin{array}{}\text{(4)}& f_{\mathrm{w}}\left(V\right)={\displaystyle \frac{k}{c}}{\left({\displaystyle \frac{V}{C}}\right)}^{k-\mathrm{1}}\exp \left[-{\left({\displaystyle \frac{V}{C}}\right)}^k\right],\end{array}$$

where k is a shape parameter and C is a scale parameter. The shape parameter k was set to 2 to match the Rayleigh distribution advised in the DNV GL guidelines for wind turbine certification (Germanischer-Lloyd, 2010). The scale parameter C was set to 11.2838 m s⁻¹ so that the average wind speed is V_avg=10 m s⁻¹, as prescribed for a Class I wind turbine. The two quantities are linked by the relation

$$\begin{array}{}\text{(5)}& V_{\mathrm{avg}}(C,k)=C\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{e}^{-t}{t}^{\mathrm{1}/k}\mathrm{d}t,\end{array}$$

where the integral term is the gamma function $\mathrm{\Gamma }\left(x\right)=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{e}^{-t}{t}^{x-\mathrm{1}}\mathrm{d}t$ evaluated in $x=(\mathrm{1}+\mathrm{1}/k)$.

2.1 Law of the wall with roughness

Two main factors change the performance of an iced airfoil: the modification of the airfoil shape and the increase in surface roughness. While the former is computed numerically, the latter is usually estimated with empirical correlations. During ice accretion, the roughness height evolves in a complex way with both time and space (Steiner and Bansmer, 2016; McClain et al., 2020). It is common practice in CFD simulations to map the real roughness distribution into the ideal roughness studied by Nikuradse (1950), known as equivalent sand-grain roughness. For this type of roughness, the shift in velocity profile in the logarithmic region of the boundary layer is known as a function of the roughness height k_s in wall units ($k_{\mathrm{s}}^{+}$). The following relation holds:

$$\begin{array}{}\text{(6)}& u^{+}={\displaystyle \frac{\mathrm{1}}{\mathit{\kappa }}}\log \left({\displaystyle \frac{y^{+}}{k_{\mathrm{s}}^{+}}}\right)+B\left(k_{\mathrm{s}}^{+}\right),\end{array}$$

where u⁺ is the non-dimensional tangential velocity in wall units, y⁺ is the non-dimensional wall distance in wall units, κ is the von Kármán constant (κ≈0.41), and $B\left(k_{\mathrm{s}}^{+}\right)$ is an additive constant. In particular, y⁺, $k_{\mathrm{s}}^{+}$, and u⁺ are defined as

$$\begin{array}{}\text{(7)}& {\displaystyle }& {\displaystyle }{y}^{+}={\displaystyle \frac{y}{{\mathit{\delta }}_{\mathit{\nu }}}},\text{(8)}& {\displaystyle }& {\displaystyle }{k}_{\mathrm{s}}^{+}={\displaystyle \frac{k_{\mathrm{s}}}{{\mathit{\delta }}_{\mathit{\nu }}}},\text{(9)}& {\displaystyle }& {\displaystyle }{u}^{+}={\displaystyle \frac{u}{u_{\mathit{\tau }}}},\end{array}$$

where ${\mathit{\delta }}_{\mathit{\nu }}=\frac{\mathit{\nu }}{u_{\mathit{\tau }}}$ is the viscous length scale, $u_{\mathit{\tau }}=\sqrt{\frac{{\mathit{\tau }}_{\mathrm{w}}}{\mathit{\rho }}}$ is the friction velocity (the velocity scale of the turbulent fluctuations at the wall), ν is the dynamic viscosity, τ_w is the wall shear stress, and ρ is the fluid density.

Depending on the value of $k_{\mathrm{s}}^{+}$, different roughness regimes are defined. The smooth regime is defined for $k_{\mathrm{s}}^{+}<\mathrm{5}$. In this case, the roughness elements are submerged within the viscous sublayer ($y^{+}<\mathrm{5}$), and its effect on the flow is negligible. In the viscous sublayer the relation $u^{+}=y^{+}$ holds. For $y^{+}\gtrsim \mathrm{30}$, the clean law of the wall is recovered:

$$\begin{array}{}\text{(10)}& u^{+}={\displaystyle \frac{\mathrm{1}}{\mathit{\kappa }}}\log \left(y^{+}\right)+\mathrm{5.1}.\end{array}$$

In the transitionally rough regime ($\mathrm{5}<k_{\mathrm{s}}^{+}\lesssim \mathrm{70}$), the additive constant $B\left(k_{\mathrm{s}}^{+}\right)$ (Eq. 6) depends on $k_{\mathrm{s}}^{+}$. In the fully rough regime ($k_{\mathrm{s}}^{+}\gtrsim \mathrm{70}$), typical of ice, $B\left(k_{\mathrm{s}}^{+}\right)$ becomes independent of $k_{\mathrm{s}}^{+}$. Its value is equal to ∼8.0, according to Schlichting and Gersten (2017). However, the estimation of a single k_s value from a time- and space-dependent roughness distribution is not trivial. Empirical relations have been developed in the aeronautic industry. Adopting them for wind turbines is common, even if icing occurs in very different environmental conditions.

Although the use of the relation for k_s developed by Shin et al. (1991) for the LEWICE code is still widespread, the code now implements a newer relation by Wright (2008). Shin's relation was specifically developed to match the ice shapes predicted by the LEWICE code to experimental ones. For this reason, it may lack generality. On the other hand, Wright's relation “was determined from experimental measurements of roughness heights as a function of the calculated freezing fraction at the stagnation point. It was not reverse-engineered in order to match ice shape predictions” (Wright, 2008). The relation is

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{k_{\mathrm{s}}}{c}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1000}}}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sqrt{\mathrm{0.15}+{\displaystyle \frac{\mathrm{0.3}}{f_{\mathrm{0}}}}},\end{array}$$

where f₀ is the freezing fraction at the stagnation point, equal to 1 for rime ice (all the impinging water freezes upon impact). Thus, according to this equation, $\frac{k_{\mathrm{s}}}{c}=\mathrm{0.34}\times {\mathrm{10}}^{-\mathrm{3}}$ for the case under analysis. For comparison, Shin's relation provides $\frac{k_{\mathrm{s}}}{c}=\mathrm{0.58}\times {\mathrm{10}}^{-\mathrm{3}}$ for convective heat transfer and $\frac{k_{\mathrm{s}}}{c}=\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{3}}$ for drag prediction.

However, these relations do not include the effect of time on k_s. They were specifically developed for the atmospheric conditions typical of aviation, where icing occurs at high wind speeds for a short time. On the other hand, ice accretion on wind turbines may last for hours. Since the roughness height increases with time, much higher roughness can be found on wind turbine blades. However, such values are currently unknown.

2.2 CFD simulation set-up

The SU2 code solves RANS equations using an edge-based finite volume discretisation in space on general unstructured grids. The convective and viscous fluxes are then evaluated at the midpoint of an edge. An upwind flux difference splitting (FDS) numerical scheme was chosen to solve the convective fluxes in the incompressible solver. Second-order accuracy of the numerical method was obtained by applying a monotonic upstream-centered scheme for conservation laws (MUSCL) for flux reconstruction. The gradients of the variables at each node were reconstructed using the Green–Gauss theorem. During reconstruction, gradients were limited using the slope limiter by Venkatakrishnan (1995) to avoid spurious oscillations of the variables.

Steady-state problems are solved with a pseudo-time-step technique, in which the solution is marched in time until the time derivative term vanishes and a steady-state solution is reached. An adaptive Courant–Friedrichs–Lewy (CFL) method was used for convergence acceleration in pseudo-time. Convergence was reached when the root mean square of the residual in the entire domain was reduced by at least 3 orders of magnitude for all variables, and the normalised relative difference between two consecutive iterations of lift and drag coefficients, averaged over 100 iterations, was smaller than 10⁻⁶.

Each airfoil was simulated approximately between the positive and negative stall angle with a step of 2^∘. The Spalart–Allmaras (SA) turbulence model was chosen. The discretisation error was assessed with the grid convergence index (GCI) method by Patrick J. Roache (Celik et al., 2008) by simulating each AoA on three different grids. The refinement factor of the grids was $\sqrt{\mathrm{2}}$; i.e. the average area of the elements was doubled each time. Then, flow transition was taken into account by applying the algebraic model by Bas and Cakmakcioglu, SA-BC (Cakmakcioglu et al., 2018) with a freestream turbulence intensity of 0.1 %. Finally, the aerodynamic coefficients were extrapolated to cover the entire range of AoA values (Viterna and Janetzke, 1982) and corrected for 3D effects (Du and Selig, 1998; Eggers et al., 2003) using NREL's tool AirfoilPrep.

Roughness was included in CFD simulations by applying the Boeing extension for the SA turbulence model. The rough-wall model was developed by Spalart (Aupoix and Spalart, 2003) and was recently implemented in SU2 (Ravishankara et al., 2020). It is a modification of the SA model to account for wall roughness so that the logarithmic law of the wall for rough walls (Eq. 6) is correctly represented. It was applied during and after ice accretion. We must point out that this is a low-Re model; i.e. the flow field is computed numerically down to the viscous sublayer, located at $y^{+}<\mathrm{5}$, using a turbulence model. This requires $y_{\mathrm{wall}}^{+}<\mathrm{1}$ to correctly capture the different boundary layer regions. High-Re models, i.e. models using wall functions, would require $y_{\mathrm{wall}}^{+}$ belonging to the log-law region; i.e. $y_{\mathrm{wall}}^{+}\approx \mathrm{50}$. Simulations with rough-wall functions on iced airfoils can be found in a recent work by Yassin et al. (2021).

An unstructured hybrid mesh was used to discretise the domain. It was generated using the code uhMesh (Dussin et al., 2009). The circular domain was made of an O grid of quadrangular elements around the airfoil surrounded by an unstructured grid of triangular elements. uhMesh generates the boundary layer grid with an advancing-front technique with the possibility of local insertion of triangles, while the outer grid is created by computing a Delaunay triangulation using the Bowyer–Watson algorithm. A first cell height of 10⁻⁶ c (chord) ensured that $y_{\mathrm{wall}}^{+}$ was lower than 1 on the entire airfoil in every simulation. A far-field distance of 240 c ensured the independence of the solution from this parameter. On the finest grid, the characteristic lengths applied were 8 c at the far field, 0.3 $c/\mathrm{1000}$ near the leading edge, $c/\mathrm{1000}$ near the trailing edge, and $c/\mathrm{100}$ elsewhere. A close-up view of the fine grid of the NACA 64₃-618 airfoil is shown in Fig. 2.

Figure 2 Fine grid (NACA 64₃-618 airfoil).

During the icing event, we assumed that the presence of roughness due to ice in the stagnation point region caused the transition to a fully turbulent flow. Thus, flow transition was neglected from the beginning of the icing event. This modelling hypothesis was necessary since no roughness-induced transition model is currently available in SU2. The fully turbulent hypothesis may affect two results, i.e. the ice accretion simulations and the aerodynamic coefficients of the iced airfoils.

Regarding flow transition on the iced airfoil, it is interesting to analyse the results of the rough-flat-plate experiment by Feindt (1957). The experiment is commonly used to verify new roughness-induced transition models with transport equations (Dassler et al., 2010; Langel et al., 2014; Min and Yee, 2021) and their boundary conditions (Son and Kim, 2022). First, let us introduce Feindt's k_s Reynolds number, $R{e}_{k_{\mathrm{s}}}=\frac{\mathit{\rho }{U}_{\mathrm{\infty }}{k}_{\mathrm{s}}}{\mathit{\mu }}$. From the experiments, the critical k_s Reynolds number after which the roughness affects transition is $R{e}_{k_{\mathrm{s},\mathrm{cr}}}=\mathrm{130}$. Moreover, as $R{e}_{k_{\mathrm{s}}}$ increases, the transition point moves upstream and the width of the transition region decreases. For $R{e}_{k_{\mathrm{s}}}\gtrsim \mathrm{300}$, the transition point is located at $R{e}_{x_{\mathrm{t}}}=\frac{\mathit{\rho }{U}_{\mathrm{\infty }}{x}_{\mathrm{t}}}{\mathit{\mu }}<\mathrm{0.1}$. On the wind turbine under analysis, considering the outer half of the blade, a rotating velocity of 11 rpm, a minimum k_s of $\mathrm{0.3}\times {\mathrm{10}}^{-\mathrm{3}}$, and a maximum k_s of $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{3}}$, $R{e}_{k_{\mathrm{s}}}$ varies from a minimum value of 750 at mid-span to a maximum value of 15 000 at the tip. Thus, the fully turbulent approximation is acceptable for the computation of aerodynamic coefficients.

Regarding ice accretion simulations, Min and Yee (2021) have recently included the effect of the roughness-induced transition in ice accretion simulations and compared the results with fully turbulent rough simulations. The transition model improved the accuracy in glaze ice simulations, while the rime case was unaffected. Indeed, water droplets freeze upon impact with rime ice, and the numerical solution mainly depends on the collection efficiency. Thus, the fully turbulent hypothesis is also appropriate to study a rime ice accretion.

2.3 The ice accretion problem

The problem of ice accretion is clearly unsteady. As ice grows on the surface, the shape of the airfoil changes, modifying the flow field, the droplet trajectories, and the ice shape as well. During this process, two timescales can be identified. One is related to the growth of ice, and the other is related to the modification of the flow field due to ice growth. The former is, in general, much larger than the latter. For this reason, a quasi-steady, multi-step approach must be adopted. The total accretion time is divided into smaller intervals. In each sub-interval, the flow field and droplets' trajectory are kept constant, and an ice accretion step is performed. The interaction between the gas and the liquid (droplet) phase can be taken into account by using an Eulerian two-fluid model (Re and Abgrall, 2020; Sotomayor-Zakharov and Bansmer, 2021). Otherwise, a one-way coupling approach can be used by computing the particle trajectories with a Lagrangian approach on the previously computed flow field. After the small ice growth, the geometry is updated and the loop is repeated until the final time is reached. For the geometry update, besides simple re-meshing, conservative mesh adaptation techniques (Cirrottola et al., 2021; Colombo and Re, 2022; Donizetti et al., 2023) or mesh-less immersed boundary methods (Lavoie et al., 2022) can be used to avoid failures in grid generation.

Each of these tasks was performed by different software. Once more, SU2 was used for the computation of the flow field. The Lagrangian particle tracking PoliDrop was used to compute the trajectories of the water droplets and the resulting collection efficiency β on the airfoil surface. The ice accretion engine PoliMIce computed the local ice thickness by solving a simplified Stefan problem. Finally, uhMesh was used to generate the grid of the new geometry. No smoothing was applied to the geometries unless grid generation failed. Roughness was applied where ice was predicted using Wright's relation (Eq. 11).

In PoliDrop, we used an iterative method to compute the collection efficiency up to an arbitrary precision. The seeding region was updated at each iteration by adding new particles where needed. A uniform seeding front was initialised as a linear grid with equally spaced elements. At the first iteration, the parcels not hitting the airfoil were identified and removed so that the seeding front was reduced in size. The first two parcels flying just above and below the object were not removed so that the impingement limits were refined as well. Then, at each iteration, elements were incrementally split, evolving the current cloud front and computing the collection efficiency β on the target surface. The simulation stopped when the difference in the L₂ norm between two consecutive iterations of computations of β was below a specified threshold (tol):

$$\begin{array}{}\text{(12)}& \Vert \mathrm{err}\mathit{\beta }{\Vert }_{\mathrm{2}}={\left(\sum _{j=\mathrm{1}}^n{\left[\left({\mathit{\beta }}_j^{[k-\mathrm{1}]}-{\mathit{\beta }}_j^{\left[k\right]}\right)\mathrm{\Delta }{s}_j\right]}^{\mathrm{2}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}<\mathrm{tol},\end{array}$$

where the index k identifies the iteration number, j is the cell on the airfoil, and Δs_j is the size of cell j. At each iteration, the number of parcels doubled, $\Vert \mathrm{err}\mathit{\beta }{\Vert }_{\mathrm{2}}$ halved, and the time to complete the iteration doubled as well. To reduce $\Vert \mathrm{err}\mathit{\beta }{\Vert }_{\mathrm{2}}$ by 1 order of magnitude, approximately three to four iterations were required, and the computational time increased by a factor of 8 or 16, respectively.

Figure 3 Blade discretisation and sections chosen for ice accretion.

It is clear that the choice of the number of time steps in a multi-step ice accretion simulation and the accuracy of the Lagrangian particle tracking are crucial to obtain an accurate solution efficiently. A proper combination of these parameters is required. These were chosen by comparing the numerical solution with three experimental test cases by Han et al. (2012) at the Adverse Environment Rotor Test Stand (AERTS) facility of Pennsylvania State University. The set-up was then tested on six additional test cases from the same series of experiments. Results are shown in Sect. 3.

Once a satisfactory set-up for icing simulations was found, the icing event on the full blade was simulated. Icing was monitored on five independent sections, as shown in Fig. 3. Each section was located at the midpoint of two BEM nodes and is representative of the ice that, on average, is accreted on the two nodes. Thus, local boundary conditions of each iced section (i.e. the relative velocity V_rel and the local AoA α) were computed as the mean value of the two adjacent nodes. Then, the aerodynamic coefficients found on the iced section were applied to the two adjacent nodes. In this way, the entire outer half of the blade was covered. Blade flexibility, wind shear, and tilt and cone angles of the rotor all contributed to generating a periodic output, with a period corresponding to that of a blade revolution, i.e. with a strong 1P component. Since a steady wind was considered, oscillations were limited. In particular, $\mathrm{\Delta }\mathit{\alpha }<\pm \mathrm{0.7}$^∘ and $\mathrm{\Delta }{V}_{\mathrm{rel}}<\pm \mathrm{1}$ m s⁻¹. Thus, the mean values of α and V_rel over one period were computed to consider a steady input for ice accretion. In OpenFAST, fully turbulent aerodynamic coefficients were used as input at this stage to include the effect of early transition due to icing.

The quasi-steady approximation was applied independently to each section, using different time steps according to the local ice accretion rate. The specific time steps used for each section are presented in Sect. 4.1. A matching time of 30 min was chosen to check if it was necessary to update the local boundary conditions due to a change in the equilibrium condition of the wind turbine, similarly to what was done by Zanon et al. (2018). An empirical relation was retrieved using OpenFAST to check the estimated variation in the AoA as ice was accreting on the blade sections. However, the difference in α and V_rel was always negligible during ice accretion.

2.4 Aerodynamics of the iced blade

Due to the high uncertainty in roughness height estimation, two values for $k_{\mathrm{s}}/c$ were considered when computing the aerodynamic coefficients of the iced sections after the icing event. The first one was estimated with Wright's formula (Eq. 11), which generally prescribes $k_{\mathrm{s}}/c=\mathrm{0.34}\times {\mathrm{10}}^{-\mathrm{3}}$ for rime ice. This roughness height is identified in the text with the letter W. Then, this value was increased by 1 order of magnitude to $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{3}}$. Since this value is close to the one prescribed by Shin's relation corrected for drag prediction, we refer to this case with the letter S.

It will be shown in Sect. 3.2 that the impingement limits are slightly under-predicted by PoliMIce and other codes during a steady ice accretion. Moreover, blade vibrations and the highly unsteady incoming wind are likely to increase the wet region of the blade surface during real wind turbine operation. Thus, in numerical studies of wind turbine icing, applying roughness on a greater area than the ice shape only is standard practice. For instance, Homola et al. (2012) used roughness on the entire surface of each airfoil, while Etemaddar et al. (2014) applied roughness to the ice shape and the first 25 % of the airfoil chord.

To quantify the effect of this modelling choice on airfoil performance and power losses, besides using two different roughness heights, we also considered two regions to which roughness was applied. In the first case, it was applied to the ice shape only. This case is denoted in the text by adding the subscript std (standard) to the letter defining the roughness height: W_std ($k_{\mathrm{s}}/c=\mathrm{0.34}\times {\mathrm{10}}^{-\mathrm{3}}$) and S_std ($k_{\mathrm{s}}/c=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{3}}$). The second roughness region included the first one and extended beyond it, covering a length corresponding to 25 % of the chord of Section A, i.e. 0.44 m, on all sections and on both the upper and the lower surface. On the other sections, from B to E, this corresponded to 18 %, 15 %, 13 %, and 11 % of the chord, respectively. This case is denoted in the text by applying the subscript ext (extended) to the letter defining the roughness height: W_ext and S_ext. The two regions are shown in Fig. 4 on Section B. An overview of the four test cases is provided in Table 3. The value of 25 % for Section A was picked to match the value chosen by Etemaddar et al. (2014). However, a different modelling choice was made, and this value was kept constant dimensionally (0.44 m) among all sections. By keeping fixed the dimensional width of the ext region, its non-dimensional width reduced as the chord of the airfoil increased, being more consistent with the physics of the problem. Indeed, the greater the chord, the greater the pressure gradient generated by the section. Thus, water droplets are deflected away, and the wet area on the section reduces, at least in non-dimensional terms.

Figure 4 Definition of the std and of the ext roughness regions on an iced airfoil. The std region corresponds to the numerically computed ice shape. The ext region includes the std region and extends for 0.44 m beyond the impingement limits of the ice shape, regardless of the chord of the airfoil. On Section B, this corresponds to 18 % of the chord.

3 Validation of the numerical set-up

3.1 CFD solver

To validate the set-up of the CFD solver, the aerodynamic coefficients of the clean airfoils of the blade were computed and compared with experimental data. The aerodynamic coefficients of DU airfoils were measured by Ruud van Rooij of the Delft University of Technology at a Reynolds number of 7×10⁶. NACA 64₃-618 coefficients were taken from Abbott et al. (1945) at a Reynolds number of 6×10⁶. All airfoil data are provided in the DOWEC report (Kooijman et al., 2003). The comparison between numerical simulations and experimental results is shown in Figs. 5–10. A correction of −0.4^∘ was applied to experimental data of the NA18 airfoil as suggested by Timmer (2009) due to a possible error in the orientation of the model in the wind tunnel. The moment coefficient was computed with respect to $\frac{c}{\mathrm{4}}$ and is positive for nose-up rotations.

Table 3 Roughness heights (a) and regions (b) tested, as well as the resulting test matrix (c).

Figure 5 Aerodynamic coefficients of the NA18 airfoil.

Figure 6 Aerodynamic coefficients of the DU21 airfoil.

Figure 7 Aerodynamic coefficients of the DU25 airfoil.

Figure 8 Aerodynamic coefficients of the DU30 airfoil.

Figure 9 Aerodynamic coefficients of the DU35 airfoil.

Figure 10 Aerodynamic coefficients of the DU40 airfoil.

First, we analyse the results for fully turbulent flows. All fully turbulent simulations showed satisfactory grid convergence for almost every aerodynamic coefficient computed. The estimate of the discretisation error was computed with the GCI method (Celik et al., 2008) and was represented through error bars. In the attached flow regime, the lift coefficient was underestimated for all airfoils, while the drag coefficient was overestimated. The error with respect to experimental data increased, together with the relative thickness of the airfoils. It was maximum at the root of the blade. The positive stall angle and lift coefficient were over-predicted for $t/c\le \mathrm{30}$ %. The maximum lift coefficient became under-predicted for $t/c\ge \mathrm{35}$ %, with the error increasing for increasing $t/c$. At negative stall, good predictions were made up to $t/c\le \mathrm{21}$ %, while minimum lift coefficients were underestimated (in absolute values) for $t/c\ge \mathrm{30}$ %.

When the algebraic BC transition model was included in the system of equations, the aerodynamic coefficients were accurately predicted for all airfoils for attached flows, regardless of their relative thickness. The absolute value of the maximum lift coefficient increased, at both positive and negative stall. This led to more accurate predictions of positive stall for $t/c\ge \mathrm{35}$ % and of negative stall for $t/c\ge \mathrm{25}$ %.

The Boeing extension for the SA model was then tested to compare the numerical results with the law of the wall for rough surfaces. The relation is presented in Eq. (6). Two cases were considered. In the first case, a roughness height $k_{\mathrm{s}}/c=\mathrm{0.5}\times {\mathrm{10}}^{-\mathrm{3}}$ was applied on the entire surface of the NA18 airfoil. In the second case, $k_{\mathrm{s}}/c$ was set to $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{3}}$. Each value of $k_{\mathrm{s}}/c$ leads to different $k_{\mathrm{s}}^{+}$ values. This occurs since the skin friction varies locally and so does the viscous length scale δ_ν. In these simulations, the Reynolds number was 6.6×10⁶ and the AoA was 0^∘. The results are shown in Fig. 11. The velocity profile in wall units u⁺ is shown as a function of the non-dimensional wall distance y⁺ at different stations along the airfoil and compared to the theoretical results obtained with the local $k_{\mathrm{s}}^{+}$. In the top row, results are for $k_{\mathrm{s}}/c=\mathrm{0.5}\times {\mathrm{10}}^{-\mathrm{3}}$ on the suction side of the airfoil. In the bottom row, $k_{\mathrm{s}}/c=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{3}}$ and results are extracted from the pressure side. The numerical solution for the smooth airfoil is also shown and compared with the theoretical behaviour (Eq. 10). For the values of $k_{\mathrm{s}}/c$ under analysis, all the resulting $k_{\mathrm{s}}^{+}$ values belonged to the fully rough regime, typical of ice. The model accurately captured the different shifts in the logarithmic region of the law of the wall.

Figure 11 Law of the wall using Boeing extension for the SA turbulence model. NA18 airfoil, α=0^∘, $Re=\mathrm{6.6}\times {\mathrm{10}}^{\mathrm{6}}$. Top row: $k_{\mathrm{s}}/c=\mathrm{0.0005}$. Bottom row: $k_{\mathrm{s}}/c=\mathrm{0.005}$.

3.2 Ice accretion simulations

Two different approaches were tested for ice accretion. These were almost equivalent in overall computational time. In the first one, the collection efficiency β of the droplets was finely computed during each time step, setting $\Vert \mathrm{err}\mathit{\beta }{\Vert }_{\mathrm{2}}<\mathrm{3}\times {\mathrm{10}}^{-\mathrm{6}}$ (Eq. 12). In the second one, the tolerance was set to $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{5}}$ and the number of icing steps was increased. Numerical results were compared with experimental rime ice accretion on a rotating blade section, consisting of an S809 airfoil with c=0.267 m. Experiments were carried out by Han et al. (2012). AERTS test cases no. 20, no. 21, and no. 22 were chosen to find the best computational set-up. Then, the chosen set-up was tested on six additional cases, i.e. AERTS no. 4 and nos. 15–19. Test cases nos. 20–22 were chosen as the primary benchmark since they are the longest available and test conditions are the most similar to those of the icing event under analysis. The three test cases only differ in the duration of the icing event, which was 30, 60, and 90 min, respectively. Test conditions are reported in Table 4.

Table 4 Test conditions of AERTS test cases.

A time step of 15 min was used when $\Vert \mathrm{err}\mathit{\beta }{\Vert }_{\mathrm{2}}=\mathrm{3}\cdot {\mathrm{10}}^{-\mathrm{6}}$, while 3 min was chosen when $\Vert \mathrm{err}\mathit{\beta }{\Vert }_{\mathrm{2}}=\mathrm{3}\cdot {\mathrm{10}}^{-\mathrm{5}}$ to try to match the computational time. Results are reported in Fig. 12. In both cases, the ice impingement limit on the lower surface was underestimated. This is a common issue in numerical ice accretion simulations. A real cloud is made by a distribution of droplet diameters, and the MVD is just an indicator of the median of this distribution. Parcels with higher diameters have a higher mass, and the pressure gradient deflects their trajectory less. Thus, a wider portion of the airfoil gets wet. This phenomenon may be overcome in numerical simulations with a multi-bin approach, i.e. by performing the weighted average of the collection efficiency computed with different droplet diameters from the diameter distribution (Sirianni et al., 2022). The effect of the uncertainty in other operating conditions was analysed in a recent work by Gori et al. (2022). On the other hand, by using a finer time discretisation, a more accurate solution at the leading edge was obtained. A higher number of ice layers led to a better representation of the physics of the problem while limiting the propagation of errors from one step to the other. This permitted the reduction in the accuracy in the computation of β without losing accuracy in the computation of the solution.

Figure 12 PoliMIce simulations of three AERTS test cases. Each column represents the same test case. The three test cases have the same atmospheric conditions and only differ in total ice accretion time (left: 30 min; centre: 60 min; right: 90 min). In each row the same numerical set-up is used for the multi-step ice accretion (top: Δt=15 min, $\Vert \mathrm{err}\mathit{\beta }{\Vert }_{\mathrm{2}}<\mathrm{3}\times {\mathrm{10}}^{-\mathrm{6}}$; bottom: Δt=3 min, $\Vert \mathrm{err}\mathit{\beta }{\Vert }_{\mathrm{2}}<\mathrm{3}\times {\mathrm{10}}^{-\mathrm{5}}$).

Figure 13 Comparison between PoliMIce simulations, LEWICE simulations, and experiments of AERTS test cases. $\Vert \mathrm{err}\mathit{\beta }{\Vert }_{\mathrm{2}}<\mathrm{3}\times {\mathrm{10}}^{-\mathrm{5}}$.

Moreover, a noticeable reduction in the elapsed real time for the entire 90 min simulation was found (approx. 13 %). For these reasons, the approach consisting of a high number of time steps with lower accuracy for β was chosen.

The set-up with $\Vert \mathrm{err}\mathit{\beta }{\Vert }_{\mathrm{2}}=\mathrm{3}\cdot {\mathrm{10}}^{-\mathrm{5}}$ was tested on the remaining AERTS test cases. The time step was adjusted according to the LWC. For test cases no. 18 and no. 19, Δt=180 s was used once more. For test cases no. 4 and nos. 15–17, Δt=112.5 s was chosen to match the accumulation parameter LWC⋅Δt of all the previous simulations. Results are shown in Fig. 13. The results obtained with the selected set-up show an almost perfect match with LEWICE simulations. Good agreement is found with the experiments in terms of ice thickness and impingement limits on the upper surface. On the lower surface, impingement limits are at times underestimated.

4 Results and discussion

4.1 Blade icing

The local boundary conditions computed at the beginning of the icing event are reported in Table 5, together with the time step chosen for each section.

Table 5 Local boundary conditions of the five sections under analysis.

There was no need to update the boundary conditions during ice accretion. For instance, the AoA of Section B increased by 0.35^∘ after the icing event due to the degradation of the aerodynamic performances of the wind turbine. However, this may not hold if greater roughness height and extension values were considered during ice accretion. The computed ice shapes are shown in Fig. 14 in non-dimensional form, while a detailed view of the multi-step process on Section B is shown in Fig. 15.

Figure 14 Non-dimensional comparison of the ice shapes on Sections A–E.

Figure 15 Multi-step ice accretion on Section B. The grid for the computation of the aerodynamic coefficients is superimposed onto the computed ice shape, highlighting the removal of highly concave regions. The region within the red box is enlarged in Fig. 16.

Figure 16 Detail of the computational grid on the final ice shape of Section B, showing the orthogonality of the grid. The region corresponds to the red box shown in Fig. 15.

The ice shapes on Sections A, B, and C (i.e. NA18 sections) were very similar. Their main difference was the length of the horn, which decreased towards the root of the blade. Some small secondary protrusions were formed on the main ice shape. These are due to some small oscillations of the collection efficiency, which eventually got amplified step after step because the geometry was not smoothed unless strictly required by the grid generator. Section E was almost unaltered. On this section, 0.42 kg m⁻¹ of ice was found. The ice mass accreted on the blade increased almost linearly, up to 3.35 kg m⁻¹ on Section A. The total accreted mass was estimated to be lower than 100 kg, i.e. less than 0.5 % of the total mass of the blade. Thus, it was chosen to neglect the additional mass during the aeroelastic simulations, although its distribution may have altered the modal response of the wind turbine.

4.2 Iced blade aerodynamics

The aerodynamic coefficients of the iced sections were then computed in the four cases defined in Sect. 2.4, i.e. W_std, S_std, W_ext, and S_ext (see Table 3 and Fig. 4). Before the computation, the size of the concave regions around the secondary protrusions in Sections A and B was reduced to improve grid quality while maintaining the same overall ice shape. The result of this process is shown in Fig. 15 for Section B by superimposing the computational grid used for the aerodynamic coefficients onto the ice shape coming from the ice accretion simulations. A red box identifies the most critical region for grid generation, which is shown in detail in Fig. 16, together with the orthogonality of the grid. uhMesh guaranteed good-quality meshing by generating a hybrid boundary layer grid, introducing triangles in highly concave and convex regions.

Figure 17 Aerodynamic coefficients of Section A.

Figure 18 Aerodynamic coefficients of Section B.

Figure 19 Aerodynamic coefficients of Section C.

Results of the CFD simulations are shown in Figs. 17–21, considering the effect of roughness height and extension. A quantitative comparison is provided in Table 6, where the percentage variation in the aerodynamic coefficients with respect to the clean case is reported for each section and each roughness at α=4^∘. The aerodynamic coefficients were non-dimensionalised with respect to the clean airfoil chord. The moment coefficient was computed with respect to the same point of the clean airfoil ($c_{\mathrm{clean}}/\mathrm{4}$). In all cases, the presence of ice caused a degradation of the aerodynamic performances due to both the ice shape and roughness. As expected, stall was anticipated, the slope of the lift coefficient decreased, the drag coefficient increased, and the moment coefficient changed significantly. The greatest difference was found when a higher roughness was applied to a wider portion of the airfoils (S_ext). The case with smaller roughness on the same region followed (W_ext). When roughness was applied where ice was predicted (W_std and S_std), the results were similar. However, the behaviour of each section was different. We may analyse the results by distinguishing between the effects of ice shapes and roughness.

Figure 20 Aerodynamic coefficients of Section D.

Figure 21 Aerodynamic coefficients of Section E.

Figure 22 Lift-to-drag ratio of Section B considering different cases. Lines represent airfoil efficiency in various cases. Shaded areas represent different contributions to loss in efficiency, providing a qualitative superposition of effects. (a) $k_{\mathrm{s}}/c=\mathrm{0.34}\times {\mathrm{10}}^{-\mathrm{3}}$. (b) $k_{\mathrm{s}}/c=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{3}}$.

Table 6 Aerodynamic penalties on Sections A–E at α=4^∘.

We start considering the W_std and S_std cases. Given the decreasingly big ice horns and the small difference between the two cases on Sections A, B, and C, we can conclude that the ice shape was mainly responsible for the aerodynamic penalty on NA18 sections. This can be clearly seen in Fig. 22, where the results for the ice shape without roughness were included. In the figure, the lift-to-drag ratio of Section B is represented as a function of the AoA. The efficiency of the section considering a smooth ice shape was almost coincident with the W_std case, and only a slight decrease was found with a higher roughness on the ice shape (S_std). For completeness, the results of the clean, smooth, fully turbulent airfoil were included to qualitatively highlight the effect of icing at the beginning of ice accretion, when the ice shape is negligible and the transition occurs earlier and earlier due to increased roughness. On the other hand, on DU sections the ice shape was small and so was the region where roughness was applied. For attached flows, results were almost identical to those of the respective clean airfoils with a fully turbulent flow. Moreover, in this range of AoA values, the performance degradation on the thicker Section E was slightly higher compared to the thinner Section D. This is coherent with the results of the fully turbulent clean airfoils. For these reasons, we can conclude that the difference between the clean and the two iced cases in this flow regime was due neither to roughness nor to the ice shape. In fact, it was due to the early transition caused by the presence of roughness, which was modelled with a fully turbulent flow. In reality, roughness height affects the transition, but it was shown at the end of Sect. 2.2 that the assumption of a fully turbulent flow is reasonable for such a long ice accretion. On these same sections, positive stall was almost unaffected compared to the fully turbulent solution. Negative stall occurred earlier, in particular on Section D. This effect was related to the small, downward-pointing ice shape. Given these results, we can conclude that the effect of ice shape becomes predominant over roughness as the horn grows in size, in accordance with previous studies (Battisti, 2015).

We now consider the two cases of extended roughness: W_ext and S_ext. For these two cases, the results were different from each other and were also different from the std cases in almost every simulation. We highlight once more that the region of extended roughness was equal in size (0.44 m) on all sections, and so it increased in terms of non-dimensional airfoil length from root to tip. By looking at NA18 sections, at the lowest AoA values, the aerodynamic coefficients coincided in all cases. Thus, at these AoA values, the aerodynamic penalties were still produced by the ice shapes. Negative stall occurred at $\mathit{\alpha }=-\mathrm{8}$^∘ when the separated regions generated from the leading and the trailing edge merged. As the AoA increased, the effect of the roughness extension became more and more important, together with the value of k_s. The slope of the C_L(α) curves decreased, while drag and moment coefficients increased. This effect is peculiar since roughness should have little effect on the aerodynamic coefficients when ice horns are well developed. The extended roughness region caused a high increase in skin friction in a geometrically smooth region of the section, increasing the viscous drag. Moreover, the flow expanded less on the suction side of the sections and was compressed less on the pressure side, causing a noticeable reduction in lift. The differences were much higher when roughness was increased by 1 order of magnitude (i.e. for case S_□). The difference between the std and the ext cases increased towards the tip of the blade since roughness was applied on a wider portion of the airfoil. On DU sections, however, the opposite occurred. For both the roughness heights tested, Section E was more sensitive to roughness than Section D in both attached flow and stall conditions. This occurred despite roughness covering a slightly shorter portion of the innermost section. Previously, it was shown that the ice shape only affects the negative stall of Section D. In the other flow regimes of Section D and on all flow regimes of Section E, we may think to have a fully turbulent airfoil, i.e. an airfoil where the transition is fixed at the stagnation point. For a transition-fixed flow, Somers (2005) found that the detrimental effect of leading-edge roughness increases with the relative thickness of the airfoil, as occurred in this case.

4.3 Effect of icing on power production

In the previous section, it was shown that the differences between W_std and S_std cases are negligible. Moreover, only small differences are found with the W_ext case. Thus, in this section, only the lower-roughness, tighter-impingement case (W_std) and the higher-roughness, wider-impingement case (S_ext) are compared with the clean case (SA-BC).

The C_P–TSR curves (Eq. 2) were computed for a pitch angle of 0^∘ using the aerodynamic module of OpenFAST, AeroDyn. Results are shown in Fig. 23. In the iced cases, the C_P values were lower for any TSR. As expected, the lowest values were found in the S_ext case. The highest decrease in C_P occurred at TSR>7. These values are used at low wind speeds when the wind turbine operates in Region 1.5. In particular, from the cut-in wind speed up to 8 m s⁻¹, the TSR decreases from 15.3 to its optimum value of ∼7. The power coefficient became negative for TSR values between 12 and 13. It is worth noticing how the TSR corresponding to the maximum C_P changed from case to case. It was approximately 7 for the clean case, while it decreased to 6.5 for W_std and increased to 7.1 for S_ext.

Figure 23 C_P–TSR curves with pitch angle β=0^∘.

Then, the power curves were computed with a steady inflow. They are shown in Fig. 24. Power losses are shown in Fig. 25 as both absolute and normalised differences with respect to the clean case. With ice, power production started at 4 m s⁻¹. The normalised power loss was maximum at cut-in wind speed and diminished as the TSR decreased from the start-up value of approximately 15 to a constant value in Region 2. In this region, the power loss is approximately 6 % for W_std and 9 % for S_ext. The TSR value obtained through the generator torque controller in Region 2 was approximately 7.4 in the clean case and 7.2 in the iced cases. This means that an almost optimum TSR was used for S_ext, while a sub-optimum one was used for both W_std and the clean case. By regulating the generator torque, their power output may increase.

Next, the power curves were computed with the turbulent inflow prescribed by the IEC in DLC 1.1. They are shown in Fig. 26, while power losses are shown in Fig. 27. The first clear effect of the increased turbulence intensity is the inflexion of the power curve close to the rated speed. On the other hand, the non-constant, non-uniform wind speed made the wind turbine produce slightly more power at low wind speeds. Regarding power losses, with a mean wind speed of 3 m s⁻¹ a slightly higher power was produced with ice with respect to the clean case. This result may differ if a different random seed was used to generate the realisation of the turbulent wind used as input. At higher wind speeds, the trend was similar to that of a steady inflow. However, due to the variability of wind, there was no clear distinction between the different controller regions, and the results of steady wind were smoothed out. In general, higher power losses were found at all wind speeds, except for the nominal rated speed (11 m s⁻¹). Power losses in turbulent wind reduced and became null at approximately 15 m s⁻¹. The actual rated speed remained unchanged after the icing event at 17 m s⁻¹.

Figure 24 Power curve with a steady inflow.

Figure 25 Power losses with a steady inflow. (a) Absolute difference. (b) Normalised difference.

Figure 26 Power curve with a turbulent inflow.

Figure 27 Power losses with a turbulent inflow. (a) Absolute difference. (b) Normalised difference.

As visible from Figs. 25 and 27, the effect of roughness on power production depended on wind speed. At low wind speeds, the power loss was similar in both cases under analysis, while differences increased with wind speed. This trend was aligned with the one found by Etemaddar et al. (2014), while it differed from those found by Homola et al. (2012) and Turkia et al. (2013). In order to give a single figure of the difference between the two roughness cases, the Weibull-averaged power P_W was computed and compared with the clean case. Its value was 2618, 2528, and 2483 kW for the clean, W_std, and S_ext cases, respectively. The average power loss of W_std was 3.44 %. For the S_ext case, it was 5.16 %, which is 50 % higher than W_std. This difference is not negligible, even though the ice shapes were well developed and the region of extended roughness was rather limited. The same quantity was computed for the power curves computed by Homola et al. (2012), Turkia et al. (2013), and Etemaddar et al. (2014) and is reported in Table 7. Once more, our results were consistent with those by Etemaddar et al. (2014), where the icing event lasted 24 h and roughness was applied on 25 % of the chord of the blade. On the other hand, Homola et al. (2012) predicted an average power loss of about 10 % for an icing event of one-third of the duration of the one analysed in the current study but in the same atmospheric conditions. In this case, roughness was applied on the entire blade surface using Shin's relation.

Homola et al. (2012)Turkia et al. (2013)Etemaddar et al. (2014)

Table 7 Weibull-averaged power loss computed for the current study and compared with other authors.

From these results, it is clear that the research on numerical simulations of icing on wind turbines should focus on water impingement limits and roughness height. Regarding the impingement limits, better results may be obtained by considering unsteady ice accretion simulations. However, the detail required for time discretisation is unknown. This is not sufficient since it is not possible to obtain reliable results by using the classical empirical correlations for k_s coming from the aeronautic field. These relations were developed for different systems operating in completely different environments. In situ roughness measurements are required to remove uncertainty of this parameter. Proper numerical predictions would allow an improvement in the design of ice protection systems and wind turbine controllers during icing events.

5 Conclusions

In this paper, we conducted a detailed numerical simulation of ice accretion on the NREL 5 MW wind turbine blade using the BEM approach. To increase the precision in the computation of ice shapes, we proposed to use independent time steps during a multi-step ice accretion simulation. Moreover, we showed that it is possible to reduce the computational time required for ice accretion simulations by increasing the error in the collection efficiency and adding a very small ice thickness during each step.

Then, we analysed the effect of roughness on the aerodynamic performances of the iced sections. Due to the uncertainty of these parameters, we considered two roughness heights and two roughness extensions on each section. We computed the aerodynamic coefficients for each case, and we assessed whether the aerodynamic penalty was due to ice, roughness, or both. It was shown that roughness can significantly affect the aerodynamics of an iced section, even when a complex ice shape is present, as long as k_s is sufficiently high.

Finally, we computed the power curves for the low-roughness (W_std) and the high-roughness (S_ext) cases and compared them with the results of the clean wind turbine. We computed a Weibull-averaged power for each case to introduce a single figure indicating the severity of the icing event. The power loss was 50 % higher for the high-roughness case.

This high variability in the prediction of power losses suggests two main areas of research for future work. The first one should be focused on the correct detection of the impingement limits of water droplets in the highly unsteady environment in which wind turbines work. The second one should be focused on the characterisation of roughness distribution and height on real wind turbine blades.