A Novel Correlation Model for Horizontal Axis Wind Turbines Operating at High-Interference Flow Regimes

Abstract

Driven by economics-of-scale factors, wind-turbine rotor sizes have increased formidably in recent years. Larger rotors with lighter blades of increased flexibility will experiment substantially higher levels of deformation. Future turbines will also incorporate advanced control strategies to widen the range of wind velocities over which energy is captured. These factors will extend turbine operational regimes, including flow states with high interference factors. In this paper we derive a new empirical relation to both improve and extend the range of Blade Element Momentum (BEM) models, when applied to high interference-factor regimes. In most BEM models, these flow regimes are modeled using empirical relations derived from experimental data. However, an empirical relation that best represents these flow states is still missing. The new relation presented in this paper is based on data from numerical experiments performed on an actuator disk model, and implemented in the context of a novel model of the BEM family called the DRD-BEM (Dynamic Rotor Deformation—BEM), recently introduced in Ponta, et al., 2016. A detailed description of the numerical experiments is presented, followed by DRD-BEM simulation results for the case of the benchmark NREL-5MW Reference Wind Turbine with this new polynomial curve incorporated.

1. Introduction

As wind energy penetration continues to grow, there is a need to improve both quality and quantity of power being generated by wind turbines [1,2]. Wind energy generators now need to efficiently capture energy from the wind over a wider range of wind velocities at a frequency compatible with the national grid. Reduction in proportion of conventional generators in power systems also means increasing participation of wind turbines in maintaining system stability making cost effective power capture extremely relevant in the current scenario. Economics of scale suggests that one way to maximize power captured would be to scale up turbine size [3] and as recent trends suggest, that is exactly how the industry has responded. State of the art turbines have reached a rotor diameter of over 120 m and are capable of generating 7.5 MW of power [4]. Larger turbines with rotor diameters ranging from 150 m to 180 m, generating up to 8 MW are already in the pipeline [5]. This radical increase in turbine size has also led to an increase in turbine loads and manufacturing costs [3], thereby, prompting the need for lighter blades with greater flexibility. These demands have led to the development of innovative and complex control mechanisms using both active and passive techniques to reduce blade loads and improve the quality of power generation.

With increasing perspectives on the use of adaptive blades (see Ponta et al. [6] for a detailed review of the topic), and the introduction of innovative control strategies, wind turbine rotors will likely undergo greater levels of deformation during operational regimes over a wide range of wind velocities. These factors can lead to rotor blade sections operating at flow states characterized by high axial induction factors. These may often occur in transient states during control actions [7,8]. Other scenarios which could drive blade sections into this regime would be when the turbine operates in part load conditions during high winds storing kinetic energy [9,10,11,12,13], sudden drop in effective or actual wind [14], and failure of control actions. Sebastian [14] explains how off-shore wind turbines can frequently encounter such flow states.

Under these conditions, the momentum balance on a streamtube control-volume, which the widely used Blade Element Momentum model (BEM) is based upon, can no longer be used (this topic will be covered in detail in Section 2.5). In most BEM schemes, the modeling of this high axial induction factor flow regime is based on an empirical relation derived from experimental data. This empirical method was introduced to the BEM model by Glauert [15], using data from Lock et al. [16], Lock and Bateman [17], Munk [18]. Several modifications to this empirical fitting have since been proposed by Buhl [19], Burton et al. [20], Wilson [21], and Ponta et al. [22], among others. However, it is still unclear as to which empirical relation provides the most accurate description of this operational state. One of the reasons could be the uncertainty in the experimental results, especially in the very high end of the interference ranges. The original experimental results show a lot of dispersion which can be attributed to sensor limitations of that era, as well as the unsteady nature of these operational regimes.

It must further be pointed out that the original experiments were performed on propellers with constant chord and incidence in closed wind tunnels. The results have not yet been validated on conventional blades in open jet wind tunnels. With the uncertainty in effect of tunnel interference, several authors, including Glauert [15], have questioned the applicability of these experiments on conventional full-size rotors in the absence of boundary layer effects of wind tunnels Johnson [23]. As mentioned earlier in this section, turbines now-a-days will transition through high axial induction factor flow regimes more frequently. As the dynamic response of the extremely flexible state-of-the-art blades gains in importance, it becomes imperative to capture blade response accurately at each time step at every section of the blade. An inaccurate representation of these transitory states makes such an analysis unreliable. In an innovative approach to this problem, flow across sections that have entered these high induction factor regimes, is instead, numerically simulated as flow across an actuator disk using turbulence models.

By using such sophisticated Computational Fluid Dynamics (CFD) tools, in this paper, we aim at a better representation of these high-interference flow states, implemented in the context of a novel model of the BEM family called the DRD-BEM (Dynamic Rotor Deformation—BEM), recently introduced by Ponta et al. [6]. For this study, flow across the actuator disk model was simulated using values of thrust coefficient similar to the above mentioned experimental data at high-interference regimes. In addition, low interference regimes were also simulated together with the classic streamtube control volume solution, as a form of validation. The next section will present a brief summary of the DRD-BEM, followed by a section analyzing the physics of the actuator-disk theory, and a further section introducing a new extended polynomial curve to fully characterize these high interference-factor flow states. A detailed description of the numerical experiments is presented, followed by DRD-BEM simulation results for the case of the benchmark NREL-5MW Reference Wind Turbine [24], with this new polynomial curve incorporated into the model. Results will be presented for characteristic scenarios where the turbine enters into flow states with high interference factors.

