An Overview of Power Loss Estimation in Wind Turbines Due to Icing

Abstract

Wind turbines are susceptible to severe meteorological conditions, which can result in power loss. Several methods have been proposed to estimate the extent of power loss in wind turbines. This chapter aims to establish a foundation for new research and investigations into the impact of icing on wind turbine power output. It provides an overview of various methodologies available for estimating power loss in wind turbines under icing conditions. One of the prominent methods utilized in the past decade is computational fluid dynamics (CFD), enabling three-dimensional numerical simulations of wind turbines. When combined with the blade element momentum theory (BEM), CFD can also facilitate two-dimensional simulations. By analyzing these methodologies, researchers can gain insights into the estimation techniques suitable for studying icing effects on wind turbine performance. Understanding the magnitude of power loss under icing conditions is crucial for optimizing wind turbine design, operation, and maintenance strategies. Overall, this chapter contributes to the body of knowledge by consolidating various methods employed for power loss estimation in wind turbines subjected to icing.

1. Introduction

Wind energy has emerged as one of the most promising and sustainable sources of electricity generation in recent years. Wind turbines, comprising large rotating blades, convert the kinetic energy of wind into electrical energy. However, wind turbines are subject to various operational challenges, one of which is the formation of ice on the turbine blades during cold and humid conditions, especially in Nordic countries characterized by higher wind potential. It is predicted that the installed wind power capacity in cold climates will increase to 224 GW in 2025 from the reported power of 156 GW in 2020 [1]. However, icing on wind turbine blades has a significant impact on their performance, leading to reduced power generation and increased maintenance costs [2, 3, 4, 5, 6]. This occurs when supercooled liquid droplets or freezing rain come into contact with the blades and freeze upon impact. This ice accumulation alters the aerodynamic properties of the blades, resulting in decreased lift and increased drag forces. As a consequence, the turbine experiences reduced rotational speed and power output.

The phenomenon of power loss in wind turbines due to icing has garnered substantial attention in recent years, as it poses a significant challenge for wind farm operators and renewable energy stakeholders. Accurate estimation of the power loss caused by icing is crucial for efficient wind turbine operation, maintenance planning, and optimizing overall power generation. The impact of icing losses on the yearly power output and aerodynamic performance of the wind turbine is reported in [2, 7] respectively.

Traditionally, the assessment of ice accretion effects on wind turbines has been limited to two-dimensional analyses, focusing on a simplified representation of an airfoil section of the blade. However, horizontal-axis wind turbines (HAWTs) usually have a twisted blade geometry comprising of multiple section forms [8]. To accurately capture the impact of icing, three-dimensional simulations are necessary to account for the intricate blade geometry. It’s important to note that icing effects extend beyond the blade itself, affecting the entire wind turbine. Various three-dimensional phenomena come into play, such as the influence of airflow along the radial direction induced by the rotation of the blades. Considering these three-dimensional effects is crucial for a comprehensive understanding of the icing’s impact on the wind turbine system.

Several factors influence the magnitude of power loss due to icing in wind turbines. The severity of icing is primarily determined by environmental conditions, including air temperature, humidity, and wind speed. Other factors such as the design and characteristics of turbine blades, the type of ice formed, and the operating conditions of the turbine also play a role in determining the extent of power loss. Estimating the power loss due to icing is a complex task that requires a multidisciplinary approach. Researchers and engineers have employed various methodologies to quantify the impact of icing on wind turbine performance. Experimental investigations, computational fluid dynamics (CFD) simulations, and field measurements have been utilized to understand the aerodynamic effects of ice accretion and its subsequent influence on turbine power production.

This chapter focuses on the current state-of-the-art in estimating power loss in wind turbines under icing conditions using numerical methods, particularly computational fluid dynamics (CFD). It aims to provide a comprehensive understanding of the modeling of icing phenomenon and the analysis of wind turbine performance. The objective is to present the available methods for quantifying power loss in wind turbines when subjected to icing conditions. The study also recognizes the significance of comprehending and addressing the consequences of icing on wind turbines, considering its considerable influence on power generation, structural stability, and operational safety. The primary objective of this research is to enhance our understanding of icing phenomena by examining the different factors involved and investigating methodologies for estimating power losses in wind turbines when subjected to icing conditions. Through an analysis of existing research and identification of research gaps, this study aims to facilitate the advancement of improved icing prediction systems and effective mitigation strategies within the wind energy industry.

The chapter is structured into several sections: Section 2 explains the CFD approach used for simulating ice accretion, Section 3 discusses the Blade Element Momentum (BEM) theory for aerodynamic analysis, Section 4 reviews the CFD-BEM and three dimensional (3D) CFD approaches, and finally, Section 5 presents the main conclusions and future research directions.

2. CFD approach used for icing conditions in wind turbines

CFD-based models for icing conditions in wind turbines involve using techniques to simulate and analyze the impact of ice accretion on wind turbine performance. These models aim to provide insights into the complex interactions between airflow, water droplets, ice formation, and the resulting effects on turbine aerodynamics and power production. CFD-based models for icing conditions in wind turbines provide valuable insights into the complex phenomena associated with ice formation and its impact on turbine performance. These models aid in understanding the icing risks, optimizing turbine design, and developing effective strategies to mitigate the adverse effects of icing in wind energy applications. As shown in [8, 9], the numerical simulation of the CFD-based model can be described by fluid flow simulation, surface thermodynamics, droplet behavior and phase changes.

2.1 Fluid flow simulation

Fluid flow simulation in CFD-based models for icing conditions in wind turbines involves simulating the airflow around the turbine components to understand the complex flow patterns and their interaction with ice formation. By simulating the fluid flow around wind turbines using CFD-based models, engineers and researchers can gain valuable insights into the complex flow phenomena associated with icing conditions. This information is crucial for optimizing turbine designs, assessing ice-related risks, and developing mitigation strategies to enhance the safety and performance of wind energy systems.

The fundamental equations that describe fluid flow are the Navier-Stokes equation [10] which are used to solve for the velocity, pressure, and temperature fields in the flow domain. In the case of steady-state simulations, the equations simplify to the Reynolds-averaged Navier-Stokes (RANS) equations. In their general form, the Navier-Stokes equations can be written as follows:

Conservation of mass (continuity equation):

$$\mathrm{\partial \rho }/\mathrm{\partial t}+\text{∇}·\left(\mathrm{\rho u}\right)=0$$E1

where ρ is the density of the fluid, t is time, u is the velocity vector, and ∇ represents the divergence operator.

Conservation of momentum (Navier-Stokes equations):

$$\partial \left(\mathrm{\rho u}\right)/\mathrm{\partial t}+\text{∇}·\left(\mathrm{\rho u}\otimes \mathrm{u}\right)=-\text{∇p}+\text{∇}·\mathrm{\tau }+\mathrm{\rho g}$$E2

where ∂(ρu)/∂t is the rate of change of momentum, ⊗ represents the outer product (tensor product), p is the pressure, τ is the stress tensor, and g is the gravitational acceleration vector.

Conservation of energy:

The energy equation accounts for the conservation of thermal energy and can be written as:

$$\partial \left(\mathrm{\rho e}\right)/\mathrm{\partial t}+\text{∇}·\left(\mathrm{\rho e}+\mathrm{p}\right)\mathrm{u}=\text{∇}·\left(\mathrm{k}\nabla \mathrm{T}\right)+\text{∇}·\left(\text{μeff}\nabla \mathrm{u}\right)+\mathrm{\rho u}·\mathrm{g}$$E3

where ∂(ρe)/∂t is the rate of change of energy, e is the internal energy per unit mass, p is the pressure, k is the thermal conductivity of the fluid, T is the temperature, μeff is the effective viscosity, ∇T is the gradient of temperature, and u · g represents the dot product of velocity and gravitational acceleration.

These equations, coupled with appropriate boundary conditions, form the basis for solving fluid flow problems using the Navier-Stokes equations. It’s worth noting that the equations can be further simplified and specialized depending on the specific assumptions and modeling techniques used in different scenarios. Appropriate boundary conditions need to be specified to simulate the icing conditions accurately. This includes specifying the inlet conditions (velocity, temperature, and turbulence properties), the airfoil/wind turbine blade surface conditions, and the icing conditions (e.g., supercooled droplet size, liquid water content). An icing model is incorporated to simulate the growth and accumulation of ice on the wind turbine blades. These models typically consider the impingement of supercooled water droplets, heat transfer, phase change, and ice accretion processes. Different icing models exist, such as the Messinger model or the Leontiev model, which can be implemented based on the specific requirements of the simulation.

