- Email: [email protected]
Efficient Nonlinear Wind-Turbine Modeling For Control Applications
Efficient Nonlinear Wind-Turbine Modeling For Control Applications Morten D. Pedersen ∗,†,§ , Thor I. Fossen ∗,‡ ∗ Department of Engineering Cybernetics, Norwegian University of Science and Technology (NTNU), O.S.Bragstadspl.2D, Trondheim NO-7491, Norway † Nowitech - Norwegian Research Centre for Offshore Wind Technology ‡ Centre for Ships and Ocean Structures (CeSOS), Norwegian University of Science and Technology, Trondheim NO-7491, Norway § Corresponding Author; [email protected]
Abstract: A simple nonlinear model for variable pitch wind turbines is described. The model is aimed specifically at nonlinear control design and analysis. We propose a parametric model and subsequently describe a method to determine the parameters for any turbine that falls into the scope of the model. The result is a model that describes a turbine rotor, which is free to move in all six degrees of freedom. This freedom is required to model the new breed of floating offshore turbines. The resulting model is compared to experimental data with good results. Cyclic forces at large yaw offsets are not modeled well, however. Limitations of the model are discussed. Keywords: Wind turbine; control; dynamic model; identification. 1. INTRODUCTION The controlled variables in a VPVS turbine are the blade pitch angles and the generator torque. A wind turbine is limited in its power output by the so-called rated power. This threshold can vary from a few kW up to several MW in the largest turbines. Structural and electromechanical considerations prohibit any additional power extraction beyond this point. The controller in a wind turbine usually acts to optimize the power absorption below rated power, and to limit the power absorption above. This is achieved by a combination of torque control and pitch control. Active load alleviation may also be applied through control. The control algorithms for modern variable-pitch wind turbines, that one may find in the literature, are usually based on some type of simplified wind turbine model. These may be in the form of lookup-tables as in [Bianchi et al., 2006, MihetPopa et al., 2004], or linear models obtained from complex numerical simulation tools [Namik and Stol, 2010, Stol and Balas, 2003]. Hybrid models blending lookup tables with mechanical models have also been used [Bottasso et al., 2007]. These and even simpler approaches predominate. However, none of these approaches yield real insight into the behavior of the system. Linear models are only valid in a small envelope around the linearization point, which requires several individual models to cover the operational domain of the turbine. On the other hand, lookup-tables are not able to reproduce the transient effects due to the turbine wake; a highly dynamic structure. The use of lookup-tables is furthermore quite inconvenient when considering effects other than those that occur in axial flow conditions; the required data-set will be too large for practical purposes. Wind turbines are highly nonlinear machines, and their response to wind and control inputs vary dramatically over their
operational domain. With the knowledge that a wind turbine is a highly non-linear system, it would seem useful to have a model that describes the turbine in this context. Nonlinear model-based control usually require a simple but accurate model. This work aims to produce such a tool. A simple nonlinear model will allow researchers to leverage the large field of nonlinear control theory towards wind-turbine applications. Improved control allows optimal use of the existing generating machinery, as well as better utilization of the local wind resource. In addition, control may extend the service-life of turbines by reducing unwanted motion. This work focuses on the turbine rotor, as good models for substructures and electrical machinery abound in the literature. 2. MODELING THE WIND TURBINE ROTOR The turbine rotor is modeled as a flat disc swept by N , uniformly spaced, rigid blades. In order to yield a tractable analytical model, the following simplifications are made: • The rotor disk is modeled as one rigid structure, that is, we do not consider blade flexibility or root hinges. • The blades have their center of mass and center of lift on the pitch axis. • We assume that the lift curve slope of the blade sections is linear throughout. • Drag is significant only in the axial degrees of freedom, surge and shaft rotation. The modeling of a wind turbine rotor may be summarized eloquently by Figure 1. The aim of our rotor modeling is then to specify the contents of each of the three components in the diagram with an appropriate theory. The dynamic equations of the rotor are given a robot-like vectorial form. This type of model has seen great success in the field of marine control, [Fossen, 2011]. The compact form and physical insight that
Inflow
Wake model
ωz
Circulations
w Lifting Theory
Flow Environment
zbi
xbi
Blade Motion
Rotor Dynamics
ybi φ + δi
Forces and Moments
zr
Ω Controls
u
xr yr
θ
Fig. 1. Components of the rotor model, adopted from [Peters et al., 1989].
v ωy
r
such a model can provide seems very useful for wind-turbine modeling.
pIr R zI
2.1 Turbine Geometry and Kinematics
xI
Consider the rotor depicted in Figure 2. Three reference frames are used: • I, Inertial frame • r, Rotor frame • bi , Blade frames, N in all.
The blade frames and the rotor frame have coincident origins. Each blade, indexed with i, is associated with a unique frame, bi . We let ra ∈ R3 denote a coordinate expressed in frame a, and ei ∈ R3 , denote a unit vector. The cross-product operator is defined as: " # " # a × b := S(a)b,
S(a) =
0 −a3 a2 a3 0 −a1 −a2 a1 0
,
a=
a1 a2 a3
The model is developed in the rotor frame. The rotor frame is nearly stationary under normal operating conditions. Furthermore, this frame relates readily to motions and forces on a support structure. A coordinate in the rotor frame relates to the inertial frame via the transformation: rI = Rr rr Where: Rr (ψ, θ) =
cos(ψ) − sin(ψ) 0 sin(ψ) cos(ψ) 0 0 0 1
·
cos(θ) 0 sin(θ) 0 1 0 − sin(θ) 0 cos(θ)
The blades are spaced by an angle δi = 2π(i − 1)/N for i = 1 . . . N . The individual blade frames are arrived at through the transformation: rri = Rb,i rb with: Rb,i (φ + δi ) =
1 0 0 0 cos(φ + δi ) − sin(φ + δi ) 0 sin(φ + δi ) cos(φ + δi )
The blade pitching axes are parameterized in the longitudinal coordinate r: T rb = [0, 0, r] = rez The matrix: Ri := Rr Rb,i describes the attitude of a turbine blade using Euler angles in a yaw-pitch-roll sequence. The position of a point on the blade longitudinal axis, in the inertial frame, may now written as: rIi = Ri ez r + pIr where the rotor frame is placed at pIr .
yI ψ
Fig. 2. Rotor kinematics. The blade and rotor frames have been separated for clarity. The angle φ describes the rotor azimuth angle. The angles θ and ψ describe the platform pitch and yaw respectively. Quasi-velocities The motion of the rotor frame is expressed in terms of quasi-velocities. In our case, these are the translational velocity of the rotor as well as it’s rotational rate, as seen from the rotor frame. Define: h i ω
r
:=
Ω ωy ωz
,
v
r
T
I
:= Rr p ˙r =
u v w
r
,
ξ :=
v r ω
˙ See Figure 2 for a graphical representation. The where Ω = φ. vector ξ contains the quasi-velocities used to describe platform motion. A Jacobian may be used to relate this vector to a point velocity expressed in any frame, [Spong et al., 2006]: va = Ja ξ (1) Furthermore, Jacobians have the useful property that point forces expressed in any coordinate system may be mapped back to the generalized forces: τ = JTa f a (2) A vector of generalized forces in 6 DOF will be denoted τ ∈ R6 throughout the text. The velocity of the blade axes, as expressed in the blade and inertial frames, will be required. The relevant Jacobians are stated as: viI = Ji ξ,
vib,i = Li ξ
Ji = Rr [I, −S(Rb,i ez r)] ,
(3) Li =
RTb,i
[I, −S(Rb,i ez r)]
2.2 Aerodynamic Model The flow on the rotor is composed of: • External flow; wind and wind shear. • Wake induced flow. • Flow from motion of the rotor. Lift Forces Let w(r) = [wx , wy , wz ]T and f (r) = T [fx , fy , fz ] denote the vector of apparent flow velocity and aerodynamic forces in the blade frame, respectively. These are functions of the radius r. The circulation on the airfoil
is lumped into a point vortex at the quarter-chord, see Figure 3. Assuming small angles, the sectional circulation may be computed as, [Katz and Plotkin, 2001]: Γ = πcwn ' πc(wx + wy (β + α)) (4) where c = c(r), β = β(r) and α are the blade chord, pitch and twist respectively. This approach requires that the blade section operates in the linear regime, i.e. not under stalled conditions at very high angles of attack. The three dimensional vortex lifting xb β+α
wn ' wx + wy β
Γ
wx
wn
yb wy
Fig. 3. Sectional airfoil. law, [Saffman, 1992]: df = ρΓS(w)dl = −ρΓS(ez )wdr (5) yields the differential force on a blade section. The apparent flow w is comprised of the flow due to the wind and the induced velocities, as well as the flow induced by own motion. The form of induced flow used in the Pitt & Peters model, eq.6, may conveniently be expressed in terms of rotor velocities. The mean inflow corresponds to a surge motion, while the cyclic modes correspond to rotor roll and pitch. The external flow may similarly be expressed in terms of rotor displacements. With eq.3, the flow on the blade section may be written: wi = Li −ζ − ξ + KT η where we have introduced the vectors ζ ∈ R6 , which represents the external flow and η ∈ R3 which represents the induced flow. The external flow is produced by the hub-height freestream velocity V∞ entering at an angle ψ∞ w.r.t. reference frame I: RTr (ψ − ψ∞ , θ)ex ζ= V RTr (ψ − ψ∞ , θ)ey (α0 /h) ∞ Wind shear is incorporated by a linearized power law model, taking the hub-height h, and shear exponent α0 as variables. See Fig.4. The constant matrix K ∈ M6×3 maps the induced flow onto the rotor: T
K
=
1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1
We collect these into the apparent flow vector: ν := −ζ − ξ + KT η, ν ∈ R6 This lets us write the differential force on blade i as: dfi = −Γi S(ez )Li νdr Individual Pitch The pitch angle α may vary between blades. A common way to express the individual pitch on the blades, is through cyclic pitch, Johnson [1994]. Here, the pitch on blade i is produced by: αi = [1, cos(φ + δi ), sin(φ + δi )] α, α := [α0 , αc , αs ]T The collective pitch is α0 , while αc and αs are the cyclic pitch. Generalized Aerodynamic Force By eq.2, we may express the distributed aerodynamic forces in terms of the quasicoordinates: dτiA = LTi dfi = −Γi LTi S(ez )Li νdr The generalized aerodynamic force is obtained by summing and integrating over the blades:
τ
A
X
−ρ
=
!
R N
Z 0
Γi L T i S(ez )Li dr
ν
i=1
= ρπA [A0 (ν, α) + cos(N φ)Ac (ν, α) + sin(N φ)As (ν, α)] ν
where the lift forces have been expressed in terms of three lift matrices, which are functions of the apparent flow ν and the individual pitch vector α. The forces are composed of a mean force as well as a cyclic component. The matrices are scaled by the rotor swept area A. The matrix A0 is conveniently the same for any blade count. The matrices Ac and As do not have this property, and must be computed for a specific blade count. The matrices are all skew symmetric: A0 = −AT0 , Ac = −ATc , As = −ATs The matrix elements are furthermore linear in the apparent flow. Note that the pitch control enters through these matrices. This is unconventional, but preserves the skew-symmetry. Model Parameters The lift matrices contain integrals in the rotor chord and twist. These are most conveniently expressed in terms of dimensionless parameters, σn and tn . For different powers of r, these appear as: Z R Z R crn dr =
N 0
π Rn+2 σn , n+1
βcrn dr =
N 0
π Rn+2 tn n+1
In terms of these parameters, the mean force matrix is: S(γ1 ) RG A0 = −RG R2 S(γ2 ) with: γ1 :=
γ2 :=
0 2ν3 (σ0 α0 + t0 ) + Rσ1 (ν6 − ν4 αs ) −2ν2 (σ0 α0 + t0 ) − Rσ1 (ν5 − ν4 αc )
1
4
1 24
0 4ν3 (σ2 α0 + t2 ) + 3Rσ3 (ν6 − ν4 αs ) −4ν2 (σ2 α0 + t2 ) − 3Rσ3 (ν5 − ν4 αc )
G :=
−
R 3
σ1 4
(σ2 α0 + t2 )ν4 −
σ1 2
ν1
1 0 1 0 2 0 0
0
0 1
2
ν2 αc + ν3 αs 0 0 0 (3ν2 αc + ν3 αs )/4 (ν3 αc + ν2 αs )/4 0 (ν3 αc + ν2 αs )/4 (ν2 αc + 3ν3 αs )/4
The matrices As and Ac do not stay similar for various blade counts, and are not presented here. Drag Forces The forces resulting from viscous drag do not have an elegant theory like lift forces. Furthermore, drag forces are usually significantly smaller than lift forces. A complete analytical model of the viscous drag forces on a turbine, in all degrees of freedom, would be prohibitively complex. Such a model would be counterproductive, as simplicity is the main aim of the present work. To avoid an overly complicated model the drag is only analyzed in the surge and roll degrees of freedom, considering purely axial flow. The drag in surge is modeled by a constant drag coefficient, whereas the drag coefficient in the rotational degree of freedom is taken as a curve-fit in the wind speed: 1 1 dx = − ρACD,x ν12 , dφ = − ρAR3 CD,φ (V∞ )Ω2 2 2 Pitt & Peters Wake Model The wake behind the rotor is produced as the rotor exerts a force on the air flowing through it. Wake models range from simple momentum theory, to advanced free vortex or CFD methods, [Bossanyi et al., 2000]. A very attractive wake model suitable for the present work is
found in [Pitt and Peters, 1981]. This model has seen extensive use in helicopter and wind-turbine engineering. Importantly, the approach allows the modeling of yaw misalignment w.r.t to the wind, a common occurrence in operation. The P&P model describes the inflow distribution as: η = η0 + r cos(φ)ηc + r sin(φ)ηs (6) The inflow is scaled by three inflow states, which is gathered in a vector η := [η0 , ηc , ηs ]T . See Figure 4 for a graphical representation. The evolution of these wake states is governed by a set of nonlinear differential equations: A ρARMw η˙ + ρAVw L−1 w (χ)η = −Kτ Here expressed in a dimensional form. We employ K to select the required axial force as well as pitch and roll moments from the generalized aerodynamic force, as required by the wake model. The mass matrix is augmented with an adjustment factor as suggested by Suzuki [2000]. The wake skew angle χ denotes the angle of the wake with respect to the rotor, see Figure 9.
M(φ) = N m
η0
Z
V∞
xI ψ∞
Y
ηs
xr
ψ ηc
yI
r CN (φ)S(ex Ω) − S(ω )CN (φ) 2 r
r
S(ωa )
cg r
cg CN (φ)S(ωa ) rg
ωar
i
S(ω )DN (φ) − DN (φ)S(ex Ω)
r
:= ω − Ωex , and: 0 − cos(φ) − sin(φ) , N =1 0 0 N cos(φ) X 0 sin(φ) 0 CN (φ) := S(Rb,i ez ) = 0 0 0 ∀ N ≥2 i=1 0 00 , 0 0 0 1 0 0 2 0 cos (φ) cos(φ) sin(φ) , N < 3 2 N 0 cos(φ)sin(φ) sin (φ) X 2 1 0 0 DN (φ) := − S (Rb,i ez ) = 1 0 , 0 N ≥3 i=1 0 20 1 where,
2
The integrated quantities are given by: Z R Z R ρ(r)dr, 0
1 cg := m
ρ(r)rdr,
s 1 m
rg :=
0
Z
R
ρ(r)r2 dr 0
Note that a blade count with N < 3 will inevitably result in a time varying system. Utilizing the above results we may write: M(φ)ξ˙ + C(φ, ξ)ξ = τe
X
yr
I −cg CN (φ) 2 cg CN (φ) rg DN (φ)
C(φ, ξ) = N m
m := χ ' ψ − ψ∞
h
Fig. 4. Inflow modes and wake skew angle in the Pitt & Peters model.
This form of the equations of motion avoids the need to consider each blade separately. Furthermore, the equations are in the rotor-frame, thus allowing easy integration with the aerodynamic model. The dynamic model have several notable properties. The mass matrix is symmetric. The coriolis matrix is T skew for N ≥ 2. It may be verified that h symmetric C = −C i T ˙ ξ M(φ) − 2C(φ, ξ) ξ = 0, which implies that the system is conservative, [Siciliano et al., 2009]. 2.4 Integrated Model The model equations may now be posed on the compact form, a set of first order nonlinear ODE’s:
2.3 Dynamic Model A simple dynamic model is derived to complement the aerodynamic model. The acceleration of a particle along the blade axis is: aIi = v˙ iI = J˙ i ξ + Ji ξ˙ The force balance of our particle along the longitudinal axis is: ρ(r)aIi dr = fiI dr where ρ(r) is the mass density along the blade and fi an external forcing. Using eq.2 one can write: ρ(r)JTi aIi dr = JTi fiI dr ⇒ ρ(r)JTi Ji drξ˙ + ρ(r)JTi J˙ i drξ = JTi fiI dr
It is necessary to sum up the contribution for the particles from all blades, N in all. A subsequent integration over r will yield the final dynamic equations. Define:
M(φ) :=
N Z X
N Z X i
T
ρ(r)Ji Ji dr,
C(φ, ξ) :=
0
i
τe :=
R
i
R
Ji
N Z X
T I fi dr
0
Integration and summation yields:
0
R T ˙ ρ(r)Ji J i dr
M(φ)ξ˙ + C(φ, ξ)ξ − ρπA [A0 + cos(N φ)Ac + sin(N φ)As ] ν = τg + τd + τ −1
RMw η˙ + Vw Lw η + πK [A0 + cos(N φ)Ac + sin(N φ)As ] ν = 0 T
ν = −ζ − ξ + K η
where τd denotes the forces due to drag. The rotor system will be connected to a substructure exerting a force on the rotor: τ . Gravity forces enter through τg , but are trivial, as a balanced flat rotor will act as point mass at the point of connection. This model bears great resemblance to vectorial models of marine craft, Fossen [2011]. A particular note of interest is the skew symmetry of the lift matrices, which has not been noted previously. Skew symmetry may play an important role in stability proofs for control systems, Siciliano et al. [2009]. 3. PARAMETER ESTIMATION The parameters in the lift matrices serve as good first estimates for the model. However, tuning these parameters against a more advanced approach, or empirical data, will yield a better model. To this end, each parameter is augmented with an adjustment factor, e.g. σ ˆ1 = σ1 . The drag model is also fitted. The Pitt & Peters model is based on an actuator disc model of the rotor, Pitt and Peters [1981]. An actuator disc does
not model the effects of a finite number of blades, rather, it considers the aerodynamic force to be smeared over an infinitely thin disc, Bossanyi et al. [2000]. The effects of having a finite number of blades is that the induced flow on the blades, will be larger than the disc average. This is often referred to as tip-loss. The disc averaged uniform induction, η0 , is used in the Pitt & Peters model. An adjustment factor is introduced to capture the effects of the tip loss, as well as non-uniform inflow, η˜0 = kη η0 , where η0 denotes the induction as it appears on the blades. The aim of the parameter estimation is to optimize the adjustment factors, fit the drag model, as well as setting kη . To achieve this, a numerical BEM tool is used to generate a large set of data points, covering the operational envelope of the turbine. BEM is the industry standard means of obtaining wind turbine performance estimates, and can be relied on for accurate data. Such a program yields the steady state axial forces and torques on the turbine, as it operates in axial flow. The steady axial subsystem of the present aerodynamic model is extracted; rows and columns 1 and 4 in the mean lift matrix A0 . The present aerodynamic model is unsteady, so the equilibrium points of the model are fitted against the data-set. The equilibrium of the model in axial flow is when: 1 4(1 − kη a)kη a = πλ [3(1 − a)σ1 − 2λ (t2 + α0 σ2 )] (7) 3 where dimensionless variables have been used: a = η0 /V∞ , λ = ΩR/V∞ Let a? denote the physical solution to eq. 7, (the other solution is discarded): q 12kη + 3πλσ1 −
?
a
2
9 4kη + πλσ1
=
Fig. 5. A typical behavior for pitch-controlled, variable-speed turbines. The pitch is activated as the maximum power is reached, to limit power, in this case 5MW. The tip speed ratio is decreased as the rotational rate is kept constant for high wind speeds. procedure will in fact work poorly if the data-set is confined to the reference trajectory. The model has two significant limitations; the P&P model will be unstable for large induction ratios, a > 0.5. Also, departures from a linear lift slope will violate the assumptions inherent in the lift modeling. To avoid fitting the model to operating conditions that lies beyond it’s scope, data points exceeding these limitations are removed. If these provisions are not taken, gross errors resulting from limitations in the model will drive the optimization procedure towards a less desirable result. The data points are subsequently weighed. The weighing is
2 (−3σ + 2λ (t + α σ )) + 48πλkη 1 2 0 2
2 24kη
Four dimensionless quantities, are used to fit the parameters in the model. CT (α0 , λ) = CP (α0 , λ) =
T 1 ρAV 2 ∞ 2
P 1 ρAV 3 ∞ 2
'
1 3
?
πλ 3(1 − a )σ1 − 2λ (t2 + α0 σ2 ) ?
' (1 − a )CT (α0 , λ)
CB (α0 , λ) =
N mB 1 ρARV 2 ∞ 2
CS (α0 , λ) =
N fS 1 ρAV 2 ∞ 2
'
1 6
?
πλ 4(1 − a )σ2 − 3λ (t3 + α0 σ3 ) ?
?
' π(1 − a ) 2(1 − a )σ0 − λ (t1 + α0 σ1 )
Here a? has been inserted. These are the thrust, power, out of plane root bending and in-plane root shear coefficients respectively. The parameter t0 does not appear in these coefficients, and is used without an adjustment. By inserting a? into the coefficients, one fits the equilibria in the dynamic model towards the data set. One has furthermore removed the need to fit the induction itself. The induction is simply an internal state of the model, enforcing correct steady state behavior. 3.1 Selection and Weighing of Data Set To produce a data set for fitting, a BEM code is used to generate a typical reference trajectory in terms of the pitch angle and tip-speed ratio, with respect to varying wind speed. See Fig.5 for an illustration. The data points to be fitted are generated by adding normally distributed random offsets in pitch-angle and tip-speed ratio from the reference trajectory. The fitting 1
Blade Element Momentum
Fig. 6. Left: Comparison of the wind speed probability function at the site (red), and in the data set (histogram). Right: Data set, color coded with respect to weighing. done so as to render the weighted likelihood of encountering a certain wind speed in the data set, similar to the probability at the turbine site. This probability may be assumed to be on the form of a Weibull curve, Bossanyi [2003]: V k−1 f (V∞ ) = k ∞k exp −(V∞ /c)k c Figure 6 illustrates the result of the weighing. The parameters in the model, are subsequently fitted by weighted least squares. This is achieved by the Gauss-Newton method. The optimum is found within 4-10 iterations. The lift and drag forces are fitted separately, and one must use a BEM code that gives separate coefficients for the potential and viscous effects. Fitting towards the sum of these yields poor results. The four coefficients must be fitted jointly, due the linking introduced by the equilibrium induction.
4. MODEL VALIDATION
Mz @kNmD 50
The fitted model replicates the equilibrium forces from the BEM code with good results. The model has been run with
100
150
200
250
300
350
f@degD
-20 -40 -60 -80 -100 -120 -140
Fig. 9. Yaw moment due to a yaw misalignment on the 2MW Tjæreborg turbine. The yaw offset is 32◦ , wind speed is 8.5[m/s]. Model is dashed.
Fig. 7. BEM results compared to fitted model for a 5MW turbine. Model data are dashed. the same controller that were used for the BEM. Note the very slight deviation in the pitch angle at high wind speeds, which is likely caused by violation of the small angle assumption in eq. 4. The point of optimum power extraction CP? = CP (α0? , λ? ) is also predicted very well. An accurate estimate of the power Model BEM
? CP
α?0
λ?
0.481 0.482
-0.219◦
7.520 7.418
-0.746◦
optimum is important if one desires to develop optimizing control methods. The dynamic wake implementation is also
Fig. 8. The nonlinear model compared to experiment and a lookup-table approach at wind speeds 12.5[m/s](top) and 7.4[m/s](bottom). The response to a pitch step is depicted. The experimental data are from the 2MW Tjæreborg turbine. quite successful as seen in Figure 8 2 . Note that the proposed model captures the transients in the experiment, whereas the lookup approach fails in this respect. A final test is the response to a large yaw misalignment. Although the model yields a reasonable mean value of the resulting yaw moment Mz , the cyclic component is very much underestimated. The reason for this failure is unknown at present. A likely cause is the very 2
The Authors thanks Martin O.L. Hansen for providing experimental data.
simple inflow model being unable to capture the complicated unsteady flow field in a yawed flow. Under-prediction of cyclic forces is indeed present in many inflow models, Bossanyi [2003]. It may be better to retain only the steady component in simple models such as the present one. REFERENCES F. Bianchi, H. De Battista, and R. Mantz. Wind Turbine Control Systems: Principles, Modelling and Gain-scheduling Design (Advances in Industrial Control). 2006. E. Bossanyi, T. Burton, D. Sharpe, and N. Jenkins. Wind energy handbook. John W iley, 2000. EA Bossanyi. Wind turbine control for load reduction. Wind Energy, 6(3):229–244, 2003. ISSN 1099-1824. CL Bottasso, A. Croce, and B. Savini. Performance comparison of control schemes for variable-speed wind turbines. In Journal of Physics: Conference Series, volume 75, page 012079. IOP Publishing, 2007. T.I. Fossen. Handbook of Marine Craft Hydrodynamics and Motion Control. Wiley, 2011. W. Johnson. Helicopter theory. Dover Pubns, 1994. J. Katz and A. Plotkin. Low-speed aerodynamics, volume 13. Cambridge University Press, 2001. L. Mihet-Popa, F. Blaabjerg, and I. Boldea. Wind turbine generator modeling and simulation where rotational speed is the controlled variable. Industry Applications, IEEE Transactions on, 40(1):3–10, 2004. H. Namik and K. Stol. Individual blade pitch control of floating offshore wind turbines. Wind Energy, 13(1):74–85, 2010. ISSN 1099-1824. D.A. Peters, D.D. Boyd, C.J. He, et al. Finite-State InducedFlow Model for Rotors in Hover and Forward Flight. Journal of the American Helicopter Society, 34:5, 1989. D.M. Pitt and D.A. Peters. Theoretical prediction of dynamicinflow derivatives. Vertica, 5(1):21–34, 1981. PG Saffman. Vortex force and bound vorticity. Vortex Dynamics, 1992. B. Siciliano, L. Sciavicco, and L. Villani. Robotics: modelling, planning and control. Springer Verlag, 2009. M.W. Spong, S. Hutchinson, and M. Vidyasagar. Robot modeling and control. John Wiley & Sons Hoboken, NJ, 2006. K.A. Stol and M.J. Balas. Periodic disturbance accommodating control for blade load mitigation in wind turbines. Journal of solar energy engineering, 125:379, 2003. A. Suzuki. Application of dynamic inflow theory to wind turbine rotors. 2000.
Abstract: A simple nonlinear model for variable pitch wind turbines is described. The model is aimed specifically at nonlinear control design and analysis. We propose a parametric model and subsequently describe a method to determine the parameters for any turbine that falls into the scope of the model. The result is a model that describes a turbine rotor, which is free to move in all six degrees of freedom. This freedom is required to model the new breed of floating offshore turbines. The resulting model is compared to experimental data with good results. Cyclic forces at large yaw offsets are not modeled well, however. Limitations of the model are discussed. Keywords: Wind turbine; control; dynamic model; identification. 1. INTRODUCTION The controlled variables in a VPVS turbine are the blade pitch angles and the generator torque. A wind turbine is limited in its power output by the so-called rated power. This threshold can vary from a few kW up to several MW in the largest turbines. Structural and electromechanical considerations prohibit any additional power extraction beyond this point. The controller in a wind turbine usually acts to optimize the power absorption below rated power, and to limit the power absorption above. This is achieved by a combination of torque control and pitch control. Active load alleviation may also be applied through control. The control algorithms for modern variable-pitch wind turbines, that one may find in the literature, are usually based on some type of simplified wind turbine model. These may be in the form of lookup-tables as in [Bianchi et al., 2006, MihetPopa et al., 2004], or linear models obtained from complex numerical simulation tools [Namik and Stol, 2010, Stol and Balas, 2003]. Hybrid models blending lookup tables with mechanical models have also been used [Bottasso et al., 2007]. These and even simpler approaches predominate. However, none of these approaches yield real insight into the behavior of the system. Linear models are only valid in a small envelope around the linearization point, which requires several individual models to cover the operational domain of the turbine. On the other hand, lookup-tables are not able to reproduce the transient effects due to the turbine wake; a highly dynamic structure. The use of lookup-tables is furthermore quite inconvenient when considering effects other than those that occur in axial flow conditions; the required data-set will be too large for practical purposes. Wind turbines are highly nonlinear machines, and their response to wind and control inputs vary dramatically over their
operational domain. With the knowledge that a wind turbine is a highly non-linear system, it would seem useful to have a model that describes the turbine in this context. Nonlinear model-based control usually require a simple but accurate model. This work aims to produce such a tool. A simple nonlinear model will allow researchers to leverage the large field of nonlinear control theory towards wind-turbine applications. Improved control allows optimal use of the existing generating machinery, as well as better utilization of the local wind resource. In addition, control may extend the service-life of turbines by reducing unwanted motion. This work focuses on the turbine rotor, as good models for substructures and electrical machinery abound in the literature. 2. MODELING THE WIND TURBINE ROTOR The turbine rotor is modeled as a flat disc swept by N , uniformly spaced, rigid blades. In order to yield a tractable analytical model, the following simplifications are made: • The rotor disk is modeled as one rigid structure, that is, we do not consider blade flexibility or root hinges. • The blades have their center of mass and center of lift on the pitch axis. • We assume that the lift curve slope of the blade sections is linear throughout. • Drag is significant only in the axial degrees of freedom, surge and shaft rotation. The modeling of a wind turbine rotor may be summarized eloquently by Figure 1. The aim of our rotor modeling is then to specify the contents of each of the three components in the diagram with an appropriate theory. The dynamic equations of the rotor are given a robot-like vectorial form. This type of model has seen great success in the field of marine control, [Fossen, 2011]. The compact form and physical insight that
Inflow
Wake model
ωz
Circulations
w Lifting Theory
Flow Environment
zbi
xbi
Blade Motion
Rotor Dynamics
ybi φ + δi
Forces and Moments
zr
Ω Controls
u
xr yr
θ
Fig. 1. Components of the rotor model, adopted from [Peters et al., 1989].
v ωy
r
such a model can provide seems very useful for wind-turbine modeling.
pIr R zI
2.1 Turbine Geometry and Kinematics
xI
Consider the rotor depicted in Figure 2. Three reference frames are used: • I, Inertial frame • r, Rotor frame • bi , Blade frames, N in all.
The blade frames and the rotor frame have coincident origins. Each blade, indexed with i, is associated with a unique frame, bi . We let ra ∈ R3 denote a coordinate expressed in frame a, and ei ∈ R3 , denote a unit vector. The cross-product operator is defined as: " # " # a × b := S(a)b,
S(a) =
0 −a3 a2 a3 0 −a1 −a2 a1 0
,
a=
a1 a2 a3
The model is developed in the rotor frame. The rotor frame is nearly stationary under normal operating conditions. Furthermore, this frame relates readily to motions and forces on a support structure. A coordinate in the rotor frame relates to the inertial frame via the transformation: rI = Rr rr Where: Rr (ψ, θ) =
cos(ψ) − sin(ψ) 0 sin(ψ) cos(ψ) 0 0 0 1
·
cos(θ) 0 sin(θ) 0 1 0 − sin(θ) 0 cos(θ)
The blades are spaced by an angle δi = 2π(i − 1)/N for i = 1 . . . N . The individual blade frames are arrived at through the transformation: rri = Rb,i rb with: Rb,i (φ + δi ) =
1 0 0 0 cos(φ + δi ) − sin(φ + δi ) 0 sin(φ + δi ) cos(φ + δi )
The blade pitching axes are parameterized in the longitudinal coordinate r: T rb = [0, 0, r] = rez The matrix: Ri := Rr Rb,i describes the attitude of a turbine blade using Euler angles in a yaw-pitch-roll sequence. The position of a point on the blade longitudinal axis, in the inertial frame, may now written as: rIi = Ri ez r + pIr where the rotor frame is placed at pIr .
yI ψ
Fig. 2. Rotor kinematics. The blade and rotor frames have been separated for clarity. The angle φ describes the rotor azimuth angle. The angles θ and ψ describe the platform pitch and yaw respectively. Quasi-velocities The motion of the rotor frame is expressed in terms of quasi-velocities. In our case, these are the translational velocity of the rotor as well as it’s rotational rate, as seen from the rotor frame. Define: h i ω
r
:=
Ω ωy ωz
,
v
r
T
I
:= Rr p ˙r =
u v w
r
,
ξ :=
v r ω
˙ See Figure 2 for a graphical representation. The where Ω = φ. vector ξ contains the quasi-velocities used to describe platform motion. A Jacobian may be used to relate this vector to a point velocity expressed in any frame, [Spong et al., 2006]: va = Ja ξ (1) Furthermore, Jacobians have the useful property that point forces expressed in any coordinate system may be mapped back to the generalized forces: τ = JTa f a (2) A vector of generalized forces in 6 DOF will be denoted τ ∈ R6 throughout the text. The velocity of the blade axes, as expressed in the blade and inertial frames, will be required. The relevant Jacobians are stated as: viI = Ji ξ,
vib,i = Li ξ
Ji = Rr [I, −S(Rb,i ez r)] ,
(3) Li =
RTb,i
[I, −S(Rb,i ez r)]
2.2 Aerodynamic Model The flow on the rotor is composed of: • External flow; wind and wind shear. • Wake induced flow. • Flow from motion of the rotor. Lift Forces Let w(r) = [wx , wy , wz ]T and f (r) = T [fx , fy , fz ] denote the vector of apparent flow velocity and aerodynamic forces in the blade frame, respectively. These are functions of the radius r. The circulation on the airfoil
is lumped into a point vortex at the quarter-chord, see Figure 3. Assuming small angles, the sectional circulation may be computed as, [Katz and Plotkin, 2001]: Γ = πcwn ' πc(wx + wy (β + α)) (4) where c = c(r), β = β(r) and α are the blade chord, pitch and twist respectively. This approach requires that the blade section operates in the linear regime, i.e. not under stalled conditions at very high angles of attack. The three dimensional vortex lifting xb β+α
wn ' wx + wy β
Γ
wx
wn
yb wy
Fig. 3. Sectional airfoil. law, [Saffman, 1992]: df = ρΓS(w)dl = −ρΓS(ez )wdr (5) yields the differential force on a blade section. The apparent flow w is comprised of the flow due to the wind and the induced velocities, as well as the flow induced by own motion. The form of induced flow used in the Pitt & Peters model, eq.6, may conveniently be expressed in terms of rotor velocities. The mean inflow corresponds to a surge motion, while the cyclic modes correspond to rotor roll and pitch. The external flow may similarly be expressed in terms of rotor displacements. With eq.3, the flow on the blade section may be written: wi = Li −ζ − ξ + KT η where we have introduced the vectors ζ ∈ R6 , which represents the external flow and η ∈ R3 which represents the induced flow. The external flow is produced by the hub-height freestream velocity V∞ entering at an angle ψ∞ w.r.t. reference frame I: RTr (ψ − ψ∞ , θ)ex ζ= V RTr (ψ − ψ∞ , θ)ey (α0 /h) ∞ Wind shear is incorporated by a linearized power law model, taking the hub-height h, and shear exponent α0 as variables. See Fig.4. The constant matrix K ∈ M6×3 maps the induced flow onto the rotor: T
K
=
1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1
We collect these into the apparent flow vector: ν := −ζ − ξ + KT η, ν ∈ R6 This lets us write the differential force on blade i as: dfi = −Γi S(ez )Li νdr Individual Pitch The pitch angle α may vary between blades. A common way to express the individual pitch on the blades, is through cyclic pitch, Johnson [1994]. Here, the pitch on blade i is produced by: αi = [1, cos(φ + δi ), sin(φ + δi )] α, α := [α0 , αc , αs ]T The collective pitch is α0 , while αc and αs are the cyclic pitch. Generalized Aerodynamic Force By eq.2, we may express the distributed aerodynamic forces in terms of the quasicoordinates: dτiA = LTi dfi = −Γi LTi S(ez )Li νdr The generalized aerodynamic force is obtained by summing and integrating over the blades:
τ
A
X
−ρ
=
!
R N
Z 0
Γi L T i S(ez )Li dr
ν
i=1
= ρπA [A0 (ν, α) + cos(N φ)Ac (ν, α) + sin(N φ)As (ν, α)] ν
where the lift forces have been expressed in terms of three lift matrices, which are functions of the apparent flow ν and the individual pitch vector α. The forces are composed of a mean force as well as a cyclic component. The matrices are scaled by the rotor swept area A. The matrix A0 is conveniently the same for any blade count. The matrices Ac and As do not have this property, and must be computed for a specific blade count. The matrices are all skew symmetric: A0 = −AT0 , Ac = −ATc , As = −ATs The matrix elements are furthermore linear in the apparent flow. Note that the pitch control enters through these matrices. This is unconventional, but preserves the skew-symmetry. Model Parameters The lift matrices contain integrals in the rotor chord and twist. These are most conveniently expressed in terms of dimensionless parameters, σn and tn . For different powers of r, these appear as: Z R Z R crn dr =
N 0
π Rn+2 σn , n+1
βcrn dr =
N 0
π Rn+2 tn n+1
In terms of these parameters, the mean force matrix is: S(γ1 ) RG A0 = −RG R2 S(γ2 ) with: γ1 :=
γ2 :=
0 2ν3 (σ0 α0 + t0 ) + Rσ1 (ν6 − ν4 αs ) −2ν2 (σ0 α0 + t0 ) − Rσ1 (ν5 − ν4 αc )
1
4
1 24
0 4ν3 (σ2 α0 + t2 ) + 3Rσ3 (ν6 − ν4 αs ) −4ν2 (σ2 α0 + t2 ) − 3Rσ3 (ν5 − ν4 αc )
G :=
−
R 3
σ1 4
(σ2 α0 + t2 )ν4 −
σ1 2
ν1
1 0 1 0 2 0 0
0
0 1
2
ν2 αc + ν3 αs 0 0 0 (3ν2 αc + ν3 αs )/4 (ν3 αc + ν2 αs )/4 0 (ν3 αc + ν2 αs )/4 (ν2 αc + 3ν3 αs )/4
The matrices As and Ac do not stay similar for various blade counts, and are not presented here. Drag Forces The forces resulting from viscous drag do not have an elegant theory like lift forces. Furthermore, drag forces are usually significantly smaller than lift forces. A complete analytical model of the viscous drag forces on a turbine, in all degrees of freedom, would be prohibitively complex. Such a model would be counterproductive, as simplicity is the main aim of the present work. To avoid an overly complicated model the drag is only analyzed in the surge and roll degrees of freedom, considering purely axial flow. The drag in surge is modeled by a constant drag coefficient, whereas the drag coefficient in the rotational degree of freedom is taken as a curve-fit in the wind speed: 1 1 dx = − ρACD,x ν12 , dφ = − ρAR3 CD,φ (V∞ )Ω2 2 2 Pitt & Peters Wake Model The wake behind the rotor is produced as the rotor exerts a force on the air flowing through it. Wake models range from simple momentum theory, to advanced free vortex or CFD methods, [Bossanyi et al., 2000]. A very attractive wake model suitable for the present work is
found in [Pitt and Peters, 1981]. This model has seen extensive use in helicopter and wind-turbine engineering. Importantly, the approach allows the modeling of yaw misalignment w.r.t to the wind, a common occurrence in operation. The P&P model describes the inflow distribution as: η = η0 + r cos(φ)ηc + r sin(φ)ηs (6) The inflow is scaled by three inflow states, which is gathered in a vector η := [η0 , ηc , ηs ]T . See Figure 4 for a graphical representation. The evolution of these wake states is governed by a set of nonlinear differential equations: A ρARMw η˙ + ρAVw L−1 w (χ)η = −Kτ Here expressed in a dimensional form. We employ K to select the required axial force as well as pitch and roll moments from the generalized aerodynamic force, as required by the wake model. The mass matrix is augmented with an adjustment factor as suggested by Suzuki [2000]. The wake skew angle χ denotes the angle of the wake with respect to the rotor, see Figure 9.
M(φ) = N m
η0
Z
V∞
xI ψ∞
Y
ηs
xr
ψ ηc
yI
r CN (φ)S(ex Ω) − S(ω )CN (φ) 2 r
r
S(ωa )
cg r
cg CN (φ)S(ωa ) rg
ωar
i
S(ω )DN (φ) − DN (φ)S(ex Ω)
r
:= ω − Ωex , and: 0 − cos(φ) − sin(φ) , N =1 0 0 N cos(φ) X 0 sin(φ) 0 CN (φ) := S(Rb,i ez ) = 0 0 0 ∀ N ≥2 i=1 0 00 , 0 0 0 1 0 0 2 0 cos (φ) cos(φ) sin(φ) , N < 3 2 N 0 cos(φ)sin(φ) sin (φ) X 2 1 0 0 DN (φ) := − S (Rb,i ez ) = 1 0 , 0 N ≥3 i=1 0 20 1 where,
2
The integrated quantities are given by: Z R Z R ρ(r)dr, 0
1 cg := m
ρ(r)rdr,
s 1 m
rg :=
0
Z
R
ρ(r)r2 dr 0
Note that a blade count with N < 3 will inevitably result in a time varying system. Utilizing the above results we may write: M(φ)ξ˙ + C(φ, ξ)ξ = τe
X
yr
I −cg CN (φ) 2 cg CN (φ) rg DN (φ)
C(φ, ξ) = N m
m := χ ' ψ − ψ∞
h
Fig. 4. Inflow modes and wake skew angle in the Pitt & Peters model.
This form of the equations of motion avoids the need to consider each blade separately. Furthermore, the equations are in the rotor-frame, thus allowing easy integration with the aerodynamic model. The dynamic model have several notable properties. The mass matrix is symmetric. The coriolis matrix is T skew for N ≥ 2. It may be verified that h symmetric C = −C i T ˙ ξ M(φ) − 2C(φ, ξ) ξ = 0, which implies that the system is conservative, [Siciliano et al., 2009]. 2.4 Integrated Model The model equations may now be posed on the compact form, a set of first order nonlinear ODE’s:
2.3 Dynamic Model A simple dynamic model is derived to complement the aerodynamic model. The acceleration of a particle along the blade axis is: aIi = v˙ iI = J˙ i ξ + Ji ξ˙ The force balance of our particle along the longitudinal axis is: ρ(r)aIi dr = fiI dr where ρ(r) is the mass density along the blade and fi an external forcing. Using eq.2 one can write: ρ(r)JTi aIi dr = JTi fiI dr ⇒ ρ(r)JTi Ji drξ˙ + ρ(r)JTi J˙ i drξ = JTi fiI dr
It is necessary to sum up the contribution for the particles from all blades, N in all. A subsequent integration over r will yield the final dynamic equations. Define:
M(φ) :=
N Z X
N Z X i
T
ρ(r)Ji Ji dr,
C(φ, ξ) :=
0
i
τe :=
R
i
R
Ji
N Z X
T I fi dr
0
Integration and summation yields:
0
R T ˙ ρ(r)Ji J i dr
M(φ)ξ˙ + C(φ, ξ)ξ − ρπA [A0 + cos(N φ)Ac + sin(N φ)As ] ν = τg + τd + τ −1
RMw η˙ + Vw Lw η + πK [A0 + cos(N φ)Ac + sin(N φ)As ] ν = 0 T
ν = −ζ − ξ + K η
where τd denotes the forces due to drag. The rotor system will be connected to a substructure exerting a force on the rotor: τ . Gravity forces enter through τg , but are trivial, as a balanced flat rotor will act as point mass at the point of connection. This model bears great resemblance to vectorial models of marine craft, Fossen [2011]. A particular note of interest is the skew symmetry of the lift matrices, which has not been noted previously. Skew symmetry may play an important role in stability proofs for control systems, Siciliano et al. [2009]. 3. PARAMETER ESTIMATION The parameters in the lift matrices serve as good first estimates for the model. However, tuning these parameters against a more advanced approach, or empirical data, will yield a better model. To this end, each parameter is augmented with an adjustment factor, e.g. σ ˆ1 = σ1 . The drag model is also fitted. The Pitt & Peters model is based on an actuator disc model of the rotor, Pitt and Peters [1981]. An actuator disc does
not model the effects of a finite number of blades, rather, it considers the aerodynamic force to be smeared over an infinitely thin disc, Bossanyi et al. [2000]. The effects of having a finite number of blades is that the induced flow on the blades, will be larger than the disc average. This is often referred to as tip-loss. The disc averaged uniform induction, η0 , is used in the Pitt & Peters model. An adjustment factor is introduced to capture the effects of the tip loss, as well as non-uniform inflow, η˜0 = kη η0 , where η0 denotes the induction as it appears on the blades. The aim of the parameter estimation is to optimize the adjustment factors, fit the drag model, as well as setting kη . To achieve this, a numerical BEM tool is used to generate a large set of data points, covering the operational envelope of the turbine. BEM is the industry standard means of obtaining wind turbine performance estimates, and can be relied on for accurate data. Such a program yields the steady state axial forces and torques on the turbine, as it operates in axial flow. The steady axial subsystem of the present aerodynamic model is extracted; rows and columns 1 and 4 in the mean lift matrix A0 . The present aerodynamic model is unsteady, so the equilibrium points of the model are fitted against the data-set. The equilibrium of the model in axial flow is when: 1 4(1 − kη a)kη a = πλ [3(1 − a)σ1 − 2λ (t2 + α0 σ2 )] (7) 3 where dimensionless variables have been used: a = η0 /V∞ , λ = ΩR/V∞ Let a? denote the physical solution to eq. 7, (the other solution is discarded): q 12kη + 3πλσ1 −
?
a
2
9 4kη + πλσ1
=
Fig. 5. A typical behavior for pitch-controlled, variable-speed turbines. The pitch is activated as the maximum power is reached, to limit power, in this case 5MW. The tip speed ratio is decreased as the rotational rate is kept constant for high wind speeds. procedure will in fact work poorly if the data-set is confined to the reference trajectory. The model has two significant limitations; the P&P model will be unstable for large induction ratios, a > 0.5. Also, departures from a linear lift slope will violate the assumptions inherent in the lift modeling. To avoid fitting the model to operating conditions that lies beyond it’s scope, data points exceeding these limitations are removed. If these provisions are not taken, gross errors resulting from limitations in the model will drive the optimization procedure towards a less desirable result. The data points are subsequently weighed. The weighing is
2 (−3σ + 2λ (t + α σ )) + 48πλkη 1 2 0 2
2 24kη
Four dimensionless quantities, are used to fit the parameters in the model. CT (α0 , λ) = CP (α0 , λ) =
T 1 ρAV 2 ∞ 2
P 1 ρAV 3 ∞ 2
'
1 3
?
πλ 3(1 − a )σ1 − 2λ (t2 + α0 σ2 ) ?
' (1 − a )CT (α0 , λ)
CB (α0 , λ) =
N mB 1 ρARV 2 ∞ 2
CS (α0 , λ) =
N fS 1 ρAV 2 ∞ 2
'
1 6
?
πλ 4(1 − a )σ2 − 3λ (t3 + α0 σ3 ) ?
?
' π(1 − a ) 2(1 − a )σ0 − λ (t1 + α0 σ1 )
Here a? has been inserted. These are the thrust, power, out of plane root bending and in-plane root shear coefficients respectively. The parameter t0 does not appear in these coefficients, and is used without an adjustment. By inserting a? into the coefficients, one fits the equilibria in the dynamic model towards the data set. One has furthermore removed the need to fit the induction itself. The induction is simply an internal state of the model, enforcing correct steady state behavior. 3.1 Selection and Weighing of Data Set To produce a data set for fitting, a BEM code is used to generate a typical reference trajectory in terms of the pitch angle and tip-speed ratio, with respect to varying wind speed. See Fig.5 for an illustration. The data points to be fitted are generated by adding normally distributed random offsets in pitch-angle and tip-speed ratio from the reference trajectory. The fitting 1
Blade Element Momentum
Fig. 6. Left: Comparison of the wind speed probability function at the site (red), and in the data set (histogram). Right: Data set, color coded with respect to weighing. done so as to render the weighted likelihood of encountering a certain wind speed in the data set, similar to the probability at the turbine site. This probability may be assumed to be on the form of a Weibull curve, Bossanyi [2003]: V k−1 f (V∞ ) = k ∞k exp −(V∞ /c)k c Figure 6 illustrates the result of the weighing. The parameters in the model, are subsequently fitted by weighted least squares. This is achieved by the Gauss-Newton method. The optimum is found within 4-10 iterations. The lift and drag forces are fitted separately, and one must use a BEM code that gives separate coefficients for the potential and viscous effects. Fitting towards the sum of these yields poor results. The four coefficients must be fitted jointly, due the linking introduced by the equilibrium induction.
4. MODEL VALIDATION
Mz @kNmD 50
The fitted model replicates the equilibrium forces from the BEM code with good results. The model has been run with
100
150
200
250
300
350
f@degD
-20 -40 -60 -80 -100 -120 -140
Fig. 9. Yaw moment due to a yaw misalignment on the 2MW Tjæreborg turbine. The yaw offset is 32◦ , wind speed is 8.5[m/s]. Model is dashed.
Fig. 7. BEM results compared to fitted model for a 5MW turbine. Model data are dashed. the same controller that were used for the BEM. Note the very slight deviation in the pitch angle at high wind speeds, which is likely caused by violation of the small angle assumption in eq. 4. The point of optimum power extraction CP? = CP (α0? , λ? ) is also predicted very well. An accurate estimate of the power Model BEM
? CP
α?0
λ?
0.481 0.482
-0.219◦
7.520 7.418
-0.746◦
optimum is important if one desires to develop optimizing control methods. The dynamic wake implementation is also
Fig. 8. The nonlinear model compared to experiment and a lookup-table approach at wind speeds 12.5[m/s](top) and 7.4[m/s](bottom). The response to a pitch step is depicted. The experimental data are from the 2MW Tjæreborg turbine. quite successful as seen in Figure 8 2 . Note that the proposed model captures the transients in the experiment, whereas the lookup approach fails in this respect. A final test is the response to a large yaw misalignment. Although the model yields a reasonable mean value of the resulting yaw moment Mz , the cyclic component is very much underestimated. The reason for this failure is unknown at present. A likely cause is the very 2
The Authors thanks Martin O.L. Hansen for providing experimental data.
simple inflow model being unable to capture the complicated unsteady flow field in a yawed flow. Under-prediction of cyclic forces is indeed present in many inflow models, Bossanyi [2003]. It may be better to retain only the steady component in simple models such as the present one. REFERENCES F. Bianchi, H. De Battista, and R. Mantz. Wind Turbine Control Systems: Principles, Modelling and Gain-scheduling Design (Advances in Industrial Control). 2006. E. Bossanyi, T. Burton, D. Sharpe, and N. Jenkins. Wind energy handbook. John W iley, 2000. EA Bossanyi. Wind turbine control for load reduction. Wind Energy, 6(3):229–244, 2003. ISSN 1099-1824. CL Bottasso, A. Croce, and B. Savini. Performance comparison of control schemes for variable-speed wind turbines. In Journal of Physics: Conference Series, volume 75, page 012079. IOP Publishing, 2007. T.I. Fossen. Handbook of Marine Craft Hydrodynamics and Motion Control. Wiley, 2011. W. Johnson. Helicopter theory. Dover Pubns, 1994. J. Katz and A. Plotkin. Low-speed aerodynamics, volume 13. Cambridge University Press, 2001. L. Mihet-Popa, F. Blaabjerg, and I. Boldea. Wind turbine generator modeling and simulation where rotational speed is the controlled variable. Industry Applications, IEEE Transactions on, 40(1):3–10, 2004. H. Namik and K. Stol. Individual blade pitch control of floating offshore wind turbines. Wind Energy, 13(1):74–85, 2010. ISSN 1099-1824. D.A. Peters, D.D. Boyd, C.J. He, et al. Finite-State InducedFlow Model for Rotors in Hover and Forward Flight. Journal of the American Helicopter Society, 34:5, 1989. D.M. Pitt and D.A. Peters. Theoretical prediction of dynamicinflow derivatives. Vertica, 5(1):21–34, 1981. PG Saffman. Vortex force and bound vorticity. Vortex Dynamics, 1992. B. Siciliano, L. Sciavicco, and L. Villani. Robotics: modelling, planning and control. Springer Verlag, 2009. M.W. Spong, S. Hutchinson, and M. Vidyasagar. Robot modeling and control. John Wiley & Sons Hoboken, NJ, 2006. K.A. Stol and M.J. Balas. Periodic disturbance accommodating control for blade load mitigation in wind turbines. Journal of solar energy engineering, 125:379, 2003. A. Suzuki. Application of dynamic inflow theory to wind turbine rotors. 2000.