Skip to main content

Simultaneous estimation of mass and aerodynamic rotor imbalances for wind turbines

Abstract

The safe operation of wind turbines requires a well-balanced rotor. The balancing of the rotor requires a method to determine its imbalances. We propose an algorithm for the reconstruction of two types of imbalances, i.e., mass and aerodynamic imbalances from pitch angle deviation. The algorithm is based on the inversion of the (nonlinear) operator equation that links the imbalance distribution of the rotor to its vibrations during operation of the wind turbine. The algorithm requires a simple finite element model of the wind turbine as well as the minimization of a Tikhonov functional with a nonlinear operator. We propose the use of a gradient-based minimization routine. The approach is validated for artificial vibration data from a model of a Nordwind NTK 500 wind turbine.

1 Introduction

In the growing field of wind energy extraction, the topic of rotor imbalances of wind turbines (WT) is of vital importance for the operation, safety, and durability of the turbines. Approximately 20%-50% of WT have significant rotor imbalances. This leads to damages of important components, high repair expenses, and reduced output [1]. Rotor imbalances can be categorized into mass and aerodynamic imbalances. Mass imbalances arise from inhomogeneous mass distributions of the rotor including blades and hub. They induce centrifugal forces that depend on the square of the rotational frequency Ω rot . The consequences are harmonic vibrations with Ω rot in the rotor plane (lateral). The vibration amplitude depends on the ratio between Ω rot and the resonance frequency (first eigenfrequency) Ω 0 of the turbine. Mass imbalances can be eliminated or reduced using balancing counterweights. Aerodynamic imbalances evolve from deviations in the aerodynamic properties of the blades. These lead to different thrust and tangential forces and moments, and they depend not only on the rotational frequency but also on the wind speed. Aside from lateral vibrations, they induce axial and torsional vibrations of WT. The main reason for aerodynamic imbalance is a relative deviation of the pitch angles of the blades. This can only be eliminated by a correction of those angles. Additional masses are not helpful here [1]. In practice, both types of imbalances are frequently observed both separately and in combination.

Usually, the detection of imbalances is based on spectral analysis and order analysis methods [2]. Another technique is to monitor the power characteristic [3], whereby power mean values are observed, and deviations from faultless conditions are used for the calculation of alarm limits. Unfortunately, this requires a learning phase under faultless conditions. This also holds for other signal processing methods, cf. [2]. Another disadvantage is the fact that the amount of the imbalance, e.g. absolute value and location of a mass imbalance or the deviation of one or more pitch angles, is not computable.

The state of the art in rotor imbalance determination still requires an expert team on-site. Initially, aerodynamic imbalances must be detected using optical methods. After the elimination of the aerodynamic imbalances, the mass imbalance is determined using vibration measurements with and without reference weights.

In the last few years, the authors have developed methods which allow for imbalance determination using data that can be collected by a condition monitoring system (CMS). As a consequence, the amount and location of rotor imbalance can be computed without the expensive measurements using reference weights on-site. The methods are model-based, but they use a very simple model of the WT that can be obtained easily for each turbine type from elementary geometrical and physical parameters. The general basics of our approach are summarized in Section 2. Essentially, the method consists of the following steps: 1. modeling of the WT; 2. description of the loads (forces, moments) arising from rotor imbalance; 3. solution of the forward problem, i.e. the computation of the resulting vibrations in every node of the WT model for a given imbalance situation and 4. solution of the inverse problem where the unknown imbalance is determined from vibration data obtained in the nacelle. In our new method we began by determining the mass imbalance [4]. Later aerodynamic imbalance from pitch angle deviations was included [5].

In case of pure mass imbalances, only lateral vibrations need to be considered. Accordingly, in the model of the WT only degrees of freedom (DOF) in the lateral direction needed to be taken into account. Let us call this restricted approach the mass imbalance approach (MIA). The presence of aerodynamic imbalances in addition to mass imbalances would lead to additional lateral forces and therefore to incorrect reconstructions using MIA since the assumption in this model is that aerodynamic imbalances are not present. The model and the description of the forces were expanded in [5] in order to include the effect of pitch angle deviations as a main reason for aerodynamic imbalances. Nevertheless, MIA can serve to provide an initial guess for the solution of the problem for combined mass and aerodynamic imbalances. In the combined imbalance case, the inverse problem is nonlinear. As the problem is ill-posed, Tikhonov regularization is used. This includes the minimization of the Tikhonov functional which was done in [5] using a direct search method. Although we could obtain first test results, the routine was far to slow and not stable enough for extensive test computations. The acceleration of the minimization routines requires the computation of the Frechet derivative of the forward problem operator A. As most of the theory for Tikhonov regularization is formulated for real valued vector spaces, the forward operator A has to be formulated in a real setting.

This paper focuses on a stable reconstruction method of mass imbalance and pitch angle deviation based on the minimization of the Tikhonov functional. To this end, the nonlinear forward operator A is derived in detail and its Frechet derivative is computed which is new compared to [5]. Now the minimization of the Tikhonov functional is performed using the steepest decent method. The resulting reconstruction algorithm is stable and reasonable fast. Additionally, a suitable metric is defined to measure the reconstruction quality. We present test results obtained with the gradient-based algorithm for artificial data. The paper is organized as follows: In Section 3 we derive the forward problem operator A. We will be very brief if there are no changes to the preceding paper [5] and more detailed in the description of the new parts. In particular, the Frechet derivative of A is computed. Section 4 is concerned with the inverse problem. Regularization and minimization techniques using gradient-based methods are presented as well as the final algorithm that reconstructs the imbalances from measurement data. In Section 5 we present numerical test computations for the reconstruction of pitch angle deviations only and a combination of mass imbalances and pitch angle deviations as well as results for the reconstruction quality in the presented metric.