2. Theoretical Background: The DRD-BEM Model

This section gives a brief account of the DRD-BEM, the simulation model used in this study. A more detailed description can be found in Ponta et al. [22]. For details of the structural model used in the DRD-BEM model, the reader is referred to Otero and Ponta [25]. The level of description provided in this section should be sufficient to make this paper self-contained.

2.1. Historical Context

Among the various approaches used to model the aerodynamics of a horizontal axis wind turbine, the stream tube approach known as the Blade element momentum (BEM) theory is the most widely used aerodynamic model (a detailed description of the classical BEM can be found in Manwell et al. [26] and Burton et al. [20]). However, to account for certain limitations of the BEM theory when applied to real world conditions, there have been several improvements and corrections to the classical BEM model. Among the examples of research carried out on the BEM model over the past decade, there are studies dealing with analysis of coning rotors Crawford [27], inclusion of centrifugal pumping Lanzafame and Messina [28], pressure drop across the rotor and wake expansion Madsen et al. [29], modifications to the dynamic stall model Dai et al. [30], and effect of wake influence Vaz et al. [31]. However, there are still certain assumptions that need to be accounted for in the BEM model.

The motivation behind the DRD-BEM model are the limitations imposed by the underlying assumption in classical BEM that all blade sections are normal to the rotor’s radial direction. As a result of this, the blade section misalignments caused by rotor deformation are not properly represented. This misrepresentation becomes even more pronounced in case of large scale deformations, something that can be expected given the flexibility of modern blades. This issue is tackled by the Dynamic Rotor Deformation (DRD) part of the DRD-BEM. It ensures that in the DRD-BEM approach the rotor deformation is accurately represented owing to the action of a set of orthogonal matrices on the incident velocity and aerodynamic-force vectors. These matrix transformations account for all misalignments caused by various factors like deformation of blade sections, actions of mechanical devices controlling pitching, yawing or azimuthal rotation, changing winds, and several design features like coning or tilting of the rotor. Figure 1 shows schematically the instantaneous position of a blade element at distance $r_h$ from the hub. The radial thickness of the annular actuator swept by that blade element, in a coordinate system aligned with the hub, is given by $\delta r_h$. Also seen in that figure is the blade element’s spanwise length $\delta l$. As will be seen later, the DRD-BEM model judiciously uses these variables to ensure that all instantaneous misalignments are accounted for by constantly updating the area of the annular actuator corresponding to the respective blade element.

2.2. Structural Model: The Blade as a Complex Beam

The DRD-BEM models the structural response by using the Generalized Timoshenko Beam Model (GTBM). The GTBM is a generalized implementation of the Timoshenko beam theory using a variational asymptotic method. It accounts for warping of beam sections by interpolating sectional warping using a 2-D Finite element model to obtain 3-D warping functions. These warping functions are used to asymptotically approximate the 3-D strain energy at each section in terms of the six strain measures including those of the classical 1-D Timoshenko theory and the transverse shear in two directions. The complex geometry and material properties are represented by stiffness matrices obtained a priori for each section. The generalized implementation of the Timoshenko theory ensures a fully populated 6x6 stiffness matrix at each section of the beam coupling all 6 modes of deformation. The 3-D beam problem is thus reduced to a nonlinear 1-D beam. This unsteady nonlinear 1-D beam problem, together with the aerodynamic model, discussed later comprises the full aeroelastic scheme.

This aeroelastic scheme is solved along the beam’s reference axis, L, for every time step using an advanced ODE algorithm. From each time step of the ODE solution the 3-D fields can be recovered using the previously derived 3-D warping functions. As seen in Figure 1, the solution variables first represented in a coordinate system along L are transformed into the instantaneous coordinate system along the new configuration l. This is done with the help of the $C_{Ll}$ matrices, updated at every time step of the ODE solution of the structural model. Thus the instantaneous position of each blade section is accurately tracked by this intrinsic coordinate system. As a result the aerodynamic part of the aeroelastic analysis gets this updated information at every time step for each blade section. At the same time the results of the aerodynamic solution constantly update the aerodynamic loads acting at each blade section as inputs to the structural model for every time step.

2.3. The DRD-BEM Model

As mentioned earlier, a detailed description of the DRD-BEM procedure is presented in Ponta et al. [22]. The first step in the DRD-BEM procedure is to account for the effect of the annular actuator on the undisturbed wind velocity in the hub coordinate system, ${{\mathbf{W}}_{\infty }}_h$. First, the wind velocity which is originally given in a coordinate system aligned with the wind itself as ${{\mathbf{W}}_{\infty }}_{wind}$ is transformed to ${{\mathbf{W}}_{\infty }}_h$ using a set of orthogonal matrices as follows

