Coupled flap and edge wise blade motion due to a quadratic wind force definition

The wind force on turbine blades, consisting of a drag and lift component, depends nonlinearly on the relative wind velocity. This relative velocity comprises mean wind speed, wind speed fluctuations and the structural response velocity. The nonlinear wind excitation couples the flap wise and edge wise response of a turbine blade. To analyze this motion coupling, an isolated blade is modelled as a continuous cantilever beam and corresponding nonlinear expressions for the drag and lift force are defined. After reduction of the model on the basis of its principal modes, the nonlinear response up to the second order is derived with the help of a Volterra series expansion and the harmonic probing technique. This technique allows for response analysis in the frequency domain, via which the combined flap and edge wise responses can easily be visualized. As a specific case, the characteristics of the NREL5 turbine blades are adopted. For both non-operating and operating conditions, blade responses in a turbulent wave field, based on a Kaimal spectrum, are determined. The second-order responses are shown to cause additional motion coupling, and moreover, are proven not to be negligible straightforwardly.


Introduction
The design of support structures for offshore wind turbines requires detailed knowledge on the stress variations due to ongoing aerodynamic and hydrodynamic forcing.The accumulated fatigue damage that results from both force processes limits the structural life-time.The stochastic character of both forces, as well as the nonlinear flow-force-structure interaction relations render the estimation of the fatigue life-time difficult.A key concept in the prediction of the stress variations during the structural life-time is the aerodynamic damping of the turbine rotor.This damping represents, whenever positive, a resisting force as a function of the structural velocity.Its magnitude depends on the operational state of the turbine and the directionality of the forces.
Many studies have been devoted to the understanding of the damping effect of the rotor, largely because the damping may turn negative, inducing aeroelastic instabilities.An early contribution by Petersen et al. [1] presented linearized models for quasi-steady aerodynamics and investigated the coupling between blade and turbine vibrations.Further development of aeroelastic models, predominantly focussing on dynamic stall, has been published by Rasmussen et al. [2] and Hansen et al. [3].A specific topic that has been studied extensively is the damping of edge wise rotor whirling modes, since models and measurements showed different damping ratios for forward and backward edge wise whirling modes [4][5][6][7][8].Hansen [8] not only addressed stall-induced vibrations for stallregulated turbines, but also classical flutter for pitch-regulated turbines.Damping ratios for global turbine modes have been determined in [9][10][11].
A common approach in analyzing the aerodynamic rotor damping starts with the definition of a nonlinear aeroelastic wind turbine model.After linearizing about steady-state operational conditions, the model is reformulated in terms of the modal state-space, coupling structural and aerodynamic degrees of freedom.Application of the Coleman transformation enables the derivation of collective non-periodic rotor modes together with the main lateral and torsional turbine modes.Damping ratios per dynamic mode are identified with the help of an eigenanalysis of the resulting model.
The Morison equation and the Kutta-Joukowski theorem relate the drag and lift force on a turbine blade nonlinearly to the effective wind flow.This effective wind flow comprises the unsteady wind velocity, which can be defined in terms of mean wind velocity and wake turbulence, and the structural response velocity, via which the aerodynamic damping can be expressed.This definition of the drag and lift force reveals the nonlinear coupling between the flap and edge wise blade excitations and the structural response velocities.By linearizing the force formulations, nonlinear coupling and crosscoupling terms are neglected.
In this paper, the nonlinear blade motion, coupled via the drag and lift force, is analyzed for an isolated rotating blade.A continuous blade model is reduced on the basis of its primary undamped modes.After adoption of the characteristics of the pitch-regulated NREL5 5 MW turbine blade characteristics [12] and the definition of a turbulence field, the dynamic response is determined in the frequency domain by means of Volterra series expansion and the harmonic probing technique.Besides the linear forcing terms, also the quadratic contributions are accounted for.The analysis in the frequency domain demonstrates the coupled response associated with the excited frequencies, while the excitation-related damping can easily be identified.

