Mass and Aerodynamic Imbalance Estimates of Wind Turbines

Due to its effect on the operation time of wind turbines, rotor imbalances of a wind turbine have to be detected early enough. We present a method that determines inhomogeneous mass distributions of the rotor as well as deviations in the pitch angles of the rotor blades from vibrational data only. To this end, a mathematical model connecting the load caused by the imbalances to the resulting vibrations was developed. After discretization, the resulting vibration equation was solved analytically. The inverse problem, i.e., the calculation of the mass and aerodynamic imbalance from vibrational data, was solved by using nonlinear regularization theory. Numerical simulations were performed using artificial vibration data.


Introduction
The efficient operation of wind parks calls for a long life time of each wind turbine.As many defects of a turbine can be related to vibrations of the system, the investigation of vibrations caused by rotor imbalances gets more and more attention.The growing size of new wind turbines leads to a more flexible structure and therefore even bigger vibrational amplitudes.Additionally, the operation of large off-shore wind parks requires a careful remote monitoring of imbalances.A well balanced rotor will prevent early fatigue and ensure a safe and economic operation of the wind turbine.Imbalances affect the drive train components as well as the structural health of the turbine [1].Hence the detection and elimination of imbalances are of vital importance.
There are two main imbalance causes: the so-called mass imbalance arising from inhomogeneous mass distributions caused by, e.g., manufacturing inaccuracies or water inclusions in the blades' texture, and aerodynamic imbalances, arising, e.g., from errors in the pitch angles or profile changes of the blades.In wind turbines with no Condition Monitoring System (CMS), imbalances are only detected if the vibrations of the turbine are clearly visible, which is only the case if the turbine rotates with a frequency close to the bending eigenfrequency of the turbine, or if there are already damages in the drive train components.If the turbine is equipped with a standard CMS, an increased presence of the rotating frequency and its multiples in the Fourier spectrum can indicate rotor imbalances.In both cases an expert team has to be employed to detect the imbalance location and quantity.In a time consuming process, the team first tries to detect aerodynamic imbalances by using mainly optical methods.If the aerodynamic imbalances have been removed, the remaining vibrations are measured.With a known test mass placed on one blade, a second set of measurements is obtained, and from both measurements together the mass imbalances can be calculated.This process requires not only a lot of time but also expensive manpower.Therefore, it is advantageous to use automatic detection methods.Presently, such methods to detect mass and aerodynamic imbalances by using signal processing methods as well as trend analysis to generate an alarm system require a learning phase performed under faultless condition of the rotor [2].However, they do not allow to compute the actual value and position of the imbalance.To this aim, vibration measurements with and without test weights as described above are still necessary.To replace this expensive procedure by a new CMS feature is the main goal of the authors work.
First results in that direction were achieved for the observation and reconstruction of mass imbalances only [3].The mass imbalance reconstruction from noisy vibration measurements was treated as an inverse problem with methods from the regularization theory [4].To solve an inverse problem, first we have to handle the direct problem of computing the turbine vibrations for a given imbalance.This can be done by using a mathematical model of the turbine that solves the vibration equation with the Finite Element Method (FEM).Although there already exist several models for wind turbines, in [3] we developed a new and comparably simple model for two reasons: Firstly, in order to find a solution for the inverse problem the direct problem has to be solved several times.Thus the computations should be very fast.Secondly, the imbalance identification algorithm should be applicable to a large variety of wind turbine types, which means the model should be adaptable without much effort.In short, in our approach the information from the test run with reference imbalances is replaced by the model information, the additional measurements are superficial due to the mathematical model.However, the model proposed in [3] was restricted to mass imbalances and did not consider aerodynamic imbalances.Moreover, we dealt with a one dimensional model only, i.e., we considered lateral vibrations.Beside axial and torsional vibrations, aerodynamic imbalances also cause vibrations in lateral direction.Thus, neglecting aerodynamic imbalances leads to wrong reconstructions.The resulting error depends on the location of the aerodynamic imbalance and on the rotational frequency, as loads from aerodynamic imbalances depend on the frequency.If the aerodynamic imbalance is located close to the mass imbalance the forces in lateral direction add and we will end up with a reconstructed mass imbalance that is too big.On the other hand, if the aerodynamic imbalance is located opposite to the mass imbalance, we will get a mass imbalance reconstruction that is too small.
In the present paper we consider a model that allows for a reconstruction of both mass imbalances and aerodynamic imbalance caused by deviations in pitch angles at the same time.Since aerodynamic imbalances also cause vibrations in axial directions and torsional vibration around the tower axis, we had to expand our turbine model in these dimensions.To describe the forces and moments from aerodynamic imbalances we have used the blade element momentum (BEM) method.Due to the BEM method, the direct problem of relating the imbalance cause (here the pitch angle deviation and a mass imbalance) to the vibrations of the turbine becomes non-linear.Nevertheless, the regularization techniques for solving the inverse problem also apply to non-linear problems.To ensure a good imbalance reconstruction we observed that the initial value should already be a fairly good guess for the mass imbalance, which we obtained by our old method described above, assuming the absence of aerodynamic imbalances.We reconstructed a good guess for the mass imbalance from the lateral vibrations neglecting the influence from the pitch error.The result is surely not the correct mass imbalance but serves very well as an initial value for the reconstruction of both imbalance causes.
The paper is organized as follows: In Section 2, we will present the system of equations connecting vibrations and imbalances.The mass and stiffness matrices are calculated by a model of the turbine that considers vibration in lateral and axial direction as well as torsional vibrations around the tower axis.The description of forces and moments from mass imbalances is relatively easy, the aerodynamic forces and moments are derived with the BEM method.In Section 3, the vibration equation is solved explicitly, whereas in Section 4, the inverse problem of deriving the imbalances from the vibrations is considered.It results in a nonlinear optimization problem.Section 5 contains examples for the numerical calculation of mass and aerodynamic imbalances from artificial vibration data at the top of the tower.