2 General approach

The procedure to determine imbalances based on a model of the WT can be summarized as follows: The load p(x) arising from an imbalance x induces displacements or vibrations u of the WT in different directions. Here we deal with displacements in axial direction (along the y-axis), lateral displacements (along the z-axis), and torsion around the tower center axis (x-axis), cf. Figure 1. The displacements and the imbalance loads are coupled via a Partial Differential Equation (PDE) in time and space which is seldom explicitly solvable. Therefore, we use the Finite Element Method (FEM) to transform the PDE into a system of Ordinary Differential Equations (ODE) in time. To this end we consider the rotor and nacelle as point masses and the WT tower as a flexible shaft that is divided into n elements with n+1 nodes (see Figure 1). The ODE system reads

M u ¨ (t)+Su(t)=p(t),
(1)

where M and S represent the mass and stiffness matrices, u the vector of displacements or degrees of freedom (DOF) of each node, and p the vector of loads for each DOF. There could be an additional damping term D u Ë™ (t) on the left hand side of (1), but so far we have achieved sufficiently good results with damping neglected.

Figure 1
figure 1

Model of a wind turbine and real wind turbine. Top: Finite-Element model of a WT with coordinate system. Bottom: NTK500/41 in Denmark.

The matrices M and S need to be derived for each different type of WT separately. This is the first step in our procedure. The second step pertains to the description of the load p(x) as a function of the imbalances x. In the third step, the now fully described ODE system (1) is solved. The solution operator A Ëœ maps the imbalances x to the displacement u at each node of the FE model: A Ëœ x=u. Now we restrict the displacement vector u that contains the displacement of each model DOF to the vector y of the DOF that are accessible for measurements. The operator A Ëœ is accordingly restricted to A. The resulting operator equation is Ax=y where the operator A maps the imbalance x of the rotor in terms of mass imbalance and pitch angle deviation to the displacement y in the nodes that can be observed by measurements. This operator equation represents the forward problem. Hence step 1-3 are connected with the solution of the forward problem. Now in a fourth step we solve the inverse problem, i.e., we determine the imbalance x from measurements y. Since the inverse problem is ill-posed, we use regularization techniques, [6]. This approach has been applied successfully in the more simple case of mass imbalances only [4]. First steps in the direction of reconstructing mass imbalances and pitch angle deviation simultaneously were done in [5].

3 Forward problem

3.1 Mass and stiffness matrices

The derivation of M and S for the simple model of MIA was taken from [7] and presented in [4]. In the situation of mass and aerodynamic imbalances the model needs to be expanded to displacements in axial direction and torsion around the tower central axis. Still, the derivation of M and S follows exactly the Ritz method described in [[7], p.155-160]. Once the model matrices are derived for a certain WT type they need to be optimized for each individual turbine of this type with respect to the first eigenfrequency Ω 0 . The general optimization procedure is described in [4].

Thus far we have modeled six different types of WT in this way: Vestas V80, V82, V90, Südwind S77, NTK500/41 and GE1.5. For the test computations presented in this paper, we have used the model of a Nordtank turbine NTK500/41, see Figure 1 (bottom), with 35 m hub height because we are provided with aerodynamic data for the blades of this WT.

3.2 Forces and moments from imbalances – the load vector p

The imbalance of the rotor creates forces and moments that affect the top node of the FE model. The forces and moments are ordered in accordance to the DOF in the displacement vector u:

p T (t)=(0,…,0, F y , F z , M x , M y , M z ),
(2)

where F y and F z denote forces in y- and z-direction, M x , M y , and M z are moments of forces or torques w.r.t. the x-, y- and z-axis, respectively. In particular, we have M y =0 as we assume that mechanical energy related to M y is transformed into electrical energy. The forces and moments are the sum of the forces and moments from mass imbalance and pitch angle deviation:

F y = F 1 + F 2 + F 3 , F z = T z + ( F c ) z , M x = M x ( m ) + M x ( F ) + M x ( T ) , M z = M z ( m ) + M z ( F ) + M z ( T ) .
(3)

Here, the F 1 , F 2 , F 3 are the thrust forces at each blade. If there is no deviation in the pitch angle of the blades all three forces are equal. T z is the projection of the sum of the tangential forces at each blade T= T 1 + T 2 + T 3 onto the z-axis. It is zero if there is no deviation in the pitch angles. F c denotes the centrifugal force induced by the mass imbalance. Its projection onto the z-axis ( F c ) z adds to F z . The moments M x ( m ) , M z ( m ) are induced by F c of the mass imbalance, M x ( F ) , M z ( F ) , M x ( T ) , M z ( T ) correspond to the thrust force and tangential force. Whereas the mass imbalance force can be directly computed, the aerodynamic imbalance forces need to be determined via, e.g., the Blade-Element-Momentum (BEM) method, a nonlinear procedure. In the BEM method, the blades are divided into a finite number of elements which are annulus segments with the center at the root of the blades. The cross section of each element is called an ‘airfoil’. If the airfoil data, the angle of attack of the wind, and the relative wind velocity as well as a lift and drag coefficient table are given, the normal or thrust forces F and the tangential forces T can be calculated according to the BEM method, see, e.g., [8]. The BEM method is used in a lot of programs for the simulation of aerodynamic effects, cf. [9], which were tested on real WT.

