Aerodynamic interaction of rain and wind turbine blades: the significance of droplet slowdown and deformation for leading-edge erosion

Leading-edge rain erosion is a severe problem in the wind energy community since it leads to blade damage and a reduction in annual energy production by up to a few percent. The impact speed of rain droplets is a key driver of the erosion rate; therefore, its precise computation is essential. This study investigates the aerodynamic interaction of rain droplets and wind turbine blades. Based on findings from the literature and an analysis of the relevant parameter space, it is found that the aerodynamic interaction leads to a reduction in the impact speed. Additionally, the rain droplets deform and break up as they approach the wind turbine blade. An existing Lagrangian particle model, developed for research in aircraft icing, is adapted, extended, and validated for leading-edge rain erosion to study the process in more detail. Results show that the droplet slowdown reduces predicted damage toward the tip of the blade by over 50 %. The model indicates that the aerodynamic blade interaction affects small droplets significantly more than large droplets. Due to this drop-size dependency, the damage accumulation is shifted toward higher-rain-intensity events. Additionally, the droplet impact speed is sensitive to the aerodynamic nose radius of the airfoil. Due to this sensitivity and its drop-size dependency, the slowdown effect provides interesting levers for erosion mitigation via blade design or operational adjustments. To conclude, the aerodynamic interaction between droplet and blade is non-negligible and needs to be taken into account in erosion lifetime models.

1 Introduction

Leading-edge rain erosion is a severe problem for wind turbines. During precipitation events, hydrometeors impact the blade and, over time, create roughening. This damage can grow to large pits that can reach deep into the different structural layers of a wind turbine blade. The damage must be repaired periodically to prevent wind turbine blade failure. Roughening also disturbs the airflow over the blade and leads to a loss in annual energy production (AEP) (Barfknecht et al., 2022). Forecasting and understanding the mechanisms of erosion are important not only for maintenance, but also for operational adjustments of the turbine, such as the erosion-safe mode (Bech et al., 2018).

A key parameter that influences the rain erosion lifetime is the impact speed of the rain droplets. A common approach is to relate the droplet impact speed V_impact via a power law to an incubation metric N, for example, the incubation time, i.e., the operational time until visible erosion damage occurs:

$$\begin{array}{}\text{(1)}& N\propto {\displaystyle \frac{\mathrm{1}}{V_{\mathrm{impact}}^{\mathit{\beta }}}},\end{array}$$

where β is a constant. The variable N can have various meanings depending on the damage model, such as the number of impacts or the impingement. Common to all models is that the magnitude of the parameter β is significant. Parameters for β reported in the literature are 5.7 in Hoksbergen et al. (2022), 16.92 in Shankar Verma et al. (2021), and 7.2–10.5 in Bech et al. (2022). While the reported values in the literature differ significantly based on the test apparatus used and exact definition of N, they all preserve the character of the equation, namely that small changes in V_impact will yield vastly different erosion lifetimes. It is, therefore, important to accurately determine the impact speed. The impact speed is the surface-normal component of the impact vector, which is calculated as the difference between the blade section and rain droplet velocity vector; i.e.,

$$\begin{array}{}\text{(2)}& V_{\mathrm{impact}}=\left({\mathit{V}}_{\sec }-{\mathit{V}}_{\mathrm{rain}}\right)\cdot {\mathit{n}}_{\mathrm{LE}},\end{array}$$

where n_LE is the surface-normal component of the leading edge. The droplet's velocity vector is usually considered to be comprised of components of the droplet's terminal velocity and its advection velocity with the wind (Barfknecht et al., 2022; Verma et al., 2020). However, in aeronautics, it has already been known for a long time that rain droplets and wings can interact aerodynamically (Nicholson, 1968). This leads to rain droplet deformation and slowdown when observed from the wing (Vargas and Feo, 2011), thus adding an extra velocity component to the problem. The potential slowdown of rain droplets has so far received no attention in the wind energy community. One exception is the work of Prieto and Karlsson (2021), although only limited results for spherical droplets were obtained. No droplet deformation was included in their analysis.

Sor et al. (2019) performed measurements in which water droplets were seeded in a rotating-arm test rig. A blunt airfoil was mounted on the arm. High-speed photographs were taken that showed the droplets close to impact with the wing. Figure 1 shows part of their results. As the airfoil approaches the droplets, they start to deform from a spherical to an oblate shape. Shortly before impact, the droplets undergo breakup. While the experiments were performed for aircraft icing research, the parameter space fits the one encountered in leading-edge erosion of wind turbines very well. These findings stand contrary to current practice in leading-edge erosion research, where it is assumed that the droplets are spherical on impact (Hoksbergen et al., 2023; Fæster et al., 2021; Keegan et al., 2012; Verma et al., 2020). The measurements of Sor et al. (2019) imply that rain droplets can undergo breakup, and, therefore, the rain droplets' appearance on impact can be considered complex in shape.

Figure 1 High-speed photography of falling water droplets of 1.75 mm diameter approaching an airfoil in a rotating-arm test rig. The different frames show the temporal progression; the shadow that can be observed in the last two frames is the approaching airfoil. The droplets are first round; then become oblate; and, before impact, break up with smaller droplets being emitted from the rim. The free-stream velocity is 60 m s⁻¹ and airfoil chord is 1.05 m. The photographs are reproduced with permission from Sor et al. (2019).

Figure 2 shows the results of a similar experiment performed by Vargas and Feo (2011). It can be observed how the rain droplets' relative speed changes in front of the airfoil. The water droplets slow down as they approach the leading edge of the rotating airfoil. Droplets with a free-stream velocity of 90 m s⁻¹ experience a velocity reduction of almost 12 m s⁻¹. Considering the exponent of the damage law in Eq. (1), this effect is highly relevant. It appears, therefore, that the effect of droplet slowdown and deformation cannot be neglected when studying leading-edge rain erosion and needs to be further understood.

Figure 2 Relative velocity between a 0.49 mm diameter water droplet and an airfoil; relative velocity is given as a function of the distance to the leading edge; blunt airfoil with a chord length of 0.47 m; five different free-stream (f.s.) velocities; data were collected in a rotating-arm test rig and are reproduced from Vargas and Feo (2011).

The research presented here investigates the impact of rain droplet slowdown and deformation on the erosion lifetime prediction of wind turbine blades. It is important to note that this study assumes that the problem is observed in the reference frame of the airfoil. From an airfoil's perspective, the incoming droplet's speed reduces; hence the term slowdown is used. An observer located on the ground will see the droplets gain speed. Since the effect reduces the impact speed, the term slowdown seems appropriate. The paper is organized as follows: Sect. 2 gives a brief introduction to droplet deformation and breakup. Additionally, the parameter space of the problem is investigated. In Sect. 3, an existing droplet model, developed for research in aircraft icing, is adapted, extended, calibrated, and validated to study the slowdown and deformation process. Subsequently, computational fluid dynamics (CFD) simulations are performed on airfoils of reference turbines to obtain their background velocity field to determine parameters needed in the model. Finally, the precipitation data used and the damage model are discussed in further detail. In Sect. 4, our proposed slowdown and deformation model is employed to analyze the sensitivity of the droplet model with respect to the droplet diameter and the airfoil's aerodynamic nose radius. This is followed by combining the model with the precipitation data and then computing the impact of the droplet slowdown and deformation on the lifetime of two reference turbines. Finally, in Sect. 5 a summary is provided, conclusions are drawn, and recommendations are given.

2 Discussion of the underlying physics

An understanding of the physics of the droplet along its trajectory is necessary before a suitable approach can be chosen to model droplet slowdown and deformation. A complete review of the known processes encountered during aerodynamic droplet deformation and breakup is out of the scope of this work. However, since droplet deformation and breakup comprise a rather new phenomenon for the leading-edge erosion community, a brief summary with a discussion of the parameter space seems appropriate.

For aerodynamic droplet deformation and breakup, the important non-dimensional numbers are the Weber number (We) and the Ohnesorge (Oh) number (Jackiw and Ashgriz, 2021). They read

$$\begin{array}{}\text{(3)}& {\displaystyle }& {\displaystyle }\mathit{We}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{air}}{V}_{\mathrm{slip}}^{\mathrm{2}}{\mathit{\varphi }}_{\mathrm{0}}}{{\mathit{\sigma }}_{\mathrm{water}}}},\text{(4)}& {\displaystyle }& {\displaystyle }\mathit{Oh}={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{water}}}{\sqrt{{\mathit{\rho }}_{\mathrm{water}}{\mathit{\sigma }}_{\mathrm{water}}{\mathit{\varphi }}_{\mathrm{0}}}}},\end{array}$$

with density ρ, surface tension σ, and dynamic viscosity μ, where the subscripts “air” and “water” indicate the corresponding medium. ϕ₀ represents the droplet diameter and V_slip is the slip velocity, i.e., the difference between the velocity of the air (V_air; see Sect. 3.2) and the droplet (V_x; see Sect. 3.1). The Weber number relates the inertial forces to the surface tension forces, whereas the Ohnesorge number relates the viscous to the inertial and surface tension forces. Depending on the Weber number, droplets subject to aerodynamic forces can first undergo deformation and subsequently also break up. Figure 3 shows an often-cited graph taken from Hsiang and Faeth (1995). It depicts how droplets are expected to behave depending on the Weber and Ohnesorge number. From the figure, it is evident that for Oh<0.1, the expected behavior is a function of the Weber number only.

Figure 3 Droplet deformation and breakup modes as a function of the Weber and the Ohnesorge numbers; the figure is reproduced with permission from Hsiang and Faeth (1995).

Aerodynamic droplet breakup consists of two phases, the initiation, also called the deformation phase, and the breakup phase (Jackiw and Ashgriz, 2021). During the deformation phase, the droplet's shape flattens. At some point the droplet breaks up into smaller droplets. This process is also called secondary breakup. Different breakup modes exist such as bag, bag and stamen, multimode, and shear breakup. Some of these modes are shown in Table 1. After the breakup stage is complete, the original droplet will have decayed into a series of small droplets that can be characterized by a drop-size distribution. Subsequently, the resultant droplets might deform and break up again, forming a decay cascade. For more information about the fundamental mechanics of droplet dynamics, the reader is recommended to read Jackiw and Ashgriz (2021, 2022).

García-Magariño (2016)

Table 1 Examples of droplet deformation and breakup in the measurement campaigns of García-Magariño (2016); estimated Weber numbers on impact calculated with the model from Sect. 3.1; photographs are reproduced with permission.