System Equation
Before we can think about reconstructing imbalances from vibration measurements we have to be able to handle the other direction, i.e., for a given imbalance cause we should compute the resulting vibrations of the wind turbine tower.In order to derive appropriate formulas we have to make simplifying assumptions.First we assume that the wind turbine is a flexible shaft (the tower) with an additional mass (the nacelle and the rotor) at the top point, where one part of that mass (the rotor) rotates with a rotational frequency Ω.The movement of such a shaft is described by a Partial Differential Equation (PDE) in time and space.An explicit solution of that equation is seldom possible.Using the Finite Element Method (FEM), the problem can be discretized in space, which results in an Ordinary Differential Equation (ODE) in time.The solution of the ODE is much easier and is a good approximation to the solution of the PDE.The final ODE has the form where M denotes the mass or inertia matrix, S the stiffness matrix, u the vector of the degrees of freedom and p the load vector.
For the FEM approach, we divide the tower into elements that are in our case hollow cylinders.We consider movements in lateral or z-direction, in axial or y-direction and torsional movements around the x-axis, see Figure 1 for the coordinate system.There are no movements in x-direction as the tower is supposed to be rigid against shear forces.

Establishment of the Mass and Stiffness Matrix
The mass matrix M and the stiffness matrix S are derived based on the Ritz method following exactly the way described in [5].We have divided the tower into N elements with one node at the bottom and one at the top of each element.Each node has five degrees of freedom (DOF), • v i , w i -displacement in yand z-direction, • β xi -cross section rotation (or torsion angle), and • β yi , β zi , -cross section slope in yand z-direction, c.f. Figure 2. We have to take into account the boundary and transition conditions between the elements.The global system displacement vector or vector of DOF is defined by where the index number denotes the node number.We assume that the base of the tower does not move.Therefore all DOF of the first node are zero and we have 5N entries.The system stiffness matrix and mass matrix are build from the element stiffness and element mass matrices according to the order in the global displacement vector and have the dimension 5N × 5N .The mass characteristics from rotor and nacelle have to be added in the mass matrix at the corresponding DOF.

