Viscous Flow and Dynamic Stall Effects on Vertical-Axis Wind Turbines

The present paper describes a numerical method, aimed to simulate the flow field of vertical-axis wind turbines, based on the solution of the steady, incompressible, laminar Navier-Stokes equations in cylindrical coordinates. The flow equations, written in conservation law form, are discretized using a control volume approach on a staggered grid. The effect of the spinning blades is simulated by distributing a time-averaged source terms in the ring of control volumes that lie in the path of turbine blades. The numerical procedure used here, based on the control volume approach, is the widely known 
“SIMPLER” algorithm. The resulting algebraic equations are solved by the TriDiagonal Matrix Algorithm (TDMA) in the r- and z-directions and the Cyclic TDMA in the 0-direction. The indicial model is used to simulate the effect of dynamic stall at low tip-speed ratio values. The viscous model, developed here, is used to predict aerodynamic loads and performance for the Sandia 17-m wind turbine. Predictions of the viscous model are compared with both experimental data and results from the CARDAAV aerodynamic code based on the Double-Multiple Streamtube Model. According to the experimental 
results, the analysis of local and global performance predictions by the 3D viscous model including dynamic stall effects shows a good improvement with respect to previous 2D models.

ind energy is considered to be a promising renew- able energy source, particularly for remote areas, like islands and mountainous regions.Advances in largescale wind energy systems, significant cost reductions, more efficient manufacturing techniques, increasing technical experience and the increasing environmental awareness have contributed greatly to the use of wind energy.Although invented in 1931, the Darrieus turbine did not see extensive development until the 1970s during the world wide oil crisis where a considerable amount of work has been done with regard to the aerodynamics of vertical-axis wind turbines.Aerodynamic analysis of vertical-axis wind turbines is based on several methods which can be classified into two major categories: momentum models and vortex models.Updated reviews of the state of the art of such methods, including the appropriate references are provided by Strickland [1986], Turyan et al. [1987], and Paraschivoiu [1988].
The major contribution to the unsteady aerodynamics of the Darrieus rotor is dynamic stall, which occurs at low tip-speed ratios and its effects have a significant role in drive-train calculation, generator sizing and overall system design.Therefore, the ability to predict dynamic stall is of crucial importance for optimizing the Darrieus wind turbine.The theoretical methods of dynamic stall are still limited in their scope and prohibitively expen- sive in term of computer CPU time.Although Navier- Stokes solvers have been used to simulate dynamic stall (Shida et al. [1987], Daube et al. [1989], Tuncer et al. [1990], and Visbal [1990]), most of the studies were concerned with helicopter retreating-blade stall or agile fighter aircraft and have thus considered pitching airfoils and not the rotation motion encountered by the blades of a Darrieus turbine.That is why our research group, attached to the J.-A.Bombardier Aeronautical Chair, is working on a Navier-Stokes solver able to simulate the flow around an airfoil in Darrieus motion under dynamic stall conditions (Tclaon et al. [1993]).Our research group has also investigated several semi-empirical methods such as Gormont, MIT and Indicial models (Allet and  Paraschivoiu [1988]).For practical reasons, the indicial model, which is one of the most recent semi-empirical models, is used in this study to simulate the effects of dynamic stall.

THREE-DIMENSIONAL VISCOUS MODEL
The steady, incompressible, laminar Navier-Stokes equations in cylindrical coordinates are solved for the flow field and Performance characteristics of a vertical-axis wind turbine.This approach was first developed for the Euler equations (Rajagopalan [1985]) and extended to Navier-Stokes equations (Rajagopalan [1986]).The flow equations, written in conservation law form, are dis- cretized using a control volume approach and the result- ing elliptic equations are solved by a line-by-line method.The numerical procedure used for solving the flow governing equations is based on the widely known "SIMPLER" algorithm developed by Patankar [1980].The turbine blades are modeled through the source terms in the momentum equations used to solve the flow field.The essentials of this analysis are presented in the following subsections.
Flow Field Modeling Some general features of the Sandia 17-m vertical-axis wind turbine are presented in Fig. 1.The top view of the rotor according to section (A-A) shown by Fig. 2 give: 0: or) (1) where V (ucos + wsin) V (v-or) (2) Since or, the linear velocity of the turbine blade_z is known, the only unknown in the above equation is Vabs.
Solution of the velocity field would yield the knowledge of all the quantities of interest.
With the assumptions stated previously, the governing equations for the flow field of a vertical-axis wind turbine are: In the simulation of the incompressible flow, the basic equations are generally expressed as the conservation of physical values such as mass or momentum.According to the Finite Volume Method (FVM), the governing equations are integrated over each control volume cell to derive the discretized equations from them.For velocity- pressure coupled flows, a staggered grid system is known to give more realistic solutions and is adopted in the present analysis.