$${{\mathbf{W}}_{\infty }}_h=\left({\mathbf{C}}_{{\theta }_{az}}{\mathbf{C}}_{{\theta }_{tlt}}{\mathbf{C}}_{\Delta {\theta }_{yaw}}{{\mathbf{W}}_{\infty }}_{wind}\right).$$

(1)

The first 2 coordinate transformations in Equation (1) represented by transformation matrices ${\mathbf{C}}_{\Delta {\theta }_{yaw}}$ and ${\mathbf{C}}_{{\theta }_{tlt}}$ account, respectively, for the misalignment at the nacelle due to yaw and the misalignment brought about by rotor tilt. The third orthogonal matrix, ${\mathbf{C}}_{{\theta }_{az}}$, rotates the blade around the main shaft to its new position at that instant transforming the wind velocity into the hub coordinate system h (see Figure 2a).

As this undisturbed wind flows across an annular actuator, it is deflected by the forces exerted on it by the actuator. Therefore, the wind velocity at the annular actuator, represented by ${\mathbf{W}}_h$, is given by,

$${\mathbf{W}}_h=\left[\begin{array}{c}{W}_{\infty h_x}(1-a)\hfill \\ W_{\infty h_y}+\mathsf{\Omega }{r}_h{a}^{\prime }\hfill \\ W_{\infty h_z}\hfill \end{array}\right],$$

(2)

where a and $a^{\prime }$ are the induction factors representing, respectively, the fraction of induced velocity in the axial and tangential directions of the hub coordinate system. $r_h$, as shown in Figure 1, gives the instantaneous radial distance between the hub and the blade section in the hub coordinate and $\mathsf{\Omega }$ represents angular velocity of the turbine rotor.

${\mathbf{W}}_h$ is projected onto a coordinate system aligned with the blade section through a series of coordinate transformations. The first orthogonal transformation matrix is ${\mathbf{C}}_{{\theta }_{cn}}$ which accounts for coning and ${\mathbf{C}}_{{\theta }_p}$ which performs a rotation about the blade’s pitch axis. This projects the wind onto blade coordinate system defined (see Figure 2b). ${\mathbf{C}}_{Lb}$ projects the velocity vector onto the intrinsic coordinate system of the blade sections along the undeformed reference line L. Finally, ${\mathbf{C}}_{lL}$ is used to represent the undisturbed wind velocity in a coordinate system aligned with the instantaneous deformed configuration along reference line l. To get the final wind velocity as seen by the blade section we need to add two more components. First is the vibration velocities ${\mathbf{v}}_{str}$ of that section as given by the structural model. The other component is the mechanical velocities, ${\mathbf{v}}_{mech}$, accounting for yawing, pitching, azimuthal rotation or any other mechanical motion associated with that blade section. Since both these components are expressed in the l coordinate system, they can be added directly to the wind velocity obtained after all of the above mentioned coordinate transformations to give

$${\mathbf{W}}_l=\left({\mathbf{C}}_{lL}{\mathbf{C}}_{Lb}{\mathbf{C}}_{{\theta }_p}{\mathbf{C}}_{{\theta }_{cn}}{\mathbf{W}}_h\right)+{\mathbf{v}}_{str}+{\mathbf{v}}_{mech}.$$

(3)

The relative wind velocity magnitude in the plane of the blade section is given by $W_{rel}=\sqrt{{W_{l_x}}^2+{W_{l_y}}^2}$. $W_{rel}$ along with angle of attack $\alpha$ are used to compute the aerodynamic lift and drag per unit length of that section giving the aerodynamic force vector

$$\mathit{\delta }{\mathbf{F}}_{rel}=\left[\begin{array}{c}d{F}_{lift}\\ d{F}_{drag}\\ 0\end{array}\right]\delta l,$$

(4)

The aerodynamic load is then obtained in the hub h coordinate system by applying the inverse of the orthogonal matrices used earlier, this time, however, in the reverse order

$$\mathit{\delta }{\mathbf{F}}_h={\mathbf{C}}_{{\theta }_{cn}}^T{\mathbf{C}}_{{\theta }_p}^T{\mathbf{C}}_{Lb}^T{\mathbf{C}}_{lL}^T{\mathbf{C}}_{Lthal}\mathit{\delta }{\mathbf{F}}_{rel}\delta l,$$

(5)

Here ${\mathbf{C}}_{Lthal}$ projects the lift and drag forces onto the chord-normal and chord-wise directions in the l coordinate system.