Blade model
A turbine blade is modelled as a cantilever beam of length R, rigidly fixed to the hub (see figure 1).The beam is described within a rotating frame of reference.The r axis represents the longitudinal axis of the undeformed blade, where the origin coincides with the point of rotation.2 presents a cross section of the blade, in which the x and y axes coincide with its principal axes.The global XY coordinate system is adopted to describe flap and edge wise blade motion; the X axis coincides with the plane of rotation and the Y axis is directed normal to this plane.The angle β describes the angle between the local x and the global X axis and is composed of both fixed blade twist and varying blade pitch.The blade twist generally varies along its longitudinal axis.Here, the blade twist is taken constant, implying that the angle β can be adopted to refer to blade pitch only.
Based on Meirovitch [13] and in accordance with Burton et al. [14], the blade is described as a geometrically linear Euler-Bernoulli beam.The initially uncoupled equations of motion for deformation in x and y direction -combined in the vector u -read: Both the displacement vector u and the force vector F comprise an x and a y component.The coefficient matrices are all diagonal; m describes the distributed mass, c s represents the structural damping, which is assumed to be proportional to mass and stiffness, EI consists of the bending stiffness components with respect to the principal axes, and T describes the tension, resulting from rotation of the blade.All matrices are r dependent.

Force definition: drag and lift
The blade cross-section is subjected to an air flow field W (see figure 2).The local width of the aerofoil is given by c, and α represents the angle between the flow vector and the local x axis, the socalled angle of attack.The exposure of the aerofoil to the air flow results in a drag force D and a lift force L, see figure 3. The Morison equation is adopted to express the relation between the air flow and the drag force, which in this case is restricted to viscous drag only.
where ρ the air density and C d is the Reynolds-number-dependent drag coefficient, accounting for the frictional stresses that can develop as a result of the air flow.
The lift force for attached flow can be determined from the Kutta-Joukowski theorem, which defines the force as a function of the air circulation strength Γ and the ambient air flow velocity: The air circulation strength can be defined as: where e r is the unit vector along the local r axis.For relatively small angles α, the lift coefficient C l can be expressed as: The resulting force vector F is a function of r and time t.The air flow vector W is both r and t dependent, as is automatically the angle of attack α.The chord width c varies in space only.The same applies for the drag coefficient C d if its Reynolds number dependency is neglected.

Constituents of the air flow vector
Both drag and lift are defined as a function of the air flow vector W, which is positioned with an angle α to the local x axis of the aerofoil (see figure 3).The air flow can be expressed as a summation of vectors: where W represents the mean air flow velocity and w the velocity fluctuations around the mean.The time derivative of u represents the structural response velocity, which is responsible for the added damping.
In terms of the global coordinate system, the mean air flow velocity consists of the mean wind velocity Y W and the in-plane rotational velocity calculated from the rotational speed of the rotor Ω and the r coordinate of the aerofoil cross section under consideration.The velocity fluctuation vector defines the turbulence field in both X and Y direction.Moreover, for rotating turbines, the blade experiences the mean in-plane wind velocity as harmonically varying.The variation of Y W along the rotor azimuth is thought to be expressed in w Y .
In order to incorporate these vectors in the equations of motion, equation ( 1), the vector components need to be transformed to the local reference system by rotation over the angle β.After transformation to the local reference system, the vector components are referred to as W x and W y , and follow from: )

Modal decomposition
The motion of the turbine blade is described by a system of coupled nonlinear partial differential equations, allowing for motions in both x and y direction, see equation (1).Coupling takes place via the forcing term, due to the structural response velocity.Likewise, the nonlinearity is pronounced in the forcing term, since eventually in both drag and lift forcing terms the multiplication W|W| appears.Due to the presence of the damping and the centrifugal stiffening term, the system of differential equations contains operators that are not self-adjoint, implying that no classical, i.e., real, dynamic modes with fixed nodes exist.In order to identify structural modes of vibration, a phase shift within each mode should be accounted for.Nevertheless, the system of partial differential equations is reduced to a system of ordinary differential equations by means of modal decomposition.Both u x and u y are expressed as an infinite series of generalized coordinates q and shape functions ψ, . The adopted shapes functions, or modes, correspond to those of an undamped and untensioned blade.Since these modes do not fulfil the orthogonality conditions with respect to all differential operators, the initial partial differential equations cannot be decomposed into a system of uncoupled ordinary differential equations.By restricting the decomposition to the first modes for both x and y directions, a system of two ordinary differential equations can be obtained, only coupled via the structural response velocity.After multiplication of both differential equations by (1)   x ψ and (1)   y ψ , respectively, and integration over r, the following ordinary differential equations remain: , and .