Since turbulence plays a crucial role in wind turbine flows, the choice of turbulence model depends on factors such as the level of accuracy required, computational resources available, and the specific flow characteristics being simulated. The Reynolds-Averaged Navier-Stokes (RANS) models and Large Eddy Simulation (LES) models remain the two commonly used turbulence models [9]. In the RANS model, the governing equations of fluid motion, known as the Navier-Stokes equations described above, are averaged in time. This averaging process results in the set of equations that describe the mean flow behavior. The time-averaged equations are then solved numerically to obtain the mean flow field.

The RANS equations are expressed in terms of the mean velocity components, pressure, and other flow properties. These equations include the conservation of mass, momentum, and energy, similar to the original Navier-Stokes equations. However, in the RANS model, additional terms called Reynolds stresses appear in the momentum equations to account for the effect of turbulent fluctuations on the mean flow.

To close the RANS equations, a turbulence model is needed to approximate the Reynolds stresses. The turbulence model provides a closure relation for the Reynolds stresses based on the mean flow variables. The most commonly used turbulence model in RANS simulations is the eddy viscosity model, where the Reynolds stresses are related to the mean velocity gradients through an eddy viscosity coefficient.

The eddy viscosity coefficient is usually determined by solving an additional transport equation called the turbulence model equation. This equation accounts for the evolution of turbulence and provides a link between the mean flow and the turbulent fluctuations. Different turbulence models, such as the k-epsilon model [2, 11], the Spalart-Allmaras model [12], k-omega models [13] and the Reynolds Stress models (RSM) [14] have all been developed with varying levels of complexity and accuracy to capture different aspects of turbulence.

The k-epsilon models are based on the two-equation approach, which includes transport equations for turbulent kinetic energy (k) and its dissipation rate (epsilon). The eddy viscosity is computed using the turbulent kinetic energy and the turbulent length scale. Popular variations include the Standard k-epsilon model and the RNG k-epsilon model, which employ different assumptions and closure coefficients. This model is quite popular because of its relatively steady, good, and rapid convergence rate as well as low computational cost. However, in some cases, it is not adequate in flows with high-pressure gradients.

Similar to k-epsilon models, k-omega models are also two-equation models. They involve transport equations for turbulent kinetic energy (k) and the specific dissipation rate (omega). The eddy viscosity is calculated using the turbulent kinetic energy and the turbulent frequency, which is related to the specific dissipation rate.

RSMs consider the full Reynolds stress tensor in the turbulent flow equations. These models involve solving transport equations for the six components of the Reynolds stress tensor, and the eddy viscosity is computed from the resolved Reynolds stresses.

The Spalart-Allmaras model is a one-equation model which employs a transport equation for the eddy viscosity, which is determined based on the turbulent kinetic energy and a modified turbulent length scale. The model introduces a turbulent viscosity that adapts to the flow’s nature, including near-wall effects and low-Reynolds-number flows. Due to its ability to strike a balance between computational efficiency and desired precision in turbulent flow simulations, this turbulence model finds widespread use in the field of aerodynamics. Moreover, numerous researchers have acclaimed this model as the top-performing choice for accurately simulating ice accumulation on electrical cables, wind turbines, and aircraft.

It’s worth noting that these models have their own advantages and limitations, and their performance varies depending on the flow characteristics and specific applications. The choice of a suitable turbulence model depends on factors such as flow regime, geometrical complexity, and available experimental or reference data for validation. Once the RANS equations and the turbulence model are defined, they can be solved numerically using appropriate numerical methods, such as finite difference, finite volume, or finite element methods. The resulting simulations provide information about the mean flow field, including velocity profiles, pressure distribution, and other important flow characteristics.

The LES models are more computationally demanding but provide a higher level of turbulence resolution. Unlike the Reynolds-averaged Navier-Stokes (RANS) model that relies on averaging turbulent fluctuations, LES focuses on resolving large-scale turbulent structures while modeling the effect of small-scale fluctuations. In LES, the flow is divided into large eddies and small-scale turbulence. The large eddies, which contain the most significant energy-containing structures, are resolved directly in the simulation, while the small-scale turbulence is modeled. This approach allows for a more accurate representation of the flow physics compared to RANS models. The LES equations consist of the filtered Navier-Stokes equations, where a filter operation is applied to the flow variables to separate the resolved scales from the unresolved ones. This filter operation effectively removes the small-scale fluctuations, and the resulting equations are then solved numerically to simulate the resolved flow field. To model the unresolved scales, a subgrid-scale (SGS) model is employed. The SGS model represents the effect of small-scale turbulence on the resolved flow by providing closure relations for the filtered variables. These closure relations are typically based on assumptions and modeling approaches, such as the Smagorinsky model, which relates the SGS stress tensor to the strain rate of the resolved flow.

The LES approach requires resolving a sufficient range of turbulent scales to capture the important features of the flow. This is achieved by using a fine computational mesh that can resolve the larger eddies while employing suitable numerical schemes to capture the details of the flow dynamics. LES is particularly useful for studying flows with significant turbulence and large-scale structures, such as turbulent boundary layers, swirling flows, and turbulent wakes. It provides more accurate predictions compared to RANS models and is often used in engineering applications where capturing detailed turbulence dynamics is critical. However, LES simulations are computationally expensive due to the high resolution required to resolve the large eddies. Therefore, LES is generally used for research purposes or in cases where the detailed understanding of the turbulence physics is essential. In practical engineering simulations, hybrid approaches combining RANS and LES, such as detached eddy simulation (DES), are often employed to strike a balance between accuracy and computational cost.

Both the RANS and LES models introduce additional equations to account for turbulence quantities such as turbulent kinetic energy and turbulent dissipation rate. The computational domain is discretized into a grid or mesh of cells or elements. This discretization can be structured (e.g., Cartesian) or unstructured (e.g., tetrahedral or hexahedral elements). The governing equations are then solved numerically on this grid using methods such as finite volume, finite element, or finite difference techniques. The discretized equations are solved iteratively using numerical algorithms. The solution procedure involves applying appropriate solvers to compute the flow field variables, such as pressure, velocity, and temperature. The solution is obtained by iterating until convergence is achieved. Once the simulation is complete, post-processing is carried out to analyze the results. This may involve visualizing the flow patterns, calculating performance metrics (e.g., power output, thrust), and evaluating the ice accretion on the wind turbine blades.

Different turbulence models can lead to discrepancies in predicting turbulent flows due to their inherent assumptions and limitations. Some of the discrepancies that can arise from turbulence models include but not limited to boundary layer flows, complex geometries, compressible flows and high Reynolds number flows.

The k-epsilon model is known to struggle with accurately predicting near-wall flows, such as boundary layers. It can overpredict the turbulent kinetic energy near walls, leading to incorrect velocity profiles and flow separation predictions. Reynolds stress models, on the other hand, offer improved predictions in boundary layers by explicitly considering the Reynolds stress terms. They capture the anisotropic nature of turbulence and can provide more accurate velocity profiles and boundary layer characteristics.

Turbulence models like the Spalart-Allmaras model may have limitations in predicting separated flows accurately. They might fail to capture the complex behavior near separation regions, resulting in inaccurate pressure distributions and flow reattachment lengths. Reynolds stress models, with their ability to account for the anisotropy of turbulence, can offer better predictions of separated flows, including reattachment locations and flow structures.

Turbulence models, including k-epsilon, k-omega, and Spalart-Allmaras, may struggle in simulating flows around complex geometries, such as bluff bodies or turbomachinery. They might provide less accurate predictions of flow separation, pressure losses, and flow recirculation regions. Reynolds stress models, with their explicit modeling of Reynolds stress terms, have the potential to better capture complex flow features and flow physics in such scenarios.

The k-epsilon model and Spalart-Allmaras model can be limited in accurately predicting compressible flows, especially in regions with shock waves or strong density gradients. They might fail to capture the complex interactions between turbulence and compressibility effects. Reynolds stress models and some k-omega models, which explicitly account for the Reynolds stress terms, can better handle compressible flows and provide improved predictions of shock-boundary layer interactions.

2.2 Surface thermodynamics