Establishment of the Load Vector
To complete the description of the system Equation ( 1), it remains to describe the load vector p.Within the scope of this paper, we assume the wind turbine is only affected at the top node of our FE model, and neglect all loads created, e.g., by wind, in other nodes.The forces and moments in the load vector are ordered according to the order of the DOF in the displacement vector u.Hence, the load vector has the form where F y and F z denote forces in yand z-direction, M x , M y , and M z are moments of forces or torques w.r.t. the x-, yand z-axis, respectively.Here M y = 0 because we assume mechanical energy related to M y is transformed into electrical energy.
In the following, we want to consider loads that result from • inhomogeneous mass distribution of the rotor, and • aerodynamic forces arising from incorrect pitch angles settings of the rotor blades.

Mass Imbalance
The inhomogeneity of the rotor mass distribution is modeled as point mass m located in the rotor plane at a distance r under an angle θ m w.r.t. the zero mark at blade A, see Figure 3.The forces caused by this mass are gravity and centrifugal force.With the angular frequency ω = 2πΩ, the absolute value for the centrifugal force is given by Its projections onto the zand x-axis are Here, t is time variable, and φ is the angle between blade A and the x-axis.Because the rotational plane has a distance L to the tower, the forces (F c ) z and (F c ) x also produce moments around the xand the z-axes with respect to the tower Besides, the gravity force of the point mass also creates a small moment around the z-axis.However, this moment is neglected.

Aerodynamic Imbalance
To describe the aerodynamic loads on a wind turbine, we have employed the well known Blade Element Momentum (BEM) theory, c.f. [7,10].In the BEM theory, we divide the blades into a finite number of elements which are sections of the blades into annulus segments with the center at the root of the blades.The cross section of each element is called "airfoil", see Figure 4.If we are given the airfoil data, the angle of attack of the wind, and the relative wind velocity, as well as a lift and drag coefficient table, we can calculate the normal or thrust forces F and the tangential forces T according to the BEM method, see, e.g., [10], and Figures 5 and 6.
The local pitch angle θ for each blade element, i.e., the angle between chord and the plane of rotation, is the sum of the adjusted pitch angle θ p at the blades root and the twist of the blade β: The result of this procedure are forces on each of the blades distributed over all blade elements.However, for manipulation it is more convenient to change the distributed forces on a blade into one equivalent force, see Figure 7.
The distributed forces and the corresponding lengths are calculated by Then, the aerodynamic forces can be modeled as three different thrust forces (normal to the rotor plane) and three different tangential forces as shown in Figures 5 and 6.Observing Figure 5 we see that torques caused by F 1 , F 2 and F 3 are given by where we denote by M 2 x and M 2 z the torques around the xand the z-axis on the rotor of the wind turbine and by ϕ = 2π 3 (120 • ) the angle between the rotor blades.Note that if all blades have the same pitch angle, we have F 1 = F 2 = F 3 and l 1 = l 2 = l 3 .This means that the torques M 2 x and M 2 z vanish.In addition, it is easy to see that the sum of F 1 , F 2 and F 3 is the force applying to the rotor in y direction Observing the Figure 5, we know the total force T is Projecting T onto z axis and x axis yields As mentioned before, we have a small distance L between rotor plane and tower, thus T z and T x also produces moments around the xand the z-axes: Load Vector With the Formulas (3)-( 8) the forces and moments from aerodynamic and mass imbalances in (2), summarize to:

