Rotor aerodynamic power limits at low tip speed ratio using CFD

When investigating limits of rotor aerodynamic models, the Betz limit serves as a solid marker of an upper limit which no model should be able to exceed. A century ago Joukowsky (1912) proposed a rotor aerodynamic model utilizing a rotating actuator disc with a constant circulation. This model has since then been the subject of much controversy as it predicts a power performance that for all tip speed ratios exceeds the Betz limit and which goes to infinity when the tip speed ratio goes to zero. Recently, it was demonstrated that the Joukowsky model is fully consistent with the inviscid Euler equations and that the apparent inconsistency partly can be explained by the lack of viscous effects (Sprensen and van Kuik [4]). However, even including a term to account for the effect of viscosity at small tip speed ratios, the model still predicts a power yield that exceeds the Betz limit. In the present work we study in detail, using a CFD actuator line model, the flow behavior for rotors at small tip speed ratios. It is shown that the excessive swirl appearing towards the rotor center at small tip speed ratios generates vortex breakdown, causing a recirculating zone in the wake that limits the power yield of the rotor. The appearance of vortex breakdown has a similar effect on the flow behavior as the vortex ring state that usually appears at higher tip speed ratios. Limits to where vortex breakdown might occur with tip speed ratio and rotor loading as parameter are investigated and presented in the paper. The limits found correspond to well-known criterion for vortex breakdown onset for swirling flows in general. By applying a criterion for vortex breakdown in combination with the general momentum theory, the power performance always stays below the Betz limit.


Introduction
Independent of the various theories, the upper limit for the aerodynamic power performance of wind turbine rotors tends to converge towards the Betz limit through the normal range of operating tip speed ratios. Here tip speed ratio is defined as λ = ΩR/U o where R is the rotor radii, Ω the angular velocity and U o the freestream velocity. At low and very low λ-values the model of Joukowsky [1] predicts a power coefficient, C P , which always is above the Betz limit and which increases to infinity for λ → 0. Inherent in the model is also an excessive swirl and infinite axial loading towards the rotor center, which is clearly non-physical. However, as shown by Sørensen and van Kuik [4], the issue is easily solved by replacing the inner part by a vortex core. CFD analysis of the model limits permits investigation of the flow field associated with the loading prescribed by the model e.g. separation or other flow phenomenon. In a recent study (see Sørensen and Mikkelsen [2]) an axi-symmetric CFD actuator disc model formulated in streamfunction-vorticity variables, was used to investigate model limits associated with the onset of vortex breakdown, which appear in the wake of an actuator disc operating at low tip speed ratio. In parallel, the general momentum theory combined with known criteria for limits in swirling flows was compared to the axi-symmetric CFD prediction of rotor-induced vortex breakdown in the wake. The result showed that vortex breakdown can be predicted numerically as part of the wake of low tip speed ratio operating rotors. It was further shown that a simple criterion for vortex breakdown used in combination with the general momentum theory, by limit, the power so that it always stays below the Betz limit. The present investigation seeks to extend the axi-symmetric CFD to a fully 3D analysis representing the rotor loading using the actuator line technique rather than an axi-symmetric actuator disc. Thus, the wake is represented by a system of distinct helical tip and root vorticies instead of an axi-symmetric vortex sheet.

Rotor modeling
The governing equations for the rotor modeling is here shortly recaptured, leading to the applied line/disc loadings used in the CFD simulations. For a more detailed derivation, see Sharpe [3] or Sørensen and van Kuik [4]. Applying usual actuator disc terminology [5], let U o , u D and u 1 be the axial freestream, disc and far wake velocity, u θ the swirl velocity, r and R the local and tip radii, Ω the angular velocity and λ = ΩR/U o the tip speed ratio, then the interference factors a and a ′ are introduced as Considering the general momentum theory with the rotor represented by an axi-symmetric actuator disc loaded with constant circulation Γ = 2πru θ , the dimensionsless circulation q is introduced as where a ′ t refers to the tip or edge of the AD. From mass, momentum and energy conservation the following two equations can be derived: where limits for b ∈ [0, 1]. From the definition of the power coefficient and conservation of angular moment of momentum, we get with A = πR 2 as the disc area. Eqs. (3) and (4) constitute a closed set of equations that may be solved for a and q. Combining the three equations and eliminating either a or q gives the explicit expressions for a and q in terms of b and λ as and a = 3 In order to apply disc or line loading in the CFD domain the surface forces f r , f θ , f z are given as (see Sørensen and Kuik, 2011) For the Joukowsky rotor this leeds to u θ ∼ 1/r thus requiring that a limiting function should be introduced with a cut-off radius, r c , such that u θ = 0 at r → 0, corresponding to an inner non-lifting part of the rotor. Defining a linear azimuthal velocity distribution for r < r c and letting it be inversely proportional to the radius for r ≥ r c , we get the following dimensionless azimuthal velocity distribution, with δ = r c /R. However, a smoother variant of eq.(8) of the form (see Delery [6]) is applied. Inserting eq.(7) into eq.(9) results in the following distribution of surface forces, These are the forces employed in the CFD computations. The computations may be carried out by using λ and b as independent variables and employing eqs. (5) and (6) to compute the associated a and q values. However, it shall be emphasized that a and q in the numerical simulations need not to be coupled through the parameter b. Indeed, the parameters in eq.(10) can be chosen arbitrarily, e.g. the azimuthal force distribution could have been specified through a.

