Joukowsky actuator disc momentum theory

Abstract. Actuator disc theory is the basis for most rotor design methods, albeit with many extensions and engineering rules added to make it a well-established method. However, the off-design condition of a very low rotational speed Ω of the disc is still a topic for scientific discussions. Several authors have presented solutions of the associated momentum theory for actuator discs with a constant circulation, the so-called Joukowsky discs, showing the efficiency Cp → ∞ for the tip speed ratio λ → 0. The momentum theory is very sensitive to the choice of the radius δ of the core of the centreline vortex as the pressure and velocity gradients become infinite for δ → 0. Usually the vortex core area is not included in the momentum balance, as it vanishes for δ → 0. However, the pressure in the vortex core behaves as a Delta function and so contributes to the balance, thereby cancelling the singular behaviour. Applying this in the momentum balance results in Cp → 0 for λ → 0, instead of Cp → ∞. The Joukowsky actuator disc theory is confirmed by a very good match with numerically obtained results. At the disc the velocity in the meridian plane is shown to be constant. The Joukowsky calculations give higher Cp values than corresponding solutions for discs with a Goldstein-based wake circulation published in literature.


Introduction
Although the concept of the actuator disc is more than 100 years old, it is still the basis for rotor design codes using the blade element momentum theory developed over these 100 years (see van Kuik et al., 2015).In recent years the behaviour of actuator disc flows with a low rotational speed has been studied by several authors, providing several solutions depending on the type of load that is applied (see, e.g., Sørensen, 2015).Research has focussed on rotors and discs with a Joukowsky distribution (Joukowsky, 1918), having a constant circulation in the wake, or with a Betz distribution (Betz, 1927), yielding a helicoidal wake structure moving with a uniform axial velocity.Goldstein (1929) was the first to find a solution for this wake for lightly loaded propellers (see Okulov et al., 2015, for an overview).Both distributions were assumed to represent the circulation distribution of an ideal rotor.The present paper considers the Joukowsky distribution and compares the results with solutions of the Betz-Goldstein distribution modified for heavily loaded actuator discs reported in Okulov and Sørensen (2008), Okulov (2014), and Wood (2015).
For a Joukowsky actuator disc, the swirl of the wake is induced by a discrete vortex at the wake centre line, leading to an infinite azimuthal velocity and pressure for the radius r → 0. The question of how to model the discrete vortex and how this impacts the momentum balance has been studied by, e.g., de Vries (1979), Sharpe (2004), Xiros and Xiros (2007), Wood (2007), Sørensen and van Kuik (2011), and van Kuik (2016).Apart from the last reference, the reported performance predictions show a remarkable result: in the limit to zero rotational speed, the efficiency of the disc increases to infinity, which is highly non-physical.Within the inviscid flow regime, the analysis in Sørensen and van Kuik (2011) is considered to be exact apart from the choice of the vortex core at the axis of the wake.The centreline vortex is assumed to be a Rankine vortex of which the core diameter is proportional to the wake radius.The analysis of Sørensen and van Kuik (2011) shows that adding a disturbance parameter to the momentum balance removes the non-physical result of infinite efficiency for zero rotational speed, no matter how small this disturbance is.This is an indication that the mo-mentum balance is very sensitive to small deviations in the flow parameters.
A failed attempt to reproduce the results of Sørensen and van Kuik (2011) by the potential flow actuator disc code described in van Kuik and Lignarolo (2016) initiated a reanalysis of the vortex core model and its impact on the momentum theory.In van Kuik (2016) the results of the modified momentum theory were published, showing a confirmation by the potential flow calculations.Still an interpretation question regarding the choice of the vortex core model was left unanswered.The present paper includes the content of this conference paper extended with an answer to the interpretation question and with several detailed results.In Sect. 2 the equations of motion for Joukowsky actuator disc flows are given as well as the properties of the disc load, far wake and vortex core.Herewith, the general mass, momentum and energy balances are derived in Sect. 3 and combined to allow performance predictions of Joukowsky actuator discs.Section 4 describes the numerical approach and its results, which are compared with the momentum theory results in Sect. 5.This section also includes the comparison with the Betz-Goldstein solutions reported in literature.Section 6 presents the conclusions.