Two figures were created to analyze the parameter space for the leading-edge erosion problem in more detail. Figures 4 and 5 show the Ohnesorge and the Weber numbers, respectively, for droplets of varying sizes and free-stream velocities.¹ This set of simulations spans the parameter space in terms of non-dimensional numbers for the erosion problem. The Ohnesorge number is not dependent on the flow velocity and is, apart from the physical constants, a function of the droplet diameter only. Figure 4 shows that the Ohnesorge numbers stay below 0.1, indicating that the droplet breakup is governed by the Weber number only. The Weber numbers lie in a very broad range of 1 to 800. The wide range of Weber numbers encountered in this problem leads to very different droplet behaviors. The droplet behavior is expected to range from simple deformation for small and slow droplets to shear breakup for larger and faster droplets. In Table 1, example images of the different breakup modes in a rotating-arm test rig are given, together with an approximate Weber number close to impact.

Figure 4 Ohnesorge number plotted against the rain droplet diameter; the values of the Ohnesorge number stay below 0.1.

Figure 5 Contour plot of Weber numbers at impact time as a function of droplet diameter and free-stream velocity; the slip velocity required for the Weber number computation was calculated with the model from Sect. 3.1; the model parameters are R_c=0.07 m and n=1.1.

Most fundamental research in the literature about droplets is based on experiments in shock tubes and steady disturbances (Hsiang and Faeth, 1995). However, in the present problem, the droplets traverse a velocity field that changes depending on the distance to the airfoil. Figure 6 gives an example velocity field. Therefore, it is not possible to directly translate the graph of Fig. 3 to, for example, the outcomes in Table 1. For this problem, the shape and extent of the background velocity field must also be considered. It is intuitive to assume that a larger airfoil will have more influence on the behavior of the droplet than a small airfoil, even though the Weber number of the droplet is similar for both airfoils close to impact. Therefore, one needs to conclude that, while the general body of droplet breakup and deformation is extensive, only very limited knowledge exists that is applicable to the wind turbine rain erosion problem.

Figure 6 Non-dimensional velocity field along the stagnation streamline vs. the dimensional distance to the leading edge; the FFA-W3-211 airfoil of the IEA 15 MW turbine at 0° angle of attack was chosen; calculated with the methodology from Sect. 3.2.

Since current erosion research treats droplets as spherical and thus as a coherent structure when impacting on a blade, it is also assumed that the entire water mass of a single droplet possesses the same impact velocity. The question is whether droplet breakup invalidates this assumption. If the velocity that describes the droplet deformation is on the same order of magnitude as the droplet slowdown itself, then, with Eq. (1) in mind, the damage potential of a droplet might be significantly influenced. To understand this aspect further, additional frames of the 0.191 mm droplet from Table 1 are shown in Fig. 7. In six distinct frames, a purple and an orange arrow indicate the maximum extent of the bag that forms during the breakup. With the timestamp and indicated length scale, the growth velocity of this bag can be obtained using a simple backward Euler finite-difference scheme. The velocities obtained are indicated next to the frame number in the figure. It can be seen that the velocity is fairly low in Frames 2 and 3 when the bag is just beginning to form. However, as soon as the bag starts rapidly growing, the velocity quickly increases to a peak of 42.66 m s⁻¹. Close to impact, this velocity reduces to a still significant value of 22 m s⁻¹. This example shows that the water's velocity inside a droplet that undergoes breakup (close to impact) is not constant in space and time. The exact velocity distribution inside the droplet is probably breakup-mode-dependent, and droplets that only undergo gradual deformation will preserve a reasonably constant velocity throughout the droplet. To further elaborate on this argument, if droplets fracture into sub-droplets during a breakup, each resulting droplet will have a distinct impact velocity. To conclude the findings, experiments suggest that droplets approaching wind-turbine-sized airfoils either are deformed or will show breakup shortly before impact. Additionally, droplets that undergo breakup will show a non-homogeneous impact velocity distribution across their water mass.

Figure 7 A 0.191 mm diameter droplet approaching an airfoil and showing a bag breakup mode; the scale, time, frame identifier, and bag expansion speed are also indicated; the photographs are reproduced with permission from Fig. 5.6 of García-Magariño (2016).

3 Methodology

3.1 One-way coupled Lagrangian particle model

The influence of droplet deformation and breakup on the blade lifetime under erosion is investigated with a model that adequately describes the relevant physical processes. Various Lagrangian droplet deformation models exist in the literature, such as the TAB, NLTAB3, DDB, and Droplet Ratio Deformation (DRD) models (Sor and García-Magariño, 2015; Schmehl, 2004). However, to the authors' knowledge, to date, no single Lagrangian model can describe the full range of complex phenomena of droplet slowdown and breakup in sufficient detail. Some advanced models attempt to model particular regimes, such as in Sichani and Emami (2015) for a droplet under deformation and up to the onset of bag rupture. Direct numerical simulation (DNS) of the air and water phase could capture all relevant physics and phenomena, especially when higher-order numerical schemes are applied. However, its computational expense makes it prohibitive when a large parameter space is supposed to be studied. Thus a gap exists with computationally affordable but low-accuracy Lagrangian particle models on one side and highly accurate but extremely costly DNS codes on the other.

This dilemma is resolved by simplifying the problem based on educated assumptions. In particular, it is argued that the model's foremost aim must be the accurate prediction of the droplet slowdown velocity. As shown in Eq. (1), a small error in the impact velocity leads to a large error in the computed erosion lifetime. The second central simplification is that, for the conclusions of this study, the exact shape of the droplet on impact does not need to be predicted very accurately. This simplification is based on the assumption that an error in the droplet's shape during impact has a smaller influence on the erosion lifetime than an error in the impact velocity. However, at the same time, the prediction of the droplet's shape prior to impact needs to be accurate enough to minimize the error in the impact velocity. It is noteworthy that the droplet's shape on impact can be an input for a damage metric that is required to calculate an erosion lifetime. This aspect is discussed in Sect. 3.4.

Additional simplifying assumptions are made to model the problem in a Lagrangian one-particle setting. It is assumed that droplets will preserve a coherent shape during the entire approach toward the airfoil, i.e., not fracture, and thus can be represented as a single particle. This assumption neglects the potential effect of the non-homogeneous impact velocity of rain droplets during and after breakup. Based on the reference measurements from the literature that were presented before, it is also assumed that the cascade breakup does not occur.

Considering these requirements, the DRD model from Sor and García-Magariño (2015) and Sor et al. (2016) was chosen. It was specifically developed to compute the trajectory of water droplets in the vicinity of approaching airfoils and stems from the research group that has also published the measurements on droplet breakup discussed before. It has shown superior performance compared to other droplet models and is based on a one-way coupled Lagrangian approach. The original method uses three equations. One equation models the rain droplet's deformation from a sphere to the shape of an oblate spheroid. The other two equations model the movement of the droplet in a two-dimensional space. For the present study, the model was modified in such a way that the movement of the droplet can be considered one-dimensional only. It is important to note that the DRD model neither accounts for droplet breakup nor imposes any limit on the maximum deformation of a droplet. As a remedy, a heuristic modification is proposed in the following.

Two fundamental equations of motion (EOMs) are at the model's core. They read

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle }m{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}x}{\mathrm{d}{t}^{\mathrm{2}}}}=F_{\mathrm{drag}},\text{(6)}& {\displaystyle }& {\displaystyle \frac{\mathrm{3}}{\mathrm{16}}}m{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}a}{\mathrm{d}{t}^{\mathrm{2}}}}=F_{\mathit{\sigma }}+F_{\mathrm{p}}.\end{array}$$

Equation (5) represents the EOM along the droplet trajectory, whereas Eq. (6) is the EOM that represents the deformation of the droplet from a spheroid to an oblate spheroid. $m=\mathrm{4}/\mathrm{3}\mathit{\pi }{R}_{\mathrm{0}}^{\mathrm{3}}{\mathit{\rho }}_{\mathrm{water}}$ is the mass of the droplet, and x is the position of the droplet along its trajectory. The possible candidates for the droplet trajectory will later be assessed in Sect. 3.2 together with Fig. 10. a is the semi-major axis of an oblate spheroid, as shown in Fig. 8. b is the semi-minor axis and can be calculated as $b=R_{\mathrm{0}}^{\mathrm{3}}/a^{\mathrm{2}}$, where R₀ is the starting radius of a spherical rain droplet.

Figure 8 Representation of the droplet shapes used in the model with the relevant geometrical parameters: (a) spheroid and (b) oblate spheroid.

The drag force acting on the droplet is computed using

$$\begin{array}{}\text{(7)}& F_{\mathrm{drag}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\rho }}_{\mathrm{air}}{V}_{\mathrm{slip}}^{\mathrm{2}}{C}_{\mathrm{D}}{A}_{\mathrm{a}}.\end{array}$$

V_slip is the velocity difference between the air and the droplet; it reads

$$\begin{array}{}\text{(8)}& V_{\mathrm{slip}}=V_{\mathrm{air}}-{\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}}.\end{array}$$

The calculation of the background velocity V_air at a particular value of x is explained in Sect. 3.2. The droplet's instantaneous frontal area A_a is calculated by simply taking A_a=πa². The drag coefficient is composed of a static and a dynamic component:

$$\begin{array}{}\text{(9)}& C_{\mathrm{D}}=C_{\mathrm{static}}+C_{\mathrm{dynamic}}.\end{array}$$

They read

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }{C}_{\mathrm{static}}=C_{{\mathrm{D}}_{\mathrm{sphere}}}^{b/a}{C}_{{\mathrm{D}}_{\mathrm{disk}}}^{\mathrm{1}-b/a},\text{(11)}& {\displaystyle }& {\displaystyle }{C}_{\mathrm{dynamic}}=k{\displaystyle \frac{b}{V_{\mathrm{slip}}^{\mathrm{2}}}}{\displaystyle \frac{\mathrm{d}{V}_{\mathrm{slip}}}{\mathrm{d}t}},\end{array}$$

where k is a calibration constant. The static component represents an interpolation between the drag coefficient of a sphere and a disk. In Eq. (6), two forces are acting against each other. The surface pressure term drives deformation, whereas the surface tension term counteracts deformation. The pressure term is calculated as

$$\begin{array}{}\text{(12)}& F_{\mathrm{p}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\rho }}_{\mathrm{air}}{V}_{\mathrm{slip}}^{\mathrm{2}}{C}_{\mathrm{p}}{A}_{\mathrm{0}}.\end{array}$$

C_p is again a calibration constant. Additionally, note the constant frontal area that is calculated with the initial droplet radius R₀; that is $A_{\mathrm{0}}=\mathit{\pi }{R}_{\mathrm{0}}^{\mathrm{2}}$. This choice is motivated in more detail in the original paper of the model. The surface tension force is written as