For the presentation of the forces and moments, we will use the following notations and abbreviations (cf. Figures 2 and 3):

  • mr, Ï• m – absolute value and angle w.r.t. blade A of the mass imbalance,

  • θ i

    – deviation of the pitch angle of i th blade to an adjusted pitch angle θ,

  • ω=2πΩ

    – angular frequency (note: WT rotates clockwise seen from the wind),

  • φ=2/3Ï€

    – angle between blades,

  • Ï• – initial location of blade A w.r.t. a zero mark (equals 0 if the mark can be chosen at blade A),

  • L – distance from hub to nacelle midpoint (x-axis),

  • F i :=F( x i )

    – thrust force at blade i with pitch angle deviation θ i , i=1,2,3,

  • T i :=T( x i )

    – tangential force at blade i with pitch angle deviation θ i , i=1,2,3,

  • l i :=l( x i )

    – corresponding length (distance of resulting force vector from nacelle) for pitch angle θ i , i=1,2,3.

Figure 2
figure 2

Notations and force from mass imbalance.

Figure 3
figure 3

Thrust and tangential forces. Left: Thrust forces F 1 , F 2 , F 3 and corresponding lengths. Right: Tangential forces T 1 , T 2 , T 3 .

3.2.1 Forces and moments from mass imbalance

The centrifugal force of a mass m with the distance r from the hub rotating with the frequency ω has the absolute value

| F c |= ω 2 mr.

The projections of F c onto the z- and x-axis are

( F c ) z = | F c | sin ( ω t + ϕ + ϕ m ) , ( F c ) x = | F c | cos ( ω t + ϕ + ϕ m ) ,
(4)

cf. Figure 2. The rotor has a small distance L to the tower center, the x-axis. Thus the forces induce moments given by

M x ( m ) = ( F c ) z â‹…L, M z ( m ) = ( F c ) x â‹…L.
(5)

3.2.2 Forces and moments from pitch angle deviation

The thrust force F and the tangential force T on a blade depend on the geometry and aerodynamics of the blade, the wind speed, the rotational speed, and the pitch angle. If the pitch angles are not the same for every blade the forces are also different. The derivation of the thrust forces F 1 , F 2 , F 3 , the tangential forces T 1 , T 2 , T 3 and the corresponding lever arms l 1 , l 2 , l 3 (cf. Figure 3) uses the nonlinear Blade Element Momentum (BEM) method. This is described in detail in [5, 8]. The BEM method is used to compute the forces for each blade element. For our purpose the distributed forces are converted into one equivalent force at a corresponding lever arm l i , cf. Figure 3. Also from Figure 3 we observe that the thrust forces F i act in y direction. Thus

F y = F 1 + F 2 + F 3 .
(6)

The forces create moments around the x- and the z-axis:

M x ( F ) = F 1 l 1 sin ( ω t + ϕ ) + F 2 l 2 sin ( ω t + ϕ + φ ) + F 3 l 3 sin ( ω t + ϕ + 2 φ ) , M z ( F ) = F 1 l 1 cos ( ω t + ϕ ) + F 2 l 2 cos ( ω t + ϕ + φ ) + F 3 l 3 cos ( ω t + ϕ + 2 φ ) .
(7)

The total tangential force

T= T 1 + T 2 + T 3

has components in x- and z-direction:

T z = T 1 sin ( ω t + ϕ ) + T 2 sin ( ω t + ϕ + φ ) + T 3 sin ( ω t + ϕ + 2 φ ) , T x = T 1 cos ( ω t + ϕ ) + T 2 cos ( ω t + ϕ + φ ) + T 3 cos ( ω t + ϕ + 2 φ ) .
(8)

As mentioned before we have a small distance L between rotor plane and tower center, thus T z and T x also induce moments around the x- and the z-axes:

M x ( T ) = T z â‹…L, M z ( T ) = T x â‹…L.
(9)

3.2.3 Combination and rewriting of forces and moments – the final load vector

We observe that F y is the only constant force, whereas all other forces and moments are harmonic with the rotating angular frequency ω. The system’s answer to a constant force are vibrations with the natural eigenfrequencies ω 0 , ω 1 ,… of the WT centered around a constant displacement:

u 1 (t)= [ I − cos ( M − 1 S t ) ] S − 1 p 1 (x),

where p 1 T :=(0,…,0, F y ,0,0,0,0), and the ω i are the roots of the eigenvalues of M − 1 S. For WT the only eigenfrequency close to the operating frequency ω is ω 0 . For safety reasons, the operation of the WT with ω= ω 0 is avoided. As only vibrations with rotating frequency ω will be extracted from the measurement, the influence of F y can be neglected. Therefore we can omit F y by setting it to zero.

In detail we have

F z = T 1 sin ( ω t + ϕ ) + T 2 sin ( ω t + ϕ + φ ) + T 3 sin ( ω t + ϕ + 2 φ ) + ω 2 mr sin ( ω t + ϕ + ϕ m ) , M x = F 1 l 1 sin ( ω t + ϕ ) + F 2 l 2 sin ( ω t + ϕ + φ ) + F 3 l 3 sin ( ω t + ϕ + 2 φ ) + L F z , M z = F 1 l 1 cos ( ω t + ϕ ) + F 2 l 2 cos ( ω t + ϕ + φ ) + F 3 l 3 cos ( ω t + ϕ + 2 φ ) + L [ T ( x 1 ) cos ( ω t + ϕ ) + T ( x 2 ) cos ( ω t + ϕ + φ ) + T ( x 3 ) cos ( ω t + ϕ + 2 φ ) + ω 2 mr cos ( ω t + ϕ + ϕ m ) ] .