The equations for a disc with constant circulation
The flow is governed by the Euler equation: in which ρ is the fluid density (kg m −3 ), f the force density (N m −3 ), p the static pressure (N m −2 ), v the velocity vector (m s −1 ) and H = p + 1 2 ρv × v the total pressure (N m −2 ).A cylindrical reference system (x, r, ϕ) is applied, with the positive x coinciding with the downwind wake axis, and with r and ϕ the radial and azimuthal coordinate (see Fig. 1).For the special case of a disc flow with constant circulation induced by a free potential flow vortex at the axis of the wake with a vortex core having radius δ(x), the azimuthal velocity in the wake is The functions δ(x) and F(r/δ(x)) remain unspecified.Figure 1 shows (half of) the cross section through the stream tube in the meridian plane, with the disc and fully developed wake indicated.The shaded area is the vortex core with an increasing radius towards the far wake due to the flow deceleration.The fully developed far wake is indicated by the index 1.If there is no index, the variables are taken at the position of the actuator disc.The index 0 is used for flow variables in the undisturbed, upstream flow.The disc has radius R and area A, while A 1 is the area of the far wake with radius R 1 .The analysis starts with δ being non-zero, after which the limit of δ → 0 is taken.The only assumption made is that lim δ→0 δ 1 = 0 while δ 1 δ > 1. (3)

The disc load
Only the pressure and the azimuthal velocity will be discontinuous across the infinitely thin disc, so integration of the axial and azimuthal component of Eq. ( 1) gives where F denotes a surface load (N m −2 ), the difference between the down-and upwind side of the disc and e the unit vector.As v ϕ = 0 at the upwind side of the disc, v ϕ = v ϕ .In Eq. ( 5) the Bernoulli equation integrated across the disc thickness has been used: The local power converted by the force field f is f ×v, which has to be equal to the local contribution to the torque, rf ϕ , times the rotational speed .The converted power f × v becomes Integration of Eq. ( 7) across the disc combined with the azimuthal component of Eq. ( 4) gives the general expression This shows that the work done by the force field is expressed in a change in the total pressure or Bernoulli constant H .With Eq. ( 2), for the Joukowsky disc, It follows that outside the core, H is constant, by which Eq. ( 6) shows that any non-uniformity in the pressure jump is due to the creation of swirl across the disc.The swirlpressure jump does not change H , so does not contribute to the conversion of power, by which Eq. ( 6) may be interpreted as p = p converting-H + p conserving-H .The sign conventions are that the rotational speed > 0 and < 0 so H < 0, implying that energy is extracted from the flow.The thrust T is obtained by integration of Eq. ( 6) on the disc area.In dimensionless form the thrust coefficient is C T = T /( 1 2 ρU 2 0 πR 2 ), containing both terms on the righthand side of Eq. ( 6), here denoted as H and ϕ.For δ → 0, where λ is the non-dimensional tip speed ratio λ = R U o .In the same way the power coefficient is defined as C p = P /( 1 2 ρU 3 0 πR 2 ), where P denotes the absorbed power.Evaluation of C p is done in Sect.3.1.

The far wake for r ≥ δ 1
With the conservation of circulation, the Bernoulli Eq. ( 9) in the wake is written as Differentiating with respect to r and combining it with the radial pressure equilibrium in the far wake, it is clear that v x,1 is constant.By this Eq.( 12) can be written as At the wake boundary the pressure has to be undisturbed (p 0 ), so p * = 1 2 ρv 2 ϕ,R 1 and, with Eq. ( 2), This shows that the pressure variation in the far wake is caused only by the swirl.By merging Eq. ( 15) with Eqs. ( 9) and ( 12), the second term on the right-hand side appears as a loss in H due to swirl: This is consistent with the optimization of rotors according to Glauert's theory which involves minimization of the swirl (see, e.g., Sørensen, 2015).