Surface thermodynamics deals with the heat transfer processes occurring at the surface of the wind turbine blades. In the presence of icing conditions, the transfer of heat between the blade surface and the surrounding air is crucial. The surface temperature distribution affects ice accretion and subsequent ice behavior. Modeling surface thermodynamics involves considering factors such as conduction, convection, and radiation heat transfer mechanisms. By accurately modeling surface thermodynamics, the simulation can provide insights into the temperature distribution on the wind turbine blades.

Heat conduction occurs when there is a temperature gradient within a solid material. In the case of wind turbine blades, heat is conducted from the interior of the blade to the surface. The thermal conductivity of the blade material determines how efficiently heat is transferred. It is important to accurately model the heat conduction within the blade, considering its composition and any internal temperature variations that may arise due to operational conditions or heating elements. The heat conduction equation, based on Fourier’s law, governs the conduction process:

$$\mathrm{q}=-\mathrm{k}\text{∇T}$$E4

where q is the heat flux vector, k is the thermal conductivity, and ∇T is the temperature gradient. By solving this equation, the temperature distribution within the blade can be determined, which directly affects the surface temperature.

Convection refers to the transfer of heat between a solid surface and a moving fluid (in this case, air). In icing simulations, convective heat transfer occurs between the wind turbine blade surface and the surrounding air, affecting the blade temperature and ice formation. The convective heat transfer coefficient, h, characterizes the rate of heat transfer through convection. It depends on several factors, including the local flow velocity, air properties (density, viscosity, and specific heat capacity), blade surface geometry, and roughness. Determining the convective heat transfer coefficient can be challenging as it relies on accurately capturing the complex flow patterns around the blade surface. Empirical correlations, such as those based on Nusselt number relationships, or more advanced CFD methods can be employed to estimate the convective heat transfer coefficient in the icing simulation.

Radiation heat transfer occurs through electromagnetic waves emitted by a surface. In the context of wind turbine blades, both incoming solar radiation and thermal radiation from the surrounding environment contribute to the heat exchange. The net radiation heat transfer between the blade surface and its surroundings depends on various factors, including the surface emissivity, absorptivity, reflectivity, and the temperature difference between the blade and the surroundings. It’s important to consider the radiative properties of the blade surface and account for both solar and thermal radiation in the icing simulation.

Surface roughness influences the formation and behavior of the boundary layer, which affects the convective heat transfer. Rough surfaces can disrupt the laminar flow and induce turbulence, enhancing heat transfer and potentially altering ice formation patterns. Accurate representation of the blade surface roughness is crucial for realistic icing simulations. Wettability refers to the ability of the surface to repel or retain water. Hydrophobic surfaces can reduce ice adhesion and promote ice shedding. Surface coatings or modifications can be incorporated in the simulation to represent different surface conditions.

There are various physical models available that plays a crucial role in examining the thermodynamic properties of the icing process on wind turbine blades when water droplets impact the blade surface. One pioneering and significant contribution in this area is the Messinger model, dating back to 1953. The Messinger model is a numerical approach based on the first-order differential equations of mass conservation and heat transfer that relates to icing and is based on the temperature of unheated surfaces. The Messinger model serves as the fundamental basis for simulating icing on aircraft wings, however, its application has been extended to wind turbine systems employing a straightforward approach that relies on mass and energy balance principles as shown in Eqs. (5) and (6) respectively.

$${\dot{m}}_{\mathit{\text{in}}}+{\dot{m}}_{\mathit{imp}}={\dot{m}}_{\mathit{ice}}+{\dot{m}}_{e,s}+{\dot{m}}_{\mathit{\text{out}}}$$E5

$${\dot{Q}}_{\mathit{cov}}+{\dot{Q}}_{\mathit{es}}+{\dot{Q}}_{\mathit{imp}}={\dot{Q}}_{\mathit{kin}}+{\dot{Q}}_{\mathit{ice}}+{\dot{Q}}_{\mathit{\text{in}}}$$E6

${\dot{m}}_{\mathit{\text{in}}}$ is the rate of water inflow, ${\dot{m}}_{\mathit{\text{out}}}$ is the rate of water outflow, ${\dot{m}}_{\mathit{imp}}$ is the rate of droplets impingement, ${\dot{m}}_{\mathit{ice}}$ is the ice accretion rate and ${\dot{m}}_{e,s}$ is either the rate of evaporation or the rate of sublimation depending on the type of ice.

The Messinger model performs well for dry icing conditions, where there is low liquid water content (LWC) in the air, and the air temperature is considerably below freezing. However, the model encounters difficulties when dealing with conditions that involve higher LWC and air temperatures closer to freezing. To address this limitation and improve the model’s predictive capabilities under conditions of higher LWC and air temperatures near freezing, an extension to the Messinger model was proposed. This extension involves the incorporation of conduction terms, which helps enhance the model’s accuracy and reliability in such scenarios. By including these conduction terms, the improved model can better handle conditions with higher LWC and air temperatures close to freezing, providing more accurate predictions of ice formation on structures in these challenging conditions. Ice formation relies on the extended Messinger model, which is a conventional approach utilizing differential equations to describe phase change phenomena, known as the Stefan problem. This thermodynamic model relies on a phase transition phenomenon, commonly referred to as the Stefan problem. The model is controlled by the conduction of heat through both ice and water, as well as considerations of mass balance and the conditions for phase change as shown in the equations below.

$$\frac{\partial T}{\partial t}=\frac{k_i}{{\rho }_i{C}_{\mathit{pi}}}\frac{{\partial }^{2}T}{\partial y^2}$$E7

$$\frac{\partial \theta }{\partial t}=\frac{k_w}{{\rho }_w{C}_{\mathit{pw}}}\frac{{\partial }^2\theta }{\partial y^2}$$E8

$${\rho }_i\frac{\partial B}{\partial t}+{\rho }_w\frac{\partial h}{\partial t}={\rho }_a\beta V_{\infty }+{\dot{m}}_{\mathit{\text{in}}}-{\dot{m}}_{e,s}$$E9

$${\rho }_i{L}_F\frac{\partial B}{\partial t}=k_i\frac{\partial T}{\partial y}-k_w\frac{\partial \theta }{\partial y}$$E10

where θ and T are the temperatures, k_w and k_i are the thermal conductivities, C_pw and C_pi are the specific heats, h and B are the thicknesses of water and ice layers, respectively.

To obtain the ice and water thicknesses and the temperature distribution at each layer, specific boundary and initial conditions need to be defined. These conditions are established on the basis of the following assumptions:

1. The ice is assumed to be in perfect contact with the airfoil surface, and its temperature is considered equal to the air temperature.

2. The temperature remains continuous at the boundary between ice and water, where it is set to the freezing temperature.

3. At the interface between air and water (glaze ice) or air and ice (rime ice), the heat flux is influenced by various factors, including convection, radiation, latent heat release, cooling caused by incoming droplets, heat carried by runback water, evaporation or sublimation, aerodynamic heating, and the kinetic energy of incoming droplets.

4. Airfoil surface is initially clean

Another frequently utilized model is the Makkonen model, which formulates a differential equation representing the gradual accumulation of ice over time. This model incorporates the LWC in the equation to characterize the ice accretion process. The Makkonen model readily facilitates the simulation of ice load and ice accretion rate. Typically, ice load is expressed in kg/m, while the ice accretion rate is given in kg/m/h. The latter measurement pertains to each meter of the cylinder employed in the model. Alternatively, the ice accretion rate may be presented in terms of active ice hours, denoting the hours when the ice accretion rate exceeds a predefined threshold, often set at 10 g/m/h. This model relies on three fractions: collision efficiency, sticking efficiency, and accretion efficiency. Makkonen model variants in combination with numerical weather prediction (NWP) can predict icing and production losses on wind turbines.

The challenge associated with employing physical models is the absence of measured values for parameters like liquid water content (LWC) and median volume diameter (MVD). As a result, Numerical Weather Prediction (NWP) is utilized to estimate these values. However, this approach has a limitation in that these estimated values are typically available for the entire wind park, and not at a per wind turbine level. Though the Messinger and Makkonen models can both used to simulate ice accretion in wind turbines, they may have discrepancies in their approaches and predictive capabilities. While the Messinger model is primarily based on first-order differential equations of mass conservation and heat transfer. It focuses on the energy required for icing protection and the freezing fraction of impinging water. On the other hand, the Makkonen model involves a differential equation that describes the ice accretion over time, considering the LWC as a key parameter.