We use the trigonometric identities to decompose F z , M x , and M z into cosine and sine parts of the argument ωt:

F z = F z c cos ( ω t ) + F z s sin ( ω t ) , M x = M x c cos ( ω t ) + M x s sin ( ω t ) , M z = M z c cos ( ω t ) + M z s sin ( ω t ) .

With the definitions

c ϕ : = cos ( ϕ ) , c 1 : = cos ( φ ) = cos ( 2 π / 3 ) = cos ( 4 π / 3 ) = cos ( 2 φ ) = − 0.5 , s ϕ : = sin ( ϕ ) , s 1 : = sin ( φ ) = sin ( 2 π / 3 ) = − sin ( 4 π / 3 ) = − sin ( 2 φ ) ≈ 0.866 ,

we get

F z c = T 1 s ϕ + T 2 [ s ϕ c 1 + c ϕ s 1 ] + T 3 [ s ϕ c 1 − c ϕ s 1 ] + ω 2 [ s ϕ x 4 + c ϕ x 5 ] , F z s = T 1 c ϕ + T 2 [ c ϕ c 1 − s ϕ s 1 ] + T 3 [ c ϕ c 1 + s ϕ s 1 ] + ω 2 [ c ϕ x 4 − s ϕ x 5 ] , M x c = F 1 l 1 s ϕ + F 2 l 2 [ s ϕ c 1 + c ϕ s 1 ] + F 3 l 3 [ s ϕ c 1 − c ϕ s 1 ] + L F z c , M x s = F 1 l 1 c ϕ + F 2 l 2 [ c ϕ c 1 − s ϕ s 1 ] + F 3 l 3 [ c ϕ c 1 + s ϕ s 1 ] + L F z s , M z c = ( F 1 l 1 + L T 1 ) c ϕ + ( F 2 l 2 + L T 2 ) [ c ϕ c 1 − s ϕ s 1 ] + ( F 3 l 3 + L T 3 ) [ c ϕ c 1 + s ϕ s 1 ] + L ω 2 [ c ϕ x 4 − s ϕ x 5 ] , M z s = − ( F 1 l 1 + L T 1 ) s ϕ − ( F 2 l 2 + L T 2 ) [ s ϕ c 1 + c ϕ s 1 ] − ( F 3 l 3 + L T 3 ) [ s ϕ c 1 − c ϕ s 1 ] − L ω 2 [ s ϕ x 4 + c ϕ x 5 ] .
(10)

According to our model and the corresponding nodes the load vector from aerodynamic and mass imbalances is given by

p ( t ) = p c cos ( ω t ) + p s sin ( ω t ) with p c : = ( 0 , … , 0 , F z c , M x c , 0 , M z c ) T , p s : = ( 0 , … , 0 , F z s , M x s , 0 , M z s ) T .
(11)

3.3 Derivation of the forward operator A

3.3.1 The imbalance vector x and the data vector y

The imbalance vector x contains the pitch angle deviations as well as the information on mass imbalances. It is given by

x=( x 1 , x 2 , x 3 , x 4 , x 5 )= ( θ 1 , θ 2 , θ 3 , mr cos ( ϕ m ) , mr sin ( ϕ m ) ) ,
(12)

where the θ i are the pitch angle deviations of each blade, and the mass imbalance is split into real part x 4 =ℜ(mrexp(i ϕ m )) and imaginary part x 5 =ℑ(mrexp(i ϕ m )). The forward operator A maps the imbalance x onto the data y. We obtain the data by measuring the amplitude c and phase shift γ of the vibration with frequency ω at each sensor position. We have at most 3 sensors located at the last (upper) node of the model. The corresponding DOF at this node have the numbers {N−4,N−3,N−2}. We can measure

  1. 1.

    the vibration c y (x)cos(ωt− γ y (x)) at DOF (N−4) (displacement in y-direction),

  2. 2.

    the vibration c z (x)cos(ωt− γ z (x)) at DOF (N−3) (displacement in z-direction),

  3. 3.

    the vibration c β x (x)cos(ωt− γ β x (x)) at DOF (N−2) (torsion around the x-axis).

For all vibrations we record the cosine and sine part of the vibration, cf. (13). Using

ccos(ωt−γ)=ccos(γ)cos(ωt)+csin(γ)sin(ωt),
(13)

we can arrange the data vector as

y= ( c y ( x ) cos ( γ y ( x ) ) c z ( x ) cos ( γ z ( x ) ) c β x ( x ) cos ( γ β x ( x ) ) c y ( x ) sin ( γ y ( x ) ) c z ( x ) sin ( γ z ( x ) ) c β x ( x ) sin ( γ β x ( x ) ) )

for all three sensors.

3.3.2 Solution of the ODE system

In order to describe the operator A that maps the vector of imbalance x onto the data y we solve (1) for the cosine and the sine part of p, cf. (11), separately. The solution u is the sum of the solutions u c and u s of the two inhomogeneous systems

M u ¨ c ( t ) + S u c ( t ) = p c ( x ) cos ( ω t ) , M u ¨ s ( t ) + S u s ( t ) = p s ( x ) sin ( ω t ) .

Although we have neglected the damping in the model of the WT, there will always be some small damping in every system. Therefore the homogeneous solution is decaying. As the vibration data are obtained over a long time period, we can assume that the homogeneous solution has already died away. With the abbreviation

