Identification of ice mass accumulated on wind turbine blades using its natural frequencies

Abstract

This work demonstrates a technique to identify information about the ice mass accumulation on wind turbine blades using its natural frequencies, and these frequencies reduce differently depending on the spatial distribution of ice mass along the blade length. An explicit relation to the natural frequencies of a 1-kW wind turbine blade is defined in terms of the location and quantity of ice mass using experimental modal analyses. An artificial neural network model is trained with a data set (natural frequencies and ice masses) generated using that explicit relation. After training, this artificial neural network model is given an input of natural frequencies of the iced blade (identified from experimental modal analysis) corresponding to 18 test cases, and it identified ice masses’ location and quantity with a weighted average percentage error value of 17.53%. The proposed technique is also demonstrated on the NREL 5-MW wind turbine blade data.

Keywords

Wind turbine blade ice detection natural frequency experimental modal analysis artificial neural network

Introduction

Wind turbines are increasingly installed in the northeastern and the mid-Atlantic United States, Canada, and Northern Europe due to good wind resources and land availability. In these regions, humidity along with low temperatures in the winter increases the risk of ice accumulation on wind turbine components. The global wind energy installations in cold climate regions reached a capacity of 127 GW at the end of 2015, and the forecast is that it would reach a capacity of 186 GW by the end of 2020 (Emerging from the Cold, 2016). Icing of the rotor blades results in reduced turbine power output as it reduces the lift force and increases the drag force acting on blade’s airfoil sections (Rindeskär, 2010; Turkia et al., 2013). Additionally, the nacelle vibration amplitudes increase during the icing conditions (Skrimpas et al., 2015), so icing can be more reliably detected using the power curve analysis along with nacelle vibration analysis. Etemaddar et al. (2014) investigated the effects of atmospheric ice accumulation on the aerodynamic performance and structural response of the wind turbines, and they predicted that the relative change in mean value is higher than the change in standard deviation for most of the response quantities (rotor speed, torque, power, thrust, and structural loads) of the iced blade. Gantasala et al. (2016) analyzed the modal behavior of a 2-MW wind turbine blade under icing conditions. Their analysis showed that the aeroelastic damping factors of the vibration modes reduce with ice accumulation on the blade. Rissanen et al. (2016) proposed simulation parameters for predicting the dynamic behavior of iced wind turbines to generate icing design load cases for the new IEC 61400-1 ed4.

Ice accumulation on the wind turbine blades is not uniform and its mass distribution changes under different stages of icing and ice shedding during operation. The natural frequencies of a blade reduce according to the quantity and location of ice mass along its length. Lorenzo et al. (2013) investigated the influence of ice mass on the natural frequencies of the NREL 5-MW turbine by extracting model parameters from the operational modal analysis. Alsabagh et al. (2015) considered different ice mass distributions as defined in the ISO 12494:2001 (2001) standard on a multi-megawatt wind turbine blade. They analyzed the influence of ice mass on natural frequencies and dynamic magnification factors (ratio of the dynamic deflection to static deflection), by considering only the mass changes in the blade. Wind turbines are usually exposed to turbulent wind conditions which will excite natural frequencies of the blade, and these frequencies can be extracted from the blade’s vibration measurements. The ice detection systems like BLADEcontrol (Brenner, 2016), fos4blade IceDetection (Cattin and Heikkilä, 2016; fos4X data sheet, n.d.), and Wolfel IDD.Blade (Wölfel, 2016) detect ice based on the deviations in blade’s natural frequencies. These systems can detect icing while the turbine is in operation or at a standstill condition, which enables automatic restart of the turbine once the blades are ice-free after de-icing. These systems indicate the state of icing on the blades such as ice-free, non-critical ice, and critical ice, but cannot identify the location and quantity of ice mass. The most widely used anti-icing and de-icing systems work on the basis of heating resistance or hot air blowing techniques (Ilinca, 2011). Both these systems require an external power source (with a capacity of 150–225 kW for a typical 3-MW range wind turbine (Nielsen, 2017; Roloff, 2017)) and consume 1%–4% of annual energy production, depending on the icing severity (Peltola et al., 2003). The required power to remove ice can be minimized if the location and quantity of ice mass are determined approximately so that the relevant energy needed at the appropriate locations of the blade can be supplied to remove ice. In addition to that, the de-icing process can be initiated and monitored remotely without the need of turbine maintenance personnel’s physical presence at the site. This has motivated the authors to pursue an idea to determine the location and quantity of ice mass accumulated on wind turbine blades.

Fan and Qiao (2011) summarized various damage identification methods for beam- or plate-type structures that use different vibration features (natural frequency–based methods, mode shape–based methods, curvature mode shape–based methods, and methods using both mode shapes and frequencies). These methods rely on the effect of damage induced changes in the physical properties of the structure (mass, damping, and stiffness) on its modal properties. If the changes in physical properties are known, their effect on modal properties can be estimated using either mathematical models of the structural behavior (like finite element method (FEM) models) or experiments (like modal analysis), which is depicted as a forward problem in Figure 1(a). Damage identification methods rely on the inverse problems (as shown in Figure 1(b)) where the location and severity of a damage in the structure are determined based on its modal properties. The authors of this article used an artificial neural network (ANN) model in Gantasala et al. (2017) to solve such inverse problem for identifying added masses on a cantilever beam for any given set of natural frequencies of the corresponding beam. Initially, a data set of natural frequencies of the beam is created using its FEM model where different quantities of added masses are considered at different locations on the beam. A neural network model is trained with this data set considering natural frequencies of the beam as an input and corresponding added masses used in the calculations as an output. It approximates a nonlinear relationship between these inputs and outputs which can be used to identify added masses for any given set of natural frequencies of the beam consisting added masses. This method has a limitation that it requires an updated FEM model (eigenvalues of the FEM model closely match with those estimated from experimental modal analysis of the structure) to create the data set required for training the neural network. The requirement of an FEM model of the structure makes it practically difficult to use such technique for ice detection on wind turbine blades as its FEM model is a proprietary information of the manufacturer. In this work, the authors propose a new method to create the data set required for training the neural network. Experimental modal analyses (output only) are conducted on a small-scale wind turbine blade considering two different quantities of ice mass at three different locations along the blade (ice mass is considered in one location at a time). Later, the relations between natural frequencies of the wind turbine blade and the quantity of ice mass used in each location are obtained using polynomial fitting (natural frequencies are considered as dependent variables and ice masses are considered as independent variables). These relations are used to create a data set of blade’s natural frequencies considering different quantities of ice masses at three different locations along the blade. A neural network model is trained with that data set and used to identify ice masses for a given set of natural frequencies of the blade estimated from experiments (different ice masses are considered randomly along the three different locations defined on the blade). Those ice masses are compared with the actual ice masses considered in the experiments.