The vortex core
The momentum theory results are very sensitive to the choice of δ and δ 1 because of the logarithmic singularity resulting from the integration of the pressure due to the azimuthal velocity: at the disc −ρπ R δ v 2 ϕ rdr = −ρ 2 4π ln R δ , and similarly in the far wake −ρ 2 4π ln R 1 δ 1 .Previous solutions have dealt with the singularity in different ways.Sørensen and van Kuik (2011) have adopted δ δ 1 = R R 1 , assuming that the vortex core grows with the stream-tube radius.This removes the singularity as −ρ Kuik (2016) assumes δ = δ 1 leading to the power coefficient C p → 0 for λ → 0. However, as discussed in van Kuik (2016), both core models do not comply with the inviscid flow equations, so the impact of the vortex core model to the momentum balance merits an additional investigation.
Both analyses used the vortex core boundary as a lower limit in the integration of momentum and energy on the control volume used in momentum theory.This implies that the vortex core is excluded, motivated by its vanishing dimension in the limit δ, δ 1 → 0. Here, it will be included, while the same limit is taken.
With δ(x) denoting the local core radius, with δ ≤ δ(x) ≤ δ 1 , the Bernoulli Eq. ( 9) in the vortex core region becomes www.wind-energ-sci.net/2/307/2017/Wind Energ.Sci., 2, 307-316, 2017 where v s is the velocity in the meridian plane.As v s , U 0 and the last term on the right-hand side remain finite for δ(x) → 0, they do not contribute in this limit to the axial momentum balance drawn on the core volume.Consequently, this balance reduces to a balance of pressures acting on the control volume boundaries, integrated as a load in x direction: where the path of integration of the third integral is the core boundary δ(x) with 0 ≤ x ≤ x 1 .This integral is evaluated with Eq. ( 17), of which the finite terms do not contribute as the area of the core boundary projected in x direction is π (δ 2 1 − δ 2 ), which vanishes for δ(x) → 0. Only the pressure terms originating from v ϕ contribute as these behave like a Delta function: The combination of Eqs. ( 18) and ( 19) gives irrespective of the choice of core model δ(x), F (r/δ(x)) .This result will be used in Sect.3.1, where the momentum balance for the entire stream tube is studied.Unless specifically indicated all equations in the forthcoming sections are derived for the flow region outside of the vortex core.

The momentum, mass and energy balance
The momentum equation drawn on the stream tube as control volume (see Fig. 1) is written as where T is the thrust (N), being the integrated pressure jump across the disc.The boundaries of the momentum balance volume are the stream-tube boundary and those of the cross sections A 0 and A 1 , far up-and downstream.As discussed in many references, amongst others in van Kuik and Lignarolo (2016), the pressure at the stream-tube boundary does not contribute to the momentum balance and so is not included in Eq. ( 21).
Figure 1 shows the pressure distributions appearing on the left-hand side of Eq. ( 21), including the thrust: 1. constant pressure jump across the disc giving the jump in Bernoulli parameter H according to the first term on the right-hand side of Eq. ( 6).
2. pressure distribution due to the jump in v ϕ for r ≥ δ according to the second term on the right-hand side of Eq. ( 6).This term conserves H .
3. the same pressure distribution in the far wake due to the v ϕ distribution for r ≥ δ 1 according to the first term on the right-hand side of Eq. ( 15), conserving H .
4. constant pressure to achieve p 1 − p 0 = 0 according to the second term on the right-hand side of Eq. ( 15) or Eq. ( 16).
5. the contribution by the vortex core cross sections (Eq.20).
When all contributions are expressed in by Eqs. ( 2) and ( 9), integrated, subjected to limδ → 0, substituted in Eq. ( 21) and divided by the disc surface π R 2 , the result is where the terms on the left-hand side have been named in accordance with Fig. 1.The term between square brackets simplifies to ln(R/R 1 ).In other words: only the wake expansion area c >R contributes to this term.The mass balance is with the bar above v x indicating that it is the average value.
The energy balance follows from Eq. ( 16): Mixing Eqs. ( 22) and ( 23) simplifies the momentum balance, yielding Wind Energ.Sci., 2, 307-316, 2017 www.wind-energ-sci.net/2/307/2017/The non-dimensional vortex q = − 2πRU o is introduced.As < 0 q > 0. Furthermore, from here on v x and v x,1 indicate the dimensionless value v x U 0 and v x,1 U 0 .Herewith Eq. ( 9) becomes and the momentum balance becomes and the energy balance becomes An analytical solution of Eqs. ( 27) and ( 28) is not found.An implicit expression of v x,1 in the independent variables λ, q is obtained by writing Eq. ( 28) as an expression for v x with the help of Eq. ( 23) and substituting this in Eq. ( 27): This can be solved numerically for v x,1 .The wake expansion follows from Eq. ( 28) and the velocity at the disc from Eq. ( 27).The power coefficient follows by integration of Eq. ( 7) on the disc area: By mixing Eqs. ( 27) and ( 28), the velocity at the disc can be written as As (1 + ln (R/R 1 ) 2 ) < (R/R 1 ) 2 for R < R 1 the ratio is < 1. Consequently v x < 0.5 v x,1 + 1 .The ratio in Eq. ( 31) is the ratio between the left-hand sides of the momentum balance Eq. ( 27) and energy balance Eq. ( 28) or, in other words, between the total load exerted on the flow in the streamtube control volume and the non-conservative load which is the load performing work.By this, Eq. ( 31) is equivalent to Eq. ( 6) of van Kuik and Lignarolo (2016) where the distinction between conservative and non-conservative loads is used www.wind-energ-sci.net/2/307/2017/Wind Energ.Sci., 2, 307-316, 2017 to explain the results of the momentum theory applied to an annulus of the stream tube.This analysis, without swirl, is done again with swirl in Sect.4.3.