Direct Problem
So far we have derived the matrices M, S, and the load vector p = p(θ 1 , θ 2 , θ 3 , mr, φ m ).The load p is a function of the absolute value mr of the mass imbalances, its phase angle φ m , and the pitch angles θ 1 , θ 2 and θ 3 .These five variables determine the forces and moments of forces in (9).The solution of ( 1) is the sum of the general solution of the homogeneous equation and a particular solution of the inhomogeneous equation with right hand side p.Since p can be written as a sum of a vector p 1 that contains constant entries only, and vectors p 2 and p 3 that contain the amplitudes of the cosand sinterms, respectively, we can solve the inhomogeneous equation analytically.Using standard methods, c.f. [11], we get where Therefore, given a mass imbalance (mr, φ m ) and a set of pitch angles (θ 1 , θ 2 , θ 3 ) the corresponding vibrations u are calculated by first determining p(θ 1 , θ 2 , θ 3 , mr, φ m ), and applying Formulas (10) and (11) afterwards.The forward problem can be solved for every node in our FE-model.However, in view of the inverse problem, where we rely on measured vibration data, we have to consider that vibration sensors are only mounted at a few positions.In practice there will be an acceleration sensor in y-direction, one in z-direction in the middle of the tower top, and, in the best case, one in z-direction at the end of the tower top.From the difference of the last two the torsional movement around the x-axis can be derived.Thus only data for the DOF (v N , w N , β xN ) will be available.Therefore it is necessary to restrict our solution to these positions or DOF.We denote the restricted vibration by g = u sensor .The final forward operator A that relates the imbalances causes to the vibrations is the mapping which performs the calculation of p summarized in ( 9), the solution calculation (10) and the restriction step.Unfortunately, this operator is nonlinear due to the BEM method.This property will determine the possible solution methods for the inverse problem of reconstructing (θ 1 , θ 2 , θ 3 , mr, φ m ) from measurements of g.

Inverse Problem Description
In recent years the theory of treating (non-)linear ill-posed problems has been well established.For an overview see [4].For the following, let A be a linear or nonlinear operator between real Hilbert spaces.With given data g, we want to find a solution f of the equation The problem is called ill-posed, if the solution f does not depend continuously on the data g.In fact, if we only have noisy data with noise level δ, i.e., g δ − g ≤ δ, then (provided the inverse A −1 exists) f δ = A −1 g δ might be an arbitrarily bad approximation to a solution of (13).To obtain a stable solution, one has to use so called regularization methods.The general idea is to approximate the discontinuous inverse operator by a family of continuous operators T α .The regularization parameter α has to be chosen such that lim δ→0 α(δ, g δ ) = 0 holds.For nonlinear operators, equation ( 13) might have several solutions.Thus we choose the concept of a f -minimum-norm solution, i.e., we are looking for a solution closest to an a priory given function f .A prominent example for a regularization method is Tikhonov regularization.The regularization operator is defined by with the Tikhonov functional given by For the determination of the regularization parameter α, we can use the so called Morozov's discrepancy principle where α is chosen s.t.
holds [12,13].For linear operators, the minimizer of the Tikhonov functional can be computed by solving a linear system.In case of a nonlinear operator, optimization methods have to be used additionally to compute a minimizer of (15).A classical approach to minimize the functional J α (f ) is the use of gradient methods.The gradient of the Tikhonov functional is given by where the linear operator A (f ) * is the adjoint of the derivative of A at the point f [14].For other minimization methods we refer to [15,16].Please note that the above mentioned methods require the derivative of A, which is in our case not available at the moment.Instead we used the Matlab routine, fminsearch.fminsearch uses the simplex search method which is a direct search method.Hence, we can avoid to use numerical of analytic gradients in case that A does not have Frechet derivative.The weak point of this method is the slowness in running time.

Computational Results
We have tested the performance of the reconstruction technique for a turbine of the type VESTAS V80-2MW with 100 m tower height and artificial data.We achieved the artificial data by employing the forward operator (12) for a given imbalances situation.The exact vibration data were disturbed by an additive and multiplicative error in order to simulate the noise that arises in measurement.This required the following steps: 1. Construction of mass and stiffness matrix using the technical parameters of the V80-2MW (geometry, material properties, first (bending) eigenfrequency 0.255 Hz) 2. Setting of a 2 degree pitch angle deviation at the blade C and a mass imbalance of 500 kgm at angle φ m = 4π/3 = 4.19 rad: [θ 1 , θ 2 , θ 3 , mr, φ m ] = [0, 0, 2, 500, 4.19] 3. Building the load vector p using formulas (3) − (9) (airfoil code NACA63-421, wind speed 6 m/s, rotational speed 18 rpm, i.e., Ω = 0.3 Hz) 4. Solving equation ( 1) using (10) 5. Adding 10% noise to the data to simulate the measurements 6. Calculate an approximate solution (θ 1 , θ 2 , θ 3 , ∆mr , φ m ) by minimizing the functional (15) with an appropriate regularization parameter