Figure 1.

(a) Forward problem and (b) inverse problem.

Experimental setup and modal analysis

A 1-kW wind turbine blade as shown in Figure 2 is considered in this study and its root is rigidly fixed to a stationary support. Geometric details of the blade are given in Table 1. Four uni-axial accelerometers (Bruel & Kjær Type 4507) are placed on the blade to measure vibration accelerations in the flapwise direction using a National Instruments NI-9234 data acquisition card along with the DEWESoft7 data acquisition software.

Figure 2.

Experimental modal analysis setup.

Table 1.

Details of the 1-kW wind turbine blade.

Length (m)	0.97
Weight (kg)	1.263
Maximum chord length (m)	0.11
Minimum chord length (m)	0.055

Two different output only modal identification techniques are used in this section to extract modal parameters of the blade from experimentally measured vibration accelerations (excitation forces are not measured). In the first method, the blade is excited with a low-pass filtered (up to 200 Hz, as first three flapwise modes of the blade are within this frequency range) Gaussian white noise force using a vibration shaker (TIRAvib S-50018) as shown in Figure 2, and its vibration accelerations are measured using four accelerometers with a sampling frequency of 2048 Hz for 200 s. The measured signals for a short duration of 4 s are shown in Figure 3. Modal parameters of the blade (natural frequencies and mode shapes) are extracted using a frequency-domain decomposition (FDD) technique on the measured vibration accelerations. It identifies modal parameters in the frequency domain using a non-parametric model. The measured vibration signals are divided into segments to calculate power spectral density (PSD) estimates using Welch’s overlapped segment averaging estimator. The PSD matrix calculated at discrete frequencies within the sampling frequency is decomposed by means of a singular value decomposition (SVD). These singular values are the estimates of the auto spectral densities of the single-degree-of-freedom (SDOF) systems, whereas singular vectors are estimates of the structural mode shapes. The PSD matrix at any frequency is mainly dominated by the vibration modes around that frequency. Near the resonance frequencies of the structure, the SVD of the PSD matrix will return a predominant singular value signifying a single vibration mode contributes to the PSD matrix at that frequency and the corresponding mode shape is characterized by its singular vector. The authors used an open-source MATLAB code written by Cheynet (2016) for modal parameter identification using the FDD, and further theoretical details on the FDD technique can be found in Brincker et al. (2001). The singular values calculated in the current case for the PSD matrix at various frequencies are shown in Figure 4, where first three flapwise natural frequencies of the blade can be identified at the peaks in the first singular value.

Figure 3.

Vibration accelerations of the blade for a low-pass filtered Gaussian white noise excitation.

Figure 4.

Singular values of the PSD matrix.

In the second modal identification method, the blade is excited with an impact hammer and its free vibration responses are measured. Impact hammer is only used to excite free vibrations of the blade, and its forces are not measured. After an impact hit on the blade, data acquisition system is triggered to acquire free vibration responses with a sampling frequency of 2048 Hz for 6 s (free vibration responses are decaying within this duration). The measured vibration accelerations for a short duration are shown in Figure 5. These measured vibration accelerations are used in an open-source Structural Modal Identification Toolsuite (SMIT) written by Chang et al. (2012). That toolsuite consists of five different time-domain-based output only system identification algorithms (their theoretical background is presented in Chang and Pakzad (2014)). In this study, subspace state space system identification (N4SID-OO) algorithm available in that toolsuite is chosen for extracting modal parameters of the blade. A state space model of certain order is assumed to model the vibration behavior of the test structure in the subspace method where it identifies the state space matrices corresponding to the test structure based on its discrete time-sampled vibration signals. The measured vibration signals are arranged in the form of a Hankel matrix with an even number of rows, and it is further divided into two sub-matrices (consisting equal number of rows), where the upper sub-matrix is called as past output and the lower sub-matrix is called as future output. Using an SVD on the projection of row space of future outputs onto the row space of the past outputs, state space matrices corresponding to the test structure are obtained, from which modal parameters are extracted. This algorithm is not described here in detail as it is out of the scope of this article, more details on this algorithm can be found in Van Overschee and De Moor (1994).

Figure 5.

Free vibration accelerations of the blade measured using four sensors.

A stability diagram is plotted in Figure 6 to identify true poles of the system for different model orders. Three dominant peaks can be identified at 10.1698, 46.1688, and 125.96 Hz in the PSD plot which matches with the stable poles (for all model orders above 10) identified in the subspace identification algorithm, and they correspond to the first three flapwise modes of the blade.

Figure 6.

Stability diagram for the modes of the system.

The natural frequencies identified using both modal identification methods are compared in Table 2, and their mode shapes are plotted in Figure 7. The modal parameters identified using two modal identification methods match very well with each other. The difference between these identified frequencies is marginal, and the mode shapes are only deviating at the location of sensor 2 in second and third modes. The authors are only interested in the natural frequencies of the blade, and mode shapes are extracted for confirming those frequencies correspond to first three vibration modes of the blade. One of these two identification methods can be used for further experiments with ice on the blade, the authors chose time-domain-based modal identification technique compared to the other method as it is more convenient and quicker to measure free vibration responses (for a duration of 6 s only) of the blade.

Table 2.

Comparison of natural frequencies of the 1-kW wind turbine blade identified using FDD and N4SID-OO modal identification methods.

Mode	Natural frequency (Hz)
Mode	FDD	N4SID-OO
1st flapwise	10.16	10.1698
2nd flapwise	46.13	46.1688
3rd flapwise	126.1	125.9600

FDD: frequency-domain decomposition.

Figure 7.

