- Email: [email protected]
A phenomenological model for the dynamic response of wind turbines to turbulent wind
ARTICLE IN PRESS
Journal of Wind Engineering and Industrial Aerodynamics 92 (2004) 159–183
A phenomenological model for the dynamic response of wind turbines to turbulent wind Alexander Rauh*, Joachim Peinke Institut fur . Physik, Universitat . Oldenburg, D-26111 Oldenburg, Germany Received 11 December 2002; received in revised form 10 November 2003; accepted 12 November 2003
Abstract To predict the average power output of a wind turbine, a response model is proposed which takes into account: (i) the delayed response to the longitudinal wind speed fluctuations; (ii) a response function of the turbine with arbitrary frequency dependence; (iii) wind fields of arbitrary turbulence intensity. In the limit of low turbulence intensity, the dynamical ansatz as proposed in 1992 by Rosen and Sheinman is reproduced. It is shown, how the response function of the turbine can be obtained from simulation experiments of a specific wind turbine. For two idealized situations the dynamic effect of fluctuating wind is estimated at turbulence intensities 0pIu p0:5: At the special mean wind speed V ¼ 8 m=s; the turbine response function is determined from simulation data published by Sheinman and Rosen in 1992 and 1994. r 2003 Elsevier Ltd. All rights reserved. Keywords: Wind turbine; Power prediction; Dynamic response
1. Introduction The prediction of the average power output of a wind turbine, which is accurate enough for an economic assessment, is still a challenging problem. We will focus in this article on the question of how longitudinal wind speed fluctuations affect the power output through a delayed, frequency-dependent response by the wind turbine. The phenomenological model to be proposed includes an arbitrary, time-dependent response function, which, normally, has to be determined from specific aerodynamic *Corresponding author. Tel.: +49-441-798-3460; fax: +49-441-798-19-3460. E-mail address: [email protected] (A. Rauh). 0167-6105/$ - see front matter r 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.jweia.2003.11.001
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
160
and electric–mechanical modelling of the wind turbine. Our model allows for considering arbitrary turbulence intensities. We have not included a possible feedback interaction to the turbine by the power grid, or interferences of turbines within a wind farm, which were considered in a recent paper by Hansen et al. [1]. In another recent study of dynamical effects by Pedersen et al. [2], the electric– mechanical response of a wind farm within an islanding experiment was emphasized rather than the effect of wind speed fluctuations. Our model generalizes a dynamic ansatz proposed by Rosen and Sheinman [3,4] for the case of weak turbulence intensity. Let us briefly sketch their ansatz. They start from the stationary power curve Ls ðuÞ; where u is the axis parallel component of the wind velocity. The instantaneous speed uðtÞ is decomposed as Z 1 T uðtÞ ¼ V þ vðtÞ with V uðtÞ ¼ dtuðtÞ; ð1Þ T 0 where T is a short-term averaging time, e.g. 10 min: The turbulent part, vðtÞ; is a stochastic variable described, according to Shinozuka and Jan [5], as follows: vðtÞ ¼
N pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X 2Sðfi ÞDf sin ð2pfi t þ yi Þ;
ð2Þ
i¼1
where the frequencies fi ¼ iDf ; i ¼ 1; 2; y; are equidistant, yi is a random phase equally distributed over the interval from 0 to 2p; and for iak the phases yi and yk are statistically independent. We will denote in this article the phase average of a function H as Z 2p 1 dy1 y dyN H: ð3Þ /HS ¼ ð2pÞN 0 Clearly, the mean value /vS ¼ 0; and /v2 S ¼
N X
Sðfi ÞDf :
ð4Þ
i¼1
The function Sðf Þ is the spectral density of the turbulent part v of the longitudinal velocity component, and depends, in general, on the short-term mean speed V and on the hub height of the turbine. The corresponding turbulence intensity, Iu ; is defined through ffi 1 pffiffiffiffiffiffiffiffiffiffiffi /v2 S; Iu ¼ ð5Þ V where the index u refers to the longitudinal component u of the wind velocity. In two applications in Section 8, we will adopt for Sðf Þ a distribution function from boundary layer meteorology. In a next step, Rosen and Sheinman [3,4] consider the case of weak turbulence by the following Taylor series up to second order in v: Ls ðuÞ ¼ Ls ðV Þ þ l1 ðV Þv þ l2 ðV Þv2 :
ð6Þ
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
161
In order to include dynamic effects, they modify the stationary average /Ls ðuÞS ¼ Ls ðV Þ þ l2 /v2 S as follows: Ldyn ðV Þ :¼ Ls ðV Þ þ l2
N X
Df Sðfi ÞG1 ðfi ; V Þ:
ð7Þ
i¼1
The function G1 was derived in [3] from simulations by means of a specific rotorgearbox-generator model. The final power prediction was then established by the ensemble average with respect to the mean velocities V : The above dynamic modification corresponds to a coupling between the turbine, through G1 ; and the wind field, through S: As a matter of fact, our generalization was partly guided by the condition to obtain the form (7) in the limit of weak turbulence. The model, to be defined in Section 2, is a linear differential equation for the instantaneous power output LðtÞ of the turbine, and is characterized by a time-dependent relaxation function rðtÞ: The solution of this equation depends on turbine response and wind field in a nonlinear way. Nevertheless, it turns out that phase averaging with respect to the stochastic wind field (2) can rigorously be reduced to a one-dimensional integral containing a Gaussian distribution function. From the general solution, we will derive the weak turbulence limit and thus achieve a connection of the turbine function gðtÞ of our model to the spectral function G1 ðf Þ defined in (7) according to Sheinman and Rosen [3]. As a nontrivial result, we show that gðtÞ can be uniquely determined from the Rosen–Sheinman function G1 : To this end, we invoke causality for gðtÞ; i.e. gðtÞ ¼ 0 for to0 and exploit a Kramers–Kronig relation between the real and imaginary parts of gðtÞ in frequency space. The necessary proofs and, in particular, the solution of a singular integral equation will be outlined in an appendix together with a simple example for demonstration. In two rather idealized applications, we model the timely shutdown behaviour of a wind turbine and calculate the average power output for the special stationary power curve (6) in dependence of the turbulence intensity 0pIu p0:5: To this end, the turbine response function is estimated from graphical data published in [3,4]. Thereby, we restrict ourselves to the data for the mean wind speed V ¼ 8 m=s: In general, the turbine function gðtÞ should be determined by simulation runs of a specific turbine model. A simple way, within our model, consists in simulating the response to velocity jumps. We will discuss this briefly by an example in Section 3. In frequency space, gðtÞ can be obtained, in principle, by simulating the response to periodic wind gusts of given frequency and deterministic phase. As is shown in Appendix D, one has to choose both time and phase in a proper way, in order to avoid involved relations. The partly involved algebraic manipulations of this study and the numerical evaluations were performed with the aid of Mathematica [6].
ARTICLE IN PRESS 162
A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
2. Response model We denote by LðtÞ the instantaneous power output of the wind turbine at time t: It is assumed to be uniquely related to the instantaneous longitudinal wind speed: t-uðtÞ measured, e.g. at the nacelle, or at a certain distance upwind from the turbine as recommended in [7]. The function Ls ðuðtÞÞ corresponds to the stationary power to which the turbine would relax, if the particular speed uðtÞ was frozen in. The power curve u-Ls ðuÞ is taken to be known. The following relaxation model is adopted: dLðtÞ ¼ rðtÞ½Ls ðuðtÞÞ LðtÞ: ð8Þ dt Apart from mathematical simplicity, the form of this model was motivated by two requirements: (1) it should reproduce the dynamical power response (7) in the limit of weak turbulence, as postulated in [3]; (2) it should be more general in order to include arbitrary turbulence intensities. The relaxation function rðtÞ; in general should depend both on the wind speed uðtÞ ¼ V þ vðtÞ and the time derivative u’ ¼ v’: We assume the form rðtÞ ¼ r1 ðuðtÞÞ þ r2 ð’vÞ:
ð9Þ
The part r1 describes then, e.g. relaxation, at constant wind speed, towards the stationary power Ls ðuÞ from a nonstationary value Lð0ÞaLs ðuÞ: For simplicity we will assume r1 as a constant, i.e. depending on the mean wind speed V only with r1 ðuÞ ¼ r1 ðV Þ ¼ r0 : The function r2 ð’vÞ is chosen within the scheme of linear response theory as follows: Z t r2 ðtÞ ¼ dt0 v’ðt0 Þgðt t0 Þ; ð10Þ N
where g describes the response of the turbine. The response function g; in general, depends, on V and Iu ; as control parameters. Within a particular time interval tA½0; T these parameters are available, in principle, from the preceding short-time period. The function r2 ¼ r2 ð’vÞ is chosen to model the instantaneous interaction of the turbine with the turbulent wind. The turbine response has to be causal, i.e. gðt t0 Þ ¼ 0 if t0 > t: As a consequence, in frequency space we have the product of the Fourier transforms of v’ðt0 Þ and gðt t0 Þ: N X pffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi 2pfi ri Df cos ð2pfi t þ yi þ gðfi ÞÞ; ri ¼ Gðfi Þ 2Sðfi Þ; ð11Þ r2 ðtÞ ¼ i¼0
where the Fourier transform of g is written in terms of modulus Gðf Þ ¼ Gðf Þ and phase gðf Þ ¼ gðf Þ: Z pffiffiffiffiffiffiffi 1 N dt gðtÞ exp ð2pIf tÞ ¼ Gðf Þ exp ðIgðf ÞÞ; I ¼ 1: ð12Þ 2p 0 As a first application, we will study in a later section the effect of large turbulence intensity in a simplified shutdown model. To this end we assume that the stationary
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
power curve depends on the turbulence intensity as follows: ( Ls ðu; 0Þ if Iu pI ; Ls ðuÞ-Ls ðu; Iu Þ ¼ 0 if Iu > I :
163
ð13Þ
3. Some simple examples To get some insight into the model, we consider three simple cases: (i) constant speed u0 with v’ 0 and an initial power Lð0ÞaLs ðu0 Þ; (ii) a constant speed jump from u1 to u2 at time t ¼ 0; i.e. v’ ¼ ðu2 u1 ÞdðtÞ; (iii) an idealized shutdown at time t ¼ 0 and instantaneous power output Lð0Þ in the turbulent wind field (2). In example (i) both r ¼ r0 and Ls ðu0 Þ are constant and the solution of (8) exponentially approaches Ls ðu0 Þ as LðtÞ ¼ Ls ðu0 Þ þ ðLð0Þ Ls ðu0 ÞÞ exp ½r0 t:
ð14Þ
The relaxation time 1=r0 is, e.g., of the order of a few seconds in the case of response measurements of blade root bending moments to sudden changes of the blade pitch angle on the Tjaereborg Turbine [8]. In case (ii) we have rðtÞ ¼ r0 þ ðu2 u1 ÞgðtÞ; Lð0Þ ¼ Ls ðu1 Þ; Ls ðtÞ ¼ Ls ðu2 Þ ¼ const: for t > 0: With the definition Z t # ¼ RðtÞ dt0 rðt0 Þ ¼ r0 t þ RðtÞ;
ð15Þ
RðtÞ ¼ ðu2 u1 Þ
0
Z
t
dt0 gðt0 Þ
ð16Þ
0
the solution of (8) can be written # # LðtÞ ¼ exp½RðtÞLð0Þ þ Ls ðu2 Þ exp½RðtÞ
Z Z
t
# 0 Þ dt0 rðt0 Þ exp ½Rðt
0 t
d # 0 Þ exp ½Rðt dt0 0 # # # 0 Þ 1Þ ¼ exp ½RðtÞLð0Þ þ Ls ðu2 Þ exp ½RðtÞðexp ½Rðt # ¼ Ls ðu2 Þ þ ðLð0Þ Ls ðu2 ÞÞ exp ½RðtÞ: # # ¼ exp½RðtÞLð0Þ þ Ls ðu2 Þ exp½RðtÞ
dt0
ð17Þ
If RðtÞ as defined in (16) is bounded in time, which we assume, then the power LðtÞ exponentially converges towards the stationary power Ls ðu2 Þ once more. In case (iii) we implement the stochastic wind field (2) and redefine the function RðtÞ as follows: Z t N pffiffiffiffiffiffi X Df ri sin ð2pfi t þ yi þ gðfi ÞÞ: dt0 rðt0 Þ ¼ r0 t þ RðtÞ Rð0Þ; RðtÞ ¼ 0
i¼0
ð18Þ
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
164
In our shutdown model the stationary power Ls 0 for t > 0: The solution of (8) therefore reads as follows: LðtÞ ¼ exp ½r0 t RðtÞ þ Rð0Þ Lð0Þ:
ð19Þ
The average with respect to the phases yi can rigorously be carried out. As is shown in Appendix A, the result of phase averaging leads to /LðtÞS ¼ exp ½r0 t þ w2 ðtÞ Lð0Þ;
ð20Þ
where w2 ðtÞ ¼
N X
Df r2i sin2 ðpfi tÞ:
ð21Þ
i¼1
As is seen, the stochastic wind field gives rise to an effectively larger relaxation time, i.e. shutdown is slowed down to some extent. Nevertheless, the shutdown state of the turbine eventually is reached by an exponential law, since w2 ðtÞ is bounded in time.
4. General integration and phase averaging In terms of RðtÞ as defined in (18), the solution of the linear differential equation (8) can be written as follows: Z t LðtÞ ¼ exp ½r0 t RðtÞ þ Rð0Þ Lð0Þ þ dt0 exp ½r0 ðt t0 ÞPðt; t0 Þ; ð22Þ 0
Pðt; t0 Þ ¼ rðt0 Þ exp ½RðtÞ þ Rðt0 Þ Ls ðV þ vðt0 ÞÞ:
ð23Þ
The phase averaging of the initial value term was carried out in example (iii), see (20) and (21), with the result that this term exponentially decays to zero. Therefore, we will omit this term henceforth, and consider Z t LðtÞ ¼ dt0 exp ½r0 ðt t0 ÞPðt; t0 Þ: ð24Þ 0
As is shown in Appendix B, we obtain the following phase average (we set t ¼ t t0 )): sffiffiffiffiffiffiffiffiffiffiffi 1 /PðtÞS ¼ exp ½C2 ðtÞ 4pC3 Z N C13 1 2
dz r0 2C12 ðtÞ þ z exp z 4C3 C3 N
Ls ðV 2C23 ðtÞ þ zÞ;
ð25Þ
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
165
where C1 ¼
N 1X Df Z2i ; 4 i¼1
C2 ðtÞ ¼
N X
N 1X Df Zi ri sin ð2pfi tÞ; C12 ðtÞ ¼ 4 i¼1
C23 ðtÞ ¼
Df r2i sin2 ðpfi tÞ;
C3 ¼
i¼1
C13 ¼
N 1X Df z2i ; 4 i¼1
N 1X Df Zi zi sin ðgðfi ÞÞ; 4 i¼1
N 1X Df zi ri ðcosðgðfi ÞÞ cos ð2pfi t þ gðfi ÞÞÞ 4 i¼1
ð26Þ
with Zi ¼ 2pfi ri ;
zi ¼
pffiffiffiffiffiffiffiffiffiffiffiffi 2Sðfi Þ;
ð27Þ
the magnitude C1 appears in an intermediate step of Appendix B. Since the phase averaged function Pðt; t0 Þ depends on the time difference only, from (24) we obtain Z t /LðtÞS ¼ dt exp ½r0 t /PðtÞS: ð28Þ 0
It should be noticed that the half-width of the Gaussian exponent in the integrand of (25), 1=ð4C3 Þ ¼ V 2 =ð2Iu2 Þ; depends on mean values of the wind speed distribution only. Furthermore, C2 ðtÞ ¼ w2 ðtÞ which occurred in the shutdown example, see (21). Moreover, it is seen that the ensemble average of the power output depends in a nonlinear way on the spectrum both of the wind gusts and the turbine response function in general.
5. The limit of weak turbulence In this section we validate the model, to some extent, by deriving from (25) and (28) the ansatz (7) of Rosen and Sheinman. Thereby we will relate the adhoc function G1 ðf ; V Þ; proposed in [3], to the Fourier transform of the turbine function gðtÞ: The strength of the turbulent wind field will be marked by multiplying v with the nondimensional parameter a; which eventually will be set back to the value 1. pffiffiffiffiffiffiffiffiffi Equivalently, Sðf Þ is multiplied by a in expression (25). In the limit a-0; the C coefficients defined in (26), are proportional to a2 : In view of the exponential term of the z-integrand of (25), z is effectively of order a and is scaled by this factor. Referring to (6), we can write for part of the integrand in (25) hðzÞ :¼ exp ½a2 C2 ðr0 2a2 C12 þ azC13 =C3 ÞLs ðV 2a2 C23 þ azÞ ¼ r0 Ls ðV Þ þ azðC13 Ls ðV Þ þ C3 r0 l1 Þ=C3 þ a2 fLs ðV Þð2C12 þ C2 r0 Þ 2C23 r0 l1 þ z2 ðr0 l2 þ C13 =C3 l1 Þg þ Oða3 Þ:
ð29Þ
ARTICLE IN PRESS 166
A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
The z-integration leads to the result, up to terms of order a3 : sffiffiffiffiffiffiffiffiffiffiffi Z N 1 1 2 dz hðzÞ exp z ¼ r0 Ls ðV Þ /PðtÞS ¼ 4pC3 N 4C3 þ a2 f2C12 ðtÞLs ðV ÞÞ C2 ðtÞr0 Ls ðV Þ þ 2C13 l1 2C23 ðtÞr0 l1 þ 2C3 r0 l2 g:
ð30Þ
The time integration, according to (28), amounts to integrals of the type Z t m cos ðnÞ þ r0 sin ðnÞ dt exp ½r0 t sin ðmt þ nÞ ¼ m2 þ r20 0 m cos ðmt þ nÞ þ r0 sin ðmt þ nÞ exp½r0 t : m2 þ r20 ð31Þ In view of the time averaging of LðtÞ we will neglect the exponentially decaying terms, which will cause errors of the order 1=ðr0 TÞ; where T is the averaging time. Restoring the ri ; Zi ; zi from (11) and (27) and restoring a ¼ 1; we arrive at the following weak turbulence limit (up to terms proportional to 1=ðr0 TÞ): /LðtÞS ¼ Ls ðV Þ þ l2
N X
Df Sðfi Þ þ DLdyn
ð32Þ
i¼1
with DLdyn ¼ l1
N X
Df
i¼1
4p2 fi2 Sðfi ÞGðfi Þ 2pfi cos ðgðf ÞÞ sin ðgðf ÞÞ : i i r0 4p2 fi2 þ r20
ð33Þ
The symbol DLdyn denotes the average excess power as compared to the quasistationary average /Ls ðV þ vðtÞÞS ¼ Ls ðV Þ þ l2 /vðtÞ2 S þ Oðv3 Þ; it can be positive or negative in principle. When expressions (32) and (33) are compared with (7), we obtain the function G1 of Rosen and Sheinman [3] (up to normalization by Ls ðV Þ) in the form G* 1 ðf ; V Þ l2 G1 ðf ; V Þ ¼ l2 þ DLdyn : ð34Þ A slightly more compact form can be obtained by introducing the angle Yðf Þ defined by tan ðYÞ ¼ 2pf =r0 : ð2pf Þ2 G* 1 ðf ; V Þ ¼ l2 þ l1 Gðf Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cos ðgðf Þ þ Yðf ÞÞ: r0 4p2 f 2 þ r20
ð35Þ
The coefficients l1;2 depend on the mean wind speed V : Since l2 ðV Þ may have zeros, we prefer to replace G1 by G* 1 ; and write the phase and time averaged power output in the limit of weak turbulence intensity as follows: /LðtÞS ¼ Ls ðV Þ þ
N X i¼1
Df Sðfi ÞG* l ðfi ; V Þ:
ð36Þ
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
167
As is seen from (35), one has optimal phase adaption of the turbine for gðf Þ ¼ Yðf Þ which gives rise to a positive excess power DLdyn in the weak turbulence limit.
6. Kramers–Kronig analysis For later applications we will derive a specific turbine function gðtÞ from the spectrum G1 ðf Þ; which is graphically represented in Fig. 2 of Ref. [3]. For this purpose, we will give the connection between a general response spectrum G1 ðf Þ and the underlying turbine function gðtÞ in the following. For convenience, we will replace f by the angular frequency o ¼ 2pf : In frequency space we denote the real and imaginary parts of the turbine function by GR and GI ; respectively: Z N 1 GR ðoÞ þ IGI ðoÞ GðoÞ ðcosðoÞ þ I sinðoÞÞ ¼ dtgðtÞ exp ½Iot: ð37Þ 2p N Furthermore, in view of (33) and (34), we introduce the function F1 ; which is related to G1 ; as follows: G1 ðoÞ ¼ 1 þ
l1 o2 F1 ðoÞ: l2 o2 þ r20
ð38Þ
From (33) and (34) we obtain the following connection of F1 to the turbine function: o ð39Þ GR ðoÞ GI ðoÞ ¼ F1 ðoÞ: r0 Because the turbine function gðtÞ is causal with gðtÞ ¼ 0 for to0; the imaginary part GI is related to GR by the Kramers–Kronig relation [9] Z N 1 GR ðxÞ ; ð40Þ GI ðoÞ ¼ þ P dx p xo N where P denotes the principal value; there is a ‘‘þ’’—sign in (40) since, as compared to [9], our definition (37) of the Fourier transform differs through o- o: We thus arrive at the following integral equation for the real part GR of the turbine function: Z N o1 GR ðxÞ P ¼ F1 ðoÞ: GR ðoÞ dx ð41Þ r0 p xo N The above integral equation is solved in Appendix C.1. using time space. The corresponding magnitudes in time space are denoted by lower case letters: GR ðoÞ-gR ðtÞ; F1 ðoÞ-f1 ðtÞ: Because the turbine function gðtÞ is real, we have the symmetries: GR ðoÞ ¼ GR ðoÞ;
GI ðoÞ ¼ GI ðoÞ;
gR ðtÞ ¼ gR ðtÞ;
f1 ðtÞ ¼ f1 ðtÞ; ð42Þ
ARTICLE IN PRESS 168
A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
where the last relation is a consequence of (41). Since gðtÞ is not analytic at t ¼ 0; the solution of (41) has to be written separately for t > 0 and to0: Z t gR ðtÞ ¼ exp ½r0 tgR ð0Þ þ r0 dt0 exp ½r0 ðt t0 Þf1 ðt0 Þ for t > 0; 0 Z t dt0 exp ½þr0 ðt t0 Þf1 ðt0 Þ for to0: ð43Þ gR ðtÞ ¼ exp ½þr0 tgR ð0Þ r0 0
It is easily seen, that the above solution obeys the symmetry stated in (42). The turbine function g; we are looking for, is now simply obtained from the following relation, which is generally proved in Appendix C.2.: ( 2gR ðtÞ if t > 0; gðtÞ ¼ ð44Þ 0 if to0: The still open constant gR ð0Þ in solution (43) remains to be discussed. Clearly, if the turbine function is continuous at t ¼ 0; then gðt ¼ 0Þ ¼ 0; and because of (44) the constant gR ð0Þ ¼ 0: On the other hand, a jump of gðtÞ at t ¼ 0 is reflected in the behaviour of the Fourier transform for o-N: As is shown in Appendix C.3, the constant is uniquely determined by F1 as follows: gR ð0Þ ¼ pr0 lim F1 ðoÞ: o-N
ð45Þ
7. Deriving a turbine function In the following we attempt to infer a concrete turbine function gðtÞ from graphical data of Fig. 2 in [3]; we restrict ourself to the case of mean wind speed V ¼ 8 m=s: To derive the function F1 ðoÞ from G1 by means of (38), we have to estimate the parameters l1;2 ðV Þ; for which, once more, we have only graphical data available as presented in [3,4]. The two l parameters, essentially, are the first and second derivatives of the stationary power curve Ls ðV Þ: As should be remarked, a measured power curve is based on 10-min time averages, in general, and differs more or less from the ideal power curve inferred from spatiotemporally constant winds. In our group we have a new method of stochastic analysis to extract the deterministic law of a physical magnitude in a noisy environment [10]. A first application to extract the ideal power curve from given data has been published in [11]. To obtain Ls ðV Þ underlying the Refs. [3,4], we multiply the power coefficient Cp0 ; shown in Fig. 2, of [4], by K V 3 ; to obtain the stationary power Ls ðV Þ ¼ K Cp0 V 3 ; where K is a constant. Within the plot range of Cp0 of [4], namely 6:75pV p12 ms1 ; the power curve can be fitted by a straight line: Ls E constant þð87 m2 =s2 Þ KV ; In this way we infer l1 Eð87 m2 =s2 Þ K: The parameter l2 ; which according to (6) is proportional to the second derivative of Ls ðV Þ; is inferred from the nondimensional function dðV Þ introduced in [3] and plotted there in Fig. 1. The
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
169
5
4
G1
3 2 1
0.01
0.03
0.1
0.3
1.0
f Fig. 1. Approximated Rosen–Sheinman [3] response function G1 ðf Þ for mean wind speed V ¼ 8 m=s: The discrete points have been read off from the original. Apart from the straight line, the curve is an interpolating polynomial of order six in the variable x ¼ xmax ð1 þ 12 Log10 ðf ÞÞ:
connection to l2 is 3 l2 ¼ dðV Þ 2 Ls ðV Þ ¼ 3dðV ÞKCp0 V : V Thus, from the data of [3,4], we obtain approximately the following ratio: l1 87 m2 =s2 E11:5 m=s; E l2 3Cp0 V
ð46Þ
ð47Þ
where at V ¼ 8 m=s we read off dðV ÞE34 and Cp0 ¼ 0:42 from Fig. 1 in [3] and Fig. 2 in [4], respectively. In the next step we read off from Fig. 2 in [3] a set of 12 discrete curve points taking care of the logarithmic scale of the abscissa: ðxðfi Þ; G1 ðfi ÞÞ; i ¼ 1; 2; y; 12 with xðfi Þ ¼ xmax ð1 þ 1=2 Log10 ðfi ÞÞ: The logarithmic plot was reasonably approximated by a sixth order polynomial in the variable xðf Þ: The corresponding interpolation curve, is plotted in Fig. 1 together with the dicrete points read off from the original plot. We now interpolate the function F1 as defined in (38) by a polynomial of order 10 in the variable f ; which replaces the logarithmic variable xðf Þ: As a function of P i o ¼ 2pf the interpolating polynomial 10 a o is characterized by the following i¼0 i coefficients: i ai
0 6:33639
1 2 45:3742 126:113
3 192:997 9
10
0:222442
0:009655:
i
6
7
ai
47:3118
12:9315 2:23907
8
4 183:530
5 114:204
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
170
As is seen from Fig. 1, we set approximately G1 ðf Þ ¼ 1 in the interval 0pf p0:05075: There is, of course, no data basis to extrapolate G1 ðf Þ to large frequencies f : We will set G1 ¼ 1 for large f ; which implies the asymptotic value F1 ¼ 0: According to Fig. 1, the frequency support of G1 and thus of F1 is in the domain f p1 so that time space resolution is not below, say 1 s: It is therefore compatible with the data to assume the turbine function gðt ¼ 0Þ ¼ 0; and as a consequence of (45) F1 ðoÞ ¼ 0 for o-N: We will therefore approximate F1 ðoÞ further by setting F1 ðoÞ ¼ 0 for oXo2 ; where F1 ðo2 Þ ¼ 0: 8 for 0popo1 and o > o2 ; > <0 P10 i ð48Þ F1 ðoÞ ¼ for o1 popo2 ; o1 ¼ 0:3188; o2 ¼ 3:970; i¼0 ai o > : F1 ðoÞ for oo0: With the aid of (43) and (44), and setting gR ð0Þ ¼ 0; we now calculate the turbine functions gðtÞ as follows: Z t gðtÞ ¼ 2r0 dt0 exp ½r0 ðt t0 Þ f1 ðt0 Þ 0 Z t Z N 0 0 dt exp ½r0 ðt t Þ do F1 ðoÞ exp ½Iot0 ; tX0: ð49Þ ¼ 2r0 N
0
When time integration is carried out first and the symmetry F1 ðoÞ ¼ F1 ðoÞ is regarded, the following integral remains: Z o2 exp½r0 t þ r0 cosðotÞ þ o sinðotÞ gðtÞ ¼ 4r0 doF1 ðoÞ : ð50Þ r20 þ o2 o1 The result of the numerical evaluation is plotted in Fig. 2. In the given case, the area under the gðtÞ curve vanishes, which actually serves as a check of the numerics. This is a consequence of the property F1 ð0Þ ¼ 0; which by (39) implies
0.2
g
0.1
t 10
20
30
40
50
-0.1
Fig. 2. Turbine function gðtÞ corresponding to Fig. 1 for mean wind speed V ¼ 8 m=s: The dimensions of t and g are s and 1=m; respectively.
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
GR ð0Þ ¼ 0; and because of GR ð0Þ ¼ 1=ð2pÞ 2 statement follows.
RN 0
dt gR ðtÞ ¼ 1=ð2pÞ
171
RN 0
dt gðtÞ; the
8. Applications In this section, we will present first numerical results for two rather idealized cases: (1) The shutdown model of Section 3, which offers an exponential dependence on the turbulence intensity Iu : (2) An application of the general theory of Section 4, to the special stationary power curve Ls ðuÞ as given in (6). From meteorological boundary layer theory, we adopt the following spectral density of the wind speed fluctuations [12]: SðoÞ do ¼ ðVIu Þ2
zV 11 do; p 1 þ 33z o 5=3 2pV
Z
N
do SðoÞ ¼ ðVIu Þ2 ;
ð51Þ
0
where z denotes the hub height of the wind turbine. As compared to [10], definition (2) requires an additional factor V 2 : For notational economy, we have denoted Sðf Þ Sðo=ð2pÞÞ-SðoÞ by the same symbol. We choose the values z ¼ 15 m and V ¼ 8 m=s: To work out the explicit dependence on Iu ; we introduce the spectral density at Iu ¼ 1 : S1 ðoÞ :¼ SðoÞIu ¼1 : Among the C-functions defined in (26), the following two are independent on time t: 1 C3 ¼ Iu2 V 2 ¼ 32Iu2 ; 2 Z 1 2 N C13 ¼ Iu do S1 ðoÞoGI ðoÞ ¼ 0:246403Iu2 ; 2 0
ð52Þ
where the latter result is obtained by numerical Simpson integration from discrete data fields for S1 and the imaginary part, GI ; of the turbine function in frequency space. Also the remaining functions of (26), C2 ðtÞ; C12 ðtÞ; C23 ðtÞ have the form C ðt; Iu Þ ¼ Iu2 C ðt; Iu ¼ 1Þ: The three fields are numerically worked out at Iu ¼ 1 in time steps of Dt ¼ 0:05 s and are available as stored data fields. To evaluate these functions, we need the real and imaginary parts, GR ðoÞ and GI ðoÞ; of the Fourier transform of the turbine function gðtÞ: The two functions were obtained in frequency steps of Do ¼ 0:01 by numerical integration from gðtÞ as given in (50). Relation (39) served as an efficient test of the numerics. Once more, GR and GI are available as stored data fields. For illustration see Fig. 3. To indicate the relevant frequency domain (support), we plot in Fig. 4 the product Sðf ÞG2 ðf Þ; which occurs in the integrand of the function C2 ðtÞ: We now numerically evaluate the shutdown model of Section 3. The time decay to Ls ¼ 0 is given in (20) with w2 ðtÞ Iu2 C2 ðt; Iu ¼ 1Þ: For Iu ¼ 0 and after time Te ¼ 1=r0 the initial power output Lð0Þ is decayed by exp½1: This decay is slown down
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
172
GR
GR,I
0.1
0.2
0.1
f
0.3
-0.1 Fig. 3. Real and imaginary parts GR and GI ; respectively, of the Fourier transform of the turbine function gðtÞ as a function of the frequency f at mean wind speed V ¼ 8 m=s: The units of f and GR;I are 1=s and s=m; respectively.
0.3
SG
0.2
0.1
0.1
0.05
0.2
f Fig. 4. Product of the spectral density of the wind field and the square of the modulus of the turbine function, SG :¼ Sðf ÞG2 ðf Þ; as a function of the frequency f : Mean wind speed V ¼ 8 m=s and turbulence intensity Iu ¼ 0:5: Ordinate and abscissa are in units of seconds and reciprocal seconds, respectively.
by the factor DðIu ÞX1 defined as follows: DðIu Þ :¼ exp½Iu2 C2 ðTe ; Iu ¼ 1Þ;
Te ¼ 1=r0 ¼ 5 s:
ð53Þ
As is seen from Fig. 5, this factor amounts to a dynamic effect of about 11% at turbulence intensity Iu ¼ 0:5: We now turn, as our second application, to the evaluation of (25) for the stationary power (6) giving rise to the approximation Ls ðV 2C23 ðtÞ þ zÞ ¼ Ls ðV Þ þ l1 ðV Þð2C23 ðtÞ þ zÞ þ l2 ðV Þð2C23 ðtÞ þ zÞ2 :
ð54Þ
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
173
1.1
D
1.08 1.06 1.04 1.02
0.1
0.2
0.3
0.4
0.5
Iu Fig. 5. Slowingdown factor D ¼ exp½w2 ðtÞ of the shutdown model at time t ¼ 1=r0 ¼ 5 s as a function of turbulence intensity Iu : Mean wind speed is V ¼ 8 m=s:
R The z integrals in (25) are of the type dz zn exp½z2 =ð4C3 Þ with n ¼ 0; 1 2 3: Since only the even powers survive, we obtain /PðtÞS ¼ exp½C2 ð2C12 þ r0 ÞLs ðV Þ þ 2l1 exp½C2 ðC13 þ 2C12 C23 r0 Þ 2 2 þ 2C12 C3 þ 2C23 r0 C3 r0 Þ: þ 2l2 exp½C2 ð4C13 C23 4C12 C23
ð55Þ We substitute C -Iu2 Cð1Þ ; where the latter refer to Iu ¼ 1: It is observed that jC2ð1Þ ðtÞjoE0:5: Considering turbulence intensities Iu p0:5; we safely can expand exp½C2 in powers of Iu : Going up to Iu6 and using henceforth the previous notation Cð1Þ -C for simplicity, we obtain /PðtÞS ¼ P0 ðtÞ þ P2 ðtÞIu2 þ P4 ðtÞIu4 þ P6 ðtÞIu6 þ OðIu8 Þ ð56Þ with P0 ðtÞ ¼ r0 Ls ðV Þ;
ð57Þ
P2 ðtÞ ¼ Ls ðV Þð2C12 C2 r0 Þ þ 2l1 ðC13 C23 r0 Þ 2l2 C3 r0 ;
ð58Þ
C22 r0 P4 ðtÞ ¼ Ls ðV Þ 2C12 C2 þ þ 2l1 ðC13 C2 þ 2C12 C23 þ C2 C23 r0 Þ 2 2 r0 þ C2 C3 r0 Þ; þ 2l2 ð4C13 C23 þ 2C12 C3 þ 2C23
ð59Þ
P6 ðtÞ ¼ Ls ðV ÞðC12 C22 C23 r0 Þ þ l1 ðC13 C22 4C12 C2 C23 C22 C23 r0 Þ 2 2 4C12 C2 C3 4C2 C23 r0 C22 C3 r0 Þ: þ l2 ð8C13 C2 C23 8C12 C23
ð60Þ
ARTICLE IN PRESS 174
A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
According to (24) the time-averaged power output has to be calculated from Z t Z T 1 /LðtÞSV ¼ dt dt exp½r0 t/PðtÞSV ; ð61Þ T 0 0 where we indicated by the subscript V that our phase average refers to a given mean wind speed V : We have in mind an averaging time of T ¼ 10 min and assume the relaxation time 1=r0 ¼ 5 s: Within the relatively short time T; the mean speeds V are approximately realized by a Weibull distribution, W ðV Þ;Rin general. This means, the long-term mean power has to be estimated from dV W ðV Þ /LðtÞSV : This application is restricted to the calculation of /LðtÞSV at V ¼ 8 m=s: Henceforth, we will omit the subscript V : Due to the exponential factor in the t integral above, this integral levels off to a constant value for times t > t ; where t is of the order of a couple of relaxation times 1=r0 : We assume Tbt and approximate the average power output as follows: Z t /LðtÞSE dt exp½r0 t/PðtÞS; t > t : ð62Þ 0
In the numerical integrations we take t ¼ 35 s corresponding to 700 discrete time steps. In order to single out the dynamical effect, we have to identify the quasi-stationary contribution /Ls ðV þ vÞS ¼ Ls ðV Þ þ l2 /v2 S ¼ Ls 1=2C3 l2 ;
ð63Þ
which corresponds to the contributions P0 and the last term of P2 ; as given in (57) and (58), respectively. As a check, the t integral from P2 leads to the weak turbulence limit (32) and (33) for Tb1=r0 : For convenience, we normalize the power output by the stationary power Ls ðV Þ: This amounts to setting Ls ¼ 1 in (57)–(60) and to replace l1;2 by l* 1;2 :¼ l1;2 =Ls ðV Þ: Numerical time integration up to t ¼ 35 s leads to the following result for the dynamical effect: /LðtÞS ¼ 1 þ l* 2 V 2 Iu2 þ DLdyn ; DLdyn ¼ 3:73l* 1 Iu2 þ ð1:10l* 1 þ 10:24l* 2 ÞIu4 þ ð0:199l* 1 þ 4:43l* 2 ÞIu6 þ Oð1=ðr0 TÞÞ: ð64Þ Whereas relative numerical errors are of the order 103 ; the main approximation consists in neglecting terms of the order 1=ðr0 TÞ: From Section 7, we infer l* 1 ¼ 0:40 s=m and l* 2 ¼ 0:035 s2 =m2 : In Fig. 6 we plot the relative dynamical effect in percent: DLdyn :¼ 100 DLdyn =/LðtÞS:
ð65Þ
A comparison of our results with those reported, e.g. in Fig. 4 of Ref. [4] appears very difficult by several reasons, in particular in view of our idealized example and the fact that our model and the parameters used are not sufficiently validated. We remark that at given mean speed V ; the quasi-stationary power should increase quadratically with higher turbulence intensity Iu which seems to be not the case in Figs. 4(d) and (e) of [4].
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
175
20
Ddyn
15
10
5
0.1
0.2
0.3
0.4
0.5
Iu Fig. 6. Relative dynamical power output, Ddyn :¼ DLdyn ; in percent as a function of the turbulence intensity Iu for mean wind speed V ¼ 8 m=s:
9. Conclusions The turbine model we have introduced, gives rise to relatively compact analytic expressions for the average power output in stochastic wind fields. It may be an efficient tool to assess the effect of wind speed fluctuations on the average power output of a wind turbine. The model is particularly suitable to examine the dynamic effect of turbulent wind fields which can be described by standard models of boundary layer meteorology. On the other hand, through the time-dependent relaxation function rðtÞ ¼ r1 ðtÞ þ r2 ðtÞ there is ample freedom to implement specific response behaviour, possibly including feedback interaction with the power grid which, however, is not considered in this article. To examine the response to turbulent wind fields, we emphasized the coupling between turbine and change rates of wind speed described through r2 : The other part of the relaxation function, r1 ; was set constant for simplicity. In a simplified shutdown model, wind fluctuations tend to weaken the relaxation towards standstill. Moreover, in the case of a special stationary power curve, we find a dynamic enhancement of the average power output with increasing turbulence intensity. As a partial validation of our model, we see the limit of weak turbulence intensity, which reproduces the dynamic ansatz by Rosen and Sheinman [3]. Verification of our model beyond the weak turbulence limit constitutes an ambitious future task.
Acknowledgements The authors are indebted to the referees for their instructive comments, in particular for bringing our attention to literature which became crucial for the final
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
176
version of the article. One of us, A.R., is thankful to Dr. Nikolaus von der Heydt for a stimulating correspondence.
Appendix A. Method of phase averaging Be f ðAÞ an arbitrary function of the stochastic variable A defined as A¼
N pffiffiffiffiffiffi X Df Ai ;
Ai ¼ ri ðsinð2pfi t þ yi þ gðfi ÞÞ sinðyi þ gðfi ÞÞÞ;
ðA:1Þ
i¼1
where ri has been defined in (11). With the aid of the Fourier representation of the Dirac delta function d; we can write Z N Z N 1 /f ðAÞS ¼ dx f ðxÞ dk exp½Ikx/exp½IkAS: ðA:2Þ 2p N N Furthermore, by the statistical independence of the phases yi i ¼ 1; 2; y; N; and taking into account the smallness of the frequency increment Df ; we obtain in turn /exp½IkAS ¼
N D Y
h pffiffiffiffiffiffi iE exp Ik Df Ai
ðA:3Þ
i¼1
and, because /Ai S ¼ 0; D h pffiffiffiffiffiffi iE pffiffiffiffiffiffi exp Ik Df Ai ¼ 1 Ik Df /Ai S 12 k2 Df /A2i S þ OððDf Þ3=2 Þ ¼ exp½ 12 k2 Df /A2i Sð1 þ OððDf Þ3=2 ÞÞ; /A2i S ¼ r2i ð1 cosð2pfi tÞÞ ¼ 2r2i sin2 ðpfi tÞ;
ðA:4Þ ðA:5Þ
which with f ðxÞ ¼ exp½x leads to the result stated in (21).
Appendix B. Phase averaging of Pðt; t0 Þ The method of Appendix A is extended as follows: Z N Z N Z N 0 /Pðt; t ÞS ¼ dx dy dzðr0 þ xÞ exp½yLs ðV þ zÞ N
N
N
/dðx r2 ðt0 ÞÞdðy Rðt0 Þ þ RðtÞÞdðz vðt0 ÞÞS; /dðx r2 ðt0 ÞÞdðy Rðt0 Þ þ RðtÞÞdðz vðt0 ÞÞS Z N Z N Z N 1 ¼ dk dk dk3 exp½Iðk1 x þ k2 y þ k3 zÞ 1 2 ð2pÞ3 N N N
/exp½Iðk1 r2 ðt0 Þ þ k2 ðRðt0 Þ RðtÞÞ þ k3 vðt0 ÞÞS:
ðB:1Þ
ðB:2Þ
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
The latter average has the factorizing structure " # N D N h Y pffiffiffiffiffiffi iE 1X 2 ¼ exp exp IBi Df /Bi SDf 2 i¼1 i¼1
177
ðB:3Þ
with Bi ¼ k1 Zi cosð2pfi t0 þ yi þ gðfi ÞÞ þ k2 ri ðsinð2pfi t0 þ yi þ gðfi ÞÞ sinð2pfi t þ yi þ gðfi ÞÞÞ þ k3 zi sinð2pfi t0 þ yi Þ;
ðB:4Þ
where the symbols ri ; Zi ; and zi are defined in (11) and (27). It can straightforwardly be shown that (we write t ¼ t t0 ) /B2i S ¼ 12 ðZ2i k12 þ 4r2i sin2 ðpfi tÞk22 þ z2i k32 þ 2Zi ri sinð2pfi tÞk1 k2 r2i cosð2pfi tÞ þ Zi zi sinðgðfi ÞÞk1 k3 sinð2pfi t þ gðfi ÞÞ þ k1 k4 zi ri cosð2pfi t þ gðfi ÞÞ þ 2zi ri fcosðgðfi ÞÞ cosð2pfi t þ gðfi ÞÞgk2 k3 Þ:
ðB:5Þ
With the abbreviations introduced in (26) we obtain /exp½Iðk1 r2 ðt0 Þ þ k2 ðRðt0 Þ RðtÞÞ þ k3 vðt0 ÞÞS ¼ exp½C1 k12 þ C2 k22 þ C3 k32 þ C12 k1 k2 þ C13 k1 k3 þ C23 k2 k3 :
ðB:6Þ
Since / B2i S cannot be positive, the exponent in (B.6) is a negative definite quadratic form in the variables k1 to k3 ; In principle, there could exist submanifolds where the form is zero, but these should be unlikely (of zero probability) in view of the dependence on t and on ri ; Zi ; zi ; gðfi Þ i ¼ 1; 2; y; N: We thus assume that the integrals in (B.2) exist. Rather than diagonalizing the quadratic form (B.6), it turns out to be more convenient to carry out in turn the integrations in the original variables, beginning with k1 and proceeding with k2 to k3 and then with x; y and z in (B.2). All integrals are Gaussian giving rise to more or less involved algebraic expressions which have been manipulated by Mathematica [6] without any complication. The result is stated in (25) and (26).
Appendix C. Proofs to Kramers–Kronig section C.1. Solving integral equation (41) The integral equation (41) is transformed into time space, where the convolution integral becomes the integral of a product. The singularity of the kernel is regularized as follows: 1 xo : 2 xo e þ ðx oÞ2
ðC:1Þ
ARTICLE IN PRESS 178
A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
In the same sense we introduce the following modification of the Heaviside step function: ( exp½et if t > 0; e > 0: ðC:2Þ ye ðtÞ ¼ 0 if to0: Its Fourier transform reads Z N 1 1 e Io Ye ðoÞ ¼ dt exp½et Iot ¼ 2p 0 2p e2 þ o2
ðC:3Þ
with the consequences o ¼ IpfYe ðoÞ Ye ðoÞg 2 e þ o2 and o 1 ¼ 2 2 e þo 2p
Z
ðC:4Þ
N
dt qe ðtÞ exp½Iot;
qe ðtÞ ¼ Ipfye ðtÞ ye ðtÞg:
ðC:5Þ
N
With the notation introduced in Section 6, we obtain in time space Z N Z N o GR ðxÞðx oÞ o dx ¼ dt gR ðtÞqe ðtÞ exp½Iot pr0 N e2 þ ðx oÞ2 2p2 r0 N Z N I d ¼ 2 dt gR ðtÞqe ðtÞ exp½Iot 2p r0 N dt Z N I d ¼ 2 dt ½gR ðtÞqe ðtÞ exp½Iot; 2p r0 N dt
ðC:6Þ
where we assume that gR ðtÞ is bounded for t-7N: The integral equation (41) now becomes a differential equation in time space I d ½gR ðtÞðtÞqe ðtÞ ¼ f1 ðtÞ: gR ðtÞ þ ðC:7Þ pr0 dt Making use of I 0 q ðtÞ ¼ p e
(
e exp½et
if t > 0;
e exp½et
if to0;
ðC:8Þ
we obtain the following linear differential equations: g0R ðtÞ þ sþ ðtÞgR ðtÞ ¼ þr0 exp½et f1 ðtÞ; sþ ¼ e þ r0 exp½et g0R ðtÞ þ s ðtÞgR ðtÞ ¼ r0 exp½et f1 ðtÞ;
for t > 0;
ðC:9Þ
s ¼ e r0 exp½et for to0: ðC:10Þ
With the aid of the abbreviations, Z t r0 S7 :¼ dt0 s7 ¼ 8et þ ðexp½7et 1Þ-7r0 t e 0
ðC:11Þ
the solution of the integral equation in time space becomes in the limit e-0; Z t dt0 exp½r0 ðt t0 Þ f1 ðt0 Þ for t > 0; ðC:12Þ gR ðtÞ ¼ exp½r0 tgR ðþ0Þ þ r0 þ0
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
gR ðtÞ ¼ exp½þr0 tgR ð0Þ r0
Z
179
t
dt0 exp½þr0 ðt t0 Þ f1 ðt0 Þ for to0:
ðC:13Þ
0
The symmetry gR ðtÞ ¼ gR ðtÞ implies the equality of the integration constants gR ðþ0Þ ¼ gR ð0Þ: Moreover, we can replace the lower limits of the integrals unambiguously by zero, because the point t0 ¼ 0 is of measure zero. This leads to the form of the solution as stated in (43). For demonstration, we check the solution gR ðtÞ for the following simple turbine function: gðtÞ ¼ g0 sinðOtÞ exp½kt
if t > 0; gðtÞ ¼ 0 else:
ðC:14Þ
Fourier transform reads g0 O k2 o2 þ O2 ; 2p ðk2 þ ðo þ OÞ2 Þðk2 þ ðo OÞ2 Þ g0 O 2ko GI ðoÞ ¼ : 2p ðk2 þ ðo þ OÞ2 Þðk2 þ ðo OÞ2 Þ
GR ðoÞ ¼
ðC:15Þ
From (39) we immediately obtain F1 ðoÞ ¼
g0 O ðr0 ðk2 o2 þ O2 Þ þ 2ko2 : 2pr0 ðk2 þ ðo OÞ2 Þðk2 þ ðo þ OÞ2 Þ
ðC:16Þ
Now we transform the special functions GR and F1 back into time space using standard contour integration with the result gR ðtÞ ¼
f1 ðtÞ ¼
g0 exp½kjtj sin½Ojtj; 2
ðC:17Þ
g0 exp½kjtj ðO cosðOtÞ k sinðOjtjÞ þ r sinðOjtjÞÞ: 2
In this special case gR ð0Þ ¼ 0: The integral in (C.13) gives Z t g0 r0 dt0 exp½r0 ðt t0 Þ f1 ðt0 Þ ¼ exp½kt sinðOtÞ gR ðtÞ; 2 0
ðC:18Þ
t > 0;
ðC:19Þ
which confirms solution (43) identically for t > 0: The case to0 in (43) follows analogously. C.2. Proof of Eq. (44) Let us consider Z N gðtÞ ¼ doðGR ðoÞ þ IGI ðoÞÞ exp½Iot N
ðC:20Þ
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
180
and replace GI by the Kramers–Kronig relation (40). Carrying out the o integration first, we Z N Z find N I exp½Iot dxGR ðxÞ P do p xo N N ( RN þ N dxGR ðxÞ exp½xt ¼ þgR ðtÞ if t > 0; ¼ ðC:21Þ RN N dxGR ðxÞ exp½xt ¼ gR ðtÞ if to0; which gives the result (44). By comparing (C.15) and (C.16), this general result is confirmed once more by the special example. C.3. Determining the integration constant gR ð0Þ Let us consider the causal function gðtÞ with gðtÞ ¼ 0 for to0; and a jump g0 at t ¼ 0 with g0 ¼ lim gðeÞ: e-þ0
ðC:22Þ
For large o and a small e > 0; we have Z N Z N I d 2p GðoÞ ¼ dt gðtÞ exp½Iot ¼ dt gðtÞ exp½Iot o dt e e Z I N d gðtÞ exp½Iot ¼ dt o e dt Z Z I e I N d gðtÞ exp½Iot: ¼ dt g0 dðtÞ exp½Iot þ dt o e o þe dt ðC:23Þ Applying once more partial integration to the last integral shows that it is of order 1=o2 for large o: Of course, we assume that gðtÞ is sufficiently regular for t > 0: We thus obtain in the limit e-0 g0 lim GðoÞ lim ðGR ðoÞ þ IGI ðoÞÞ ¼ I ; ðC:24Þ o-N o-N 2po which implies that the real part, GR ðoÞ; vanishes faster than with 1=o; and GI ðoÞ- g0 =ð2poÞ: Inserting this asymptotic expression into (39), we find lim F1 ðoÞ ¼ g0 =ð2pr0 Þ:
o-N
ðC:25Þ
On the other hand, we consider the offset of solutions (43), g0R ðtÞ :¼ exp½r0 jtjgR ð0Þ; in frequency space: Z N 1 r0 0 GR ðoÞ :¼ dt g0R ðtÞ exp½Iot ¼ gR ð0Þ 2 ; ðC:26Þ 2p N pðr0 þ o2 Þ which by the Kramers–Kronig relation (40) leads to the corresponding imaginary part o GI0 ðoÞ ¼ gR ð0Þ 2 : ðC:27Þ pðr0 þ o2 Þ
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
181
The left-hand side of (39) then gives rise to 0 GR ðoÞ
o 0 gR ð0Þ G ðoÞ ¼ : r0 I pr0
ðC:28Þ
A comparison with (C.25) gives the result claimed in (45): gR ð0Þ ¼ pr0 lim F1 ðoÞ ¼ g0 =2;
ðC:29Þ
o-N
which, because gðtÞ ¼ 2gR ðtÞ for t > 0; proves that solution (43) for the turbine function gðtÞ is uniquely determined by the Rosen–Sheinman response function G1 ðoÞ:
Appendix D. Response to periodic wind gust If a mechanical model of the wind turbine is available, one can simulate the response to the following wind field vðtÞ ¼ a sinðot þ yÞ:
ðD:1Þ
As will be shown, both the modulus GðoÞ and the phase gðoÞ of the turbine function can be measured, in principle, by special choices of time t and phase y: We examine the limit a-0; up to second order, of LðtÞ as defined in the deterministic equations (22) and (23). Using (8) and (10), we have rðtÞ ¼ r0 þ a s o cosðot þ gðoÞ þ yÞ;
s ¼ 2pGðoÞ;
RðtÞ ¼ a s sinðot þ gðoÞ þ yÞ:
ðD:2Þ
We are interested in the response for times tb1=r0 and thus neglect terms which are of order exp½r0 t as compared to terms of order 1. As a consequence, we can disregard the initial value term proportional to Lð0Þ in (22). From (22), (23), and (50) we obtain up to quadratic order in a; Z t LðtÞ ¼ dt0 exp½r0 ðt t0 Þ RðtÞ þ Rðt0 Þ 0
fLs ðV Þ þ l1 vðt0 Þ þ l2 vðt0 Þ2 þ ?g
¼
Z
ðD:3Þ
t
dt0 exp½r0 ðt t0 Þfr0 Ls ðV Þ þ aM1 ðt0 Þ þ a2 M2 ðt0 Þgð1 þ Oða3 ÞÞ
0
ðD:4Þ with M1 ¼ Ls ðV Þs ðo cosðot0 þ g þ yÞ r0 sinðot þ g þ yÞ þ r0 sinðot0 þ g þ yÞÞ þ l1 r0 sinðot0 þ yÞ;
ðD:5Þ
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
182
M2 ¼
r0 s l2 ð1 cosð2ot0 þ 2yÞ þ l1 ðr0 cosðgÞ o sinðgÞ r0 cosðg þ oðt t0 ÞÞ 2 2 þ r0 cosðg þ 2y þ oðt þ t0 ÞÞ r0 cosðg þ 2y þ 2ot0 Þ þ o sinðg þ 2y s2 Ls ðV Þð2r0 þ r0 cosð2g þ 2y þ 2otÞ þ 2r0 cosðoðt t0 ÞÞ 4 2r0 cosð2g þ 2y þ oðt þ t0 ÞÞ þ r0 cosð2g þ 2y þ 2ot0 Þ þ 2o sinðoðt t0 ÞÞ þ 2ot0 ÞÞ
þ 2o sinð2g þ 2y þ oðt þ t0 ÞÞ 2o sinð2g þ 2y þ 2ot0 ÞÞ:
ðD:6Þ
0
The t -integration in (D.5) is elementary. In the result we choose t ¼ 2pN=o; where N is a sufficiently large integer such that terms proportional to exp½r0 t can be neglected. This means, we run the simulations up to this special time. The result, which now depends on the phases g and y only, can be written with LðyÞ :¼ Lðt ¼ 2pN=o; yÞ: LðyÞ ¼ Ls ðV Þ þ aL1 ðyÞ þ a2 L2 ðyÞ þ Oða3 Þ;
ðD:7Þ
r0 ðo cosðyÞ þ r0 sinðyÞÞ; o2 þ r20
ðD:8Þ
L1 ðyÞ ¼ l1 L2 ðyÞ ¼ l2
2ðr20
þ l1
1 ðr2 þ 4o2 r20 cosð2yÞÞ 2r0 o sinð2yÞ þ o2 Þ 0
2r0 ðr20
so2 ðcosðgÞðr30 þ 4r0 o2 Þ þ o2 Þðr20 þ 4o2 Þ
þ cosðg þ 2yÞðr30 2r0 o2 Þ þ sinðgÞðr20 o 4o3 Þ þ 3r20 o sinðg þ 2yÞÞ:
ðD:9Þ
In the next step, we make runs with phase angle y þ p to eliminate the first order term L1 : D1 ðyÞ :¼ 12 ðLðyÞ þ Lðy þ pÞÞ ¼ Ls ðV Þ þ a2 L2 ðyÞ:
ðD:10Þ
To get rid of terms which are independent of y; we form D0 ðyÞ :¼ 12 ðD1 ðyÞ þ D1 ðy þ p=2ÞÞ
ðD:11Þ
and, after restoring s ¼ 2pGðoÞ; we obtain D2 ðyÞ :¼ D1 ðyÞ D0 ðyÞ r0 ¼ l2 2 ðr0 cosð2yÞ þ 2o sinð2yÞÞ 2ðr0 þ 4o2 Þ þ l1
po2 G ððr2 2o2 Þ cosðg þ 2yÞ ðr20 þ o2 Þ ðr20 þ 4o2 Þ 0
þ 3r0 o sinðg þ 2yÞÞ:
ðD:12Þ
As final step, we choose the special values y1 ¼ 0 and y2 ¼ p=4 and arrive at the following linear equations for the two unknowns x and y: x :¼ GðoÞ cosðgðoÞÞ
and
y :¼ GðoÞ sinðgðoÞÞ
ðD:13Þ
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
183
a1 x þ a2 y ¼ 2ðr20 þ o2 Þðr20 þ 4o2 ÞD2 ð0Þ þ a0 ; b1 x þ b2 y ¼ ðr20 þ o2 Þðr20 þ 4r20 ÞD2 ðp=4Þ þ b0
ðD:14Þ
a1 ¼ 2po2 l1 ðr20 2o2 Þ;
a0 ¼ r20 l2 ðr20 þ o2 Þ;
ðD:15Þ
b0 ¼ r0 ol2 ðr20 þ o2 Þ:
ðD:16Þ
with
b1 ¼ 3pr0 o3 l1 ;
a2 ¼ 6pr0 o2 l1 ;
b2 ¼ po2 l1 ðr20 þ 2o2 Þ;
The determinant D of this equation system is D¼
p2 o4 l21 a0 if l1 ðV Þa0; ðr20 þ o2 Þðr20 þ 4o2 Þ
ðD:17Þ
which shows that the above method to determine both modulus G and phase g of the turbine function in Fourier space is feasible in principle.
References [1] A.D. Hansen, P. S^rensen, F. Blaabjerg, J. Becho, Dynamic modelling of wind farm grid interaction, Wind Eng. 26 (2002) 191–208. [2] J.K. Pedersen, K.O. Helgelsen-Pedersen, N. Kj^lstad Poulsen, V. Akhmatov, A. Hejde Nielsen, Contribution to a dynamic wind turbine model validation from wind farm islanding experiment, Electr. Power Syst. Res. 64 (2003) 41–51. [3] Y. Sheinman, A. Rosen, A dynamic model of the influence of turbulence on the power output of a wind turbine, J. Wind Eng. Ind. Aerodyn. 39 (1992) 329–341. [4] A. Rosen, Y. Sheinman, The average power output of a wind turbine in a turbulent wind, J. Wind Eng. Ind. Aerodyn. 51 (1994) 287–302. [5] M. Shinozuka, C.-M. Jan, Digital simulation of random processes and its applications, J. Sound Vib. 25 (1972) 111–128. [6] S. Wolfram, The Mathematica Book, 4th edition, Wolfram Media, Cambridge University Press, Cambridge, 1999. [7] IEC, Wind turbine generator systems, Part 12: wind turbine power performance testing, IEC 6140012, International Standard, 1st Edition, 1998, pp. 61400–12. [8] T. Burton, D. Sharpe, N. Jenkins, E. Bossanyi, Wind Energy Handbook, Wiley, New York, 2001. [9] L.D. Landau, E.M. Lifschitz, Statistische Physik, Akademie-Verlag, Berlin, 1966. [10] M. Siefert, A. Kittel, R. Friedrich, J. Peinke, On a quantitative method to analyse dynamical and measurement noise, Europhys. Lett. 61 (4) (2003) 466–472. . [11] E. Bottcher, E. Anahua, J. Peinke, The Dynamical Behaviour of Wind Gusts, Proceedings of Global Windpower Conference in Paris 2002, European Wind Energy Association, Brussels, as CD. [12] J.C. Kaimal, J.J. Finnigan, Atmospheric Boundary Layer Flows, Oxford University Press, New York, 1994.
Journal of Wind Engineering and Industrial Aerodynamics 92 (2004) 159–183
A phenomenological model for the dynamic response of wind turbines to turbulent wind Alexander Rauh*, Joachim Peinke Institut fur . Physik, Universitat . Oldenburg, D-26111 Oldenburg, Germany Received 11 December 2002; received in revised form 10 November 2003; accepted 12 November 2003
Abstract To predict the average power output of a wind turbine, a response model is proposed which takes into account: (i) the delayed response to the longitudinal wind speed fluctuations; (ii) a response function of the turbine with arbitrary frequency dependence; (iii) wind fields of arbitrary turbulence intensity. In the limit of low turbulence intensity, the dynamical ansatz as proposed in 1992 by Rosen and Sheinman is reproduced. It is shown, how the response function of the turbine can be obtained from simulation experiments of a specific wind turbine. For two idealized situations the dynamic effect of fluctuating wind is estimated at turbulence intensities 0pIu p0:5: At the special mean wind speed V ¼ 8 m=s; the turbine response function is determined from simulation data published by Sheinman and Rosen in 1992 and 1994. r 2003 Elsevier Ltd. All rights reserved. Keywords: Wind turbine; Power prediction; Dynamic response
1. Introduction The prediction of the average power output of a wind turbine, which is accurate enough for an economic assessment, is still a challenging problem. We will focus in this article on the question of how longitudinal wind speed fluctuations affect the power output through a delayed, frequency-dependent response by the wind turbine. The phenomenological model to be proposed includes an arbitrary, time-dependent response function, which, normally, has to be determined from specific aerodynamic *Corresponding author. Tel.: +49-441-798-3460; fax: +49-441-798-19-3460. E-mail address: [email protected] (A. Rauh). 0167-6105/$ - see front matter r 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.jweia.2003.11.001
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
160
and electric–mechanical modelling of the wind turbine. Our model allows for considering arbitrary turbulence intensities. We have not included a possible feedback interaction to the turbine by the power grid, or interferences of turbines within a wind farm, which were considered in a recent paper by Hansen et al. [1]. In another recent study of dynamical effects by Pedersen et al. [2], the electric– mechanical response of a wind farm within an islanding experiment was emphasized rather than the effect of wind speed fluctuations. Our model generalizes a dynamic ansatz proposed by Rosen and Sheinman [3,4] for the case of weak turbulence intensity. Let us briefly sketch their ansatz. They start from the stationary power curve Ls ðuÞ; where u is the axis parallel component of the wind velocity. The instantaneous speed uðtÞ is decomposed as Z 1 T uðtÞ ¼ V þ vðtÞ with V uðtÞ ¼ dtuðtÞ; ð1Þ T 0 where T is a short-term averaging time, e.g. 10 min: The turbulent part, vðtÞ; is a stochastic variable described, according to Shinozuka and Jan [5], as follows: vðtÞ ¼
N pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X 2Sðfi ÞDf sin ð2pfi t þ yi Þ;
ð2Þ
i¼1
where the frequencies fi ¼ iDf ; i ¼ 1; 2; y; are equidistant, yi is a random phase equally distributed over the interval from 0 to 2p; and for iak the phases yi and yk are statistically independent. We will denote in this article the phase average of a function H as Z 2p 1 dy1 y dyN H: ð3Þ /HS ¼ ð2pÞN 0 Clearly, the mean value /vS ¼ 0; and /v2 S ¼
N X
Sðfi ÞDf :
ð4Þ
i¼1
The function Sðf Þ is the spectral density of the turbulent part v of the longitudinal velocity component, and depends, in general, on the short-term mean speed V and on the hub height of the turbine. The corresponding turbulence intensity, Iu ; is defined through ffi 1 pffiffiffiffiffiffiffiffiffiffiffi /v2 S; Iu ¼ ð5Þ V where the index u refers to the longitudinal component u of the wind velocity. In two applications in Section 8, we will adopt for Sðf Þ a distribution function from boundary layer meteorology. In a next step, Rosen and Sheinman [3,4] consider the case of weak turbulence by the following Taylor series up to second order in v: Ls ðuÞ ¼ Ls ðV Þ þ l1 ðV Þv þ l2 ðV Þv2 :
ð6Þ
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
161
In order to include dynamic effects, they modify the stationary average /Ls ðuÞS ¼ Ls ðV Þ þ l2 /v2 S as follows: Ldyn ðV Þ :¼ Ls ðV Þ þ l2
N X
Df Sðfi ÞG1 ðfi ; V Þ:
ð7Þ
i¼1
The function G1 was derived in [3] from simulations by means of a specific rotorgearbox-generator model. The final power prediction was then established by the ensemble average with respect to the mean velocities V : The above dynamic modification corresponds to a coupling between the turbine, through G1 ; and the wind field, through S: As a matter of fact, our generalization was partly guided by the condition to obtain the form (7) in the limit of weak turbulence. The model, to be defined in Section 2, is a linear differential equation for the instantaneous power output LðtÞ of the turbine, and is characterized by a time-dependent relaxation function rðtÞ: The solution of this equation depends on turbine response and wind field in a nonlinear way. Nevertheless, it turns out that phase averaging with respect to the stochastic wind field (2) can rigorously be reduced to a one-dimensional integral containing a Gaussian distribution function. From the general solution, we will derive the weak turbulence limit and thus achieve a connection of the turbine function gðtÞ of our model to the spectral function G1 ðf Þ defined in (7) according to Sheinman and Rosen [3]. As a nontrivial result, we show that gðtÞ can be uniquely determined from the Rosen–Sheinman function G1 : To this end, we invoke causality for gðtÞ; i.e. gðtÞ ¼ 0 for to0 and exploit a Kramers–Kronig relation between the real and imaginary parts of gðtÞ in frequency space. The necessary proofs and, in particular, the solution of a singular integral equation will be outlined in an appendix together with a simple example for demonstration. In two rather idealized applications, we model the timely shutdown behaviour of a wind turbine and calculate the average power output for the special stationary power curve (6) in dependence of the turbulence intensity 0pIu p0:5: To this end, the turbine response function is estimated from graphical data published in [3,4]. Thereby, we restrict ourselves to the data for the mean wind speed V ¼ 8 m=s: In general, the turbine function gðtÞ should be determined by simulation runs of a specific turbine model. A simple way, within our model, consists in simulating the response to velocity jumps. We will discuss this briefly by an example in Section 3. In frequency space, gðtÞ can be obtained, in principle, by simulating the response to periodic wind gusts of given frequency and deterministic phase. As is shown in Appendix D, one has to choose both time and phase in a proper way, in order to avoid involved relations. The partly involved algebraic manipulations of this study and the numerical evaluations were performed with the aid of Mathematica [6].
ARTICLE IN PRESS 162
A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
2. Response model We denote by LðtÞ the instantaneous power output of the wind turbine at time t: It is assumed to be uniquely related to the instantaneous longitudinal wind speed: t-uðtÞ measured, e.g. at the nacelle, or at a certain distance upwind from the turbine as recommended in [7]. The function Ls ðuðtÞÞ corresponds to the stationary power to which the turbine would relax, if the particular speed uðtÞ was frozen in. The power curve u-Ls ðuÞ is taken to be known. The following relaxation model is adopted: dLðtÞ ¼ rðtÞ½Ls ðuðtÞÞ LðtÞ: ð8Þ dt Apart from mathematical simplicity, the form of this model was motivated by two requirements: (1) it should reproduce the dynamical power response (7) in the limit of weak turbulence, as postulated in [3]; (2) it should be more general in order to include arbitrary turbulence intensities. The relaxation function rðtÞ; in general should depend both on the wind speed uðtÞ ¼ V þ vðtÞ and the time derivative u’ ¼ v’: We assume the form rðtÞ ¼ r1 ðuðtÞÞ þ r2 ð’vÞ:
ð9Þ
The part r1 describes then, e.g. relaxation, at constant wind speed, towards the stationary power Ls ðuÞ from a nonstationary value Lð0ÞaLs ðuÞ: For simplicity we will assume r1 as a constant, i.e. depending on the mean wind speed V only with r1 ðuÞ ¼ r1 ðV Þ ¼ r0 : The function r2 ð’vÞ is chosen within the scheme of linear response theory as follows: Z t r2 ðtÞ ¼ dt0 v’ðt0 Þgðt t0 Þ; ð10Þ N
where g describes the response of the turbine. The response function g; in general, depends, on V and Iu ; as control parameters. Within a particular time interval tA½0; T these parameters are available, in principle, from the preceding short-time period. The function r2 ¼ r2 ð’vÞ is chosen to model the instantaneous interaction of the turbine with the turbulent wind. The turbine response has to be causal, i.e. gðt t0 Þ ¼ 0 if t0 > t: As a consequence, in frequency space we have the product of the Fourier transforms of v’ðt0 Þ and gðt t0 Þ: N X pffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi 2pfi ri Df cos ð2pfi t þ yi þ gðfi ÞÞ; ri ¼ Gðfi Þ 2Sðfi Þ; ð11Þ r2 ðtÞ ¼ i¼0
where the Fourier transform of g is written in terms of modulus Gðf Þ ¼ Gðf Þ and phase gðf Þ ¼ gðf Þ: Z pffiffiffiffiffiffiffi 1 N dt gðtÞ exp ð2pIf tÞ ¼ Gðf Þ exp ðIgðf ÞÞ; I ¼ 1: ð12Þ 2p 0 As a first application, we will study in a later section the effect of large turbulence intensity in a simplified shutdown model. To this end we assume that the stationary
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
power curve depends on the turbulence intensity as follows: ( Ls ðu; 0Þ if Iu pI ; Ls ðuÞ-Ls ðu; Iu Þ ¼ 0 if Iu > I :
163
ð13Þ
3. Some simple examples To get some insight into the model, we consider three simple cases: (i) constant speed u0 with v’ 0 and an initial power Lð0ÞaLs ðu0 Þ; (ii) a constant speed jump from u1 to u2 at time t ¼ 0; i.e. v’ ¼ ðu2 u1 ÞdðtÞ; (iii) an idealized shutdown at time t ¼ 0 and instantaneous power output Lð0Þ in the turbulent wind field (2). In example (i) both r ¼ r0 and Ls ðu0 Þ are constant and the solution of (8) exponentially approaches Ls ðu0 Þ as LðtÞ ¼ Ls ðu0 Þ þ ðLð0Þ Ls ðu0 ÞÞ exp ½r0 t:
ð14Þ
The relaxation time 1=r0 is, e.g., of the order of a few seconds in the case of response measurements of blade root bending moments to sudden changes of the blade pitch angle on the Tjaereborg Turbine [8]. In case (ii) we have rðtÞ ¼ r0 þ ðu2 u1 ÞgðtÞ; Lð0Þ ¼ Ls ðu1 Þ; Ls ðtÞ ¼ Ls ðu2 Þ ¼ const: for t > 0: With the definition Z t # ¼ RðtÞ dt0 rðt0 Þ ¼ r0 t þ RðtÞ;
ð15Þ
RðtÞ ¼ ðu2 u1 Þ
0
Z
t
dt0 gðt0 Þ
ð16Þ
0
the solution of (8) can be written # # LðtÞ ¼ exp½RðtÞLð0Þ þ Ls ðu2 Þ exp½RðtÞ
Z Z
t
# 0 Þ dt0 rðt0 Þ exp ½Rðt
0 t
d # 0 Þ exp ½Rðt dt0 0 # # # 0 Þ 1Þ ¼ exp ½RðtÞLð0Þ þ Ls ðu2 Þ exp ½RðtÞðexp ½Rðt # ¼ Ls ðu2 Þ þ ðLð0Þ Ls ðu2 ÞÞ exp ½RðtÞ: # # ¼ exp½RðtÞLð0Þ þ Ls ðu2 Þ exp½RðtÞ
dt0
ð17Þ
If RðtÞ as defined in (16) is bounded in time, which we assume, then the power LðtÞ exponentially converges towards the stationary power Ls ðu2 Þ once more. In case (iii) we implement the stochastic wind field (2) and redefine the function RðtÞ as follows: Z t N pffiffiffiffiffiffi X Df ri sin ð2pfi t þ yi þ gðfi ÞÞ: dt0 rðt0 Þ ¼ r0 t þ RðtÞ Rð0Þ; RðtÞ ¼ 0
i¼0
ð18Þ
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
164
In our shutdown model the stationary power Ls 0 for t > 0: The solution of (8) therefore reads as follows: LðtÞ ¼ exp ½r0 t RðtÞ þ Rð0Þ Lð0Þ:
ð19Þ
The average with respect to the phases yi can rigorously be carried out. As is shown in Appendix A, the result of phase averaging leads to /LðtÞS ¼ exp ½r0 t þ w2 ðtÞ Lð0Þ;
ð20Þ
where w2 ðtÞ ¼
N X
Df r2i sin2 ðpfi tÞ:
ð21Þ
i¼1
As is seen, the stochastic wind field gives rise to an effectively larger relaxation time, i.e. shutdown is slowed down to some extent. Nevertheless, the shutdown state of the turbine eventually is reached by an exponential law, since w2 ðtÞ is bounded in time.
4. General integration and phase averaging In terms of RðtÞ as defined in (18), the solution of the linear differential equation (8) can be written as follows: Z t LðtÞ ¼ exp ½r0 t RðtÞ þ Rð0Þ Lð0Þ þ dt0 exp ½r0 ðt t0 ÞPðt; t0 Þ; ð22Þ 0
Pðt; t0 Þ ¼ rðt0 Þ exp ½RðtÞ þ Rðt0 Þ Ls ðV þ vðt0 ÞÞ:
ð23Þ
The phase averaging of the initial value term was carried out in example (iii), see (20) and (21), with the result that this term exponentially decays to zero. Therefore, we will omit this term henceforth, and consider Z t LðtÞ ¼ dt0 exp ½r0 ðt t0 ÞPðt; t0 Þ: ð24Þ 0
As is shown in Appendix B, we obtain the following phase average (we set t ¼ t t0 )): sffiffiffiffiffiffiffiffiffiffiffi 1 /PðtÞS ¼ exp ½C2 ðtÞ 4pC3 Z N C13 1 2
dz r0 2C12 ðtÞ þ z exp z 4C3 C3 N
Ls ðV 2C23 ðtÞ þ zÞ;
ð25Þ
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
165
where C1 ¼
N 1X Df Z2i ; 4 i¼1
C2 ðtÞ ¼
N X
N 1X Df Zi ri sin ð2pfi tÞ; C12 ðtÞ ¼ 4 i¼1
C23 ðtÞ ¼
Df r2i sin2 ðpfi tÞ;
C3 ¼
i¼1
C13 ¼
N 1X Df z2i ; 4 i¼1
N 1X Df Zi zi sin ðgðfi ÞÞ; 4 i¼1
N 1X Df zi ri ðcosðgðfi ÞÞ cos ð2pfi t þ gðfi ÞÞÞ 4 i¼1
ð26Þ
with Zi ¼ 2pfi ri ;
zi ¼
pffiffiffiffiffiffiffiffiffiffiffiffi 2Sðfi Þ;
ð27Þ
the magnitude C1 appears in an intermediate step of Appendix B. Since the phase averaged function Pðt; t0 Þ depends on the time difference only, from (24) we obtain Z t /LðtÞS ¼ dt exp ½r0 t /PðtÞS: ð28Þ 0
It should be noticed that the half-width of the Gaussian exponent in the integrand of (25), 1=ð4C3 Þ ¼ V 2 =ð2Iu2 Þ; depends on mean values of the wind speed distribution only. Furthermore, C2 ðtÞ ¼ w2 ðtÞ which occurred in the shutdown example, see (21). Moreover, it is seen that the ensemble average of the power output depends in a nonlinear way on the spectrum both of the wind gusts and the turbine response function in general.
5. The limit of weak turbulence In this section we validate the model, to some extent, by deriving from (25) and (28) the ansatz (7) of Rosen and Sheinman. Thereby we will relate the adhoc function G1 ðf ; V Þ; proposed in [3], to the Fourier transform of the turbine function gðtÞ: The strength of the turbulent wind field will be marked by multiplying v with the nondimensional parameter a; which eventually will be set back to the value 1. pffiffiffiffiffiffiffiffiffi Equivalently, Sðf Þ is multiplied by a in expression (25). In the limit a-0; the C coefficients defined in (26), are proportional to a2 : In view of the exponential term of the z-integrand of (25), z is effectively of order a and is scaled by this factor. Referring to (6), we can write for part of the integrand in (25) hðzÞ :¼ exp ½a2 C2 ðr0 2a2 C12 þ azC13 =C3 ÞLs ðV 2a2 C23 þ azÞ ¼ r0 Ls ðV Þ þ azðC13 Ls ðV Þ þ C3 r0 l1 Þ=C3 þ a2 fLs ðV Þð2C12 þ C2 r0 Þ 2C23 r0 l1 þ z2 ðr0 l2 þ C13 =C3 l1 Þg þ Oða3 Þ:
ð29Þ
ARTICLE IN PRESS 166
A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
The z-integration leads to the result, up to terms of order a3 : sffiffiffiffiffiffiffiffiffiffiffi Z N 1 1 2 dz hðzÞ exp z ¼ r0 Ls ðV Þ /PðtÞS ¼ 4pC3 N 4C3 þ a2 f2C12 ðtÞLs ðV ÞÞ C2 ðtÞr0 Ls ðV Þ þ 2C13 l1 2C23 ðtÞr0 l1 þ 2C3 r0 l2 g:
ð30Þ
The time integration, according to (28), amounts to integrals of the type Z t m cos ðnÞ þ r0 sin ðnÞ dt exp ½r0 t sin ðmt þ nÞ ¼ m2 þ r20 0 m cos ðmt þ nÞ þ r0 sin ðmt þ nÞ exp½r0 t : m2 þ r20 ð31Þ In view of the time averaging of LðtÞ we will neglect the exponentially decaying terms, which will cause errors of the order 1=ðr0 TÞ; where T is the averaging time. Restoring the ri ; Zi ; zi from (11) and (27) and restoring a ¼ 1; we arrive at the following weak turbulence limit (up to terms proportional to 1=ðr0 TÞ): /LðtÞS ¼ Ls ðV Þ þ l2
N X
Df Sðfi Þ þ DLdyn
ð32Þ
i¼1
with DLdyn ¼ l1
N X
Df
i¼1
4p2 fi2 Sðfi ÞGðfi Þ 2pfi cos ðgðf ÞÞ sin ðgðf ÞÞ : i i r0 4p2 fi2 þ r20
ð33Þ
The symbol DLdyn denotes the average excess power as compared to the quasistationary average /Ls ðV þ vðtÞÞS ¼ Ls ðV Þ þ l2 /vðtÞ2 S þ Oðv3 Þ; it can be positive or negative in principle. When expressions (32) and (33) are compared with (7), we obtain the function G1 of Rosen and Sheinman [3] (up to normalization by Ls ðV Þ) in the form G* 1 ðf ; V Þ l2 G1 ðf ; V Þ ¼ l2 þ DLdyn : ð34Þ A slightly more compact form can be obtained by introducing the angle Yðf Þ defined by tan ðYÞ ¼ 2pf =r0 : ð2pf Þ2 G* 1 ðf ; V Þ ¼ l2 þ l1 Gðf Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cos ðgðf Þ þ Yðf ÞÞ: r0 4p2 f 2 þ r20
ð35Þ
The coefficients l1;2 depend on the mean wind speed V : Since l2 ðV Þ may have zeros, we prefer to replace G1 by G* 1 ; and write the phase and time averaged power output in the limit of weak turbulence intensity as follows: /LðtÞS ¼ Ls ðV Þ þ
N X i¼1
Df Sðfi ÞG* l ðfi ; V Þ:
ð36Þ
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
167
As is seen from (35), one has optimal phase adaption of the turbine for gðf Þ ¼ Yðf Þ which gives rise to a positive excess power DLdyn in the weak turbulence limit.
6. Kramers–Kronig analysis For later applications we will derive a specific turbine function gðtÞ from the spectrum G1 ðf Þ; which is graphically represented in Fig. 2 of Ref. [3]. For this purpose, we will give the connection between a general response spectrum G1 ðf Þ and the underlying turbine function gðtÞ in the following. For convenience, we will replace f by the angular frequency o ¼ 2pf : In frequency space we denote the real and imaginary parts of the turbine function by GR and GI ; respectively: Z N 1 GR ðoÞ þ IGI ðoÞ GðoÞ ðcosðoÞ þ I sinðoÞÞ ¼ dtgðtÞ exp ½Iot: ð37Þ 2p N Furthermore, in view of (33) and (34), we introduce the function F1 ; which is related to G1 ; as follows: G1 ðoÞ ¼ 1 þ
l1 o2 F1 ðoÞ: l2 o2 þ r20
ð38Þ
From (33) and (34) we obtain the following connection of F1 to the turbine function: o ð39Þ GR ðoÞ GI ðoÞ ¼ F1 ðoÞ: r0 Because the turbine function gðtÞ is causal with gðtÞ ¼ 0 for to0; the imaginary part GI is related to GR by the Kramers–Kronig relation [9] Z N 1 GR ðxÞ ; ð40Þ GI ðoÞ ¼ þ P dx p xo N where P denotes the principal value; there is a ‘‘þ’’—sign in (40) since, as compared to [9], our definition (37) of the Fourier transform differs through o- o: We thus arrive at the following integral equation for the real part GR of the turbine function: Z N o1 GR ðxÞ P ¼ F1 ðoÞ: GR ðoÞ dx ð41Þ r0 p xo N The above integral equation is solved in Appendix C.1. using time space. The corresponding magnitudes in time space are denoted by lower case letters: GR ðoÞ-gR ðtÞ; F1 ðoÞ-f1 ðtÞ: Because the turbine function gðtÞ is real, we have the symmetries: GR ðoÞ ¼ GR ðoÞ;
GI ðoÞ ¼ GI ðoÞ;
gR ðtÞ ¼ gR ðtÞ;
f1 ðtÞ ¼ f1 ðtÞ; ð42Þ
ARTICLE IN PRESS 168
A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
where the last relation is a consequence of (41). Since gðtÞ is not analytic at t ¼ 0; the solution of (41) has to be written separately for t > 0 and to0: Z t gR ðtÞ ¼ exp ½r0 tgR ð0Þ þ r0 dt0 exp ½r0 ðt t0 Þf1 ðt0 Þ for t > 0; 0 Z t dt0 exp ½þr0 ðt t0 Þf1 ðt0 Þ for to0: ð43Þ gR ðtÞ ¼ exp ½þr0 tgR ð0Þ r0 0
It is easily seen, that the above solution obeys the symmetry stated in (42). The turbine function g; we are looking for, is now simply obtained from the following relation, which is generally proved in Appendix C.2.: ( 2gR ðtÞ if t > 0; gðtÞ ¼ ð44Þ 0 if to0: The still open constant gR ð0Þ in solution (43) remains to be discussed. Clearly, if the turbine function is continuous at t ¼ 0; then gðt ¼ 0Þ ¼ 0; and because of (44) the constant gR ð0Þ ¼ 0: On the other hand, a jump of gðtÞ at t ¼ 0 is reflected in the behaviour of the Fourier transform for o-N: As is shown in Appendix C.3, the constant is uniquely determined by F1 as follows: gR ð0Þ ¼ pr0 lim F1 ðoÞ: o-N
ð45Þ
7. Deriving a turbine function In the following we attempt to infer a concrete turbine function gðtÞ from graphical data of Fig. 2 in [3]; we restrict ourself to the case of mean wind speed V ¼ 8 m=s: To derive the function F1 ðoÞ from G1 by means of (38), we have to estimate the parameters l1;2 ðV Þ; for which, once more, we have only graphical data available as presented in [3,4]. The two l parameters, essentially, are the first and second derivatives of the stationary power curve Ls ðV Þ: As should be remarked, a measured power curve is based on 10-min time averages, in general, and differs more or less from the ideal power curve inferred from spatiotemporally constant winds. In our group we have a new method of stochastic analysis to extract the deterministic law of a physical magnitude in a noisy environment [10]. A first application to extract the ideal power curve from given data has been published in [11]. To obtain Ls ðV Þ underlying the Refs. [3,4], we multiply the power coefficient Cp0 ; shown in Fig. 2, of [4], by K V 3 ; to obtain the stationary power Ls ðV Þ ¼ K Cp0 V 3 ; where K is a constant. Within the plot range of Cp0 of [4], namely 6:75pV p12 ms1 ; the power curve can be fitted by a straight line: Ls E constant þð87 m2 =s2 Þ KV ; In this way we infer l1 Eð87 m2 =s2 Þ K: The parameter l2 ; which according to (6) is proportional to the second derivative of Ls ðV Þ; is inferred from the nondimensional function dðV Þ introduced in [3] and plotted there in Fig. 1. The
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
169
5
4
G1
3 2 1
0.01
0.03
0.1
0.3
1.0
f Fig. 1. Approximated Rosen–Sheinman [3] response function G1 ðf Þ for mean wind speed V ¼ 8 m=s: The discrete points have been read off from the original. Apart from the straight line, the curve is an interpolating polynomial of order six in the variable x ¼ xmax ð1 þ 12 Log10 ðf ÞÞ:
connection to l2 is 3 l2 ¼ dðV Þ 2 Ls ðV Þ ¼ 3dðV ÞKCp0 V : V Thus, from the data of [3,4], we obtain approximately the following ratio: l1 87 m2 =s2 E11:5 m=s; E l2 3Cp0 V
ð46Þ
ð47Þ
where at V ¼ 8 m=s we read off dðV ÞE34 and Cp0 ¼ 0:42 from Fig. 1 in [3] and Fig. 2 in [4], respectively. In the next step we read off from Fig. 2 in [3] a set of 12 discrete curve points taking care of the logarithmic scale of the abscissa: ðxðfi Þ; G1 ðfi ÞÞ; i ¼ 1; 2; y; 12 with xðfi Þ ¼ xmax ð1 þ 1=2 Log10 ðfi ÞÞ: The logarithmic plot was reasonably approximated by a sixth order polynomial in the variable xðf Þ: The corresponding interpolation curve, is plotted in Fig. 1 together with the dicrete points read off from the original plot. We now interpolate the function F1 as defined in (38) by a polynomial of order 10 in the variable f ; which replaces the logarithmic variable xðf Þ: As a function of P i o ¼ 2pf the interpolating polynomial 10 a o is characterized by the following i¼0 i coefficients: i ai
0 6:33639
1 2 45:3742 126:113
3 192:997 9
10
0:222442
0:009655:
i
6
7
ai
47:3118
12:9315 2:23907
8
4 183:530
5 114:204
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
170
As is seen from Fig. 1, we set approximately G1 ðf Þ ¼ 1 in the interval 0pf p0:05075: There is, of course, no data basis to extrapolate G1 ðf Þ to large frequencies f : We will set G1 ¼ 1 for large f ; which implies the asymptotic value F1 ¼ 0: According to Fig. 1, the frequency support of G1 and thus of F1 is in the domain f p1 so that time space resolution is not below, say 1 s: It is therefore compatible with the data to assume the turbine function gðt ¼ 0Þ ¼ 0; and as a consequence of (45) F1 ðoÞ ¼ 0 for o-N: We will therefore approximate F1 ðoÞ further by setting F1 ðoÞ ¼ 0 for oXo2 ; where F1 ðo2 Þ ¼ 0: 8 for 0popo1 and o > o2 ; > <0 P10 i ð48Þ F1 ðoÞ ¼ for o1 popo2 ; o1 ¼ 0:3188; o2 ¼ 3:970; i¼0 ai o > : F1 ðoÞ for oo0: With the aid of (43) and (44), and setting gR ð0Þ ¼ 0; we now calculate the turbine functions gðtÞ as follows: Z t gðtÞ ¼ 2r0 dt0 exp ½r0 ðt t0 Þ f1 ðt0 Þ 0 Z t Z N 0 0 dt exp ½r0 ðt t Þ do F1 ðoÞ exp ½Iot0 ; tX0: ð49Þ ¼ 2r0 N
0
When time integration is carried out first and the symmetry F1 ðoÞ ¼ F1 ðoÞ is regarded, the following integral remains: Z o2 exp½r0 t þ r0 cosðotÞ þ o sinðotÞ gðtÞ ¼ 4r0 doF1 ðoÞ : ð50Þ r20 þ o2 o1 The result of the numerical evaluation is plotted in Fig. 2. In the given case, the area under the gðtÞ curve vanishes, which actually serves as a check of the numerics. This is a consequence of the property F1 ð0Þ ¼ 0; which by (39) implies
0.2
g
0.1
t 10
20
30
40
50
-0.1
Fig. 2. Turbine function gðtÞ corresponding to Fig. 1 for mean wind speed V ¼ 8 m=s: The dimensions of t and g are s and 1=m; respectively.
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
GR ð0Þ ¼ 0; and because of GR ð0Þ ¼ 1=ð2pÞ 2 statement follows.
RN 0
dt gR ðtÞ ¼ 1=ð2pÞ
171
RN 0
dt gðtÞ; the
8. Applications In this section, we will present first numerical results for two rather idealized cases: (1) The shutdown model of Section 3, which offers an exponential dependence on the turbulence intensity Iu : (2) An application of the general theory of Section 4, to the special stationary power curve Ls ðuÞ as given in (6). From meteorological boundary layer theory, we adopt the following spectral density of the wind speed fluctuations [12]: SðoÞ do ¼ ðVIu Þ2
zV 11 do; p 1 þ 33z o 5=3 2pV
Z
N
do SðoÞ ¼ ðVIu Þ2 ;
ð51Þ
0
where z denotes the hub height of the wind turbine. As compared to [10], definition (2) requires an additional factor V 2 : For notational economy, we have denoted Sðf Þ Sðo=ð2pÞÞ-SðoÞ by the same symbol. We choose the values z ¼ 15 m and V ¼ 8 m=s: To work out the explicit dependence on Iu ; we introduce the spectral density at Iu ¼ 1 : S1 ðoÞ :¼ SðoÞIu ¼1 : Among the C-functions defined in (26), the following two are independent on time t: 1 C3 ¼ Iu2 V 2 ¼ 32Iu2 ; 2 Z 1 2 N C13 ¼ Iu do S1 ðoÞoGI ðoÞ ¼ 0:246403Iu2 ; 2 0
ð52Þ
where the latter result is obtained by numerical Simpson integration from discrete data fields for S1 and the imaginary part, GI ; of the turbine function in frequency space. Also the remaining functions of (26), C2 ðtÞ; C12 ðtÞ; C23 ðtÞ have the form C ðt; Iu Þ ¼ Iu2 C ðt; Iu ¼ 1Þ: The three fields are numerically worked out at Iu ¼ 1 in time steps of Dt ¼ 0:05 s and are available as stored data fields. To evaluate these functions, we need the real and imaginary parts, GR ðoÞ and GI ðoÞ; of the Fourier transform of the turbine function gðtÞ: The two functions were obtained in frequency steps of Do ¼ 0:01 by numerical integration from gðtÞ as given in (50). Relation (39) served as an efficient test of the numerics. Once more, GR and GI are available as stored data fields. For illustration see Fig. 3. To indicate the relevant frequency domain (support), we plot in Fig. 4 the product Sðf ÞG2 ðf Þ; which occurs in the integrand of the function C2 ðtÞ: We now numerically evaluate the shutdown model of Section 3. The time decay to Ls ¼ 0 is given in (20) with w2 ðtÞ Iu2 C2 ðt; Iu ¼ 1Þ: For Iu ¼ 0 and after time Te ¼ 1=r0 the initial power output Lð0Þ is decayed by exp½1: This decay is slown down
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
172
GR
GR,I
0.1
0.2
0.1
f
0.3
-0.1 Fig. 3. Real and imaginary parts GR and GI ; respectively, of the Fourier transform of the turbine function gðtÞ as a function of the frequency f at mean wind speed V ¼ 8 m=s: The units of f and GR;I are 1=s and s=m; respectively.
0.3
SG
0.2
0.1
0.1
0.05
0.2
f Fig. 4. Product of the spectral density of the wind field and the square of the modulus of the turbine function, SG :¼ Sðf ÞG2 ðf Þ; as a function of the frequency f : Mean wind speed V ¼ 8 m=s and turbulence intensity Iu ¼ 0:5: Ordinate and abscissa are in units of seconds and reciprocal seconds, respectively.
by the factor DðIu ÞX1 defined as follows: DðIu Þ :¼ exp½Iu2 C2 ðTe ; Iu ¼ 1Þ;
Te ¼ 1=r0 ¼ 5 s:
ð53Þ
As is seen from Fig. 5, this factor amounts to a dynamic effect of about 11% at turbulence intensity Iu ¼ 0:5: We now turn, as our second application, to the evaluation of (25) for the stationary power (6) giving rise to the approximation Ls ðV 2C23 ðtÞ þ zÞ ¼ Ls ðV Þ þ l1 ðV Þð2C23 ðtÞ þ zÞ þ l2 ðV Þð2C23 ðtÞ þ zÞ2 :
ð54Þ
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
173
1.1
D
1.08 1.06 1.04 1.02
0.1
0.2
0.3
0.4
0.5
Iu Fig. 5. Slowingdown factor D ¼ exp½w2 ðtÞ of the shutdown model at time t ¼ 1=r0 ¼ 5 s as a function of turbulence intensity Iu : Mean wind speed is V ¼ 8 m=s:
R The z integrals in (25) are of the type dz zn exp½z2 =ð4C3 Þ with n ¼ 0; 1 2 3: Since only the even powers survive, we obtain /PðtÞS ¼ exp½C2 ð2C12 þ r0 ÞLs ðV Þ þ 2l1 exp½C2 ðC13 þ 2C12 C23 r0 Þ 2 2 þ 2C12 C3 þ 2C23 r0 C3 r0 Þ: þ 2l2 exp½C2 ð4C13 C23 4C12 C23
ð55Þ We substitute C -Iu2 Cð1Þ ; where the latter refer to Iu ¼ 1: It is observed that jC2ð1Þ ðtÞjoE0:5: Considering turbulence intensities Iu p0:5; we safely can expand exp½C2 in powers of Iu : Going up to Iu6 and using henceforth the previous notation Cð1Þ -C for simplicity, we obtain /PðtÞS ¼ P0 ðtÞ þ P2 ðtÞIu2 þ P4 ðtÞIu4 þ P6 ðtÞIu6 þ OðIu8 Þ ð56Þ with P0 ðtÞ ¼ r0 Ls ðV Þ;
ð57Þ
P2 ðtÞ ¼ Ls ðV Þð2C12 C2 r0 Þ þ 2l1 ðC13 C23 r0 Þ 2l2 C3 r0 ;
ð58Þ
C22 r0 P4 ðtÞ ¼ Ls ðV Þ 2C12 C2 þ þ 2l1 ðC13 C2 þ 2C12 C23 þ C2 C23 r0 Þ 2 2 r0 þ C2 C3 r0 Þ; þ 2l2 ð4C13 C23 þ 2C12 C3 þ 2C23
ð59Þ
P6 ðtÞ ¼ Ls ðV ÞðC12 C22 C23 r0 Þ þ l1 ðC13 C22 4C12 C2 C23 C22 C23 r0 Þ 2 2 4C12 C2 C3 4C2 C23 r0 C22 C3 r0 Þ: þ l2 ð8C13 C2 C23 8C12 C23
ð60Þ
ARTICLE IN PRESS 174
A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
According to (24) the time-averaged power output has to be calculated from Z t Z T 1 /LðtÞSV ¼ dt dt exp½r0 t/PðtÞSV ; ð61Þ T 0 0 where we indicated by the subscript V that our phase average refers to a given mean wind speed V : We have in mind an averaging time of T ¼ 10 min and assume the relaxation time 1=r0 ¼ 5 s: Within the relatively short time T; the mean speeds V are approximately realized by a Weibull distribution, W ðV Þ;Rin general. This means, the long-term mean power has to be estimated from dV W ðV Þ /LðtÞSV : This application is restricted to the calculation of /LðtÞSV at V ¼ 8 m=s: Henceforth, we will omit the subscript V : Due to the exponential factor in the t integral above, this integral levels off to a constant value for times t > t ; where t is of the order of a couple of relaxation times 1=r0 : We assume Tbt and approximate the average power output as follows: Z t /LðtÞSE dt exp½r0 t/PðtÞS; t > t : ð62Þ 0
In the numerical integrations we take t ¼ 35 s corresponding to 700 discrete time steps. In order to single out the dynamical effect, we have to identify the quasi-stationary contribution /Ls ðV þ vÞS ¼ Ls ðV Þ þ l2 /v2 S ¼ Ls 1=2C3 l2 ;
ð63Þ
which corresponds to the contributions P0 and the last term of P2 ; as given in (57) and (58), respectively. As a check, the t integral from P2 leads to the weak turbulence limit (32) and (33) for Tb1=r0 : For convenience, we normalize the power output by the stationary power Ls ðV Þ: This amounts to setting Ls ¼ 1 in (57)–(60) and to replace l1;2 by l* 1;2 :¼ l1;2 =Ls ðV Þ: Numerical time integration up to t ¼ 35 s leads to the following result for the dynamical effect: /LðtÞS ¼ 1 þ l* 2 V 2 Iu2 þ DLdyn ; DLdyn ¼ 3:73l* 1 Iu2 þ ð1:10l* 1 þ 10:24l* 2 ÞIu4 þ ð0:199l* 1 þ 4:43l* 2 ÞIu6 þ Oð1=ðr0 TÞÞ: ð64Þ Whereas relative numerical errors are of the order 103 ; the main approximation consists in neglecting terms of the order 1=ðr0 TÞ: From Section 7, we infer l* 1 ¼ 0:40 s=m and l* 2 ¼ 0:035 s2 =m2 : In Fig. 6 we plot the relative dynamical effect in percent: DLdyn :¼ 100 DLdyn =/LðtÞS:
ð65Þ
A comparison of our results with those reported, e.g. in Fig. 4 of Ref. [4] appears very difficult by several reasons, in particular in view of our idealized example and the fact that our model and the parameters used are not sufficiently validated. We remark that at given mean speed V ; the quasi-stationary power should increase quadratically with higher turbulence intensity Iu which seems to be not the case in Figs. 4(d) and (e) of [4].
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
175
20
Ddyn
15
10
5
0.1
0.2
0.3
0.4
0.5
Iu Fig. 6. Relative dynamical power output, Ddyn :¼ DLdyn ; in percent as a function of the turbulence intensity Iu for mean wind speed V ¼ 8 m=s:
9. Conclusions The turbine model we have introduced, gives rise to relatively compact analytic expressions for the average power output in stochastic wind fields. It may be an efficient tool to assess the effect of wind speed fluctuations on the average power output of a wind turbine. The model is particularly suitable to examine the dynamic effect of turbulent wind fields which can be described by standard models of boundary layer meteorology. On the other hand, through the time-dependent relaxation function rðtÞ ¼ r1 ðtÞ þ r2 ðtÞ there is ample freedom to implement specific response behaviour, possibly including feedback interaction with the power grid which, however, is not considered in this article. To examine the response to turbulent wind fields, we emphasized the coupling between turbine and change rates of wind speed described through r2 : The other part of the relaxation function, r1 ; was set constant for simplicity. In a simplified shutdown model, wind fluctuations tend to weaken the relaxation towards standstill. Moreover, in the case of a special stationary power curve, we find a dynamic enhancement of the average power output with increasing turbulence intensity. As a partial validation of our model, we see the limit of weak turbulence intensity, which reproduces the dynamic ansatz by Rosen and Sheinman [3]. Verification of our model beyond the weak turbulence limit constitutes an ambitious future task.
Acknowledgements The authors are indebted to the referees for their instructive comments, in particular for bringing our attention to literature which became crucial for the final
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
176
version of the article. One of us, A.R., is thankful to Dr. Nikolaus von der Heydt for a stimulating correspondence.
Appendix A. Method of phase averaging Be f ðAÞ an arbitrary function of the stochastic variable A defined as A¼
N pffiffiffiffiffiffi X Df Ai ;
Ai ¼ ri ðsinð2pfi t þ yi þ gðfi ÞÞ sinðyi þ gðfi ÞÞÞ;
ðA:1Þ
i¼1
where ri has been defined in (11). With the aid of the Fourier representation of the Dirac delta function d; we can write Z N Z N 1 /f ðAÞS ¼ dx f ðxÞ dk exp½Ikx/exp½IkAS: ðA:2Þ 2p N N Furthermore, by the statistical independence of the phases yi i ¼ 1; 2; y; N; and taking into account the smallness of the frequency increment Df ; we obtain in turn /exp½IkAS ¼
N D Y
h pffiffiffiffiffiffi iE exp Ik Df Ai
ðA:3Þ
i¼1
and, because /Ai S ¼ 0; D h pffiffiffiffiffiffi iE pffiffiffiffiffiffi exp Ik Df Ai ¼ 1 Ik Df /Ai S 12 k2 Df /A2i S þ OððDf Þ3=2 Þ ¼ exp½ 12 k2 Df /A2i Sð1 þ OððDf Þ3=2 ÞÞ; /A2i S ¼ r2i ð1 cosð2pfi tÞÞ ¼ 2r2i sin2 ðpfi tÞ;
ðA:4Þ ðA:5Þ
which with f ðxÞ ¼ exp½x leads to the result stated in (21).
Appendix B. Phase averaging of Pðt; t0 Þ The method of Appendix A is extended as follows: Z N Z N Z N 0 /Pðt; t ÞS ¼ dx dy dzðr0 þ xÞ exp½yLs ðV þ zÞ N
N
N
/dðx r2 ðt0 ÞÞdðy Rðt0 Þ þ RðtÞÞdðz vðt0 ÞÞS; /dðx r2 ðt0 ÞÞdðy Rðt0 Þ þ RðtÞÞdðz vðt0 ÞÞS Z N Z N Z N 1 ¼ dk dk dk3 exp½Iðk1 x þ k2 y þ k3 zÞ 1 2 ð2pÞ3 N N N
/exp½Iðk1 r2 ðt0 Þ þ k2 ðRðt0 Þ RðtÞÞ þ k3 vðt0 ÞÞS:
ðB:1Þ
ðB:2Þ
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
The latter average has the factorizing structure " # N D N h Y pffiffiffiffiffiffi iE 1X 2 ¼ exp exp IBi Df /Bi SDf 2 i¼1 i¼1
177
ðB:3Þ
with Bi ¼ k1 Zi cosð2pfi t0 þ yi þ gðfi ÞÞ þ k2 ri ðsinð2pfi t0 þ yi þ gðfi ÞÞ sinð2pfi t þ yi þ gðfi ÞÞÞ þ k3 zi sinð2pfi t0 þ yi Þ;
ðB:4Þ
where the symbols ri ; Zi ; and zi are defined in (11) and (27). It can straightforwardly be shown that (we write t ¼ t t0 ) /B2i S ¼ 12 ðZ2i k12 þ 4r2i sin2 ðpfi tÞk22 þ z2i k32 þ 2Zi ri sinð2pfi tÞk1 k2 r2i cosð2pfi tÞ þ Zi zi sinðgðfi ÞÞk1 k3 sinð2pfi t þ gðfi ÞÞ þ k1 k4 zi ri cosð2pfi t þ gðfi ÞÞ þ 2zi ri fcosðgðfi ÞÞ cosð2pfi t þ gðfi ÞÞgk2 k3 Þ:
ðB:5Þ
With the abbreviations introduced in (26) we obtain /exp½Iðk1 r2 ðt0 Þ þ k2 ðRðt0 Þ RðtÞÞ þ k3 vðt0 ÞÞS ¼ exp½C1 k12 þ C2 k22 þ C3 k32 þ C12 k1 k2 þ C13 k1 k3 þ C23 k2 k3 :
ðB:6Þ
Since / B2i S cannot be positive, the exponent in (B.6) is a negative definite quadratic form in the variables k1 to k3 ; In principle, there could exist submanifolds where the form is zero, but these should be unlikely (of zero probability) in view of the dependence on t and on ri ; Zi ; zi ; gðfi Þ i ¼ 1; 2; y; N: We thus assume that the integrals in (B.2) exist. Rather than diagonalizing the quadratic form (B.6), it turns out to be more convenient to carry out in turn the integrations in the original variables, beginning with k1 and proceeding with k2 to k3 and then with x; y and z in (B.2). All integrals are Gaussian giving rise to more or less involved algebraic expressions which have been manipulated by Mathematica [6] without any complication. The result is stated in (25) and (26).
Appendix C. Proofs to Kramers–Kronig section C.1. Solving integral equation (41) The integral equation (41) is transformed into time space, where the convolution integral becomes the integral of a product. The singularity of the kernel is regularized as follows: 1 xo : 2 xo e þ ðx oÞ2
ðC:1Þ
ARTICLE IN PRESS 178
A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
In the same sense we introduce the following modification of the Heaviside step function: ( exp½et if t > 0; e > 0: ðC:2Þ ye ðtÞ ¼ 0 if to0: Its Fourier transform reads Z N 1 1 e Io Ye ðoÞ ¼ dt exp½et Iot ¼ 2p 0 2p e2 þ o2
ðC:3Þ
with the consequences o ¼ IpfYe ðoÞ Ye ðoÞg 2 e þ o2 and o 1 ¼ 2 2 e þo 2p
Z
ðC:4Þ
N
dt qe ðtÞ exp½Iot;
qe ðtÞ ¼ Ipfye ðtÞ ye ðtÞg:
ðC:5Þ
N
With the notation introduced in Section 6, we obtain in time space Z N Z N o GR ðxÞðx oÞ o dx ¼ dt gR ðtÞqe ðtÞ exp½Iot pr0 N e2 þ ðx oÞ2 2p2 r0 N Z N I d ¼ 2 dt gR ðtÞqe ðtÞ exp½Iot 2p r0 N dt Z N I d ¼ 2 dt ½gR ðtÞqe ðtÞ exp½Iot; 2p r0 N dt
ðC:6Þ
where we assume that gR ðtÞ is bounded for t-7N: The integral equation (41) now becomes a differential equation in time space I d ½gR ðtÞðtÞqe ðtÞ ¼ f1 ðtÞ: gR ðtÞ þ ðC:7Þ pr0 dt Making use of I 0 q ðtÞ ¼ p e
(
e exp½et
if t > 0;
e exp½et
if to0;
ðC:8Þ
we obtain the following linear differential equations: g0R ðtÞ þ sþ ðtÞgR ðtÞ ¼ þr0 exp½et f1 ðtÞ; sþ ¼ e þ r0 exp½et g0R ðtÞ þ s ðtÞgR ðtÞ ¼ r0 exp½et f1 ðtÞ;
for t > 0;
ðC:9Þ
s ¼ e r0 exp½et for to0: ðC:10Þ
With the aid of the abbreviations, Z t r0 S7 :¼ dt0 s7 ¼ 8et þ ðexp½7et 1Þ-7r0 t e 0
ðC:11Þ
the solution of the integral equation in time space becomes in the limit e-0; Z t dt0 exp½r0 ðt t0 Þ f1 ðt0 Þ for t > 0; ðC:12Þ gR ðtÞ ¼ exp½r0 tgR ðþ0Þ þ r0 þ0
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
gR ðtÞ ¼ exp½þr0 tgR ð0Þ r0
Z
179
t
dt0 exp½þr0 ðt t0 Þ f1 ðt0 Þ for to0:
ðC:13Þ
0
The symmetry gR ðtÞ ¼ gR ðtÞ implies the equality of the integration constants gR ðþ0Þ ¼ gR ð0Þ: Moreover, we can replace the lower limits of the integrals unambiguously by zero, because the point t0 ¼ 0 is of measure zero. This leads to the form of the solution as stated in (43). For demonstration, we check the solution gR ðtÞ for the following simple turbine function: gðtÞ ¼ g0 sinðOtÞ exp½kt
if t > 0; gðtÞ ¼ 0 else:
ðC:14Þ
Fourier transform reads g0 O k2 o2 þ O2 ; 2p ðk2 þ ðo þ OÞ2 Þðk2 þ ðo OÞ2 Þ g0 O 2ko GI ðoÞ ¼ : 2p ðk2 þ ðo þ OÞ2 Þðk2 þ ðo OÞ2 Þ
GR ðoÞ ¼
ðC:15Þ
From (39) we immediately obtain F1 ðoÞ ¼
g0 O ðr0 ðk2 o2 þ O2 Þ þ 2ko2 : 2pr0 ðk2 þ ðo OÞ2 Þðk2 þ ðo þ OÞ2 Þ
ðC:16Þ
Now we transform the special functions GR and F1 back into time space using standard contour integration with the result gR ðtÞ ¼
f1 ðtÞ ¼
g0 exp½kjtj sin½Ojtj; 2
ðC:17Þ
g0 exp½kjtj ðO cosðOtÞ k sinðOjtjÞ þ r sinðOjtjÞÞ: 2
In this special case gR ð0Þ ¼ 0: The integral in (C.13) gives Z t g0 r0 dt0 exp½r0 ðt t0 Þ f1 ðt0 Þ ¼ exp½kt sinðOtÞ gR ðtÞ; 2 0
ðC:18Þ
t > 0;
ðC:19Þ
which confirms solution (43) identically for t > 0: The case to0 in (43) follows analogously. C.2. Proof of Eq. (44) Let us consider Z N gðtÞ ¼ doðGR ðoÞ þ IGI ðoÞÞ exp½Iot N
ðC:20Þ
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
180
and replace GI by the Kramers–Kronig relation (40). Carrying out the o integration first, we Z N Z find N I exp½Iot dxGR ðxÞ P do p xo N N ( RN þ N dxGR ðxÞ exp½xt ¼ þgR ðtÞ if t > 0; ¼ ðC:21Þ RN N dxGR ðxÞ exp½xt ¼ gR ðtÞ if to0; which gives the result (44). By comparing (C.15) and (C.16), this general result is confirmed once more by the special example. C.3. Determining the integration constant gR ð0Þ Let us consider the causal function gðtÞ with gðtÞ ¼ 0 for to0; and a jump g0 at t ¼ 0 with g0 ¼ lim gðeÞ: e-þ0
ðC:22Þ
For large o and a small e > 0; we have Z N Z N I d 2p GðoÞ ¼ dt gðtÞ exp½Iot ¼ dt gðtÞ exp½Iot o dt e e Z I N d gðtÞ exp½Iot ¼ dt o e dt Z Z I e I N d gðtÞ exp½Iot: ¼ dt g0 dðtÞ exp½Iot þ dt o e o þe dt ðC:23Þ Applying once more partial integration to the last integral shows that it is of order 1=o2 for large o: Of course, we assume that gðtÞ is sufficiently regular for t > 0: We thus obtain in the limit e-0 g0 lim GðoÞ lim ðGR ðoÞ þ IGI ðoÞÞ ¼ I ; ðC:24Þ o-N o-N 2po which implies that the real part, GR ðoÞ; vanishes faster than with 1=o; and GI ðoÞ- g0 =ð2poÞ: Inserting this asymptotic expression into (39), we find lim F1 ðoÞ ¼ g0 =ð2pr0 Þ:
o-N
ðC:25Þ
On the other hand, we consider the offset of solutions (43), g0R ðtÞ :¼ exp½r0 jtjgR ð0Þ; in frequency space: Z N 1 r0 0 GR ðoÞ :¼ dt g0R ðtÞ exp½Iot ¼ gR ð0Þ 2 ; ðC:26Þ 2p N pðr0 þ o2 Þ which by the Kramers–Kronig relation (40) leads to the corresponding imaginary part o GI0 ðoÞ ¼ gR ð0Þ 2 : ðC:27Þ pðr0 þ o2 Þ
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
181
The left-hand side of (39) then gives rise to 0 GR ðoÞ
o 0 gR ð0Þ G ðoÞ ¼ : r0 I pr0
ðC:28Þ
A comparison with (C.25) gives the result claimed in (45): gR ð0Þ ¼ pr0 lim F1 ðoÞ ¼ g0 =2;
ðC:29Þ
o-N
which, because gðtÞ ¼ 2gR ðtÞ for t > 0; proves that solution (43) for the turbine function gðtÞ is uniquely determined by the Rosen–Sheinman response function G1 ðoÞ:
Appendix D. Response to periodic wind gust If a mechanical model of the wind turbine is available, one can simulate the response to the following wind field vðtÞ ¼ a sinðot þ yÞ:
ðD:1Þ
As will be shown, both the modulus GðoÞ and the phase gðoÞ of the turbine function can be measured, in principle, by special choices of time t and phase y: We examine the limit a-0; up to second order, of LðtÞ as defined in the deterministic equations (22) and (23). Using (8) and (10), we have rðtÞ ¼ r0 þ a s o cosðot þ gðoÞ þ yÞ;
s ¼ 2pGðoÞ;
RðtÞ ¼ a s sinðot þ gðoÞ þ yÞ:
ðD:2Þ
We are interested in the response for times tb1=r0 and thus neglect terms which are of order exp½r0 t as compared to terms of order 1. As a consequence, we can disregard the initial value term proportional to Lð0Þ in (22). From (22), (23), and (50) we obtain up to quadratic order in a; Z t LðtÞ ¼ dt0 exp½r0 ðt t0 Þ RðtÞ þ Rðt0 Þ 0
fLs ðV Þ þ l1 vðt0 Þ þ l2 vðt0 Þ2 þ ?g
¼
Z
ðD:3Þ
t
dt0 exp½r0 ðt t0 Þfr0 Ls ðV Þ þ aM1 ðt0 Þ þ a2 M2 ðt0 Þgð1 þ Oða3 ÞÞ
0
ðD:4Þ with M1 ¼ Ls ðV Þs ðo cosðot0 þ g þ yÞ r0 sinðot þ g þ yÞ þ r0 sinðot0 þ g þ yÞÞ þ l1 r0 sinðot0 þ yÞ;
ðD:5Þ
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
182
M2 ¼
r0 s l2 ð1 cosð2ot0 þ 2yÞ þ l1 ðr0 cosðgÞ o sinðgÞ r0 cosðg þ oðt t0 ÞÞ 2 2 þ r0 cosðg þ 2y þ oðt þ t0 ÞÞ r0 cosðg þ 2y þ 2ot0 Þ þ o sinðg þ 2y s2 Ls ðV Þð2r0 þ r0 cosð2g þ 2y þ 2otÞ þ 2r0 cosðoðt t0 ÞÞ 4 2r0 cosð2g þ 2y þ oðt þ t0 ÞÞ þ r0 cosð2g þ 2y þ 2ot0 Þ þ 2o sinðoðt t0 ÞÞ þ 2ot0 ÞÞ
þ 2o sinð2g þ 2y þ oðt þ t0 ÞÞ 2o sinð2g þ 2y þ 2ot0 ÞÞ:
ðD:6Þ
0
The t -integration in (D.5) is elementary. In the result we choose t ¼ 2pN=o; where N is a sufficiently large integer such that terms proportional to exp½r0 t can be neglected. This means, we run the simulations up to this special time. The result, which now depends on the phases g and y only, can be written with LðyÞ :¼ Lðt ¼ 2pN=o; yÞ: LðyÞ ¼ Ls ðV Þ þ aL1 ðyÞ þ a2 L2 ðyÞ þ Oða3 Þ;
ðD:7Þ
r0 ðo cosðyÞ þ r0 sinðyÞÞ; o2 þ r20
ðD:8Þ
L1 ðyÞ ¼ l1 L2 ðyÞ ¼ l2
2ðr20
þ l1
1 ðr2 þ 4o2 r20 cosð2yÞÞ 2r0 o sinð2yÞ þ o2 Þ 0
2r0 ðr20
so2 ðcosðgÞðr30 þ 4r0 o2 Þ þ o2 Þðr20 þ 4o2 Þ
þ cosðg þ 2yÞðr30 2r0 o2 Þ þ sinðgÞðr20 o 4o3 Þ þ 3r20 o sinðg þ 2yÞÞ:
ðD:9Þ
In the next step, we make runs with phase angle y þ p to eliminate the first order term L1 : D1 ðyÞ :¼ 12 ðLðyÞ þ Lðy þ pÞÞ ¼ Ls ðV Þ þ a2 L2 ðyÞ:
ðD:10Þ
To get rid of terms which are independent of y; we form D0 ðyÞ :¼ 12 ðD1 ðyÞ þ D1 ðy þ p=2ÞÞ
ðD:11Þ
and, after restoring s ¼ 2pGðoÞ; we obtain D2 ðyÞ :¼ D1 ðyÞ D0 ðyÞ r0 ¼ l2 2 ðr0 cosð2yÞ þ 2o sinð2yÞÞ 2ðr0 þ 4o2 Þ þ l1
po2 G ððr2 2o2 Þ cosðg þ 2yÞ ðr20 þ o2 Þ ðr20 þ 4o2 Þ 0
þ 3r0 o sinðg þ 2yÞÞ:
ðD:12Þ
As final step, we choose the special values y1 ¼ 0 and y2 ¼ p=4 and arrive at the following linear equations for the two unknowns x and y: x :¼ GðoÞ cosðgðoÞÞ
and
y :¼ GðoÞ sinðgðoÞÞ
ðD:13Þ
ARTICLE IN PRESS A. Rauh, J. Peinke / J. Wind Eng. Ind. Aerodyn. 92 (2004) 159–183
183
a1 x þ a2 y ¼ 2ðr20 þ o2 Þðr20 þ 4o2 ÞD2 ð0Þ þ a0 ; b1 x þ b2 y ¼ ðr20 þ o2 Þðr20 þ 4r20 ÞD2 ðp=4Þ þ b0
ðD:14Þ
a1 ¼ 2po2 l1 ðr20 2o2 Þ;
a0 ¼ r20 l2 ðr20 þ o2 Þ;
ðD:15Þ
b0 ¼ r0 ol2 ðr20 þ o2 Þ:
ðD:16Þ
with
b1 ¼ 3pr0 o3 l1 ;
a2 ¼ 6pr0 o2 l1 ;
b2 ¼ po2 l1 ðr20 þ 2o2 Þ;
The determinant D of this equation system is D¼
p2 o4 l21 a0 if l1 ðV Þa0; ðr20 þ o2 Þðr20 þ 4o2 Þ
ðD:17Þ
which shows that the above method to determine both modulus G and phase g of the turbine function in Fourier space is feasible in principle.
References [1] A.D. Hansen, P. S^rensen, F. Blaabjerg, J. Becho, Dynamic modelling of wind farm grid interaction, Wind Eng. 26 (2002) 191–208. [2] J.K. Pedersen, K.O. Helgelsen-Pedersen, N. Kj^lstad Poulsen, V. Akhmatov, A. Hejde Nielsen, Contribution to a dynamic wind turbine model validation from wind farm islanding experiment, Electr. Power Syst. Res. 64 (2003) 41–51. [3] Y. Sheinman, A. Rosen, A dynamic model of the influence of turbulence on the power output of a wind turbine, J. Wind Eng. Ind. Aerodyn. 39 (1992) 329–341. [4] A. Rosen, Y. Sheinman, The average power output of a wind turbine in a turbulent wind, J. Wind Eng. Ind. Aerodyn. 51 (1994) 287–302. [5] M. Shinozuka, C.-M. Jan, Digital simulation of random processes and its applications, J. Sound Vib. 25 (1972) 111–128. [6] S. Wolfram, The Mathematica Book, 4th edition, Wolfram Media, Cambridge University Press, Cambridge, 1999. [7] IEC, Wind turbine generator systems, Part 12: wind turbine power performance testing, IEC 6140012, International Standard, 1st Edition, 1998, pp. 61400–12. [8] T. Burton, D. Sharpe, N. Jenkins, E. Bossanyi, Wind Energy Handbook, Wiley, New York, 2001. [9] L.D. Landau, E.M. Lifschitz, Statistische Physik, Akademie-Verlag, Berlin, 1966. [10] M. Siefert, A. Kittel, R. Friedrich, J. Peinke, On a quantitative method to analyse dynamical and measurement noise, Europhys. Lett. 61 (4) (2003) 466–472. . [11] E. Bottcher, E. Anahua, J. Peinke, The Dynamical Behaviour of Wind Gusts, Proceedings of Global Windpower Conference in Paris 2002, European Wind Energy Association, Brussels, as CD. [12] J.C. Kaimal, J.J. Finnigan, Atmospheric Boundary Layer Flows, Oxford University Press, New York, 1994.