- Email: [email protected]
Enhancement of flutter stability in wind turbines with a new type of passive damper of torsional rotation of blades
Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 171–179
Contents lists available at ScienceDirect
Journal of Wind Engineering & Industrial Aerodynamics journal homepage: www.elsevier.com/locate/jweia
Enhancement of flutter stability in wind turbines with a new type of passive damper of torsional rotation of blades Bei Chen a, Zili Zhang b, Xugang Hua a, *, Søren R.K. Nielsen c, Biswajit Basu d a
Key Laboratory for Wind and Bridge Engineering, Hunan University, 410082 Changsha, China Department of Engineering, Aarhus University, 8000 Aarhus, Denmark Department of Civil Engineering, Aalborg University, 9220 Aalborg, Denmark d Department of Civil, Structural and Environmental Engineering, School of Engineering, Trinity College Dublin, Dublin 2, Ireland b c
A R T I C L E I N F O
A B S T R A C T
Keywords: Horizontal axis wind turbines Flutter analysis Viscous dampers Vibration control
Classical flutter of wind turbine blades is defined as an aeroelastic instability where a torsional blade mode couples to a flapwise bending mode to result in a mutual rapid growth of the amplitude of the flapwise and torsional motions beyond certain rotational speed of the rotor. Till date, flutter instability of wind turbine rotor has not been considered to be a serious problem. However, with the advent of larger turbines fitted with relatively softer blades, classical flutter may become a more important design consideration. As an economically viable alternative to increasing the torsional rigidity (which is more expensive), the use of a type of viscous damper mounted inside the wind turbine blades is proposed to promote the flutter critical rotational speed of the rotor. The flutter critical rotational speed is determined by using a time-domain aeroelastic simulation technique based on the quasi-steady aerodynamic forces and also considering turbulence intensity of incoming wind. The optimal damper location and the optimal damping constant values are investigated for the DTU 10 MW wind turbine. It is shown that at an optimal tuning of the damper, the flutter critical rotational speed is only marginally dependent on the turbulence intensity of the incoming wind field, on the mean wind velocity and on the structural damping of the flutter torsional eigenmode. The novel tuned damper device proposed in this paper may increase the flutter critical rotational speed with more than 100%.
1. Introduction It is known that three typical aeroelastic instabilities may occur in modern commercial wind turbines: divergence, stall-induced vibration and classical flutter. Divergence is a static phenomenon wherein the resulting wind velocity (consist of both the ambient wind speed and the rotational speed of the rotor) becomes large enough so that the load produced for an incremental angle of attack change due to blade twisting is greater than the reaction load produced by the elastic restoring forces for the same amount of twist (Lobitz and Veers, 1998). Stall-induced vibration is a single mode phenomenon characterized by predominantly flapwise blade oscillations, and normally occurs at high angle of attack causing separation of the boundary layer on the backside of the profile. The stall induces negative aerodynamic damping, related to a l negative gradient ∂c ∂α of the lift coefficient with the angle of attack (Hansen, 2007). The separation point may be stationary, which is referred to as static stall. At the core of dynamic stall the separation point
moves on the backside of the profile, the boundary layer may even vary between fully attached or fully detached during part of the oscillation (Larsen et al., 2007). Pitch-regulated variable speed wind turbines normally do not operate in stall and the risk of stall-induced vibration is not as serious as for stall-regulated wind turbines. Flutter is a well-known destructive instability from aircraft industry but historically has not been an issue in modern commercial wind turbines’ design. Commensurately, among the utility-scale wind turbines that have been built, rarely has classical flutter ever been observed. In short, classical flutter of a beam-like structure such as wind turbine blades is explained as an unfavourable aerodynamic coupling between flapwise and torsional vibrations of blades in the presence of incident aerodynamic forces, and appears in a rapid growth of the amplitudes of the flapwise and torsional motions, which may leads to major structural failure of the blades. In a steady homogeneous wind field the onset of flutter takes place, if the average power supplied to the structure by the self-induced loads exceeds the average power dissipated by the structural
* Corresponding author. E-mail addresses: [email protected] (B. Chen), [email protected] (Z. Zhang), [email protected] (X. Hua), [email protected] (S.R.K. Nielsen), [email protected] (B. Basu). https://doi.org/10.1016/j.jweia.2017.12.011 Received 22 August 2016; Received in revised form 7 December 2017; Accepted 7 December 2017 0167-6105/© 2017 Elsevier Ltd. All rights reserved.
B. Chen et al.
Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 171–179
earlier work on a 20 KW HAWT with stiff, aeroelastically tailored 5 m blades, Lobitz and Veers (1998) predicted flutter critical rotational speed that were six times the nominal rotor speed (ΩΩcr ¼ 6), rendering flutter a moot issue. The predicted flutter limits of a SNL 9-meter CX-100 experimental blade were shown to be four times higher than the nominal rotor speed (ΩΩcr ¼ 5) (Resor and Paquette, 2011). Study on 1.5 MW-sized WindPACT wind turbine with 34 m blades has a safety margin of 2–2.5 (ΩΩcr ¼ 2:0 2:5) (Lobitz, 2004). Similar result has also been found by Hansen for a 5 MW NREL wind turbine (Hansen, 2007). Further, for the 10 MW NOWITECH wind turbine with 89 m blades, the critical rotational speed for flutter is found to be about 1.6 times the nominal rotor speed (ΩΩcr ¼ 1:6) (Vatne, 2011). However, result from a 13.2 MW wind turbine with 100 m blades in length showed that there is little or no margin on flutter speed (ΩΩcr ¼ 1:0 1:1) reported by Griffith and Resor (Griffith and Ashwill; Resor et al., 2012). The trend shows that the ratio of critical rotational speed for flutter to the nominal rotor speed ΩΩcr drops significantly as the blades grows in length from 5 m to 100 m. Therefore, it is believed that with a likely reducing blade torsional stiffness with increasing blade radius (considering current design trends) the critical rotational speed for flutter Ωcr will be approach the nominal rotational speed Ω of the rotor (Politakis et al., 2008). This also means that aeroelastic instabilities, especially flutter, would then become one of the principal design drivers (Hansen, 2007). Great efforts have been devoted to mitigating loads and improving the aeroelastic stability of wind turbines (Belamadi et al., 2016; Bottasso et al., 2016; Moshfeghi et al., 2017; Bernhammer et al., 2016). Due to the limitation of space in wind turbine blades, few works have been carried out on suppressing flutter of wind turbine blade, especially using mechanical countermeasures, such as tuned mass damper (TMD). Hence, the adaptive blades that involving the use of aeroelastic tailoring (also known as bend-twist coupling (BTC)), wherein the blade twists as it bends under the action of aerodynamic loads to shed wind-induced loads and increase the aeroelastic stability of wind turbine has been investigated by a lot of researchers (De Goeij et al., 1999; Veers et al., 1998; Rafiee et al., 2016). Even though the BTC proved to be effective in mitigating the fatigue loads for large-scale wind turbine blades, it may increases the blades proclivity for flutter at the same time (Bottasso et al., 2013; Hayat and Ha, 2015; Hayat et al., 2016). Lobitz's (Lobitz and Veers, 1998) research shown that the flutter speed of the rotor can be increased by around 11% at the optimal values of the coupling coefficient for 20 kW HAWT, since the blade is less stable as the extreme values of the coupling coefficient are approached due to the divergency limit. Further, Hayat (Hayat et al., 2016) extended Lobitz's study using a beam element model to simulate the composite laminate material structure. Results shown the flutter speed of the rotor can be increased by about 7:6 9:5% using lighter and stiffer carbon fibers. As mentioned above, there have been several researches reported on increasing flutter limits using BTC. However, their effectiveness is very limited. Inspired from this, a new ’electromagnetic’ type of linear viscous damping device is proposed in this paper in order to mitigate the torsional and flapwise vibration of the blades, hence increasing the flutter critical rotational speed of the rotor. Firstly, a wind turbine model is presented, which considers the coupling between aerodynamic, elastic and inertial loads of the blades. The parameters of the model have been calibrated to the DTU 10 MW wind turbine (Bak et al., 2013; Bak et al.,). The elastic deformations of the blades are coupled via the tower motion, which is modeled by a reduced five degree-of-freedom model, and via the drive train, which is modeled by a two degree-of-freedom model. The pitch of the blades are controlled by a collective PI controller with feedback from the rotor speed. Furthermore, a linear viscous damper model is provided. The optimal parameters of the viscous damper are calculated. Finally, the sensitivity study of this damping system to the turbulence, the mean wind velocity, the structural damping and the parameters of this damper is performed.
damping and possibly attached active or passive damping mechanism mounted on the structure (Politakis et al., 2008). Based on the studies of Bisplinghoff, Fung and Hansen (Hansen, 2004a, 2007; Bisplinghoff et al., 2013; Fung, 2004), the risk of classical flutter of wind turbine blades with attached boundary layer may rise if the following three main criteria are satisfied: the frequency ratio between a flapwise bending mode and a torsional mode is close to one, the rotational speed of the rotor is sufficiently high and the position of the shear center of the airfoil is placed between the center of mass and aerodynamic center. Considering the first criterion, the risk of flutter increases as the ratio between the first torsional mode and one of the lower flapwise modes approaches one, such that they can couple in a flutter mode. As a reference point, for the 1.5 MW Horizontal Axis Wind Turbine (HAWT) with 35 m blades, it was the second flapwise bending mode and the first torsional blade mode which are coupled in flutter (Lobitz, 2004). For the 5 MW National Renewable Energy Laboratory (NREL) turbine, used as a benchmark in numerous analyses of wind turbines, the flutter instability took place between the third flapwise bending mode and the first torsional blade mode (Hansen, 2007). Similar result has also been found by Bir for a 5 MW NREL wind turbine (Bir and Jonkman, 2007). For the second criterion, the rotational speed of the rotor Ω must be sufficiently high, resulting in higher energy in the aerodynamic forces which ’feed’ the flutter mode. Finally, the onset of flutter requires that the shear center is placed between the aerodynamic center and the center of mass. The inertial load is acting at the center of mass, whereas the aerodynamic load is acting at the aerodynamic center. Since the aerodynamic and the inertial loads approximately are in counter phase, they provide a maximum of elastic torsion with the indicated position of the shear center. Lobitz's work shown the flutter speed limit decreases when the center of mass is moved towards the trailing edge of the blade. Furthermore, when the center of mass is ahead of the shear center of the blade profile, flutter is unlikely at any rotational speed of the rotor (Lobitz, 2005). Since the flutter mode is not known beforehand, thus the flutter can only be predicted by aeroelastic eigenvalue analysis or time domain simulation. Several studies have been carried out on wind turbine stability analysis, especially for flutter. Hansen (2003; 2004a; 2004b) proposed a multi-blade coordinate transformation (MBC, (Bir, 2008)) based eigenvalue analysis. Firstly the time-independent aeroelastic motion equations which couple the structural motion with the dynamics of the aerodynamic forces were derived using Lagrange's equations and MBC transformation. Then the eigenvalue problem which is ill-conditioned based on these autonomous equations was solved by using a modal expansion of the structural DOFs based on the undamped modes of the turbine at standstill. This enabled the computation of the natural frequencies, logarithmic decrements, and mode shapes of the aeroelastic turbine modes. The lowest nominal rotational speed of the rotor which corresponding to a negative damping coefficient exhibits in eigenmode is designated as the flutter critical rotational speed. Similar method has been used by Lobitz (2004). Alternatively, the aeroelastic stability may be checked by a pure mathematical approach of using Floquet theory on the non-transformed linearized equations (Nayfeh and Mook, 1995). However, this approach should be avoided since it doesn't have any physical meaning (Hansen, 2003). Furthermore, the stability analysis may also be performed by observing a certain envelope of the flapwise or torsional motions, indicating the evolution with time of the amplitude of the involved response quantities. Flutter is defined to take place, where the envelope of the torsional angle exceeds a certain critical value θmax . This approach is useful if the instability is dominated by nonlinear mechanisms, either structural or aeroelastic (Chen et al., 2017a). According to previous studies, the critical rotational speed for flutter turns out to be higher than the nominal rotor speed. Smaller wind turbines do not experience flutter due to their high torsional angular frequency in comparison to the angular eigen frequency in the critical flapwise bending mode of the blades (Lobitz, 2004; Hansen, 2008). In
172
B. Chen et al.
Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 171–179
2. Modelling of flutter induced vibrations of wind turbines
Generator shaft kg
2.1. Wind turbine model
q7 (t)
The motion of the tower is described in a fixed, global ðX1 ; X2 ; X3 Þ-coordinate system with origin O at the interface to the foundation, see Fig. 1. The displacement of each blade is described relative to a moving, local ðx1 ; x2 ; x3 Þ-coordinate system with origin O' at the center of the hub. The position of the local coordinate system attached to blade j is merely specified by the azimuthal angle Ψj ðtÞ, which is considered positive when rotating clockwise seen from an upwind position. Then, the azimuthal angles Ψj ðtÞ may be given as:
2π 2π Ψj ðtÞ ¼ Ωt þ q6 ðtÞ þ ðj 1Þ≃Ωt þ ðj 1Þ; 3 3
j ¼ 1; 2; 3
Mr (t)
Rotor of generator
q6 (t)
Jr
Mg (t)
Mg (t) kr
Stator of generator
Rotor shaft
Fig. 2. 2-DOF model of flexible drive train with odd number of gear stages. Definition of degrees of freedom q6 ðtÞ and q7 ðtÞ.
vector q0 ðtÞ. Based on q0 ðtÞ and q_ 0 ðtÞ the kinetic and strain energy of the tower and the drive train minus the rotor may be calculated. The kinetic energy related to the rotor is accounted for at the analysis of the blades. Each blades is modeled as a Bernoulli-Euler beam with variable mass per unit length μðx3 Þ and variable bending stiffness in the flapwise and edgewise directions. Correspondingly, deformations from the shear forces are ignored. The elastic (bending) center E of all cross-sections along the blade are assumed to be placed along the x3 -axis. For each cross section a principal axis ðy1 ; y2 ; y3 Þ-coordinate is defined as shown in Fig. 3. The y3 -axis is unidirectional to the x3 -axis, and the orientation of the y2 -axis is defined by the angle δðx3 Þ, indicating the angle from the y2 -axis into the x2 -axis in the positive x3 -direction. δðx3 Þ represents the pretwist of the blade. The aerodynamic center A, the shear center S and the mass center G are also shown on Fig. 3. The position of the aerodynamic center is assumed to be fixed in time, either corresponding to the situation of attached boundary layers on the back-side of the profile, or to deep stall, where the boundary layers are totally detached. The external aerodynamic lift force per unit length pl and drag force per unit length pd acting in the aerodynamic center is transferred to the shear center in order to decouple the equations of motion for the elastic torsion and bending displacement, when formulated in the ðy1 ; y2 ; y3 Þ-coordinate system. The inertial loads, including the gravity load, are acting in the mass center G. The gravity loads are referred to the shear center for the same reason as for the aerodynamic loads. The components of velocity vector and angular velocity vector of the mass center need to be interpolated from the time derivates of the degrees-of-freedom of the element. Based on these interpolation functions the kinetic energy, as well as the mass matrix and the gyroscopic damping and stiffness matrices of the element can be calculated. The fact that the inertial loads are acting in the mass center G, and not in the elastic center E, introduces couplings in the equation of motion of importance for the flutter stability. Fig. 4 illustrates the discretization of a blade into a number n of finite elements, including the global node numbering. x3;j indicates the position of node j of the element on the x3 -axis. Also shown is the beam element j between the nodes j and j þ 1 and the definitions of the twelve degrees of freedom r1;j ; …; r12;jþ1 in the principal axis coordinate system. lj ¼ x3;jþ1 x3;j indicates the element length. The degrees-of-freedom specifying the elastic deformation of blade j relative to the hub in the
(1)
where Ψj ðtÞ indicates the azimuthal angle of blade j ¼ 1; 2; 3, Ω is the nominal rotational speed of the rotor. q6 ðtÞ is a degree of freedom specifying the deviatoric rotation of the hub on top of the nominal rotational angle Ωt. q6 ðtÞ is supposed to be numerical small quantities compared to Ωt, hence q6 ðtÞ can be ignored in Eq. (1). A five degree-of-freedom attached modal model is applied for the description of the tower motion defined by the displacements q1 ðtÞ and q2 ðtÞ of the tower top in the global X1 - and X2 -directions, and the rotations q3 ðtÞ, q4 ðtÞ, q5 ðtÞ. L indicates the length of the blade from the hub, h is the height of the tower from the base to the nacelle, and the horizontal distance from the center of the tower top to the origin O' of the moving coordinate systems is denoted s. The drive train is modeled by the degrees of freedom q6 ðtÞ and q7 ðtÞ as shown in Fig. 2, indicating the rotational angle of the rotor and generator rotor, respectively. Mr ðtÞ and Mg ðtÞ denote the workconjugated rotor torque and generator torque. Jr and Jg denote the mass moment of inertia of the rotor and the generator, and kr and kg denote the St. Venant torsional stiffness of the rotor and generator shafts. For more details about the tower and the drive train modelling, the reader is referred to Zhang (Zhang et al., 2014; Zhang, 2015). The degrees of freedom q1 ðtÞ; …; q7 ðtÞ are assembled in the column
Fig. 1. Definition of the fixed and the moving frames of reference and the degrees of freedom q1 ðtÞ; …; q5 ðtÞ. Fig. 3. Definition of principal axis coordinate system.
173
B. Chen et al.
Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 171–179
Fig. 4. Nodal numbering and definition of degrees of freedom of a beam element.
2
3 qðtÞ _ 5 zðtÞ ¼ 4 qðtÞ βðtÞ
local ðx1 ; x2 ; x3 Þ-coordinate system, are assembled in the column vector qj ðtÞ; j ¼ 1; 2; 3 of dimension 6n. Small elastic vibrations are assumed, so that linear beam theory may be assumed with due consideration of the geometrical stiffening in the x3 -direction, and the geometrical softening in transverse direction. For a more thorough elaboration, including the mathematical derivation of the displacement and velocity vector of the blades, the reader is referred to Zhang et al. (2014), Zhang (2015). The degrees of freedom of the system are assemble in the column vectors qðtÞ of dimension 7 þ 18n:
The parameter βðtÞ indicates a full-span collective pitch angle with a time delay determined by the first order filter:
_ ¼ 1 ðβ0 ðq6 ; q_6 Þ βðtÞÞ βðtÞ τ
(2)
2.2. Aeroelastic load model
Based on the applied discretization of the blades the kinetic energy _ and the potential energy UðqÞ of the tower, drive train and blades Tðq; qÞ _ may be calculated as a function of q and qðtÞ. Then, the equations of motion of the tower, the drive train and the blades follow from the EulerLagrange stationarity condition of analytical dynamics (Pars, 1979):
_ _ d ∂Tðq; qÞ ∂Tðq; qÞ ∂UðqÞ _ tÞ þ ¼ fðq; q; dt ∂q_ ∂q ∂q
(5)
βðtÞ is considered positive, when rotating in the positive x3 -direction, see Fig. 3. β0 ðq6 ; q_6 Þ indicates the pitch demand, which is modeled as PI controller with feed-back from q6 ðtÞ and q_6 ðtÞ (Ogata, 2010). The parameter τ specifies the reaction time of the pitch actuators.
2
3 q0 ðtÞ 6 q1 ðtÞ 7 7 qðtÞ ¼ 6 4 q2 ðtÞ 5 q3 ðtÞ
(4)
In accordance with (Commission et al, ), the turbulence model is constructed based on Taylors hypothesis of frozen turbulence, corresponding to a frozen turbulence field that is convected into the rotor in the global X1 direction with a mean velocity V0 . The frozen field is
€ þ CðtÞqðtÞ _ þ KðtÞqðtÞ ¼ fðz; tÞ ⇒ MðtÞqðtÞ
(3)
assumed to be a zero mean homogeneous and isotropic stochastic field, with a spatial covariance structure as elaborated in (Batchelor, 1953). Calibrated from the theoretical covariance structure of the turbulence, the first order AR model (Krenk et al., 2012) performs a first-order filtering of the white noise input, resulting in continuous, non-differentiable sample curves of the turbulence field at the rotor plane in the fixed frame of reference (Zhang et al., 2014). The present turbulence model is in good agreement with one of TurbSim's turbulence models (i.e. the Von Karman's turbulence model) (Jonkman, 2009; Krenk et al., 2012), but is highly efficient numerically and easy to incorporate into the 907 DOF model. The model generates correlated time histories corresponding to such frozen turbulence field at the points of the mesh shown in Fig. 5 consisting of 8 radial lines with 6 points in each. As seen, the fixed frame components of the convected turbulence are generated at nn ¼ na ⋅nr þ 1 nodal points at each time step t ¼ 0; △t; 2△t; …, where na denotes the
Because the elastic deformations have been assumed to be linear, the mass matrix MðtÞ, the damping matrix CðtÞ and the stiffness matrix KðtÞ merely depend q0 ðtÞ and q_ 0 ðtÞ, and more specific on q6 ðtÞ and q_6 ðtÞ. As can be seen in Eq. (1), the q6 ðtÞ and q_6 ðtÞ are supposed to be numerically small quantities compared to Ωt and Ω, and can be ignored. Hence, the mass matrix MðtÞ, the damping matrix CðtÞ and the stiffness matrix KðtÞ can be expressed as a function of t. Furthermore, the mass matrix MðtÞ and the damping matrix CðtÞ are symmetric and periodic, but due to the gyroscopic stiffness contributions the stiffness matrix KðtÞ is asymmetric periodic matrix. fðz; tÞ on the right hand side of Eq. (3) indicates the load vector conjugated to qðtÞ from self-induced aeroelastic load, pitch control load and wind load. zðtÞ indicates the state vector of the system defined as:
174
B. Chen et al.
Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 171–179
is quasi-steady, i.e. the unsteady effects are not pronounced; if k > 0:05, the unsteady effects cannot be neglected. Based on our calculation, the reduced frequency k at the outer part of the blade is much smaller than 0.05, as shown in Table 1. Looking at the values of the reduced frequencies in Table 1, it is possible to see that the reduced frequencies k in the range 0.001–0.03. Actually, Lobitz's research on a MW-sized wind turbine blade demonstrated that the computed flutter limit of a blade using quasi-steady theory is conservative than that estimated using unsteady (Theodorsen) theory. Furthermore, when flutter speed obtained from quasi-steady theory is similar to the nominal rotational speed of the rotor, the results are no more reliable and may adversely affect blade design (Lobitz, 2004). In other words, for those wind turbines fitting with relatively softer blades (e.g. 13 MW HAWT, 20 MW HAWT or HAWTs with bigger rated capacity), the quasi-steady theory may be unapplicable. In the present study, we haven't considered the unsteady aerodynamics that can be derived from Theodorsen's theory (Lobitz, 2004; Theodorsen and Mutchler, 1935; Fung, 2002; Dowell et al., 1989), because the main goal of current research is to investigate the effectiveness of the linear viscous damping device on mitigating flutter. However, the unsteady aerodynamics will be considered in our future studies. Assuming locally 2D-flow around the profile the lift and drag load per unit length pl ðx3 ; tÞ and pd ðx3 ; tÞ at the position x3 at the time t become:
9 1 > pl ðx3 ; tÞ ¼ cl ðαÞ ρV 2 ðx3 ; tÞcðx3 Þ > = 2 > 1 > pd ðx3 ; tÞ ¼ cd ðαÞ ρV 2 ðx3 ; tÞcðx3 Þ ; 2
Fig. 5. Nodal numbering and definition of degrees of freedom of a beam element (Zhang et al., 2014).
number of radial lines from the center and nr denotes the number of equidistantly placed nodes along a given radial line. Then, the fixed frame components of the rotational sampled turbulence vector on each blade with the azimuthal angle Ψj are calculated through the linear interpolation of the turbulence of the adjacent radial lines according to the position of the blade at each time step. Finally, the moving frame components of the rotational sampled turbulence can be determined through the following coordinate transformation equation:
3 2 1 0 v1;j ðx3 ; tÞ 4 v2;j ðx3 ; tÞ 5 ¼ 4 0 cosΨj 0 sinΨj v3;j ðx3 ; tÞ 2
3 32 0 v1;j ðX; tÞ sinΨj 54 v2;j ðX; tÞ 5 cosΨj v3;j ðX; tÞ
where cd ðαÞ and cl ðαÞ indicate the lift and drag coefficients obtained from static 2D wind tunnel test data, αðx3 ; tÞ is the instantaneous angle of attack, ρ is the mass density of air, Vðx3 ; tÞ is the resulting wind velocity and cðx3 Þ is the chord length of the profile at the coordinate x3 . The mean wind velocity V0 is directed in the global X1 -direction, and is assumed to be constant over the rotor area. Further, the nominal rotational speed of the rotor Ω is assumed to be constant in time. Then, the resulting wind velocity on the profile becomes:
(6) Vðx3 ; tÞ ¼
ence with X ¼ ½0; x3 sinΨj ; x3 cosΨj T . Due to the longitudinal correlation of the incoming turbulence, a certain periodicity is present as spectral peaks at Ω; 2Ω; 3Ω… in the frequency domain representation of the rotational sampled turbulence. This simple AR model used here does not represent the low-frequency, large-scale turbulent structures very well, given the homogeneity and isotropy assumption. Furthermore, the dynamics of the tower are more related to the frequency component of turbulence in the vicinity of the tower frequency. In this paper, the rotational sampled effect seems to be more important and is well accounted for by the present turbulence model (Zhang et al., 2014). Here the quasi-steady aerodynamic is employed which means the boundary layers on the profile are assumed to be fully attached and a change of the angle of attack is causing a comparable change of the aerodynamic loads without time delay. According to the criteria proposed by Leishman (2006), if the reduced frequency k < 0:05 the problem
V1 ðx3 ; tÞ ¼ ð1 aÞV0 þ v1 ðx3 ; tÞ u_1 ðx3 ; tÞ V2 ðx3 ; tÞ ¼ ð1 þ a0 ÞΩx3 þ v2 ðx3 ; tÞ u_2 ðx3 ; tÞ
0.7 1
0.0040.03 0.0010.01
(8)
(9)
where aðx3 Þ and a0 ðx3 Þ indicate the axial and tangential induction coefficients, determined by Blade Element Momentum (BEM) (Hansen, 2008). vj ðx3 ; tÞ indicate the components in the moving frame of reference of the rotational sampled turbulence at the position x3 , at the time t. u_j ðx3 ; tÞ denote the corresponding moving frame components of the velocity vector of the blade at the abscissa x3 , due to elastic deformations of the tower and the blade. The instantaneous angle of attack is given as:
αðx3 ; tÞ ¼ φðx3 ; tÞ þ βðtÞ þ δðx3 Þ
(10)
where φðx3 ; tÞ indicates the flow angle defined as, see Fig. 3:
tanφðx3 ; tÞ ¼
Table 1 Typical reduced frequencies for wind turbine operations. k
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi V 21 ðx3 ; tÞ þ V 22 ðx3 ; tÞ
Fig. 3 shows the wind velocities in the rotor plane. V1 ðx3 ; tÞ and V2 ðx3 ; tÞ are the instantaneous wind velocities in local x1 -direction and x2 -direction, respectively, seen by an observer fixed to the moving (x1 , x2 , x3 )-coordinate system. These are given as:
where v1;j ðx3 ; tÞ; v2;j ðx3 ; tÞ; v3;j ðx3 ; tÞ are rotational sampled turbulence components at the position x3 in jth blade, in the moving frames of reference. v1;j ðX; tÞ; v2;j ðX; tÞ; v3;j ðX; tÞ are rotational sampled turbulence components at the same position for blade j in the fixed frame of refer-
x3 =L
(7)
V1 ðx3 ; tÞ V2 ðx3 ; tÞ
(11)
The lift coefficient cl ðαÞ and drag coefficient cd ðαÞ are calculated based on the instantaneous angle of attack together with 2D profile data, where the delay effect on the loads due to change of α has not been considered. 175
B. Chen et al.
Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 171–179
The damping forces of the eddy current damper can be presented by:
The aerodynamic force along the blade is calculated through the BEM method with Prandtls tip loss factor and Glauert correction (Hansen, 2008). Non-linear aeroelasticity is also considered by introducing the local deformation velocities of the blade into the calculation of the flow angle and the angle of attack, as expressed in Eq. (9). As a result, high aerodynamic damping is introduced in the blade flap-wise and the fore-aft tower vibrations, but relatively low aerodynamic damping in the blade edgewise and the lateral tower vibrations (Chen et al., 2017b).
where q_ r ðtÞ denotes the torsional degree of freedom velocity vector of the section i and section j, Cr denotes the eddy current damping matrix in the torsional degree of freedom of these two nodes. cr denotes the equivalent damping coefficient.
3. Viscous damper design
Cr ¼ c r
f r ðtÞ ¼ Cr q_ r ðtÞ
Fig. 6 shows schematic diagram of the proposed eddy current damper, which is placed at the inner part of the blade. It is well known that when a conductive material is subjected to a time-varying magnetic flux, eddy currents are generated in the conductor. These eddy currents circulate inside the conductor generating a magnetic field of opposite polarity as the applied magnetic field (Sodano and Bae, 2004). The interaction of the two magnetic fields causes a torque that resists the change in magnetic flux. The torque is proportional to the slip r_6;j r_6;i , i.e. the difference between the angular velocity of the adjacent shafts, corresponding to a linear viscous damper with the damper constant cr . The value of the damper constant can easily be changed by adjusting the magnetic field, making an optimal semi-active control to changing wind load conditions possible. The damper device can be separated into two parts: the magnetic part and the conducting part. The magnetic part is attached to the inner shaft of the length a fixed at node i. The conducting part consist of two parts: the conducting tube and the transmission shaft of the length b fixed at node j. The damping torque to be transmitted to node j is much smaller than the elastic torque transmitted by the blade structure. Hence, the dimension of the shaft is correspondingly small. In order to make maintenance of the damper device, the support node i should be placed close to the nacelle. The length b of the shaft corresponding to the support node j is determined, so the performance of the damper is optimized. If the damper is not activated corresponding to cr ¼ 0, the blade renders into flutter in the fundamental torsional mode. At the other extreme corresponding to cr ¼ ∞, the rotations are locked at node j, and no mechanical energy is extracted from the blade. In this case the blade renders into flutter with a mode fixed at node j and with a critical rotational speed Ωcr of the rotor somewhat larger than for cr ¼ 0. It turns out that an optimal value of cr exists, where a maximum value of Ωcr is achieved and a maximum. The optimal value of cr depends both on the position of node j as indicated by the parameter b and on the turbulence in the incoming flow to the rotor as defined by the mean wind velocity V0 and the turbulence intensity I. V0 will change the pitch angle βðtÞ, and hence may influence the flutter critical rotational speed of the rotor.
1 1 ; 1 1
(12)
q_ r ðtÞ ¼
r_6;i r_6;j
(13)
4. Numerical example Here the 907-DOF model has been used for analyzing flutter stability. According to Hansen's previous study (Hansen, 2002), the dynamics of supporting structure may affect the aerodynamic damping of blade and also the risk of flutter. However, results from his later researches (Hansen, 2004a, 2004b) shown that the single-blade model and full-wind-turbine model analyses for a particular turbine predicts similar flutter limits, suggesting that the single-blade flutter analysis is sufficient for predicting the flutter limits. It was also demonstrated that the frequency characteristics of blade are less affected by the supporting structure dynamics. Nevertheless, we used this full-wind-turbine model in the present study, partly because this model has been fully developed and widely used in our previous works (Zhang, 2015), and partly because we want to mimic the real situation. 4.1. Parameters of DTU 10 MW wind turbine Data from the DTU 10 MW reference wind turbine (Bak et al., 2013) have been used to calibrate the structural model. The constant parameters employed in the 907-DOF wind turbine model are calculated and provided in Table 2. Table 2 Parameters in the 907-DOF wind turbine model. Parameter
Value
Unit
Parameter
Value
Unit
L s h Jr Jg kr kg
89.166 5 119 1.5696 108 2328.24 3.48 109 1.392 106
m m m kgm2 kgm2 Nm/rad Nm/rad
Ω g a θmax ρ
1.0053 9.81 4.815 0.04 1.0 103
rad/s m/s2 m rad kg/m3
Fig. 6. The eddy current damper.
176
B. Chen et al.
Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 171–179
Table 3 Natural frequencies of the 907-DOF wind turbine model. Item
Value
Unit
Description
ωf 1 ωf 2 ωf 3 ωe1 ωe2 ωe3 ωθ1 ωfa ωl
4.15 11.56 23.86 6.09 18.50 39.18 37.77 2.75 2.75
rad/s rad/s rad/s rad/s rad/s rad/s rad/s rad/s rad/s
1st flapwise frequency of blade 2nd flapwise frequency of blade 3rd flapwise frequency of blade 1st edgewise frequency of blade 2nd edgewise frequency of blade 3rd edgewise frequency of blade 1st torsional frequency of blade 1st tower fore-aft frequency 1st tower lateral frequency
The main natural frequencies for the 907-DOF wind turbine model, relevant for the analysis below, are listed in Table 3. In this model, the 1st flapwise frequency, 1st edgewise frequency and the 1st torsional frequency have been chosen to generate the structural damping matrix Cst based on the Cauchey's damping model (Vatne, 2011; Bak et al., 2013). Furthermore, we also assumed that the structural damping ratios corresponding to these three modes (the 1st flapwise mode, the 1st edgewise mode and the 1st torsional mode) are the same. Since the aeroelastic mode and structural mode of wind turbine blades are different when they operate on the rotor, as Campbell diagram shows (Hansen, 2004a), the blades of wind turbine do not vibrate in pure blade mode, i.e. flapwise, edgewise and torsional modes may couple in the rotor whirling modes of the wind turbine (Hansen, 2003). The aeroelastic modes coupling through the aerodynamic forces may lead to strong interactions, such as classic flutter. In this paper, the numerical simulation result shows it is the third flapwise bending mode that couples with the first torsional mode in flutter. Actually, according to flutter criteria in classic flutter theory (Hansen, 2007), the natural frequencies of a torsional mode and a flapwise bending mode must be relatively close for them to couple in a flutter mode. In the present case, the first torsional mode frequency is indeed very close to the third flapwise bending mode with a frequency ratio around 1.58. In the following sections, the effects of different parameters on the optimal value of the flutter critical rotational speed Ωcr are discussed in detail. These parameters include the level of structural damping ratio ζ0 , the level of the turbulence as measured by the turbulence intensity I, the transmission shaft length of the damper b, and the mean wind velocity V0 . Since the main load components of wind turbine blades comes from the rotational speed of the rotor and the turbulence component is numerically too small in comparison with other load components, further discussion about turbulence effect on flutter limits is not elaborated in this paper but will be the topic of follow-on studies.
Fig. 7. Time series of torsional angle of the blade tip at ccstr ¼ 0, V0 ¼ 15m=s, Lb ¼ 0:4268, ζ 0 ¼ 0:005, I ¼ 0:1.
Fig. 8. Relative flutter critical rotational speed as a function of the relative damping ratio at V0 ¼ 15m=s, Lb ¼ 0:4268. ____ , turbulence intensity I ¼ 0; , turbulence intensity I ¼ 0:1; ⋯, turbulence intensity I ¼ 0:3.
deviations from the mean value θmax ¼ μθ þ 3 σ θ has been chosen as a compromised to limit the probability of false alarms due to elastic exceedances of the threshold level, and to prevent structural damages at the onset of flutter (Chen et al., 2017a). Fig. 8 shows the influence of the flutter critical rotational speed Ωcr as a function of cr for different values of turbulence intensity I and structural damping ratio ζ0 , with V0 ¼ 15m=s and b ¼ 0:4268L. cr has been normalized with respected to cst , where cst indicates the component at the main diagonal of the structural damping matrix Cst at node j in the rotational degree of freedom. Ωcr has been nondimensionalized with respect to Ω. As seen, each curve has a peak at a specific relative damping ratio, which represents the maximum value of the critical rotational speed for flutter as discussed above. Furthermore, the flutter critical rotational speed is only marginally dependent on the turbulence intensity and will decrease with the increase of turbulence intensity. With the increasing value of structural damping ratio ζ0 , the peak will move in the direction of lower relative damping ratio. Fig. 9 shows the variation of Ωcr as a function of Lb for fixed value of ccstr , V0 , ζ0 and I. As can be seen, there is an optimal value of Lb, which maxi-
4.2. Result and discussion Here the rotor was purposefully operated at ever increasing rpm until the flutter boundary was reached and dramatic classical flutter oscillations were observed. Further, the simulation time at each increased rotational speed of rotor is long enough to obtain the steady-state vibration of blade tip. Fig. 7 shows the time series of the torsional degree of the blade tip under flutter critical conditions for fixed values of V0 , Lb, ccstr , ζ 0 and I. For convenience, the time history prior to that shown has been omitted, and only the initial segment of the unstable behaviour is displayed where flutter is evident. As can be viewed, as the envelope of the torsional angle exceeds the critical value θmax , the vibration starts to increase beyond limits, indicating instability. Here, the threshold value θmax has been set as the mean value of the time series of torsional angle of the blade tip plus three standard deviations under the stability situation for the following reason. Let μθ and σ θ indicate the mean value and the standard deviation of elastic deformations of the torsional rotation θðtÞ, caused by the turbulence alone. The indicated statistical moments and determined by ergodic sampling of a suitable line series before the onset of flutter. The value of 3 standard
mizes the value of
177
Ωcr Ω
(i.e. Ωcr ¼ 3.6 Ω), and it is almost 2.3 times the
B. Chen et al.
Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 171–179
5. Concluding remarks
Fig. 9. Variation of the flutter critical rotational speed as a function of V0 ¼ 15m=s, ζ 0 ¼ 0:005, I ¼ 0:1.
b L
at
cr cst
A new type of eddy current based ’electromagnetic’ linear viscous damper attached to the inner part of wind turbine blades has been proposed in the paper for enhancing the flutter stability. An investigation is carried out to evaluate the effectiveness of the proposed viscous damper in increasing the flutter critical rotational speed through a sophisticated 907-DOF aeroelastic model. The influence of the modal damping ratio of the structure, the turbulence intensity, the mean wind velocity and the transmission shaft length of the damper device on the flutter critical rotational speed have been discussed. It is concluded that if the transmission shaft length of the damper device b is almost half of the length of the blade L, the linear viscous damper at optimal tuning can increase the flutter critical rotational speed more than 100%. Too small or too large values of the equivalent damping coefficient can decrease the flutter critical rotational speed compared to the optimal tuned value. It is demonstrated that the mean wind velocity, the turbulence intensity and the structural damping of the torsional mode particularly in the flutter mode have little effect on the achieved enhanced flutter stability. All these results show that the proposed viscous damper is a promising damping device with the possibility of semi-active mode of operation and can be utilized in wind turbine blades for flutter control.
¼ 9,
Acknowledgements The first author gratefully acknowledges the financial support from the State's Key Project of Research and Development Plan (No. 2016YFE0127900) and the National Science Foundation of China (No. 51422806, No. 51278189). References Bak, C., Zahle, F., Bitsche, R., Kim, T., Yde, A., Henriksen, L., Natarajan, A., Hansen, M., 2013. Description of the DTU 10 MW Reference Wind Turbine. DTU Wind Energy Report-I-0092. C. Bak, F. Zahle, R. Bitsche, T. Kim, A. Yde, L. Henriksen, P. B. Andersen, A. Natarajan, M. H. Hansen, Design and performance of a 10 MW wind turbine, Wind Energy 124. Batchelor, G.K., 1953. The Theory of Homogeneous Turbulence. Cambridge university press. Belamadi, R., Djemili, A., Ilinca, A., Mdouki, R., 2016. Aerodynamic performance analysis of slotted airfoils for application to wind turbine blades. J. Wind Eng. Ind. Aerod. 151, 79–99. Bernhammer, L.O., van Kuik, G.A., De Breuker, R., 2016. Fatigue and extreme load reduction of wind turbine components using smart rotors. J. Wind Eng. Ind. Aerod. 154, 84–95. Bir, G., 2008. Multiblade coordinate transformation and its application to wind turbine analysis. In: ASME Wind Energy Symposium, pp. 1–15. Bir, G., Jonkman, J., 2007. Aeroelastic instabilities of large offshore and onshore wind turbines. J. Phys. Conf. Ser. 75, 012069. IOP Publishing. Bisplinghoff, R.L., Ashley, H., Halfman, R.L., 2013. Aeroelasticity. Courier Corporation. Bottasso, C., Campagnolo, F., Croce, A., Tibaldi, C., 2013. Optimization-based study of bend–twist coupled rotor blades for passive and integrated passive/active load alleviation. Wind Energy 16 (8), 1149–1166. Bottasso, C.L., Croce, A., Gualdoni, F., Montinari, P., 2016. Load mitigation for wind turbines by a passive aeroelastic device. J. Wind Eng. Ind. Aerod. 148, 57–69. Chen, B., Hua, X.G., Zhang, Z., Basu, B., Nielsen, S.R.K., 2017. Monitoring of wind turbine blades for flutter instability. Struct. Monitor. Mainten. 4 (2), 115–131. Chen, B., Zhang, Z., Hua, X., Basu, B., Nielsen, S.R.K., 2017. Identification of aerodynamic damping in wind turbines using time-frequency analysis. Mech. Syst. Signal Process. 91, 198–214. I. E. Commission, et al., Wind Turbine Part 1: Design Requirements, IEC 61400–61401, International Electrotechnical Commission, Geneva, Switzerland. De Goeij, W., Van Tooren, M., Beukers, A., 1999. Implementation of bending-torsion coupling in the design of a wind-turbine rotor-blade. Appl. Energy 63 (3), 191–207. Dowell, E.H., Curtiss, H.C., Scanlan, R.H., Sisto, F., 1989. A Modern Course in Aeroelasticity, vol. 3. Springer. Fung, Y.C., 2002. An Introduction to the Theory of Aeroelasticity. Courier Corporation. Fung, Y.C., 2004. An Introduction to the Theory of Aeroelasticity. Dower. D. T. Griffith, T. D. Ashwill, reportThe Sandia 100-meter All-glass Baseline Wind Turbine Blade: SNL100-00, Sandia National Laboratories, Albuquerque, Report No. SAND2011–3779. Hansen, M., 2002. Vibrations of a three-bladed wind turbine rotor due to classical flutter. In: 40th AIAA Aerospace Sciences Meeting and Exhibit, 2002 ASME Wind Energy Symposium.
Fig. 10. Variation of the flutter critical rotational speed with the mean wind velocity V0 at cr b cst ¼ 9, L ¼ 0:4268, ζ 0 ¼ 0:005, I ¼ 0:1.
value, without the damper devise (i.e. Ωcr ¼ 1.6 Ω). The maximum value for b ¼ 0:4268L is Ωcr ≃ 3:3Ω. As seen on Fig. 9, slight increase to Ωcr ≃ 3:6Ω can be achieved by fixing the damper device at b ≃ 0:58L. However, this increase is at the expense of a longer torque shaft to the damper device. Fig. 10 shows the variation of Ωcr as a function of the mean wind velocity V0 for fixed values of ccstr , Lb, ζ0 and I. Evidently, Ωcr is found to decrease monotonically for increasing mean wind velocity V0 . However, Ωcr only decreases slightly with V0 , which is similar to Politakis et al.’s result (Politakis et al., 2008). Furthermore, as seen from Eq. (8) and (9), the resulting wind velocity V will increases with V0 , which in turn will increases the proneness to instability. In conclusion, Figs. 8 and 10 show that the optimal damper constant cr for optimal flutter stability only depends slightly on the incoming wind field to the rotor as characterized by the mean wind velocity V0 and the turbulence intensity I. 178
B. Chen et al.
Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 171–179 Moshfeghi, M., Shams, S., Hur, N., 2017. Aerodynamic performance enhancement analysis of horizontal axis wind turbines using a passive flow control method via split blade. J. Wind Eng. Ind. Aerod. 167, 148–159. Nayfeh, A.H., Mook, D.T., 1995. Nonlinear Oscillations. John Wiley and Sons. Ogata, K., 2010. Modern Control Engineering, fifth ed. Pearson. Pars, L.A., 1979. A Treatise on Analytical Dynamics. OX Bow Press. Politakis, G., Haans, W., Van Bussel, G., 2008. Suppression of Classical Flutter Using a ’Smart Blade’. American Institute of Aeronautics and Astronautics AIAA. Rafiee, R., Tahani, M., Moradi, M., 2016. Simulation of aeroelastic behavior in a composite wind turbine blade. J. Wind Eng. Ind. Aerod. 151, 60–69. Resor, B., Paquette, J., 2011. Uncertainties in prediction of wind turbine blade flutter. In: 49 th AIAA Aerospace Sciences Meeting and Exhibit. Resor, B.R., Owens, B.C., Griffith, D.T., 2012. Aeroelastic instability of very large wind turbine blades. In: European Wind Energy Conference Annual Event, Copenhagen, Denmark. Sodano, H.A., Bae, J.-S., 2004. Eddy current damping in structures. Shock Vib. Digest 36 (6), 469. Theodorsen, T., Mutchler, W., 1935. General Theory of Aerodynamic Instability and the Mechanism of Flutter. National Advisory Committee for Aeronautics, Washington, DC, USA. Vatne, S.R., 2011. Aeroelastic Instability and Flutter for a 10 MW Wind Turbine. Master’s thesis. Institutt for energi-og prosessteknikk. Veers, P., Bir, G., Lobitz, D., 1998. Aeroelastic tailoring in wind-turbine blade applications. In: Proceedings, Windpower, 98, pp. 291–304. Zhang, Z., 2015. Passive and Active Vibration Control of Renewable Energy Structures. Ph.D. thesis. Aalborg University, The Faculty of Engineering and Science, Division of Reliability, Dynamics and Marine Engineering. Zhang, Z., Nielsen, S.R.K., Blaabjerg, F., Zhou, D., 2014. Dynamics and control of lateral tower vibrations in offshore wind turbines by means of active generator torque. Energies 7 (11), 7746–7772.
Hansen, M., 2003. Improved modal dynamics of wind turbines to avoid stall-induced vibrations. Wind Energy 6 (2), 179–195. Hansen, M., 2004. Stability analysis of three-bladed turbines using an eigenvalue approach. In: 2004 ASME Wind Energy Symposium. American Institute of Aeronautics and Astronautics; American Society of Mechanical Engineers, pp. 192–202. Hansen, M.H., 2004. Aeroelastic stability analysis for for wind turbines using an eigen value approach. Wind Energy 7 (2), 133–143. Hansen, M.H., 2007. Aeroelastic instability problems for wind turbines. Wind Energy 10 (6), 551–577. Hansen, M.O.L., 2008. Aerodynamics of Wind Turbines, second ed. Earthscan. Hayat, K., Ha, S.K., 2015. Load mitigation of wind turbine blade by aeroelastic tailoring via unbalanced laminates composites. Compos. Struct. 128, 122–133. Hayat, K., de Lecea, A.G.M., Moriones, C.D., Ha, S.K., 2016. Flutter performance of bend–twist coupled large-scale wind turbine blades. J. Sound Vib. 370, 149–162. Jonkman, B.J., 2009. Turbsim User's Guide: Version 1.50. National Renewable Energy Laboratory Colorado. Krenk, S., Svendsen, M.N., Høgsberg, J., 2012. Resonant vibration control of three-bladed wind turbine rotors. AIAA J. 50 (1), 148–161. Larsen, J.W., Nielsen, S.R.K., Krenk, S., 2007. Dynamic stall model for wind turbine airfoils. J. Fluid Struct. 23 (7), 959–982. Leishman, G.J., 2006. Principles of Helicopter Aerodynamics. Cambridge university press. Lobitz, D.W., 2004. Aeroelastic stability predictions for a MW-sized blade. Wind Energy 7 (3), 211–224. Lobitz, D.W., 2005. Parameter sensitivities affecting the flutter speed of a MW-sized blade. J. Sol. Energy Eng. 127 (4), 538–543. Lobitz, D.W., Veers, P.S., 1998. Aeroelastic behavior of twist-coupled HAWT blades. In: ASME/AIAA Wind Energy Symposium, Reno, NV, pp. 75–83.
179
Contents lists available at ScienceDirect
Journal of Wind Engineering & Industrial Aerodynamics journal homepage: www.elsevier.com/locate/jweia
Enhancement of flutter stability in wind turbines with a new type of passive damper of torsional rotation of blades Bei Chen a, Zili Zhang b, Xugang Hua a, *, Søren R.K. Nielsen c, Biswajit Basu d a
Key Laboratory for Wind and Bridge Engineering, Hunan University, 410082 Changsha, China Department of Engineering, Aarhus University, 8000 Aarhus, Denmark Department of Civil Engineering, Aalborg University, 9220 Aalborg, Denmark d Department of Civil, Structural and Environmental Engineering, School of Engineering, Trinity College Dublin, Dublin 2, Ireland b c
A R T I C L E I N F O
A B S T R A C T
Keywords: Horizontal axis wind turbines Flutter analysis Viscous dampers Vibration control
Classical flutter of wind turbine blades is defined as an aeroelastic instability where a torsional blade mode couples to a flapwise bending mode to result in a mutual rapid growth of the amplitude of the flapwise and torsional motions beyond certain rotational speed of the rotor. Till date, flutter instability of wind turbine rotor has not been considered to be a serious problem. However, with the advent of larger turbines fitted with relatively softer blades, classical flutter may become a more important design consideration. As an economically viable alternative to increasing the torsional rigidity (which is more expensive), the use of a type of viscous damper mounted inside the wind turbine blades is proposed to promote the flutter critical rotational speed of the rotor. The flutter critical rotational speed is determined by using a time-domain aeroelastic simulation technique based on the quasi-steady aerodynamic forces and also considering turbulence intensity of incoming wind. The optimal damper location and the optimal damping constant values are investigated for the DTU 10 MW wind turbine. It is shown that at an optimal tuning of the damper, the flutter critical rotational speed is only marginally dependent on the turbulence intensity of the incoming wind field, on the mean wind velocity and on the structural damping of the flutter torsional eigenmode. The novel tuned damper device proposed in this paper may increase the flutter critical rotational speed with more than 100%.
1. Introduction It is known that three typical aeroelastic instabilities may occur in modern commercial wind turbines: divergence, stall-induced vibration and classical flutter. Divergence is a static phenomenon wherein the resulting wind velocity (consist of both the ambient wind speed and the rotational speed of the rotor) becomes large enough so that the load produced for an incremental angle of attack change due to blade twisting is greater than the reaction load produced by the elastic restoring forces for the same amount of twist (Lobitz and Veers, 1998). Stall-induced vibration is a single mode phenomenon characterized by predominantly flapwise blade oscillations, and normally occurs at high angle of attack causing separation of the boundary layer on the backside of the profile. The stall induces negative aerodynamic damping, related to a l negative gradient ∂c ∂α of the lift coefficient with the angle of attack (Hansen, 2007). The separation point may be stationary, which is referred to as static stall. At the core of dynamic stall the separation point
moves on the backside of the profile, the boundary layer may even vary between fully attached or fully detached during part of the oscillation (Larsen et al., 2007). Pitch-regulated variable speed wind turbines normally do not operate in stall and the risk of stall-induced vibration is not as serious as for stall-regulated wind turbines. Flutter is a well-known destructive instability from aircraft industry but historically has not been an issue in modern commercial wind turbines’ design. Commensurately, among the utility-scale wind turbines that have been built, rarely has classical flutter ever been observed. In short, classical flutter of a beam-like structure such as wind turbine blades is explained as an unfavourable aerodynamic coupling between flapwise and torsional vibrations of blades in the presence of incident aerodynamic forces, and appears in a rapid growth of the amplitudes of the flapwise and torsional motions, which may leads to major structural failure of the blades. In a steady homogeneous wind field the onset of flutter takes place, if the average power supplied to the structure by the self-induced loads exceeds the average power dissipated by the structural
* Corresponding author. E-mail addresses: [email protected] (B. Chen), [email protected] (Z. Zhang), [email protected] (X. Hua), [email protected] (S.R.K. Nielsen), [email protected] (B. Basu). https://doi.org/10.1016/j.jweia.2017.12.011 Received 22 August 2016; Received in revised form 7 December 2017; Accepted 7 December 2017 0167-6105/© 2017 Elsevier Ltd. All rights reserved.
B. Chen et al.
Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 171–179
earlier work on a 20 KW HAWT with stiff, aeroelastically tailored 5 m blades, Lobitz and Veers (1998) predicted flutter critical rotational speed that were six times the nominal rotor speed (ΩΩcr ¼ 6), rendering flutter a moot issue. The predicted flutter limits of a SNL 9-meter CX-100 experimental blade were shown to be four times higher than the nominal rotor speed (ΩΩcr ¼ 5) (Resor and Paquette, 2011). Study on 1.5 MW-sized WindPACT wind turbine with 34 m blades has a safety margin of 2–2.5 (ΩΩcr ¼ 2:0 2:5) (Lobitz, 2004). Similar result has also been found by Hansen for a 5 MW NREL wind turbine (Hansen, 2007). Further, for the 10 MW NOWITECH wind turbine with 89 m blades, the critical rotational speed for flutter is found to be about 1.6 times the nominal rotor speed (ΩΩcr ¼ 1:6) (Vatne, 2011). However, result from a 13.2 MW wind turbine with 100 m blades in length showed that there is little or no margin on flutter speed (ΩΩcr ¼ 1:0 1:1) reported by Griffith and Resor (Griffith and Ashwill; Resor et al., 2012). The trend shows that the ratio of critical rotational speed for flutter to the nominal rotor speed ΩΩcr drops significantly as the blades grows in length from 5 m to 100 m. Therefore, it is believed that with a likely reducing blade torsional stiffness with increasing blade radius (considering current design trends) the critical rotational speed for flutter Ωcr will be approach the nominal rotational speed Ω of the rotor (Politakis et al., 2008). This also means that aeroelastic instabilities, especially flutter, would then become one of the principal design drivers (Hansen, 2007). Great efforts have been devoted to mitigating loads and improving the aeroelastic stability of wind turbines (Belamadi et al., 2016; Bottasso et al., 2016; Moshfeghi et al., 2017; Bernhammer et al., 2016). Due to the limitation of space in wind turbine blades, few works have been carried out on suppressing flutter of wind turbine blade, especially using mechanical countermeasures, such as tuned mass damper (TMD). Hence, the adaptive blades that involving the use of aeroelastic tailoring (also known as bend-twist coupling (BTC)), wherein the blade twists as it bends under the action of aerodynamic loads to shed wind-induced loads and increase the aeroelastic stability of wind turbine has been investigated by a lot of researchers (De Goeij et al., 1999; Veers et al., 1998; Rafiee et al., 2016). Even though the BTC proved to be effective in mitigating the fatigue loads for large-scale wind turbine blades, it may increases the blades proclivity for flutter at the same time (Bottasso et al., 2013; Hayat and Ha, 2015; Hayat et al., 2016). Lobitz's (Lobitz and Veers, 1998) research shown that the flutter speed of the rotor can be increased by around 11% at the optimal values of the coupling coefficient for 20 kW HAWT, since the blade is less stable as the extreme values of the coupling coefficient are approached due to the divergency limit. Further, Hayat (Hayat et al., 2016) extended Lobitz's study using a beam element model to simulate the composite laminate material structure. Results shown the flutter speed of the rotor can be increased by about 7:6 9:5% using lighter and stiffer carbon fibers. As mentioned above, there have been several researches reported on increasing flutter limits using BTC. However, their effectiveness is very limited. Inspired from this, a new ’electromagnetic’ type of linear viscous damping device is proposed in this paper in order to mitigate the torsional and flapwise vibration of the blades, hence increasing the flutter critical rotational speed of the rotor. Firstly, a wind turbine model is presented, which considers the coupling between aerodynamic, elastic and inertial loads of the blades. The parameters of the model have been calibrated to the DTU 10 MW wind turbine (Bak et al., 2013; Bak et al.,). The elastic deformations of the blades are coupled via the tower motion, which is modeled by a reduced five degree-of-freedom model, and via the drive train, which is modeled by a two degree-of-freedom model. The pitch of the blades are controlled by a collective PI controller with feedback from the rotor speed. Furthermore, a linear viscous damper model is provided. The optimal parameters of the viscous damper are calculated. Finally, the sensitivity study of this damping system to the turbulence, the mean wind velocity, the structural damping and the parameters of this damper is performed.
damping and possibly attached active or passive damping mechanism mounted on the structure (Politakis et al., 2008). Based on the studies of Bisplinghoff, Fung and Hansen (Hansen, 2004a, 2007; Bisplinghoff et al., 2013; Fung, 2004), the risk of classical flutter of wind turbine blades with attached boundary layer may rise if the following three main criteria are satisfied: the frequency ratio between a flapwise bending mode and a torsional mode is close to one, the rotational speed of the rotor is sufficiently high and the position of the shear center of the airfoil is placed between the center of mass and aerodynamic center. Considering the first criterion, the risk of flutter increases as the ratio between the first torsional mode and one of the lower flapwise modes approaches one, such that they can couple in a flutter mode. As a reference point, for the 1.5 MW Horizontal Axis Wind Turbine (HAWT) with 35 m blades, it was the second flapwise bending mode and the first torsional blade mode which are coupled in flutter (Lobitz, 2004). For the 5 MW National Renewable Energy Laboratory (NREL) turbine, used as a benchmark in numerous analyses of wind turbines, the flutter instability took place between the third flapwise bending mode and the first torsional blade mode (Hansen, 2007). Similar result has also been found by Bir for a 5 MW NREL wind turbine (Bir and Jonkman, 2007). For the second criterion, the rotational speed of the rotor Ω must be sufficiently high, resulting in higher energy in the aerodynamic forces which ’feed’ the flutter mode. Finally, the onset of flutter requires that the shear center is placed between the aerodynamic center and the center of mass. The inertial load is acting at the center of mass, whereas the aerodynamic load is acting at the aerodynamic center. Since the aerodynamic and the inertial loads approximately are in counter phase, they provide a maximum of elastic torsion with the indicated position of the shear center. Lobitz's work shown the flutter speed limit decreases when the center of mass is moved towards the trailing edge of the blade. Furthermore, when the center of mass is ahead of the shear center of the blade profile, flutter is unlikely at any rotational speed of the rotor (Lobitz, 2005). Since the flutter mode is not known beforehand, thus the flutter can only be predicted by aeroelastic eigenvalue analysis or time domain simulation. Several studies have been carried out on wind turbine stability analysis, especially for flutter. Hansen (2003; 2004a; 2004b) proposed a multi-blade coordinate transformation (MBC, (Bir, 2008)) based eigenvalue analysis. Firstly the time-independent aeroelastic motion equations which couple the structural motion with the dynamics of the aerodynamic forces were derived using Lagrange's equations and MBC transformation. Then the eigenvalue problem which is ill-conditioned based on these autonomous equations was solved by using a modal expansion of the structural DOFs based on the undamped modes of the turbine at standstill. This enabled the computation of the natural frequencies, logarithmic decrements, and mode shapes of the aeroelastic turbine modes. The lowest nominal rotational speed of the rotor which corresponding to a negative damping coefficient exhibits in eigenmode is designated as the flutter critical rotational speed. Similar method has been used by Lobitz (2004). Alternatively, the aeroelastic stability may be checked by a pure mathematical approach of using Floquet theory on the non-transformed linearized equations (Nayfeh and Mook, 1995). However, this approach should be avoided since it doesn't have any physical meaning (Hansen, 2003). Furthermore, the stability analysis may also be performed by observing a certain envelope of the flapwise or torsional motions, indicating the evolution with time of the amplitude of the involved response quantities. Flutter is defined to take place, where the envelope of the torsional angle exceeds a certain critical value θmax . This approach is useful if the instability is dominated by nonlinear mechanisms, either structural or aeroelastic (Chen et al., 2017a). According to previous studies, the critical rotational speed for flutter turns out to be higher than the nominal rotor speed. Smaller wind turbines do not experience flutter due to their high torsional angular frequency in comparison to the angular eigen frequency in the critical flapwise bending mode of the blades (Lobitz, 2004; Hansen, 2008). In
172
B. Chen et al.
Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 171–179
2. Modelling of flutter induced vibrations of wind turbines
Generator shaft kg
2.1. Wind turbine model
q7 (t)
The motion of the tower is described in a fixed, global ðX1 ; X2 ; X3 Þ-coordinate system with origin O at the interface to the foundation, see Fig. 1. The displacement of each blade is described relative to a moving, local ðx1 ; x2 ; x3 Þ-coordinate system with origin O' at the center of the hub. The position of the local coordinate system attached to blade j is merely specified by the azimuthal angle Ψj ðtÞ, which is considered positive when rotating clockwise seen from an upwind position. Then, the azimuthal angles Ψj ðtÞ may be given as:
2π 2π Ψj ðtÞ ¼ Ωt þ q6 ðtÞ þ ðj 1Þ≃Ωt þ ðj 1Þ; 3 3
j ¼ 1; 2; 3
Mr (t)
Rotor of generator
q6 (t)
Jr
Mg (t)
Mg (t) kr
Stator of generator
Rotor shaft
Fig. 2. 2-DOF model of flexible drive train with odd number of gear stages. Definition of degrees of freedom q6 ðtÞ and q7 ðtÞ.
vector q0 ðtÞ. Based on q0 ðtÞ and q_ 0 ðtÞ the kinetic and strain energy of the tower and the drive train minus the rotor may be calculated. The kinetic energy related to the rotor is accounted for at the analysis of the blades. Each blades is modeled as a Bernoulli-Euler beam with variable mass per unit length μðx3 Þ and variable bending stiffness in the flapwise and edgewise directions. Correspondingly, deformations from the shear forces are ignored. The elastic (bending) center E of all cross-sections along the blade are assumed to be placed along the x3 -axis. For each cross section a principal axis ðy1 ; y2 ; y3 Þ-coordinate is defined as shown in Fig. 3. The y3 -axis is unidirectional to the x3 -axis, and the orientation of the y2 -axis is defined by the angle δðx3 Þ, indicating the angle from the y2 -axis into the x2 -axis in the positive x3 -direction. δðx3 Þ represents the pretwist of the blade. The aerodynamic center A, the shear center S and the mass center G are also shown on Fig. 3. The position of the aerodynamic center is assumed to be fixed in time, either corresponding to the situation of attached boundary layers on the back-side of the profile, or to deep stall, where the boundary layers are totally detached. The external aerodynamic lift force per unit length pl and drag force per unit length pd acting in the aerodynamic center is transferred to the shear center in order to decouple the equations of motion for the elastic torsion and bending displacement, when formulated in the ðy1 ; y2 ; y3 Þ-coordinate system. The inertial loads, including the gravity load, are acting in the mass center G. The gravity loads are referred to the shear center for the same reason as for the aerodynamic loads. The components of velocity vector and angular velocity vector of the mass center need to be interpolated from the time derivates of the degrees-of-freedom of the element. Based on these interpolation functions the kinetic energy, as well as the mass matrix and the gyroscopic damping and stiffness matrices of the element can be calculated. The fact that the inertial loads are acting in the mass center G, and not in the elastic center E, introduces couplings in the equation of motion of importance for the flutter stability. Fig. 4 illustrates the discretization of a blade into a number n of finite elements, including the global node numbering. x3;j indicates the position of node j of the element on the x3 -axis. Also shown is the beam element j between the nodes j and j þ 1 and the definitions of the twelve degrees of freedom r1;j ; …; r12;jþ1 in the principal axis coordinate system. lj ¼ x3;jþ1 x3;j indicates the element length. The degrees-of-freedom specifying the elastic deformation of blade j relative to the hub in the
(1)
where Ψj ðtÞ indicates the azimuthal angle of blade j ¼ 1; 2; 3, Ω is the nominal rotational speed of the rotor. q6 ðtÞ is a degree of freedom specifying the deviatoric rotation of the hub on top of the nominal rotational angle Ωt. q6 ðtÞ is supposed to be numerical small quantities compared to Ωt, hence q6 ðtÞ can be ignored in Eq. (1). A five degree-of-freedom attached modal model is applied for the description of the tower motion defined by the displacements q1 ðtÞ and q2 ðtÞ of the tower top in the global X1 - and X2 -directions, and the rotations q3 ðtÞ, q4 ðtÞ, q5 ðtÞ. L indicates the length of the blade from the hub, h is the height of the tower from the base to the nacelle, and the horizontal distance from the center of the tower top to the origin O' of the moving coordinate systems is denoted s. The drive train is modeled by the degrees of freedom q6 ðtÞ and q7 ðtÞ as shown in Fig. 2, indicating the rotational angle of the rotor and generator rotor, respectively. Mr ðtÞ and Mg ðtÞ denote the workconjugated rotor torque and generator torque. Jr and Jg denote the mass moment of inertia of the rotor and the generator, and kr and kg denote the St. Venant torsional stiffness of the rotor and generator shafts. For more details about the tower and the drive train modelling, the reader is referred to Zhang (Zhang et al., 2014; Zhang, 2015). The degrees of freedom q1 ðtÞ; …; q7 ðtÞ are assembled in the column
Fig. 1. Definition of the fixed and the moving frames of reference and the degrees of freedom q1 ðtÞ; …; q5 ðtÞ. Fig. 3. Definition of principal axis coordinate system.
173
B. Chen et al.
Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 171–179
Fig. 4. Nodal numbering and definition of degrees of freedom of a beam element.
2
3 qðtÞ _ 5 zðtÞ ¼ 4 qðtÞ βðtÞ
local ðx1 ; x2 ; x3 Þ-coordinate system, are assembled in the column vector qj ðtÞ; j ¼ 1; 2; 3 of dimension 6n. Small elastic vibrations are assumed, so that linear beam theory may be assumed with due consideration of the geometrical stiffening in the x3 -direction, and the geometrical softening in transverse direction. For a more thorough elaboration, including the mathematical derivation of the displacement and velocity vector of the blades, the reader is referred to Zhang et al. (2014), Zhang (2015). The degrees of freedom of the system are assemble in the column vectors qðtÞ of dimension 7 þ 18n:
The parameter βðtÞ indicates a full-span collective pitch angle with a time delay determined by the first order filter:
_ ¼ 1 ðβ0 ðq6 ; q_6 Þ βðtÞÞ βðtÞ τ
(2)
2.2. Aeroelastic load model
Based on the applied discretization of the blades the kinetic energy _ and the potential energy UðqÞ of the tower, drive train and blades Tðq; qÞ _ may be calculated as a function of q and qðtÞ. Then, the equations of motion of the tower, the drive train and the blades follow from the EulerLagrange stationarity condition of analytical dynamics (Pars, 1979):
_ _ d ∂Tðq; qÞ ∂Tðq; qÞ ∂UðqÞ _ tÞ þ ¼ fðq; q; dt ∂q_ ∂q ∂q
(5)
βðtÞ is considered positive, when rotating in the positive x3 -direction, see Fig. 3. β0 ðq6 ; q_6 Þ indicates the pitch demand, which is modeled as PI controller with feed-back from q6 ðtÞ and q_6 ðtÞ (Ogata, 2010). The parameter τ specifies the reaction time of the pitch actuators.
2
3 q0 ðtÞ 6 q1 ðtÞ 7 7 qðtÞ ¼ 6 4 q2 ðtÞ 5 q3 ðtÞ
(4)
In accordance with (Commission et al, ), the turbulence model is constructed based on Taylors hypothesis of frozen turbulence, corresponding to a frozen turbulence field that is convected into the rotor in the global X1 direction with a mean velocity V0 . The frozen field is
€ þ CðtÞqðtÞ _ þ KðtÞqðtÞ ¼ fðz; tÞ ⇒ MðtÞqðtÞ
(3)
assumed to be a zero mean homogeneous and isotropic stochastic field, with a spatial covariance structure as elaborated in (Batchelor, 1953). Calibrated from the theoretical covariance structure of the turbulence, the first order AR model (Krenk et al., 2012) performs a first-order filtering of the white noise input, resulting in continuous, non-differentiable sample curves of the turbulence field at the rotor plane in the fixed frame of reference (Zhang et al., 2014). The present turbulence model is in good agreement with one of TurbSim's turbulence models (i.e. the Von Karman's turbulence model) (Jonkman, 2009; Krenk et al., 2012), but is highly efficient numerically and easy to incorporate into the 907 DOF model. The model generates correlated time histories corresponding to such frozen turbulence field at the points of the mesh shown in Fig. 5 consisting of 8 radial lines with 6 points in each. As seen, the fixed frame components of the convected turbulence are generated at nn ¼ na ⋅nr þ 1 nodal points at each time step t ¼ 0; △t; 2△t; …, where na denotes the
Because the elastic deformations have been assumed to be linear, the mass matrix MðtÞ, the damping matrix CðtÞ and the stiffness matrix KðtÞ merely depend q0 ðtÞ and q_ 0 ðtÞ, and more specific on q6 ðtÞ and q_6 ðtÞ. As can be seen in Eq. (1), the q6 ðtÞ and q_6 ðtÞ are supposed to be numerically small quantities compared to Ωt and Ω, and can be ignored. Hence, the mass matrix MðtÞ, the damping matrix CðtÞ and the stiffness matrix KðtÞ can be expressed as a function of t. Furthermore, the mass matrix MðtÞ and the damping matrix CðtÞ are symmetric and periodic, but due to the gyroscopic stiffness contributions the stiffness matrix KðtÞ is asymmetric periodic matrix. fðz; tÞ on the right hand side of Eq. (3) indicates the load vector conjugated to qðtÞ from self-induced aeroelastic load, pitch control load and wind load. zðtÞ indicates the state vector of the system defined as:
174
B. Chen et al.
Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 171–179
is quasi-steady, i.e. the unsteady effects are not pronounced; if k > 0:05, the unsteady effects cannot be neglected. Based on our calculation, the reduced frequency k at the outer part of the blade is much smaller than 0.05, as shown in Table 1. Looking at the values of the reduced frequencies in Table 1, it is possible to see that the reduced frequencies k in the range 0.001–0.03. Actually, Lobitz's research on a MW-sized wind turbine blade demonstrated that the computed flutter limit of a blade using quasi-steady theory is conservative than that estimated using unsteady (Theodorsen) theory. Furthermore, when flutter speed obtained from quasi-steady theory is similar to the nominal rotational speed of the rotor, the results are no more reliable and may adversely affect blade design (Lobitz, 2004). In other words, for those wind turbines fitting with relatively softer blades (e.g. 13 MW HAWT, 20 MW HAWT or HAWTs with bigger rated capacity), the quasi-steady theory may be unapplicable. In the present study, we haven't considered the unsteady aerodynamics that can be derived from Theodorsen's theory (Lobitz, 2004; Theodorsen and Mutchler, 1935; Fung, 2002; Dowell et al., 1989), because the main goal of current research is to investigate the effectiveness of the linear viscous damping device on mitigating flutter. However, the unsteady aerodynamics will be considered in our future studies. Assuming locally 2D-flow around the profile the lift and drag load per unit length pl ðx3 ; tÞ and pd ðx3 ; tÞ at the position x3 at the time t become:
9 1 > pl ðx3 ; tÞ ¼ cl ðαÞ ρV 2 ðx3 ; tÞcðx3 Þ > = 2 > 1 > pd ðx3 ; tÞ ¼ cd ðαÞ ρV 2 ðx3 ; tÞcðx3 Þ ; 2
Fig. 5. Nodal numbering and definition of degrees of freedom of a beam element (Zhang et al., 2014).
number of radial lines from the center and nr denotes the number of equidistantly placed nodes along a given radial line. Then, the fixed frame components of the rotational sampled turbulence vector on each blade with the azimuthal angle Ψj are calculated through the linear interpolation of the turbulence of the adjacent radial lines according to the position of the blade at each time step. Finally, the moving frame components of the rotational sampled turbulence can be determined through the following coordinate transformation equation:
3 2 1 0 v1;j ðx3 ; tÞ 4 v2;j ðx3 ; tÞ 5 ¼ 4 0 cosΨj 0 sinΨj v3;j ðx3 ; tÞ 2
3 32 0 v1;j ðX; tÞ sinΨj 54 v2;j ðX; tÞ 5 cosΨj v3;j ðX; tÞ
where cd ðαÞ and cl ðαÞ indicate the lift and drag coefficients obtained from static 2D wind tunnel test data, αðx3 ; tÞ is the instantaneous angle of attack, ρ is the mass density of air, Vðx3 ; tÞ is the resulting wind velocity and cðx3 Þ is the chord length of the profile at the coordinate x3 . The mean wind velocity V0 is directed in the global X1 -direction, and is assumed to be constant over the rotor area. Further, the nominal rotational speed of the rotor Ω is assumed to be constant in time. Then, the resulting wind velocity on the profile becomes:
(6) Vðx3 ; tÞ ¼
ence with X ¼ ½0; x3 sinΨj ; x3 cosΨj T . Due to the longitudinal correlation of the incoming turbulence, a certain periodicity is present as spectral peaks at Ω; 2Ω; 3Ω… in the frequency domain representation of the rotational sampled turbulence. This simple AR model used here does not represent the low-frequency, large-scale turbulent structures very well, given the homogeneity and isotropy assumption. Furthermore, the dynamics of the tower are more related to the frequency component of turbulence in the vicinity of the tower frequency. In this paper, the rotational sampled effect seems to be more important and is well accounted for by the present turbulence model (Zhang et al., 2014). Here the quasi-steady aerodynamic is employed which means the boundary layers on the profile are assumed to be fully attached and a change of the angle of attack is causing a comparable change of the aerodynamic loads without time delay. According to the criteria proposed by Leishman (2006), if the reduced frequency k < 0:05 the problem
V1 ðx3 ; tÞ ¼ ð1 aÞV0 þ v1 ðx3 ; tÞ u_1 ðx3 ; tÞ V2 ðx3 ; tÞ ¼ ð1 þ a0 ÞΩx3 þ v2 ðx3 ; tÞ u_2 ðx3 ; tÞ
0.7 1
0.0040.03 0.0010.01
(8)
(9)
where aðx3 Þ and a0 ðx3 Þ indicate the axial and tangential induction coefficients, determined by Blade Element Momentum (BEM) (Hansen, 2008). vj ðx3 ; tÞ indicate the components in the moving frame of reference of the rotational sampled turbulence at the position x3 , at the time t. u_j ðx3 ; tÞ denote the corresponding moving frame components of the velocity vector of the blade at the abscissa x3 , due to elastic deformations of the tower and the blade. The instantaneous angle of attack is given as:
αðx3 ; tÞ ¼ φðx3 ; tÞ þ βðtÞ þ δðx3 Þ
(10)
where φðx3 ; tÞ indicates the flow angle defined as, see Fig. 3:
tanφðx3 ; tÞ ¼
Table 1 Typical reduced frequencies for wind turbine operations. k
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi V 21 ðx3 ; tÞ þ V 22 ðx3 ; tÞ
Fig. 3 shows the wind velocities in the rotor plane. V1 ðx3 ; tÞ and V2 ðx3 ; tÞ are the instantaneous wind velocities in local x1 -direction and x2 -direction, respectively, seen by an observer fixed to the moving (x1 , x2 , x3 )-coordinate system. These are given as:
where v1;j ðx3 ; tÞ; v2;j ðx3 ; tÞ; v3;j ðx3 ; tÞ are rotational sampled turbulence components at the position x3 in jth blade, in the moving frames of reference. v1;j ðX; tÞ; v2;j ðX; tÞ; v3;j ðX; tÞ are rotational sampled turbulence components at the same position for blade j in the fixed frame of refer-
x3 =L
(7)
V1 ðx3 ; tÞ V2 ðx3 ; tÞ
(11)
The lift coefficient cl ðαÞ and drag coefficient cd ðαÞ are calculated based on the instantaneous angle of attack together with 2D profile data, where the delay effect on the loads due to change of α has not been considered. 175
B. Chen et al.
Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 171–179
The damping forces of the eddy current damper can be presented by:
The aerodynamic force along the blade is calculated through the BEM method with Prandtls tip loss factor and Glauert correction (Hansen, 2008). Non-linear aeroelasticity is also considered by introducing the local deformation velocities of the blade into the calculation of the flow angle and the angle of attack, as expressed in Eq. (9). As a result, high aerodynamic damping is introduced in the blade flap-wise and the fore-aft tower vibrations, but relatively low aerodynamic damping in the blade edgewise and the lateral tower vibrations (Chen et al., 2017b).
where q_ r ðtÞ denotes the torsional degree of freedom velocity vector of the section i and section j, Cr denotes the eddy current damping matrix in the torsional degree of freedom of these two nodes. cr denotes the equivalent damping coefficient.
3. Viscous damper design
Cr ¼ c r
f r ðtÞ ¼ Cr q_ r ðtÞ
Fig. 6 shows schematic diagram of the proposed eddy current damper, which is placed at the inner part of the blade. It is well known that when a conductive material is subjected to a time-varying magnetic flux, eddy currents are generated in the conductor. These eddy currents circulate inside the conductor generating a magnetic field of opposite polarity as the applied magnetic field (Sodano and Bae, 2004). The interaction of the two magnetic fields causes a torque that resists the change in magnetic flux. The torque is proportional to the slip r_6;j r_6;i , i.e. the difference between the angular velocity of the adjacent shafts, corresponding to a linear viscous damper with the damper constant cr . The value of the damper constant can easily be changed by adjusting the magnetic field, making an optimal semi-active control to changing wind load conditions possible. The damper device can be separated into two parts: the magnetic part and the conducting part. The magnetic part is attached to the inner shaft of the length a fixed at node i. The conducting part consist of two parts: the conducting tube and the transmission shaft of the length b fixed at node j. The damping torque to be transmitted to node j is much smaller than the elastic torque transmitted by the blade structure. Hence, the dimension of the shaft is correspondingly small. In order to make maintenance of the damper device, the support node i should be placed close to the nacelle. The length b of the shaft corresponding to the support node j is determined, so the performance of the damper is optimized. If the damper is not activated corresponding to cr ¼ 0, the blade renders into flutter in the fundamental torsional mode. At the other extreme corresponding to cr ¼ ∞, the rotations are locked at node j, and no mechanical energy is extracted from the blade. In this case the blade renders into flutter with a mode fixed at node j and with a critical rotational speed Ωcr of the rotor somewhat larger than for cr ¼ 0. It turns out that an optimal value of cr exists, where a maximum value of Ωcr is achieved and a maximum. The optimal value of cr depends both on the position of node j as indicated by the parameter b and on the turbulence in the incoming flow to the rotor as defined by the mean wind velocity V0 and the turbulence intensity I. V0 will change the pitch angle βðtÞ, and hence may influence the flutter critical rotational speed of the rotor.
1 1 ; 1 1
(12)
q_ r ðtÞ ¼
r_6;i r_6;j
(13)
4. Numerical example Here the 907-DOF model has been used for analyzing flutter stability. According to Hansen's previous study (Hansen, 2002), the dynamics of supporting structure may affect the aerodynamic damping of blade and also the risk of flutter. However, results from his later researches (Hansen, 2004a, 2004b) shown that the single-blade model and full-wind-turbine model analyses for a particular turbine predicts similar flutter limits, suggesting that the single-blade flutter analysis is sufficient for predicting the flutter limits. It was also demonstrated that the frequency characteristics of blade are less affected by the supporting structure dynamics. Nevertheless, we used this full-wind-turbine model in the present study, partly because this model has been fully developed and widely used in our previous works (Zhang, 2015), and partly because we want to mimic the real situation. 4.1. Parameters of DTU 10 MW wind turbine Data from the DTU 10 MW reference wind turbine (Bak et al., 2013) have been used to calibrate the structural model. The constant parameters employed in the 907-DOF wind turbine model are calculated and provided in Table 2. Table 2 Parameters in the 907-DOF wind turbine model. Parameter
Value
Unit
Parameter
Value
Unit
L s h Jr Jg kr kg
89.166 5 119 1.5696 108 2328.24 3.48 109 1.392 106
m m m kgm2 kgm2 Nm/rad Nm/rad
Ω g a θmax ρ
1.0053 9.81 4.815 0.04 1.0 103
rad/s m/s2 m rad kg/m3
Fig. 6. The eddy current damper.
176
B. Chen et al.
Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 171–179
Table 3 Natural frequencies of the 907-DOF wind turbine model. Item
Value
Unit
Description
ωf 1 ωf 2 ωf 3 ωe1 ωe2 ωe3 ωθ1 ωfa ωl
4.15 11.56 23.86 6.09 18.50 39.18 37.77 2.75 2.75
rad/s rad/s rad/s rad/s rad/s rad/s rad/s rad/s rad/s
1st flapwise frequency of blade 2nd flapwise frequency of blade 3rd flapwise frequency of blade 1st edgewise frequency of blade 2nd edgewise frequency of blade 3rd edgewise frequency of blade 1st torsional frequency of blade 1st tower fore-aft frequency 1st tower lateral frequency
The main natural frequencies for the 907-DOF wind turbine model, relevant for the analysis below, are listed in Table 3. In this model, the 1st flapwise frequency, 1st edgewise frequency and the 1st torsional frequency have been chosen to generate the structural damping matrix Cst based on the Cauchey's damping model (Vatne, 2011; Bak et al., 2013). Furthermore, we also assumed that the structural damping ratios corresponding to these three modes (the 1st flapwise mode, the 1st edgewise mode and the 1st torsional mode) are the same. Since the aeroelastic mode and structural mode of wind turbine blades are different when they operate on the rotor, as Campbell diagram shows (Hansen, 2004a), the blades of wind turbine do not vibrate in pure blade mode, i.e. flapwise, edgewise and torsional modes may couple in the rotor whirling modes of the wind turbine (Hansen, 2003). The aeroelastic modes coupling through the aerodynamic forces may lead to strong interactions, such as classic flutter. In this paper, the numerical simulation result shows it is the third flapwise bending mode that couples with the first torsional mode in flutter. Actually, according to flutter criteria in classic flutter theory (Hansen, 2007), the natural frequencies of a torsional mode and a flapwise bending mode must be relatively close for them to couple in a flutter mode. In the present case, the first torsional mode frequency is indeed very close to the third flapwise bending mode with a frequency ratio around 1.58. In the following sections, the effects of different parameters on the optimal value of the flutter critical rotational speed Ωcr are discussed in detail. These parameters include the level of structural damping ratio ζ0 , the level of the turbulence as measured by the turbulence intensity I, the transmission shaft length of the damper b, and the mean wind velocity V0 . Since the main load components of wind turbine blades comes from the rotational speed of the rotor and the turbulence component is numerically too small in comparison with other load components, further discussion about turbulence effect on flutter limits is not elaborated in this paper but will be the topic of follow-on studies.
Fig. 7. Time series of torsional angle of the blade tip at ccstr ¼ 0, V0 ¼ 15m=s, Lb ¼ 0:4268, ζ 0 ¼ 0:005, I ¼ 0:1.
Fig. 8. Relative flutter critical rotational speed as a function of the relative damping ratio at V0 ¼ 15m=s, Lb ¼ 0:4268. ____ , turbulence intensity I ¼ 0; , turbulence intensity I ¼ 0:1; ⋯, turbulence intensity I ¼ 0:3.
deviations from the mean value θmax ¼ μθ þ 3 σ θ has been chosen as a compromised to limit the probability of false alarms due to elastic exceedances of the threshold level, and to prevent structural damages at the onset of flutter (Chen et al., 2017a). Fig. 8 shows the influence of the flutter critical rotational speed Ωcr as a function of cr for different values of turbulence intensity I and structural damping ratio ζ0 , with V0 ¼ 15m=s and b ¼ 0:4268L. cr has been normalized with respected to cst , where cst indicates the component at the main diagonal of the structural damping matrix Cst at node j in the rotational degree of freedom. Ωcr has been nondimensionalized with respect to Ω. As seen, each curve has a peak at a specific relative damping ratio, which represents the maximum value of the critical rotational speed for flutter as discussed above. Furthermore, the flutter critical rotational speed is only marginally dependent on the turbulence intensity and will decrease with the increase of turbulence intensity. With the increasing value of structural damping ratio ζ0 , the peak will move in the direction of lower relative damping ratio. Fig. 9 shows the variation of Ωcr as a function of Lb for fixed value of ccstr , V0 , ζ0 and I. As can be seen, there is an optimal value of Lb, which maxi-
4.2. Result and discussion Here the rotor was purposefully operated at ever increasing rpm until the flutter boundary was reached and dramatic classical flutter oscillations were observed. Further, the simulation time at each increased rotational speed of rotor is long enough to obtain the steady-state vibration of blade tip. Fig. 7 shows the time series of the torsional degree of the blade tip under flutter critical conditions for fixed values of V0 , Lb, ccstr , ζ 0 and I. For convenience, the time history prior to that shown has been omitted, and only the initial segment of the unstable behaviour is displayed where flutter is evident. As can be viewed, as the envelope of the torsional angle exceeds the critical value θmax , the vibration starts to increase beyond limits, indicating instability. Here, the threshold value θmax has been set as the mean value of the time series of torsional angle of the blade tip plus three standard deviations under the stability situation for the following reason. Let μθ and σ θ indicate the mean value and the standard deviation of elastic deformations of the torsional rotation θðtÞ, caused by the turbulence alone. The indicated statistical moments and determined by ergodic sampling of a suitable line series before the onset of flutter. The value of 3 standard
mizes the value of
177
Ωcr Ω
(i.e. Ωcr ¼ 3.6 Ω), and it is almost 2.3 times the
B. Chen et al.
Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 171–179
5. Concluding remarks
Fig. 9. Variation of the flutter critical rotational speed as a function of V0 ¼ 15m=s, ζ 0 ¼ 0:005, I ¼ 0:1.
b L
at
cr cst
A new type of eddy current based ’electromagnetic’ linear viscous damper attached to the inner part of wind turbine blades has been proposed in the paper for enhancing the flutter stability. An investigation is carried out to evaluate the effectiveness of the proposed viscous damper in increasing the flutter critical rotational speed through a sophisticated 907-DOF aeroelastic model. The influence of the modal damping ratio of the structure, the turbulence intensity, the mean wind velocity and the transmission shaft length of the damper device on the flutter critical rotational speed have been discussed. It is concluded that if the transmission shaft length of the damper device b is almost half of the length of the blade L, the linear viscous damper at optimal tuning can increase the flutter critical rotational speed more than 100%. Too small or too large values of the equivalent damping coefficient can decrease the flutter critical rotational speed compared to the optimal tuned value. It is demonstrated that the mean wind velocity, the turbulence intensity and the structural damping of the torsional mode particularly in the flutter mode have little effect on the achieved enhanced flutter stability. All these results show that the proposed viscous damper is a promising damping device with the possibility of semi-active mode of operation and can be utilized in wind turbine blades for flutter control.
¼ 9,
Acknowledgements The first author gratefully acknowledges the financial support from the State's Key Project of Research and Development Plan (No. 2016YFE0127900) and the National Science Foundation of China (No. 51422806, No. 51278189). References Bak, C., Zahle, F., Bitsche, R., Kim, T., Yde, A., Henriksen, L., Natarajan, A., Hansen, M., 2013. Description of the DTU 10 MW Reference Wind Turbine. DTU Wind Energy Report-I-0092. C. Bak, F. Zahle, R. Bitsche, T. Kim, A. Yde, L. Henriksen, P. B. Andersen, A. Natarajan, M. H. Hansen, Design and performance of a 10 MW wind turbine, Wind Energy 124. Batchelor, G.K., 1953. The Theory of Homogeneous Turbulence. Cambridge university press. Belamadi, R., Djemili, A., Ilinca, A., Mdouki, R., 2016. Aerodynamic performance analysis of slotted airfoils for application to wind turbine blades. J. Wind Eng. Ind. Aerod. 151, 79–99. Bernhammer, L.O., van Kuik, G.A., De Breuker, R., 2016. Fatigue and extreme load reduction of wind turbine components using smart rotors. J. Wind Eng. Ind. Aerod. 154, 84–95. Bir, G., 2008. Multiblade coordinate transformation and its application to wind turbine analysis. In: ASME Wind Energy Symposium, pp. 1–15. Bir, G., Jonkman, J., 2007. Aeroelastic instabilities of large offshore and onshore wind turbines. J. Phys. Conf. Ser. 75, 012069. IOP Publishing. Bisplinghoff, R.L., Ashley, H., Halfman, R.L., 2013. Aeroelasticity. Courier Corporation. Bottasso, C., Campagnolo, F., Croce, A., Tibaldi, C., 2013. Optimization-based study of bend–twist coupled rotor blades for passive and integrated passive/active load alleviation. Wind Energy 16 (8), 1149–1166. Bottasso, C.L., Croce, A., Gualdoni, F., Montinari, P., 2016. Load mitigation for wind turbines by a passive aeroelastic device. J. Wind Eng. Ind. Aerod. 148, 57–69. Chen, B., Hua, X.G., Zhang, Z., Basu, B., Nielsen, S.R.K., 2017. Monitoring of wind turbine blades for flutter instability. Struct. Monitor. Mainten. 4 (2), 115–131. Chen, B., Zhang, Z., Hua, X., Basu, B., Nielsen, S.R.K., 2017. Identification of aerodynamic damping in wind turbines using time-frequency analysis. Mech. Syst. Signal Process. 91, 198–214. I. E. Commission, et al., Wind Turbine Part 1: Design Requirements, IEC 61400–61401, International Electrotechnical Commission, Geneva, Switzerland. De Goeij, W., Van Tooren, M., Beukers, A., 1999. Implementation of bending-torsion coupling in the design of a wind-turbine rotor-blade. Appl. Energy 63 (3), 191–207. Dowell, E.H., Curtiss, H.C., Scanlan, R.H., Sisto, F., 1989. A Modern Course in Aeroelasticity, vol. 3. Springer. Fung, Y.C., 2002. An Introduction to the Theory of Aeroelasticity. Courier Corporation. Fung, Y.C., 2004. An Introduction to the Theory of Aeroelasticity. Dower. D. T. Griffith, T. D. Ashwill, reportThe Sandia 100-meter All-glass Baseline Wind Turbine Blade: SNL100-00, Sandia National Laboratories, Albuquerque, Report No. SAND2011–3779. Hansen, M., 2002. Vibrations of a three-bladed wind turbine rotor due to classical flutter. In: 40th AIAA Aerospace Sciences Meeting and Exhibit, 2002 ASME Wind Energy Symposium.
Fig. 10. Variation of the flutter critical rotational speed with the mean wind velocity V0 at cr b cst ¼ 9, L ¼ 0:4268, ζ 0 ¼ 0:005, I ¼ 0:1.
value, without the damper devise (i.e. Ωcr ¼ 1.6 Ω). The maximum value for b ¼ 0:4268L is Ωcr ≃ 3:3Ω. As seen on Fig. 9, slight increase to Ωcr ≃ 3:6Ω can be achieved by fixing the damper device at b ≃ 0:58L. However, this increase is at the expense of a longer torque shaft to the damper device. Fig. 10 shows the variation of Ωcr as a function of the mean wind velocity V0 for fixed values of ccstr , Lb, ζ0 and I. Evidently, Ωcr is found to decrease monotonically for increasing mean wind velocity V0 . However, Ωcr only decreases slightly with V0 , which is similar to Politakis et al.’s result (Politakis et al., 2008). Furthermore, as seen from Eq. (8) and (9), the resulting wind velocity V will increases with V0 , which in turn will increases the proneness to instability. In conclusion, Figs. 8 and 10 show that the optimal damper constant cr for optimal flutter stability only depends slightly on the incoming wind field to the rotor as characterized by the mean wind velocity V0 and the turbulence intensity I. 178
B. Chen et al.
Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 171–179 Moshfeghi, M., Shams, S., Hur, N., 2017. Aerodynamic performance enhancement analysis of horizontal axis wind turbines using a passive flow control method via split blade. J. Wind Eng. Ind. Aerod. 167, 148–159. Nayfeh, A.H., Mook, D.T., 1995. Nonlinear Oscillations. John Wiley and Sons. Ogata, K., 2010. Modern Control Engineering, fifth ed. Pearson. Pars, L.A., 1979. A Treatise on Analytical Dynamics. OX Bow Press. Politakis, G., Haans, W., Van Bussel, G., 2008. Suppression of Classical Flutter Using a ’Smart Blade’. American Institute of Aeronautics and Astronautics AIAA. Rafiee, R., Tahani, M., Moradi, M., 2016. Simulation of aeroelastic behavior in a composite wind turbine blade. J. Wind Eng. Ind. Aerod. 151, 60–69. Resor, B., Paquette, J., 2011. Uncertainties in prediction of wind turbine blade flutter. In: 49 th AIAA Aerospace Sciences Meeting and Exhibit. Resor, B.R., Owens, B.C., Griffith, D.T., 2012. Aeroelastic instability of very large wind turbine blades. In: European Wind Energy Conference Annual Event, Copenhagen, Denmark. Sodano, H.A., Bae, J.-S., 2004. Eddy current damping in structures. Shock Vib. Digest 36 (6), 469. Theodorsen, T., Mutchler, W., 1935. General Theory of Aerodynamic Instability and the Mechanism of Flutter. National Advisory Committee for Aeronautics, Washington, DC, USA. Vatne, S.R., 2011. Aeroelastic Instability and Flutter for a 10 MW Wind Turbine. Master’s thesis. Institutt for energi-og prosessteknikk. Veers, P., Bir, G., Lobitz, D., 1998. Aeroelastic tailoring in wind-turbine blade applications. In: Proceedings, Windpower, 98, pp. 291–304. Zhang, Z., 2015. Passive and Active Vibration Control of Renewable Energy Structures. Ph.D. thesis. Aalborg University, The Faculty of Engineering and Science, Division of Reliability, Dynamics and Marine Engineering. Zhang, Z., Nielsen, S.R.K., Blaabjerg, F., Zhou, D., 2014. Dynamics and control of lateral tower vibrations in offshore wind turbines by means of active generator torque. Energies 7 (11), 7746–7772.
Hansen, M., 2003. Improved modal dynamics of wind turbines to avoid stall-induced vibrations. Wind Energy 6 (2), 179–195. Hansen, M., 2004. Stability analysis of three-bladed turbines using an eigenvalue approach. In: 2004 ASME Wind Energy Symposium. American Institute of Aeronautics and Astronautics; American Society of Mechanical Engineers, pp. 192–202. Hansen, M.H., 2004. Aeroelastic stability analysis for for wind turbines using an eigen value approach. Wind Energy 7 (2), 133–143. Hansen, M.H., 2007. Aeroelastic instability problems for wind turbines. Wind Energy 10 (6), 551–577. Hansen, M.O.L., 2008. Aerodynamics of Wind Turbines, second ed. Earthscan. Hayat, K., Ha, S.K., 2015. Load mitigation of wind turbine blade by aeroelastic tailoring via unbalanced laminates composites. Compos. Struct. 128, 122–133. Hayat, K., de Lecea, A.G.M., Moriones, C.D., Ha, S.K., 2016. Flutter performance of bend–twist coupled large-scale wind turbine blades. J. Sound Vib. 370, 149–162. Jonkman, B.J., 2009. Turbsim User's Guide: Version 1.50. National Renewable Energy Laboratory Colorado. Krenk, S., Svendsen, M.N., Høgsberg, J., 2012. Resonant vibration control of three-bladed wind turbine rotors. AIAA J. 50 (1), 148–161. Larsen, J.W., Nielsen, S.R.K., Krenk, S., 2007. Dynamic stall model for wind turbine airfoils. J. Fluid Struct. 23 (7), 959–982. Leishman, G.J., 2006. Principles of Helicopter Aerodynamics. Cambridge university press. Lobitz, D.W., 2004. Aeroelastic stability predictions for a MW-sized blade. Wind Energy 7 (3), 211–224. Lobitz, D.W., 2005. Parameter sensitivities affecting the flutter speed of a MW-sized blade. J. Sol. Energy Eng. 127 (4), 538–543. Lobitz, D.W., Veers, P.S., 1998. Aeroelastic behavior of twist-coupled HAWT blades. In: ASME/AIAA Wind Energy Symposium, Reno, NV, pp. 75–83.
179