Force decomposition
In order to take a closer look at the right hand sides of the sets of differential equations, first the scalar equations of the forcing vector F are presented, using equation ( 2), ( 3) and ( 4): It can immediately be recognized that the decomposition concerns the multiplications W .This expression reveals the time dependency of the angle of attack and introduces an additional coupling between the motions in x and y direction.It should be noted that blade torsion, mainly resulting from the twisted shape of the blade, may bring an important contribution to the time variation of the angle of attack α.Since the twist angle is neglected in the current analysis, no α variation due to torsion takes place.The drag coefficient C d is assumed to be α independent.This assumption is valid for relatively small values of α, a condition that was already adopted for the lift coefficient C l .
In accordance with equation ( 6), the force expressions ( 7) and ( 8) need to be multiplied with 1 x ψ and 1 y ψ , respectively, and integrated with respect to r over the interval 0 ≤ r ≤ R. It should be noted that in general the chord width c and the drag coefficient C d are r dependent.The same applies for the tangential velocity of the aerofoil X W , which is calculated from the rotational speed of the rotor Ω and r.
Both F x and F y depend on the wind fluctuations in x and y direction.Moreover, both expressions contain structural response terms related to motion in both directions.A closer inspection of these equations can be facilitated by the application of Taylor series expansion.If only the linear terms are considered, and both the drag coefficient and the in-plane wind velocity are set to zero, the following expressions remain: ( ) (1) (1) (1) cos 2 sin 2 cos 2 sin 2 Here, the drag coefficient has been set to zero, because in the linear regime its value is small compared to the lift contribution.By setting the in-plane wind velocity to zero, the rotor is assumed to be perfectly yawed towards the direction of the mean wind velocity.For convenience the wind velocity terms are expressed in terms of the global reference system.The linear forcing terms provide some valuable insight in the coupling of the lateral blade motions.First of all, the forcing in x direction is amplified by the dynamic response in y direction.For a rotating rotor, when the pitch angle is small, this negative added damping term is proportional to the mean wind velocity.Besides, out of plane wind velocity fluctuations generate a force which is directed oppositely to in plane wind velocity fluctuations.Equation ( 9) also shows that for a non-operating rotor, with blades feathered into the direction of the mean wind velocity, the linear lift force in x direction is zero.
The linear y forcing expression (10) consists of two fluctuation terms and two structural response terms.For a rotating blade with small pitch angle, the force due to in-plane fluctuations depends linearly on the mean wind velocity.The force resulting from out-of-plane fluctuations is proportional to the rotational velocity of the blade.For small pitch angles, the structural response terms generally represent positive added damping terms.The aerodynamic damping of a rotating blade, as presented in [15], can be recognized in the fourth forcing term, relating the structural response in y direction to the force in y direction.

Volterra series expansion
Systems containing nonlinearities of the polynomial type can be analyzed in the frequency domain with the application of the Volterra series expansion [15].By means of Volterra frequency response functions, input-output relations per polynomial order can be established.The technique is based on the expansion of the nonlinear operators into a series of homogeneous operators, consisting of multiple convolutions [16].Considering the k-th generalized coordinate, the frequency response can be thought of as an infinite series, representing the first and all higher-order contributions: Since this study only considers the primary modes, the superscript k is omitted from here on.The i-thorder contribution is defined by the corresponding Volterra kernels ( ) i j H .For instance, for the system under consideration with two degrees of freedom, the first-order of the Volterra series yields: The explicit reference to the frequency dependency is adopted for convenience sake.For the secondorder term it follows: where both the Volterra kernels and the input signals are a function of ω and ω 1 .The frequency ω 1 represents an auxiliary integration variable.An elegant method to derive the Volterra kernels is the harmonic probing technique, which is based on the idea that a harmonic input results in a harmonic output [17].For linear systems, this approach is very common.The identification of higher-order kernels require the balance of multiple harmonic input signals with the system output.Worden et al. [18] described this kernel-identification technique for multi-input multi-output systems.The main limitation with respect to the application of this method is the complexity of the nonlinear system.If the system contains nonlinearities of the polynomial type, the kernels can in general be identified.
The response to finite-order polynomials can be determined exactly with the application of the Volterra series expansion in combination with the harmonic probing technique.Polynomials of infinite-order, however, require truncation in order to suit Volterra analysis.The same applies for finite-order polynomials that consist of too many terms, since for higher-order kernels the method becomes computationally expensive.
The nonlinearities of the system under consideration are not of the polynomial type.To facilitate the Volterra expansion and the harmonic balancing, a Taylor series expansion is applied to the right hand sides of the equations of motion, equation ( 7) and (8).The input and output variables w x , w y , q x and q y are taken as expansion variables.Expansion is carried out about the zero means and to the third-order, implying that linear and quadratic components are accounted for.

NREL5 turbine
To analyze the structural response of a wind turbine blade, the blade characteristics of the academic NREL5 turbine are adopted [12].A typical aspect of these blades is the relatively high Lock number, which expresses the aerodynamic lift capability of a blade in comparison to its weight.As a result of this, other researchers -among which Bir and Jonkman [11] -have found high aerodynamic damping values, which may be not in the same range as for other turbine types.The first flap wise natural frequency, for bending around the local x axis, is approximately 4.21 rad/s, whereas the first edge wise natural frequency is approximately 6.79 rad/s.

