Wind turbine ice detection using AEP loss method: A case study

This paper describes the comparison of a statistical and numerical case study of wind resource assessment and estimation of resultant Annual Energy Production due to ice of a wind park in ice prone cold region. Three years Supervisory Control and Data Acquisition data from a wind park located in arctic region is used for this study. Statistical analysis shows that the relative power loss due to icing related stops is the main issue for this wind park. While Larsen wake model is used for the CFD simulations, where results show that it is important to use the wake loss model for CFD simulations of wind resource assessment and AEP estimation of a wind park. A preliminary case study about wind park layout optimization has also been carried out which shows that AEP can be improved by optimizing the wind park layout and CFD simulations can be a good tool.


Introduction
Due to increasing demand of electrical power and needs of protecting the environment, there has been a growing need of rapid expansion and better use of renewable energy resources to reduce the carbon emissions. Cold regions have good wind resources, but environmental challenges such as icing has been recognized as hindrance in proper utilization of these good wind resources. Worldwide, installed wind energy capacity in ice prone cold regions in 2015 was 86.5 GW, which is expected to reach 123 GW in year 2020 (Davis, 2014). Therefore, wind resource assessment in cold regions is important both for the operation of wind parks and also to provide more accurate wind energy production forecasts (Bourgeois and Beckford, 2018). The assessment of wind resources at cold climate sites is challenging and include aspects as terrain topology of the area where the wind park is to be located, wind behavior, type of buildings or any other human intervention around the wind park location that may induce changes in the terrain roughness as well as the geometry of the wind turbine and its expected performance in the calculation of the tested scenarios. Low temperatures and icing climate set additional challenges for wind resource assessment in cold regions.
The International Energy Agency (IEA) T19 set out a statistical method to assess production losses due to icing based on standard SCADA data available from modern wind parks, named ''T19IceLossMethod'' and developed by international expert group IEA Wind TCP Task 19 (Bourgeois and Beckford, 2018). T19IceLossMethod is the standardized open source Python code model (IEA Wind Task 19, 2019) for assessment of wind power production loss (Lehtomaki et al., 2018) and compared different sites with a systematic analysis method, as well as evaluate the effectiveness of various blade heating systems (Lehtomaki et al., 2019). The scope of work includes a detailed historical SCADA data analysis using 10-minute averages from four different wind turbine models and analysis of icing conditions on sites (T19IceLossMethod, 2019).
CFD based wind resource assessment can provides improved agreement with the field measurements, when compared with the statistical flow models. Davis et al. (Davis et al., 2016) performed CFD analysis for wind resource assessment and wind park design layout optimization, where they found that design of a wind park layout based upon CFD simulations can help to increase the wind energy production. Wind turbine wake effects are important for the wind resource assessment and can be responsible for approximately 10%-20% losses of the AEP (Bilal, 2016). Wind turbine wake effects on neighboring turbines effects the flow behavior, which also effects the wind turbine performance and power production. The development of models to describe wind turbine wake effects on energy production, include different initiatives of simplified or linearized models such as Jensen model (Katic et al., 1986) or more complex such as the eddy viscosity together with Gaussian shape (Anderson, 2009;Renke, 2007) or second order model combined with the Prandtl's assumption of velocity been function of the distance power 21/3 (Larsen, 1988) or considering turbulence on the wake recovery as function of the distance from the wind turbine (Ishihara et al., 2004). Wake structure creates a flow disturbance behind the wind turbine as well as in front of the turbine, which differ from the flow around an idealized turbine (Gonza´lez-Longatt et al., 2012) as shown in Figure 1.
Wake losses and aerodynamic behavior of wind turbine can be better addressed once the flow over the wind park site is better understood. Wind flow over steep slopes and ridges causes flow separation which is hard to capture and simulate using mesoscale models such as WRF. Furthermore, the atmospheric stability influences the wind flow (Bilal et al., 2015). This paper describes the statistical and numerical case study of wind resource assessment and wind park design layout effects on energy production in icing climate. T19IcelossMthod based statistical analysis and CFD based numerical simulations are carried out in comparison with the 3 years (2013-2015) wind park SCADA data. To better understand the wind turbine wake effects on flow behavior and resultant power production, Larsen wake model is used to calculate the AEP of each wind turbine.