B:= ( − ω 2 M + S ) − 1

the solutions are given by

u c ( t ) = B p c cos ( ω t ) , u s ( t ) = B p s sin ( ω t )

with p c and p s from (11). Since both vectors p c and p s have only entries at the positions {(N−3),(N−2),N} and u(t) can only be measured at DOF {(N−4),(N−3),(N−2)}, we restrict B, p c , p s , and u by withdrawing the relevant rows and columns (subscript r, see (15)) and get

u r (t)=( B r p c , r )cos(ωt)+( B r p s , r )sin(ωt).

It follows that A is given by

A(x)=y= ( B r p c , r B r p c , s ) ,
(14)

where

p c , r (x):= ( F z c M x c M z c ) , p s , r (x):= ( F z s M x s M z s ) ,and B r := B { [ N − 4 : N − 2 ] ; [ N − 3 , N − 2 , N ] } .
(15)

3.4 Derivative

We recall that parts of p c , r and p s , r occurring in the operator A are derived with the nonlinear BEM method. Hence A is a nonlinear operator. The Frechet derivative A ′ of A is given by the Jacobi matrix

A ′ (x)= ( B r B r ) ( ∂ p c , r ∂ x 1 ⋯ ∂ p c , r ∂ x 5 ∂ p s , r ∂ x 1 ⋯ ∂ p s , r ∂ x 5 ) ,
(16)

with

∂ p c , r ∂ x i = ( ∂ F z c ∂ x i ∂ M x c ∂ x i ∂ M z c ∂ x i ) , ∂ p s , r ∂ x i = ( ∂ F z s ∂ x i ∂ M x s ∂ x i ∂ M z s ∂ x i ) .

We denote by ∂ T i the partial derivative ∂ T i /∂ x i and by ∂ ( F l ) i the partial derivative ∂( F i ⋅ l i )/∂ x i . Starting from (10) we have, e.g.

∂ F z c ∂ x 1 = s ϕ ∂ T 1 , ∂ F z c ∂ x 2 = ( s ϕ c 1 + c ϕ s 1 ) ∂ T 2 , ∂ F z c ∂ x 3 = ( s ϕ c 1 − c ϕ s 1 ) ∂ T 3 , ∂ F z c ∂ x 4 = ω 2 s ϕ , ∂ F z c ∂ x 5 = ω 2 c ϕ , ∂ M x c ∂ x 1 = s ϕ ( ∂ ( F l ) 1 + L ∂ T 1 ) , ∂ M x c ∂ x 2 = ( s ϕ c 1 + c ϕ s 1 ) ( ∂ ( F l ) 2 + L ∂ T 2 ) , … .

It remains to compute the partial derivatives ∂ T i and ∂ ( F l ) i , i=1,2,3. Since T, F, and l are obtained via a nonlinear numerical routine we cannot compute the partial derivatives explicitly. Instead, we use the differential quotients as an approximation provided Δ is small, e.g.

∂ T j ∂ x j ( x j ) ≈ T j ( x j + Δ ) − T j ( x j ) Δ , … , ∂ ( F j l j ) ∂ x j ( x j ) ≈ ∂ F j ∂ x j ( x j ) l j + F j ∂ l j ∂ x j ( x j ) , j = 1 , 2 , 3 .

Numerical experiments show, that Δ≈0.1 provides a sufficiently accurate approximation of the partial derivatives. Thus, the matrix

∂ p = ( ∂ p c , r ∂ x 1 ⋯ ∂ p c , r ∂ x 5 ∂ p s , r ∂ x 1 ⋯ ∂ p s , r ∂ x 5 )
(17)

can be computed.

4 Regularization

The inverse problem of (14) is to find the imbalance vector x for given data y that, in case of real measurements, are contaminated by noise. We assume that the noisy data y δ fulfill the inequality ∥y− y δ ∥≤δ, i.e., a bound on the noise is known. The operator A is nonlinear and ill-posed. To obtain a stable solution, one has to use regularization methods, [6, 10]. Theoretical results for the regularization of nonlinear problems were derived in the last decade, cf. [11–16]. We have used the most prominent Tikhonov regularization.

4.1 Tikhonov functional

In Tikhonov regularization we determine an approximation x α δ of the solution x as a minimizer of the Tikhonov functional J α (x) with a rule to choose the regularization parameter Î±. We use Morozov’s discrepancy principle, i.e., we define a sequence of α k , e.g., α k + 1 =q α k , q<1, and compute x α k δ as minimizer of J α k (x) for each k and stop the iteration if for the first time

∥ A ( x α k δ ) − y δ ∥ ≤τδ,Ï„>1 fixed,

holds. In order to minimize the Tikhonov functional for each α k we will use a gradient-based method. The Tikhonov functional is given by

J α (x)= ∥ A ( x ) − y δ ∥ 2 +α ∥ x − x ¯ ∥ 2 ,

where ∥⋅∥ is the usual l 2 (n)-norm with n=5, and x ¯ is an a priori guess of the solution. In our case, x ¯ can be the vector (0,0,0, mr ¯ cos( φ ¯ m ), mr ¯ sin( φ ¯ m )) where the mass imbalance is computed using MIA. In view of the different absolute values of the entries in the vector x it is appropriate to choose a weighted norm instead of ∥ x − x ¯ ∥ 2 . According to the common applications we chose c=(1,1,1,0.01,0.01) and define