Excitation
In reality, the unsteady wind field as experienced by a rotating turbine blade, consists of the ambient turbulence, wake turbulence and gust slicing effects.The definition of representative frequency input signals is a scientific field on its own.Nevertheless, to illustrate the frequency response to wind fluctuations, a free-field Kaimal spectral density function is adopted to describe the turbulence field: For the numerical analyses, a mean wind velocity of 11.4 m/s, a longitudinal turbulence intensity of 10% and a longitudinal turbulence length scale of 150 m are adopted.The amplitude signal is defined with frequency steps of 0.02 rad/s, within the frequency band from -10.0 rad/s to 10.0 rad/s.The phase spectra are randomly generated for phase angles from -π to π, implying that the complex frequency signal represents the Fourier transform of a random time signal.

Non-operating turbine
First, the response of a standstill blade is analyzed.The blade is assumed to be pointing upwards and feathered into the direction of the mean wind velocity, implying a pitch angle of ½π.The aeroelastic drag coefficient is set at 0.1, irrespective of the angle of attack.The lift coefficient follows from equation (5).
Figure 4 presents the linear frequency-domain responses.Due to the chosen pitch angle, the local x axis coincides with the global Y axis; the local y axis is directed oppositely with respect to the global X axis.This implies that the blade motions purely result from the wind fluctuations in the same direction.Close inspection of equation ( 9) and (10) reveals that the x component of the linear lift force for this scenario equals zero; the linear (1)   ; x Y q ɶ completely results from drag.Equation (10) indicates that in- plane wind fluctuations generates y motion due to the varying lift force.The dynamic amplification, however, is highly damped by the structural response.No linear coupling between the motions is identified.
The second-order responses, presented by figure 5 and figure 6, clearly reveal some motion coupling for the non-rotating blade.The x response results from auto-excitation in both X and Y direction, where the coupled ( 2)   ; ɶ x XX q response is of the same order of magnitude as the linear response.The responses show clear peaks at the frequency corresponding to the natural frequency of the considered motion.Whether these responses are lift of drag dominated requires a more detailed analysis.Contrarily to the linear response, the second order response in y direction also result from the combined excitations.The response amplitude, however, is an order of magnitude smaller than the linear response amplitude.To what extent the structural response damps the dynamic amplification cannot be derived instantly from these figures.The shape of the peaks, however, indicates a high amount of effective damping.

Operating turbine
In contrast to the non-operating turbine, the blade is now given a rotational speed of 12.1 rpm.This coincides with the rated rotor speed of the adopted turbine.Blade pitch is set at 0.1.All other parameters are the same as in the previous non-operating case.
The linear responses are illustrated by figure 7, which shows the combined flap and edge wise response amplitudes, responses in local y and x direction respectively, to both in-plane (X) and out-ofplane (Y) excitations.The main, however, highly damped motion is observed in flap wise direction, resulting from wind fluctuations in the global Y direction.From this motion, commonly, the fore-aft aerodynamic damping of wind turbines is derived.Besides, the excitation in global Y direction induces edge wise motion, with a peak at the corresponding natural frequency.The large (1)   ; y Y q ɶ response, compared to the (1)   ; y X q ɶ response, results from the relatively high rotational velocity of the blade with respect to the mean wind velocity (see equation (10)).According to equation ( 9), the edge wise forcing is affected by the flap wise response.Due to the absence of a pronounced peak in the flap wise response, this interaction cannot be deduced from figure 7.  Second-order y response (operating turbine).The second-order responses, depicted in figure 8 and figure 9, show most contribution in edge wise x direction, resulting from global out-of-plane excitation.This motion contributes considerably to the peak that was already observed for the linear response, see figure 7.All other response contributions are at least one order of magnitude smaller than the linear responses.For the specific case under consideration it seems reasonable to neglect these second-order effects.Noteworthy, however, is the clear peak in the (2)   ; ɶ y YY q response, which is observed at the natural frequency of the edge wise motion.