Limit values of the Joukowsky momentum theory for
λ → 0, λ → ∞ and for maximum C p For large values of λ, the wake angular momentum should go to 0, and the momentum theory should become the one-dimensional theory yielding the well-known Betz-Joukowsky maximum value for C p .According to Eq. ( 26) q is inversely proportional to λ for constant H or λq.In the balances Eqs. ( 27) and ( 28), the q 2 terms vanish for λ → ∞, with which the momentum theory without wake swirl is indeed recovered.
For the limit λ → 0, flow states with λq being constant are studied.The energy balance Eq. ( 28) shows that the highest value for q 2 (R/R 1 ) 2 is obtained for v x,1 = 0: With v x,1 = 0, the right-hand side of the momentum balance is 0 as is clear from Eq. ( 22), by which it becomes 2λq Elimination of q 2 from Eqs. ( 32) and ( 33) gives the wake expansion for the highest q−lowest λ: As an example, 2λq = 8/9 results in R 1 R = 2.77, q = 0.924 from Eq. ( 32) and λ = 0.48.Both v x and v x,1 are 0, but the ratio of v x v x,1 → 7.69.This flow state is characterized by a full blockage by the disc, creating a wake with azimuthal flow only, so there is no change in axial momentum.The associated pressure distributions in the wake and at the disc balance each other.A lower value of λ is not possible for this value of λq.For λq = 0 with λ = 0, Eq. ( 34) gives In the wake only the azimuthal velocity is non-zero, reaching qR R 1 = 1 at the far wake boundary r = R 1 .The wake expansion is close to the experimental value ≈ 1.6 of the wake expansion behind a solid disc reported in Craze (1977).
C p,max (λ) is obtained by optimizing the solutions for fixed λ varying q.