∥ x ∥ c 2 : = ∑ i = 1 5 c i x i 2 = ( c . ∗ x ) T x , 〈 x , z 〉 c : = ∑ i = 1 5 c i x i z i = ( c . ∗ x ) T z ,

where c.∗x= ( c 1 x 1 , … , c 5 x 5 ) T denotes the component-by-component multiplication. We remark that 〈 x , z 〉 c =〈c.∗x,z〉.

Hence the Tikhonov functional now reads

J α ( x ) = ∥ A ( x ) − y δ ∥ 2 + α ∥ x − x ¯ ∥ c 2 = ∥ A ( x ) − y δ ∥ 2 + α ∥ c . ∗ ( x − x ¯ ) ∥ 2 .
(18)

4.2 The gradient of Tikhonov functional

In order to minimize the functional, the gradient has to be computed. Since A is nonlinear we use the Taylor expansion of A(x+h),

A(x+h)=A(x)+ A ′ (x)h+O ( h 2 ) ,

where A ′ is the Frechet derivative of A computed in (16) and (17). We have

J α ( x + h ) = ∥ A ( x ) − y δ ∥ 2 + α ∥ x − x ¯ ∥ c 2 + 2 〈 1 . / c . ∗ { [ A ′ ( x ) ] T ( A ( x ) − y δ ) } + α ( x − x ¯ ) , h 〉 c + O ( h 2 ) = J α ( x ) + 2 〈 [ A ′ ( x ) ] T ( A ( x ) − y δ ) + α ( c . ∗ ( x − x ¯ ) ) , h 〉 + O ( h 2 ) .

It follows that

1 2 ▽ J α ( x ) = 1 . / c . ∗ { [ A ′ ( x ) ] T ( A ( x ) − y δ ) } + α ( x − x ¯ ) = [ A ′ ( x ) ] T ( A ( x ) − y δ ) + α c . ∗ ( x − x ¯ ) .

4.3 Minimization with steepest decent and algorithm

The functional (18) can be minimized for each fixed α by the steepest decent method. Starting from an initial value x 0 the minimizer is computed iteratively by updating the current iterate in the direction of steepest descent, i.e., the direction of the negative gradient:

x k + 1 = x k − β ▽ J α ( x ) = x k + β 2 { [ A ′ ( x k ) ] ∗ ( y δ − A ( x k ) ) − α c . ∗ ( x k − x ¯ ) } .

There are theoretical results on how to choose the step size β, [6, 13]. In the numerical realization, we have chosen β=1/(2 α 2 ). If α becomes very small, β is chosen fixed, e.g. β=10. The final minimization algorithm is shown in Algorithm 4.

Figure
figure 4

Algorithm 4

5 Numerical tests

5.1 Settings

5.1.1 Model data

We have tested the imbalance reconstruction algorithm for a wind turbine of the type NTK 500 with constant rotational frequency, since the aerodynamic data for this WT can be found in the literature. Namely, we use the following system parameters:

  • Masses: m r =9,030 kg (rotor), m n =15,400 kg (nacelle), m t =22,500 kg (tower).

  • Geometry: Tower height: h=33.8 m, outer diameter at tower top: d out top =1.69 m, outer diameter at tower bottom: d out bot =2.4 m, distance from rotor to nacelle midpoint L=1.5 m.

  • Eigenfrequency: Ω 0 =0.925 Hz.

We have set the zero mark at blade A. Thus blade B correspond to an angle of 240∘ and blade C corresponds to 120∘, resp. For example, assuming a pitch angle deviation of blade B of 2∘ compared to blade A and C and a mass imbalance of 50 kgm at blade A, the corresponding imbalance vector is x= ( 0 , 2 , 0 , 50 , 0 ) T ; if the imbalance is at blade C, we have x= ( 0 , 2 , 0 , 50 cos ( 2 / 3 π ) , 50 sin ( 2 / 3 π ) ) T . We obtain radial, axial and torsion vibration data y by applying the forward operator to x. To simulate real measurements, the data y is contaminated with Gaussian noise. The disturbed data y δ as well as the noise δ=∥y− y δ ∥ are used as input data for the reconstruction process.

5.1.2 Reconstruction quality

Due to the different magnitudes of mass imbalance and pitch angle deviation we have to choose an appropriate metric in order to evaluate the reconstruction quality of our method. Additionally, we observe that classes of pitch angle vectors actually represent the same pitch angle deviation, namely

( x 1 , x 2 , x 3 )∼( x 1 +s, x 2 +s, x 3 +s),s∈R.

In order to obtain accuracy results we chose from the class of the reconstructed pitch angle deviation the representative ( x α , 1 δ +s, x α , 2 δ +s, x α , 3 δ +s) with the smallest distance to the original pitch angle deviation ( x 1 , x 2 , x 3 ), i.e. we chose s s.t.

s=arg min t ∈ R ∥ ( x 1 , x 2 , x 3 ) − ( x α , 1 δ + t , x α , 2 δ + t , x α , 3 δ + t ) ∥ 2 .

As one can easily verify s is given by

s= 1 3 [ ( x 1 − x α , 1 δ ) + ( x 2 − x α , 2 δ ) + ( x 3 − x α , 3 δ ) ] .
(19)

In practice the original pitch angle deviation is not known. Here we chose the number s in a way that at least one of the x i , i=1,2,3, is zero. The residual vibrations are the same for all representatives of the class.

We split the norm of the deviation of the reconstruction x α δ and the original imbalance x into two parts. The first part is given by the pitch angle deviations,