Conclusions
This paper describes the derivation of a reduced model of an isolated wind turbine blade.The coupled nonlinear wind excitation, including structural motion was elaborated and the structural response was analyzed in the frequency domain.To this end, use was made of a Volterra series expansion and the harmonic probing technique.The estimation of the flap and edge wise response due to combined nonlinear wind excitation was the main purpose of this work.
Dynamic responses were determined for both a non-operating and an operating blade.Use was made of the blade characteristics of the NREL5 turbine.The nonlinear force formulations were expanded into a Taylor series and truncated after the quadratic terms.The linear responses clearly visualized the combined flap and edge wise motion of the blade.Second-order responses were proven to be present and in some cases shown to be of the same order of magnitude as the linear responses.
Some additional remarks with respect to the estimation of second-order effects need to be made.First of all, no linear relation between the input and output amplitudes exists and second, the response frequencies cannot one-to-one be related to the frequency content of the excitation.This implies that a turbine blade may be excited at a natural frequency, while the excitation does not contain much energy at this specific frequency.From this it follows that conclusions with respect to second-order contributions cannot be generalized on the basis of a single analysis.In this specific case, the high aerodynamic damping of the NREL5 blade may have concealed second-order effects that can be of concern for other blade types.
Unavoidably, the presented approach carries some limitations.First, the geometrically linear blade formulation does not allow for blade twist along the longitudinal axis.Its validity for larger deformations is limited too.Second, the modal reduction is restricted to the first modes only.Since the adopted modes are not adjusted for secondary effects as centrifugal stiffening and added damping, their validity may be poor.Especially for the highly damped NREL5 blade the extent to which the assumed modes are uncoupled should be doubted.Next, torsional motion is not incorporated in the analysis, implying that variations in angle of attack are not fully accounted for.Furthermore, the aeroelastic relations, which are functions of this angle of attack are assumed to be linear, or almost linear.These relations are only valid for attached flows, implying that the angle of attack should be relatively small.With respect to the aerodynamics, the adopted wind turbulence field cannot be assumed to be representative for a rotating rotor, since rotational effects and wake effects are not incorporated.Last of all, due to the coupled flap and edge wise excitation, also cubic and even higher nonlinearities appear in the forcing formulation.These terms have been neglected.

Figure 1 .
Figure 1.Blade model.Figure 2. Blade section exposed to wind velocity vector.

Figure 2 .
Figure 1.Blade model.Figure 2. Blade section exposed to wind velocity vector.

Figure 3 .
Figure 3. Lift and dragdefinitions.Figure2presents a cross section of the blade, in which the x and y axes coincide with its principal axes.The global XY coordinate system is adopted to describe flap and edge wise blade motion; the X axis coincides with the plane of rotation and the Y axis is directed normal to this plane.The angle β describes the angle between the local x and the global X axis and is composed of both fixed blade twist and varying blade pitch.The blade twist generally varies along its longitudinal axis.Here, the blade twist is taken constant, implying that the angle β can be adopted to refer to blade pitch only.Based on Meirovitch[13] and in accordance with Burton et al.[14], the blade is described as a geometrically linear Euler-Bernoulli beam.The initially uncoupled equations of motion for deformation in x and y direction -combined in the vector u -read:

Figure
Figure 3. Lift and dragdefinitions.Figure2presents a cross section of the blade, in which the x and y axes coincide with its principal axes.The global XY coordinate system is adopted to describe flap and edge wise blade motion; the X axis coincides with the plane of rotation and the Y axis is directed normal to this plane.The angle β describes the angle between the local x and the global X axis and is composed of both fixed blade twist and varying blade pitch.The blade twist generally varies along its longitudinal axis.Here, the blade twist is taken constant, implying that the angle β can be adopted to refer to blade pitch only.Based on Meirovitch[13] and in accordance with Burton et al.[14], the blade is described as a geometrically linear Euler-Bernoulli beam.The initially uncoupled equations of motion for deformation in x and y direction -combined in the vector u -read: the aeroelastic coefficients C d and C l are α dependent.Equation (5) introduced a simple relation between C l and α, which is valid for attached flows.With reference to figure 2, the angle of attack can be expressed as sin y W α = the autospectral density function for the wind fluctuation w Y .Y W is the corresponding mean wind velocity, 1 Y w L the turbulence length scale and σ Y w the standard deviation of the wind fluctuations.The wind fluctuation in X direction, which is either lateral or vertical, has been defined for a turbulence intensity of 0

Figure 9 .
Figure 9. Second-order y response (operating turbine).The second-order responses, depicted in figure8and figure9, show most contribution in edge wise x direction, resulting from global out-of-plane excitation.This motion contributes considerably to the peak that was already observed for the linear response, see figure7.All other response contributions are at least one order of magnitude smaller than the linear responses.For the specific case under consideration it seems reasonable to neglect these second-order effects.Noteworthy, however, is the clear peak in the(2)