Computational Grid
The computational domain is subdivided into control volumes by series of grid lines orthogonal to the f'r,0, z coordinate directions.The size of the control volumes is decided based on the accuracy required in the region.Patankar 1980] suggested to use a staggered grid for the storage of each variable to avoid a wavy solution and this strategy is adopted in the present study.In the staggered grid, the velocity components are calculated for the points that lie on the faces of the control volumes.Thus, the r-direction velocity u is calculated at the faces that are normal to the r-direction.It is easy to see how the locations for the velocity components v and w are to be defined.Typical cells for each variable (u, v, w, p) are shown on Figures 3 and 4.

Turbine Modeling
The motion of vertical-axis wind turbine blades intro- duces primarily a change in the local momentum of the fluid due to the forces generated by the rotating blades.Thus, an obvious place to introduce the action of the blades in the governing equations of the flow field is through the source terms in the momentum equations, namely Sr, So and Sz.These source terms are valid only for the computational cells that lie in the path of the turbine blades.Following the detailed procedure illus- trated by Allet [1993], only the final form of the time-averaged source terms is given here: S WCOS(CDUCOS CLV + Czwsin) (7) S W(CDV + CLV CDtOr  VISCOUS FLOW AND DYNAMIC STALL These are forces felt by the blades per unit of span, and are thus subtracted in the discretized momentum equations to effect an equal but opposite reaction on the fluid at a specific cell.

DYNAMIC STALL ANALYSIS
Dynamic stall is an unsteady flow phenomenon which refers to the stalling behavior of an airfoil when the angle of attack is changing with time.It is characterized by dynamic delay of stall to angles significantly beyond the static stall angle and by massive recirculating regions moving downstream over the airfoil surface.In the case of the Darrieus turbine, when the operational wind speed approaches its maximum, all blades sections exceed the static stall angle, the angle of attack changes rapidly, and the whole blade works in dynamic stall conditions.This phenomenon may cause structural fatigue, and even stall flutter leading to catastrophic failure, and thus can become, in many cases, the primary limiting factor in the performance and structure of the turbine.Therefore, in order to predict the aerodynamic performances of such structures and provide aerodynamic loading information to structural dynamic codes, a numerical model must be able to capture the complex aerodynamics encountered by a vertical-axis rotor blade, and more particularly the dynamic stall phenomenon.
To provide some representation of dynamic stall effects, there are several methods which can be classified in two categories: the theoretical and semi-empirical methods.Theoretical methods are attractive and accurate in studying local insight of the dynamic stall phenomenon (boundary layer evolution, pressure distribution) even if the required computational resources are prohibitive.However for practical cases, semi-empirical meth- ods are still used to simulate dynamic stall effects within engineering accuracy.A general review of some dynamic stall prediction models is provided by Reddy and Kaza  [1987].

Indicial Model
The methodology of semi-empirical models often relies on the reduction and resynthesis of aerodynamic test data from unsteady airfoil tests.In the interest of computa- tional simplicity, some models sacrifice physical realism and so may have limited generality in their application.
The aim of this section is to introduce the indicial method used in our aerodynamic analysis to simulate dynamic stall.This model, originally developed by Beddoes 1983, 1984] in the early 80' s, has continued to improve since then (Beddoes and Leishman [1989]).The theoretical approach of the indicial method is quite different from the other methods.The indicial method is based on the fact that dynamic stall is a superposition of distinct effects that may be studied separately by using indicial function (response of a system to a disturbance).These effects consists of three distinct parts: Potential flow solution (linear solution); Separated flow solution for the non-linear airloads; Deep stall solution for vortex induced airloads.This model is then represented by a relatively simple set of equations and, most importantly, allows closer identification of the interacting phenomena involved in dynamic stall.Potential flow calculation.Any unsteady aerody- namic model must be able to represent correctly the attached flow behavior.This flow induces two aerodynamic components.The first component is a circulatory one (CNc) due to the change in the angle of attack and has a lag behavior.Cuc(t CIaE coso e (10) where e E, the effective angle of attack is defined by: with OE--" On_ + 0 (t) (11) c(t) Ao-A exp(-blS' -A 2 exp (-b2s' (12)   The constants of Eq. ( 12) are determined on the basis of theoretical considerations and optimization to fit experi- mental results (Beddoes [1984]), so the constants are given as A o 1.0, A 0.3, A 2 0.7, bl 0.35 and b 2 0.68.The second aerodynamic response is an impulsive component resulting from an induced velocity normal to the airfoil surface generated by the airfoil movement and is given by: where (4(t)AoO CNI( The normal force coefficient CN resulting from potential flow is obtained by a simple summation:

CN(t CNc(t + CNI(t (15)
The tangential force coefficient is given by: CT(t CLa12E sin12 e (16) Non-linear effects.The unsteady non-linear effects are simulated through empirical correlations.Leading edge and trailing edge separations must be predicted with accuracy because of their great influence on the aerodynamic coefficients.To model trailing edge separation, the position of the separation point f is required.This position will be corrected later for unsteady cases.f= 0.38 exp(12 12)/S 12<-12 (17) f= 0.62 + 0.38(12 12)S (exp(-1/S)-1) 121 < 12 <--t22 (18) f= 0.025 + 0.175 exp(122 c)/S 2 12 >c 2 (19) A typical curve of the separation point f vs. the angle of attack 12 is shown in Fig. 5.The parameters 121, o2, S, $2, $3 are found from the curve-fitting of the static lift and drag curves.The non-linear value of the normal force coefficient is determined by" CNN-L + f0.5)2 CuC (20) For unsteady flows, temporal effects influence the posi- tion of the separation point.Two first order lags are therefore introduced to take these effects into account.
The first function simulates the airfoil unsteady pressure response.
C N, +p(t) CNN_ L (21) where The onset of leading edge separation is related to the following criterion: CN, > CN1 with CN1 being the critical normal force coefficient illustrated in Fig. 6.A second linear lag function is introduced to simulate additional effects of the unsteady boundary layer response: where ( ,--2Vt A new angle of attack, 12F is defined from the value of CNF," CNF' (25)   This angle is used instead of ot in Eqs. 17-19to compute the position of the unsteady separation.Deep dynamic stall.For most dynamic stall cases, a vortex appears near the leading edge and moves down- stream on the upper surface of the airfoil.Until the detachment of the vortex, it seems there is no significant change in the airfoil pressure distribution.Thus, the induced vortex lift contribution is equal to the loss of lift due to the separation (Fig. 7) From experimental observation, it has been proven that the reattachment of the unsteady flow occurs later relatively to the steady case.This delay is introduced by angle of attack FIGURE 7 Dynamic stall vortex contribution.
replacing o from eqs. 17-19 by a modified angle o given by: (28)

RESULTS AND DISCUSSION
This section presents a selection of results obtained by both the three-dimensional viscous model developed previously and the CARDAAV code based on streamtube theory.The corresponding results are compared with the experimental data on the Sandia 17-m wind turbine configuration provided by Akins [1989].The dynamic stall is simulated by the indicial model for both aerodynamic codes.First, to validate the viscous code, preliminary tests dealing with the determination of the computational grid and the convergence of the numerical model are presented.

Considerations
To validate this solver, preliminary results for validating the computational procedure are presented.Thus, for the wind turbine problem, it was necessary to first establish the size of the computational domain which is referred to as the world size.The domain was considered to be large or accurate enough when the power coefficient reached an asymptotic value.This exercise was carried out on the Sandia 17-m wind turbine rotating at 42.2rpm and a tip-speed ratio Xeq 6.12.This value was used instead of other experimental values to avoid the effects of dynamic stall phenomena on the computed results.
As shown in Fig. 8, the power coefficient Cp reached its asymptotic value at a world size of about 60Req.The same exercise was carried out for the 0-grid size which is illustrated by Fig. 9.According to this figure, a minimum of 40 points are needed to get a sufficient computational grid in the 0-direction.For the z-direction, Fig. 10 shows that 26 points at least were necessary to get an asymptotic value of Cp.Thus, for all the following computa- tions, the fineness of the grid was maintained at (66 48 26) points.A last feature about Fig. 10 is that at the end of the optimization exercise, the calculated value of the power coefficient Cp is almost equal to the corre- sponding experimental value (Cp(exp) 0.387).A conver- gence history of Cp with the number of iterations is presented in Fig. 11 for the same conditions used in the determination of the grid size.As seen in this figure, a converged solution is obtained as early as 100 iterations World size in diameters Effect of world size in r-direction on power coefficient.
and even before that.In general, the calculations were taken to 250 iterations to assure the convergence behav- ior of the computed solution for different operational regimes of the Darrieus wind turbine.Following this optimization exercise, the resulting computational mesh (66 48 26) used by the viscous solver is shown by Fig. 12.

Local Aerodynamic Coefficients
The distribution of local normal force coefficient C N at the equator level is plotted as a function of the azimuth angle to for values of tip-speed ratio Xeq varying from 2.33 to 4.60.for the Sandia 17-m wind turbine rotating at 38.7rpm according to the experimental data (see Figs. 13-15).Figure 13 shows that the viscous model represents quite accurately the distribution of the normal force coefficient for almost all azimuth angle values.However, near 0 -60 deg. in the downwind zone, predictions are different from experimental data.Akins 1989], in his experimental data analysis, has compared the experimen- tal data to values predicted by a vortex model coupled with Gormont model (as dynamic stall model) as modi- fied by Mass6 [1981]  that the wake interaction in the downwind zone doesn't allow the establishment of dynamic stall.However this effect is not predicted by the vortex model (used by Akins in his analysis) nor by the viscous model.In general, values predicted by CARDAAV code are in good agreement with experimental data even if the maximum value in the upwind part is slightly overesti- mated.With increasing tip-speed ratio (Xeq 3.09 and 4.60), this feature (wake interaction) tends to disappear and the values predicted by both models agree quite well with the corresponding experimental data.This first comparison of the normal force coefficient tends to demonstrate a certain superiority of the indicial model as a dynamic stall model when compared to other semi-empirical models (Gormont model for instance).Thus, CARDAAV code coupled with Gormont model predicted a maximum value of CN in the upwind part much greater than the corresponding experimental value CNmax 1.6 (Allet and Paraschivoiu [1988]).To better visualize the dynamic stall phenomena characterized by an hysteresis loop on normal force coefficient Cy curves vs. angle of attack o, Fig. 16   normal force coefficient vs the angle of attack for a tip-speed ratio Xeq 2.49 predicted by both models.This figure shows that for both models dynamic stall does exist in the upwind part (hysteresis loop) and that the reattachment of the boundary layer does happen at the same time for both models.
The values of local tangential force coefficient Cw predicted by both aerodynamic codes are presented as a function of the azimuth angle in Figs.(17)(18)(19).Generally speaking, the viscous code reproduces quite well the experimental values of the tangential force coefficient Cw except near 0 30 deg. in the upwind zone where the predicted values are slightly overestimated.The onset of dynamic stall tends to give a kind of irregularity to the curves of the tangential force coefficient which tends to disappear with increasing tip-speed ratio (Figs. 18and   19).CARDAAV predictions of Caare in general in good agreement with experimental data, though the part de- fined by 150 deg.< 0 deg. in the upwind zone presents predictions completely different from the corresponding experimental data.Moreover, for almost all cases pre- sented, this model overestimates the maximum value of the tangential force coefficient in the upwind and down- wind zones of the Darrieus wind turbine.

Global Performance
The evaluation of the power coefficient Cp was computed on the Sandia 17-m wind turbine operating at 42.2 and 50.6rpm.Even if the effect of dynamic stall is not really visible on power coefficient curves, however this kind of curves allows one to have a better idea about the accuracy of the results at low wind velocities while the  power curves deal with greater wind velocities.Predicted and experimental values of the power coefficient Cp vs. tip speed ratio Xeq are compared in Figs. 20and 21.In the dynamic stall zone (1 < Xeq < 4), both models seem to reproduce quite accurately this zone however in the transition zone (4 < Xeq < 6), both models underestimate the maximum value of the power coefficient.For the unstalled region, it is obvious that the viscous code offers better predictions than CARDAAV code. Figure 21 was reproduced mainly from the reference of Touryan et al.
[1987] for a rotational speed of 50.6rpm.This figure offers a more interesting comparison between several aerodynamic models and the viscous code.These aero- dynamic codes are: double-multiple streamtube model (CARDAAV, CARDAAX)of Paraschivoiu [1988], the model based on vortex theory (VDART3) developed by Strickland [1979] and the model based on local circula- tion theory (MCL) developed by Mass6 [1986].First, we   have to notice that this figure shows results of two versions of the CARDAAV code: one using the indicial model and the other using Gormont model as dynamic stall model.Moreover, the second version of CARDAAV takes into account secondary effects which are important at high tip-speed ratio values.Consequently, differences are noticed between predicted results of both versions of CARDAAV.In the dynamic stall regime (Xeq < 4), the coupling CARDAAV/Indicial code seems to give better results than CARDAAV/Gormont code.According to experimental data, the viscous model provides a good representation of the power coefficient in all ranges of tip-speed ratio Xeq even if the Local Circulation method seems to be the best model that reproduces quite accu- rately the distribution of the performance coefficient C P in almost all operational ranges of Xeq.
The dynamic stall phenomenon is clearly visible on the power curves characterized by a plateau at high speed (Figs.22 and 23).Comparison of rotor power values predicted by both models and experimental data is presented in Figures 22 and 23.For low wind velocities (Veq < 10m/s), both models seem to give a good approximation of the experimental values of rotor power, though in the dynamic stall regime (Veq > 12.5m/s), the CARDAAV code using the indicial model does not predict an almost constant value and keeps on increasing.The viscous model predicts a plateau close to the experimental one for a rotational speed of 42.2rpm.However, for a rotational speed of 50.6 rpm, Fig. 23 shows that predicted values of the power are underesti- mated in the dynamic stall regime even if the plateau still exists.

CONCLUSION
The main objective of this study was to develop a new computational procedure to analyze the global perfor- mance and blade loads of vertical-axis wind turbines.This has been accomplished with a numerical solver based on the solution of the steady, incompressible, laminar Navier-Stokes equations using Finite Vol- ume Method (FVM).Dynamic stall effects were simu- lated by an indicial model (as dynamic stall model) which tackles this problem at a more physical level of approximation but still in a sufficiently simple manner.Analysis of the computed results presented in the previous section has demonstrated certain success of the viscous solver coupled with the indicial model regarding its ability to accurately predict the aerodynamic charac- teristics of Darrieus wind turbine.The results also tend to support the fact that, among other semi-empirical mod- els, the indicial model offers the best representation of dynamic stall for the prediction of wind turbine performance.
As a follow-up, considerable improvement could also be achieved for future predictions.For instance this numerical solver could be easily extended to include the turbulent nature of the wind by using a stochastic model (see reference of Brahimi and Paraschivoiu [1992]).Tip-Speed Ratio (toReq/U) velocity components in the r,to,z directions, m/s absolute velocity of the wind, m/s component of Vab. in the n-direction, rn/s relative velocity, rn/s freestream wind velocity, m/s time averaging factor (NCpVrelA0/4'rr) tip speed ratio at the equator (toReq/U) FIGUREGeneral configuration of Sandia 17-m VAWT.

FIGURE 4
FIGURE 4 Typical control volume cells for u, v, w and p(3-D view).

FIGURE 5
FIGURE 5 Typical curve of the position of the flow separation point function of o.

4MachFIGURE 6
FIGURE 6 Critical normal force coefficient (CN1) for the onset of leading edge separation function of the Mach number.
FIGURE 9 Effect of grid points in 0-direction on power coefficient.
FIGURE 15 Normal force coefficient vs. azimuth angle at 38.7 rpm and Xeq 4.60.
blade airfoil section lift coefficient of blade airfoil section lift static curve slope normal force coefficient of blade airfoil section power coefficient of the turbine tangential force coefficient of blade airfoil section flow separation point (% c) blade length, m Mach number coordinate system attached to the blade number of blades static pressure, Pa cylindrical coordinate system attached'to the blade rotor radius at a certain level, m rotor radius at the equator, m nondimensional time (t(1-M2)2V/c) rotor swept area, m source terms in the momentum equations time, .
Akins concluded in his analysis local angle of attack, deg.static stall angle, deg.local normal angle, deg.dynamic stall reattachment factor, deg.