ε pitch ( x α δ , x ) := ∥ ( x 1 , x 2 , x 3 ) − ( x α , 1 δ + s , x α , 2 δ + s , x α , 3 δ + s ) ∥ 2 ,

where s is given by (19). The second part is the difference for the mass imbalance in the usual l 2 -norm:

ε mass ( x α δ , x ) := ∥ ( x α , 4 δ , x α , 5 δ ) − ( x 4 , x 5 ) ∥ 2 .

In this way we can give a sophisticated account of the reconstruction accuracy of both imbalances types separately. The relative errors are defined by

ε pitch rel ( x α δ , x ) : = ε pitch ( x α δ , x ) ∥ ( x α , 1 δ + s , x α , 2 δ + s , x α , 3 δ + s ) ∥ 2 , ε mass rel ( x α δ , x ) : = ε mass ( x α δ , x ) ∥ ( x α , 4 δ , x α , 5 δ ) ∥ 2 ,

provided the denominator is not zero.

5.2 Aerodynamic imbalance caused by pitch angle deviation

The first test was concerned with the reconstruction of pitch angle deviation imbalances in the absence of mass imbalances. Here we present two examples.

Example 1 We have chosen a pitch angle deviation of 2∘ at blade B, i.e., the imbalance vector is given by x= ( 0 , 2 , 0 , 0 , 0 ) T . Figure 5, top picture, shows the original and the reconstructed pitch angle deviations for a relative noise of δ rel =10%. In Table 1 the reconstructed pitch angle deviation as well as the reconstruction error is given for several noise levels. Since the mass imbalance is zero, the absolute mass imbalance error is given. In Figure 6 the reconstruction error as a function of the noise level and the development of the residual ∥A( x n )− y δ ∥ during the iteration process are presented. The original vibrations in the imbalanced case and the residual vibrations after balancing with the reconstructed balancing setting are shown in Figure 7. The original vibrations are relatively small due to the fact that the turbine rotates with a fixed rotor frequency which is only about 70% of the tower eigenfrequency. The remaining vibrations after balancing are almost invisible.

Figure 5
figure 5

Comparison of reconstruction and original pitch angle deviation. Top picture: (0,2,0,0,0) with 10% noise level; bottom picture: (1,0,−0.5,0,0) with 7% noise level.

Figure 6
figure 6

Reconstruction performance. Reconstruction error versus noise level (top picture), residual during iteration (bottom picture).

Figure 7
figure 7

Imbalance vibrations and residual vibrations after balancing.

Table 1 Reconstruction results for imbalance x= ( 0 , 2 , 0 , 0 , 0 ) T .

Example 2 The reconstruction result for an imbalance x= ( 1 , 0 , − 0.5 , 0 , 0 ) T for a reconstruction with noisy data and noise level of δ rel =7% is shown in Figure 5, bottom picture. Here, the parameter s is computed as s=0.167 and the relative pitch angle reconstruction error is ε pitch rel ( x α δ ,x)=7%.

The initial approximation x 0 for both examples was chosen as x 0 = ( 0 , 0 , 0 , 0 , 0 ) T .

5.3 Combined mass and aerodynamic imbalance

As a second test we have considered a combined imbalance setting, a pitch angle deviation of 2∘ of blade B and a mass imbalance of 50 kgm at blade A. Blade A corresponds to an initial angle of 0∘, hence x= ( 0 , 2 , 0 , 50 , 0 ) T . The data were disturbed with a noise level of 7%. To test the performance regarding the initial value for the iteration process we have chosen two different initial approximations: firstly the zero vector x 0 = ( 0 , 0 , 0 , 0 , 0 ) T and secondly x 0 = ( 0 , 0 , 0 , 51.5 , − 4.8 ) T that represents zero pitch angle deviation and a mass imbalance of 51.8 kgm at an angle of −5∘. This mass imbalance value was reconstructed from the noisy data with MIA. The reconstruction results for both start vectors are presented in Figure 8 and Table 2.

Figure 8
figure 8

Pitch angle (left) and mass imbalance reconstruction for different initial values.

Table 2 Reconstruction results for imbalance x= ( 0 , 2 , 0 , 50 , 0 ) T .

We observe that the reconstruction process with zero vector as initial imbalance vector converges more slowly but results in a better pitch angle reconstruction. The reconstruction with the precomputed mass imbalance is faster and gives good results for both imbalance types. The original imbalance vibrations and residual vibrations after balancing are plotted in Figure 9 for all vibration directions. The amplitude of the vibrations is reduced to less than 10% in case of torsional and axial vibrations and even 2% fin case of radial vibration.

Figure 9
figure 9

Imbalance vibrations and residual vibrations after balancing. Imbalance ( 0 , 2 , 0 , 50 , 0 ) T , start vector ( 0 , 0 , 0 , 51.5 , − 4.8 ) T .

As mentioned in the Introduction, the problem of simultaneous reconstruction of mass imbalances and pitch angle deviation was already treated in [5]. In this first attempt, a direct search method for minimizing the nonlinear Tikhonov functional was used, which was far too slow for extensive test computations. Also convergence was not always be achieved. With our gradient-based reconstruction method the computation time was reduced significantly to a more reasonable time of about 2 minutes on 2.8 GHz Intel Core2 Duo processor. Convergence could also be observed in all tests.

6 Summary and outlook