Wind park terrain and layout
Wind park used for this study consists of 14 wind turbines located at 400 m a.s.l in the arctic region of Norway. Three years (2013-2015) Supervisory Control and Data Acquisition (SCADA) from all wind turbines is collected and sorted for this study. This data mainly includes wind speed, wind direction, temperature and power production. Wind park terrain is hilly and is mainly covered with grass in summer and snow in winter. Figure 2 shows the wind park layout, terrain, and wind rose map.
Wind direction is divided in 12 sectors, where first sector is oriented directly north. Wind climatology is presented as a wind rose, giving the average wind speed distribution divided in the velocity intervals (bins) and wind directions (sectors). The frequency distribution has been fitted to a Weibull distribution. Results of wind climatology for all 12 sectors are given in Table 1.

Statistical model
SCADA based T19IceLossMethod is the statistical model developed to minimize the required inputs and to create a solution that would not require any external measurements, such as ice detectors. Iced wind turbine power losses are defined by comparing the performance to the calculated power curve using heated anemometers from nacelle and the measured reference, expected power curve. In general, the method can be divided into three main steps (Bourgeois and Beckford, 2018) and site air density and air pressure are used to calibrate the nacelle wind speed as equation (1): where v site is calibrated nacelle wind speed, v std measured nacelle wind speed, T site is nacelle temperature and T std is standard temperature of 15°C, which resulting to air density of 1.225 kg/m 3 at sea level P std = 101325 Pa ambient air pressure, h is site elevation in meters above sea level (Engineering ToolBox, 2003). In order to detect production loss, a baseline needs to be established. Icing event is defined as a period of time when the output power of wind turbine is below a set limit in cold temperature. These limits are calculated form a subset of the data where temperature is high enough so that it is reasonable to assume that there is very little icing. The detailed setup is described in Table 2. The statistical model uses the 10th percentile of the power curve to trigger an icing alarm. The percentile method was chosen mainly due to robustness and simplicity of the definition (Engineering ToolBox, 2003) This method  classifies icing events into three main categories: Class a means the reduced production due to icing, Class b means the stop caused by icing and Class c is apparent over production (T19IceLossMethod, 2019).
Most common icing effect on turbine operation is production loss, and the production losses are caused by degraded aerodynamics of the blade caused by ice accretion along the blade surface. If turbine stops during an icing event, it is assumed that the icing event is considered to be a contributing factor to the stop, where an icing induced stop is assumed to be ended once turbine output power is above the icing alarm limit again. In some cases turbine seems to be producing more than expected, this can happen if for example an anemometer is experiencing icing. An iced cup anemometer will show lower wind speeds than what the actual wind speeds are. This in turn will result in apparent overproduction in SCADA.

Numerical model
CFD based numerical simulations are carried out using WindSIM software, which is a modern Wind Park Design Tool (WPDT) that helps to optimize the wind park energy production by using non-linear mathematical methods. 3D terrain model of the wind park is generated where the domain is discretized using hexahedral numerical grid. Iterative numerical simulations of air flow behavior are performed by solving Reynolds Averaged Navier-Stokes equations (RANS). The energy equation is neglected during this study, as temperature is assumed constant in the region close to the ground surface. So, the exchange of heat and water vapor at the Earth's surface is neglected. In order to account for surface roughness in the numerical simulations, Wieringa's classification is used. RNG-based k-e turbulence model is used due to its better agreement with the flow profiles and length of the separated flow region. The detailed setup is shown in Table 3.
The default value of terrain surface roughness heights is read from the grid file in WindSIM, or, alternatively a constant roughness height can also be imposed in the model by specifying a non-zero value for the roughness height. The roughness height is used in the velocity profile log-law, given in equation (2): where u is the wind velocity; u Ã is the friction velocity, k is von Ka´rma´n constant (k = 0.435); z is the coordinate in vertical direction; z 0 is the roughness height. Wind park digital terrain layout used for the CFD simulations is shown in Figure 3.
In this study, the wind park terrain corresponds to mesoscale topology of the region around the wind park. The areas selected for the wind flow simulations is the section that involves different sizes in the direction of the flow to determine the influence of the natural formation of the region on the wind profile. In the wind park terrain, the presence of a natural channel along the mountain line west-east oriented at the north side and water body at the south represent a challenge for the selection of the domain to do the numerical simulations. Estimation of latitudinal and longitudinal extension of the domain is quite important because the wind behavior will be directly affected by the surface shape, following the mass and momentum conservation equations. For a fixed latitude range that cover the top of the mountain and a more flatter area to the south, three different territorial domains based on the variation of the starting point of the domain in the longitudinal direction are tested in this numerical study. Mesh sensitivity analysis are carried out, where a final refined mesh by default placed in the center of the domain and the cell distribution is uniform within the refined area with an increasing cell size toward the borders with cell resolution of 28m is selected for the numerical study. Figure 4 shows the numerical mesh used for this study, the grid spacing, and cells is described in Table 4. Larsen et al. wake loss model (Anderson, 2009) is used for this study. This wake loss model is based on calculating the normalized velocity deficit and rotational axisymmetric along the x-axis, dV, as described in equation (3):   where, U is the free stream velocity, and V is air velocity at some point after the turbine rotor. This model is derived from the turbulent boundary layer equations and a similarity assumption. Shows as equation (4): where, C T is the thrust coefficient, and the constant A r and variable C 1 based on rotor diameter (m) D, are described in equations (5) and (6): where, x 0 is described in equation (7), and referred to R 95 which described by using R nb that correlated with the ambient turbulence intensity at hub height, I a , as shown in equations (8) and (9):