Flow and pressure field
The computer code described in van Kuik and Lignarolo (2016) has been adapted to include wakes with swirl.Axial and radial velocities are calculated by summation of the induction by each of the several thousand vortex rings which constitute the wake boundary.The azimuthal velocities are calculated from Eq. ( 2).The shape and strength of the vortex rings are adapted in the convergence scheme to satisfy the two boundary conditions: zero pressure jump across the wake boundary, and zero cross flow.The first boundary condition p wake-boundary = 0 is expressed in |v| and input parameter H : ( 1 2 ρ|v| 2 ) − H = 0.In van Kuik and Lignarolo (2016), v only had an axial and radial component, but now the azimuthal component also enters this boundary condition.The strength of the vortex at the axis follows from Eq. ( 26) expressed in H and the second input parameter λ : q = − H /(ρU 2 0 λ).Apart from these changes, the code and the numerical parameters are unmodified.The results satisfy the same accuracy requirements as described in van Kuik and Lignarolo (2016).Figures 2 and 3 show the streamlines, expressed in the stream-function and isobars of the disc flow with H /( 1 2 ρU 2 0 ) = −0.8888and λ = 0.731 and 1.018.The isobars in the wake show the pressure gradient due to the swirl.

Constant meridian velocity at the disc
As shown in Figs. 2 and 3 the pressure at the upstream side of the disc is constant, which implies, by the Bernoulli equation, that the absolute velocity |v| upstream of the disc is constant.Figure 4 shows the values of the axial, radial and azimuthal velocity component at the disc as well as the absolute value |v| meridian = v x 2 + v r 2 .Similar distributions of the axial velocity have been calculated by several others, e.g.Madsen (1996), Crawford (2006), Madsen et al. (2007), Mikkelsen et al. (2009) and Madsen et al. (2010).The fact that |v| meridian is constant has been reported by van Kuik and Lignarolo (2016) for actuator disc flows without swirl.The explanation given in van Kuik and Lignarolo (2016) is now extended to include discs with swirl.
The radial component of Eq. ( 1) just upstream of the disc is with s being the coordinate along the streamline and r the radial coordinate.The pressure does not depend on r when it is shown that the radial velocity reaches a maximum at the disc when following a streamline.Along any streamline passing the disc, v r increases when the position of observation s 0 travels from far upstream to the disc s disc , due to the decreasing distance to the vorticity γ in the wake boundary, so ∂v r /∂s > 0. Following the streamline in the wake, so with s 0 > s disc , two regions can be distinguished: the vorticity between s disc and s 0 induces a negative v r and so contributes to ∂v r /∂s < 0, while the induction by the vorticity at s > s 0 is approximately constant as the wake downstream of s 0 remains semi-infinite, with γ depending only slightly on s for s > s 0 .The result is that ∂v r /∂s = 0 at the disc position, so from Eq. ( 35) the pressure upstream of the disc is constant and, from the Bernoulli, equation |v| meridian is constant; QED.
Figure 5 shows the calculated strength of the vortex sheet for the load case of Fig. 3, confirming the reasoning.With |γ | having a maximum at its leading edge, the non-uniformity of γ contributes to a negative induction of v r at streamline positions s 0 > s disc .It should be noted that the distribution in Fig. 5 does not show the irregular behaviour at the leading edge as shown in Fig. 9 of van Kuik and Lignarolo (2016).The explanation is that the distance between the first vortex rings in this previous paper is smaller than the radius of the vortex ring core, leading to this irregularity.Calculations with a smaller core size, not yet reported, have removed this irregularity, thereby not having any observable impact on the flow pattern and integrated numerical results.In the present calculation the distance between rings is always larger than the radius of the core of the vortex ring.Now that the pressure at the upstream side of the disc is known to be constant, the radial derivative of Eq. ( 6) becomes with p − being the pressure at the downstream side of the disc.This is the radial equilibrium expression Eq. ( 13) for the flow in the wake.Apparently the radial distribution of p is only linked to v ϕ , not to the other velocity components.