The paper deals with the mathematical determination of two types of imbalances in a WT: aerodynamic imbalances from pitch angle deviation of the blades and mass imbalances of the rotor. Presently, both types have to be determined on-site by an expert team performing time consuming measurement procedures. The mathematical approach provides a reconstruction method that can be implemented into a condition monitoring system. Instead of the expert team, the method requires the knowledge of geometrical and physical parameters of the WT and aerodynamic airfoil data of the blades. Whereas the former usually can be obtained from inspection reports, the latter represent very sensible data and are hardly available for newer WT. For that reason, so far the method could only be tested for an older NTK500 turbine. Unfortunately, our industrial partners could not provide the necessary vibration data for this type. Thus we had to restrict all tests to artificial data. For future applicability, the WT owners and manufacturers need to be persuaded that tit in beneficial for the safe operation of the turbine to provide the sensible airfoil data in a confined way.

The problem of imbalance determination was addressed by the authors in former papers. In this paper we have presented an improved algorithm to reconstruct mass imbalance and pitch angle deviation of the rotor of a WT at the same time. Including pitch angle deviations into the reconstruction leads to a nonlinear forward operator. The regularization of the inverse problem requires the minimization of the Tikhonov functional. Straight forward direct search method proved to be to slow and not stable enough for our purpose. Hence we employed gradient based minimization methods. To that end, we have derived the Frechet derivative of the forward operator that occurs in the gradient of the Tikhonov functional. The accelerated algorithm was successfully applied to several test examples with artificial vibration data. In case of real vibration data we expect equally good reconstruction results as in the test computations, provided that the necessary aerodynamic data of the blades are available.

References

  1. Donth A, Grunwald A, Heilmann C, Melsheimer M: Improving performance of wind turbines through blade angle optimization and rotor balancing. EWEA O/M strategies 350; 2011.

  2. Caselitz P, Giebhardt J: Rotor condition monitoring for improved operational safety of offshore wind energy converters. ASME J Sol Energy Eng 2005, 127: 253–261. 10.1115/1.1850485

    Article  Google Scholar 

  3. Hameed Z, Hong YS, Cho YM, Ahn SH, Song CK: Condition monitoring and fault detection of wind turbines and related algorithms: a review. Renew Sustain Energy Rev 2009, 13: 1–39. 10.1016/j.rser.2007.05.008

    Article  Google Scholar 

  4. Ramlau R, Niebsch J: Imbalance estimation without test masses for wind turbines. ASME J Sol Energy Eng 2009.,131(1): Article ID 011010 Article ID 011010

  5. Niebsch J, Ramlau R, Nguyen TT: Mass and aerodynamic imbalance estimates of wind turbines. Energies 2010, 3: 696–710. 10.3390/en3040696

    Article  Google Scholar 

  6. Engl HW, Hanke M, Neubauer A: Regularization of Inverse Problems. Kluwer Academic, Dordrecht; 1996.

    Book  MATH  Google Scholar 

  7. Gasch R, Knothe K: Strukturdynamik 2. Springer, Berlin; 1989.

    Book  MATH  Google Scholar 

  8. Ingram G: Wind turbine blade analysis using the blade element momentum method. Durham University, Durham; 2005.

  9. Ahlström A: Aeroelastic simulations of wind turbine dynamics. PhD thesis. Royal Institute of Technology, Department of Mechanics, Stockholm; 2005. Ahlström A: Aeroelastic simulations of wind turbine dynamics. PhD thesis. Royal Institute of Technology, Department of Mechanics, Stockholm; 2005.

  10. Louis AK: Inverse und Schlecht Gestellte Probleme. Teubner, Stuttgart; 1989.

    Book  MATH  Google Scholar 

  11. Kaltenbacher B, Neubauer A, Scherzer O: Iterative Regularization Methods for Nonlinear Ill-Posed Problems. de Gruyter, Berlin; 2008.

    Book  MATH  Google Scholar 

  12. Ramlau R, Teschke G: Tikhonov replacement functionals for iteratively solving nonlinear operator equations. Inverse Probl 2005,21(5):1571–1592. 10.1088/0266-5611/21/5/005

    Article  MATH  MathSciNet  Google Scholar 

  13. Ramlau R: TIGRA – an iterative algorithm for regularizing nonlinear ill-posed problems. Inverse Probl 2003,19(2):433–4670. 10.1088/0266-5611/19/2/312

    Article  MATH  MathSciNet  Google Scholar 

  14. Ramlau R: Morozov’s discrepancy principle for Tikhonov regularization of nonlinear operators. Numer Funct Anal Optim 2002,23(1,2):147–172.

    Article  MATH  MathSciNet  Google Scholar 

  15. Scherzer O: The use of Morozov’s discrepancy principle for Tikhonov regularization for solving nonlinear ill-posed problems. Computing 1993, 51: 45–60. 10.1007/BF02243828

    Article  MATH  MathSciNet  Google Scholar 

  16. Ramlau R: On the use of fixed point iterations for the regularization of nonlinear ill-posed problems. J Inverse Ill-Posed Probl 2005,13(2):175–200.

    MATH  MathSciNet  Google Scholar 

Download references

Acknowledgements

The authors work supported by the Austrian Research Promotion Agency (FFG).

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Jenny Niebsch.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

Both authors contributed equally to the writing of this paper. All authors read and approved the final manuscript.

Authors’ original submitted files for images

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 2.0 International License (https://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Niebsch, J., Ramlau, R. Simultaneous estimation of mass and aerodynamic rotor imbalances for wind turbines. J.Math.Industry 4, 12 (2014). https://doi.org/10.1186/2190-5983-4-12

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/2190-5983-4-12

Keywords