Ice detection
Result from T19IceLossMethod are compared with the field SCADA data from the wind park. All three categories of icing events: Class a, Class b, and Class c are used as relative production losses due to icing (% of kWh), relative losses due to icing related stops (% of kWh) and over production hours (% of total hours) respectively. The results from this statistical model are shown in Table 5, and the calculation of losses time of all icing events is illustrated in Table 6. The maximum values are shown with pink color and the minimum values are shown with green color. According to the statistical method, the bold and underlined data zones are where the maximum power losses in each icing event categories, which means annual power production is not ideal in these turbines, and need to optimize the wind turbine sites or operation methods to improve the actual wind power productions. Whereas the bold and italicized data zones are the minimum power losses where the actual wind energy production is ideal in these turbines. To consider the operating hours of each turbines during year 2013 to 2015, turbine 01 are the maximum loss of total hours in 2013 and 2015. However, turbine 01 in 2014 has the minimum and turbine 13 has the maximum loss of total operating hours. Results of maximum ice loss detection on wind power production for turbine 01 and turbine 13 are shown in Figure 5 by using 10th (P10) and 90th (P90) percentile of power in the bin compared with standard power production with median output power in the bin. Analysis shows that the maximum reduced power production due to icing with T01 is 2.1% in 2013, 0.3% in 2014, and T08 is 0.5% in 2015, whereas the maximum ice loss with class a for all 14 wind turbines. Analysis about the maximum stop caused by icing with T02 is 9.1% in 2013, T13 is 29.6% in 2014, and T01 is 13.4% in 2015 respectively, where the maximum ice loss with class b for all 14 wind turbines. For the ice loss with class c for all 14 turbines, analysis shows that T04 is 3.0% in 2013, T09 is 4.2% in 2014, and T04 is 4.2% in 2015 with the maximum apparent overproduction. In total, added up ice losses in different categories of all 14 wind turbines, class b is the maximum ice detection category for all 3 years of the wind park, with 57.1%, 54.6%, and 55.5% respectively, whereas class a is the minimum ice detection category during 2013 to 2015 with ice loss 4.7%, 2.8%, and  3.0% separately. In other words, the relative power losses due to icing related stops are the main issue for this wind park.

Wind resource assessment
Results from CFD simulations are compared with the field SCADA data from the wind park. Both simulated and measured results of average wind conditions are used for AEP estimation. For the CFD simulations, three hub heights (35, 80, and 125 m) which equal the lower and upper tip height and the hub height are used for the comparison. Wind resource maps for all three hub heights are shown in Figure 6, where simulation results show that the wind park locates and operates in rich wind resources area. The results from wind resource map also shows that wind speed at the hub height of 125 m is the highest and 35 m is the lowest in all 3 years. The main reason for it is the planetary boundary layer effects, that is, the friction near the surface, due to no-slip condition, which results in some mass flux of the air being displaced into layers at higher altitudes, which by mass-conservation principle will result in adjacent layers slowing down. In addition, according to the pioneer research work (Liu et al., 2010) and the topographic characteristics of actual wind park with specific surface roughness height study in this paper, when the airflow through the mountain and lake, the flow will form a recirculation zone. If the wind speed is high or the atmosphere structure is unstable, a strong vortex motion will appear on the leeward side of the mountain. This vortex motion may cause the wind speed to rotate at a certain height, thus showing the phenomenon that wind speed decreases with height. This affects the velocity deficit in near-surface layers that result in change of wake and associated power losses.