Parameter Identification
For our test computations, we set the exact parameters [θ 1 , θ 2 , θ 3 , ∆mr, φ m ] = [0, 0, 2, 500, 4.1888].We have employed Tikhonov regularization for different kind of data.In a first attempt we used vibration data from all nodes of the model, leading to almost perfect reconstructions.However, in practice only the displacements in yand zdirection at the top of the tower are available.Reconstructions with these data were done in a second attempt.Those data and their noisy version are presented in Figure 8.
The data was contaminated with 10% noise.Tikhonov regularization with α = 10 −6 was used for all reconstructions.The results are displayed in Table 1.The error is split into two parts related to aerodynamical and mass imbalances:  The results in Table 1 show that the quality of the reconstruction depends on the initial value for the reconstruction.Good initial values, in particular for the mass imbalance, lead to good reconstructions.The main reason for this behavior is the fact that the Tikhonov functional with nonlinear operator might have several, even local minima, where the iteration might get stuck.Therefore, it is crucial to find a good initial value for the minimization routine.As mentioned in the Introduction, the algorithm proposed in [3], reconstructing mass imbalances only, is inaccurate if we neglect existing aerodynamic imbalances.However, it could provide us with an approximate mass imbalance that serves well as an initial guess for our new algorithm.In Table 1, the last result is calculated by using the initial value from the above mentioned algorithm.The performance of the algorithm was tested on several other examples, see Table 2.In these examples, we have used vibrations of the top of the tower as artificial data.The data were disturbed by 10% noise.Again, the initial values were computed with the algorithm from [3].

Balancing
For the balancing process, we changed the mass and aerodynamic imbalance in the first example according to the reconstruction (using the last reconstruction in Table 1).Hence, the new pitch angles are given by [0, 0, 2] − [0.37, −0.92, 1.75] = [−0.37,0.92, 0.25], and the residual mass imbalance is computed by 500e i4π/3 + 491.7e i (4.1929−π) .Thus, it is 8.5 kgm located at the phase angle −2.33 rad or 226.5 • .Note that in practice this is achieved by placing extra weights on the blades.Using these new parameters we have calculated the residual vibration which is plotted in Figure 9. Clearly, the remaining vibrations are much less after balancing.

Summary
In this paper we developed an algorithm that can reconstruct both mass imbalances and aerodynamic imbalances arising from pitch angle errors in the rotor of a wind turbine.For the first time it is possible to calculate those values from vibration measurements at the top of the tower only and at the same time.To this end we have developed a turbine model of a Vestas V80-2MW that provides us with the system mass and stiffness matrix.Furthermore, we calculated the load vector that arises in the presence of imbalances using the BEM method.We have solved the vibration equation that connects the vibrations of the system and the load explicitly.This has enabled us to handle the direct problem of computing vibrations in every node of the model for given imbalances, which is necessary to solve the inverse problem by using regularization techniques.With numerical examples we have confirmed that our method works well if the initial guess for the mass imbalance is good enough.This value can be obtained from a previously developed algorithm that works for mass imbalances only.In the future the algorithm has to be tested with data from field measurements.Its successful application will lead to a a new condition monitoring tool that will provide a continuous monitoring of imbalances and a correct distinction between mass and aerodynamic imbalances from pitch angle deviations without using expensive optometrical measurements and vibration measurements with additional test weights.

Figure 2 .
Figure 2. Displacement of a beam.

Figure 7 .
Figure 7. Force equivalent to distributed force.

Figure 8 .
Figure 8. Exact vibration and disturbed vibration in the top of the tower.

Figure 9 .
Figure 9. Vibrations in z direction before and after balancing with the computed parameters.