$$\begin{array}{}\text{(13)}& F_{\mathit{\sigma }}=-{\displaystyle \frac{\mathrm{4}}{\mathrm{3}}}{\mathit{\sigma }}_{\mathrm{water}}{\displaystyle \frac{\mathrm{d}{A}_{\mathrm{s}}}{\mathrm{d}a}},\end{array}$$

where σ_water is the surface tension of water and $\frac{\mathrm{d}{A}_{\mathrm{s}}}{\mathrm{d}a}$ is the derivative of the surface area of an oblate spheroid with respect to a. Following the approach of Sor and García-Magariño (2015), Sor et al. (2016), Sor and García-Magariño (2021), and Schmehl (2004), the surface area of an oblate spheroid A_s reads

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{A_{\mathrm{s}}}{A_{\mathrm{s},\mathrm{0}}}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{a}{R_{\mathrm{0}}}}\right)}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{R_{\mathrm{0}}}{a}}\right)}^{\mathrm{4}}{\displaystyle \frac{\mathrm{arctanh}\mathit{ϵ}}{\mathit{ϵ}}},\end{array}$$

where $A_{\mathrm{s},\mathrm{0}}=\mathrm{4}\mathit{\pi }{R}_{\mathrm{0}}^{\mathrm{2}}$ is the surface area of a sphere. The derivative becomes

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{A_{\mathrm{s},\mathrm{0}}}}{\displaystyle \frac{\mathrm{d}{A}_{\mathrm{s}}}{\mathrm{d}\stackrel{\mathrm{‾}}{a}}}=& \stackrel{\mathrm{‾}}{a}-{\displaystyle \frac{\mathrm{2}}{{\stackrel{\mathrm{‾}}{a}}^{\mathrm{5}}}}{\displaystyle \frac{\mathrm{arctanh}\mathit{ϵ}}{\mathit{ϵ}}}+{\displaystyle \frac{\mathrm{3}}{\mathrm{2}{\stackrel{\mathrm{‾}}{a}}^{\mathrm{5}}\left({\stackrel{\mathrm{‾}}{a}}^{\mathrm{6}}-\mathrm{1}\right)}}\\ & \left({\stackrel{\mathrm{‾}}{a}}^{\mathrm{6}}-{\displaystyle \frac{\mathrm{arctanh}\mathit{ϵ}}{\mathit{ϵ}}}\right),\end{array}\end{array}$$

where

$$\begin{array}{}\text{(16)}& \mathit{ϵ}=\sqrt{\mathrm{1}-{\left({\displaystyle \frac{b}{a}}\right)}^{\mathrm{2}}}=\sqrt{\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{{\stackrel{\mathrm{‾}}{a}}^{\mathrm{6}}}}}\end{array}$$

and

$$\begin{array}{}\text{(17)}& \stackrel{\mathrm{‾}}{a}={\displaystyle \frac{a}{R_{\mathrm{0}}}}.\end{array}$$

Finally $\frac{\mathrm{d}{A}_{\mathrm{s}}}{\mathrm{d}a}$ is obtained by

$$\begin{array}{}\text{(18)}& {\displaystyle \frac{\mathrm{d}{A}_{\mathrm{s}}}{\mathrm{d}a}}={\displaystyle \frac{\mathrm{1}}{R_{\mathrm{0}}}}{\displaystyle \frac{\mathrm{d}{A}_{\mathrm{s}}}{\mathrm{d}\stackrel{\mathrm{‾}}{a}}}.\end{array}$$

$C_{{\mathrm{D}}_{\mathrm{sphere}}}$ has been calculated with the Schiller–Naumann relation as given in Sommerfeld et al. (2008):

$$\begin{array}{}\text{(19)}& C_{{\mathrm{D}}_{\mathrm{sphere}}}=\left\{\begin{array}{ll}\mathrm{27.6}& Re\le \mathrm{1},\\ \frac{\mathrm{24}}{Re}\left(\mathrm{1}+\mathrm{0.15}R{e}^{\mathrm{0.687}}\right)& \mathrm{1}<Re<\mathrm{1000},\\ \mathrm{0.4383}& Re\ge \mathrm{1000}.\end{array}\right.\end{array}$$

Note that the drag coefficient was clamped for Re≤1 and Re≥1000. The Reynolds number Re reads

$$\begin{array}{}\text{(20)}& Re={\displaystyle \frac{V_{\mathrm{slip}}{\mathit{\rho }}_{\mathrm{air}}\mathrm{2}{R}_{\mathrm{0}}}{{\mathit{\mu }}_{\mathrm{air}}}}.\end{array}$$

In the original form, the model does not account for the influence of droplet breakup. The model permits the droplet to grow without restriction. From the literature, such as Jackiw and Ashgriz (2021, 2022), Hsiang and Faeth (1995), and Schmehl (2004), it is known that, depending on the Weber number, there exists a maximum diameter at which droplets start to break up. It usually lies in the range of 1.5 to $\mathrm{2}a/R_{\mathrm{0}}$. In this study, the assumption is made that when the droplets reach a specific maximum a, they will stop growing, and the value of a will be fixed for the remainder of the simulation. In particular, the following formula is used:

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{a_{\max }}{R_{\mathrm{0}}}}=\min \left(\mathrm{2.2},\mathrm{3.4966}{\mathit{We}}_{\mathrm{\infty }}^{-\mathrm{0.1391}}\right).\end{array}$$

The formula was obtained by fitting a power function to a set of reference data shown in Fig. 9. Further, the limit of 2.2 was chosen based on the data in Fig. 3 of García-Magariño et al. (2021). It is important to note that

$$\begin{array}{}\text{(22)}& {\mathit{We}}_{\mathrm{\infty }}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{air}}{V}_{\mathrm{\infty }}^{\mathrm{2}}{\mathit{\varphi }}_{\mathrm{0}}}{{\mathit{\sigma }}_{\mathrm{water}}}},\end{array}$$

which is different to Eq. (3) since the free-stream velocity (Eq. 27) is used instead of the slip velocity (Eq. 8). The motivation for this is the fact that for the limited sets of published data on $a_{\max }/R_{\mathrm{0}}$, the corresponding impact velocity is not always given. Additionally, We_impact is not known a priori but rather as a result of the simulation, therefore necessitating an iterative approach for solving the set of equations. The assumption can be justified by realizing that V_impact≈V_∞ represents a conservative estimate. Since $a_{\max }/R_{\mathrm{0}}$ should be decreasing with increasing We_impact, assuming We_impact≈We_∞ will lead to a higher estimated Weber number. Thus, droplet slowdown will be underpredicted due to an underprediction in $a_{\max }/R_{\mathrm{0}}$. Section 3.3 shows that the limiter introduced here deals with the droplet breakup satisfactorily.

Figure 9 The limit of maximum droplet dimension $a_{\max }/R_{\mathrm{0}}$ is a function of the free-stream Weber number; the sources of data points are given in Table 2.

Vargas et al. (2012)García-Magariño (2016)Vargas et al. (2012)Feo et al. (2012)García-Magariño (2016)García-Magariño (2016)Sor (2017)Sor (2017)Sor (2017)Feo et al. (2012)

Table 2 Sources of reference data for the maximum droplet dimension limiter $a_{\max }/R_{\mathrm{0}}$ that is shown in Fig. 9.

The resulting set of differential equations describing the droplet model is

$$\begin{array}{}\text{(23)}& {\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}}=V_x,{\displaystyle \frac{\mathrm{d}{V}_x}{\mathrm{d}t}}={\displaystyle \frac{F_{\mathrm{drag}}}{m}},\end{array}$$

$$\begin{array}{}\text{(24)}& {\displaystyle \frac{\mathrm{d}a}{\mathrm{d}t}}=V_{\mathrm{a}},{\displaystyle \frac{\mathrm{d}{V}_{\mathrm{a}}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{16}}{\mathrm{3}}}{\displaystyle \frac{F_{\mathit{\sigma }}+F_{\mathrm{p}}}{m}}.\end{array}$$

The initial conditions for the droplet equations are set as

$$\begin{array}{}\text{(25)}& {\displaystyle }& {\displaystyle }{x}_{\mathrm{0}}=\mathrm{0},V_{x,\mathrm{0}}=\mathrm{0},\text{(26)}& {\displaystyle }& {\displaystyle }{a}_{\mathrm{0}}=R_{\mathrm{0}}+\mathrm{eps},V_{\mathrm{a},\mathrm{0}}=\mathrm{0}.\end{array}$$

“eps” is a very small number, e.g., $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{12}}$. This is necessary since Eq. (15) is not defined for a=R₀. V_x is the velocity of the droplet, and V_a is the expansion velocity of the semi-major axis. An additional differential equation is needed to describe the movement of the blade. It reads

$$\begin{array}{}\text{(27)}& {\displaystyle \frac{\mathrm{d}{x}_{\mathrm{blade}}}{\mathrm{d}t}}=V_{\mathrm{blade}}=-V_{\mathrm{\infty }}=\mathrm{const},\end{array}$$

with the initial conditions sufficiently far away from the droplet:

$$\begin{array}{}\text{(28)}& {\displaystyle }{x}_{\mathrm{blade},\mathrm{0}}\gg R_c.\end{array}$$

Here a sufficiently far distance can, for example, be 20R_c. The definition of R_c is explained in Sect. 3.2. In this study, the differential equations were solved using a simple Runge–Kutta method. The simulation is stopped when the distance between the airfoil and the droplet falls below a certain threshold; i.e.,

$$\begin{array}{}\text{(29)}& \mathrm{\Delta }x=|x-x_{\mathrm{blade}}|<\mathrm{eps},\end{array}$$

where “eps” is once again a small number. The relative velocity between the droplet and the airfoil can be defined as

$$\begin{array}{}\text{(30)}& \mathrm{\Delta }V=|V_x-V_{\mathrm{blade}}|;\end{array}$$

see Figs. 2 and 18. In addition, the slowdown velocity can be defined as

$$\begin{array}{}\text{(31)}& V_{\mathrm{slowdown}}={\left({\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}}\right)}_{\mathrm{at}\mathrm{impact}},\end{array}$$

which is just the velocity gained by the droplet, since V_blade=const. The slowdown velocity can be interpreted as the reduction in impact velocity.

Table 3 summarizes the physical and calibrations constants used in the model. In the original method of Sor et al. (2016) C_p was given as C_p=0.93. However, in this study, it was found that setting $C_{\mathrm{p}}=C_{{\mathrm{D}}_{\mathrm{disk}}}$ provided results that matched more closely the impact velocities of the validation cases in Fig. 18.

Sor et al. (2016)Sor et al. (2016)