AEP estimation and comparison
Numerical analysis are carried out based on the geographical long-term conditions in order to calculate the AEP of the wind park using Larsen wake loss model with actual surface roughness height. Each year SCADA data is further analyzed and compared by T19IceLossMethod statistical model and CFD simulation. This approach is used to get the better estimate of energy based availability within field SCADA data, statistical model and numerical simulations of the wind park. In 2013 and 2015, both statistical and CFD results showed overestimated AEP values when compared with the SCADA data, whereas the underestimated AEP values show mostly in 2014 for CFD simulations. In other words, the wind park have better estimation of energy power production compare in year 2013 and 2015. Then further CFD analysis are carried out by taking into account the wind turbine wake losses. Results of this parametric study are shown in Table 7 and Figure 7. Analysis show that CFD results without use of wake model overestimated the AEP in comparison to SCADA data, whereas CFD results with the use of Larsen wake loss model is found in better agreement with the SCADA data in 2013 and 2015, whereas for 2014 are in reverse. This case study highlights that it is important to use the proper wind turbine wake model for CFD analysis while performing the wind resource assessment.

Wind park layout optimization
To further improve the AEP and reduce wake loss effects of the wind park, a preliminary numerical case study is carried out to optimize the existing wind park layout (wind turbine locations). WindSIM Park Optimizer is used for this purpose. Park optimizer is a numerical tool coupled with the WindSIM that helps to maximize the wind farm profitability by optimizing the wind park layout. Park optimizer uses two type of optimization approaches, (1) basic optimization, (2) WFD (wind farm design) optimization. The basic optimization is based on heuristic algorithm which gives near optimal results. Algorithms such as genetic algorithm, simulated annealing and similar trial-and error algorithms for the wind park layout design also provide good results but they are slow in process and do not provide information about the quality. The WFD optimizer is highly innovation approach based on formal research operation methods and is derived from state of art optimization solver. For this case study, we have used WFD optimizer. Results from CFD simulations of year 2014 presented in previous section are used as input to the WFD optimizer, where WFD optimizer optimally relocated the wind turbines with respect to the highest average wind speed and wind energy. Figure 8 shows the optimized layout of the wind park in comparison to existing layout.
To verify that how much AEP is increased after optimization of the wind park layout, CFD simulations for optimized wind park layout are carried out. Larsen wake model is used for this study. Table 8 shows the AEP and total wake loss results of the Wind Park. Before optimization of the wind park layout design, gross AEP from CFD analysis has been calculated as 95.61 GWh/y, whereas after optimizing the wind park layout the gross AEP calculations from CFD simulations shows a value of 112 GWh/y. Comparison of AEP for optimized design layout   of the wind park shows that 16.39 GWh/y energy increased as compared to the existing layout design. Moreover, the CFD based on Larsen wake model have decreased 7% in total wake loss compare before and after optimization used WFD optimizer. This shows that CFD based numerical optimization approach can be a good way to better simulate the wind resources and optimize the wind park layout and the power production.

Conclusion
SCADA data analysis, T19IceLossMehtod and CFD simulations are carried out for a wind park located in arctic region. AEP from SCADA data is compared with statistical and CFD simulations, where a good agreement is found. Statistical analysis shows the relative power loss due to icing related stops (Class b) is the main issue for this wind park from 2013 to 2015. While the CFD analysis shows that it is important to use the proper wake loss model for better estimation of the wind resources using numerical simulations. Analysis also show that CFD results without use of wake loss model can over/under-estimate the AEP in comparison to SCADA data, whereas CFD results with the use of wake loss models are in good agreement with the SCADA data. Analysis shows the phenomenon that wind speed decreases with height increase due to the wind park have a strong vortex motion appear on leeward side of the mountain and affects the velocity deficit in near-surface layers, which result in change of wake and associated power losses. Moreover, study about wind park layout optimization shows that AEP can be optimized by optimizing the wind park layout and CFD simulations can be used as a tool in this regards.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The work reported in this paper is supported by the University of Tromsø PhD project no (381100/74104) and the publication charges for this article have been funded by a grant from the publication fund of UiT-The Arctic University of Norway.

Code and data availability
T19IceLossMethod is the standardized open source Python code model open-source simulation tools developed by international expert group IEA Wind TCP Task 19. The GitHub repository https://github.com/IEAWind-Task19/T19IceLossMethod. The SCADA data and other farm-specific information for the commercial wind farm discussed in Section ''Wind park terrain and layout'' are not publicly available due to proprietary privacy restrictions.