The actuator disc theory for wind turbines including losses and non-uniform velocity distribution

In a paper published in 1989 the actuator disc model proposed by Betz to represent the flow across a wind turbine has been set right by taking losses into account. In the present paper this model is improved by introducing global correction factors at the entry and the exit of the machine. The resulting extent of the actuator disc concept shows that management of the velocity field may lead to a satisfactory adjustment between the power output and the velocity deficit range. A significant shift in the characteristics of the machine, associated to a slight enlargement of the wake, is made possible in the frame of the new model. Comparison reveals a better fit with experimental data.


Introduction
The concept of actuator disc, extended to horizontal wind turbines by Betz [1,2], is a by-product of the propeller theory.The criterion, usually referred to Betz, has been discovered independently, nearly at the same time, by Betz, Lanchester and Joukowsky [3].The Betz approach (Appendix A), gives rise to some inconsistent results originating from the use of the Bernoulli equation on both sides of the machine, an unfounded assumption which does not take physics into consideration.The error has been repeated in the literature day after day, with the result that the Betz' model still remains a reference for the first step study in the field of wind turbines, even in recent articles.However, this model has been modified 25 years ago by including losses [4].Afterwards, following a statement by Theodorsen [5], Prado has proposed a middle term adjustment between Betz' model and our own, by superimposing a fictitious source, located for the best inside the wake, in order to get a better fit with experimental results [6].In the present contribution the simplified flow field from [4] is improved by taking the non uniformity of the velocity profiles into account.Contrary to what happens generally in the frame of quasi one dimensional flow, here non uniformity has not trivial consequences.The usual corrections act through a global factor, the range of variation of which has the same order of magnitude that the mean velocity deficit across the disc.This fact is essential in the new model as it has significant effects on the losses and on the power output delivered by the machine.
The actuator disc is a vertical porous flat plate of area S facing an incompressible fluid of density ρ which moves at the uniform speed U (Fig. 1).The retarded flow induced by the disc is characterized by the volumetric discharge q and the head H absorbed by the disc, the extracted power being qH.Making use of the primary quantities ρ, U and S, the parameters a and C p introduced by dimensional analysis are defined by q = (1 − a)SU and qH = (ρU 3 S/2)C p : a is the velocity interference factor and C p the power output coefficient.
Comparison with the working of a turbine looks a priori complicated, because of the complementary data relating to the rotor, i.e. the angular speed of rotation Ω, the radius of the rotor R, the number of blades and the geometric characteristics of the blades.For a family of turbines having the same number of geometrically similar blades, the sole supplementary dimensionless parameter is the tip speed factor ΩR/U .Rotation has little influence in front of a turbine, contrary to the exit where the flow cannot be considered as almost one dimensional as the tip speed is generally higher than U .So, it seems that comparison will fail at first because no counterpart can be associated to the tip speed.Moreover, comparison is to be made between an averaged in time periodical flow and a steady one.Finally, depending on fewer parameters, the actuator disc can bring no more that a partial representation of a turbine.
Nowadays the actuator disc is still used for studying meandering and turbulence decay in turbine wakes, the objective being to maximize the performance of wind farms [7].As for wind turbines, they are investigated by approaches deriving from the lifting line theory.In the model put forward by Glauert [8] the Prandtl horseshoe vortices located along the span of a wing are replaced by a rotating vortical system attached to the blades composing the rotor.Every blade behaves as a rotating wing along which the velocity and angle of attack vary versus the distance at the axis of rotation.As the flow is considered free of losses, referring to the Betz limit seems pertinent.Really, the main reason of the Betz criterion success lies in its relationship to theoretical approaches within the frame of inviscid fluid.The Glauert model has accompanied the progresses in the wing theory throughout the last decades [9][10][11], during which the Betz limit still remained an out of reach unquestionable criterion when efficiency was raised.Later on, computational fluid dynamics modeling has been used to address wind turbines similarly to all branches of fluid mechanics, especially during the periods when fuel crises arose.A recent general review devoted to wind turbines [11] contains more than 130 references!Our purpose is to put forward a reformed theory of the actuator disc in case of retarded flow.The present paper is directly issued from a previous one [4].The approach is the same.The actuator disc is encased in a coaxial conduit (Fig. 1) which is useful only temporary in order to simplify the calculations.The disc is loaded with a sudden jump in pressure.How this purpose is fulfilled lies beyond the scope of the analysis.Generally the disc consists of a wire netting or a metallic grid [7].When passing through the disc the viscous stresses developed in the fluid are negligible with regard to the momentum variations issued from sudden changes in direction and enlargements.Inertia effects prevail provided that the Reynolds number characterizing the mesh is high, a circumstance which is always met in practice.The lattice acts similarly to honey combs in wind tunnels, by producing turbulence in the wake.Stating that the velocity distribution is almost one-dimensional is consistent with constant values of the pressure on each face.However, the fact that the velocities are slightly different on both sides implies that some changes occur during crossing.The exact location of the surfaces S − and S + marking the entering and the leaving velocity distributions are not known, except that these boundaries are situated very close to the area S occupied by the disk.A fictitious thickness must be introduced: that is the distance between the flow separations in front of the lattice and the reattachments behind it.In fact, this thickness is the actual one when the disc is regarded as an aerodynamic slice shaped device, able to generate the entry and exit values of the pressure.