Table 3 The physical properties are given for an ambient temperature 288.15 K and ambient pressure of 101 325 Pa.

3.2 Calculation of the background velocity

A necessary input to the model is the background velocity field V_air. The droplet traverses this field while approaching the airfoil (see Eq. 8). It is dependent on the size and shape of the wind turbine's airfoil. This study treats the problem as one-dimensional. From this assumption, a range of possibilities for the implied trajectory of the droplet emerge. Figure 10 illustrates these possibilities. In the limit, there are two possible trajectories for small and large droplets, respectively. Very small droplets are expected to follow the streamline of the flow, while large droplets are expected to follow a ballistic trajectory. In practice, the rain droplets will follow a trajectory that lies in the region between these two. To find a characteristic velocity field that can be used for further study, two popular reference turbine designs were chosen, the NREL 5 MW and the IEA 15 MW (Jonkman et al., 2009; Gaertner et al., 2020).

Figure 10 Ballistic and streamline trajectories of a droplet approaching an airfoil operated under an angle of attack α; the blade twist angle and pitch are set to zero.

The first step taken here toward obtaining V_air is to perform CFD calculations of the flow field surrounding the reference turbines' airfoils using OpenFOAM. The simulations were carried out using the simpleFoam solver with the k–ω shear stress transport (SST) turbulence model. A free-stream velocity of 90 m s⁻¹ was chosen. A 2D mesh of around 100 000 cells has been used for the computations. In this application, a rather coarse computational grid is satisfactory since the area of interest is located ahead of the leading edge. In this area, the solutions are well-behaved, and problematic areas with flow separation are located far downstream. Subsequently, the one-dimensional velocity field was extracted from the solution using ParaView. Two fields were extracted, one for the ballistic trajectory and one for the streamline trajectory. The latter was obtained by seeding an upstream streamline from the leading edge in ParaView and subsequently extracting the velocity vector along this line.

Instead of directly using the extracted fields as a model input, they were parameterized, which allows us to better compare the different airfoils by looking at the model parameters. As in Lopez-Gavilan et al. (2020), the underlying parametrization model is the potential flow solution of a cylinder representing the nose of the airfoil. The horizontal velocity component for the potential flow in the stagnation streamline, when viewed from the airfoil, reads

$$\begin{array}{}\text{(32)}& {\displaystyle \frac{V_{\mathrm{air}}}{V_{\mathrm{\infty }}}}=\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{{\left(\mathrm{1}+\frac{\mathrm{\Delta }x}{R_c\left(\mathit{\alpha }\right)}\right)}^n}}.\end{array}$$

For Sect. 3.1 it is advised to transform this equation to the reference frame of an external observer.

R_c(α) is the radius of the cylinder, and Δx is the distance (from the droplet) to the cylinder. However, here n and R_c(α) are free parameters that are fitted to the extracted field from the CFD simulation. Therefore, R_c(α) and n should not be regarded as geometric but rather as aerodynamic parameters; i.e., R_c(α) is the aerodynamic nose radius. It is also a function of the angle of attack. It is heuristically found that it is possible to collapse the one-dimensional velocity field for different angles of attack; i.e., the solution is self-similar to an (arbitrary) scaling value. In this case, the self-similar variable is taken to be the distance from the leading edge at which the velocity has dropped to the 95 % value of the free-stream velocity (x_95 %). This self-similar property is shown in Fig. 11. The left plot shows the velocity field against the dimensional distance to the leading edge. In the right plot the velocity field is collapsed by scaling with x_95 %.

Figure 11 Non-dimensional velocity field along the stagnation streamline vs. the dimensional (a) and non-dimensional (b) distance to the leading edge; airfoil: FFA-W3-211.

The self-similarity allows us to represent the velocity field at different angles of attack by scaling R_c(α). The influence of the angle-of-attack variation on the self-similar parameters can be sufficiently represented by considering a second-order polynomial. Therefore,

$$\begin{array}{}\text{(33)}& x_{\mathrm{95}\mathit{\%}}\left(\mathit{\alpha }\right)=\left(C_{\mathrm{1}}{\mathit{\alpha }}^{\mathrm{2}}+C_{\mathrm{2}}\mathit{\alpha }+\mathrm{1}\right)x_{\mathrm{95}\mathit{\%},\mathrm{0}},\end{array}$$

and thus also

$$\begin{array}{}\text{(34)}& R_c\left(\mathit{\alpha }\right)=\left(C_{\mathrm{1}}{\mathit{\alpha }}^{\mathrm{2}}+C_{\mathrm{2}}\mathit{\alpha }+\mathrm{1}\right)R_{c,\mathrm{0}},\end{array}$$

where R_c,0 is the aerodynamic nose radius at a 0° angle of attack. In Fig. 12 the FFA-W3-211 airfoil's variation of x_95 %(α) is shown in conjunction with the polynomial fit. The first step in the parametrization process is to find x_95 %(α) for every angle of attack. The 0° angle-of-attack field is then used to find the parameters n and R_c,0. Last but not least, the parameters C₁ and C₂ are found by fitting the polynomial to x_95 %(α). All best-fit parameters were found using MATLAB's fmincon function and a least-squares minimization function.

Figure 12 x_95 % location of the FFA-W3-211 airfoil as a function of the angle of attack; the velocity field is along the streamline trajectory.

The final parameters are given in Table 4. The values of R_c(α) were made dimensionless with the airfoil chord c. The table shows a general trend when comparing thicker airfoils to thinner airfoils. Thicker airfoils have a higher aerodynamic nose radius and exponent. Therefore, x_95 % is also higher, meaning thicker airfoils influence droplets farther upstream. Two diverging behaviors can be noticed regarding the parameters for the angle-of-attack correction. For the flow that follows a ballistic trajectory, an increasing angle of attack leads to a decreasing R_c(α), whereas for the flow along the stagnation streamline, an increasing R_c(α) can be noticed. Therefore, in comparison to the 0° angle of attack, small droplets are expected to be influenced more, whereas large droplets are expected to be influenced less when the angle of attack is increased.

Table 4 Best-fit parameters of V_air for the NREL 5 MW and IEA 15 MW turbine airfoils; subscripts “b” and “s” stand for ballistic and streamline path, respectively; the coefficients C₁ and C₂ are dimensional; they are given for the angle of attack α in degrees – see Eq. (34); their units are given in the parentheses of the column header.

Figures 13 and 14 show the dimensional aerodynamic nose radius R_c and the exponent n along the blade of the NREL 5 MW and the IEA 15 MW wind turbines. For consistency with the official definitions of both reference turbines, the airfoils were kept constant between the officially defined stations of the NREL turbine but were linearly interpolated between the stations of the IEA turbine. This explains the sawtooth pattern in the results of the NREL turbine in both plots. Figure 15 gives the angle-of-attack distributions along the blade that are used in Eq. (34). The influence on the angle of attack on the aerodynamic nose radius is less than 5 % for the ballistic trajectory and around 30 % for the streamline trajectory. Both R_c and n are larger for the IEA reference turbine than for the NREL design. This has three main reasons. The IEA turbine has a higher dimensionless aerodynamic nose radius $R_c/c$ for its airfoils. It also has a larger chord, and the angles of attack are higher. Due to the similarity of the IEA 15 MW turbine to current state-of-the-art offshore turbines, it is argued that the values of R_c=0.07 m and n=1.1, as they can be found at around $r/R_{\mathrm{blade}}=\mathrm{0.9}$, represent a good baseline for the remainder of this study. With these findings in mind, it is worth noting that with R_c=0.071 m and n=1.2, one obtains a very good fit of the reference velocity field of Case F and Case G (see Table 5). Hence, the parameter space of the reference data is close to the parameter space encountered in leading-edge erosion.

Figure 13 Dimensional aerodynamic nose radius R_c along the dimensionless blade distance for the IEA 15 MW and NREL 5 MW turbines.

3.3 Validation of the model

Two tests are performed to validate the model. First, the model is compared against well-known relations for the terminal falling conditions of water droplets. To this end, V_air was set to zero and a gravity term was added to Eq. (5). Secondly, for comparison, a set of rotating-arm test-rig reference data are compiled from different sources.

Figure 14 Aerodynamic exponent n along the dimensionless blade distance for the IEA 15 MW and NREL 5 MW turbines.

Figure 15 Element-wise angle of attack α along the dimensionless blade distance. The data were obtained from vortex method simulations in Barfknecht et al. (2022) with a tip-speed ratio of 9 for the IEA 15 MW turbine and of 7.55 for the NREL 5 MW turbine.

Best (1950b) gives a relation for the terminal velocity of falling water droplets that reads

$$\begin{array}{}\text{(35)}& V_{\mathit{\varphi }}=\mathrm{9.32}{e}^{\mathrm{0.0405}h}\left(\mathrm{1}-e^{-{\left(\mathrm{0.565}{\mathit{\varphi }}_{\mathrm{0}}\right)}^{\mathrm{1.147}}}\right),\end{array}$$

where h is the altitude in kilometers, which was set to zero for this study. In this equation, the droplet diameter ϕ₀ must be given in millimeters. A relation for the shape ($a_{\mathit{\varphi }}/R_{\mathrm{0}}$) of droplets in terminal conditions is given by Brandes et al. (2002). It reads

$$\begin{array}{}\text{(36)}& \begin{array}{rl}{\displaystyle \frac{a_{\mathit{\varphi }}}{R_{\mathrm{0}}}}& =\left(\mathrm{0.9951}+\mathrm{0.02510}{\mathit{\varphi }}_{\mathrm{0}}-\mathrm{0.03644}{\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{2}}\right.\\ & {\left.+\mathrm{0.005030}{\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{3}}-\mathrm{0.0002492}{\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{4}}\right)}^{-\mathrm{1}/\mathrm{3}}.\end{array}\end{array}$$

Note, in the original formulation, Eq. (36) was given as the ratio $a/b$. This has been converted here to $a_{\mathit{\varphi }}/R_{\mathrm{0}}$ by assuming the shape of an oblate spheroid. In addition, ϕ₀ must be given in millimeters for this equation. Figures 16 and 17 compare both formulas with the results from the model. Excellent agreement is achieved for the shape of the droplet. For droplet diameters of up to 2.5 mm, the terminal velocity is almost identical to the reference value. Afterwards, a slight deviation can be noticed.

Figure 16 Validation of terminal velocity for falling water droplets as a function of the droplet diameter.

Figure 17 Validation of the terminal dimensionless semi-major axis for falling water droplets as a function of the droplet diameter.

A set of reference data from the literature have been compiled for the second test. Unfortunately, the data quality differs based on whether they were directly available or had to be derived by, for example, measuring distances from published images of high-speed photography. Table 5 summarizes the reference cases.