This force is then equated with the forces obtained by performing the momentum balance in the hub coordinate system. The component normal to the annular actuator, $\delta F_{h_x}$, equated to the momentum change on $W_{\infty h_x}$ and the tangential component $\delta F_{h_y}$, is equated to the momentum change on $W_{\infty h_y}$. An interference pattern is obtained, enabling us to use an iterative scheme to obtain the two induction factors a and $a^{\prime }$. These induction factors are then used to obtain ${\mathbf{W}}_h$ and then following all of the steps explained above compute the distributed aerodynamic forces on the entire blade per unit span length. $\mathit{\delta }{\mathbf{F}}_l={\mathbf{C}}_{Lthal}\mathit{\delta }{\mathbf{F}}_{rel}$ gives the chord-normal and chord-wise aerodynamic loads and $d{M}_{aer}=\frac{1}{2}\rho C_m{W}_{rel}^2{c}^2$ computes the aerodynamic moment around the first axis of l, per unit span length at each blade section. Here $C_m$ is the aerodynamic pitch coefficient of the airfoil at the corresponding angle of attack, $\rho$ is the density of air and c represents chord length of that section. The effect of gravity on the distributed forces and moments is also accounted for with the help of inertia properties of the equivalent beam represented by inertia matrices obtained as part of the GTBM process of dimensional reduction described in Section 2.2. These forces provide the input to the GTBM structural model.

2.4. Dynamic Update of Corrective Factors

To account for the various assumptions made in the blade element momentum theory certain correction factors have to be used. In DRD-BEM these factors, unlike the classical implementation of BEM, are updated dynamically. These include the rotational-augmentation effects as proposed by Du and Selig [33] and Eggers [34], dynamic-stall effects based on the works of Leishman and Beddoes [35,36,37], extrapolation of aerodynamic-coefficient data proposed by Viterna and Janetzke [38], tip and hub loss corrections (see Section 3.8.3 of Manwell et al. [26], and Ponta et al. [22]), and the tower-interference models proposed by Bak et al. [39] and Powles [40]. The DRD-BEM also incorporates the flexibility of multiple data sets for aerodynamic coefficients of different airfoils associated with factors varying the aerodynamic situation at the rotor. These factors could include varying Reynolds number, effects of flow control devices, etc.

2.5. The Glauert Correction

Another important modification needed in any model based on momentum theory is the correction of the thrust on the annular actuator when operating in the so-called turbulent-wake state. This correction, is the focus of the present paper, and plays a key role when the turbine operates at axial induction factors a greater than 0.5 (in practical implementations, this limit is lowered to about 0.3 to 0.45, depending on the corrective curve adopted).

The classic streamtube-momentum solution for the actuator-disk representing the thrust coefficient $CT$ as a function of a gives a parabolic curve of the form $CT=4a(1-a)$, which is depicted in Figure 3. At $a=0.5$, that parabola reaches its vertex. Beyond that the basic assumptions of momentum theory on a streamtube control-volume become invalid, as the flow at the pseudo-infinite position, in the far wake downstream, would theoretically stagnate and start to propagate upstream. Physically, this flow reversal cannot occur, and what actually happens is that flow re-enters from outside of the wake, creating vortex structures and increasing the turbulent activity behind the rotor’s disk. This slows down the flow passing through the rotor, while the thrust continues to increase.

This work focuses on finding out a new extended empirical relation to compute the thrust for flow states with a high induction factor in the turbulent-wake state and the vortex-ring state. As explained in Section 1, instead of using momentum balance, Glauert [15] proposed finding $CT$ from an empirical fitting to data from experiments on propellers provided by Lock and Townend [41]. As seen in Figure 3, Glauert’s empirical fit is tangent to the $CT$ curve at $a=0.4$. Alternative fitting functions to the same set of data were proposed by other researcher works like Wilson [21] and Burton et al. [20]. As noted by Buhl [19], a discontinuity between the $CT$ function and the fitting curve can be observed upon the application of tip and hub losses. This makes convergence difficult to achieve for models where the induction factors are calculated iteratively. Buhl [19] solved the problem by using an empirical relation where the tangential match was enforced regardless of the tip and hub loss corrections. Another solution was proposed by Ponta et al. [22] in the form of a Power-Law fitting. The latter fitting curve, while significantly reducing the error of approximation, also avoids the discontinuity problem, since it always intercepts the stream-tube $CT$ function whether or not the tip and hub losses are taken into account. The DRD-BEM model has the capacity to switch between any of these correction methods.

As mentioned in Section 1, the set of experiments reported by Lock and Townend [41] was recreated by numerically solving an actuator-disk set-up in ANSYS-Fluent. This numerical set-up also allowed to extend those experiments to the low interference conditions where $a<0.4$, and regions of very-high values of interference factors in the vortex-ring state, characterized by $a>0.9$, which go beyond the range of Lock and Townend [41] experiments. This new set of numerical data was then used to define a new empirical curve which can be used to find the thrust at high induction factor regimes. This new curve covers any induction factor between $a=0.4$ and $a=1.4$, and merges seamlessly with the $CT$ versus a curve from the streamtube theory for values of $a<0.4$. The following sections discuss the set of numerical experiments performed, together with the analysis of physical aspects associated with the actuator-disk model operating at flow states with high induction factors.

3. The Physics of the Actuator Disk Model