Rotor induced vortex breakdown
It is well-known that strongly swirling flows may encounter vortex breakdown, which is manifested as an abrupt change in flow topology associated with reversal of the axial flow component. Vortex breakdown may appear on delta wings, in combustion chambers and in cyclone separators. For flows about rotors vortex breakdown has the same effect on the flow state as the so-called vortex ring state, which appears at high loadings (see e.g. Sørensen et al., [7]). Most theories for developing criteria for the onset of vortex breakdown are based on stability analysis of velocity distributions of the form of eqs. (1) and (2), i.e. a constant axial flow associated with a swirling flow with constant circulation. Since vortex breakdown is associated with the presence of reversal flow, it is clear that this, independent of loading and torque, limits the possibility of the rotor to enhance further the power production. The common way of characterizing the strength of a vortex is through the definition of a swirl number. Defining the swirl number as the angular momentum divided by axial momentum, and introducing the velocity profiles given in eqs. (1) and (2), we get As a general rule, vortex breakdown occurs when the swirl number reaches some critical value.
One of the very first theoretical investigations of vortex breakdown is due to Squire [8], who derived the simple criterion that vortex breakdown may appear when the maximum angular velocity is greater than 1.2 multiplied by the axial velocity. Later experiments on vortices have shown that the phenomenon in practice appears when the angular velocity is about 1.4 greater than the axial velocity, Delery [6]. Assuming that the lifting part of the rotor disc extends from r = r c to r = R, we get the following criterion for vortex breakdown, The evaluation of this criterion is investigated using CFD with rotor loading applied according to the above equations.

Numerical modeling
The Navier-Stokes computations are carried out using the 3D flow solver EllipSys3D [9; 10]. The code solves the discretized incompressible Navier-Stokes equations using a finite volume approach. In the computations, the flow field around the wind turbine rotor is simulated using the actuator line (ACL) technique [11], where the presence of rotor in the computation is replaced by the body forces. The forces along the blade are computed using eq.(10), which are strictly disc loading expressions. However, considering a rotor with a finite number of blades, N b , the forces are distributed and scaled by π/N b to give an equivalent loading to a disc with infinitely many blades. Tip correction is not considered in the present modeling. The cut-off radius should in general be as small as possible, although the region needs to be numerically represented by a sufficient amount of grid points. Initial tests suggest that δ ≈ 15% gives good numerical convergence with the chosen grid and line resolution.

Result and discussions
Computed numerical results are presented in the following. The tip speed ratios investigated range from λ = 0.5 to λ = 2.0, and the value of b is adjusted until the onset of vortex breakdown is present in the wake.   Figure 5 depicts the normalized local thrust coefficient C T (r) for the cases considered, showing that in general the thrust exhibits a C T ∼ 1/r behavior. The C T = 8/9 level of the Betz optimum is also indicated in the figure. The impact of the cut-off function is also evident as it is reduced to zero at the center. The equivalent local C P (r) (not shown) is constant outboard and likevise decreasing towards the center, not exceeding the Betz limit at anytime for the computed cases.

Rotor loading
The finite value at the tip has been maintained in the present modeling although the loading at the tip for a line loading should decrease to zero. Further work is needed on this part. With

Limits to the thrust and power coefficients
The impact of vortex breakdown in the wake of a rotor is a limiting factor for the extraction of kinetic energy of the incoming flow, as the entire wake has shifted regime or state. In figure 8 the C T and C P values for onset of vortex breakdown associated with b and λ values are inserted, which also include axi-symmetric prediction from [2]. The lines included are based on applying the limiting swirl criteria of eq. (12) combined with the general momentum theory as previously presented. Remarkably, the different models agree extremely well, all agreeing on limiting the power coefficient towards zero for λ → 0. The δ-value used for the full CFD is chosen to δ = 0.15  whereas for the axi-symmetric results, δ = 0.2 was used in addition to a much higher resolution. The fit between the high resolution axi-symmetric results and the general momentum theory with a swirl limiter for δ = 0.2 suggest that the axi-symmetric conditions assumed by both methods is important. Higher resolution computations should be carried out to investigate further the numerical dependencies of resolution and different values of δ.

Conclusions
The vortex breakdown phenomenon associated with a rotor including a swirling wake has been investigated, showing firstly the presence of vortex breakdown at low tip speed ratios and next, the limiting impact on the maximum achievable power for low tip speed rotors. The onset of vortex breakdown compares elegantly, using either full 3D or axi-symmetric CFD, to the general momentum theory including limiting criteria for swirling flows. Further investigations are to be conducted on the detailed behavior of the wake for a more detailed characterization of rotor induced vortex breakdown.