Comparison of mode shapes identified using FDD and N4SID-OO modal identification methods.

Experimental modal analysis of the blade with ice mass

Three different zones of length 0.2 m are defined along the blade as shown in Figure 8, where two quantities of ice masses 40 and 80 g are placed (randomly distributed) in each of these three zones separately (refer to Figure 9(a) to (c)). These ice masses are not naturally formed on the blade, instead ice formed during the winter season is collected and placed on the upper side of the blade’s flapwise surface. Experimental modal analysis is carried out in the six cases described before (two quantities of ice masses are considered in three zones), and the measured free vibration responses are used in the Subspace State Space System Identification (N4SID-OO) algorithm (Chang et al., 2012) to extract corresponding modal parameters of the blade. First three flapwise natural frequencies of the blade identified in these six cases are given in Table 3. These natural frequencies are reduced differently depending on the location and the quantity of ice mass. A quadratic polynomial as shown in Figure 10 is fitted separately between each one of the three natural frequencies of the iced blade corresponding to the cases of ice mass considered in a specific zone and respective ice masses (0, 40, and 80 g) used in that zone. In total, nine polynomials are fitted between three natural frequencies of the iced blade and respective ice masses used in the three zones. These polynomials can be expressed in the form of equation (1), and they can be used to interpolate natural frequencies of the iced blade corresponding to any intermediate ice mass value between 0 and 80 g. The reduction in natural frequencies of the blade for any ice mass value between 0 and 80 g in any one of the three zones exclusively can be calculated using the polynomial coefficients (their values in the current case are given in Appendix 1) given in equations (2) to (4)

f_{i}^{M_{j}} = A_{i j} M_{j}^{2} + B_{i j} M_{j} + f_{i}^{0}

(1)

where $f_{i}^{M_{j}}$ refers to the $i th$ flapwise natural frequency of the blade corresponding to the cases of ice mass values $M = 0, 40, and 80 g$ in $j th$ zone, $A_{i j}$ and $B_{i j}$ are the coefficients of quadratic polynomial fitted between natural frequencies $f_{i}^{M_{j}}$ and ice mass $M_{j}$ , $f_{i}^{0}$ refers to the $i th$ flapwise natural frequency of the blade without any ice mass

\begin{matrix} {\begin{matrix} Δ f_{1} \\ Δ f_{2} \\ Δ f_{3} \end{matrix}}^{(m_{1}, 0, 0)} = {\begin{matrix} A_{11} \\ A_{21} \\ A_{31} \end{matrix}} m_{1}^{2} + {\begin{matrix} B_{11} \\ B_{21} \\ B_{31} \end{matrix}} m_{1}, & \forall m_{1} \in [0, 80] refers to ice mass (in g) in Zone 1 \end{matrix}

(2)

\begin{matrix} {\begin{matrix} Δ f_{1} \\ Δ f_{2} \\ Δ f_{3} \end{matrix}}^{(0, m_{2}, 0)} = {\begin{matrix} A_{12} \\ A_{22} \\ A_{32} \end{matrix}} m_{2}^{2} + {\begin{matrix} B_{12} \\ B_{22} \\ B_{32} \end{matrix}} m_{2}, & \forall m_{2} \in [0, 80] refers to ice mass (in g) in Zone 2 \end{matrix}

(3)

\begin{matrix} {\begin{matrix} Δ f_{1} \\ Δ f_{2} \\ Δ f_{3} \end{matrix}}^{(0, 0, m_{3})} = {\begin{matrix} A_{13} \\ A_{23} \\ A_{33} \end{matrix}} m_{3}^{2} + {\begin{matrix} B_{13} \\ B_{23} \\ B_{33} \end{matrix}} m_{3}, & \forall m_{3} \in [0, 80] refers to ice mass (in g) in Zone 3 \end{matrix}

(4)

where ${\begin{cases} Δ f_{1} \\ Δ f_{2} \\ Δ f_{3} \end{cases}}^{(m_{1}, 0, 0)}$ , ${\begin{cases} Δ f_{1} \\ Δ f_{2} \\ Δ f_{3} \end{cases}}^{(0, m_{2}, 0)}$ , and ${\begin{cases} Δ f_{1} \\ Δ f_{2} \\ Δ f_{3} \end{cases}}^{(0, 0, m_{3})}$ refer to the reduction in first three flapwise natural frequencies of the blade when ice masses $m_{1}, m_{2}, and m_{3}$ are exclusively present in one of the three zones.

Figure 8.

Definition of three zones along the blade length.

Figure 9.

Ice masses considered along the blade in (a) Zone 1, (b) Zone 2, (c) Zone 3, and (d) Zones 1, 2, and 3.

Table 3.

Influence of ice mass on first three flapwise natural frequencies of the 1-kW wind turbine blade.

Reference ice mass case no.	Ice mass		Natural frequency (Hz)
Reference ice mass case no.	Location	Value (g)	Mode 1	Mode 2	Mode 3
Without ice mass			10.1698	46.1688	125.9600

1	Zone 1 Zone 2 Zone 3	4000	10.1010 (−0.68%)	44.9420 (−2.66%)	125.0370 (−0.73%)

2	Zone 1 Zone 2 Zone 3	8000	10.0650 (−1.03%)	43.6090 (−5.54%)	124.1370 (−1.45%)

3	Zone 1 Zone 2 Zone 3	0400	9.9540 (−2.12%)	45.2750 (−1.94%)	120.7140 (−4.16%)

4	Zone 1 Zone 2 Zone 3	0800	9.6640 (−4.97%)	44.4490 (−3.73%)	116.4130 (−7.58%)

5	Zone 1 Zone 2 Zone 3	0040	9.5570 (−6.03%)	45.5020 (−1.44%)	124.2480 (−1.36%)

6	Zone 1 Zone 2 Zone 3	0080	9.0038 (−11.47%)	45.0259 (−2.48%)	122.4027 (−2.82%)

The percentage change in the natural frequencies of the iced blade is calculated with respect to the case of without any ice mass on the blade.

Figure 10.

Quadratic polynomial curve fit for the $i th$ natural frequency with ice mass in $j th$ zone.