As it was mentioned in previous sections, future turbines are expected to undergo larger levels of deformations, which in turn will induce higher levels of misalignment of the airfoil blade sections which respect to the incoming flow (see Otero and Ponta [42] for a detailed description of blade-section misalignment effects). These misalignments lead to changes in the effective interference at the rotor annular areas swept by those misaligned blade sections, often leading to higher axial-induction factors at various portions of the rotor. The nature of the DRD-BEM model is such that it captures these misalignments with great accuracy. This means that the simulation will detect flow states not lying in the standard operating regimes more precisely. As an example, we shall analyse such a possible scenario based on a DRD-BEM simulation of the benchmark NREL-5MW RTW turbine. Figure 4 shows simulation results of a de-loaded NREL-5MW turbine operating in the control region at a wind speed of 8 m/s, where the turbine is optimized for power extraction [24]. The turbine is de-loaded by increasing the rotor RPM from approximately 9.2 rpm to 15 rpm to obtain a 10% margin of reserve power from the turbine. Jonkman et al. [24] gives the range of optimal rpms for an NREL-5MW turbine with varying wind speeds. A power reserve of 10% is based on Ahmadyar and Verbič [43] which analyses de-loading an NREL-5MW reference turbine and also on Venne et al. [44] and Chang-Chien et al. [45] where control strategies based on de-loading the turbine along with realistic limits of rpm variation used for such strategies are analyzed. An axial induction factor greater than $0.4$ can be observed and this, as it was mentioned before, lies in the turbulent-wake state, and requires the application of the Glauert correction.

3.1. The Experimental Data

As the axial induction factor varies, the rotor goes through several interference states. Figure 5 depicts the different states as compared to the theoretical parabolic curve obtained by the momentum theory. The experimental results used by Glauert [15] are also plotted.

For values of $a<0$, the rotor acts as a regular propeller, transferring energy to the flow instead of extracting energy from it, with a negative thrust (e.g., with the thrust pointing in the upwind direction).

The normal operating state of the wind-turbine rotor, known as the windmill-brake state, is where the rotor maximizes energy extraction from the flow. In particular, the ideal condition occurs where the axial induction factor a is equal to $1/3$, which is characterized by a power coefficient $C_p=16/27$ (or approximately 59.26%). This represents an optimum value of energy-conversion efficiency for an ideal wind-energy harvester, usually referred to as the Betz’s limit (in honour of German physicist Albert Betz who proposed the analysis in 1919 [20,26]). Wind-turbine designs target this ideal as their nominal (or “rated”) operational conditions. Betz’s $a=1/3$ ideal condition lies within the windmill-brake state.

With an increase in tip speed ratio either due to a drop in wind velocity, a change in the turbine’s rotational speed, or a change in rotor orientation with respect to the wind, the axial load on the rotor disk increases in proportion to the wind velocity. This leads to an increase in the axial induction factor. Initially there is incipient recirculation in the rotor’s wake with re-entrance of flow from outside of the wake, which mixes with the slow moving wake flow. This is the so called turbulent-wake state that we described before. In terms of the actuator-disk analogy, in this state the rotor disk begins to act like a semi-porous imaginary surface with decreasing permeability. As seen in Figure 5, the turbine is said to be in this state when the axial induction factor a lies approximately between 0.5 and 1. The experiments reported by Lock et al. [16] are one of the earliest studies of high interference states where the author goes on to propose a possible empirical correction. Several others have provided a detailed analysis of the high axial induction factor regimes. These include Burton et al. [20], Stoddard [48,49], Eggleston and Stoddard [46], and Wilson and Lissaman [50]. In particular Glauert [15] forms the basis of the correction model used in BEM. More recent are the works of Sørensen et al. [51] presenting a numerical analysis of an actuator disk model in this regime, and Johnson [23] presenting a detailed analysis of several noteworthy experiments that have provided data for these flow states.

For values of a higher than 1, the rotor enters the vortex-ring state, where the recirculation starts to assume a toroidal pattern approximately centered at the rotor’s disk. As a increases beyond 1, the flow direction across the rotor starts to reverse. In this condition the turbine acts like a propeller operating in brake mode (that is the case of a propeller aeroplane in the landing run after touching down), and begins to transfer energy into the flow with an increased thrust in a downwind direction. This flow state is referred to as the propeller-brake state.

3.2. The Numerical Approach

As mentioned previously, this experimental approach has been the most widely accepted method to account for rotor interference in high a regimes, with relatively slight modifications introduced to improve Glauert’s original empirical relation. One factor to point out is the substantial amount of dispersion seen in the experimental data used by Glauert. This could be attributed to the limitations of the sensors from that era, as well as, to the unsteady nature of the flow in the close surroundings of a real rotor. On the other hand, the interference computation in BEM theory is based on an idealization steady-in-the-average evaluation of the momentum change in a streamtube control-volume, which differs from the situation across a real rotor.