2.3 Droplet behavior

The concept of droplet behavior is important to understand the formation and growth of ice on turbine surfaces. Droplet behavior refers to the motion and interaction of water droplets in the air, which can lead to ice accretion on turbine blades. To simulate droplet behavior in CFD models, several key phenomena need to be considered. These include droplet transport, evaporation, collision/coalescence, and impingement.

Droplet transport involves modeling the motion of individual water droplets in the airflow around the turbine. The transport of droplets can be described using the equations of motion, taking into account the forces acting on the droplets, such as gravity, drag, and lift forces. The Navier–Stokes equations are typically solved to determine the airflow, and droplet trajectories are computed based on the local flow field.

Water droplets can evaporate as they travel through the air due to the difference in vapor pressure between the droplet and its surroundings. The rate of droplet evaporation can be estimated using various empirical or semi-empirical correlations, such as the Reynolds analogy or the Hertz-Knudsen equation. These correlations relate the evaporation rate to factors such as droplet size, temperature, relative humidity, and air velocity.

In regions with a high concentration of droplets, collisions between droplets can occur, leading to coalescence (merging) or breakup. The probability of collision and the resulting droplet behavior depend on factors like droplet size, relative velocity, and the concentration of droplets. These effects are typically modeled using stochastic or statistical methods, such as the Smoluchowski coagulation equation or the kinetic theory of coalescence.

When droplets come into contact with a solid surface, such as wind turbine blades, they can impinge and freeze, leading to ice accretion. The impingement process depends on the droplet size, impact velocity, surface temperature, and surface characteristics. Empirical models or experimental data are often used to determine the droplet impingement efficiency, which represents the fraction of impinging droplets that freeze upon impact.

The droplets’ behavior in wind turbine icing can be described by considering the icing phenomenon as a multiphase flow because the icing rate may be affected by a mixture of ice, air, and liquid water since icing on wind turbine surfaces involves the simultaneous presence of multiple phases [15]. By treating the system as a multiphase flow, it becomes possible to model and understand the behavior of these phases, their interactions, and their effects on the icing process. Multiphase flow analysis allows for the consideration of interactions between the different phases. In the case of wind turbine icing, this includes the interaction between liquid water droplets and the turbine surfaces, the interaction between droplets and the surrounding air, as well as the growth and interaction of ice with the other phases. These interactions play a significant role in determining the rate and characteristics of ice formation. Several models are used to study droplet behavior in wind turbine icing. These models aim to simulate and understand the complex interactions between droplets, ice, and the surrounding air. The Lagrangian Particle Tracking model (LPT) and the Eulerian model remain the most used models for simulating disperse phase present in multiphase flows [16].

The Lagrangian Particle Tracking Model is a computational approach that tracks individual droplets as Lagrangian particles in the flow field surrounding a wind turbine. It is commonly used to analyze the behavior of liquid water droplets, including their impingement, adhesion, and freezing. The Lagrangian Particle Tracking Model employs Newton’s equations of motion to solve for the movement of each individual droplet or particle. It incorporates a collision model to handle interactions between particles. The trajectory calculation of the particles involves two steps: (1) determining the particle motion based on Newton’s second law and the forces acting on the particle, excluding collisions, such as gravity and the fluid phase’s traction force on the particle. (2) Considering the particle’s collision with another particle. The utilization of this method is commonly seen in 2D icing analysis. To assess the icing rate, the velocities of water droplets are determined by numerically tracking a few droplets near the object experiencing icing. However, tracking a sufficient number of droplets in 3D to create the shape of the ice layer becomes computationally expensive [2].

In the context of wind turbine icing, an essential aspect of the LPT model is the consideration of droplet deposition on the turbine surfaces. Deposition models are used to determine whether droplets adhere to the surface based on criteria such as impact velocity and surface conditions. Additionally, droplet evaporation models may be incorporated to account for the reduction in droplet size due to the heat and mass transfer between the droplets and the surrounding air. To determine the forces acting on the droplets, the LPT model requires information about the fluid flow field. Typically, this information is obtained from a separate Computational Fluid Dynamics (CFD) simulation or experimental data. The fluid flow parameters, such as velocity, pressure, and turbulence characteristics, are interpolated from the CFD data to the droplet locations and times during the Lagrangian integration.

The Eulerian model, on the other hand, focuses on describing the behavior of the continuous fluid phase [17]. In this model, the fluid flow field is divided into a grid or mesh, and the governing equations for mass, momentum, and energy conservation are solved for each grid cell. The Eulerian model provides a detailed understanding of the fluid flow characteristics, such as velocity, pressure, and turbulence, which are critical for the behavior of the dispersed phase. Within the Eulerian framework, different sub-models are used to represent the dispersed phase, such as the volume fraction or concentration fields, and additional transport equations are solved to account for the interactions between the dispersed phase and the continuous phase.

In the Eulerian model, the flow domain is discretized into a structured or unstructured grid composed of cells or control volumes. The grid is a spatial representation of the domain where the fluid flow is simulated. Each grid cell represents a small portion of the domain, and the flow properties are computed at the cell centers or vertices. The Navier-Stokes equations, are solved for each grid cell in the Eulerian model. These equations include the conservation of mass (continuity equation) and the conservation of momentum (momentum equations). Additional equations, such as the energy equation, may be solved if heat transfer is significant in the system. These equations are solved numerically using appropriate discretization schemes, such as finite difference, finite volume, or finite element methods.

The Eulerian model calculates the flow variables, such as velocity, pressure, temperature, and turbulence characteristics, at each grid cell in the computational domain. These variables are represented as fields, meaning they are functions of space and time. The values of these variables are determined by solving the governing equations and considering the boundary conditions.

The conservation equations, derived from the Navier-Stokes equations, are solved for each grid cell. The continuity equation represents the conservation of mass, ensuring that mass is conserved within each control volume. The momentum equations represent the conservation of momentum and account for the various forces acting on the fluid, such as pressure gradients, viscous forces, and external body forces. Boundary conditions are essential in the Eulerian model to specify the flow behavior at the domain boundaries. Different types of boundary conditions are applied, such as velocity inlets, pressure outlets, wall boundaries, and symmetry conditions. These boundary conditions provide the necessary information to define the flow characteristics and interactions with the surroundings.

By solving the governing equations on the grid using the Eulerian model, detailed information about the flow field, including velocity profiles, pressure distributions, and temperature variations, can be obtained. The Eulerian model is computationally efficient and well-suited for large-scale simulations of fluid flows, making it a valuable tool in various engineering applications, including aerodynamics, hydrodynamics, and industrial fluid dynamics.

The liquid water content (LWC) and droplet volume diameter have significant effects on the droplet collection efficiency. The collection efficiency refers to the fraction of droplets that are captured or collected by a specific surface or device. The LWC represents the amount of liquid water present in a unit volume of air. It has a direct impact on the droplet collection efficiency. Generally, as the LWC increases, the collection efficiency tends to increase as well. This is because a higher LWC implies a greater number of droplets available for capture, increasing the chances of collisions between the droplets and the capturing surface. The increased droplet concentration enhances the overall collection efficiency. The droplet volume diameter, often referred to as the droplet size, is the measure of the droplet’s size or diameter. It plays a crucial role in determining the droplet collection efficiency. The relationship between droplet size and collection efficiency is complex and depends on various factors, including the capturing mechanism and the size of the capturing surface.

The combined effects of LWC and droplet volume diameter can be complex and interdependent. Higher LWC provides a higher number of droplets, which can compensate for the lower collection efficiency of smaller droplets due to factors like Brownian diffusion. In some cases, an optimal droplet size range may exist where the collection efficiency is maximized.

It’s important to note that the specific capturing mechanism, the characteristics of the capturing surface, and the airflow conditions also influence the droplet collection efficiency. Additionally, the interaction between droplets, such as coalescence or breakup, can further complicate the behavior of droplet collection. Understanding the effects of LWC and droplet volume diameter on droplet collection efficiency is crucial for various applications, including aerosol sampling, wet scrubbers, and precipitation estimation in atmospheric science. Experimental measurements, numerical simulations, and theoretical models are often employed to investigate and quantify these effects under specific conditions and configurations.

2.4 Phase changes