Using equations (2) to (4), the natural frequencies of blade with ice masses $m_{1}, m_{2}, and m_{3}$ altogether in the Zones 1, 2, and 3, respectively, are approximated in equation (5) using a linear superposition of the reduction in natural frequencies when those ice masses are present individually on the blade. The validity of this approximation is verified by conducting an experimental modal analysis of the blade considering 40-g ice mass in all three zones as shown in Figure 9(d). The natural frequencies identified in this case using N4SID-OO system identification algorithm on the vibration measurements of the blade are shown in Table 4 and also compared with the frequency values calculated using equation (5) with the following ice mass values (in g) $m_{1} = 40, m_{2} = 40, and m_{3} = 40$ . The natural frequencies identified from experimental modal analysis approximately matches with those calculated using equation (5).

{\begin{matrix} f_{1} \\ f_{2} \\ f_{3} \end{matrix}}^{(m_{1}, m_{2}, m_{3})} \approx {\begin{matrix} f_{1}^{0} \\ f_{2}^{0} \\ f_{3}^{0} \end{matrix}} + {\begin{matrix} Δ f_{1} \\ Δ f_{2} \\ Δ f_{3} \end{matrix}}^{(m_{1}, 0, 0)} + {\begin{matrix} Δ f_{1} \\ Δ f_{2} \\ Δ f_{3} \end{matrix}}^{(0, m_{2}, 0)} + {\begin{matrix} Δ f_{1} \\ Δ f_{2} \\ Δ f_{3} \end{matrix}}^{(0, 0, m_{3})}

(5)

where ${\begin{matrix} f_{1} \\ f_{2} \\ f_{3} \end{matrix}}^{(m_{1}, m_{2}, m_{3})}$ refers to the first three flapwise natural frequencies of the blade with the ice masses $m_{1}, m_{2}, and m_{3}$ in the Zones 1, 2, and 3, respectively.

Table 4.

Comparison of natural frequencies of the blade with 40-g ice mass in all three zones calculated using equation (5) and identified from the experimental modal analysis.

Ice mass		Flapwise mode	Natural frequency (Hz)
Location	Value (g)	Flapwise mode	Experiment (Figure 9(d))	Using equation (5)
Zone 1	40	Mode 1	9.2580	9.2724
Zone 2	40	Mode 2	43.4300	43.3814
Zone 3	40	Mode 3	118.5000	118.0790

A valid forward relation between ice masses (as independent variables) and natural frequencies (as dependent variables) of the 1-kW wind turbine blade is defined in this section using experimental modal analyses with ice masses. So, first three natural frequencies of the iced blade can be determined using such relation without in need of a structural model like FEM when the quantities of ice masses are known. Similarly, an inverse relation with ice masses as the dependent variables and natural frequencies as the independent variables is required to find the quantities of ice masses whenever natural frequencies of the iced blade are known. It is difficult to express these inverse relations explicitly. However, optimization techniques can be used to serve the purpose of finding the quantities of ice masses which will reduce the natural frequencies in the similar way as that of measured frequencies of an iced blade. A cost function can be defined using the forward relation and target natural frequencies where ice masses are chosen as variables in the optimization techniques and tuned to match the frequencies calculated using the forward relation with target frequencies. This process needs to be repeated for every new set of target frequencies to find corresponding ice masses which makes it cumbersome to use for ice mass detection. Another viable option is to define the inverse relation between ice masses and natural frequencies using an artificial neural network (ANN) model which expresses that relation explicitly, so it is straightforward to use that relation later for detecting ice masses for any new input of natural frequencies of the iced blade.

ANN model

ANNs are inspired from the biological neural networks, where artificial nodes known as neurons are used as computational elements and connected together to form a network. The ANN can be used to estimate or approximate functions that depend on multiple variables. These networks receive input from the neurons in the input layer, and their output are given by the neurons in the output layer. There may be one or more intermediate hidden layers. The neuron behaves as an activation or mapping function (typically a nonlinear function) producing an output corresponding to the input to the neuron. The neurons in the network are connected with each other and associated with some weights and bias which either increase or decrease input signals to the neurons. These weights and bias are determined in the training process where the relation between input and output variables is established (Ata, 2015). The backpropagation algorithm is one of the most widely used algorithms to determine the weights and bias of the neurons organized in layers. In this algorithm, the network signal travels in the forward direction and the errors are propagated backward. This algorithm uses supervised learning in which the network is trained using data for which inputs as well as desired outputs are known. The backpropagation algorithm is a generalization of the least mean square algorithm that modifies network weights to minimize the mean square error $(m s e)$ between the desired and actual outputs of the network. Once trained, the network weights are frozen and can be used to predict output values for any new input samples. The advantage of ANN is that it only needs inputs and outputs of the process to approximate underlying complex nonlinear functions. Other conventional methods like regression analysis require an explicit function definition between the inputs and outputs, whereas the ANN does not require such function definition.

A simple network with two neurons in a hidden layer is shown in Figure 11. The input to the network is $x$ and output from the network is $y$ . The connections between neurons from the input to hidden layer are associated with the weights $w_{1} and w_{2}$ , and the neurons in the hidden layer are associated with the biases $b_{1} and b_{2}$ . The neurons from hidden to output layer are associated with the weights $w_{3} and w_{4}$ and the neuron in the output layer is associated with a bias $b_{3}$ . The weights of neurons either increase or decrease input signals to the neuron whereas bias has the effect of increasing or lowering the net input to the neuron. An input pattern is presented to the input layer neurons that pass them onto the hidden layer, and the neurons in the input layer do not perform any computations. Each of the hidden layer neurons computes a sum of the weighted inputs and biases $(v_{i}, i = 1, 2)$ , passes them through their activation function, and presents the result to the output layer. These will be the inputs to neuron in the output layer which compute their weighted sum and then feed to the activation function to compute the output of the network $(y)$ . The backpropagation algorithm determines these weights $w_{i}, i = 1, \dots, 4$ and biases $b_{i}, i = 1, \dots, 3$ such that the difference between network output $(y)$ and the expected value of training data is the minimum, which will be problem specific. The network shown in Figure 11 is a simple network that maps the relation between $x$ and $y$ using two neurons; in general, ANNs are used for function approximation for a data with multiple inputs and outputs. Such applications require multiple neurons organized in layers between the input and output layers. The number of input and output neurons in the network is dictated by the problem under consideration. There is no hard rule to follow for choosing the number of hidden layers, number of neurons in those layers and the number of samples required in the training data set. These parameters are determined by trial and error method for the problem under consideration.