Momentum balance per annulus
In van Kuik and Lignarolo (2016), the non-uniformity of the axial velocity at the disc is explained by applying the momentum theory to annuli instead of the entire stream tube.An annulus is defined as the volume of flow between the stream-tube values n and n−1 .With 1 < n ≤ 10 and n − n−1 = 0.1 stream tube , Figs. 2 and 3 shows 10 annuli as the volume between the plotted streamlines passing the disc.For the flow case shown in Fig. 3, the momentum balance per annulus is evaluated as follows.
The balance for the entire stream tube is defined by Eq. ( 21).In this equation it is implicitly assumed that the pressure acting at the stream-tube boundary does not contribute, as discussed in many publications, amongst which  are van Kuik and Lignarolo (2016).This is confirmed by the calculations: the force in x direction resulting from the pressure integrated along the wake boundary for −25 < x/R < 25 is 0.2 % of the non-conservative disc load H π R 2 .However, when applying Eq. ( 21) per annulus, the pressure acting at the boundaries of the annulus has to be added, so the second term in Eq. ( 21) is expanded to become S (p − p 0 ) 2π rdr, where S is the contour in the meridian plane of the annulus control volume: The pressure integral is calculated with x/R = ±25R as upand downstream limits.The momentum balance Eq. ( 37) yields the wake velocity v x,1,annulus which, combined with the mass conservation  2014) and Wood (2015) for rotors with an infinite number of blades.
results in v x,annulus at the disc.Figure 6 shows the distribution of v x resulting from the flow field calculation, the associated average value per annulus, and the value resulting from Eqs. ( 37) and (38).Apart from the outer annulus, the calculated average per annulus coincides with the momentum balance average values.This confirms results as published by Sørensen and Mikkelsen (2001) and van Kuik and Lignarolo (2016), both for disc flows without swirl: the annuli cannot be assumed to be independent, as the pressure field contributes to the axial momentum exchange leading to the non-uniform distribution of v x .

Comparison of calculated and momentum-theory-predicted performance
Figure 7 shows the comparison of the Joukowsky momentum theory and the potential flow results.The correspondence between both is excellent.A comparison with the C p,max − λ curve for discs having a modified Betz-Goldstein distribution of the circulation is shown in Fig. 8.As shown by Okulov and Sørensen (2008) and Okulov (2014) the original Betz-Goldstein solution for a rotor with a finite number of blades resulted in C p,max = 1, as the pitch of the helicoidal wake was based on the undisturbed velocity.With the pitch based on the velocity in the rotor plane, Okulov (2014) showed that C p,max reaches the well-known Betz-Joukowsky maximum 16/27 for high λ.The C p,max −λ curve of this corrected solution expanded to a rotor with an infinite number of blades is shown in Fig. 3 of Okulov (2014).An alternative solution is published in Wood (2015), where the Goldstein formulation is adapted to allow for non-zero torque when λ → 0. A comparison of the Joukowsky maximum C p curve and corresponding Betz-Goldstein-Okulov/Wood curves is given in Fig. 8.The Joukowsky distribution gives higher C p,max than the Betz-Goldstein-based distributions, with the difference vanishing for higher λ.This is confirmed by Okulov and Sørensen (2010), where rotors with a finite number of blades having a Joukowsky and Betz-Goldsteinbased distribution have been compared.

Conclusions
An actuator disc momentum theory including wake swirl has been developed resulting in the physically plausible result that C p → 0 in the limit λ → 0. For high λ the theory reproduces the results of the classical momentum theory without swirl.
The novelty in the momentum theory is the removal from the momentum balance of the singular behaviour of the pressure near the wake centreline vortex, giving rise to nonphysical results in several previously published methods.This removal is done by including the vortex core in the momentum balance.
The momentum theory results are very accurately confirmed by potential flow field calculations.
At the actuator disc the velocity in the meridian plane is constant.
The Joukowsky momentum theory results are higher than the equivalent results for rotors with an infinite number of blades optimized for Betz-Goldstein solutions.

Figure 1 .
Figure 1.Pressure distributions acting in the momentum balance.The arrows give the direction of the pressure fields acting on the flow.The meaning of a, b, c, d and e is given in Sect.3.1.

Figure 6 .Figure 7 .
Figure 6.v x at the disc for the load case H /( 1 2 ρU 2 0 ) = −0.8888and λ = 1.018: the calculated distribution, the calculated average per annulus and the result from the momentum balance per annulus.The two annuli lines coincide except in the outboard annulus.

Figure 8 .
Figure8.The Joukowsky actuator disc C p,max compared with the Betz-Goldstein C p,max solutions ofOkulov (2014) andWood (2015) for rotors with an infinite number of blades.