However, it can be argued that an actuator-disk model is still a very suitable representation of the rotor interaction with the flow, and that is the main reason the method has been widely employed for so many years. Therefore, Glauert’s basic concept of replacing the momentum balance in the streamtube by a mathematical expression giving a functional relation of $CT$ versus a (that represents a more complex turbulent-flow interference) is still valid. Moreover, several recent experimental studies on actuator-disk models have been carried by Castro [52], Hong [53], Lignarolo et al. [54], Batten et al. [55], and Harrison et al. [56]. These studies, however, concentrate mostly on wake analysis instead of the analysis of conditions at, or in the close surroundings of the disk, which is what is needed for a BEM model.

Wake measuring experiments benefit from the fact that non-steady effects created in the near surroundings of the rotor get attenuated by the distance. On the other hand, measuring flow properties in the close surroundings of the rotor’s disk is much more challenging. Thus, for conditions at the actuator itself, a numerical analysis would be better suited, as the actuator-disk is essentially an imaginary semi-porous surface that induces a pressure drop. These conditions are relatively easy to set-up in a numerical simulation.

The numerical model has been developed in ANSYS-Fluent. It consists of a 1m radius actuator disk in a 2-D axisymmetric domain. The disk is essentially an infinitesimally thin boundary lying at an interface across cells. This infinitesimal interface represents an actuator disk by virtue of a discontinuous pressure drop acting across it. For this experiment, the $CT$ acting on the disk was varied by varying this pressure drop. A relatively large domain size was selected to ensure that the external boundaries do not affect the solution at the actuator disk, especially when the flow enters the extremely unstable vortex ring state. The overall domain size was 250 m upstream and 350 m downstream of the actuator disk and extended 400 m above the horizontal axis of symmetry. Figure 6 shows a close-up view of the the local domain surrounding the actuator-disk region in one of the meshes used in this study.

The typical boundary conditions of Velocity inlet and pressure outlet were used for the external boundaries. A free-stream velocity, U, of 1 m/s was used. An unstructured triangular mesh using the ANSYS meshing tool was used to discretize the domain. The gradient adaption algorithm of Fluent was used to dynamically refine the mesh at each iteration using the curvature of velocity as the refinement criteria.The Navier-Stokes equations for the turbulent flow across this disk were solved in this discretized domain by adopting a turbulence model of the $k-ϵ$ family called the realizable$k-ϵ$. As explained in ANSYS [57] and Menter [58], the $k-ϵ$ family of models are based on the evaluation of two properties of the turbulent flow: the turbulence kinetic energy k, and its dissipation rate $ϵ$. The $k-ϵ$ models are more suitable to simulate turbulent properties on the far-field flow regions than other turbulence models like the ones in the $k-\omega$ family. That is why the $k-\omega$ SST model gradually changes from a $k-\omega$ approach in the boundary layer to a $k-ϵ$ formulation in the free-shear regions. It can be reasoned that the actuator disk, because of its semi-porous surface nature, mostly deals with free-shear effects at the zone of interaction between the disk and the incident flow than with near-wall type conditions, therefore, a $k-ϵ$ model was chosen.

In particular for the case of boundary-free shear flows with high mean shear rate, and also for flows with secondary flow features, the realizable version of the $k-ϵ$ model performs better than other models of the $k-ϵ$ family [57,59]. These are the flow conditions expected from the actuator-disk model, which are characterised by a sudden pressure drop. A thorough comparison was carried out between the $k-ϵ$ realizable model with the other turbulence models available in Fluent. The flow topology observed in the simulations supports the use of a $k-ϵ$ realizable model. Figure 7 shows an expanding wake behind the actuator disk for an axial induction factor $a=0.41$ with $CT=1.01$. The first signs of vortical wakes with recirculating patterns and re-entry flow can be seen in the zoomed-up view in Figure 7b. This pattern is absent in the simulation results for $CT=1$. As explained earlier, this flow pattern signals the beginning of the turbulent-wake state, on which the streamtube control-volume theory is no longer applicable to the actuator disk model. Significantly, that is precisely the region around which the experimental data compiled in Glauert [15] start to enter the turbulent wake state (see Figure 3).

Furthermore, it was seen that results from the $k-ϵ$ realizable model best matched the data used by Glauert for flow in the turbulent-wake state. It also closely follows the trend of the power-law fitting proposed in Ponta et al. [22] in the region lying between the turbulent-wake state and the vortex ring state. This power-law fitting has previously been shown to minimize the approximation error to the experimental data used by Glauert [15]. Another indication of the the $k-ϵ$ realizable model’s suitability is the fact that the simulation results at low axial interference factors (i.e., $a<0.5$) matches perfectly with the theoretical solution for the streamtube control-volume analysis in classic BEM theory. The streamtube analysis for such low-interference flow regimes has been validated in numerous studies involving both, controlled experiments on propellers and scaled-down turbine models, as well as field measurements on full-scale turbines. All these observations validate the use of the $k-ϵ$ realizable model to solve the actuator disk flow in our numerical experiments.

The values for the pressure drop $\Delta P$ applied across the actuator disk in the simulations were computed from the coefficient of thrust $CT$ expected for the operational range under study. $CT$ is given by:

$$CT=\frac{T}{\frac{1}{2}\rho U^{2}A}$$

(6)

where A is the area of the actuator disk, $\rho$ the air density, U is the free-stream velocity equal to 1 m/s and T the thrust on the disk, which could be written in terms of the pressure drop as $T=A\Delta P$. Therefore, the relation for pressure drop $\Delta P$ across the actuator disk is given by:

$$\Delta P=\frac{1}{2}\rho CT$$

(7)

Before applying the the numerical model to obtain the new set of data, we conducted a thorough test to confirm that the final mesh obtained by using Fluent’s gradient adaption algorithm consistently resulted in an accurate solution. To this end, a mesh convergence analysis was performed on three points selected from the $CT$ curve shown in Figure 3. The $CT$ values selected for the mesh convergence analysis were $0.8$, 1 and $1.4$. This can be seen in Figure 8, which shows a plot of axial velocity at the disk against the mesh density. These cases cover flow states ranging from low axial induction factor windmill brake states to the point were the flow in the turbine wake is well into the turbulent wake state.

The mesh density was varied using two meshing parameters: The Global mesh refinement and the element growth rate around the disk. To vary the growth rate in the proximity of the disk a smaller local domain around the disk was selected within the larger global domain. This domain extended 21 m downstream and 9 m upstream of the disk and 10 m above the axis of symmetry. As a result the increase in element count was observed primarily in the local domain, thereby, avoiding an unreasonably dense mesh in the entire domain. This would increase solution times exponentially without any contribution to solution accuracy, given that outside of the local domain, the flow can be considered as free stream. The criteria for mesh convergence was set at a velocity residual lower than 0.01 m/s between successive mesh refinements. The number of elements was increased from approximately 1000 elements to 40,000 elements, within the local domain of the disk, in steps of approximately 5000 elements. As an example, the local domain highlighted in green in Figure 6 shows a particular refinement of approximately 20,000 elements, with the global domain consisting of 35,000 elements.

It was observed that solution reached the above mentioned convergence criteria within the first couple of refinements for all three $CT$ values. Finally, these axial velocities were compared against those obtained using the Gradient Adaption algorithm in Fluent. The difference in velocities was well within the 0.01 m/s criteria mentioned above. It was, therefore, decided to go ahead with Fluent’s Gradient Adaption method.

Finally, the axial induction factor a was computed from the average flow velocity across the actuator disk obtained from the simulations using the relation $U_{dsk}=U(1-a)$, where U is the free-stream velocity, and $U_{dsk}$ is the area-weighted average of the axial velocities at the interface representing the actuator disk. The axial induction factors were plotted against their respective $CT$ values to obtain the new correction curve.

4. The New Polynomial Correction Curve

Figure 9 shows the numerical results for the $CT$ versus a relation obtained from the numerical experiments previously described, together with a best-fit curve to the data obtained by a least-squares regression of a 4th order polynomial, and the parabolic relation given by the classic streamtube momentum solution. As it could be seen, the numerical experiments not only merge smoothly with the classic streamtube-momentum solution, but also shows a very close agreement with the classic curve all over the range of values where the classic theory is valid.

A very important issue that needs to be considered for any $CT$ versus a correction method that is proposed, is the merging of the correction curve with the streamtube-momentum parabola. In particular, as previously mentioned in Section 2.5, how this merging may produce discontinuities that may affect convergence of the iterative process in the BEM solution when the tip and hub loss factors are taken into account (for a detailed description on subject of tip and hub loss factors see Section 3.8.3 of Manwell et al. [26], and Ponta et al. [22]). In this regard, determining the point of departure of the polynomial fitting curve from the theoretical parabolic curve is worth of further consideration.

As it was discussed in Section 2.5, Glauert’s original empirical fit was designed to be tangent to the $CT$ curve at $a=0.4$, which as noted by Buhl [19], creates a discontinuity between the theoretical parabolic function and the correction curve when tip and hub loss factors are applied. Buhl [19] solved the problem by proposing an empirical relation where the tangential match was enforced at $a=0.4$ regardless of the tip and hub loss corrections. In that process, the initial data points from Glauert’s experimental data set are essentially ignored, and the adapted correction curve is set as a parabola that fits the tangential merging at $a=0.4$, and it is anchored in its other end at $CT$ curve at $a=1$. This procedure is illustrated in Figure 10, where Buhl’s fitting without tip-hub losses is compared with a fitting incorporating a tip-hub loss factor $FF$ of $0.8$ (another example with $FF=0.9$ can be seen in Figure 3). Buhl [19] method is successful in avoiding the convergence problem associated with the discontinuity in the $CT$ versus a curve, but it has the caveat that as $FF$ decreases, the resulting curve does not match the trend of the experimental data (which could be seen in Figure 10), and the error of approximation increases.