Figure 11.

A simple artificial neural network model.

In this study, a nonlinear function that considers natural frequencies of the iced blade as an input and the ice masses accumulated in the three zones defined on the blade in corresponding cases as an output is approximated using ANN. As discussed before, the network parameters need to be found by trial and error method, which makes the process of finding network parameters cumbersome. However, based on the training data and problem under consideration, these parameters can be chosen, which is outlined as follows:

Number of outputs $(n_{o / p})$ : this value is fixed at 3 in this study, which are the ice masses in the three zones defined on the blade as shown in Figure 8. Increasing this number further makes the function approximation problem more complex, which requires more natural frequencies as an input to find a valid relation between inputs and outputs of the neural network.

Number of inputs $(n_{i / p})$ : for the current problem, natural frequencies of the blade are used as an input to the neural network. The accuracy and consistency of the ice mass detection are influenced by the number of natural frequencies used in the detection technique. While designing the neural network for the current problem, it is found that a minimum of three natural frequencies is required as an input to the network to establish a valid relation between inputs and outputs of the network. Also, the uncertainty associated with multiple possibilities (in contrast to local and global minima in the optimization techniques) of the ice mass combinations that have a similar reduction in blade natural frequencies is minimized when three frequencies are used as an input. So, first three flapwise natural frequencies of the blade are chosen as an input to the neural network.

Number of hidden layers $(n_{h i d d e n})$ : neural networks with one hidden layer and sigmoid function as an activation function can approximate continuous functions with arbitrary precision (Mehrotra et al., 1996). So, one hidden layer of neurons is chosen in this study.

Number of neurons $(n_{n e u})$ : when neural networks are used for classification tasks with d-dimensional input data, neurons in the hidden layer function as hyperplanes that partition d-dimensional space into various regions. For classifying d-dimensional input data into $M$ clusters using a single-layer network, a minimum of the $M$ number of hidden layer neurons is required (Mehrotra et al., 1996). If two classes are defined for the ice mass like light and heavy with respect to the mass of the blade in each zone, then eight classes of combinations of ice masses are possible in three zones. For such definition of ice mass classes, a minimum of eight neurons are required in the hidden layer to approximate the nonlinear function using ANN.

Number of samples in the data set used for training $(n_{s a m p l e})$ : a rule of thumb, obtained from related statistical problems, is to have at least 5–10 times as many training samples as the number of weights and bias variables used in the network to be trained (Mehrotra et al., 1996). If $n_{n e u}$ number of neurons are considered in the hidden layer, then the number of weights and bias variables to be determined in the training process are given by $n_{v a r i a b l e s} = (n_{i / p} + 1) n_{n e u} + n_{o / p} (n_{n e u} + 1)$ . In this study, the number of samples considered in the training data is such that they are greater than 5 or 10 times the number of variables to be determined in the network. Each sample in the data set consists of $n_{i / p}$ inputs and $n_{o / p}$ outputs.

Gantasala et al. (2017) found $n_{n e u} = 8$ and $n_{s a m p l e} = 343$ as suitable parameters for a similar problem, where a nonlinear function is approximated for the data set consisting three inputs (first three natural frequencies of a cantilever beam) and three outputs (added masses at three different locations on the beam). Same values are used for $n_{n e u}$ and $n_{s a m p l e}$ in this study, and a data set consisting 343 samples is created for training the neural network model. A total eight different cases are possible as shown in Figure 12 when two possibilities of ice mass (0 and some value $m$ ) are chosen in each zone. So, seven values between 0 and 80 g are chosen in each of the three zones to create the required data set of 343 samples (equal to $7^{3}$ ). For these ice mass values, natural frequencies of the blade are calculated using equation (5). A data set is created with the calculated natural frequencies and ice mass values used in the corresponding cases. A neural network model is designed using the MATLAB toolbox (MATLAB Neural Network Toolbox, n.d.) with the network parameters shown in Figure 13 by training with a data set of 343 samples considering first three flapwise natural frequencies of the iced blade as an input and corresponding ice masses used in the three zones as an output. An approximate inverse relation between these inputs and outputs is identified in the training process. The designed neural network model is given natural frequencies of the iced blade identified from experimental modal analyses as an input in the next section to identify ice masses.

Figure 12.

Different cases possible when two values are chosen for an ice mass in each one of the three zones.

Figure 13.

Neural network designed using the MATLAB toolbox.

Results and discussion

In this section, 18 cases of ice masses (randomly chosen between 0 and 80 g in all three zones) are considered along the length of the three zones defined on the 1-kW wind turbine blade, and its first three flapwise natural frequencies are identified from the experimental modal analysis. These ice masses are distributed randomly (i.e., linear mass density varies randomly) along the length of the zones. The natural frequencies identified in the above cases (given in Appendix 2) are given as an input to the neural network model designed in the last section, and the identified ice masses are compared against the actual ice mass values used in the experiments in Table 5. As the actual ice mass values used in the experiment and the corresponding identified mass values in all 18 test cases are known, an error metric called as weighted absolute percentage error $(WAPE)$ as defined in equation (6) is calculated for each test case. It is defined as the ratio of the sum of the absolute error in the identified masses in each zone to the sum of actual ice masses used in that test case

WAP E_{i} = \frac{\sum_{k = 1}^{3} | m_{a c t u a l}^{i k} - m_{i d e n t i f i e d}^{i k} |}{\sum_{k = 1}^{3} m_{a c t u a l}^{i k}}

(6)

where $i$ refers to the test case number; $k$ refers to the zone number; $m_{a c t u a l} and m_{i d e n t i f i e d}$ refers to ice masses actually used in the experiment and the values identified by neural network model, respectively.

Table 5.

Ice masses identified in the 18 test cases: 1-kW wind turbine blade.