Phase changes play a vital role as a main component in CFD-based models for simulating wind turbine icing. These models incorporate phase change phenomena to accurately predict the formation, growth, and behavior of ice on turbine surfaces. Understanding phase changes is essential for optimizing turbine design, developing anti-icing strategies, and ensuring the safe and efficient operation of wind turbines in icy conditions.

In wind turbine icing, phase changes primarily involve the freezing of liquid water droplets and the subsequent growth of ice. When liquid water droplets come into contact with cold turbine surfaces or experience a decrease in temperature due to environmental conditions, they undergo a phase change from a liquid to a solid state, known as freezing. Freezing occurs when the droplets’ temperature reaches the freezing point, causing a transformation into solid ice.

One crucial aspect of phase changes in wind turbine icing is supercooling. Supercooling refers to the phenomenon where a liquid is cooled below its freezing point without solidifying. Some liquid water droplets may experience supercooling before freezing due to factors such as a lack of ice nucleation sites or rapid cooling rates. Supercooling delays the phase change from liquid to solid until a critical temperature or external stimulus triggers the nucleation and growth of ice crystals.

Nucleation is the initial formation of ice crystals or nuclei within a supercooled liquid. In wind turbine icing, nucleation can occur on the surfaces of turbine blades or on the droplets themselves. Nucleation serves as a starting point for ice growth, and it can be either heterogeneous or homogeneous. Heterogeneous nucleation happens when ice crystals form on foreign particles or surface irregularities, while homogeneous nucleation occurs spontaneously in a homogeneous liquid.

Following nucleation, ice crystals continue to grow by attracting water molecules from the surrounding liquid or by condensing water vapor from the air. Ice growth is a crucial phase change mechanism in wind turbine icing as it leads to the formation of a solid ice layer. The rate of ice growth depends on various factors, including temperature, humidity, liquid water availability, and the presence of impurities. Understanding these factors and their interactions is essential for accurately modeling ice growth in CFD-based models.

Sublimation is another phase change phenomenon of significance in wind turbine icing. Sublimation occurs when ice directly changes from a solid to a gas without going through the liquid phase. In icing conditions, sublimation can cause ice on the turbine surfaces to evaporate, reducing ice accretion and mitigating the effects of icing. Proper consideration of sublimation in CFD-based models helps assess the loss of ice and its impact on turbine performance.

Incorporating phase changes into CFD-based models involves several considerations and modeling techniques. Numerical simulations and experimental studies utilize phase change models to simulate the freezing process accurately, account for supercooling effects, predict ice growth rates, and estimate ice accretion patterns on turbine surfaces. These models consider factors such as heat transfer, energy exchange, nucleation kinetics, and ice growth mechanisms. Numerous numerical methods for simulating vapor/liquid phase change in CFD have been developed. The commonly utilized methods include the volume-of-fluid (VOF) method [18, 19], level set (LS) method [20], front-tracking method (FT) [21], and the hybrid VOF/LS method [22].

The volume-of-fluid (VOF) method also known as the Enthalpy-Porosity method, is a widely used approach for simulating phase change. It tracks the volume fractions of different phases (e.g., liquid water, ice, and vapor) within a computational domain. The method solves conservation equations for mass, momentum, energy, and species transport, and uses a source term to account for phase change effects. Notwithstanding its slightly lower accuracy in capturing interface curvature, the volume-of-fluid (VOF) method has gained popularity due to its good conservative properties and relatively high computational efficiency [23]. The Level Set Method is another popular approach for tracking phase interfaces and capturing phase change phenomena. It represents the interface between different phases as a level set function. By solving an advection equation for the level set function, the method tracks the evolution of the phase interface and enables the simulation of phase change processes. The Front Tracking Method is a technique that tracks the position and characteristics of sharp phase interfaces, such as the solid-liquid interface during phase change. It involves discretizing the interface into discrete points or markers and solving equations for the motion and deformation of these markers. This method is particularly useful when simulating complex and irregular phase change interfaces.

The various phase change models used in CFD simulations each have their strengths and limitations. The VOF method is known for its good conservative properties. It accurately preserves mass and volume fractions of the different phases during phase change simulations. The FT method generally maintains good conservation properties as it explicitly tracks the phase interface. However, proper implementation is crucial to ensure accurate conservation. The LS method can face challenges in accurately conserving mass and volume fractions due to numerical diffusion and interface smearing, especially for large-scale simulations.

The LS method is generally computationally more expensive than the VOF method. Solving the advection equation for the level set function requires additional computational effort, making it comparatively slower. The FT method can be computationally expensive due to the need to explicitly track individual interface markers. The computational cost increases as the complexity and number of markers increase.

The LS method offers better accuracy in capturing interface curvature compared to the VOF method. It represents the phase interface as a level set function, allowing for more accurate representation of sharp and complex interfaces. The FT method typically provides excellent accuracy in capturing interface curvature. By explicitly tracking the movement and deformation of interface markers, it can accurately represent sharp phase interfaces. The VOF method may exhibit slightly inferior accuracy in capturing interface curvature. This is because it uses a piecewise-constant representation of the phase interface, which can lead to some numerical diffusion and smearing of sharp features.

3. Blade Element Momentum (BEM) theory

Blade Element Momentum (BEM) theory is a mathematical approach used to analyze the aerodynamic performance of rotating blades, particularly in wind turbine applications. It provides a detailed understanding of how the blades interact with the oncoming wind and generate lift and drag forces. The BEM theory provides a comprehensive understanding of the aerodynamic behavior of wind turbine blades, allowing for the optimization of rotor design and performance prediction [12]. While it simplifies the analysis by treating the blade as a series of individual sections, it provides valuable insights into the complex interactions between the blades and the surrounding air flow. The accuracy of BEM theory depends on the assumptions made, such as the airfoil characteristics, wake modeling, and tip loss corrections, which can be refined to improve the accuracy of the predictions.

The BEM theory divides the rotor blade into a number of small elements along its span, typically referred to as blade sections or blade elements. Each element is treated as an isolated airfoil and analyzed individually based on local conditions, such as local wind speed, local angle of attack, and airfoil characteristics. The main concept of BEM theory is to break down the complex aerodynamic interaction between the blade and the flow field into simpler, manageable calculations for each element. By considering the local conditions at each element, such as the local wind speed and the local angle of attack, the theory provides a detailed understanding of how the blade generates lift and drag forces.

Each blade element is typically represented by its own airfoil, which is characterized by specific lift and drag coefficients that depend on the angle of attack. These coefficients are typically obtained from experimental data or airfoil databases. The angle of attack at each element is determined by the relative wind direction, which is the vector sum of the inflow wind speed and the induced velocity caused by the rotor’s rotation and wake effects from neighboring blades. Using the lift and drag coefficients and the local conditions, the theory calculates the lift and drag forces acting on each blade element. The resultant force is then resolved into axial (thrust) and tangential (torque) components. The axial force represents the thrust generated by the rotor, while the tangential force contributes to the torque and rotational motion of the rotor.

To account for the induced effects caused by the rotor’s rotation, BEM theory considers the changes in flow direction and the additional velocity induced by the blades. The induced velocity is typically modeled as a combination of axial and tangential components. It represents the downwash effect caused by the generation of lift. At the blade tips, BEM theory applies Prandtl’s tip loss correction [24, 25]. The correction accounts for the reduced efficiency at the tips due to the high-pressure gradient and vortex shedding. The correction is typically applied by reducing the effective angle of attack at the blade tips compared to the angle of attack at the mid-span sections. This reduction in effective angle of attack is proportional to the strength of the tip vortices and is determined based on empirical or analytical models. By incorporating the tip loss correction, the BEM theory provides more accurate predictions of the aerodynamic forces and power production of the wind turbine. It improves the estimation of thrust and torque at the blade tips, which are crucial parameters for evaluating the performance and efficiency of the rotor. It reduces the calculated forces at the blade tips to improve the accuracy of the analysis. By integrating the calculated forces for each blade element along the rotor span, the theory provides an estimation of the overall thrust and torque generated by the rotor. These parameters are crucial for assessing the aerodynamic performance and power production of the wind turbine.

There are several different tip loss models used in the analysis of wind turbine aerodynamics to account for the reduced efficiency and increased losses at the blade tips. These models aim to capture the impact of tip vortices on the performance of the rotor. The Prandtl’s Classical Tip Loss Model [26, 27] is the original tip loss model proposed by Ludwig Prandtl. It assumes that the tip vortices are concentrated at the blade tips and that the flow behaves as if the blade ends at a virtual point called the “ideal tip.” The model applies an idealized downwash distribution along the spanwise direction to account for the induced flow caused by the tip vortices.