Vargas and Feo (2011)Vargas and Feo (2011)Vargas and Feo (2011)Vargas and Feo (2011)Vargas and Feo (2011)Vargas et al. (2012)Vargas et al. (2012)García-Magariño (2016)

Table 5 Summary of rotating-arm test-rig reference data used in the validation of the proposed model.

Figure 18 compares the model and measurements. It can be seen that there is good agreement between them. The model overpredicts the slowdown for Case F. Interestingly, the slight discrepancy starts already at a distance of about 0.05 m from the leading edge, a region where the other cases show excellent agreement. Cases D and E suffer from a slight underprediction of the slowdown close to the leading edge. Arguably, Case C overpredicts the slowdown. Data extraction of Case H was challenging and had to be done manually from a small series of published photographs. Therefore, the data can only be considered fair. Nevertheless, the simulation and measurements still agree reasonably well.

Figure 18 Validation of the trajectory model with reference data summarized in Table 5; markers indicate reference data and solid lines the results of the model; note that the y axis contains repeated ticks for better visualization of cases with equal free-stream velocity.

To summarize, the model agrees well with reference data for both validation cases. Recall that even slight differences in the impact speed will lead to very different lifetime predictions due to the large exponent in Eq. (1). Nevertheless, with the available data and the simple reduced-order Lagrangian model in mind, the validation results are considered adequate for lifetime predictions.

3.4 Damage model and relevant rain droplet diameters

A damage model is required to evaluate the magnitude of the slowdown effect on the lifetime of a blade. The damage model proposed in this paper is described in the following. Additionally, the equations developed here are used to compute the relevant rain droplet diameter range for the present study.

Several damage metrics have been proposed to calculate an erosion lifetime: the water hammer pressure metric, which is often used in conjunction with the Springer model (Hoksbergen et al., 2022); impingement (Bech et al., 2022); kinetic energy (Bech et al., 2018); or the material's strain (Verma et al., 2020). Arguably, the two most common models are currently the Springer model and the impingement metric. This study uses the impingement metric to calculate an erosion lifetime. The choice is explained in the following.

The Springer model (as described in Hoksbergen et al., 2022) gives an equation for the erosion lifetime by considering the number of allowable repeated impacts on one location $N_i^{*}$. The model is derived by computing the impact force F=p_whA_projected, where p_wh is the (modified) water hammer pressure and A_projected is the projected area of the droplet onto the impact target. In the case of an oblate spheroid, this would be $A_{\mathrm{projected}}=A_{\mathrm{a}}=\mathit{\pi }{a}^{\mathrm{2}}$. The assumption is made that the water hammer pressure is constant for the entire projected area of the droplet. Subsequently, a stress field within the target is computed using an analytical equation of the form $\mathit{\sigma }(F,r,\mathrm{\dots })$, where r is the distance to the impact location. Further, σ∝F. After some steps, $N_i^{*}$ is obtained. The entire derivation for the (uncoated) Springer model is given in Springer and Baxi (1972). A problematic assumption within the Springer model is the calculation of the impact force. If, for example, a single droplet is infinitely stretched, that is A_projected→∞, then σ→∞ and therefore $N_i^{*}\to \mathrm{0}$. Alternatively, a droplet that is significantly squeezed, i.e., A_projected→0, will have a lifetime of $N_i^{*}\to \mathrm{\infty }$. Both results seem unphysical and thus question the validity of the Springer model. Since the rain droplets deform significantly and, therefore, grow in the projected area, the Springer model does not seem to be an adequate choice for the present study.

Impingement is a damage metric representing the total water column that the blade intercepts until coating failure. Since impingement only considers the amount of water, it is, at least conceptually, agnostic to the impacting droplet's shape, a property that seems advantageous considering the complex shape of droplets during impact. Due to this property and its recent gain in popularity, as shown in Bech et al. (2022), Visbech et al. (2023), and Badger et al. (2022), it was chosen as the damage metric for this study.

The general formula for the accumulated impingement H during operation is

$$\begin{array}{}\text{(37)}& H\left(t\right)=W{V}_{\mathrm{collection}}t.\end{array}$$

W is the accumulated water column in meters per meter of swept air, t is the time, and V_collection is the speed at which water is collected. Here it is assumed that V_collection=V_∞.

The impingement until end of incubation, dubbed allowed impingement, is also required. H_allowed reads

$$\begin{array}{}\text{(38)}& H_{\mathrm{allowed}}={\displaystyle \frac{\mathit{\alpha }}{V_{\mathrm{impact}}^{\mathit{\beta }}}}.\end{array}$$

The equation has the form of Eq. (1). The parameters α and β were found using the measurements of Bech et al. (2022). They performed measurements in a rotating-arm erosion test rig, where they recorded H_allowed of a generic blade coating with respect to V_impact. Tests with four distinct droplet sizes ranging from 0.76 to 3.5 mm were performed. The measurements are shown in Fig. 19. Their raw data were used in this study to fit a function through the data points, leading to the best-fit parameters of $\mathit{\alpha }=\mathrm{3.4860}\times {\mathrm{10}}^{\mathrm{20}}$ and β=9.5774. Figure 19 shows that the measurements collapse well. It should be noted that the authors of the study argue that the data show some drop-size dependency with $\mathrm{7.2}\le \mathit{\beta }\le \mathrm{10.5}$. This range is found when best-fit functions are created for every droplet size individually. Nevertheless, the assumption made here is that this dependency can be neglected for the conclusions drawn in this study.

Since, as indicated earlier, there is a wide spread of reported values for β in the literature, two other exponents were considered to ensure the robustness of the drawn conclusions with respect to β. The other two exponents that were chosen are 5.7 and 7. The exponent of 5.7 originates from the Springer model (Hoksbergen et al., 2022). Even though Springer does not measure impingement, but rather impacts (per surface area), it is still considered to be worth showing. The exponent 7 represents an arbitrary value between 5.7 and 9.5774.

Figure 19 Rotating-arm erosion test-rig results by Bech et al. (2022) that relate impact velocity to impingement.

Equations (37) and (38) can be used in a Palmgren–Miner damage rule, yielding the total damage that reads

$$\begin{array}{}\text{(39)}& D=T_{\mathrm{rain}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{{\partial }_t{H}_{\mathrm{I},\mathit{\varphi }}}{H_{\mathrm{allowed}}}}\mathrm{d}\mathit{\varphi }\mathrm{d}I.\end{array}$$

T_rain is the time of operation during rain. ${\partial }_t{H}_{\mathrm{I},\mathit{\varphi }}$ is the impingement collection rate as a function of the rain intensity I and droplet diameter ϕ₀. It is defined as

$$\begin{array}{}\text{(40)}& {\partial }_t{H}_{\mathrm{I},\mathit{\varphi }}=W_{\mathit{\varphi }}{V}_{\mathrm{collection}}{f}_{\mathrm{I}}\end{array}$$

and is analogous to Eq. (37) but with W_ϕ, which depends on the droplet diameter ϕ₀, given by

$$\begin{array}{}\text{(41)}& W_{\mathit{\varphi }}={\displaystyle \frac{f_{\mathit{\varphi },\mathrm{plane}}I}{V_{\mathit{\varphi }}}},\end{array}$$

where f_ϕ,plane is a distribution that describes the amount of water associated with every droplet diameter that passes through an imaginary plane in the air. One popular model that can be used to obtain f_ϕ,plane is the Best model (Best, 1950a). It gives a probability density function (pdf) of the water mass associated with every droplet diameter in a control volume in air and is given as

$$\begin{array}{}\text{(42)}& f_{\mathit{\varphi },\mathrm{air}}=\mathrm{2.25}{\left({\displaystyle \frac{\mathrm{1}}{\mathrm{1.3}{I}^{\mathrm{0.232}}}}\right)}^{\mathrm{2.25}}{\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{2.25}-\mathrm{1}}{e}^{-{\left(\frac{{\mathit{\varphi }}_{\mathrm{0}}}{\mathrm{1.3}{I}^{\mathrm{0.232}}}\right)}^{\mathrm{2.25}}},\end{array}$$

where here, notice the units, I is the rain intensity in millimeters per hour and ϕ₀ is the droplet diameter in millimeters! f_ϕ,air can be converted into f_ϕ,plane using

$$\begin{array}{}\text{(43)}& f_{\mathit{\varphi },\mathrm{plane}}={\displaystyle \frac{f_{\mathit{\varphi },\mathrm{air}}{V}_{\mathit{\varphi }}}{\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{f}_{\mathit{\varphi },\mathrm{air}}{V}_{\mathit{\varphi }}\mathrm{d}\mathit{\varphi }}}.\end{array}$$

Note that if f_ϕ,plane is supposed to be obtained for droplet diameters in meters, then the integral in the denominator should be computed with ϕ in meters. f_ϕ,plane is plotted for five different rain intensities in Fig. 20. One can see that the water volume of lighter rain events is mainly composed of droplets with smaller diameters on the order of 0.5 to 1 mm. With increasing rain intensity, the amount of water contained in larger droplets is increasing.

Figure 20 Best's (Eq. 43) distribution over a plane as a function of droplet diameter in millimeters for five different rain intensities.

The time of operation during rain, T_rain, i.e., the rain collection time, over 1 year of operation is given by

$$\begin{array}{}\text{(44)}& T_{\mathrm{rain}}=p_{\mathrm{rain}}{T}_{\mathrm{year},\mathrm{spinning}},\end{array}$$

where T_{year,spinning} is the number of seconds in a year that the turbine spins and p_rain is the probability of rain at a particular site. It should be noted that in the Results section of this study, the damage is presented in its non-dimensional form. p_rain and T_{year,spinning} cancel out during non-dimensionalization since they are both constant. f_I is the probability density function for the various rain intensities. To find f_I, in this study, we consider the coastal site De Kooy located in the Netherlands at coordinates 52.924°, 4.780°. Hourly precipitation data from a 10-year window from 2011 to 2020 were used (KNMI, 2020). The probability density function was determined using the same approach as in Shankar Verma et al. (2021), where a lognormal distribution that reads