Test case no.	Actual ice mass values (g)			Identified ice mass values (g)
Test case no.	Zone 1	Zone 2	Zone 3	Zone 1	Zone 2	Zone 3
1	0.0	0.0	15.0	0.0	0.0	18.5
2	0.0	0.0	30.0	0.0	0.0	32.1
3	0.0	0.0	48.0	0.0	0.0	45.7
4	0.0	10.0	0.0	0.4	8.5	3.7
5	0.0	25.0	20.0	0.0	19.5	22.6
6	0.0	36.0	0.0	0.0	31.6	5.0
7	0.0	36.0	15.0	0.0	29.3	21.8
8	10.0	0.0	30.0	10.5	0.0	32.8
9	10.0	15.0	30.0	8.1	8.7	33.9
10	12.0	0.0	0.0	15.0	0.0	0.0
11	12.0	0.0	26.0	11.0	0.0	25.8
12	12.0	22.0	0.0	12.2	17.8	3.0
13	15.0	25.0	20.0	17.2	19.6	23.0
14	25.0	0.0	0.0	28.0	0.0	0.0
15	25.0	0.0	48.0	21.3	0.0	45.4
16	25.0	32.0	0.0	17.9	31.9	4.6
17	25.0	32.0	48.0	18.6	29.6	47.7
18	40.0	40.0	40.0	40.3	35.4	42.9

The WAPE is calculated for all 18 test cases and shown in Figure 14. This error changes with respect to the location and the quantity of ice mass, and the mean value of WAPE for all 18 cases is found to be 17.53%. Ice masses in the test case 4 are identified with a maximum WAPE value of 55.6%, where only 10 g of ice mass is considered in Zone 2 in the experiment. It can be observed that the WAPE value is higher in test cases where smaller ice masses are considered on the blade and it decreases with increasing ice mass values. The randomness in the way ice masses are distributed in the experiments resulted randomness in the WAPE values of 18 test cases. The proposed ice mass detection algorithm is able to approximately identify ice masses in the three zones defined on the blade based on its first three flapwise natural frequencies. Its application to a large wind turbine blade is demonstrated in the next section.

Figure 14.

Weighted absolute percentage error $(WAPE)$ calculated for the 18 test cases.

Application of the proposed technique to detect ice masses on a multi-megawatt wind turbine blade

The detection technique proposed in the earlier sections is approximately able to identify ice masses on a 1-kW wind turbine blade, the same technique is used in this section for detecting ice masses on the NREL 5-MW wind turbine blade. In the case of 1-kW wind turbine blade, an explicit relation between natural frequencies and ice masses is defined using experimental modal analyses of the iced blade; here, in this section, such relation for the NREL 5-MW wind turbine blade is defined using an eigenvalue analysis of its FEM model with ice mass. The structural details of the NREL 5-MW wind turbine blade can be found in Jonkman et al. (2009). The length of the wind turbine blade is divided into three zones as shown in Figure 15. The inner half of the blade is considered as Zone 1 and the outer half of the blade is divided into two zones, Zones 2 and 3. Two quantities of ice masses 250 and 500 kg are considered with a constant linear mass density (mass per unit length) along the length of the three zones separately as shown in Figure 15, and eigenvalues of the blade (first three flapwise natural frequencies are of interest in this study) in the corresponding cases are calculated using the BModes tool (Bir, 2005). A quadratic polynomial is fitted separately between the ice masses (0, 250, and 500 kg) used in a particular zone and respective flapwise natural frequencies of the blade (one of the three flapwise natural frequencies is considered at a time). The polynomial fitted between third flapwise natural frequency of the iced blade when ice masses (0, 250, and 500 kg) are considered exclusively in Zone 3 is shown in a figure in Appendix 3. This particular natural frequency calculated with different intermediate ice mass values between 0 and 500 kg in Zone 3 is also shown in the same figure in Appendix 3, these frequencies are following the quadratic polynomial fitted before. This emphasizes that the natural frequencies of a wind turbine blade vary quadratically with the quantity of ice mass in a specific location (one of the three zones). So, the polynomial fitted in Appendix 3 can be used to calculate third flapwise natural frequency of the blade corresponding to any ice mass value between 0 and 500 kg in Zone 3. Similarly, polynomials are fitted for other two natural frequencies of the blade with respect to ice mass in Zone 3. Likewise, polynomials are fitted for ice masses in other two zones of the blade. Using these polynomial coefficients, an explicit relation like equation (5) is defined between natural frequencies of the 5-MW wind turbine blade and ice masses in the three zones $(m_{1}, m_{2}, and m_{3})$ . The validity of this relation for the NREL 5-MW wind turbine blade is verified using an eigenvalue analysis of its FEM model with ice masses in Appendix 4. Later, this explicit relation is used to create a data set of natural frequencies of the blade corresponding to different ice masses in the three zones.

Figure 15.

Division of the NREL 5-MW wind turbine blade length into different zones.

Seven different ice mass values are considered in the range of 0–500 kg in each zone, which will generate 343 (equal to $7^{3}$ ) combinations of ice masses in the three zones, and natural frequencies in all these cases are calculated using the explicit relation defined earlier to create the data set required for training a neural network model. An ANN model (with parameters $n_{n e u} = 8 and n_{h i d d e n} = 1$ ) is trained with this data set considering first three flapwise natural frequencies of the iced blade as an input and three ice masses used in the corresponding cases as an output. ANN approximates a nonlinear relation between these inputs and outputs. Ice masses are identified for the six test cases shown in Figure 16 by giving natural frequencies of the blade (calculated using BModes tool) in the respective cases as an input to the trained neural network model. The ice masses identified from the neural network model are compared with actual ice mass values used in the test cases in Table 6. The detection technique proposed in this work is approximately able to identify ice masses on the 5-MW wind turbine blade based on its natural frequencies, and these ice masses are identified with a mean WAPE value of 14.77% in six test cases.

Figure 16.

Ice mass distributions considered along the NREL 5-MW wind turbine blade in the six test cases.

Table 6.

Ice masses identified in the 6 test cases: NREL 5 MW wind turbine blade.

Test case no.	Ice mass location	Ice mass (kg)		WAPE
Test case no.	Ice mass location	Actual	Identified	WAPE
1	Zone 1 Zone 2 Zone 3	785048	955644	15.56%

2	Zone 1 Zone 2 Zone 3	15510096	14411584	10.59%