The Glauert’s Tip Loss Model, developed by Friedrich Glauert is an extension of Prandtl’s classical tip loss model. It considers the spanwise distribution of downwash and introduces a correction factor, known as the Glauert correction, to adjust the effective angle of attack at the blade tips. The correction factor is derived based on assumptions about the circulation distribution along the span. It considers that the circulation decreases linearly from the root to the tip of the blades. This assumption implies that the strength of the tip vortices decreases towards the tip, leading to a reduction in the downwash effect on the downstream blades. The Glauert correction factor is defined as the ratio between the local induced velocity at the blade tip and the axial velocity of the inflow wind. It is typically denoted by the symbol “F” and is a function of the local radial position along the blade span. By incorporating the Glauert correction factor, the effective angle of attack at the blade tips is adjusted, which in turn affects the lift and drag forces acting on the blades. This correction helps to improve the accuracy of the analysis by accounting for the downwash effects and the resulting reduction in lift production near the blade tips. It is important to note that the Glauert model assumes an idealized spanwise circulation distribution and does not account for some complex phenomena, such as tip leakage flow or three-dimensional effects. While it provides a significant improvement over the classical model, it still involves simplifications and assumptions to make the analysis tractable.

Goldstein’s Tip Loss Model is a comprehensive approach that accounts for the complex behavior of tip vortices and their impact on the aerodynamic performance of wind turbine blades. It was developed by Alan Goldstein and provides a more accurate estimation of the downwash effect at the blade tips compared to simpler models. Goldstein’s model considers the tip vortices as a rotating vortex sheet along the blade tips, taking into account the variation in the strength of the vortices and the associated downwash distribution. The model utilizes a complex variable method to calculate the induced velocities and adjust the effective angle of attack at the blade tips. In Goldstein’s model, the blade tip is treated as a line source of vorticity, and the flow around the blades is represented by a complex potential function. This complex potential function enables the representation of both the velocity and vorticity fields in a unified manner. By solving the complex potential equation, Goldstein’s model determines the induced velocities caused by the tip vortices at each point along the span of the blades. These induced velocities are then used to adjust the effective angle of attack, taking into account the downwash effect on the downstream blades. The strength of the tip vortices is determined by the circulation distribution along the blades, which can be estimated based on the rotor’s operating conditions and the geometric characteristics of the blades. The model incorporates the variation in circulation along the span to capture the changes in tip vortex strength. The downwash distribution resulting from the tip vortices is determined by the superposition of the individual downwash velocities caused by each vortex. The induced velocities and downwash distribution are then used to adjust the effective angle of attack at the blade tips, influencing the lift and drag forces acting on the blades. Goldstein’s Tip Loss Model provides a more accurate representation of the complex behavior of the tip vortices compared to simpler models. It takes into account the spatial variation of the vortices and provides a more refined estimation of the downwash effect on the blades. This allows for improved predictions of the aerodynamic performance and power production of wind turbines, especially at high tip speed ratios and under highly loaded conditions. It is worth noting that Goldstein’s model is computationally more demanding than simpler tip loss models due to the need to solve the complex potential equation and consider the detailed vorticity distribution. However, with the advancements in computational tools and resources, this model has gained popularity in research and advanced engineering analyses to provide more accurate insights into the behavior of wind turbine blades near the tips.

4. Estimation of power loss in wind turbines under icing conditions

Estimating power loss in wind turbines under icing conditions is a complex task due to the varied and dynamic nature of ice accretion on the blades. Several factors contribute to power loss, including changes in aerodynamic performance, additional weight on the blades, and altered turbine operation.

Icing alters the smooth surface of wind turbine blades, leading to an increase in surface roughness. This roughness causes more air resistance (drag) as the blades rotate through the wind. Consequently, the blades experience reduced lift and a lower ability to convert wind energy into rotational motion, resulting in decreased power generation. As ice accumulates on the turbine blades, it adds significant weight to the rotating components. The extra weight affects the blades’ rotational speed, requiring more force from the wind to maintain the same rotation rate. This increased resistance decreases the turbine’s ability to capture wind energy efficiently and reduces power output. To prevent damage and excessive stress on the blades, wind turbines may adjust their operation during icing events. Blade pitch control is one method used to regulate blade angles and reduce the impact of icing. Additionally, turbines may partially shut down or enter a “de-icing mode” to allow accumulated ice to shed off, thereby protecting the turbine’s structural integrity. These operational adjustments during icing events lead to reduced power generation.

The severity of the icing event plays a crucial role in the power loss experienced by wind turbines. Light icing may have a minimal impact on power output, whereas severe icing can cause substantial reductions. Moreover, different types of ice, such as rime ice or glaze ice, have distinct physical properties that affect aerodynamics differently, leading to varying levels of power loss. Wind turbine blades are designed to optimize aerodynamic efficiency under normal conditions. Some turbines incorporate specialized coatings or materials on the blade surfaces to reduce ice adherence or facilitate ice shedding. Additionally, active anti-icing systems, like blade heating, may be employed to prevent or reduce ice formation. The effectiveness of these design features and anti-icing systems influences the extent of ice accretion and power loss during icing events. The wind speed and direction during icing events play a significant role in the ice distribution on the turbine blades. Icing patterns can be non-uniform, with varying ice thickness along the blade length, depending on the wind characteristics. Turbines may experience higher power losses when encountering turbulent and gusty wind conditions during icing events. Local climate conditions, including ambient temperature and humidity, impact ice formation and adhesion on the turbine blades. Colder temperatures generally lead to more severe icing, while humidity levels affect the water content in the air, influencing the potential for ice accretion.

Estimating power loss accurately requires a comprehensive understanding of these factors, along with access to meteorological data, turbine design specifications, and icing event observations. Developers and investors in wind farms rely on a thorough understanding and accurate assessment of the losses caused by icing. This is essential for estimating the expected energy generation of a project, assessing associated risks, and evaluating its financial feasibility. This section provides an overview of the CFD-BEM and full three-dimensional (3D) CFD methodologies for power loss estimation in wind turbines as a result of icing conditions.

4.1 CFD-BEM approach

The CFD-BEM (Computational Fluid Dynamics - Blade Element Momentum) approach is a computational technique that combines the strengths of CFD and BEM theory to analyze and design wind turbines [28, 29]. This approach provides a comprehensive understanding of the aerodynamic behavior and performance of wind turbine blades, allowing for more accurate predictions and optimizations [30]. While BEM theory is computationally efficient and provides a good estimation of the overall aerodynamic performance of wind turbine rotors, it lacks the ability to capture three-dimensional effects, flow unsteadiness, and complex flow phenomena. CFD, on the other hand, can simulate the detailed flow behavior but is computationally expensive. The CFD-BEM approach bridges this gap by utilizing CFD simulations to generate inputs for the BEM analysis, enabling more accurate predictions of the turbine performance.

In the CFD-BEM approach, CFD simulations are performed to model the flow field around the wind turbine. This can be achieved by incorporating appropriate ice accretion models based on the environmental conditions, such as temperature, humidity, and liquid water content. The CFD component of the approach simulates the airflow and the movement of water droplets in the vicinity of the turbine blades. These simulations help to predict the ice accretion pattern and its effects on the blade surfaces. The CFD solver solves the governing equations of fluid flow, typically the Navier-Stokes equations, in a computational domain that includes the turbine geometry and the surrounding flow domain. The simulations consider the wind inflow conditions, the rotor geometry, and the rotational motion of the blades. The CFD solver discretizes the domain into a grid or mesh and solves the equations iteratively to obtain the flow velocities, pressures, and other flow characteristics while in the BEM analysis, the rotor blades are divided into small sections or elements along the span. Each element is treated as an isolated airfoil and analyzed based on local conditions such as the local wind speed, angle of attack, and airfoil characteristics.

