- Email: [email protected]
Boundary element analysis of unsteady aerodynamics of windmill rotors in the presence of yaw
Journal of Wind Engineering and Industrial Aerodynamics, 45 (1993) 153-173 Elsevier
153
Boundary element analysis of unsteady aerodynamics of windmill rotors in the presence of yaw G. Arsuffi ENEA
C.R.E. Casaccia, Rome, Italy
G. Guj a n d L. M o r i n o Terza Universitd di Roma "La Sapienza", Rome, Italy
(Received April 12, 1991; revised version received October 27, 1992)
Summary In this paper a method for the evaluation of pressure and velocity fields for a horizontalaxis wind turbine in the presence of yaw is presented. The flow is assumed to be potential and incompressible. The formulation used is general and is applicable to unsteady incompressible potential flows around wind turbines, as well as aircraft and ship propellers. In order to assess the importance of the wake geometry a free-wake analysis is included in this paper. Comparisons with experimental results are also included.
1. I n t r o d u c t i o n T h e i n t e r e s t in u n s t e a d y a e r o d y n a m i c s of h o r i z o n t a l - a x i s wind t u r b i n e s ( H A W T ) is i n c r e a s i n g b e c a u s e of two different issues: the t r e n d t o w a r d s m o r e flexible s t r u c t u r e s [1] (in w h i c h u n s t e a d y a e r o d y n a m i c p r o b l e m s are magnified) a n d the use of m o r e s o p h i s t i c a t e d c o n t r o l s t r a t e g i e s [2] (e.g., p o w e r c o n t r o l by stall or y a w angle). In addition, f a t i g u e failures in c u r r e n t - g e n e r a t i o n m a c h i n e s s h o w the need for i m p r o v i n g the m e t h o d s of blade-load p r e d i c t i o n by i n c l u d i n g u n s t e a d y a e r o d y n a m i c effects [3]. T h e o b j e c t i v e of this p a p e r is to p r e s e n t a m e t h o d for the e v a l u a t i o n of H A W T p r e s s u r e and v e l o c i t y fields in u n s t e a d y flow conditions. Specifically, the e m p h a s i s h e r e is on the flow in the p r e s e n c e of y a w b e c a u s e of the p r a c t i c a l i n t e r e s t in H A W T design a n d of the s i m p l i c i t y of s i m u l a t i n g the p h e n o m e n o n
Correspondence to: G. Guj, Dipartimento di Meccanica e Aeronautica, Universita degli Studi di Roma, Via Eudossiana 18, 00184 Rome, Italy.
0167-6105/93/$06.00 © 1993 Elsevier Science Publishers B.V. All rights reserved.
154
(;. Arsuffi et al./ W i n d m i l l rotors itz the presence of y a w
by numerical and experimental techniques. Accordingly, a numerical model for the analysis of three-dimensional unsteady aerodynamic flows is proposed and tested in several situations. The work is divided into two parts. In the first part the boundary integral formulation for three-dimensional, unsteady, incompressible, potential flows [4] is extended to the problem described above. In the second, the formulation is applied to specific problems for assessment of the validity and the accuracy of the method by comparison with experimental results. Two types of comparisons are considered: steady flows in the absence of yaw and unsteady flows in the presence of yaw. In this second case, particular attention is given to the issue of the wake transport because of its dependence upon the yaw angle: the definition of the wake geometry is one of more critical aspects in the aerodynamic analysis of HAWT, because no systematic studies (such as those for helicopter rotors) exist in this field. In the case of axial flow, both numerical and experimental approaches have been used to study this problem [5 7]. On the other hand, in the non-axial flow analysis the approaches used are typically experimental [8 10]. Therefore, a free-wake analysis is included in this paper; the approach used is closely related to that discussed by Morino and Tseng [11]. 2. F o r m u l a t i o n
Here we present the formulation for potential incompressible flows used for the calculations. The details of the derivation are given in Appendix, where the role of the wake in potential flows is emphasized. The equation for the velocity potential, obtained from the continuity equation, is V2qS=0.
(1)
The boundary conditions, in the air frame of reference (i.e., a reference frame fixed with the undisturbed air), are: ( i ) 5¢9/5n = ~ = n " vR on the rotor surface SR (where VRis the velocity of a point on the rotor surface and n the outward normal); (ii) qS=0 at infinity; (iii) potential discontinuity Aq~= constant following a material wake point; (iv) A(Sgb/Sn)=0 on the wake surface Sw; (v) at the trailing edge, A~b on the wake is equal to the potential discontinuity on the body. The corresponding boundary integral representation is (see Appendix) qS(x.,t.)=
.
~nnG-~nn
dS-
wAq~5~ndS'
(2)
where G = - 1/47r IIx - x . II, with x denoting the point of integration. If the point x . is on the surface SR, (2) is an integral equation that allows to solve for ~b on SR, given 5 ¢ / 5 n (from the body boundary condition) and Ag5 (from the wake boundary condition).
G. Arsuffi et al./WindmiU rotors in the presence of yaw
155
This equation is discretized by dividing SR and Sw in quadrilateral elements (a~ and an, respectively), assuming ~b, 0, and A~b to be constant within each element, and using the collocation method. Using M elements on SR and N elements on Sw, one obtains (Morino and Tseng [11]) M
M
C~k= ~
Bk~O~ + ~
m=l
m=l
N
Ck~¢, + Z FknA~bn,
(3)
n=l
where Bk~
j J ~ Gk da,
(4)
Ck,, = -- f f ~ ~ ~Gk -~n da,
(5)
ffo
(6)
Fkn =
--
. ~
da,
with Gk= --1/4x IIX--Xk II, and Ckk is obtained in the limit as x , approaches Xk: this yields (see Morino and Tseng [11]) r~ ~'~kk = 21L 7 - p(0~ ~ ' k k , where p(0~ ~.Fkk is evaluated with the collocation point on the element. When the collocation point lies on the element, the integrals in Eqs. (4) and (5) have a singularity of order [Ix - Xk II- 1 (in particular, ~Gk/~n ----(X -- Xk)" n/4~ [Ix - - X k H3, and ( x - Xk)" n = 0 (]lX--Xkl] 2) since it approaches the distance from the tangent plane). This singularity is integrable, as it may be seen using polar coordinates. Explicit expressions used for the integrals are given in [12] for arbitrary non-planar quadrilateral elements (portions of hyperboloidal paraboloids), which may be used to approximate any surface without introducing gaps. Two types of analysis are considered: prescribed- and free-wake analysis. In the prescribed-wake analysis the wake is assumed to be a helicoidal surface (without skew), with pitch evaluated according to the velocity at the rotor disk as given by the blade-element theory. The unskewed-wake assumption is acceptable for small yaw angles, but, as the yaw angle increases, it becomes questionable. In order to assess the importance of the wake geometry in the case of high yaw angles, we performed a free-wake analysis. In this case, using the air frame of reference, each wake node is moved according to Xk(t + At)= Xk (t)+ Vk(t)ht, where Vk is the velocity of the node XR. For the first few spirals Vk is evaluated by taking the gradient of Eq. (2) and discretizing (near-wake model). For the remaining spirals, Vk is prescribed (intermediatewake model). A source disk is used for the far-wake model. Details are similar to those given in [11] for helicopter rotors. Next, consider the boundary conditions for non-axial flows. The normalwash may be expressed as ~ = n" VR = n" (-- v~ + ¢o x x + ~ x x), where v~ is the freestream velocity, ~ is the angular velocity of the rotor, and ~ is the angular velocity due to the yaw angle ~ (t~ = d6/dt).
(;. Arsuffi et al./Windmill rotors in the presence of yaw
156
In the p r a c t i c a l i m p l e m e n t a t i o n , it is c o n v e n i e n t to express n and x in t e r m s of t h e i r c o m p o n e n t s in the body frame, a n d v,_, ¢0, and ~ in the air frame. In o r d e r to a c c o m p l i s h this, we i n t r o d u c e the following n o t a t i o n s : a lower-case s a n s e r i f boldface letter will d e n o t e c o l u m n m a t r i c e s (e.g., the c o m p o n e n t s of a v e c t o r expressed in a specified f r a m e of reference, e i t h e r body or air). Thus, the s a m e v e c t o r will be identified by different symbols d e p e n d i n g upon the f r a m e of r e f e r e n c e in w h i c h it is considered. In p a r t i c u l a r , the c o m p o n e n t s of x in the body f r a m e will be identified with x, those in the air f r a m e with z. C a p i t a l s a n s e r i f boldface letters will d e n o t e s q u a r e m a t r i c e s (e.g., t r a n s f o r m a tion m a t r i c e s r e l a t i n g the c o l u m n matrices). In o r d e r to o b t a i n the m a t r i x r e l a t i n g x to z, c o n s i d e r the c o o r d i n a t e s y s t e m s (xx, x2, x3) and (zl, z2, z3), as well as an i n t e r m e d i a t e s y s t e m (y~, Y2, Y3), which is r o t a t i n g with a n g u l a r v e l o c i t y Ft with r e s p e c t to (z~, z2, z3), and w i t h r e s p e c t to w h i c h the s y s t e m (x~, x2, x3) r o t a t e s with the a n g u l a r v e l o c i t y ~. C h o o s i n g the Zl -axis to be p a r a l l e l to v , , the y3-axis and the z3-axis to be p a r a l l e l to D, and the x l - a x i s and the Yl-axis to be p a r a l l e l to ~, one o b t a i n s z = U1 y and y=U2x, or
U1 U2x,
z=
with
(7)
[cos_sin,i] [lo
U 1:
sin6
cos6
0
0
,
U2=
o
0
cos~t
-sinogt
0
sin (gt
cos o t
(8)
N o t e t h a t the d e r i v a t i v e s of the m a t r i c e s U1 a n d U2 h a v e the form
dUl dt
dU2 dt
=AIUI=UIAI'
-A2 U2~ U2A2,
(9)
[0 0] i00 01
where
AI=
D
0
0
0
0
0
,
A2=
0
0
--(9
0
vJ
0
(10)
.
Then, s e t t i n g v~ = [v~, 0, 0] and u s i n g Eqs. (7) a n d (9), one o b t a i n s for the normalwash
¢,=nTU]U~ ~-+v~ =nT(U]A1U2x+A2x+U~U~v~)
=
n2
n3
0
x~coscot
- xl sin ~ot
+co
(° t t c°s t) -x3 x2
+v~,
- sin 6 cos cot
.
sin 6 sin ogt
(11)
G. Arsuffi et al./Windmill rotors in the presence of yaw
157
3. N u m e r i c a l r e s u l t s In order to assess its validity and accuracy the proposed formulation has been used for the numerical evaluation of velocity and pressure fields of windmill rotors. First, the case of axial flow is considered in order to assess the validity and the accuracy of the formulation. Next, the flow in the presence of a yaw angle is analyzed and discussed. Finally, a case of nonconstant yaw (i.e., yaw dependent upon time) is examined. The following notations are used: R denotes the rotor radius, )o=~R/V~ the tip speed ratio, 0 the blade pitch angle, AR the aspect ratio. Also, N1 and N2 denote the number of elements on one side of one blade, in chordwise and spanwise directions respectively, Nw the number of wake elements in the azimuth direction, Ns the number of elements per spiral, and finally N~ the number of the free spirals in the free-wake analysis. The number of steps per revolution is always equal to Ns (Ns = 12 was used in all the computations). Two rotors were used in the calculations. The first, used in the experimental work of Guj et al. [14], is a single-bladed rotor, having rectangular planform, with R=0.14 m, AR=4.67, and 0= 12°, with a symmetrical profile having 12% maximum thickness. The second, used in the experimental work of Wagner and Deppe [15,16], is a two-bladed rotor, having tapered and twisted blades, with R=2.0 m, AR =9.97, ).= 10.5, and 0.75 =2 °, with a variable profile. They will be referred to as "rotor A" and "rotor B", respectively. A preliminary study was performed, using "rotor A", in order to study the convergence of the algorithm. The results indicate that a grid with N1 = 6 and N2 = 10 is sufficient for the evaluation of the velocity in the field, whereas N1 = 12 and Nz = 18 were used for the evaluation of the pressure on the blade surface. Next, the cases without yaw and in the presence of yaw are considered. It should be noted that the velocity data of the experimental works of Wagner and Deppe [15,16] and Guj et al. [14] are given as a function of the azimuth angle and are measured at a point fixed in space and located on a plane parallel to the disk: the distance from the disk x3/R and the radial position r/R are specified later. Let us consider first the case without yaw. A comparison with the experimental results of Wagner and Deppe [15,16] is performed (i.e., the results are obtained using "rotor B", described above). Figures 1 and 2 show the dimensionless axial and tangential velocities, respectively, at a point which has radial position r/R= 0.8 and is located: (i) on the plane of rotation for the axial velocity; (ii) upstream of the rotor (x3/R= 0.0925) for the tangential velocity. Next, consider a comparison with the experimental results of Guj et al. [14] ("rotor A", defined above). These give the velocity at a point with r/R=0.357 and located upstream of the rotor plane, with x3/R = 0.0572. The comparison is presented in Figs. 3 and 4, which show the dimensionless axial and tangential velocities, respectively, for 2 = 4.42. The agreement appears to be satisfactory. In both cases the results have the same trend, but there appears to be a vertical shift between experimental and numerical results (this is stronger in the case
G. Arsuffi et al./ Windmill rotors in the presence of yaw
158 13
I
12
8
>.
r
o o 0 0 0 3
11
1 0
0
r~ >
09
~
t.l--
x
08
07
06 Numer ical Exp,er i n ~ k a t 0
5
~
6]
28
i
4@
i
6@
,
i
80
8@
hz*mu[h
(deg)
1 12~
i
T 149
[] 0
i
ii
II
T, 168
18@
Fig. 1. C o m p a r i s o n w i t h [15]: a x i a l v e l o c i t y v e r s u s a z i m u t h a t rlR=0.80, x~lR=O; 5 = 0 < (Nw = 60, N s = 12).
o f " r o t o r B"). Specifically, the axial velocity exhibits a small relative shift t hat could be due to the u n c e r t a i n t y of the blockage effect of the wind tunnel. Another source of disagreement between numerical and experimental results is due to the presence of viscous effects which are not considered in the numerical case. Other possible explanations for the disagreement are discussed in Ref. [16] for the case of " r o t o r B". Next, consider the results in the presence of yaw. First, consider a comparison with the experimental results of Wagner and Deppe [15]. The results are obtained using " r o t o r B", described above. Figure 5 presents the power coefficient as a function of the yaw angle. The trend is captured correctly and the agreement between numerical and experimental results is satisfactory. Next, in order to assess the method on the basis of "local" quantities, consider
159
G. Arsuffi et a l . / W i n d m i l l rotors in the presence of y a w 8 15
8
L
OlO
-iJ
G
o >
-~
j
O
)
c @
[
o ~-
Oi
-o
-oo'5
d
o
O!
0
0 0
\
0 0
o
0
! ~EXl::,~r _ ,mental ,~ol .
-0 1S o
.
20
.
.
40
.
.
C:~
.
0 0
1
- T - T -
80
18e
120
140
1C~
18E~
Azimuth (deg)
Fig. 2. Comparison with [15]: tangential velocity versus azimuth at r/R =0.80, x l / R = 0,09; 5=0:" (Nw=60, Ns=12). a c o m p a r i s o n with the e x p e r i m e n t a l results of Guj et al. [14]. These were o b t a i n e d w i t h o u t a c o n t r o l on the a n g u l a r speed (the r o t o r is allowed to assume the " e q u i l i b r i u m " a n g u l a r speed a l t h o u g h a flywheel was used to minimize a n g u l a r v e l o c i t y fluctuations). Therefore, the e x p e r i m e n t a l results exhibit the s i m u l t a n e o u s effects of two p a r a m e t e r s : the yaw angle, 6, and the tip speed ratio, ;.. These effects are p r e s e n t e d in Fig. 6, which shows the dimensionless axial velocity as a f u n c t i o n of the azimuth angle, for t h r e e different cases: (i) 5 = 0 ° and )~=4.42 (as in Fig. 3), (ii) 5 = 0 ° and 2=3.90, and (iii) 5 = 1 5 ° and 2--3.901 . We m a y observe that, for the case u n d e r consideration, the two
' the values of case (iii) are those adopted in the experimental results, shown in Figs. 7 and 8 discussed later.
16ll
(;. Arsuffi et al.: WipTdmill rotors in the preset~ce r~/.v~tc 1 3
J
I I
i
f
1 2
r 8
{ 1
!
I
J b.
1 (3
0 0 ,g >
( 99 1 1i
0 X <~
of
08
I
o 7
@6
I
I
NLI~EP,IC AL E (PEPZPIENTAL
i i i i i
OB 0
--
0 0
i
I L
~ J i i
i
I 2g
J
i
i
i
]
~
i
2-40
18(O
Az, mutb,
300
(de,_-j ]'
Fig. 3. Comparison with [14]: axial velocity versus azimuth, at r/R=0.36, xl/R=0.057; ~5=0 (N~v=60, Ns= 12). p a r a m e t e r s h a v e c o n t r a s t i n g effects on the p e a k values: as ). decreases, the p e a k v a l u e s b e c o m e m o r e p r o n o u n c e d , w h e r e a s , as 6 increases, the u p p e r p e a k d e c r e a s e s (along w i t h m o s t of the curve), while the l o w e r p e a k r e m a i n s a l m o s t u n c h a n g e d . Thus, the c o m b i n e d effect of 6 and )o is not strong. This is confirmed by the e x p e r i m e n t a l results, p r e s e n t e d in Figs. 7 and 8, w h i c h show the dimensionless axial and t a n g e n t i a l v e l o c i t y as a f u n c t i o n of the a z i m u t h a n g l e for ~f= 1 5 and ), = 3.90. The c o m p a r i s o n with the e x p e r i m e n t a l r e s u l t s is satisfactory. The c o m p a r i s o n b e t w e e n Cp d i s t r i b u t i o n s at different a z i m u t h angles is s h o w n in T a b l e 1. V a r i a t i o n s with the a z i m u t h a l a n g l e are in the r a n g e of some percent, In o r d e r to verify t h a t the code is bag-free we c o m p a r e d the results with those o b t a i n e d using an i n d e p e n d e n t l y - d e v e l o p e d f r e q u e n c y - d o m a i n code based on the s a m e f o r m u l a t i o n [13]. T h e c o m p a r i s o n is p r e s e n t e d in Fig. 9, which shows the v a l u e of d) at the blade e l e m e n t 11 =6. L = 5 (i.e., at the
161
G. Arsuffi et a l . / W i n d m i l l rotors in the presence of yaw 06
ml
I__
NUMERICAL EXPERIHBNTAL
0 0
i
@ @4
3
5 0
@
>
@2
( I]
@ @
J
C O I.--
@0
--0
Q
,
2
O
J
J
i
i
L
6@
i
i
i
i
,
120
i
t
3
J
,
L
l
,
180 AzimuLh
,
i
,
240
i
,
,
J
i
J
J
i
,
31~
(deg)
Fig. 4. Comparison with [14]: tangential velocity versus azimuth at r/R = 0.36, x l / R = 0.057; 5 = 0 (Nw=60, Ns=12).
trailing-edge e l e m e n t w h e r e the section lift coefficient is higher), for 3 = 15 ° and )~= 4.5: not surprisingly, the r e s u l t s are in e x c e l l e n t a g r e e m e n t w i t h e a c h o t h e r (the difference is due to the zero initial c o n d i t i o n s for the t i m e - d o m a i n solution). Next, in o r d e r to verify t h a t the d i s c r e p a n c y b e t w e e n t h e o r e t i c a l and e x p e r i m e n t a l r e s u l t s is not due to the a s s u m p t i o n on the wake, we p e r f o r m e d a f r e e - w a k e analysis. In o r d e r to assess b e t t e r the i m p a c t of the f r e e - w a k e analysis, we c o n s i d e r e d a case t h a t is p a r t i c u l a r l y s e n s i t i v e to the w a k e g e o m e t r y , i.e., we used a h i g h e r v a l u e for the tip speed r a t i o ()~= 14 w i t h 0 = 2°). The results for the w a k e g e o m e t r y are p r e s e n t e d in Fig. 10 w h e r e the sections b e t w e e n a v e r t i c a l p l a n e and the w a k e surface are s h o w n for the first two spirals. T h e cases w i t h 3 = 0 ~ and 3 = 15" are c o m p a r e d . T h o s e for the p r e s s u r e coefficient as a f u n c t i o n of xl/c (for r/R=0.8 and a z i m u t h equal to zero) are p r e s e n t e d in T a b l e 2, for p r e s c r i b e d and free wake, respectively. T h e s e r e s u l t s
162
(;. A r s u / ] i et (ll. ; W i n d m i l l
r o t o r s lt? t h e p r e s e n c e o f y a w
06
I -iJ
04
c ¢
o i ¢
o c) L
¢
o 02
,col
Expe- i ,'t,~--,t:a I
I
09
0
10
29 Yaw a n g l e
0
t 3@
(deg)
Fig. 5. C o m p a r i s o n w i t h [15]: p o w e r c o e f f i c i e n t v e r s u s y a w a n g l e (Nw = 60, Ns = 12).
indicate that e v e n for this case, w h i c h is d e s i g n e d to e m p h a s i z e the free-wake effect, the differences are r e l a t i v e l y small. Therefore, we believe that the n u m e r i c a l results w i t h a prescribed-wake model are an a c c u r a t e r e p r e s e n t a t i o n of the potential-flow s o l u t i o n . As ment i o n e d a b o v e the d i s c r e p a n c y b e t w e e n n u m e r i c a l and e x p e r i m e n t a l results m a y be due to (1) the e x p e r i m e n t a l error (estimated to be of the order of 4% [14]), to (2) the effects of v i s c o s i t y (not i n c l u d e d in the n u m e r i c a l formulation), and to (3) the t r u n c a t i o n error (see, h o w e v e r , the c o n v e r g e n c e c o n s i d e r a t i o n ment i o n e d above). Finally, in order to s h o w the g e n e r a l i t y of the proposed method, c o n s i d e r a case w i t h n o n c o n s t a n t y a w (i.e., w i t h [2 # 0). R o t o r A w a s used w i t h 2 = 14 and 0 = 2 . We a s s u m e d 3 = 0 for t < to (where to c o r r e s p o n d s to 15 r e v o l u t i o n s , in order to r e a c h steady state), ~ = 5 1 = 1 5 ° for t > t l (where tl corresponds to 16
G. Arsuffi et al./Windmill rotors in the presence of yaw I 3
l /i /~
163
] w
$ = 0 °, ~ = 4.42 6 : 0 °, ,~ : 3.90
1 2
6 = 15°~A = 3 . 9 0
.i
8
°
/ 2"/
2
I 0
>
09
-
"~
/"
/
-o x
<
08
I -~'~'~ I"
/
/
~/,/
;/ 07
0 5
J I I L L
0
, L J , I
6Q
i
~
,
12~3
,
i
,
J
~
1~
,
,
,
240
,
J
,
L
,
J
,
i
i
3~D
h z imuLh ( d e g )
Fig. 6. Axial v e l o c i t y versus azimuth at r/R=0.36 and x l / R = 0 . 0 5 7 for different values of 6 and ). (Nw =60, Ns=12).
r e v o l u t i o n s ) , and
6=61sin2
tt-----to)'
for
to
(12)
Three extra r e v o l u t i o n s after tl w e r e i n c l u d e d in the c a l c u l a t i o n s in order to better identify the t r a n s i e n t r e g i o n in the s o l u t i o n . Figure 11 s h o w s ¢ at a g i v e n panel v e r s u s a z i m u t h angle. Figure 12 s h o w s the Cp distributions at t = to, t = to plus o n e r e v o l u t i o n , and after t w e n t y r e v o l u t i o n s w h i c h are suffic i e n t to r e a c h the periodic solution.
(L A r s u ~ et a l . / W i n d m i l l rotors in the presence of y a w
1(;4 13
12
@
I 1
>
I @
0 >
09
× "~
08
@7
@6 MJ'~_RICAL E)~ERI~AL
--
[3
0
85 8
60
1;_>0
! 80
240
380
36@
Az, rau tF, ( d e g )
F i g . 7. Comparison with [14] : axial velocity versus azimuth at r/R = 0.36, x ~/R = 0,057 ; d = 1 5 ( N w = 6 0 , N s = 12).
5. C o n c l u d i n g r e m a r k s A method for prescribed- and free-wake unsteady aerodynamic analysis of rotors in the presence of yaw has been presented. The proposed method is very general and it seems suited for a wide range of unsteady aerodynamics applications. The numerical results obtained for HAWT rotors are in satisfactory agreement with existing experimental results. The free-wake results indicate that wake distortion due to yaw should be taken into account only for high yaw angles: in this case, however, the potential-aerodynamic assumption is questionable due to the presence of significant viscous-flow effects. Because of the importance of the evaluation of blade loads in such high yaw angle conditions, we recommend that future work be directed toward the development of a viscous-flow capability.
165
G. Arsuffi et a l . / W i n d m i l l rotors in the presence of yaw
06
__.1_
l__
~[~ME~IC~_ £XP£R]]'ENTAL
r] 0
r
8 04 ).
(t
3 0
C
>O
-8
Ilnl O, _
02
-'t
0
o
@
0
c o
(
O0
--~
2
i
9
J
L
J
i
i
68
i
i
*
i
i
128
i
i
i
i
I
i
18e
i
,
*
248
,
,
i
,
,
300
,
,
,
J
36~
Az,muLh (deg) Fig. 8. Comparison with [14] : tangential velocity versus azimuth at r/R = 0.36, x~/R = 0.057; = 15° (Nw = 60, Ns = 12).
Acknowledgements This work was partially supported by the Italian Commission for Nuclear a n d A l t e r n a t i v e E n e r g y S o u r c e s ( E N E A ) u n d e r c o n t r a c t No. 13984 t o t h e University of Rome "La Sapienza".
Appendix. Integral formulation for exterior Neumann problem for Laplace equation In this appendix we present an outline of the boundary element method for i n c o m p r e s s i b l e p o t e n t i a l flows. 2 F o r t h e s a k e o f s i m p l i c i t y , w e a l w a y s a s s u m e 2 The presentation of this section follows closely that of Morino and Tseng [12], to which the reader is referred for details.
166
(;. Arsu/fi et el.~ Windmill rotors in the presence o[y¢~l~
Table 1 C~, distribution versus azimuth position for suction (pressure) side
x/c
Reference value
Relative difference (%)
azimuth = 0
90
180
270
0.49 (0.45)
- 0.881 (0.0173)
- 0.3 (-58)
0.8 (125)
0.5 (.-62)
-. 0.39 (-0.29)
0.534 (0.268)
1.5 (3.7)
3.0 (7.5)
- 0.9 (--3.7)
-0.18 (-0.06)
0.401 (0.233)
- 1.7 (-3.0)
3.7 (-6.0)
- 1.7 (3.0)
0.06 (0.18)
0.246 (0.122)
- 2.4 ( - 2.5)
- 6.5 ( - 4.9)
- 2.4 (-2.5)
0.29 (0.38)
0.0712 ( - 0.0384)
4.1 ( - 3.4)
6.2 ( - 6.0)
2.1 (-2.6)
0.45 (0.49)
-0.146 (--0.293)
0.2 (-5.1)
1.4 (-7.5)
1.6 (3.1)
that the fluid and the body are initially at rest. In addition, the frame of r e f e r e n c e is a s s u m e d to b e c o n n e c t e d w i t h t h e a i r , u n l e s s o t h e r w i s e s t a t e d . T h e g o v e r n i n g e q u a t i o n s f o r a n i n c o m p r e s s i b l e i n v i s c i d flow a r e t h e c o n t i n u ity equation and the Euler equations. From these equations one obtains K e l v i n ' s t h e o r e m w h i c h s t a t e s t h a t , t h e c i r c u l a t i o n o v e r a m a t e r i a l c o n t o u r is c o n s t a n t in t i m e . A s m e n t i o n e d a b o v e , t h e flow field is i n i t i a l l y i r r o t a t i o n a l . T h e n , a t t = 0 w e h a v e F = 0 for a l l t h e c o n t o u r s in t h e field, a n d K e l v i n ' s t h e o r e m y i e l d s t h a t I ' = 0 a t a l l t i m e s for a l l t h e m a t e r i a l c o n t o u r s t h a t rem a i n e d i n t h e field b e t w e e n t = 0 a n d t h e c u r r e n t t i m e . T h e r e f o r e , s u c h a field remains irrotational at all times, except for the points that come in contact with the boundary surface, as for these points there do not exist a contour that remains in the fluid region at all times. Indeed, they form a surface called the wake. C o n s i d e r f i r s t t h e c a s e w i t h o u t a w a k e . S i n c e t h e flow is i r r o t a t i o n a l , t h e r e e x i s t s a f u n c t i o n ~b ( v e l o c i t y p o t e n t i a l ) s u c h t h a t
v=V~,
(1)
Combining the continuity equation, V'v=0, V2~b=0
(x o u t s i d e S , ) .
w i t h Eq. (1) o n e o b t a i n s (2)
N e x t , c o n s i d e r t h e b o u n d a r y c o n d i t i o n s . A s t y p i c a l , t h e s u r f a c e o f t h e b o d y , SB, is h e r e a s s u m e d to b e i m p e r m e a b l e ; a c c o r d i n g l y , t h e b o u n d a r y c o n d i t i o n o n SB is (v v , ) ' n = 0 , w h e r e vs is t h e v e l o c i t y o f t h e p o i n t o n SR. C o m b i n i n g w i t h
G. Arsuffi et al./Windmill rotors in the presence of yaw
167
0 010
0 008 0 000 0 004
S ~
0 002
¢ o n
0
-0 000
-0 002
> -0 004 -0 0OG -0 008 - 0 010 8
30{3
600
900
~28e
~S ~
1800
Az irnuEh ( d e g )
Fig. 9. Time and frequency domain solutions [13]: ¢30 (r/R=0.76) versus azimuth ("Rotor A", Nw = 60, Ns = 12). Eq. (1), one o b t a i n s
On = vB' n .
(3)
In addition, in the air f r a m e of r e f e r e n c e used here, the b o u n d a r y c o n d i t i o n at infinity is g i v e n by v = O( IIx II-2), or, in t e r m s of the v e l o c i t y p o t e n t i a l , = O< IIx If - ~).
¢4)
T h e p r o b l e m described by Eqs. (2), (3), a n d (4) is k n o w n as the e x t e r i o r N e u m a n n p r o b l e m for the L a p l a c e equation. Once the a b o v e p r o b l e m h a s b e e n solved, the p r e s s u r e m a y be e v a l u a t e d u s i n g B e r n o u l l i ' s t h e o r e m for p o t e n t i a l i n c o m p r e s s i b l e flows (a first i n t e g r a l of
16~
(;. A r s u f i et ~1.' W i n d m i l l rotors it; the presence of y¢tw
I
8.16
~=0
'
,
'
l
°
5=15 ° 0. t4
f
"
I
0.12
i I
r~
/
t
E
F
k9
/
L,.T,
g ,a.lo
"o o re
i
I
I
0.06
! II
0.04
x ~
0.02 -0.~
I II k I
.
0.00
0.02
x.
.
0.04 Axial
.
.
0.06 po~i~.ion
.
0.88
O. 10
0.12
(m)
Fig. 10. F r e e - w a k e a n a l y s i s : w a k e g e o m e t r y (first t w o spirals) for 6 = 0 A " . N w = 8 4 , N s = 1 2 , Nf.=2).
a n d 6 = 15 ( " R o t o r
Euler's equations) which, in the air frame, is given by 5q5 ~ , 1 1 __ p= p~, 5t + ~ v - + p P
(5)
where co indicates evaluation at infinity. Next, consider the corresponding boundary integral equation formulation, the basis for which is Green's second identity , (fV2g-gV2f)d't'=
--
,/
kl~n-g~n] d,9°,
(6)
where n is the inwardly directed normal. The essence of the method consists in identifying f with ~b, and in choosing for g the so-called fundamental solution of
G. Arsuffi et al./Windmill rotors in the presence of yaw
169
Table 2 C 0 distribution:relative difference (%) b e t w e e n p r e s c r i b e d a n d f r e e - w a k e r e s u l t s
x/c suction side x/c p r e s s u r e side
0.49 -0.6 0.45 24.
- 0.39 0.7 -0.29 1.5
- 0.18 1.0 -0.06 2.1
0.06 1.2 0.18 4.9
0.29 2.8 0.38 - 20.
0.45 0.4 0.49 - 7.2
-8 04
- 0 86
-8 88
¢
-0
18
-8
12
-8
14
-8
16
I
o & >.
0 >
~0
~1
-@ 18 ' , I , , I . . 5@48 5228 54@6 5588 5760 5948 Az,mu(½
. . 6120
. 638@
6480
666@
6848
(deg)
Fig. 11. F r e e - w a k e a n a l y s i s : ~b~o v e r s u s a z i m u t h (r/R=0.76) ( " R o t o r A", N w = 8 4 , N s = 1 2 ,
NF=2).
the Laplace equation, G, which is defined by V 2 G = 5 ( x - x . ) (where 5 is the Dirac delta function), with boundary condition G= O( ]1x II- 1) at infinity. This is the well k n o w n potential for the unit source, G= -1/4~ ]1x - x . ]j. For exterior problems, we apply Eq. (6) to a volume ~/- ouside ,9~s, bound by a spherical surface, ~ R , having radius R and center x . . Taking the limit as
17()
(;. Arsuffi et el.,, W i n d m i l l rotors in the presence of y a w 0
)
Q6
02 n © I
-02
-06
t=tO t:tg * t: Lint-
i-I
I
rov
0 A
-70
-OS
00
OS
CHORD
Fig. 12. Free-wake analysis: cp distribution at r / R = 0.80 for different time steps ("Rotor A", Ns=12, N~ = 2).
Nw=84,
R tends to infinity, using the boundary condition at infinity for ~b(Eq. (4)), and the corresponding one for G, one obtains, using Eq. (2), ~b(x,)=
8~n a - ¢
~nn d,~,
(7)
where n is the outwardly directed normal to ,~B. Equation (7) may be used to obtain the values of ¢ at any point in the field if the values of ¢ and ~ ¢ / 8 n on the boundary surface &~. are known. We will refer to Eq. (7) as the boundary integral representation for the velocity potential. Note that, if x , approaches a point on the surface ,90,, the value of ~b(x,) approaches the value of ~b at that point of 9~B. This yields a compatibility condition, between the values of ¢ and 8¢len on 5%, which is satisfied by any
G. Arsuffi et al./Windmill rotors in the presence of yaw
171
function that satisfies the Laplace equation. In our case ~c~/Snon 5zB is known from Eq. (3) (Neumann boundary condition), and hence this compatibility condition yields a boundary integral equation that may be used to evaluate on 5%. This integral equation and its extensions are the key to the methodology reviewed here (as well as to BIE methods of direct type in general). Once, ¢ on 5fB is known, we may evaluate ¢ anywhere in :t~, using Eq. (7), and hence the pressure p, using Bernoulli's theorem, Eq. (5). Next, we recall that, according to d'Alembert's paradox, if the flow around an object is irrotational everywhere, then the resultant of the forces acting on the body equals zero. For this reason, the case without wake has little or no interest in aeronautics. In this following, we present a potential-flow formulation for flows around bodies capable of producing lift, i.e., a potential-flow formulation with vorticity present in the field. This may be accomplished if we allow for the potential to be discontinuous in the field. Thus, in order to complete the problem, we need boundary conditions for the wake. These are obtained from the principles of conservation of mass and momentum across a surface of discontinuity, and use (V-Vw)'n=O (no penetration) and Ap = 0 (no pressure discontinuity). For potential flows, the first yields
A ~
=o,
(s)
On the other hand, combining the condition Ap = 0 with Bernoulli's theorem, Eq. (5), yields
~(~)1 2 S t"~-2Vl-1 ~t2--1V22=A ( ~ t ) -~-½(vl'~-D2)°(D1-D2)=0'
(9)
or Dw A¢ = O, Dt
(10)
with Dw/Dt=~/~t+ Vw" V, where Vw=½(Vl + v2). Equation (10) implies that the value of A~b remains constant following a wake point, Xw, which is defined as a point having velocity Vw, and equal to the value it had when Xw left the trailing edge. Numerically, this value is obtained from the trailing-edge condition that, at the trailing edge, A¢ on the wake is equal to the potential discontinuity on the body. Next, consider the integral formulation. Let £fBw denote a surface that surrounds the body, ~B, and a thin layer that includes 5fw. Note that Eq. (2) is is valid outside this surface. Thus, Eq. (7) is still valid, as long as 5~Bis replaced with 5fBw. Next, let the two sides of the portion of ~Bw that surrounds 5fw
(;. Arsuffi et al./Windmill rotors in the presence of yaw
172
become infinitesimally close t o / / ' w . In the limit, recalling Eq. (8), one obtains
~'/l~
~nG-~bOn]
d>/-
'tv.
Adp--d'(/'"
(1l)
~,i2
This is the b o u n d a r y integral r e p r e s e n t a t i o n for the potential for flows a r o u n d lifting objects. As before, in the limit as x , tends to SR, the above representation yields a compatibility c o n d i t i o n between ~b and 5¢/5n on ,(/~ and A~b on Y~w. Note t h a t 5~/5n is k n o w n from the impermeability b o u n d a r y c o n d i t i o n on ,~h (Eq. (3)), w h e r e a s A¢ and the wake g e o m e t r y are k n o w n from the preceding time history. Thus, the compatibility c o n d i t i o n yields a b o u n d a r y integral e q u a t i o n for the u n k n o w n ~ on ,~/~. Then, Eq. (11) yields ~, and from this the velocity, v = V ¢ , as well as the pressure (from B e r n o u l l i ' s theorem, Eq. (5)) a n y w h e r e in the field. In addition, from the velocity of the wake points one o b t a i n s the wake g e o m e t r y at t + d t (with A¢ c o n s t a n t following the wake points and equal to the value it had when Xw left the trailing edge) and the process m a y be repeated.
References [1] B. Montgomerie and A. Zdunek, Dynamic response in horizontal axis wind turbines including instationary effects in the stream tube, Proc. European Wind Energy Conference, Hamburg, FRG, 1984. [2] S.E. Mattsson, Modelling and control of large horizontal axis wind power plants, Doctoral Thesis, Dept. of Automatic Control, Lund Institute of Technology, Sweden, 1984. [3[ A. Hemon, A. Zervos and S. Huberson, Numerical computation of unsteady forces on HAWT, Proc. European Wind Energy Conference, Glasgow, UK, 1989. [4] L. Morino, A general theory of unsteady compressible potential aerodynamics, NASA CR-2464. [5] B. Montgomerie, Performance methods (again!), Proc. 2nd Symposium on Aerodynamics of Wind Turbines, IEA R&D WECS Joint Action, Lyngby, DK, 1988. [6] J.C. Gohard, Free wake analysis of wind turbine aerodynamics, ASRL TR-184-14, 1978. [7] P.H. Alfredsson, J.A. Dahlberg and F.A. Bark, Some properties of the wake behind horizontal axis wind turbines, Proc. 3rd International Symposium on Wind Energy Systems, Lyngby, Denmark, pp. 469 84, Aug. 26-29, 1980. [8] M.B. Anderson, D.J. Milborrow and J.N. Ross, Performance and wake measurements on a 3 m horizontal axis wind turbine Comparison of theory, wind tunnel and field test data, Proc. 4th International Symposium on Wind Energy Systems, Stockholm, Sweden, pp. 113 135, Sep. 21-24, 1982. [9] B.R. Clayton and P. Filby, Measured effects of oblique flows and change in pitch angle on performance and wake development of model wind turbine, Proc. 4th BWEA Wind Energy Workshop, Cranfield, Bedford, United Kingdom, Mar. 24 26, 1982. [10] G.T. Atkinson, S.M.B. Wilmshurst and D.M. Wilson, The aerodynamics of rotors in yaw, Proc. 9th BWEA Conference, Edinbourgh, United Kingdom, pp. 71 78, Apr. t 3, 1987. [11] L. Morino and K. Tseng, A general theory of unsteady compressible potential flows with application to airplanes and rotors, Boundary Element Method in Nonlinear Fluid Dynamics, Ch. 6, eds. P.K. Banerjee and L. Morino (Elsevier Applied Science., London, 1990).
G. Arsuffi et al./Windmill rotors in the presence of yaw
173
[12] L. Morino, L.T. Chen and E.O. Suciu, Steady and oscillatory subsonic and supersonic aerodynamics around complex configurations, AIAA J., 13(3) (1975) 368-374. [13] F. Mastroddi and L. Morino, Time and frequency domain aerodynamics for flutter of helicopter rotors in hover, International Symposium on Boundary Element Method, Hartford, CT, USA, 1989. [14] G. Guj, M. Terzitta and G. Arsuffi, Velocity measurements upstream of a windmill rotor model, Wind Eng., 15(5) (1991) 248 260. [15] W.J. Wagner, L. Deppe, Messung des Momentanen Geschwindigkeitsfeldes einer Modellwindturbine sowie der Auftriebsverteilungen der rotierenden B1/itter, DFVLRIB 232-88 J 06, DLR, Institut fur Aeroelastik, GSttingen, FRG. [16] W.J. Wagner, private communication (1990).
153
Boundary element analysis of unsteady aerodynamics of windmill rotors in the presence of yaw G. Arsuffi ENEA
C.R.E. Casaccia, Rome, Italy
G. Guj a n d L. M o r i n o Terza Universitd di Roma "La Sapienza", Rome, Italy
(Received April 12, 1991; revised version received October 27, 1992)
Summary In this paper a method for the evaluation of pressure and velocity fields for a horizontalaxis wind turbine in the presence of yaw is presented. The flow is assumed to be potential and incompressible. The formulation used is general and is applicable to unsteady incompressible potential flows around wind turbines, as well as aircraft and ship propellers. In order to assess the importance of the wake geometry a free-wake analysis is included in this paper. Comparisons with experimental results are also included.
1. I n t r o d u c t i o n T h e i n t e r e s t in u n s t e a d y a e r o d y n a m i c s of h o r i z o n t a l - a x i s wind t u r b i n e s ( H A W T ) is i n c r e a s i n g b e c a u s e of two different issues: the t r e n d t o w a r d s m o r e flexible s t r u c t u r e s [1] (in w h i c h u n s t e a d y a e r o d y n a m i c p r o b l e m s are magnified) a n d the use of m o r e s o p h i s t i c a t e d c o n t r o l s t r a t e g i e s [2] (e.g., p o w e r c o n t r o l by stall or y a w angle). In addition, f a t i g u e failures in c u r r e n t - g e n e r a t i o n m a c h i n e s s h o w the need for i m p r o v i n g the m e t h o d s of blade-load p r e d i c t i o n by i n c l u d i n g u n s t e a d y a e r o d y n a m i c effects [3]. T h e o b j e c t i v e of this p a p e r is to p r e s e n t a m e t h o d for the e v a l u a t i o n of H A W T p r e s s u r e and v e l o c i t y fields in u n s t e a d y flow conditions. Specifically, the e m p h a s i s h e r e is on the flow in the p r e s e n c e of y a w b e c a u s e of the p r a c t i c a l i n t e r e s t in H A W T design a n d of the s i m p l i c i t y of s i m u l a t i n g the p h e n o m e n o n
Correspondence to: G. Guj, Dipartimento di Meccanica e Aeronautica, Universita degli Studi di Roma, Via Eudossiana 18, 00184 Rome, Italy.
0167-6105/93/$06.00 © 1993 Elsevier Science Publishers B.V. All rights reserved.
154
(;. Arsuffi et al./ W i n d m i l l rotors itz the presence of y a w
by numerical and experimental techniques. Accordingly, a numerical model for the analysis of three-dimensional unsteady aerodynamic flows is proposed and tested in several situations. The work is divided into two parts. In the first part the boundary integral formulation for three-dimensional, unsteady, incompressible, potential flows [4] is extended to the problem described above. In the second, the formulation is applied to specific problems for assessment of the validity and the accuracy of the method by comparison with experimental results. Two types of comparisons are considered: steady flows in the absence of yaw and unsteady flows in the presence of yaw. In this second case, particular attention is given to the issue of the wake transport because of its dependence upon the yaw angle: the definition of the wake geometry is one of more critical aspects in the aerodynamic analysis of HAWT, because no systematic studies (such as those for helicopter rotors) exist in this field. In the case of axial flow, both numerical and experimental approaches have been used to study this problem [5 7]. On the other hand, in the non-axial flow analysis the approaches used are typically experimental [8 10]. Therefore, a free-wake analysis is included in this paper; the approach used is closely related to that discussed by Morino and Tseng [11]. 2. F o r m u l a t i o n
Here we present the formulation for potential incompressible flows used for the calculations. The details of the derivation are given in Appendix, where the role of the wake in potential flows is emphasized. The equation for the velocity potential, obtained from the continuity equation, is V2qS=0.
(1)
The boundary conditions, in the air frame of reference (i.e., a reference frame fixed with the undisturbed air), are: ( i ) 5¢9/5n = ~ = n " vR on the rotor surface SR (where VRis the velocity of a point on the rotor surface and n the outward normal); (ii) qS=0 at infinity; (iii) potential discontinuity Aq~= constant following a material wake point; (iv) A(Sgb/Sn)=0 on the wake surface Sw; (v) at the trailing edge, A~b on the wake is equal to the potential discontinuity on the body. The corresponding boundary integral representation is (see Appendix) qS(x.,t.)=
.
~nnG-~nn
dS-
wAq~5~ndS'
(2)
where G = - 1/47r IIx - x . II, with x denoting the point of integration. If the point x . is on the surface SR, (2) is an integral equation that allows to solve for ~b on SR, given 5 ¢ / 5 n (from the body boundary condition) and Ag5 (from the wake boundary condition).
G. Arsuffi et al./WindmiU rotors in the presence of yaw
155
This equation is discretized by dividing SR and Sw in quadrilateral elements (a~ and an, respectively), assuming ~b, 0, and A~b to be constant within each element, and using the collocation method. Using M elements on SR and N elements on Sw, one obtains (Morino and Tseng [11]) M
M
C~k= ~
Bk~O~ + ~
m=l
m=l
N
Ck~¢, + Z FknA~bn,
(3)
n=l
where Bk~
j J ~ Gk da,
(4)
Ck,, = -- f f ~ ~ ~Gk -~n da,
(5)
ffo
(6)
Fkn =
--
. ~
da,
with Gk= --1/4x IIX--Xk II, and Ckk is obtained in the limit as x , approaches Xk: this yields (see Morino and Tseng [11]) r~ ~'~kk = 21L 7 - p(0~ ~ ' k k , where p(0~ ~.Fkk is evaluated with the collocation point on the element. When the collocation point lies on the element, the integrals in Eqs. (4) and (5) have a singularity of order [Ix - Xk II- 1 (in particular, ~Gk/~n ----(X -- Xk)" n/4~ [Ix - - X k H3, and ( x - Xk)" n = 0 (]lX--Xkl] 2) since it approaches the distance from the tangent plane). This singularity is integrable, as it may be seen using polar coordinates. Explicit expressions used for the integrals are given in [12] for arbitrary non-planar quadrilateral elements (portions of hyperboloidal paraboloids), which may be used to approximate any surface without introducing gaps. Two types of analysis are considered: prescribed- and free-wake analysis. In the prescribed-wake analysis the wake is assumed to be a helicoidal surface (without skew), with pitch evaluated according to the velocity at the rotor disk as given by the blade-element theory. The unskewed-wake assumption is acceptable for small yaw angles, but, as the yaw angle increases, it becomes questionable. In order to assess the importance of the wake geometry in the case of high yaw angles, we performed a free-wake analysis. In this case, using the air frame of reference, each wake node is moved according to Xk(t + At)= Xk (t)+ Vk(t)ht, where Vk is the velocity of the node XR. For the first few spirals Vk is evaluated by taking the gradient of Eq. (2) and discretizing (near-wake model). For the remaining spirals, Vk is prescribed (intermediatewake model). A source disk is used for the far-wake model. Details are similar to those given in [11] for helicopter rotors. Next, consider the boundary conditions for non-axial flows. The normalwash may be expressed as ~ = n" VR = n" (-- v~ + ¢o x x + ~ x x), where v~ is the freestream velocity, ~ is the angular velocity of the rotor, and ~ is the angular velocity due to the yaw angle ~ (t~ = d6/dt).
(;. Arsuffi et al./Windmill rotors in the presence of yaw
156
In the p r a c t i c a l i m p l e m e n t a t i o n , it is c o n v e n i e n t to express n and x in t e r m s of t h e i r c o m p o n e n t s in the body frame, a n d v,_, ¢0, and ~ in the air frame. In o r d e r to a c c o m p l i s h this, we i n t r o d u c e the following n o t a t i o n s : a lower-case s a n s e r i f boldface letter will d e n o t e c o l u m n m a t r i c e s (e.g., the c o m p o n e n t s of a v e c t o r expressed in a specified f r a m e of reference, e i t h e r body or air). Thus, the s a m e v e c t o r will be identified by different symbols d e p e n d i n g upon the f r a m e of r e f e r e n c e in w h i c h it is considered. In p a r t i c u l a r , the c o m p o n e n t s of x in the body f r a m e will be identified with x, those in the air f r a m e with z. C a p i t a l s a n s e r i f boldface letters will d e n o t e s q u a r e m a t r i c e s (e.g., t r a n s f o r m a tion m a t r i c e s r e l a t i n g the c o l u m n matrices). In o r d e r to o b t a i n the m a t r i x r e l a t i n g x to z, c o n s i d e r the c o o r d i n a t e s y s t e m s (xx, x2, x3) and (zl, z2, z3), as well as an i n t e r m e d i a t e s y s t e m (y~, Y2, Y3), which is r o t a t i n g with a n g u l a r v e l o c i t y Ft with r e s p e c t to (z~, z2, z3), and w i t h r e s p e c t to w h i c h the s y s t e m (x~, x2, x3) r o t a t e s with the a n g u l a r v e l o c i t y ~. C h o o s i n g the Zl -axis to be p a r a l l e l to v , , the y3-axis and the z3-axis to be p a r a l l e l to D, and the x l - a x i s and the Yl-axis to be p a r a l l e l to ~, one o b t a i n s z = U1 y and y=U2x, or
U1 U2x,
z=
with
(7)
[cos_sin,i] [lo
U 1:
sin6
cos6
0
0
,
U2=
o
0
cos~t
-sinogt
0
sin (gt
cos o t
(8)
N o t e t h a t the d e r i v a t i v e s of the m a t r i c e s U1 a n d U2 h a v e the form
dUl dt
dU2 dt
=AIUI=UIAI'
-A2 U2~ U2A2,
(9)
[0 0] i00 01
where
AI=
D
0
0
0
0
0
,
A2=
0
0
--(9
0
vJ
0
(10)
.
Then, s e t t i n g v~ = [v~, 0, 0] and u s i n g Eqs. (7) a n d (9), one o b t a i n s for the normalwash
¢,=nTU]U~ ~-+v~ =nT(U]A1U2x+A2x+U~U~v~)
=
n2
n3
0
x~coscot
- xl sin ~ot
+co
(° t t c°s t) -x3 x2
+v~,
- sin 6 cos cot
.
sin 6 sin ogt
(11)
G. Arsuffi et al./Windmill rotors in the presence of yaw
157
3. N u m e r i c a l r e s u l t s In order to assess its validity and accuracy the proposed formulation has been used for the numerical evaluation of velocity and pressure fields of windmill rotors. First, the case of axial flow is considered in order to assess the validity and the accuracy of the formulation. Next, the flow in the presence of a yaw angle is analyzed and discussed. Finally, a case of nonconstant yaw (i.e., yaw dependent upon time) is examined. The following notations are used: R denotes the rotor radius, )o=~R/V~ the tip speed ratio, 0 the blade pitch angle, AR the aspect ratio. Also, N1 and N2 denote the number of elements on one side of one blade, in chordwise and spanwise directions respectively, Nw the number of wake elements in the azimuth direction, Ns the number of elements per spiral, and finally N~ the number of the free spirals in the free-wake analysis. The number of steps per revolution is always equal to Ns (Ns = 12 was used in all the computations). Two rotors were used in the calculations. The first, used in the experimental work of Guj et al. [14], is a single-bladed rotor, having rectangular planform, with R=0.14 m, AR=4.67, and 0= 12°, with a symmetrical profile having 12% maximum thickness. The second, used in the experimental work of Wagner and Deppe [15,16], is a two-bladed rotor, having tapered and twisted blades, with R=2.0 m, AR =9.97, ).= 10.5, and 0.75 =2 °, with a variable profile. They will be referred to as "rotor A" and "rotor B", respectively. A preliminary study was performed, using "rotor A", in order to study the convergence of the algorithm. The results indicate that a grid with N1 = 6 and N2 = 10 is sufficient for the evaluation of the velocity in the field, whereas N1 = 12 and Nz = 18 were used for the evaluation of the pressure on the blade surface. Next, the cases without yaw and in the presence of yaw are considered. It should be noted that the velocity data of the experimental works of Wagner and Deppe [15,16] and Guj et al. [14] are given as a function of the azimuth angle and are measured at a point fixed in space and located on a plane parallel to the disk: the distance from the disk x3/R and the radial position r/R are specified later. Let us consider first the case without yaw. A comparison with the experimental results of Wagner and Deppe [15,16] is performed (i.e., the results are obtained using "rotor B", described above). Figures 1 and 2 show the dimensionless axial and tangential velocities, respectively, at a point which has radial position r/R= 0.8 and is located: (i) on the plane of rotation for the axial velocity; (ii) upstream of the rotor (x3/R= 0.0925) for the tangential velocity. Next, consider a comparison with the experimental results of Guj et al. [14] ("rotor A", defined above). These give the velocity at a point with r/R=0.357 and located upstream of the rotor plane, with x3/R = 0.0572. The comparison is presented in Figs. 3 and 4, which show the dimensionless axial and tangential velocities, respectively, for 2 = 4.42. The agreement appears to be satisfactory. In both cases the results have the same trend, but there appears to be a vertical shift between experimental and numerical results (this is stronger in the case
G. Arsuffi et al./ Windmill rotors in the presence of yaw
158 13
I
12
8
>.
r
o o 0 0 0 3
11
1 0
0
r~ >
09
~
t.l--
x
08
07
06 Numer ical Exp,er i n ~ k a t 0
5
~
6]
28
i
4@
i
6@
,
i
80
8@
hz*mu[h
(deg)
1 12~
i
T 149
[] 0
i
ii
II
T, 168
18@
Fig. 1. C o m p a r i s o n w i t h [15]: a x i a l v e l o c i t y v e r s u s a z i m u t h a t rlR=0.80, x~lR=O; 5 = 0 < (Nw = 60, N s = 12).
o f " r o t o r B"). Specifically, the axial velocity exhibits a small relative shift t hat could be due to the u n c e r t a i n t y of the blockage effect of the wind tunnel. Another source of disagreement between numerical and experimental results is due to the presence of viscous effects which are not considered in the numerical case. Other possible explanations for the disagreement are discussed in Ref. [16] for the case of " r o t o r B". Next, consider the results in the presence of yaw. First, consider a comparison with the experimental results of Wagner and Deppe [15]. The results are obtained using " r o t o r B", described above. Figure 5 presents the power coefficient as a function of the yaw angle. The trend is captured correctly and the agreement between numerical and experimental results is satisfactory. Next, in order to assess the method on the basis of "local" quantities, consider
159
G. Arsuffi et a l . / W i n d m i l l rotors in the presence of y a w 8 15
8
L
OlO
-iJ
G
o >
-~
j
O
)
c @
[
o ~-
Oi
-o
-oo'5
d
o
O!
0
0 0
\
0 0
o
0
! ~EXl::,~r _ ,mental ,~ol .
-0 1S o
.
20
.
.
40
.
.
C:~
.
0 0
1
- T - T -
80
18e
120
140
1C~
18E~
Azimuth (deg)
Fig. 2. Comparison with [15]: tangential velocity versus azimuth at r/R =0.80, x l / R = 0,09; 5=0:" (Nw=60, Ns=12). a c o m p a r i s o n with the e x p e r i m e n t a l results of Guj et al. [14]. These were o b t a i n e d w i t h o u t a c o n t r o l on the a n g u l a r speed (the r o t o r is allowed to assume the " e q u i l i b r i u m " a n g u l a r speed a l t h o u g h a flywheel was used to minimize a n g u l a r v e l o c i t y fluctuations). Therefore, the e x p e r i m e n t a l results exhibit the s i m u l t a n e o u s effects of two p a r a m e t e r s : the yaw angle, 6, and the tip speed ratio, ;.. These effects are p r e s e n t e d in Fig. 6, which shows the dimensionless axial velocity as a f u n c t i o n of the azimuth angle, for t h r e e different cases: (i) 5 = 0 ° and )~=4.42 (as in Fig. 3), (ii) 5 = 0 ° and 2=3.90, and (iii) 5 = 1 5 ° and 2--3.901 . We m a y observe that, for the case u n d e r consideration, the two
' the values of case (iii) are those adopted in the experimental results, shown in Figs. 7 and 8 discussed later.
16ll
(;. Arsuffi et al.: WipTdmill rotors in the preset~ce r~/.v~tc 1 3
J
I I
i
f
1 2
r 8
{ 1
!
I
J b.
1 (3
0 0 ,g >
( 99 1 1i
0 X <~
of
08
I
o 7
@6
I
I
NLI~EP,IC AL E (PEPZPIENTAL
i i i i i
OB 0
--
0 0
i
I L
~ J i i
i
I 2g
J
i
i
i
]
~
i
2-40
18(O
Az, mutb,
300
(de,_-j ]'
Fig. 3. Comparison with [14]: axial velocity versus azimuth, at r/R=0.36, xl/R=0.057; ~5=0 (N~v=60, Ns= 12). p a r a m e t e r s h a v e c o n t r a s t i n g effects on the p e a k values: as ). decreases, the p e a k v a l u e s b e c o m e m o r e p r o n o u n c e d , w h e r e a s , as 6 increases, the u p p e r p e a k d e c r e a s e s (along w i t h m o s t of the curve), while the l o w e r p e a k r e m a i n s a l m o s t u n c h a n g e d . Thus, the c o m b i n e d effect of 6 and )o is not strong. This is confirmed by the e x p e r i m e n t a l results, p r e s e n t e d in Figs. 7 and 8, w h i c h show the dimensionless axial and t a n g e n t i a l v e l o c i t y as a f u n c t i o n of the a z i m u t h a n g l e for ~f= 1 5 and ), = 3.90. The c o m p a r i s o n with the e x p e r i m e n t a l r e s u l t s is satisfactory. The c o m p a r i s o n b e t w e e n Cp d i s t r i b u t i o n s at different a z i m u t h angles is s h o w n in T a b l e 1. V a r i a t i o n s with the a z i m u t h a l a n g l e are in the r a n g e of some percent, In o r d e r to verify t h a t the code is bag-free we c o m p a r e d the results with those o b t a i n e d using an i n d e p e n d e n t l y - d e v e l o p e d f r e q u e n c y - d o m a i n code based on the s a m e f o r m u l a t i o n [13]. T h e c o m p a r i s o n is p r e s e n t e d in Fig. 9, which shows the v a l u e of d) at the blade e l e m e n t 11 =6. L = 5 (i.e., at the
161
G. Arsuffi et a l . / W i n d m i l l rotors in the presence of yaw 06
ml
I__
NUMERICAL EXPERIHBNTAL
0 0
i
@ @4
3
5 0
@
>
@2
( I]
@ @
J
C O I.--
@0
--0
Q
,
2
O
J
J
i
i
L
6@
i
i
i
i
,
120
i
t
3
J
,
L
l
,
180 AzimuLh
,
i
,
240
i
,
,
J
i
J
J
i
,
31~
(deg)
Fig. 4. Comparison with [14]: tangential velocity versus azimuth at r/R = 0.36, x l / R = 0.057; 5 = 0 (Nw=60, Ns=12).
trailing-edge e l e m e n t w h e r e the section lift coefficient is higher), for 3 = 15 ° and )~= 4.5: not surprisingly, the r e s u l t s are in e x c e l l e n t a g r e e m e n t w i t h e a c h o t h e r (the difference is due to the zero initial c o n d i t i o n s for the t i m e - d o m a i n solution). Next, in o r d e r to verify t h a t the d i s c r e p a n c y b e t w e e n t h e o r e t i c a l and e x p e r i m e n t a l r e s u l t s is not due to the a s s u m p t i o n on the wake, we p e r f o r m e d a f r e e - w a k e analysis. In o r d e r to assess b e t t e r the i m p a c t of the f r e e - w a k e analysis, we c o n s i d e r e d a case t h a t is p a r t i c u l a r l y s e n s i t i v e to the w a k e g e o m e t r y , i.e., we used a h i g h e r v a l u e for the tip speed r a t i o ()~= 14 w i t h 0 = 2°). The results for the w a k e g e o m e t r y are p r e s e n t e d in Fig. 10 w h e r e the sections b e t w e e n a v e r t i c a l p l a n e and the w a k e surface are s h o w n for the first two spirals. T h e cases w i t h 3 = 0 ~ and 3 = 15" are c o m p a r e d . T h o s e for the p r e s s u r e coefficient as a f u n c t i o n of xl/c (for r/R=0.8 and a z i m u t h equal to zero) are p r e s e n t e d in T a b l e 2, for p r e s c r i b e d and free wake, respectively. T h e s e r e s u l t s
162
(;. A r s u / ] i et (ll. ; W i n d m i l l
r o t o r s lt? t h e p r e s e n c e o f y a w
06
I -iJ
04
c ¢
o i ¢
o c) L
¢
o 02
,col
Expe- i ,'t,~--,t:a I
I
09
0
10
29 Yaw a n g l e
0
t 3@
(deg)
Fig. 5. C o m p a r i s o n w i t h [15]: p o w e r c o e f f i c i e n t v e r s u s y a w a n g l e (Nw = 60, Ns = 12).
indicate that e v e n for this case, w h i c h is d e s i g n e d to e m p h a s i z e the free-wake effect, the differences are r e l a t i v e l y small. Therefore, we believe that the n u m e r i c a l results w i t h a prescribed-wake model are an a c c u r a t e r e p r e s e n t a t i o n of the potential-flow s o l u t i o n . As ment i o n e d a b o v e the d i s c r e p a n c y b e t w e e n n u m e r i c a l and e x p e r i m e n t a l results m a y be due to (1) the e x p e r i m e n t a l error (estimated to be of the order of 4% [14]), to (2) the effects of v i s c o s i t y (not i n c l u d e d in the n u m e r i c a l formulation), and to (3) the t r u n c a t i o n error (see, h o w e v e r , the c o n v e r g e n c e c o n s i d e r a t i o n ment i o n e d above). Finally, in order to s h o w the g e n e r a l i t y of the proposed method, c o n s i d e r a case w i t h n o n c o n s t a n t y a w (i.e., w i t h [2 # 0). R o t o r A w a s used w i t h 2 = 14 and 0 = 2 . We a s s u m e d 3 = 0 for t < to (where to c o r r e s p o n d s to 15 r e v o l u t i o n s , in order to r e a c h steady state), ~ = 5 1 = 1 5 ° for t > t l (where tl corresponds to 16
G. Arsuffi et al./Windmill rotors in the presence of yaw I 3
l /i /~
163
] w
$ = 0 °, ~ = 4.42 6 : 0 °, ,~ : 3.90
1 2
6 = 15°~A = 3 . 9 0
.i
8
°
/ 2"/
2
I 0
>
09
-
"~
/"
/
-o x
<
08
I -~'~'~ I"
/
/
~/,/
;/ 07
0 5
J I I L L
0
, L J , I
6Q
i
~
,
12~3
,
i
,
J
~
1~
,
,
,
240
,
J
,
L
,
J
,
i
i
3~D
h z imuLh ( d e g )
Fig. 6. Axial v e l o c i t y versus azimuth at r/R=0.36 and x l / R = 0 . 0 5 7 for different values of 6 and ). (Nw =60, Ns=12).
r e v o l u t i o n s ) , and
6=61sin2
tt-----to)'
for
to
(12)
Three extra r e v o l u t i o n s after tl w e r e i n c l u d e d in the c a l c u l a t i o n s in order to better identify the t r a n s i e n t r e g i o n in the s o l u t i o n . Figure 11 s h o w s ¢ at a g i v e n panel v e r s u s a z i m u t h angle. Figure 12 s h o w s the Cp distributions at t = to, t = to plus o n e r e v o l u t i o n , and after t w e n t y r e v o l u t i o n s w h i c h are suffic i e n t to r e a c h the periodic solution.
(L A r s u ~ et a l . / W i n d m i l l rotors in the presence of y a w
1(;4 13
12
@
I 1
>
I @
0 >
09
× "~
08
@7
@6 MJ'~_RICAL E)~ERI~AL
--
[3
0
85 8
60
1;_>0
! 80
240
380
36@
Az, rau tF, ( d e g )
F i g . 7. Comparison with [14] : axial velocity versus azimuth at r/R = 0.36, x ~/R = 0,057 ; d = 1 5 ( N w = 6 0 , N s = 12).
5. C o n c l u d i n g r e m a r k s A method for prescribed- and free-wake unsteady aerodynamic analysis of rotors in the presence of yaw has been presented. The proposed method is very general and it seems suited for a wide range of unsteady aerodynamics applications. The numerical results obtained for HAWT rotors are in satisfactory agreement with existing experimental results. The free-wake results indicate that wake distortion due to yaw should be taken into account only for high yaw angles: in this case, however, the potential-aerodynamic assumption is questionable due to the presence of significant viscous-flow effects. Because of the importance of the evaluation of blade loads in such high yaw angle conditions, we recommend that future work be directed toward the development of a viscous-flow capability.
165
G. Arsuffi et a l . / W i n d m i l l rotors in the presence of yaw
06
__.1_
l__
~[~ME~IC~_ £XP£R]]'ENTAL
r] 0
r
8 04 ).
(t
3 0
C
>O
-8
Ilnl O, _
02
-'t
0
o
@
0
c o
(
O0
--~
2
i
9
J
L
J
i
i
68
i
i
*
i
i
128
i
i
i
i
I
i
18e
i
,
*
248
,
,
i
,
,
300
,
,
,
J
36~
Az,muLh (deg) Fig. 8. Comparison with [14] : tangential velocity versus azimuth at r/R = 0.36, x~/R = 0.057; = 15° (Nw = 60, Ns = 12).
Acknowledgements This work was partially supported by the Italian Commission for Nuclear a n d A l t e r n a t i v e E n e r g y S o u r c e s ( E N E A ) u n d e r c o n t r a c t No. 13984 t o t h e University of Rome "La Sapienza".
Appendix. Integral formulation for exterior Neumann problem for Laplace equation In this appendix we present an outline of the boundary element method for i n c o m p r e s s i b l e p o t e n t i a l flows. 2 F o r t h e s a k e o f s i m p l i c i t y , w e a l w a y s a s s u m e 2 The presentation of this section follows closely that of Morino and Tseng [12], to which the reader is referred for details.
166
(;. Arsu/fi et el.~ Windmill rotors in the presence o[y¢~l~
Table 1 C~, distribution versus azimuth position for suction (pressure) side
x/c
Reference value
Relative difference (%)
azimuth = 0
90
180
270
0.49 (0.45)
- 0.881 (0.0173)
- 0.3 (-58)
0.8 (125)
0.5 (.-62)
-. 0.39 (-0.29)
0.534 (0.268)
1.5 (3.7)
3.0 (7.5)
- 0.9 (--3.7)
-0.18 (-0.06)
0.401 (0.233)
- 1.7 (-3.0)
3.7 (-6.0)
- 1.7 (3.0)
0.06 (0.18)
0.246 (0.122)
- 2.4 ( - 2.5)
- 6.5 ( - 4.9)
- 2.4 (-2.5)
0.29 (0.38)
0.0712 ( - 0.0384)
4.1 ( - 3.4)
6.2 ( - 6.0)
2.1 (-2.6)
0.45 (0.49)
-0.146 (--0.293)
0.2 (-5.1)
1.4 (-7.5)
1.6 (3.1)
that the fluid and the body are initially at rest. In addition, the frame of r e f e r e n c e is a s s u m e d to b e c o n n e c t e d w i t h t h e a i r , u n l e s s o t h e r w i s e s t a t e d . T h e g o v e r n i n g e q u a t i o n s f o r a n i n c o m p r e s s i b l e i n v i s c i d flow a r e t h e c o n t i n u ity equation and the Euler equations. From these equations one obtains K e l v i n ' s t h e o r e m w h i c h s t a t e s t h a t , t h e c i r c u l a t i o n o v e r a m a t e r i a l c o n t o u r is c o n s t a n t in t i m e . A s m e n t i o n e d a b o v e , t h e flow field is i n i t i a l l y i r r o t a t i o n a l . T h e n , a t t = 0 w e h a v e F = 0 for a l l t h e c o n t o u r s in t h e field, a n d K e l v i n ' s t h e o r e m y i e l d s t h a t I ' = 0 a t a l l t i m e s for a l l t h e m a t e r i a l c o n t o u r s t h a t rem a i n e d i n t h e field b e t w e e n t = 0 a n d t h e c u r r e n t t i m e . T h e r e f o r e , s u c h a field remains irrotational at all times, except for the points that come in contact with the boundary surface, as for these points there do not exist a contour that remains in the fluid region at all times. Indeed, they form a surface called the wake. C o n s i d e r f i r s t t h e c a s e w i t h o u t a w a k e . S i n c e t h e flow is i r r o t a t i o n a l , t h e r e e x i s t s a f u n c t i o n ~b ( v e l o c i t y p o t e n t i a l ) s u c h t h a t
v=V~,
(1)
Combining the continuity equation, V'v=0, V2~b=0
(x o u t s i d e S , ) .
w i t h Eq. (1) o n e o b t a i n s (2)
N e x t , c o n s i d e r t h e b o u n d a r y c o n d i t i o n s . A s t y p i c a l , t h e s u r f a c e o f t h e b o d y , SB, is h e r e a s s u m e d to b e i m p e r m e a b l e ; a c c o r d i n g l y , t h e b o u n d a r y c o n d i t i o n o n SB is (v v , ) ' n = 0 , w h e r e vs is t h e v e l o c i t y o f t h e p o i n t o n SR. C o m b i n i n g w i t h
G. Arsuffi et al./Windmill rotors in the presence of yaw
167
0 010
0 008 0 000 0 004
S ~
0 002
¢ o n
0
-0 000
-0 002
> -0 004 -0 0OG -0 008 - 0 010 8
30{3
600
900
~28e
~S ~
1800
Az irnuEh ( d e g )
Fig. 9. Time and frequency domain solutions [13]: ¢30 (r/R=0.76) versus azimuth ("Rotor A", Nw = 60, Ns = 12). Eq. (1), one o b t a i n s
On = vB' n .
(3)
In addition, in the air f r a m e of r e f e r e n c e used here, the b o u n d a r y c o n d i t i o n at infinity is g i v e n by v = O( IIx II-2), or, in t e r m s of the v e l o c i t y p o t e n t i a l , = O< IIx If - ~).
¢4)
T h e p r o b l e m described by Eqs. (2), (3), a n d (4) is k n o w n as the e x t e r i o r N e u m a n n p r o b l e m for the L a p l a c e equation. Once the a b o v e p r o b l e m h a s b e e n solved, the p r e s s u r e m a y be e v a l u a t e d u s i n g B e r n o u l l i ' s t h e o r e m for p o t e n t i a l i n c o m p r e s s i b l e flows (a first i n t e g r a l of
16~
(;. A r s u f i et ~1.' W i n d m i l l rotors it; the presence of y¢tw
I
8.16
~=0
'
,
'
l
°
5=15 ° 0. t4
f
"
I
0.12
i I
r~
/
t
E
F
k9
/
L,.T,
g ,a.lo
"o o re
i
I
I
0.06
! II
0.04
x ~
0.02 -0.~
I II k I
.
0.00
0.02
x.
.
0.04 Axial
.
.
0.06 po~i~.ion
.
0.88
O. 10
0.12
(m)
Fig. 10. F r e e - w a k e a n a l y s i s : w a k e g e o m e t r y (first t w o spirals) for 6 = 0 A " . N w = 8 4 , N s = 1 2 , Nf.=2).
a n d 6 = 15 ( " R o t o r
Euler's equations) which, in the air frame, is given by 5q5 ~ , 1 1 __ p= p~, 5t + ~ v - + p P
(5)
where co indicates evaluation at infinity. Next, consider the corresponding boundary integral equation formulation, the basis for which is Green's second identity , (fV2g-gV2f)d't'=
--
,/
kl~n-g~n] d,9°,
(6)
where n is the inwardly directed normal. The essence of the method consists in identifying f with ~b, and in choosing for g the so-called fundamental solution of
G. Arsuffi et al./Windmill rotors in the presence of yaw
169
Table 2 C 0 distribution:relative difference (%) b e t w e e n p r e s c r i b e d a n d f r e e - w a k e r e s u l t s
x/c suction side x/c p r e s s u r e side
0.49 -0.6 0.45 24.
- 0.39 0.7 -0.29 1.5
- 0.18 1.0 -0.06 2.1
0.06 1.2 0.18 4.9
0.29 2.8 0.38 - 20.
0.45 0.4 0.49 - 7.2
-8 04
- 0 86
-8 88
¢
-0
18
-8
12
-8
14
-8
16
I
o & >.
0 >
~0
~1
-@ 18 ' , I , , I . . 5@48 5228 54@6 5588 5760 5948 Az,mu(½
. . 6120
. 638@
6480
666@
6848
(deg)
Fig. 11. F r e e - w a k e a n a l y s i s : ~b~o v e r s u s a z i m u t h (r/R=0.76) ( " R o t o r A", N w = 8 4 , N s = 1 2 ,
NF=2).
the Laplace equation, G, which is defined by V 2 G = 5 ( x - x . ) (where 5 is the Dirac delta function), with boundary condition G= O( ]1x II- 1) at infinity. This is the well k n o w n potential for the unit source, G= -1/4~ ]1x - x . ]j. For exterior problems, we apply Eq. (6) to a volume ~/- ouside ,9~s, bound by a spherical surface, ~ R , having radius R and center x . . Taking the limit as
17()
(;. Arsuffi et el.,, W i n d m i l l rotors in the presence of y a w 0
)
Q6
02 n © I
-02
-06
t=tO t:tg * t: Lint-
i-I
I
rov
0 A
-70
-OS
00
OS
CHORD
Fig. 12. Free-wake analysis: cp distribution at r / R = 0.80 for different time steps ("Rotor A", Ns=12, N~ = 2).
Nw=84,
R tends to infinity, using the boundary condition at infinity for ~b(Eq. (4)), and the corresponding one for G, one obtains, using Eq. (2), ~b(x,)=
8~n a - ¢
~nn d,~,
(7)
where n is the outwardly directed normal to ,~B. Equation (7) may be used to obtain the values of ¢ at any point in the field if the values of ¢ and ~ ¢ / 8 n on the boundary surface &~. are known. We will refer to Eq. (7) as the boundary integral representation for the velocity potential. Note that, if x , approaches a point on the surface ,90,, the value of ~b(x,) approaches the value of ~b at that point of 9~B. This yields a compatibility condition, between the values of ¢ and 8¢len on 5%, which is satisfied by any
G. Arsuffi et al./Windmill rotors in the presence of yaw
171
function that satisfies the Laplace equation. In our case ~c~/Snon 5zB is known from Eq. (3) (Neumann boundary condition), and hence this compatibility condition yields a boundary integral equation that may be used to evaluate on 5%. This integral equation and its extensions are the key to the methodology reviewed here (as well as to BIE methods of direct type in general). Once, ¢ on 5fB is known, we may evaluate ¢ anywhere in :t~, using Eq. (7), and hence the pressure p, using Bernoulli's theorem, Eq. (5). Next, we recall that, according to d'Alembert's paradox, if the flow around an object is irrotational everywhere, then the resultant of the forces acting on the body equals zero. For this reason, the case without wake has little or no interest in aeronautics. In this following, we present a potential-flow formulation for flows around bodies capable of producing lift, i.e., a potential-flow formulation with vorticity present in the field. This may be accomplished if we allow for the potential to be discontinuous in the field. Thus, in order to complete the problem, we need boundary conditions for the wake. These are obtained from the principles of conservation of mass and momentum across a surface of discontinuity, and use (V-Vw)'n=O (no penetration) and Ap = 0 (no pressure discontinuity). For potential flows, the first yields
A ~
=o,
(s)
On the other hand, combining the condition Ap = 0 with Bernoulli's theorem, Eq. (5), yields
~(~)1 2 S t"~-2Vl-1 ~t2--1V22=A ( ~ t ) -~-½(vl'~-D2)°(D1-D2)=0'
(9)
or Dw A¢ = O, Dt
(10)
with Dw/Dt=~/~t+ Vw" V, where Vw=½(Vl + v2). Equation (10) implies that the value of A~b remains constant following a wake point, Xw, which is defined as a point having velocity Vw, and equal to the value it had when Xw left the trailing edge. Numerically, this value is obtained from the trailing-edge condition that, at the trailing edge, A¢ on the wake is equal to the potential discontinuity on the body. Next, consider the integral formulation. Let £fBw denote a surface that surrounds the body, ~B, and a thin layer that includes 5fw. Note that Eq. (2) is is valid outside this surface. Thus, Eq. (7) is still valid, as long as 5~Bis replaced with 5fBw. Next, let the two sides of the portion of ~Bw that surrounds 5fw
(;. Arsuffi et al./Windmill rotors in the presence of yaw
172
become infinitesimally close t o / / ' w . In the limit, recalling Eq. (8), one obtains
~'/l~
~nG-~bOn]
d>/-
'tv.
Adp--d'(/'"
(1l)
~,i2
This is the b o u n d a r y integral r e p r e s e n t a t i o n for the potential for flows a r o u n d lifting objects. As before, in the limit as x , tends to SR, the above representation yields a compatibility c o n d i t i o n between ~b and 5¢/5n on ,(/~ and A~b on Y~w. Note t h a t 5~/5n is k n o w n from the impermeability b o u n d a r y c o n d i t i o n on ,~h (Eq. (3)), w h e r e a s A¢ and the wake g e o m e t r y are k n o w n from the preceding time history. Thus, the compatibility c o n d i t i o n yields a b o u n d a r y integral e q u a t i o n for the u n k n o w n ~ on ,~/~. Then, Eq. (11) yields ~, and from this the velocity, v = V ¢ , as well as the pressure (from B e r n o u l l i ' s theorem, Eq. (5)) a n y w h e r e in the field. In addition, from the velocity of the wake points one o b t a i n s the wake g e o m e t r y at t + d t (with A¢ c o n s t a n t following the wake points and equal to the value it had when Xw left the trailing edge) and the process m a y be repeated.
References [1] B. Montgomerie and A. Zdunek, Dynamic response in horizontal axis wind turbines including instationary effects in the stream tube, Proc. European Wind Energy Conference, Hamburg, FRG, 1984. [2] S.E. Mattsson, Modelling and control of large horizontal axis wind power plants, Doctoral Thesis, Dept. of Automatic Control, Lund Institute of Technology, Sweden, 1984. [3[ A. Hemon, A. Zervos and S. Huberson, Numerical computation of unsteady forces on HAWT, Proc. European Wind Energy Conference, Glasgow, UK, 1989. [4] L. Morino, A general theory of unsteady compressible potential aerodynamics, NASA CR-2464. [5] B. Montgomerie, Performance methods (again!), Proc. 2nd Symposium on Aerodynamics of Wind Turbines, IEA R&D WECS Joint Action, Lyngby, DK, 1988. [6] J.C. Gohard, Free wake analysis of wind turbine aerodynamics, ASRL TR-184-14, 1978. [7] P.H. Alfredsson, J.A. Dahlberg and F.A. Bark, Some properties of the wake behind horizontal axis wind turbines, Proc. 3rd International Symposium on Wind Energy Systems, Lyngby, Denmark, pp. 469 84, Aug. 26-29, 1980. [8] M.B. Anderson, D.J. Milborrow and J.N. Ross, Performance and wake measurements on a 3 m horizontal axis wind turbine Comparison of theory, wind tunnel and field test data, Proc. 4th International Symposium on Wind Energy Systems, Stockholm, Sweden, pp. 113 135, Sep. 21-24, 1982. [9] B.R. Clayton and P. Filby, Measured effects of oblique flows and change in pitch angle on performance and wake development of model wind turbine, Proc. 4th BWEA Wind Energy Workshop, Cranfield, Bedford, United Kingdom, Mar. 24 26, 1982. [10] G.T. Atkinson, S.M.B. Wilmshurst and D.M. Wilson, The aerodynamics of rotors in yaw, Proc. 9th BWEA Conference, Edinbourgh, United Kingdom, pp. 71 78, Apr. t 3, 1987. [11] L. Morino and K. Tseng, A general theory of unsteady compressible potential flows with application to airplanes and rotors, Boundary Element Method in Nonlinear Fluid Dynamics, Ch. 6, eds. P.K. Banerjee and L. Morino (Elsevier Applied Science., London, 1990).
G. Arsuffi et al./Windmill rotors in the presence of yaw
173
[12] L. Morino, L.T. Chen and E.O. Suciu, Steady and oscillatory subsonic and supersonic aerodynamics around complex configurations, AIAA J., 13(3) (1975) 368-374. [13] F. Mastroddi and L. Morino, Time and frequency domain aerodynamics for flutter of helicopter rotors in hover, International Symposium on Boundary Element Method, Hartford, CT, USA, 1989. [14] G. Guj, M. Terzitta and G. Arsuffi, Velocity measurements upstream of a windmill rotor model, Wind Eng., 15(5) (1991) 248 260. [15] W.J. Wagner, L. Deppe, Messung des Momentanen Geschwindigkeitsfeldes einer Modellwindturbine sowie der Auftriebsverteilungen der rotierenden B1/itter, DFVLRIB 232-88 J 06, DLR, Institut fur Aeroelastik, GSttingen, FRG. [16] W.J. Wagner, private communication (1990).