3	Zone 1 Zone 2 Zone 3	223157148	212209124	16.71%

4	Zone 1 Zone 2 Zone 3	315197193	273226166	13.95%

5	Zone 1 Zone 2 Zone 3	376252253	299318213	20.79%

6	Zone 1 Zone 2 Zone 3	464300293	468363243	11.01%

WAPE: weighted average percentage error.

Summary and conclusion

Ice accumulation on a wind turbine blade reduces its natural frequencies differently depending on the location and the quantity of ice mass. In this work, a detection technique is proposed to identify ice masses on a blade based on its first few natural frequencies. The proposed technique is initially demonstrated on a small 1-kW wind turbine blade experimental setup. Experimental modal analyses are carried out on that blade to determine the sensitivities of its natural frequencies to the quantity of ice mass and its location on the blade. An explicit relation to the natural frequencies of the blade is defined in terms of the location and quantity of ice masses using those sensitivities. Later, a data set of natural frequencies of the blade is created using that explicit relation considering different quantities of ice masses at different locations on the blade. An ANN model is trained with this data set where natural frequencies of the iced blade are used as an input and corresponding ice masses at different locations on the blade as an output of the neural network model. ANN approximates a nonlinear relation between these inputs and outputs in the training process. This model is able to approximately identify ice masses with a WAPE value of 17.53% (mean value) in the 18 test cases based on their natural frequencies.

The proposed technique can be applied to detect ice masses on large wind turbines operating in cold climate as it requires only first few natural frequencies of the blade, these frequencies are usually excited by the turbulent wind and can be estimated from the vibration measurements of the blade. The effectiveness of the proposed technique is demonstrated on the NREL 5-MW wind turbine blade where eigenvalue analysis of its FEM model is used to define an explicit relation between its natural frequencies and ice masses on the blade. A neural network model designed using that explicit relation is able to identify ice masses with a mean WAPE value of 14.77% in six test cases. This technique can be used to identify ice masses on a stationary wind turbine, which will allow to optimize and control the de-icing process remotely with the help of information on location and approximate quantities of ice masses present on the blade.

Footnotes

Appendix 1: Values of the polynomial coefficients in equations (2) to (4) for the 1-kW wind turbine blade

\begin{array}{l} [\begin{matrix} A_{11} & A_{12} & A_{13} \\ A_{21} & A_{22} & A_{23} \\ A_{31} & A_{32} & A_{33} \end{matrix}] = 10^{- 5} [\begin{matrix} 1.025 & - 2.319 & 1.862 \\ - 3.319 & 2.119 & 5.959 \\ 0.719 & 29.530 & - 4.166 \end{matrix}] \\ [\begin{matrix} B_{11} & B_{12} & B_{13} \\ B_{21} & B_{22} & B_{23} \\ B_{31} & B_{32} & B_{33} \end{matrix}] = 10^{- 2} [\begin{matrix} - 0.213 & - 0.447 & - 1.607 \\ - 2.934 & - 2.319 & - 1.905 \\ - 2.336 & - 14.300 & - 4.113 \end{matrix}] \end{array}

Appendix 2: Natural frequencies of the iced 1-kW wind turbine blade identified in the 18 test cases

Test case no.	Ice masses (g)	Natural frequencies (Hz)
Test case no.	$m_{1}, m_{2}, and m_{3}$	Mode 1	Mode 2	Mode 3
1	(0, 0, 15)	9.89	45.97	125.30
2	(0, 0, 30)	9.69	45.71	125.10
3	(0, 0, 48)	9.48	45.55	124.10
4	(0, 10, 0)	10.07	45.89	124.60
5	(0, 25, 20)	9.73	45.41	122.40
6	(0, 36, 0)	9.93	45.46	121.60
7	(0, 36, 15)	9.69	45.24	121.20
8	(10, 0, 30)	9.67	45.42	125.10
9	(10, 15, 30)	9.59	45.15	123.10
10	(12, 0, 0)	10.15	45.74	125.70
11	(12, 0, 26)	9.75	45.42	124.80
12	(12, 22, 0)	10.01	45.34	123.10
13	(15, 25, 20)	9.68	44.80	121.90
14	(25, 0, 0)	10.13	45.37	125.60
15	(25, 0, 48)	9.44	44.80	123.60
16	(25, 32, 0)	9.90	44.83	121.10
17	(25, 32, 48)	9.26	44.17	119.50
18	(40, 40, 40)	9.26	43.43	118.50

$m_{1}, m_{2}, and m_{3}$ refer to the ice masses considered in Zones 1, 2, and 3, respectively, in the test cases.

Appendix 3: Quadratic polynomial fitted between third flapwise natural frequency of the iced NREL 5-MW wind turbine blade (when ice mass is only considered in Zone 3) and respective ice masses used in the Zone 3

Appendix 4: Validation of the approximation used in equation (5) for the NREL 5-MW wind turbine blade natural frequencies with ice masses m 1,m 2,a n d m 3 in Zones 1,2,and 3,respectively

Values of the polynomial coefficients in equations (2) to (4) corresponding to the NREL 5-MW wind turbine blade (with ice masses) natural frequencies are given below

\begin{array}{l} [\begin{matrix} A_{11} & A_{12} & A_{13} \\ A_{21} & A_{22} & A_{23} \\ A_{31} & A_{32} & A_{33} \end{matrix}] = [\begin{matrix} 0.00 E + 00 & 2.29 E - 09 & 6.84 E - 08 \\ 2.72 E - 10 & 2.26 E - 08 & 2.01 E - 07 \\ 3.90 E - 09 & 2.14 E - 08 & 3.38 E - 07 \end{matrix}] \\ [\begin{matrix} B_{11} & B_{12} & B_{13} \\ B_{21} & B_{22} & B_{23} \\ B_{31} & B_{32} & B_{33} \end{matrix}] = [\begin{matrix} - 1.20 E - 06 & - 3.63 E - 05 & - 2.09 E - 04 \\ - 2.27 E - 05 & - 1.33 E - 04 & - 4.12 E - 04 \\ - 8.74 E - 05 & - 2.14 E - 04 & - 9.22 E - 04 \end{matrix}] \end{array}