Upon a closer examination of the origin of the tip-hub losses, an alternative approach can be proposed. Even though in BEM theory the $FF$ factor is applied to the momentum change in the streamtube wake, the tip-hub losses are, actually, physically associated with the trailing vortices created at the tip and root of the rotor blades, in the same way that they form at the tip and root of an airplane wing. That means that the origin of this effect is connected with the Blade Element part of the BEM theory, and as such, the effect is still going to be present even if the rotor goes beyond the range of the windmill-brake state and enters the turbulent-wake state. Therefore, shifting the entire fitted curve in accordance with the tip-hub loss factor, as it is done with the theoretical parabolic relation, seems like a logical option. By fitting the polynomial curve using the entire data set all the way from $CT=0.2$ to $CT=4$, and not just after $a=4$, and then applying the $FF$ factor, we ensure that the curves always merge gradually in the range between $a=0.28$ and $a=4$. It also ensures that the accuracy of the streamtube analysis, validated many times, is effectively used. In the practical implementation in the context of the DRD-BEM model, we adopted $a=0.28$ as the switching point between the two curves, with the theoretical curve covering the range of $a<0.28$, and the polynomial fitting the range of $a>0.28$. Figure 11 shows examples of the polynomial correction curve for the case of no tip-hub losses, and for the case of a tip-hub loss factor $FF=0.9$.

Implementation in the Context of the DRD-BEM Model

To test the performance of the new polynomial-curve correction, we incorporated it into the general scheme of the DRD-BEM model. We conducted simulations based on two operational scenarios for the NREL-5MW turbine. Both scenarios are characterized by values of the axial induction factor a higher than 0.4 on the blade section located at 90% of the span (i.e., the rotor is operating within the turbulent-wake state).

Of the two scenarios, the first one represents a typical situation during the start-up regime of an NREL-5MW, as reported by Jonkman et al. [24], with a rotor’s angular speed of 7.9 rpm and a wind velocity of 6 m/s. The mean a over 3 rotor revolutions is found to be 0.403. The second scenario represents the situation of an NREL-5MW turbine rotating at its nominal speed of 12.1 rpm, with the wind velocity dropping suddenly from the nominal velocity of 11.4 m/s to 5.8 m/s. The mean a over the same 3 cycles is 0.612. A comparison of the time evolution of the axial induction factors for the two test cases can be seen in Figure 12.

Figure 13 shows examples of the time evolution over the same 3 cycles of three relevant parameters: rotor power, thrust, and displacement of the blade axis at 90% of the span, for the first test case. Figure 14 shows the evolution of the same parameters for the second test case.

5. Conclusions

In this paper, a new evolution of the BEM theory for wind turbines operating in high axial interference regimes is presented. This analysis is based on numerical experiments on an actuator disk model. A primary reason for exploring these operational regimes is the increasing complexity of wind turbine control strategies, and the increasing complexity of the structural response of wind turbine blades. This has resulted in wind turbines often operating at regimes not usually found earlier. With state of the art and future turbine blades experiencing higher levels of deformation, flow states at higher axial induction factors are (and will be) encountered more frequently and for longer durations. Figure 4 and Figure 12 show the examples of test-case simulations where the turbine does, in fact, enter such a flow state. An important motivation for this exploration is the ability of the DRD-BEM model to capture the complex nature of dynamic rotor deformation, with the resulting misalignment of the blade airfoil sections with respect to the incident flow, and its subsequent effects on the aerodynamic forces. This feature gives the DRD-BEM a capacity to detect a wider range of induction factors, improving turbine performance predictions at a wide variety of operational regimes.

A full analysis of the flow across an actuator disk, based on a turbulent Navier-Stokes solution by the $k-ϵ$ realizable model, was conducted in ANSYS-Fluent for operational regimes covering a wide variety of axial induction factors. Results from these numerical experiments for the $CT$ versus a curve are in good agreement with the experimental data set used by Glauert [15], and are also in excellent agreement with the parabolic curve given by the classic streamtube control-volume theory all over the range where it is valid. Figure 5 and Figure 9 show how both Glauert’s data set as well as the present numerical results start departing from the theoretical curve at approximately $a=0.28$ to $a=0.4$. Also, the sudden increase in $CT$ in the experimental data, as seen in Figure 3, coincides with the appearance of wake recirculation in the numerical model (see Figure 7, and the discussion associated with it), which supports the validity of the numerical model used in this work.

A 4th order polynomial regression was proposed to fit the numerical data obtained by the numerical model. This new fitting curve was introduced in the context of the DRD-BEM model to compute $CT$ for axial induction factors a higher than $0.28$. An important observation that can be made on the Glauert [15] correction (and indeed for any other previous correction proposed since) is that it does not include the vortex-ring or the propeller-break states starting at around $a=0.9$. The new polynomial curve proposed here extends the correction all the way upto $a=1.5$ for the first time, covering regimes that were not covered previously. Hence, this study extends the capacities of BEM models in general, and the DRD-BEM model in particular to the analysis of situations beyond $a=0.9$. For the reasons mentioned before, it can be expected to see this operational flow regimes gaining in importance with the turbines of the future. Thus, it can be said that the study presented in this work fills an important gap in the available modeling resources for wind turbine performance and dynamics.