CFD is used to discretize the entire fluid domain into a grid of cells, where the governing equations (typically the Navier-Stokes equations) are solved for each cell. On the other hand, BEM is employed to discretize the boundaries or surfaces of objects in the domain, where the integral form of the governing equations is solved. In this approach, the discretization scheme involves volume discretization where the fluid domain is divided into a finite number of cells, usually using techniques like structured or unstructured grids. Structured grids have a regular arrangement of cells, such as a Cartesian or curvilinear grid. Unstructured grids have cells with arbitrary shapes, which can be more flexible for complex geometries. The boundaries or surfaces of objects in the domain are discretized into a collection of elements. These elements can be triangles, quadrilaterals, or other suitable shapes depending on the geometry and requirements of the problem. The choice of element shape affects the accuracy and efficiency of the BEM calculations. Once the domain is discretized, numerical integration methods are employed to approximate the integral terms in the governing equations. These integrals arise due to the BEM formulation, where the influence of boundary values on the interior points is determined. Various integration schemes like Gaussian quadrature or midpoint rule can be used depending on the accuracy requirements and complexity of the problem. The boundary conditions, such as velocities or pressures, are prescribed or obtained from the problem statement. These conditions are then applied to the discretized boundary elements, ensuring that the appropriate values are used during the calculations. The discretized equations, including the CFD equations for the interior cells and the BEM equations for the boundary elements, are solved iteratively or simultaneously. The solution procedure involves solving the equations for each cell or element while considering the interdependence between neighboring cells or elements. Once the solution is obtained, post-processing is performed to analyze and visualize the results. This may involve calculating derived quantities of interest, such as forces, velocities, or pressure distributions, and presenting them in a suitable format for analysis or presentation.

In the context of the CFD-BEM approach, the iterative method refers to the iterative process used to solve the discretized equations of the combined CFD and BEM formulations [31]. Since the CFD and BEM equations are coupled through the boundary conditions, an iterative approach is employed to obtain a converged solution. The process starts with an initial guess for the solution variables. These variables can include the velocity components, pressure, and other relevant quantities. The CFD equations for the interior cells are solved using the current values of the solution variables. This typically involves discretizing the governing equations, applying boundary conditions, and solving the resulting system of equations using numerical methods like finite difference, finite volume, or finite element techniques. With the updated values of the solution variables from the CFD step, the BEM equations for the boundary elements are solved. This involves evaluating the influence coefficients between boundary elements and determining the boundary values at each element. The boundary values are typically obtained by integrating the influence of neighboring elements and interior cell values using numerical integration techniques. The BEM-calculated boundary values are then used to update the boundary conditions for the CFD calculations. These updated boundary conditions are applied in the subsequent CFD step to solve the equations for the interior cells. This coupling ensures consistency between the CFD and BEM solutions. The BEM-calculated boundary values are then used to update the boundary conditions for the CFD calculations. These updated boundary conditions are applied in the subsequent CFD step to solve the equations for the interior cells. This coupling ensures consistency between the CFD and BEM solutions. Once the iterative process has converged and a satisfactory solution is obtained, post-processing is performed to analyze and visualize the results. This may include calculating derived quantities, generating plots or animations, and extracting relevant information for further analysis or interpretation.

The BEM analysis uses the inputs from the CFD simulations to estimate the lift, drag, and other aerodynamic forces on each blade element. These forces are integrated along the blade span to obtain the overall aerodynamic performance of the rotor, including thrust, power, and torque. The CFD simulations provide the necessary inputs for the BEM analysis, including the local wind flow conditions and the aerodynamic characteristics of the airfoils. The CFD results, such as the flow velocities and pressures, are interpolated and integrated into the BEM analysis to obtain more accurate estimates of the aerodynamic forces and power output of the wind turbine. This integration is typically performed at specific points along the blade span, where the BEM analysis is carried out.

Using the aerodynamic forces obtained from the BEM analysis, the power loss in the wind turbine can be estimated. The power output of a wind turbine is directly related to the aerodynamic forces acting on the blades. By comparing the power output of the iced turbine with the power output of the clean turbine (without ice accretion), the power loss can be quantified. The power loss is typically expressed as a percentage of the clean turbine power. In addition to the overall power loss estimation, the CFD-BEM approach can also provide insights into load-based power losses. Ice accretion affects the distribution of aerodynamic forces along the blade span, leading to uneven loading on different blade sections. This uneven loading can result in structural and fatigue issues. By analyzing the changes in load distribution caused by ice accretion, the CFD-BEM approach can estimate load-based power losses and identify critical areas prone to structural concerns. By considering load-based power losses, the CFD-BEM approach provides a comprehensive analysis of the impact of icing on wind turbine performance beyond the overall power loss estimation.

Power loss estimation also considers the dynamic effects induced by ice accretion. Ice adds mass to the blades, altering their natural frequencies and mode shapes. This can lead to resonant vibrations and increased fatigue damage. The CFD-BEM approach accounts for these dynamic effects by incorporating the changes in blade properties due to ice accretion, allowing for a more accurate estimation of power losses associated with blade dynamics. Ice accretion on wind turbine blades is a time-dependent process, with ice formation, growth, and shedding occurring over time. The CFD-BEM approach can perform time-dependent analyses to capture the temporal evolution of ice accretion and its impact on power losses. By simulating the icing process and updating the BEM analysis at different time steps, the approach provides a more realistic estimation of the dynamic power losses associated with changing ice conditions.

To account for uncertainties in the power loss estimation, the CFD-BEM approach can incorporate uncertainty quantification techniques. Uncertainties may arise from variations in meteorological conditions, ice accretion models, and input parameters. By performing sensitivity analyses and probabilistic simulations, the approach can provide confidence intervals or probability distributions of the estimated power losses, allowing for a better understanding of the associated uncertainties.

4.2 Full three-dimensional (3D) CFD approach

3D CFD refers to the numerical simulation of fluid flow and heat transfer in three-dimensional geometries. It is a powerful tool used to analyze and predict the behavior of fluids under various conditions, including the estimation of power loss in wind turbines during icing conditions. When wind turbines operate in cold climates, they can be exposed to icing conditions where ice accretes on the rotor blades. The presence of ice on the blades can significantly impact their aerodynamic performance, leading to reduced power output and increased structural loads. To estimate the power loss caused by icing, 3D CFD simulations can be employed. This process typically involves geometry modeling, mesh generation, applying boundary conditions, solver and simulation and post-processing.

In geometry modeling, the wind turbine’s components, including the rotor, nacelle, and blade surfaces, are accurately represented using computer-aided design (CAD) software. The rotor blades are typically modeled with high geometric fidelity, including airfoil shapes and surface roughness. The geometry model serves as the basis for constructing the computational domain for the CFD simulation. To model the wind turbine’s rotor, nacelle, and other relevant components, detailed CAD models are created. The rotor blades are of particular importance and are modeled with high geometric fidelity. This includes accurately capturing the airfoil shapes, which determine the aerodynamic performance of the blades. The CAD model should also account for other blade characteristics such as twist, chord length, and surface roughness. Surface roughness plays a crucial role in determining the airflow behavior around the blades. It affects the boundary layer development, separation points, and ultimately the aerodynamic performance. Therefore, the CAD model must incorporate realistic surface roughness profiles based on the blade’s manufacturing specifications. In addition to the blades, the nacelle, hub, and other wind turbine components are also modeled in detail. The nacelle includes the generator, gearbox, and control systems. The hub connects the rotor blades to the main shaft. Accurate representation of these components ensures that the simulation captures the overall behavior of the wind turbine system. The CAD model is typically created by experienced designers using specialized software. It requires expertise in wind turbine design and knowledge of the specific turbine being analyzed. The fidelity and accuracy of the geometry model have a direct impact on the reliability and accuracy of the subsequent CFD simulations. Once the geometry model is created, it needs to be converted into a suitable format for CFD simulations. This may involve meshing the geometry, which is the process of dividing it into small elements or cells for numerical analysis. The meshing process, discussed in the subsequent steps, is essential to accurately capture the flow behavior and aerodynamic characteristics of the wind turbine components. Overall, the geometry modeling step establishes the foundation for the subsequent stages of the 3D CFD process. It ensures that the wind turbine components are accurately represented in the computational domain, allowing for realistic simulations and analysis of the fluid flow and heat transfer during icing conditions.