$$\begin{array}{}\text{(45)}& f_{\mathrm{I}}={\displaystyle \frac{\mathrm{1}}{I\mathit{\sigma }\sqrt{\mathrm{2}\mathit{\pi }}}}{e}^{-\frac{(\ln I-\mathit{\mu }{)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}\end{array}$$

was fitted using MATLAB's lognfit function applied to the measured precipitation data of the site. μ is the mean, and σ is the standard deviation. Note here the different meanings of the symbols in comparison to before. For a rain intensity given in mm h⁻¹, the coefficients read $\mathit{\mu }=-\mathrm{0.1987}$ and σ=0.9693, whereas when I is given in m s⁻¹, the coefficients become $\mathit{\mu }=-\mathrm{15.29}$ and σ=0.9693.

By combining the previous equations, one obtains the universal Palmgren–Miner damage for an element along the blade reading

$$\begin{array}{}\text{(46)}& \begin{array}{rl}D=& p_{\mathrm{rain}}{T}_{\mathrm{year},\mathrm{spinning}}{V}_{\mathrm{collection}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}I{f}_{\mathrm{I}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\\ & {\displaystyle \frac{f_{\mathit{\varphi },\mathrm{plane}}/V_{\mathit{\varphi }}}{H_{\mathrm{allowed}}\left(V_{\mathrm{impact}}\left(\mathit{\varphi }\right)\right)}}\mathrm{d}\mathit{\varphi }\mathrm{d}I.\end{array}\end{array}$$

Here it is assumed that the turbine always spins at a constant velocity. The formula written in its cumulative form with respect to the rain intensity reads

$$\begin{array}{}\text{(47)}& \begin{array}{rl}{D}_{\mathrm{cumulative}}\left(I\right)=& p_{\mathrm{rain}}{T}_{\mathrm{year},\mathrm{spinning}}{V}_{\mathrm{collection}}\underset{\mathrm{0}}{\overset{\mathrm{I}}{\int }}{I}^{\prime }{f}_{{\mathrm{I}}^{\prime }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\\ & {\displaystyle \frac{f_{\mathit{\varphi },\mathrm{plane}}/V_{\mathit{\varphi }}}{H_{\mathrm{allowed}}\left(V_{\mathrm{impact}}\left(\mathit{\varphi }\right)\right)}}\mathrm{d}\mathit{\varphi }\mathrm{d}{I}^{\prime }.\end{array}\end{array}$$

A special version can be derived that gives the damage associated per meter of impingement at a particular rain intensity. It reads

$$\begin{array}{}\text{(48)}& {\displaystyle \frac{D\left(I\right)}{H\left(I\right)}}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{f_{\mathit{\varphi },\mathrm{air}}}{H_{\mathrm{allowed}}}}\mathrm{d}\mathit{\varphi },\end{array}$$

where H(I) is the collected impingement as a function of rain intensity,

$$\begin{array}{}\text{(49)}& H\left(I\right)=T_{\mathrm{rain}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\partial }_t{H}_{\mathrm{I},\mathit{\varphi }}\mathrm{d}\mathit{\varphi },\end{array}$$

and D(I) the accumulated damage as a function of rain intensity,

$$\begin{array}{}\text{(50)}& D\left(I\right)=T_{\mathrm{rain}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{{\partial }_t{H}_{\mathrm{I},\mathit{\varphi }}}{H_{\mathrm{allowed}}}}\mathrm{d}\mathit{\varphi }.\end{array}$$

The derivation of Eq. (48) uses the fact that

$$\begin{array}{}\text{(51)}& f_{\mathit{\varphi },\mathrm{air}}={\displaystyle \frac{f_{\mathit{\varphi },\mathrm{plane}}/V_{\mathit{\varphi }}}{\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{f}_{\mathit{\varphi },\mathrm{plane}}/V_{\mathit{\varphi }}\mathrm{d}\mathit{\varphi }}}.\end{array}$$

The distribution of water mass that is associated with every droplet diameter at a particular site can be found by combining the functions of f_ϕ,plane and f_I. The result reads

$$\begin{array}{}\text{(52)}& f_{\mathit{\varphi },\mathrm{site}}={\displaystyle \frac{\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}I{f}_{\mathit{\varphi },\mathrm{plane}}{f}_{\mathrm{I}}\mathrm{d}I}{\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}I{f}_{\mathit{\varphi },\mathrm{plane}}{f}_{\mathrm{I}}\mathrm{d}I\mathrm{d}\mathit{\varphi }}},\end{array}$$

with the corresponding cumulative density function of

$$\begin{array}{}\text{(53)}& F_{\mathit{\varphi },\mathrm{site}}=\underset{\mathrm{0}}{\overset{\mathit{\varphi }}{\int }}{f}_{{\mathit{\varphi }}^{\prime },\mathrm{site}}\mathrm{d}{\mathit{\varphi }}^{\prime }.\end{array}$$

Both functions are plotted in Fig. 21. It can be seen that the droplets in the range of 0 to 4 mm contain around 99 % of the total water content. This range needs to be studied for the slowdown effect. The droplets in the range of 0.5 to 3.0 mm account for about 92.5 % of water. $F_{\mathit{\varphi },\mathrm{site},\mathrm{50}\mathit{\%}}$ is found at a diameter of 1.54 mm.

Figure 21 Distribution (dist.) and cumulative function (cum.) of the total rain column associated with every droplet diameter.

4 Results

In this part, the slowdown and deformation model from the previous section is applied. First, in Sect. 4.1 and 4.2, the sensitivities of the droplet diameter and the aerodynamic nose radius to the slowdown and deformation are investigated. Subsequently, in Sect. 4.3, the model's influence on the erosion damage associated with rain intensities is determined. In Sect. 4.4, the distribution of the slowdown velocity along the blades of the two reference turbines is discussed. Finally, these velocities are used to determine an updated damage distribution.

4.1 The influence of droplet diameter and shape on the impact velocity

This section discusses the sensitivity of the droplet slowdown with respect to the diameter. Two droplet types are considered, spherical and oblate droplets. The results of the spherical droplets serve as a conservative bound to the problem and represent the minimum slowdown. As per Fig. 21, droplets in the range of 0.1 to 4 mm were considered for free-stream velocities ranging from 50 to 90 m s⁻¹.

Figure 22a shows the dimensional impact velocity of spherical droplets. A significant slowdown of the droplets can be observed for droplets under 0.5 mm in diameter. Larger droplets show a more gradual slowdown. The origin of this behavior can be found in the ratio of surface area to mass, which is much larger for smaller droplets, thus making them more affected by the drag force. Non-dimensionalizing the impact velocity reveals that the impact velocity for spherical droplets is self-similar – i.e., the curves collapse onto each other – as shown in Fig. 22b. The resulting curve describes the coefficient that relates the non-dimensional slowdown to the droplet diameter. For spherical droplets, this coefficient is not dependent on the free-stream velocity.

Figure 22 Impact velocity for different droplet diameters and free-stream velocities; aerodynamic nose radius R_c=0.07 m and exponent n=1.1; (a) spherical droplets, dimensional velocity; (b) spherical droplets, non-dimensional velocity; (c) oblate droplets, dimensional velocity; (d) oblate droplets, non-dimensional velocity; (e) oblate droplets, droplet semi-major axis.

The results for the deformed droplets, as shown in Fig. 22c, reveal additional effects. First, it can be observed that the impact velocities are noticeably lower. For example, droplets of 1 mm diameter and 90 m s⁻¹ free-stream velocity are slowed down by around 2.5 m s⁻¹ when kept spherical, whereas deformation leads to a slowdown of about 10 m s⁻¹. The reason for this is that the larger surface area due to the deformation leads to higher drag forces, increasing the slowdown for oblate droplets. The impact velocity graphs of the spherical droplets have a concave shape. In the graphs of the oblate droplets, a saddle point appears in the region of 0.5 mm in diameter. The prominence of this saddle point increases with increasing free-stream velocities. From 70 m s⁻¹ the impact velocity does not monotonically increase but shows a slight dip at the saddle point. It is, therefore, possible that a larger droplet has a lower impact velocity. The location of the saddle point coincides approximately with the maximum deformation of the droplet, as shown in Fig. 22e. In this figure, the deformation is shown to rise to a maximum, after which it begins to decline. The maximum corresponds to the diameter at which the limiter of Eq. (21) starts to restrict the growth of the droplets. However, the limiter is not the reason for the occurrence of the saddle points. This can be shown by simulations without a limiter, where the prominence and extent of the saddle point grow. Therefore, the saddle point must be a consequence of the non-linear coupling of the momentum and deformation equation and cannot be attributed to the limiter. It would be interesting to know whether this saddle point can also be observed in experiments. The non-dimensional impact velocity of oblate droplets is self-similar outside the region of the saddle points. In the region of the saddle point, the non-dimensional impact velocities are lower for higher free-stream velocities, indicating that an extra slowdown is obtained that is greater than the common scaling factors of the self-similar solution. It is also evident that with increasing free-stream velocities, the overlap of the curves becomes larger, meaning that, for example, the solutions of 80 and 90 m s⁻¹ are more self-similar than the ones of 50 and 90 m s⁻¹. It can be summarized that oblate droplets slow down more than their spherical peers and that the slowdown effect is sensitive with respect to the droplet diameter.

4.2 The influence of the aerodynamic nose radius on the impact velocity

The influence of the aerodynamic nose radius on the impact speed is investigated in this section for a combination of spherical and oblate droplets of 0.5 and 2.0 mm diameter. Figure 23 shows that 0.5 mm droplets are much more sensitive to a change in R_c than the larger droplets of 2.0 mm. For example, spherical droplets of 0.5 mm diameter have their normalized impact velocity reduced by about 0.1 when R_c is increased from 0.1 to 0.2 m. The impact velocity of the 2.0 mm droplets decreases in the same range by only about 0.01. In general, the curves of the spherical droplet closely overlap, indicating self-similarity. Oblate droplets show much greater sensitivity toward R_c, as seen when comparing Fig. 23c and d. Over the entire range of the investigated nose radii, the velocities of the 2.0 mm spherical droplets decrease by about 0.05, whereas a decrease of approximately 0.25 to 0.3 can be observed for the oblate droplets, i.e., 5 times larger. No self-similarity can be observed for 0.5 mm oblate droplets, as shown in Fig. 23b. The curves of the different free-stream velocities spread out as R_c increases. Here, it is interesting to note that rotating-arm test rigs will not be able to capture this effect due to their small-scale airfoils and, hence, small R_c. Consequently, the non-dimensional slowdown in the test rigs will appear similar for all free-stream velocities, while on an actual turbine it is not. The dependency on the free-stream velocity originates from the non-linear coupling of the momentum and deformation equation as discussed in the previous section; see the result for the 0.5 mm droplets of Fig. 22d. For oblate droplets of 2.0 mm, the curves again overlap closely, as is also the case in Fig. 22d. To conclude, droplets in the saddle point region are especially sensitive to a change in the nose radius. This property is interesting since it means that, especially for faster tip speeds, a higher R_c gives extra slowdown and thus reduces blade damage. Therefore, from a mitigation perspective, it appears to be attractive to utilize aerodynamically thicker airfoils; see Table 4. To summarize, the slowdown effect for oblate droplets is highly sensitive to the aerodynamic nose radius. This sensitivity provides an interesting opportunity for an erosion mitigation strategy.

Figure 23 Non-dimensional droplet impact velocity for different aerodynamic nose radii R_c and free-stream velocities; exponent n=1.1; (a) ϕ₀=0.5 mm, spherical droplets; (b) ϕ₀=0.5 mm, oblate droplets; (c) ϕ₀=2.0 mm, spherical droplets; (d) ϕ₀=2.0 mm, oblate droplets.

4.3 Sensitivity of erosion damage with respect to rain intensity

This section investigates how the droplet slowdown influences the sensitivity of the erosion damage with respect to the rain intensity. First, Eq. (48) is considered, which gives the damage associated with 1 m of impingement at a particular rain intensity. The average droplet impact speed must vary with rain intensity since every rain intensity has a distinct drop-size distribution. As a result, equal amounts of impingement originating from different rain intensities lead to varying degrees of damage. Without the slowdown effect, the impact speed of all droplets, irrespective of their diameter, is equal, and there will be no distinction in damage across the rain intensities. Note that the terminal velocity of a droplet and its dependency on the diameter are neglected here. Equation (48) was brought into a non-dimensional form with

$$\begin{array}{}\text{(54)}& \stackrel{\mathrm{‾}}{\left({\displaystyle \frac{D\left(I\right)}{H\left(I\right)}}\right)}={\displaystyle \frac{D\left(I\right)/H\left(I\right)}{{\left(D\left(I\right)/H\left(I\right)\right)}_{\mathrm{no}\mathrm{slowdown}}}}.\end{array}$$

Figure 24 Normalized erosion damage for 1 m of rain impingement at different rain intensities; default parameters of V_∞=90 m s⁻¹, R_c=0.07 m, and n=1.1; the normalization reference is with respect to no-slowdown droplets.

The three distinct damage exponents β from Sect. 3.4 were considered to establish the robustness of the results with respect to the damage metric. The results are shown in Fig. 24. Droplets without slowdown are non-dimensionalized with themselves and, thus, show damage of unity in the entire plot. The damage for spherical and oblate droplets varies with rain intensity. At low rain intensity, most water mass is contained in the smaller droplets, which experience a significant slowdown, as shown in Fig. 22a and c. Therefore, low-intensity rain shows a large reduction in its damage. As the rain intensity increases, so does the fraction of large droplets within the rain. The large droplets experience considerably less slowdown and, thus, are much more damaging.

Figure 25 Normalized cumulative damage distribution for the De Kooy weather station; R_c=0.07 m, n=1.1, and V_∞=90 m s⁻¹; (a) normalization with respect to no-slowdown droplets – see Eq. (55); (b) normalization with respect to itself – see Eq. (56) (sph.: spherical; obl.: oblate).

Even though the exponents span a wide range, the spherical and oblate droplets' curves remain close together with respect to themselves. The difference in damage between the highest and the lowest exponent is fairly constant for both droplet types across the entire range of rain intensities. This difference is approximately 0.1 and 0.175 for spherical and oblate droplets, respectively. Spherical droplets, especially for smaller rain intensities, already show so much damage reduction that the slowdown effect cannot be neglected. The difference in damage between spherical and oblate droplets is even more significant than between spherical droplets with and without slowdown. Thus it is not sufficient to assume that droplets are spherical: the deformation needs to be taken into account as well. Figure 24 also shows f_I, which is the pdf of the rain intensities. Around 80 % of all precipitation events are of the magnitude 2 mm h⁻¹ and lower. In this range, the slowdown also has the highest effect.

The non-dimensional cumulative site damage reads

$$\begin{array}{}\text{(55)}& {\left(F_{\mathrm{D}}\left(I\right)\right)}_{\mathrm{no}\mathrm{slowdown}}={\displaystyle \frac{D_{\mathrm{cumulative}}}{D_{\mathrm{no}\mathrm{slowdown}}}}.\end{array}$$

It is plotted in Fig. 25a as a function of I for all droplet types and damage exponents. Since the damage is written in its cumulative form, the damage of droplets without slowdown reaches unity for I→∞. The plot shows that for a turbine located at the De Kooy weather station, the inclusion of the droplet slowdown leads to predicted damage of 0.77 to 0.85 for spherical and 0.41 to 0.57 for oblate droplets. Alternatively, expressed in the reciprocal, the predicted lifetime is twice as long for oblate droplets. Figure 25a also shows which rain intensities contribute the most to erosion damage. For example, for droplets without slowdown, all rain events between 0 and 2 mm h⁻¹ contribute to about 55 % of the total erosion damage. From this, the question arises whether the slowdown also influences which rain intensities contribute the most toward erosion damage. To study this, different non-dimensionalization is used:

$$\begin{array}{}\text{(56)}& F_{\mathrm{D}}\left(I\right)={\displaystyle \frac{D_{\mathrm{cumulative}}}{D}}.\end{array}$$

Here every case is non-dimensionalized with itself so that the erosion damage for I→∞ is always unity. Hence, Eq. (56) can also be seen as the cumulative distribution function (cdf) of the damage with respect to the rain intensity. The results are shown in Fig. 25b. For oblate droplets, the 55 % mark of relative damage is shifted to around 2.3 mm h⁻¹ compared to 2 mm h⁻¹ for the case without slowdown. This shows that the slowdown effect not only significantly reduces the predicted erosion damage but also slightly shifts the production of erosion damage to higher rain intensities.

The shift in production of the erosion damage could also influence the viability of erosion mitigation strategies such as the erosion-safe mode. The erosion-safe mode aims at avoiding damage by either reducing the tip speed or shutting down the turbine during precipitation events. To develop this point further, the damage of Eq. (56) can be expressed as $(\mathrm{1}-F_{\mathrm{D}}(I\left)\right)\cdot \mathrm{100}$ % and be plotted against

$$\begin{array}{}\text{(57)}& \left(\mathrm{1}-F_{\mathrm{I}}\left(I\right)\right)\cdot \mathrm{100}\mathit{\%}=\left(\mathrm{1}-\underset{\mathrm{0}}{\overset{I}{\int }}{f}_{I^{\prime }}\left(I^{\prime }\right)\mathrm{d}{I}^{\prime }\right)\cdot \mathrm{100}\mathit{\%},\end{array}$$

resulting in Fig. 26. The figure should be interpreted as how much damage will be saved if X % of the highest-intensity precipitation events can be avoided. As an example, the figure shows that for droplets without slowdown, turning off a turbine during the 20 % most intense precipitation events will reduce the erosion damage by 49 %. Likewise, avoiding the 50 % most intense rain events will save 79 % of all damage. When droplet deformation and slowdown are taken into account, this curve shifts. Depending on the damage exponents avoiding the 20 % most intense rain events now avoids 53 % to 55 % of the erosion damage. Alternatively, when moving laterally, 49 % of erosion damage can be saved when 15.9 % to 17.5 % of the highest-rain-intensity events are avoided. From the figure, it is also visible that the assumption of purely spherical droplets also shifts the curve. However, this shift's magnitude is fairly low compared to oblate droplets. To conclude, the deformation and slowdown effect reduces erosion damage and impacts the viability of erosion mitigation strategies. When the erosion-safe mode is used, neglecting the slowdown effect will yield a sub-optimal utilization by reducing power production in conditions that do not contribute the most toward erosion damage.

Figure 26 Savings in the damage distribution against the X % of heaviest-rain events; R_c=0.07 m, n=1.1, and V_∞=90 m s⁻¹.

4.4 Droplet behavior for reference turbines

The impact of the droplet slowdown along the blades of two reference turbines is investigated. As previously discussed, the NREL 5 MW and IEA 15 MW turbines were chosen for this purpose. The turbines were assumed to be located at the De Kooy weather station. First, the slowdown velocities are analyzed, and the resulting normalized damage distribution is subsequently investigated. Nominal turbine operating conditions at design tip-speed ratio (TSR) were chosen as the control set point for the comparison (IEA TSR=9; NREL TSR=7.55). The parameters from Figs. 13–15 were used for the blade elements. The ballistic angle-of-attack correction coefficients of Table 4 were applied. As explained previously, the philosophies of the original reference turbines were used. This means that the airfoils of the NREL 5 MW turbine stay constant between the officially defined stations, whereas, for the IEA 15 MW turbine, airfoils are linearly interpolated between stations. Hence, a sawtooth pattern is expected in the results of the NREL turbine.

Figure 27 shows the slowdown along the blades of the reference turbines. The calculations were performed for spherical and oblate droplets and for diameters of 0.5, 1.0, 2.0, and 3.0 mm. The slowdown velocities are approximately twice as high for the IEA turbine. The reasons for this can be found in the slightly higher tip speed of the IEA turbine and the larger aerodynamic nose radius, as shown in Fig. 13. The latter, as discussed in Fig. 23, is a significant driver of the slowdown of droplets. The IEA's aerodynamic nose radius R_c is higher not only due to the larger chord, but also due to airfoils, which have, in general, a higher $R_{c,\mathrm{0}}/c$, as shown in Table 4.

Figure 27 Droplet slowdown along the non-dimensional blade distance of the NREL 5 MW and IEA 15 MW turbines; spherical and oblate droplets are considered; slowdown is shown for droplets of varying sizes; (a) NREL 5 MW, spherical droplets; (b) NREL 5 MW, oblate droplets; (c) IEA 15 MW, spherical droplets; (d) IEA 15 MW, oblate droplets.

As expected, smaller droplets show a more significant slowdown along the blade. Spherical droplets experience a decrease in the slowdown velocity from inboard to outboard. This is, at first glance, counterintuitive since the blade element speed is higher toward the tip of the blade. However, the decrease in the aerodynamic nose radius and the aerodynamic exponent, as shown in Fig. 13, offsets the increase in blade element velocity. Oblate droplets show an inverted behavior where the slowdown velocities increase to a maximum when traveling outboard. There, the slowdown effect starts to diminish again. As with the spherical droplets, there is a sharp drop at the blade's tip. In general, in the tip region, the slowdown velocities for oblate droplets are about 2 to 3 times higher than for spherical droplets. Outboard, the higher free-stream velocities must promote the deformation of droplets, leading to a larger slowdown. The sharp drop in slowdown at the very tip of the blade can be explained by the rapidly decreasing chord. To conclude, the deformation and resulting slowdown of the droplets are also critical when actual wind turbines are considered.

The point of maximum slowdown for oblate droplets shifts outboard with decreasing droplet diameters. Larger droplets see their maximum slowdown inboard, whereas the smaller droplets see their maximum outboard of the blade. This reveals another drop-size-dependent non-linearity of the slowdown effect. Larger droplets see a reduced slowdown compared to their smaller peers, and the slowdown is unevenly distributed along the blade. Inboard, large droplets see a relatively large slowdown, whereas small droplets are slowed down significantly in the erosion-prone outboard region of the blade.

An interesting observation can be made with reference to Fig. 27d, where the curves of the various droplet sizes are not only offset but also briefly overlap, e.g., at $r/R_{\mathrm{blade}}=\mathrm{0.65}$ for the 0.5 and 1.0 mm diameter droplets. Even though the droplets have different sizes, they see the same absolute slowdown. This effect was found before, as seen in Fig. 22c, where a saddle point is observed. The position of the saddle point with respect to the droplet diameter shifts for variations in R_c and n and thus leads to different overlapping points along the blade.

Figure 28 shows the non-dimensional damage along the blade. The damage was calculated using Eq. (46) with the non-dimensionalization of

$$\begin{array}{}\text{(58)}& \stackrel{\mathrm{‾}}{D}={\displaystyle \frac{D}{D_{\mathrm{no}\mathrm{slowdown}}}}.\end{array}$$

The damage was calculated for every blade section with V_collection=V_sec. As before, to investigate the sensitivity of the results, the three damage exponents of 5.7, 7, and 9.58 were considered. A damage of unity represents the damage accumulated from a turbine without any droplet slowdown.

For both turbines, the damage decreases toward the blade root, which, at first glance, seems counterintuitive. However, the slowdown velocities stay reasonably constant along the entire blade. In contrast, the blade section speeds vary linearly from close to zero to 82 and 95 m s⁻¹ for the NREL 5 MW and IEA 15 MW turbines, respectively, when moving toward the blade's tip. Hence, the ratio between slowdown and blade element speed is much higher inboard of the blade, and, therefore, the slowdown leads, inboard, to a proportionally higher damage reduction. Still, at the blade's tip, the slowdown effect is non-negligible. Large damage reductions are observed at an $r/R_{\mathrm{blade}}$ value of 0.9. Under the assumption of spherical droplets, the normalized damage is in the range of 0.82 to 0.9 for the NREL turbine. The range for oblate droplets is 0.53 to 0.7. The IEA turbine shows slightly lower non-dimensional damage. As in Fig. 24, the band formed by the damage exponents is fairly constant along the entire blade span, indicating that the results are robust with respect to the damage exponent.

To conclude, the slowdown effect significantly impacts the lifetime prediction of actual wind turbine blades. Adding droplet deformation changes the magnitude and the characteristics of the slowdown velocity along the blade. Even though the highest damage reduction can be found inboard, the slowdown effect remains significant at the blade tip. The results of Figs. 27 and 28 show how a larger R_c can effectively increase the slowdown and thus mitigate erosion damage. This lever seems especially interesting by considering the properties of the airfoils shown in Table 4, i.e., $R_{c,\mathrm{0}}/c$ and the angle-of-attack correction.

Figure 28 Damage distribution along the non-dimensional blade distance: (a) NREL 5 MW and (b) IEA 15 MW.

5 Conclusions

Based on previous findings in the literature, it can be said that experiments in a rotating-arm test-rig environment that used a parameter space relevant to current wind turbine designs have shown that droplets slow down and break up when they approach an airfoil. Hence, slowdown and deformation are also most likely occurring on actual wind turbines. Measurements have shown that the slowdown can be in excess of 10 m s⁻¹ for small droplets. The slowdown becomes less significant as the droplet diameter increases. Moreover, in the above experiments, the breakup modes of bag, bag–stamen, and shear were observed. The role of such droplet breakup in rain erosion is unknown.

From the results obtained in this study, the following main conclusions can be drawn:

• The slowdown effect leads to significant damage reductions and, consequently, should not be neglected in erosion damage modeling. On actual wind turbines, the slowdown effect varies along the blade but remains significant throughout the erosion-prone region. The conclusions regarding the slowdown in this work are robust with respect to variations in the model parameters, such as the exponents of the damage law.

• Droplet size matters! For the investigated cases, droplets under 0.25 mm diameter are slowed down so much that they contribute only marginally to the erosion damage. Large droplets are thus more damaging than their smaller peers. Furthermore, the droplet slowdown is highly sensitive toward the aerodynamics nose radius R_c. Due to an expected difference in trajectory between small and large droplets, the angle-of-attack correction of R_c(α) is projected to be more significant for smaller droplets. This correction increases the slowdown of smaller droplets.

• Droplet shape matters too! The slowdown effect is already significant for spherical droplets. However, the slowdown of oblate droplets greatly exceeds that of spherical droplets. Therefore, deformation must be taken into account. When studying the impact of droplets on blades, droplets should (at least) be modeled as being oblate. Figure 22e can be used as a suggestion for a particular shape.

• Rain intensity matters! This is due to the relationship of rain intensity and droplet size distribution. The slowdown effect is particularly significant for light-rain-intensity events. It also shifts the damage accumulation to higher precipitation intensities. Therefore, it may be beneficial to reduce the tip speed of turbines only during heavy-precipitation events to avoid erosion.

Due to the importance of the droplet slowdown effect on the erosion lifetime of the wind turbine blades, additional research is recommended:

1. Rotating-arm erosion test rigs might also encounter a slowdown effect. This effect would then need to be taken into account in order to find the true impact speed for a given free-stream velocity.

2. It is conceivable that droplets might break up in a cascade decay. Additionally, it has been shown that droplets prior to impact can represent a water mass that has a non-homogeneous velocity. The potential implications for the rain erosion damage of these two effects need to be better understood.

3. In general, more research needs to be conducted on the dynamics of droplet breakup when droplets are subjected to a transient slip velocity field, for example, when approaching an airfoil. Especially the exact conditions and non-dimensional numbers that promote the various breakup modes need to be further understood. Based on such findings, a catalog of droplet shapes just prior to impact would be beneficial, as it could be used in further studies that concern the collision of droplets with wind turbine blades as well as applications beyond wind energy.

Appendix A: List of symbols

Aa	frontal area of deformed droplet
Aprojected	projected droplet area in Springer model
As	surface area of oblate spheroid
As,0	surface area of sphere
A0	projected area of undeformed droplet
a	semi-major axis of deformed droplet
aimpact	droplet semi-major axis on impact
amax	droplet semi-major axis at onset of breakup
a0	initial droplet semi-major axis
$\stackrel{\mathrm{‾}}{a}$	non-dimensional droplet semi-major axis
aϕ	droplet semi-major axis in terminal falling
	conditions
b	droplet semi-minor axis
CD	resultant drag coefficient of droplet
$C_{{\mathrm{D}}_{\mathrm{disk}}}$	drag coefficient of disk
$C_{{\mathrm{D}}_{\mathrm{sphere}}}$	drag coefficient of sphere
Cdynamic	dynamic drag component of droplet
Cp	pressure coefficient driving droplet
	deformation
Cstatic	static drag component of droplet
C1	linear term constant of curve fit for
	angle-of-attack correction
C2	quadratic term constant of curve fit for
	angle-of-attack correction
c	airfoil chord
D	total erosion damage (during a year)
D(I)	D as a function of rain intensity
$\stackrel{\mathrm{‾}}{D}$	non-dimensional D
Dcumulative	D in cumulative form

Dno slowdown	D for when droplets do not have
	slowdown
F	impact force in Springer model
FD(I)	cdf of fD(I)
Fdrag	drag force on droplet
FI	cdf of the rain intensity
Fp	pressure force on droplet driving
	deformation
Fσ	surface tension force on droplet
Fϕ,site	cdf of water of mass associated with a
	droplet diameter at a particular site
fI	pdf of rain intensity
fϕ,air	drop-size distribution in a volume of air
fϕ,plane	drop-size distribution through plane in
	air
fϕ,site	pdf of water of mass associated with a
	droplet diameter at a particular site
H	accumulated impingement
H(I)	H as a function of rain intensity
H(t)	accumulated impingement as a function
	of time
Hallowed	allowed impingement until end of
	incubation
h	height in kilometers for droplet terminal
	velocity formula
I	rain intensity
k	calibration constant in DRD model
m	droplet mass
N	measure for incubation time
$N_i^{*}$	number of allowable impacts on one
	location as per Springer model
n	aerodynamic exponent
nLE	normal vector of the leading edge
Oh	Ohnesorge number
prain	probability of rain (at site)
pwh	water hammer pressure in Springer
	model
Rc	aerodynamic nose radius
Rc(α)	aerodynamic nose radius as a function
	of angle of attack
Rc,0	Rc(α) at 0° angle of attack
Re	droplet Reynolds number
R0	initial droplet radius
r	position along rotor blade
r	distance to impact location in Springer
	model
$r/R_{\mathrm{blade}}$	non-dimensional position along rotor
	blade
Train	duration of rain (during a year)
Tyear,spinning	time in a year the turbine spins
t	time
Va	expansion velocity of droplet
	semi-major axis
Va,0	starting expansion velocity of droplet
	semi-major axis

Vair	velocity of air (in front of the airfoil)
Vblade	blade velocity
Vcollection	water collection velocity for impingement
Vimpact	impact speed of droplet
Vrain	rain velocity vector
Vsec	blade section speed
Vsec	blade section velocity vector
Vslip	velocity difference between droplet and
	airflow
Vslowdown	slowdown velocity
Vx	droplet velocity
Vx,0	droplet velocity at start of slowdown
	simulation
V∞	ΔV far away from the airfoil
Vϕ	droplet terminal velocity
W	volume of water per volume of air
Wϕ	W for a particular drop size
We	Weber number
Weimpact	Weber number based on Vimpact
We∞	Weber number based on V∞
x	droplet position
xb	ballistic droplet trajectory position
xblade	position of the blade in slowdown
	simulation
xblade,0	starting position of the blade in slowdown
	simulation
xs	streamline droplet trajectory position
x0	starting position of droplet
x95 %	position at which ΔV reaches 95 % of V∞
x95 %(α)	x95 % as a function of the angle of attack
x95 %,0	x95 % at 0° angle of attack
α	angle of attack
α	coefficient in damage law
β	exponent in damage law
ΔV	relative velocity between droplet and airfoil
Δx	distance between droplet and airfoil
ϵ	eccentricity of droplet
μ	mean in lognormal distribution
μwater	dynamic viscosity of water
μair	dynamic viscosity of air
ρair	density of air
ρwater	density of water
σ	standard deviation in lognormal distribution
σwater	surface tension of water
σ(…)	stress field in Springer model derivation
ϕ	droplet diameter
ϕ0	initial droplet diameter
∂tH	impingement accumulation rate