Governing equations revisited by including velocity correction factors
All symbols are indicated in Figure 1.The duct sectional area is Σ.The undisturbed velocity and pressure are U and P .The sectional area of the stream tube feeding the machine, far upstream, is s.The distribution in  b, c, and d are positive as air is slowed by the turbine, with 0 < a ≤ b < 1, involving σ ≥ S. The pressures at the entry and at the exit are P and P .Along S 1 the mean pressure is P 1 and it is P 2 along S 2 and σ.Thus, there is a jump in pressure along the rim, on both sides of the disc.The case of a disc in an infinite atmosphere will be obtained by removing Σ at infinity, which makes c and d tending to zero.
All areas interfering in conservation laws are connected by: The discharge along the duct is constant: Therefore: ( Conservation of mass in the central stream tube and outside is expressed by: ( The rates of enlargement are S/s = 1/(1 − a) upstream of the disc and σ/S downstream.The usual momentum and kinetic energy correction factors, β and α, resulting from the non uniformity of the velocity distribution, are introduced in the Appendix B. At the disc β and α are defined by 3 dA, v being the axial velocity.In fact, previous integrals are taken on S − and S + , depending on the considered side of the disc.A new global correction factor λ is defined in the Appendix B, λ = (1 + α − 2β) 1/2 , the quantity between parentheses being proven positive.If the velocity is not far from uniform, α−1, β −1 and λ 2 are very close to zero, but λ cannot be considered as small in comparison with the velocity interference factor and the power output coefficient, as it will be seen in Section 4. As a consequence, terms containing λ 2 will be neglected as far as possible.
The correction factors are β , α and λ in front of the disc and β , α and λ behind it.In front of the disc the fluid velocity varies gradually against the distance from the axis, while at the exit the velocity profile is curved when the edge of the wake is approached, especially close to the rim, in relation with the velocity defect occupying the central part of the wake.The downstream distribution in velocity being much more far from uniform that upstream, λ is expected to be higher than λ .Over the area σ where the flow in the wake recovers a fully developed regime the velocity is considered as almost uniform.
As for the areas S 1 and S 2 , it is taken into account in advance that β and α tend to 1 when the casing is removed at infinity (Appendix B).
Let us recall the arguments advanced in [4] on behalf of our approach.The eddies resulting from the instability of the shear layer formed at the border of the retarded fluid, upstream of the disc, are located rather on the side where the flow is slowed, that is inside the slipstream.Consequently, the Bernoulli equation is not valid inside the progressive enlargements in front of the disc and in the wake.On the contrary it can be applied in the contracted external stream: Note that in case of a propeller the Bernoulli equation is valid on both sides of the machine as the flow in the slipstream is accelerated.
The momentum equation is written for the closed cylindrical domains indicated by dotted lines, upstream and downstream of the disc (Fig. 1).We have upstream: then, thanks to (6), Finally, Equation (4) leads to: We proceed similarly downstream: Thanks to (6), the contribution of all pressure terms is Consequently, taking (2) and (3) into account: The losses which exist at every stage of the flow inside the diverging slipstream, upstream J and downstream J , are: No losses have been taken into consideration by Betz, because he regarded the wind turbine as the reciprocal of a propeller, being probably influenced by the propeller model formerly introduced by Froude.

Effect of velocity corrections on flow characteristics
Now the enclosing duct is removed to infinity, so that c and d tend to zero.Then, according to Equations ( 6) and ( 7), we obtain As for P , the last terms in the right hand side of formula (8) tend to zero, with the result: Therefore, the pressure drop across the disc is: According to ( 9) and (12), the upstream losses are given by A similar calculation achieved downstream, by starting from ( 10) and (13), leads to: The drop in head across the disc, H, is: Energy is extracted from the machine as the last term is negligible.
Let (ρU 3 S/2)C L be the power delivered to the fluid passing in the slipstream and (ρU 3 S/2)C P the power output, that is the power extracted from the machine.With the volumetric discharge through the disc, q = (1− a)U S, these quantities are defined as usually by: Thanks to the formula (17), the power output coefficient is: The drag D is obtained from the momentum equation, without forgetting the contribution of the fluid passing through the disc.The result is: D = (P − P )S + ρU 2 S(1 − a) 2 (β − β ).Then, according to (14), When the disc is little permeable, a tending to 1, the drag becomes zero.This result comes from the fact that viscous stresses which are dominant at low Reynolds number are not taken into account, as mentioned above.The drag coefficient C D , defined by D = (ρU 2 S/2)C D , is

Consequences deriving from the existence of losses
A first approximation solution is obtained when the additional terms containing the squares of the global factors λ and λ are neglected in formulas (15) and ( 16).That leads to 2J /(ρU 2 ) = a 2 , 2J /(ρU 2 ) = −(b − a) 2 .These results, which have been published in [4], show that losses exist only upstream of the disc and that b is necessarily equal to a, making the wake cylindrical, σ = S.In this case, the upstream losses, originating in the slowing down of the fluid arriving from infinity, are exactly the ones occurring in a sudden enlargement of amplitude S/s with velocity correction factors equal to 1.
The subscript zero being reserved to this solution, it has been obtained that the velocity deficit cannot exceed 1/3 in order to make the dissipation minimum for given power output [4].The results are: The highest output C P o = 8/27 is reached at a = 1/3.This output is twice lower than that obtained by Betz.As modern wind turbines reach values of the output coefficient not far from 0.5, the model is not satisfactory.Moreover, as noticed by Prado [6], since vortices continue to be produced downstream of the disc, the wake should thicken instead of remaining cylindrical.That is exactly what happens when the full expressions of the losses are considered.
Coming back to the formulas ( 15) and ( 16), the conditions J ≥ 0 and J ≥ 0 involve The rates of spreading are subjected to S/s ≥ 1 + λ , and σ/S ≤ 1/(1 − λ ).According to the last condition the wake may spread, but little.With λ = 0.1, for example, σ/S ≤ 1.11 is necessary, so the wake cannot be larger in diameter than 1.05 times the disc.We are far from the Betz limit.Both conditions show that the assumption of one dimensional flow is founded as long as λ and λ are sufficiently low.Downstream, in order to keep J ≥ 0 we may set with the balance coefficient δ constant, 0 ≤ δ ≤ 1.That makes σ/S = 1/(1 − δλ ) and Consequently, J is negligible, which means that the Bernoulli theorem is approximately valid downstream.However, this is accompanied by an unavoidable loss upstream of the disc which may be approximated by 2J /(ρU 2 ) = a 2 , according to (15), except when the interference factor a takes values close to λ /(1 + λ ).For this particular value J cancels too, so the fluid behaves as inviscid everywhere.When J is negligible, but not J , the fluid may be considered as inviscid only downstream, with the result for a turbine associate to the actuator disc that the Glauert theory may be implemented, but with a frontal head reduced by a 2 ρU 2 /2.All dimensionless parameters C L , C P and C D can be expressed in terms of λ instead of b.Consider the power output coefficient C P .The jump in λ 2 being neglected in formula (20), C P becomes approximately: This formula shows that C P does not depend directly on λ , but on the sole parameter λ , except when a is close to its lowest value λ /(1 + λ ).Inspection of the shift produced by the corrections is made now.The following comments are illustrated in Figure 2. As stated above, no substantial modification in all previous results is issued from any jump in λ 2 , because λ 2 is negligible.On the contrary, taking into account that to the first approximation a must be lower than 1/3, variations in λ and λ cannot be neglected as they may be of the same order that a.This observation has important consequences as it will be seen below.
For given λ , all curves representing C P in terms of a, formula ( 26 has been selected, the power output is obtained by varying a beyond λ /(1 + λ ), according to (23).Then C P grows and becomes maximum at a = a * , the first value of the velocity deficit making zero the derivative dC P /da.
The derivative is zero also at the root a = 1 which is eliminated since it does not correspond to a maximum.The right root is given by and the maximum value of C P is: So, a * is lower than 1/3 while C P * is higher than 8/27.Note that a * is a decreasing function of λ .The highest possible value of λ , corresponding to a * = 0 is 1/3; it leads to C P * = 2/3.By account of continuity with the first approximation (subscript zero), the solution is unique when a is lower than the particular value a * making C P maximum.This is related to the fact that the essential of losses being proportional a 2 , when C P is given the appropriate value of a in ( 26) is the lower one because it involves the lowest losses.Consequently, the range of variation of a is restricted to The lowest value of this quantity must be positive, involving the condition: As λ and λ are positive or zero, we find again the condition λ ≤ 1/3, which ensures that a * is positive.
It must be emphasized that the increase in C P may become important, particularly at low values of a, the two contributions in the right hand side of formula (26) becoming of the same order of magnitude.Indeed, the two terms are equal if (1 − a)δλ = a.Suppose a = a * , that is a * = 1/6 according to the formula (27), and therefore δλ = 1/5.Then, the output remains almost constant.At the lowest limit, a = λ /(1 + λ ), the loss J is zero and we have C P o = 2λ /(1 + λ ) 3 and C P = 2(λ + δλ )/(1 + λ ) 3 .When δ = 1, J is zero too.So, as already signaled, there exists a particular solution without dissipation, which is far from the Betz one.This solution provides an important profit on account of C P /C P o = 1+λ /λ .It requires a searching adjustment between a, b, λ and λ .The principal condition to be met is a = λ /(1 + λ ), in order to cancel the major part of the losses, occurring upstream.Previous results show that rational values of λ are able to make the output not only higher but also approximately independent of the interference factor.Per contra, the range of variation of the deficit is reduced.Theoretical results giving rise to an increase in output beyond the limit prescribed by the first approximation model (subscript o) trend to a better fit with experimental data.In this connection, according to the literature, C P does not exceed 0.45 for horizontal wind turbines [9,10].
For the second example, keep λ = 0.08, δ = 1 and choose λ = 0.20 in order to inspect the optimum case raised above: δλ = 1/5, a * = 1/6.The condition (29) is still fulfilled and we obtain C P * = 0.463 from the formula (28), a value which is close to the maximum of C P reached experimentally.
To summarize, two facts support our approach.First, the maximum values of C p delivered by the model and those measured with real turbines are found similar thanks to the presence of λ in the expression (26) of the power output.Second, according to calculations performed for turbines in the frame of Glauert theory, the range of the velocity interference factor a delivering the highest C p lies between 1/4 and 1/3 [11], an interval which is not far from that observed in our model (Fig. 2).Do previous results mean that the global correction factor λ is an appropriate substitute of rotary effects, since the obtained efficiencies are close to the ones measured with actual turbines?Do rotary effects can be expressed only through combined influences of a and λ ?These questions are of importance as, based on the simple foundations of our approach, no direct comparison including rotary effects is possible.However that may be, the influence of velocity correction factors is not so illusory as feared at the beginning of the paper!

Conclusion
Supplementary information about the actuator disc concept has been obtained in the frame of a new approach.Surprisingly, more than 60 years have been necessary to replace the Betz' model by a rational one, and then 25 years more to complement the new model.Thanks to the introduction of global correction factors, some valuable properties are discovered about the losses, the power output and the growth of the wake, making us able to express the power output in terms of the velocity interference factor and the global velocity correction factors, The key of the approach is based upon the following hierarchy in order of magnitude: the small coefficients (α−1), (β − 1) and λ 2 having the same order, the global correction factors λ and λ may be equivalent to the interference factor, which must be lower than 1/3.Therefore, the global correction factors involve a substantial increase of the output, reaching the performances of actual turbines.Moreover, flows without losses are possible in case of special arrangement of the velocity correction.Indeed, mastering the corrections seems a right method to increase the output level and to make it independent of the velocity deficit, as in the second example of Section 4. It seems worth to push investigation somewhat further in order to inspect how the correction factors depend on the distribution of velocity.What is the efficient way to arrange the velocity distribution, at given flow rate, in order to obtain prescribed values of the correction factors?Why not experiment with grid holes different in size, shape and depth, according to their location around the center of the disc.The final aim should be to detect all correlations existing between the non-uniformity in velocity distribution and the way how a rotor operates, a challenge which is not easy.
If the velocity v is referred to v m , y = v/v m − 1, we have A ydA = 0 and the following results are easily obtained: A(β − 1) = A y 2 dA, A(α − 3β + 2) = A y 3 dA.When y is small β − 1, α − 1 and λ 2 are very close to zero.In the case of an unlimited area, y necessarily behaves at infinity in such a way that A ydA tends to zero, which requires y equivalent to 1/r n with n > 2, when the distance r at the axis tends to infinity.