The computational domain is divided into small control volumes or cells through mesh generation. The mesh captures the geometric details of the wind turbine components and discretizes the computational domain. A fine mesh is crucial to capture the flow features accurately, especially near the blade surfaces. Advanced meshing techniques, such as boundary layer resolution and adaptive mesh refinement, may be employed to enhance the simulation’s accuracy. There are various techniques for mesh generation, including structured, unstructured, and hybrid methods. In structured meshes, the cells are arranged in a regular, well-organized pattern, such as a Cartesian or cylindrical grid. Structured meshes offer good numerical accuracy and efficiency but may struggle to handle complex geometries. Unstructured meshes, on the other hand, allow for more flexibility in representing complex geometries. The cells in unstructured meshes can be of different sizes and shapes, adapting to the geometry of the wind turbine. They are typically generated using techniques like Delaunay triangulation or advancing front methods. Unstructured meshes are well-suited for capturing the flow features around complex geometries but can be computationally more expensive. Hybrid meshes combine structured and unstructured elements to leverage their respective strengths. For example, structured meshes can be used in regions where the geometry is relatively simple, while unstructured meshes can be employed in complex areas. Hybrid meshes provide a balance between accuracy and computational efficiency. Special attention should be given to the near-wall region, particularly around the blade surfaces, where the flow behavior is highly influenced by boundary layer effects. Fine mesh resolution is required in these regions to accurately capture the flow gradients and resolve the boundary layer structures. To improve the accuracy of the simulation, adaptive mesh refinement techniques can be employed. These techniques dynamically refine the mesh in areas with high flow gradients or areas of particular interest, such as regions experiencing ice accretion. By refining the mesh in these critical areas, the simulation can capture the local flow features more accurately without requiring an excessively fine mesh throughout the entire domain. Overall, mesh generation in 3D CFD involves selecting an appropriate mesh type (structured, unstructured, or hybrid) and generating a mesh with suitable resolution, especially in areas of interest. The goal is to strike a balance between accuracy and computational efficiency to ensure an effective simulation of the fluid flow and heat transfer in the wind turbine system during icing conditions.

The boundary conditions define the inflow conditions, environmental parameters, and ice accretion rates. These conditions are based on the anticipated environmental conditions during icing events. The wind inflow conditions specify the velocity, temperature, and humidity of the incoming air. Additional parameters such as ice water content and droplet size distribution are also considered to accurately model ice accretion on the blades.

The CFD solver numerically solves the governing equations, such as the Navier-Stokes equations, in the computational domain. The solver simulates the fluid flow, icing, and heat transfer processes, considering the time-dependent behavior of the icing process. The simulation progresses in discrete time steps, allowing the computation of the airflow, ice accretion, and thermal behavior of the wind turbine blades. The solver utilizes numerical algorithms and iterative methods to approximate the solutions of the governing equations.

Once the simulation is completed, post-processing techniques are employed to analyze and interpret the results. Visualization tools enable the examination of flow patterns, temperature profiles, and ice accretion distribution on the blade surfaces. Quantitative analysis provides insights into the power loss caused by icing. Power coefficients, which represent the efficiency of power conversion from the wind, are calculated and compared for different icing scenarios. These analyses help identify the regions of maximum power loss and provide valuable information for design optimization, such as anti-icing systems, blade de-icing strategies, and operational adjustments during icing events.

5. Future research directions and conclusion

Future research in power loss estimation in wind turbines due to icing can explore several directions to improve our understanding and develop more accurate models for predicting and mitigating the effects of icing. Enhancing the accuracy of icing models is crucial for more precise power loss estimation. Future research can focus on developing advanced CFD models that capture the complex interaction between ice accretion, airflow, and turbine performance. Incorporating more detailed physics, such as ice formation and shedding dynamics, can lead to more accurate predictions of power loss. Current icing models primarily focus on the aerodynamic aspects of ice accretion, but future research can expand the scope to include multi-physics modeling. This involves incorporating thermal, mechanical, and electrical aspects of ice formation and shedding. By considering the thermodynamic behavior of ice formation, mechanical stresses induced by ice loading, and the electrical effects of ice on blade performance, more comprehensive models can be developed.

To improve the accuracy of icing models, it is essential to validate them using extensive field measurements. Collaborations with wind farm operators and instrumented turbine installations can provide valuable datasets on ice accretion and turbine performance under various icing conditions. These field measurements can serve as benchmark data for validating and refining icing models, enhancing their predictive capabilities. Light Detection and Ranging (LiDAR) technology has shown promise in measuring ice thickness and shape remotely. By emitting laser pulses and analyzing the reflected signals, LiDAR systems can provide detailed information about ice accumulation on wind turbine blades. Ongoing research can focus on improving the resolution and accuracy of LiDAR-based ice measurement techniques, enabling real-time monitoring of ice conditions across large wind farms. Also, combining data from multiple sensors and measurement techniques can enhance the accuracy and reliability of ice accretion estimation. By fusing data from LiDAR, thermal imaging, and optical methods, researchers can obtain a more comprehensive understanding of ice accumulation on wind turbine components. Data fusion techniques, along with advanced signal processing and data analysis algorithms, can be explored to improve the accuracy and robustness of ice measurement systems. Advancements in ice accretion measurement techniques will not only improve power loss estimation in wind turbines but also aid in the development of effective anti-icing and de-icing strategies.

While wind tunnel tests and scaled-down turbine experiments are useful for studying ice accretion, the scaling effects between laboratory conditions and full-scale wind turbines can introduce discrepancies. Future research can investigate the scaling laws and conduct full-scale testing to validate the models under real-world conditions. Full-scale testing provides a comprehensive understanding of the impact of icing on large wind turbines and allows for the validation of model predictions at the actual operational scale. Comparative studies involving multiple wind turbine models and designs can provide insights into the influence of turbine characteristics on icing-related power loss. By testing different blade geometries, materials, and surface treatments, researchers can evaluate their performance under icing conditions. These studies help identify design features that are more resistant to ice accretion and shedding, leading to reduced power losses.

Research can focus on optimizing anti-icing and de-icing systems to improve their effectiveness and energy efficiency. This includes developing advanced coatings, heating technologies, and aerodynamic designs that minimize ice accretion or facilitate ice shedding. Exploring novel materials and techniques for preventing ice adhesion and developing adaptive de-icing systems can lead to more reliable and cost-effective solutions.

Icing has a significant impact on the performance and power loss in wind turbines. When ice accumulates on various parts of the turbine, it alters the aerodynamic characteristics, adds weight, and affects the overall operation of the turbine. Ice formation on the turbine blades disrupts the smooth flow of air and modifies the aerodynamic profile of the blades. The ice adds thickness and roughness, changing the shape of the airfoil and reducing its lift-to-drag ratio. As a result, the blades experience increased drag and reduced lift, leading to reduced turbine efficiency and power output. Icing can also cause increased turbulence in the airflow around the blades.

The irregular surface created by ice formations generates local flow separations, vortices, and eddies. This turbulent flow disrupts the ideal aerodynamic conditions and further reduces the efficiency of the blades. Turbulence also increases the loads on the turbine components, potentially leading to structural damage. Ice accumulation adds weight to the rotating turbine blades. The additional mass creates an imbalance, which can cause uneven loading on the rotor and generate vibrations. These imbalances and vibrations can lead to increased mechanical stresses and damage to the turbine structure, bearings, and other components. The presence of ice on the blades can cause erosion as the ice particles impact the leading edges. The repeated impact of ice particles can remove protective coatings and damage the blade surface, further degrading the aerodynamic performance and reducing power output.

Understanding and accurately estimating the power loss caused by icing is essential for optimizing wind farm operations, improving energy production efficiency, and ensuring the reliability of wind energy generation. In summary, power loss estimation in wind turbines due to icing is a complex and multi-faceted challenge. Ongoing research efforts, technological advancements, and collaborations between academia, industry, and policymakers are essential for improving the understanding of icing effects, refining prediction models, and implementing effective mitigation strategies. By addressing the power loss associated with icing, it is possible to enhance the reliability, efficiency, and economic viability of wind energy generation in cold climate regions, contributing to the sustainable and clean energy transition.

In conclusion, power loss estimation in wind turbines due to icing is a crucial aspect of wind farm operation in cold climate regions. Icing on wind turbine components, particularly the blades, can significantly impact performance and reduce power output. Overall, icing on wind turbines can cause a significant reduction in power output. The extent of power loss depends on the severity and duration of icing events, the turbine’s design and anti-icing measures, and the specific aerodynamic and structural characteristics of the turbine. Proper monitoring, icing detection, and mitigation strategies are crucial to minimizing power losses and ensuring safe turbine operation in icing conditions.