First three flapwise natural frequencies of the blade calculated using an eigenvalue analysis of its FEM model with ice masses (in kg) $m_{1} = 250, m_{2} = 250, and m_{3} = 250$ are compared with the frequencies calculated using equation (5) with above polynomial coefficients in the below table.

From the above table, it is clear that the natural frequencies of iced NREL 5-MW wind turbine blade can be approximately calculated using equation (5) with above polynomial coefficients for given quantities of ice masses $(m_{1}, m_{2}, and m_{3})$ .

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: This research work is carried out as a part of the project on “Wind power in cold climates” funded by the Swedish Energy Agency, which supports research and development investments in wind power.

References

Alsabagh

ASY

Virk

et al . (2015) Atmospheric ice loading and its impact on natural frequencies of wind turbines. Wind Engineering 39(1): 83–96.

Ata

(2015) Artificial neural networks applications in wind energy systems: A review. Renewable & Sustainable Energy Reviews 49: 534–562.

Bir

(2005) BModes: Software for computing rotating blade coupled modes. Available at: https://nwtc.nrel.gov/BModes (accessed 25 May 2017).

Brenner

(2016) Determination of the actual ice mass on wind turbine blades. In: Winterwind—International wind energy conference, Åre, 9–10 February.

Brincker

Zhang

Andersen

(2001) Modal identification of output-only systems using frequency domain decomposition. Smart Materials and Structures 10(3): 441.

Cattin

Heikkilä

(2016) Evaluation of ice detection systems for wind turbines. Report, VGB research project no. 392. Available at: https://www.vgb.org/vgbmultimedia/392_Final+report-p-10476.pdf (accessed 25 May 2017).

Chang

Pakzad

(2014) Observer Kalman filter identification for output-only systems using interactive Structural Modal Identification Toolsuite. Journal of Bridge Engineering 19(5). Available at: https://doi.org/10.1061/(ASCE)BE.1943-5592.0000530

Chang

Leonard

Pakzad

(2012) Structural Modal Identification Toolsuite (SMIT). Available at: http://smit.atlss.lehigh.edu (accessed 25 May 2017).

Cheynet

(2016) Automated frequency domain decomposition (AFDD): Modal parameters identification from ambient vibrations measurement using the FDD. MATLAB Central File Exchange. Available at: https://se.mathworks.com/matlabcentral/fileexchange/57153-automated-frequency-domain-decomposition–afdd- (accessed 25 May 2017).

10.

Emerging from the Cold (2016) Available at: http://www.windpowermonthly.com/article/1403504/emerging-cold (accessed 25 May 2017).

11.

Etemaddar

Hansen

MOL

Moan

(2014) Wind turbine aerodynamic response under atmospheric icing conditions. Wind Energy 17: 241–265.

12.

Fan

Qiao

(2011) Vibration-based damage identification methods: A review and comparative study. Structural Health Monitoring 10(1): 83–111.

13.

fos4X data sheet (n.d.) fos4Blade IceDetection. Available at: http://www.fos4x.de/sites/default/files/field_downloads/0307_1-2-0_4-20_01_fos4Blade_IceDetection.pdf (accessed 25 May 2017).

14.

Gantasala

Luneno

J-C

Aidanpää

J-O

(2016) Influence of icing on the modal behavior of wind turbine blades. Energies 9(11): 862.

15.

Gantasala

Luneno

J-C

Aidanpää

J-O

(2017) Investigating how an artificial neural network model can be used to detect added mass on a non-rotating beam using its natural frequencies: A possible application for wind turbine blade ice detection. Energies 10(2): 814.

16.

Ilinca

(2011) Analysis and mitigation of icing effects on wind turbines. In: Al-Bahadly

(ed.) Wind Turbines. InTech, pp. 177–215. DOI: 10.5772/15484

17.

ISO 12494:2001 (2001) Atmospheric icing of structures.

18.

Jonkman

Butterfield

Musial

et al . (2009) Definition of a 5-MW reference wind turbine for offshore system development. Technical report NREL/TP-500-38060. National Renewable Energy Laboratory. Available at: http://www.nrel.gov/docs/fy09osti/38060.pdf

19.

Lorenzo

Manzato

Peeters

et al . (2013) Virtual structural monitoring of wind turbines using operational modal analysis techniques. Key Engineering Materials 569–570: 523–530.

20.

MATLAB Neural Network Toolbox (n.d.) Available at: https://se.mathworks.com/products/neural-network/ (accessed 25 May 2017).

21.

Mehrotra

Mohan

Ranka

(1996) Elements of Artificial Neural Networks. Cambridge, MA: MIT Press.

22.

Nielsen

(2017) Vestas cold climate solutions. In: Winterwind—International wind energy conference, Skellefteå, 6–8 February.

23.

Peltola

Laakso

Antikainen

et al . (2003) Prevention of icing effects. In: Proceedings of BOREAS VI conference, wind energy production in cold climate, Pyhä, Finland, 9–11 April.

24.

Rindeskär

(2010) Modelling of icing for wind farms in cold climate. Report, Examensarbete vid Institutionen fr Geovetenskaper. Available at: http://www.elforsk.se/Global/Vindforsk/Rapporter%20VFIII/V-313/Modelling_icing_Rindesskar_2010.pdf

25.

Rissanen

Lehtomäki

Wennerkoski

et al . (2016) Modelling load and vibrations due to iced turbine operation. Wind Engineering 40(3): 293–303.

26.

Roloff

(2017) Reliable energy production during the winter months. In: Winterwind—International wind energy conference, Skellefteå, 6–8 February.

27.

Skrimpas

Kleani

Mijatovic

et al . (2015) Detection of icing on wind turbine blades by means of vibration and power curve analysis. Wind Energy 19(10): 1819–1832.

28.

Turkia

Huttunen

Wallenius

(2013) Method for estimating wind turbine production losses due to icing. Report, VTT Technology 114. Espoo: VTT.

29.

Van Overschee

De Moor

(1994) N4SID: Subspace algorithms for the identification of combined deterministic-stochastic systems. Automatica 30(1): 75–93.

30.

Wölfel

(2016) Reliable ice detection for rotor blades to increase availability and yield of wind turbines. In: Winterwind—International wind energy conference, Åre, 9–10 February.