Some effects of flow expansion on the aerodynamics of horizontal-axis wind turbines

The flow upwind of an energy-extracting horizontal-axis wind turbine expands as it approaches the rotor, and the expansion continues in the vorticity-bearing wake behind the rotor. The upwind expansion has long been known to influence the axial momentum equation through the axial component of the pressure, although the extent of the influence has not been quantified. Starting with the impulse analysis of Limacher and Wood (2020), but making no further use of impulse techniques, we derive its exact expression when the rotor is a circumferentially uniform disc. This expression, which depends on the radial velocity and the axial induction factor, is added to the thrust equation containing the pressure on the back of the disc. Removing the pressure to obtain a practically useful equation shows the axial induction in the far wake is twice the value at the rotor only at high tip speed ratio and only if the relationship between vortex pitch and axial induction in non-expanding flow carries over to the expanding case. At high tip speed ratio, we assume that the expanding wake approaches the Joukowsky model of a hub vortex on the axis of rotation and tip vortices originating from each blade. The additional assumption that the helical tip vortices have constant pitch allows a semi-analytic treatment of their effect on the rotor flow. Expansion modifies the relation between the pitch and induced axial velocity so that the far-wake area and induction are significantly less than twice the values at the rotor. There is a moderate decrease – about 6 % – in the power production, and a similar size error occurs in the familiar axial momentum equation involving the axial velocity.

1 Introduction

Conservation of axial and angular momentum are fundamental principles for wind turbine analysis. They are applied using control volumes (CVs) such as those in Fig. 1 or, more commonly, to a CV coinciding with a mean streamtube and extending into the far wake, the hypothetical region of no further wake development. For blade-element momentum theory, the CVs become expanding annular streamtubes intersecting the elements. The change in axial or angular momentum of the flow determines the net thrust or torque, respectively, acting on the rotor or blade elements (e.g. Burton et al., 2011; Hansen, 2015; Sørensen, 2016). Angular momentum is easier to analyze because in most cases it is generated only at the blades.

Figure 1 Control volumes (CVs) to be used in the present analysis. In both variants, the upstream face extends in z to −∞, where the velocity is the wind speed, and R_CV≫R. The downstream control surface is just downstream and just upstream of the rotor plane in CV₁ and CV₂, respectively, and the corresponding downstream control surfaces (CS) are labelled S_D and S_U. Taken from LW.

When a turbine extracts kinetic energy from the wind, the flow must expand both upwind and downwind of the rotor. As noted on p. 185 of Glauert (1935) and by Goorjian (1972), the axial momentum equation may receive contributions from the pressure in the expanding flow upwind of the rotor. The expansion causes the pressure force to have an axial component which alters the rotor thrust by an amount equal to the momentum flux external to the rotor. This is because the pressure forces acting on the cylindrical control surfaces at radius R_CV in Fig. 1 are entirely radial. Although the role of pressure has been recognized for a long time and is discussed by Sørensen (2016) and van Kuik (2018) among others, a satisfactory analysis of it is lacking. The first main result of the present analysis is a closed-form expression for the pressure force for a circumferentially uniform rotor.

Limacher and Wood (2020) (hereinafter “LW”) investigated steady wind turbine thrust, T, using an impulse analysis, whereby the pressure in the axial momentum equation for any CV is replaced by terms that include vorticity fluxes across the CV boundaries. This removal of pressure is achieved by the substitution of various integral identities into a standard momentum-based control volume analysis, as demonstrated by Noca (1997). We will use what we call the “impulse perspective” as explained below but not impulse techniques in this paper; the interested reader is referred to LW for a short history and more details. LW showed that by approximating a rotor as an actuator disc, T is given exactly by integration over the face S_D of CV₁, situated just downwind of the rotor on the left of Fig. 1:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{T}{\mathit{\rho }}}=\underset{S_{\mathrm{D}}}{\int }\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{w}^{\mathrm{2}}+\mathit{\lambda }wx\right)\mathrm{d}S,\end{array}$$

where ρ is the air density, and w is the circumferential velocity (in the direction of θ in Fig. 1) on S_D; in LW, w denoted the circumferential velocity at the rotor plane, which was assumed to be one-half that on S_D. This assumption is also used in the present analysis. λ is the tip speed ratio (λ>0 for clockwise rotation, as viewed from the positive z axis), and x is the radius normalized by the tip radius so that x≤1 for the rotor. The downwind face of the second CV in the figure, S_U, is just upwind of the rotor. The term “exact” will be used throughout this paper to indicate that no assumptions beyond those listed below have been invoked. Taking the “wake” to be the flow that has passed through the rotor, which rotates with a constant angular velocity, these assumptions are as follows:

1. the flow upwind of the rotor and outside the wake is inviscid, steady, and spatially uniform;

2. the total energy of the wake is reduced instantaneously at the rotor, after which it is conserved;

3. viscous and/or Reynolds stresses can be neglected on the CV surfaces;

4. the axial, u, and radial velocity, v, are continuous through the rotor disc;

5. viscous drag is negligible;

6. w is zero in the upwind flow and outside the wake;

7. the vorticity in the wake is concentrated in line vortices or vortex sheets aligned with the local streamlines in the rotating frame of reference – in other words, the wake vortices rotate rigidly with the blades and vortex lines and streamlines coincide;

8. to derive the local or differential form of Eq. (1), the vorticity piercing the lateral boundaries of the annular CVs intersecting the blade elements must have no effect on the element's thrust.

Assumption no. 7 simplifies the terms involving the trailing vorticity crossing S_D in the impulse derivation. Assumptions no. 3, no. 5, and no. 8 are likewise embedded in the equations derived by LW and are not explicitly required in the analysis to follow. As LW note, a thorough investigation of assumption no. 8 remains an important area of future research, but as yet the assumption remains necessary to recover the Kutta–Joukowsky expression for local thrust that is conventionally employed in blade-element momentum (BEM) analyses. As such, we perpetuate the use of assumption no. 8 for the time being. We also note, emphatically, that none of the eight assumptions places any restrictions on flow expansion. Since the impulse derivation of Eq. (1) is likewise unrestricted, the equation is exact in the presence of flow expansion and for any distribution of w(x).

Although Eq. (1) has been known since Glauert (1935) and appears in modern texts, such as Eq. (4.6) in van Kuik (2018), LW's analysis provides the first proof of its exactness when the trailing vortex sheets have finite thickness. LW's second main finding for circumferentially uniform, expanding flow is

$$\begin{array}{}\text{(2)}& \mathrm{0}=\underset{S_{\mathrm{U}}}{\int }\left(v^{\mathrm{2}}-a^{\mathrm{2}}\right)\mathrm{d}S=\underset{S_{\mathrm{D}}}{\int }\left(v^{\mathrm{2}}-a^{\mathrm{2}}\right)\mathrm{d}S,\end{array}$$

where $a=\mathrm{1}-u$ (when u is normalized by the wind speed, U₀) is the usual axial induction factor. To maintain consistency we use only normalized velocities from here on. The easiest was to do this is to mentally replace the density ρ by $\mathit{\rho }{U}_{\mathrm{0}}^{\mathrm{2}}$.

van Kuik (2020) found Eq. (2) was satisfied by his model of the expanding flow through a wind turbine rotor. When a and v are further assumed to be C₀ continuous on S_U and S_D, Eq. (2) tells us that $\left|a\right|=\left|v\right|$ at some radial location, and LW cite three simulations that show $\left|a\right|\approx \left|v\right|$ near the rotor tip. The vanishing of the first integral on S_U in Eq. (2) is the more general result; the vanishing of the second integral on S_D follows from assumption no. 4 above. Until the end of Sect. 3, we treat the rotor as circumferentially uniform. Since Eq. (1) contains no terms representing pressure redistribution, LW assert that its local version giving the contribution to the thrust at radius x is also exact:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}x}}=\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{w}^{\mathrm{2}}+\mathit{\lambda }wx\right)x\mathrm{d}\mathit{\theta }=\mathrm{2}\mathit{\pi }wx\left({\displaystyle \frac{w}{\mathrm{2}}}+\mathit{\lambda }x\right),\end{array}$$

where θ is the circumferential co-ordinate, defined in Fig. 1. This result is also not new: it is, for example, Eq. (4.24) of van Kuik (2018). It is often referred to as the Kutta–Joukowsky theorem for blade-element thrust because it gives the axial force, $\mathrm{d}T/\mathrm{d}x$, as the product of the circumferential velocity at the rotor, $w/\mathrm{2}+\mathit{\lambda }x$, and the sum of the circulation on all blades, 2πwx. Impulse analysis, however, can also be applied if the CV outlet is moved to the far wake to give $\mathrm{d}T/\mathrm{d}x$ in terms of the w in the far wake. The Glauert (1935) original derivation of Eq. (1) – based on the Bernoulli equation – also suggests the exactness of Eq. (3).

Equation (2) can be derived using standard CV momentum analysis, but the authors are unaware of it appearing in the literature prior to LW. It is a natural outcome of the impulse perspective which we use to investigate the effects of flow expansion on the conventional axial momentum equation. It will be shown that Eq. (2) is closely related to the effects of pressure in the upwind flow on the conventional axial momentum equation and the general relationship between a and the far-wake induction, a_∞. T is derived in Sect. 4 of Sørensen (2016) and Sect. 5.2.4 of van Kuik (2018) using a CV ending in the far wake. We take the different approach of using the CVs shown in Fig. 1 because that choice clarifies the effects of expansion. We also make further use of the impulse form of the T equation. The derivations of the remaining equations in this paper are straightforward and could have been easily done in the past if the impulse perspective had been available.

For context, we now examine the connection between the impulse- and momentum-based approaches to turbine thrust, which requires a relationship between a and w. As explained by, for example, Wood et al. (2021), the Kawada–Hardin (KH) equations for the velocity field of a constant pitch, p, constant radius helical vortex, Kawada (1936) and Hardin (1982), yield

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{p}{x}}={\displaystyle \frac{w/\mathrm{2}}{a}},\end{array}$$

as only half of the near-wake azimuthal velocity is induced by the wake (the other half is due to the blades). Pitch can also be related geometrically to a and λ by treating the wake as a non-expanding rigid helicoidal surface, as done by Okulov and Sørensen (2008). In the limit where λx≫w, we have

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{p}{x}}\approx {\displaystyle \frac{\mathrm{1}-a}{\mathit{\lambda }x}},\end{array}$$

and the preceding two equations can be combined to give $\mathit{\lambda }wx\approx \mathrm{2}a(\mathrm{1}-a)$. The high-λ limit of Eq. (3) thus becomes

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}x}}{|}_{\mathit{\lambda }\to \mathrm{\infty }}\approx \mathrm{2}\mathit{\pi }\mathit{\lambda }w{x}^{\mathrm{2}}\approx \mathrm{4}a(\mathrm{1}-a)\mathit{\pi }x,\end{array}$$

recovering the familiar 4a(1−a) integrand from classical momentum theory. At smaller λ, the relationship between the momentum- and impulse-based thrust expressions has not been fully investigated.

In the next section, we express the contribution of pressure on the expanding upstream streamtube to actuator disc thrust. The section thereafter analyzes the local form of the thrust equation. It contains our second main result about the behaviour of a – that it is negligible at λ=0 and a_∞≈2a is possible at high λ only if Eq. (4) remains valid for expanding flow. In Sect. 4, we apply the Biot–Savart law to an expanding Joukowsky wake, which contains only hub and tip vortices. On the further assumption of constant p, we show, again for the first time, that $a\le a_{\mathrm{\infty }}\le \mathrm{2}a$ depending on the extent of the vortex expansion; the larger the expansion, the closer a approaches a_∞. Not surprisingly, the far-wake radius is reduced as is the power extracted by the turbine. The final two sections contain the general discussion and conclusions, respectively.

2 Actuator disc thrust for expanding flow

Some results of the impulse analysis can be converted easily to conventional equations containing the axial velocity and the pressure on the CV surface even when the flow expands through the rotor. For example, Bernoulli's equation for P_U, the pressure on S_U, is

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\mathrm{2}{P}_{\mathrm{U}}}{\mathit{\rho }}}=\mathrm{1}-v^{\mathrm{2}}-u^{\mathrm{2}}.\end{array}$$

P_U and all pressures considered herein are gauge pressures relative to the free-stream pressure in the wind. Equation (7) allows the removal of v² from Eq. (2) to give

$$\begin{array}{}\text{(8)}& \underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{P_{\mathrm{U}}}{\mathit{\rho }}}x\mathrm{d}x=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}u(\mathrm{1}-u)x\mathrm{d}x,\end{array}$$

which is also the outcome of a conventional momentum balance on CV₂. The momentum balance on CV₁ yields

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{T}{\mathrm{2}\mathit{\pi }\mathit{\rho }}}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}u(\mathrm{1}-u)x\mathrm{d}x-\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{P_{\mathrm{D}}}{\mathit{\rho }}}x\mathrm{d}x,\end{array}$$

where P_D is the pressure on S_D. It is important to note that the effective upper limit on the integrals in Eq. (9) is outside the rotor. Nevertheless,

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{T}{\mathrm{2}\mathit{\pi }\mathit{\rho }}}=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\displaystyle \frac{P_{\mathrm{U}}-P_{\mathrm{D}}}{\mathit{\rho }}}x\mathrm{d}x=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\displaystyle \frac{\mathrm{\Delta }P}{\mathit{\rho }}}x\mathrm{d}x,\end{array}$$

since P_D=P_U for x>1. In other words, there is no pressure jump at z=0 outside the rotor. The thrust equation with integration only over the rotor can be found by rewriting Eq. (9) as

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\displaystyle \frac{T}{\mathrm{2}\mathit{\pi }\mathit{\rho }}}& =\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}a(\mathrm{1}-a)x\mathrm{d}x-\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\displaystyle \frac{P_{\mathrm{D}}}{\mathit{\rho }}}x\mathrm{d}x+\underset{\mathrm{1}}{\overset{\mathrm{\infty }}{\int }}a(\mathrm{1}-a)x\mathrm{d}x\\ & -\underset{\mathrm{1}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{P_{\mathrm{D}}}{\mathit{\rho }}}x\mathrm{d}x.\end{array}\end{array}$$

To remove the last two integrals for x≥1, we use Eq. (7) for P_D=P_U and then Eq. (2) to arrive at

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{T}{\mathit{\pi }\mathit{\rho }}}=\mathrm{2}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}a(\mathrm{1}-a)x\mathrm{d}x-\mathrm{2}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\displaystyle \frac{P_{\mathrm{D}}}{\mathit{\rho }}}x\mathrm{d}x+\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\left(a^{\mathrm{2}}-v^{\mathrm{2}}\right)x\mathrm{d}x.\end{array}$$

The first integral in Eq. (12) contributes half of the conventional thrust. The first and second integrals are components of conventional CV analysis, whereas the third integral is new. It makes Eq. (12) exact for an actuator disc when the flow expands and will be shown below to be generally positive. We now change the CV from that shown in Fig. 1 to the more commonly used one formed by the bounding streamsurface (BS) dividing the flow passing through the rotor from the external flow. BS begins at $z=-\mathrm{\infty }$, where z is the axial co-ordinate with origin at the rotor (Fig. 1). The vertical faces of the new CV are, therefore, subsets of those shown in Fig. 1. A straightforward momentum balance gives

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{T}{\mathit{\pi }\mathit{\rho }}}=\mathrm{2}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}a(\mathrm{1}-a)x\mathrm{d}x-\mathrm{2}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\displaystyle \frac{P_{\mathrm{D}}}{\mathit{\rho }}}x\mathrm{d}x+\mathrm{2}\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}{\left(P{\displaystyle \frac{\mathrm{d}x}{\mathrm{d}z}}\right)}_{\mathrm{BS}}x\mathrm{d}z,\end{array}$$

where the last integrand is evaluated on BS. $\mathrm{d}x/\mathrm{d}z$ gives the local slope of BS, so $P\mathrm{d}x/\mathrm{d}z$ is the axial component of the pressure acting on BS. It follows immediately from Eqs. (12) and (13) that

$$\begin{array}{}\text{(14)}& \underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}{\left(P{\displaystyle \frac{\mathrm{d}x}{\mathrm{d}z}}\right)}_{\mathrm{BS}}x\mathrm{d}z={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\left(a^{\mathrm{2}}-v^{\mathrm{2}}\right)x\mathrm{d}x,\end{array}$$

which gives the first quantification known to the authors of the axial force due to the expanding flow through a wind turbine rotor. It is easy to generalize this equation because there is no thrust extracted in the upwind flow. For any x and z≤0,

$$\begin{array}{}\text{(15)}& \underset{-\mathrm{\infty }}{\overset{z}{\int }}{\left(P{\displaystyle \frac{\mathrm{d}x}{\mathrm{d}z}}\right)}_{S(x,z)}x\mathrm{d}z={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{x}{\int }}\left(a^{\mathrm{2}}-v^{\mathrm{2}}\right)x\mathrm{d}x,\end{array}$$

where S(x,z) is the streamsurface passing through (x,z) so that $\mathrm{BS}=S(\mathrm{1},\mathrm{0})$. The second integral is evaluated at $z\in [-\mathrm{\infty },\mathrm{0}]$.

P_D in Eq. (12) can be evaluated in the standard manner by assuming that the unsteady Bernoulli equation is valid from immediately behind the rotor to the far wake:

$$\begin{array}{}\text{(16)}& \begin{array}{rl}-{\displaystyle \frac{\mathrm{2}{P}_{\mathrm{D}}}{\mathit{\rho }}}& =u^{\mathrm{2}}+v^{\mathrm{2}}+w^{\mathrm{2}}-u_{\mathrm{\infty }}^{\mathrm{2}}-w_{\mathrm{\infty }}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{2}{P}_{\mathrm{\infty }}}{\mathit{\rho }}}\\ & -\mathrm{2}\mathit{\lambda }{x}_{\mathrm{\infty }}{w}_{\mathrm{\infty }}+\mathrm{2}\mathit{\lambda }xw,\end{array}\end{array}$$

where the far-wake terms have the subscript “∞”. The last two terms arise from the unsteady potential terms, evaluated by assuming rigid wake rotation (see appendix B of LW). Conveniently, these terms cancel due to conservation of angular momentum, yielding

$$\begin{array}{}\text{(17)}& -{\displaystyle \frac{\mathrm{2}{P}_{\mathrm{D}}}{\mathit{\rho }}}=u^{\mathrm{2}}+v^{\mathrm{2}}+w^{\mathrm{2}}-u_{\mathrm{\infty }}^{\mathrm{2}}-w_{\mathrm{\infty }}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{2}{P}_{\mathrm{\infty }}}{\mathit{\rho }}}.\end{array}$$

Combining Eqs. (12) and (17) we get

$$\begin{array}{}\text{(18)}& {\displaystyle \frac{T}{\mathit{\pi }\mathit{\rho }}}=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\left(\mathrm{1}-u_{\mathrm{\infty }}^{\mathrm{2}}\right)x\mathrm{d}x-\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\left({\displaystyle \frac{\mathrm{2}{P}_{\mathrm{\infty }}}{\mathit{\rho }}}-w^{\mathrm{2}}+w_{\mathrm{\infty }}^{\mathrm{2}}\right)x\mathrm{d}x,\end{array}$$

where w_∞ and P_∞ are evaluated at x_∞ in the wake, connected to x at the rotor by a mean streamsurface.

In the far wake, the pressure and circumferential velocity are related by

$$\begin{array}{}\text{(19)}& {\displaystyle \frac{\mathrm{d}{P}_{\mathrm{\infty }}/\mathit{\rho }}{\mathrm{d}x}}={\displaystyle \frac{w_{\mathrm{\infty }}^{\mathrm{2}}}{x}}.\end{array}$$

The relationship between the area integrals of P and w can be found using the technique introduced by McCutchen (1985) and rediscovered by Wood (2007): multiply both sides by x² and integrate by parts for the left side. If P_∞x²→0 as x↓0, and is zero at the edge of the far wake, then

$$\begin{array}{}\text{(20)}& \underset{\mathrm{0}}{\overset{R_{\mathrm{\infty }}}{\int }}{\displaystyle \frac{P_{\mathrm{\infty }}}{\mathit{\rho }}}x\mathrm{d}x=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{R_{\mathrm{\infty }}}{\int }}{w}_{\mathrm{\infty }}^{\mathrm{2}}x\mathrm{d}x.\end{array}$$

As pointed out by van Kuik (2018) in conjunction with his Eq. (6.8), any swirl at the edge of the wake makes P_∞(x_∞)≠0. The present analysis can accommodate this behaviour, but for the present we take the simpler path of assuming P_∞(x_∞)=0. The main justification for this assumption is that we expect the magnitude of the swirl to become negligible everywhere at the edge of the wake at high λ. When P_∞(x_∞)=0, Eq. (18) reduces to

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{T}{\mathit{\pi }\mathit{\rho }}}\approx \underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\left(\mathrm{1}-u_{\mathrm{\infty }}^{\mathrm{2}}\right)x\mathrm{d}x+\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{w}^{\mathrm{2}}x\mathrm{d}x.\end{array}$$

Defining the axial induction in the far wake as $a_{\mathrm{\infty }}=\mathrm{1}-u_{\mathrm{\infty }}$, we obtain

$$\begin{array}{}\text{(22)}& {\displaystyle \frac{T}{\mathit{\pi }\mathit{\rho }}}\approx \underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{a}_{\mathrm{\infty }}\left(\mathrm{2}-a_{\mathrm{\infty }}\right)x\mathrm{d}x+\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{w}^{\mathrm{2}}x\mathrm{d}x,\end{array}$$

and the standard thrust equation is recovered if a_∞≈2a and w²≈0, which is typically the case at high λ but may not be generally correct. Note that Eq. (22) is accurate at λ=0, where the first integral is negligible but a_∞≠2a.

To recover the classical thrust equation and to provide a comparison to the analyses of Sørensen (2016) and van Kuik (2018), we now move the downwind face of the CV to the far wake and use Eq. (20). This results in

$$\begin{array}{}\text{(23)}& \begin{array}{rl}{\displaystyle \frac{T}{\mathrm{2}\mathit{\pi }\mathit{\rho }}}& =\underset{\mathrm{0}}{\overset{R_{\mathrm{\infty }}}{\int }}{a}_{\mathrm{\infty }}\left(\mathrm{1}-a_{\mathrm{\infty }}\right)x\mathrm{d}x-\underset{\mathrm{0}}{\overset{R_{\mathrm{\infty }}}{\int }}{\displaystyle \frac{P_{\mathrm{\infty }}}{\mathit{\rho }}}x\mathrm{d}x,\\ & =\underset{\mathrm{0}}{\overset{R_{\mathrm{\infty }}}{\int }}{a}_{\mathrm{\infty }}\left(\mathrm{1}-a_{\mathrm{\infty }}\right)x\mathrm{d}x+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{R_{\mathrm{\infty }}}{\int }}{w}_{\mathrm{\infty }}^{\mathrm{2}}x\mathrm{d}x.\end{array}\end{array}$$

If we ignore the second integral in Eq. (18) and the integrals in Eq. (23) containing P and w, and assume a and a_∞ are constant with x, we again recover the conventional relation a_∞≈2a by invoking conservation of mass.

In considering the local equation for $\mathrm{d}T/\mathrm{d}x$ in the next section, it is useful to have the alternative form of Eq. (23) from the impulse analysis of LW. The direct application of LW's Eq. (22),

$$\begin{array}{}\text{(24)}& {\displaystyle \frac{T}{\mathrm{2}\mathit{\pi }\mathit{\rho }}}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{w}^{\mathrm{2}}+\mathit{\lambda }wx\right)x\mathrm{d}x+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\left(v^{\mathrm{2}}-a^{\mathrm{2}}\right)x\mathrm{d}x,\end{array}$$

together with v=0 everywhere in the far wake and a_∞=0 for r>R_∞ gives

$$\begin{array}{}\text{(25)}& {\displaystyle \frac{T}{\mathrm{2}\mathit{\pi }\mathit{\rho }}}=\underset{\mathrm{0}}{\overset{R_{\mathrm{\infty }}}{\int }}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{w}_{\mathrm{\infty }}^{\mathrm{2}}+\mathit{\lambda }{w}_{\mathrm{\infty }}x\right)x\mathrm{d}x-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{R_{\mathrm{\infty }}}{\int }}{a}_{\mathrm{\infty }}^{\mathrm{2}}x\mathrm{d}x.\end{array}$$

Note that Eq. (24) holds anywhere behind the rotor, i.e. for z>0 with the second integral approaching zero as $z\downarrow +\mathrm{0}$.

The results for the far wake can be used to estimate the conventional thrust when λ=0, and the expansion is negligible by application to S_D. Equation (20) will then be approximately valid for a and P_D replacing a_∞ and P_∞, and comparison with Eq. (1) shows that the momentum flux term 2a(1−a) will be negligible. In other words, the thrust on a stationary disc occurs predominately through the pressure on its back face associated with w.

3 Local thrust in expanding flow

Having considered the thrust for the complete rotor, we now consider the local contribution at radius x. We continue to use a circumferentially uniform disc. It is easy to show that the local form of Eq. (12),

$$\begin{array}{}\text{(26)}& {\displaystyle \frac{\mathrm{1}}{\mathit{\pi }\mathit{\rho }}}{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}x}}=\left[\mathrm{2}a(\mathrm{1}-a)-{\displaystyle \frac{\mathrm{2}{P}_{\mathrm{D}}}{\mathit{\rho }}}+a^{\mathrm{2}}-v^{\mathrm{2}}\right]x,\end{array}$$

is exact for a circumferentially uniform disc in expanding flow. This can be done in at least two ways. First, using Eq. (7) and simple manipulation, the bracketed terms can be written as

$$\begin{array}{}\text{(27)}& \mathrm{1}-u^{\mathrm{2}}-v^{\mathrm{2}}-{\displaystyle \frac{\mathrm{2}{P}_{\mathrm{D}}}{\mathit{\rho }}}={\displaystyle \frac{\mathrm{2}\mathrm{\Delta }P}{\mathit{\rho }}},\end{array}$$

and the pressure difference across the annulus containing the blade elements must give the exact thrust by assumption no. 4. Secondly, starting from Eq. (15) it is easy to prove that a²−v² in Eq. (26) accounts for the difference in pressure acting on the top and bottom of the expanding annular streamtube that intersects the blade elements.

We now consider the consequences of the exact Eq. (26) for the far wake. If the w and P_∞ terms in Eq. (17) are negligible at high λ, the bracketed term in Eq. (26) becomes

$$\begin{array}{}\text{(28)}& \mathrm{2}a(\mathrm{1}-a)-{\displaystyle \frac{\mathrm{2}{P}_{\mathrm{D}}}{\mathit{\rho }}}+a^{\mathrm{2}}-v^{\mathrm{2}}\approx \mathrm{1}-u_{\mathrm{\infty }}^{\mathrm{2}}.\end{array}$$

The exactness of the local form of Eq. (23) is not easy to establish in general because all three velocity components can be important in the wake and the total pressure is not constant. This is the first reason we based our analysis on the CVs shown in Fig. 1 rather than one extending to the far wake. We note, however, that there is no interchange between pressure on BS and axial momentum in the flow outside the far wake where $v^{\mathrm{2}}=a^{\mathrm{2}}=\mathrm{0}$. In other words, the interchange is completed before the far wake is reached. This is the second reason we used the CVs in Fig. 1. The local form of Eq. (23) will have a term corresponding to the bracketed term in Eq. (26) of 2u(1−u_∞). Combining with Eq. (28), we retrieve the standard result that $u=(\mathrm{1}+u_{\mathrm{\infty }})/\mathrm{2}$ or a_∞=2a which can be accurate only at high λ; note that the discussion immediately below Eq. (22) shows the result does not hold at λ=0. Further, from Eq. (7),

$$\begin{array}{}\text{(29)}& {\displaystyle \frac{\mathrm{2}{P}_{\mathrm{U}}}{\mathit{\rho }}}=\mathrm{1}-u^{\mathrm{2}}-v^{\mathrm{2}}=\mathrm{2}a(\mathrm{1}-a)+a^{\mathrm{2}}-v^{\mathrm{2}}.\end{array}$$

If Eq. (6) is valid, then Eq. (26) becomes

$$\begin{array}{}\text{(30)}& -{\displaystyle \frac{\mathrm{2}{P}_{\mathrm{D}}}{\mathit{\rho }}}\approx \mathrm{2}a(\mathrm{1}-a)-a^{\mathrm{2}}+v^{\mathrm{2}}\end{array}$$

for x≤1.

The preceding analysis shows that $P_{\mathrm{U}}\ne -P_{\mathrm{D}}$ in general, in contrast to the familiar results of one-dimensional momentum theory. $P_{\mathrm{U}}=-P_{\mathrm{D}}$ would require a=v, which cannot hold everywhere for several reasons. First, v→0 as x↓0, whereas there is no similar constraint on a. Secondly, we argued above that the flow outside the wake has an axial momentum deficit so a≥0 but not necessarily equal to v for x>1. Equation (2) would then be violated if a=v for x≤1. Thirdly, van Kuik (2018) Sect. 5.4.4 points out that there is no theoretical requirement that $P_{\mathrm{U}}=-P_{\mathrm{D}}$. They are unlikely, however, to differ greatly in general. This suggests v→a as x→1, as argued by LW and shown by the model calculations of van Kuik (2020), who found also that v was significantly larger than a outside the wake until at least x≈1.2. If a>v over most of the rotor, then the positive a²−v² in Eq. (26) corresponds to a positive pressure exerted by the external flow on the wake.

A more definite statement about P_U and P_D can be made for stationary rotors (λ=0) following the last paragraph of the previous section. It is shown there that a is negligible at λ=0 so that P_U≈0 and P_D is associated with w behind the rotor. The inequality reduces as λ increases but is always present because of non-zero a²−v².

We now consider the far wake in more detail to determine the vortex pitch and its relation to a_∞, which are required in the next section. Equation (19) requires $P_{\mathrm{\infty }}/\mathit{\rho }=-w_{\mathrm{\infty }}^{\mathrm{2}}/\mathrm{2}$ when $w_{\mathrm{\infty }}\sim \mathrm{1}/x_{\mathrm{\infty }}$. We assume that at sufficiently high λ, the wake approximates a Joukowsky wake with the hub vortex lying along the axis of rotation and the tip vortices at radius R_∞ in the far wake, with no vorticity in between. The main justification for this assumption comes from Eq. (1). When λ=0, the first term implies that the bound vorticity, Γ, cannot be constant; Wood (2015) showed that Γ∼x². At high λ, however, the first term becomes negligible in comparison to the second for most x. The simplest wake for which the thrust remains bounded on a turbine with N blades occurs when NΓλ∼λwx is constant in x and λ; this is the Joukowsky wake in which assumption no. 7 of Sect. 1 becomes irrelevant to the flow between the tip and hub vortices. Further, the tip vortices now separate the wake and the external flow, which may have very different velocities. The vortex velocity should then be the average of these two, and the vortex lines need not align with the wake streamlines.

Outside the hub vortex core of a Joukowsky wake, $w_{\mathrm{\infty }}\sim \mathrm{1}/x_{\mathrm{\infty }}$ and, as pointed out by Sørensen (2016), the total pressure is constant for all streamsurfaces. In addition, independence will hold in the sense that the integrands in Eqs. (23) and (25) must be equal. Thus,

$$\begin{array}{}\text{(31)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{w}_{\mathrm{\infty }}^{\mathrm{2}}+\mathit{\lambda }{w}_{\mathrm{\infty }}x-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{a}_{\mathrm{\infty }}^{\mathrm{2}}=a_{\mathrm{\infty }}\left(\mathrm{1}-a_{\mathrm{\infty }}\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{w}_{\mathrm{\infty }}^{\mathrm{2}},\end{array}$$

without making any assumption about the relationship between a and a_∞. p_∞, the constant pitch of the constant radius tip vortices, is related to the velocities by the equivalent of Eq. (4): $p_{\mathrm{\infty }}/x_{\mathrm{\infty }}=w_{\mathrm{\infty }}/a_{\mathrm{\infty }}$. Equation (31) can be rewritten as

$$\begin{array}{}\text{(32)}& p_{\mathrm{\infty }}={\displaystyle \frac{\mathrm{1}-a_{\mathrm{\infty }}/\mathrm{2}}{\mathit{\lambda }}}.\end{array}$$

In the next section, λ will be calculated using Eq. (32) for a given p_∞ and the corresponding calculated value of a_∞. Equation (32) is the high−λ equivalent of Eq. (22) of Okulov and Sørensen (2008) for vortex pitch provided the convection velocity of the vortex – w in their notation but w_v here – is equal to $a_{\mathrm{\infty }}/\mathrm{2}$. Table 1 of Wood and Okulov (2017) shows that w_v→a for ideal Betz–Goldstein rotors as λ→∞, and so Eq. (32) is recovered since a_∞→2a in the same limit. Another way to view this result is that the axial velocity in the Joukowsky far wake is constant and equal to 1−a_∞ outside the vortex cores, so the tip vortices must travel downwind at a velocity of $\mathrm{1}-a_{\mathrm{\infty }}/\mathrm{2}$ to be force-free.

The KH equations for a doubly infinite helical vortex of constant radius and pitch lead to

$$\begin{array}{}\text{(33)}& a_{\mathrm{\infty }}=N\mathrm{\Gamma }/\left(\mathrm{2}\mathit{\pi }{p}_{\mathrm{\infty }}\right).\end{array}$$

If we ignore wake expansion, then $a\approx N\mathrm{\Gamma }/\left(\mathrm{4}\mathit{\pi }p\right)$ for the singly infinite near-wake, and if p_∞≈p, then a_∞≈2a. This result suggests the strategy for the next Section, where we analyze the flow associated with expanding tip vortices by assuming they have constant pitch everywhere. This allows a semi-analytic determination of their influence on the flow through the rotor. In other words, we relax one of the limitations of the KH equations, that of constant radius, which must be relaxed, but keep the limitation on p, which, hopefully, leads to results of sufficient generality.

4 The expanding Joukowsky wake with constant pitch

We assume p remains constant and use the results of the previous section and the Biot–Savart law to investigate the flow immediately behind the rotor to determine the thrust and power coefficients. The circumferentially averaged velocities are due entirely to the trailing vorticity: w is due to the hub vortex only, whereas u and v result from the expanding tip vortices only.

4.1 Biot–Savart analysis of expanding tip vortices

Without loss of generality, let the lifting line representing one blade lie instantaneously along the x axis in Fig. 1 and consider the tip vortex beginning at (1, 0, 0). We now determine the velocities induced at a point $(x,\mathit{\theta },\mathrm{0})$ in cylindrical polar co-ordinates or $(x\cos \mathit{\theta },x\sin \mathit{\theta },\mathrm{0})$ in Cartesian co-ordinates for constant p. A point on the vortex is $\left(t\right(\mathit{\beta }),\mathit{\beta },p\mathit{\beta })$ in Cartesian co-ordinates or $\left(t\right(\mathit{\beta })\cos \mathit{\beta },t(\mathit{\beta })\sin \mathit{\beta },p\mathit{\beta })$ in cylindrical polar co-ordinates, where radius t is a monotonically increasing function of the vortex angle β that asymptotes to the far-wake radius. Thus, $\mathrm{1}\le t\le R_{\mathrm{\infty }}$, and from here on, the dependence of t on β will be understood. An increment dl along the vortex is given by

$$\begin{array}{}\text{(34)}& dl=\left(-t\sin \mathit{\beta }+{\displaystyle \frac{\mathrm{d}t}{\mathrm{d}\mathit{\beta }}}\cos \mathit{\beta },t\cos \mathit{\beta }+{\displaystyle \frac{\mathrm{d}t}{\mathrm{d}\mathit{\beta }}}\sin \mathit{\beta },p\right)\mathrm{d}\mathit{\beta },\end{array}$$

and the distance d from the point to the vortex is

$$\begin{array}{}\text{(35)}& d=\left(x\cos \mathit{\theta }-t\cos \mathit{\beta },x\sin \mathit{\theta }-t\sin \mathit{\beta },-p\mathit{\beta }\right),\end{array}$$

so that

$$\begin{array}{}\text{(36)}& d^{\mathrm{2}}=x^{\mathrm{2}}+t^{\mathrm{2}}-\mathrm{2}xt\cos (\mathit{\beta }-\mathit{\theta })+p^{\mathrm{2}}{\mathit{\beta }}^{\mathrm{2}},\end{array}$$

which is an even function of β and θ. A straightforward application of the Biot–Savart law yields the three velocities associated with the trailing tip vortex as

$$\begin{array}{}\text{(37)}& \begin{array}{rl}\left(v\right(x,\mathit{\theta }),& w(x,\mathit{\theta }),a(x,\mathit{\theta }))={\displaystyle \frac{\mathrm{\Gamma }}{\mathrm{4}\mathit{\pi }}}\left(I_{\mathrm{v}},I_{\mathrm{w}},I_{\mathrm{a}}\right)=\\ & {\displaystyle \frac{\mathrm{\Gamma }}{\mathrm{4}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\left(i_{\mathrm{v}}(x,\mathit{\theta }),i_{\mathrm{w}}(x,\mathit{\theta }),i_{\mathrm{a}}(x,\mathit{\theta })\right)}{d^{\mathrm{3}}}}\mathrm{d}\mathit{\beta },\end{array}\end{array}$$

where Γ is the vortex strength,

$$\begin{array}{}\text{(38)}& {\displaystyle }& {\displaystyle }{i}_{\mathrm{v}}(x,\mathit{\theta })=-p\left[t\mathit{\beta }\cos (\mathit{\beta }-\mathit{\theta })+\left(t-\mathit{\beta }{\displaystyle \frac{\mathrm{d}t}{\mathrm{d}\mathit{\beta }}}\right)\sin (\mathit{\beta }-\mathit{\theta })\right],\text{(39)}& {\displaystyle }& {\displaystyle }{i}_{\mathrm{w}}(x,\mathit{\theta })=p\left[x+\left(\mathit{\beta }{\displaystyle \frac{\mathrm{d}t}{\mathrm{d}\mathit{\beta }}}-t\right)\cos (\mathit{\beta }-\mathit{\theta })-\mathit{\beta }t\sin (\mathit{\beta }-\mathit{\theta })\right],\mathrm{and}\text{(40)}& {\displaystyle }& {\displaystyle }{i}_{\mathrm{a}}(x,\mathit{\theta })=t^{\mathrm{2}}-xt\cos (\mathit{\beta }-\mathit{\theta })-x{\displaystyle \frac{\mathrm{d}t}{\mathrm{d}\mathit{\beta }}}\sin (\mathit{\beta }-\mathit{\theta }).\end{array}$$

In forming the circumferential averages by integrating over $\mathrm{0}\le \mathit{\theta }\le \mathrm{2}\mathit{\pi }$, all the sin (β−θ) terms will vanish as they are odd in θ. The linearity of inviscid flow leads to equal contributions to the averaged $(u,w,a)$ from the N identical and equi-spaced trailing vortices.

The simplest calculation of i_a is for x=0, for which the circumferential average $a\left(\mathrm{0}\right)=a(\mathrm{0},\mathit{\theta })$, and

$$\begin{array}{}\text{(41)}& \begin{array}{rl}a\left(\mathrm{0}\right)& ={\displaystyle \frac{N\mathrm{\Gamma }}{\mathrm{4}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{t^{\mathrm{2}}}{{\left(t^{\mathrm{2}}+p^{\mathrm{2}}{\mathit{\beta }}^{\mathrm{2}}\right)}^{\mathrm{3}/\mathrm{2}}}}\mathrm{d}\mathit{\beta },\\ & ={\displaystyle \frac{N\mathrm{\Gamma }}{\mathrm{4}\mathit{\pi }p}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{t^{\mathrm{2}}}{{\left(t^{\mathrm{2}}+z^{\mathrm{2}}\right)}^{\mathrm{3}/\mathrm{2}}}}\mathrm{d}z={\displaystyle \frac{N\mathrm{\Gamma }}{\mathrm{4}\mathit{\pi }p}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{i}_{\mathrm{a}}\left(\mathrm{0}\right)\mathrm{d}z,\\ & ={\displaystyle \frac{N\mathrm{\Gamma }}{\mathrm{4}\mathit{\pi }}}{I}_{\mathrm{a}}\left(\mathrm{0}\right).\end{array}\end{array}$$

I_a is, clearly, dependent only on the geometry of the tip vortices. For an expanding wake with constant p, Eq. (4) will underestimate a as $I_{\mathrm{a}}\left(x\right)\ge I_{\mathrm{a}}\left(\mathrm{0}\right)\ge \mathrm{1}/p$ when t is not constant. If p varied with β, then pβ in Eq. (36) would be replaced by ∫pdβ and the direct relation between $\partial /\partial z$ and $(\mathrm{1}/p)\mathrm{d}/\mathrm{d}\mathit{\beta }$ would be lost. It is likely that an analytic expression for the integrands in Eq. (37) would not be possible.

Performing the θ integration of Eq. (40) using Mathematica gives

$$\begin{array}{}\text{(42)}& {\displaystyle }\begin{array}{rl}{i}_{\mathrm{v}}\left(x\right)& ={\displaystyle \frac{p\mathit{\beta }}{\mathit{\pi }x\sqrt{p^{\mathrm{2}}{\mathit{\beta }}^{\mathrm{2}}+(x+t{)}^{\mathrm{2}}}}}\left[\left(\mathrm{1}+{\displaystyle \frac{m}{\mathrm{2}}}\right)E\left(m_{\mathrm{p}}\right)\right.\\ & \left.-K\left(m_{\mathrm{p}}\right)\right],\mathrm{and}\end{array}\text{(43)}& {\displaystyle }\begin{array}{rl}{i}_{\mathrm{a}}\left(x\right)& ={\displaystyle \frac{-\mathrm{1}}{\mathit{\pi }\sqrt{p^{\mathrm{2}}{\mathit{\beta }}^{\mathrm{2}}+(x+t{)}^{\mathrm{2}}}}}\left[\left(\mathrm{1}+{\displaystyle \frac{m}{\mathrm{2}}}-{\displaystyle \frac{mt}{\mathrm{2}r}}\right)E\left(m_{\mathrm{p}}\right)\right.\\ & \left.-K\left(m_{\mathrm{p}}\right)\right],\end{array}\end{array}$$

where $m=\mathrm{4}xt/(p^{\mathrm{2}}{\mathit{\beta }}^{\mathrm{2}}+(x-t{)}^{\mathrm{2}})$. E(.) and K(.) are the complete elliptic integrals of the second and first kind, respectively, whose argument $m_{\mathrm{p}}=m/(\mathrm{1}+m)$. Thus, v and a can be obtained by integrating Eq. (41) along the trajectory of the tip vortex, t(β) for $\mathrm{0}\le \mathit{\beta }\le \mathrm{\infty }$. This must, in general, be done numerically, but several checks are possible. In describing these, we continue to use the notation $I=\int i\mathrm{d}\mathit{\beta }$ and identify the limits to the integral if they differ from (0,∞).

If t remains constant at 1, say, and the integration is over $-\mathrm{\infty }\le \mathit{\beta }\le \mathrm{\infty }$, that is, for a doubly infinite vortex or vortices of constant radius and pitch, then $I_{\mathrm{v}}(-\mathrm{\infty },\mathrm{\infty })=\mathrm{0}$ for any x, and

$$\begin{array}{}\text{(44)}& \begin{array}{rl}{I}_{\mathrm{a}}(-\mathrm{\infty },\mathrm{\infty })& =\mathrm{2}/p,\mathrm{for}x<\mathrm{1},\\ & =\mathrm{1}/p,\mathrm{for}x=\mathrm{1},\mathrm{and}\\ & =\mathrm{0},\mathrm{otherwise}.\end{array}\end{array}$$

The interior and exterior solutions in Eq. (44) are consequences of the KH equations, derived from the velocity potential. All results in Eq. (44) follow from Eq. (37). Using NIntegrate in Mathematica and MATLAB's integral, these results were reproduced to six significant figures for a similar range of x to that used in the main text and limits of ±1000π on the integration. For a singly infinite helix, the values of I_a for z=0 when β=0 are half those in Eq. (44). These were reproduced numerically to the same accuracy. I_v is not available from the KH equations for this case.

As with any Biot–Savart analysis, the behaviour of Eqs. (42) and (43) as $x\to t\left(\mathrm{0}\right)=\mathrm{1}$ must be considered. As m_p→1, E(m_p)∼1, Eq. (19.6.1) of NIST DLMF (2021), and $K\left(m_{\mathrm{p}}\right)\sim \log (\mathrm{16}/m_{\mathrm{p}}^{\prime })/\mathrm{2}$, where $m_p^{\prime }+m_p=\mathrm{1}$, Eq. (17.3.26) of Abramowitz and Stegun (1964). The leading terms in Eqs. (42) and (43) become

$$\begin{array}{}\text{(45)}& {\displaystyle }\begin{array}{rl}{i}_{\mathrm{v}}\left(x\right)\sim & p\mathit{\beta }\left(p^{\mathrm{2}}{\mathit{\beta }}^{\mathrm{2}}+x^{\mathrm{2}}+t^{\mathrm{2}}\right)/\sqrt{p^{\mathrm{2}}{\mathit{\beta }}^{\mathrm{2}}+(x+t{)}^{\mathrm{2}}}/x\\ & \left(p^{\mathrm{2}}{\mathit{\beta }}^{\mathrm{2}}+(x-t{)}^{\mathrm{2}}\right),\mathrm{and}\end{array}\text{(46)}& {\displaystyle }\begin{array}{rl}{i}_{\mathrm{a}}\left(x\right)\sim & -\left(p^{\mathrm{2}}{\mathit{\beta }}^{\mathrm{2}}+x^{\mathrm{2}}-t^{\mathrm{2}}\right)/\sqrt{p^{\mathrm{2}}{\mathit{\beta }}^{\mathrm{2}}+(x+t{)}^{\mathrm{2}}}/\\ & \left(p^{\mathrm{2}}{\mathit{\beta }}^{\mathrm{2}}+(x-t{)}^{\mathrm{2}}\right),\end{array}\end{array}$$

showing that a logarithmic singularity occurs in i_a despite it being the integrand for the circumferentially averaged axial velocity. This is a stronger singularity than that in Chattot (2020) perturbation analysis of the flow near the edge of the rotor, which assumes a vortex cylinder wake. There is no logarithmic singularity in i_v, but the slope $\mathrm{d}{i}_{\mathrm{v}}/\mathrm{d}\mathit{\beta }$ increases without bound as β→0. These behaviours could be mitigated by using the well-known cut-off modification to the limits of the Biot–Savart integrals as was done for helical vortices by Ricca (1994); see also Sect. 11.2 of Saffman (1992). There is, however, a simpler, heuristic alternative. The upper limit on a(x) as x→1 is taken to be a_∞. A partial justification for this tactic comes from the wind tunnel measurements of a model wind turbine by Krogstad and Adaramola (2012). Their Fig. 9c shows that at λ=9.51, a≈0 at small x but rises to the extraordinary value of around 0.8 at x=1. Thus, $I_{\mathrm{a}}\le \mathrm{2}/p$ was enforced in the calculations described in the main text. Whenever this was done, I_v was assumed equal to the maximum value below the limit on I_a.

The numerical evaluation of I_a and I_v can be improved by considering the asymptotic behaviour of i_a and i_v for large β, which corresponds to small m and m_p. The leading terms are simple functions of β, allowing the infinite integrals to be approximated. For I_v, we have

$$\begin{array}{}\text{(47)}& I_{\mathrm{v}}\left(x\right)\approx I_{\mathrm{v}}(x,\widehat{\mathit{\beta }})+R_{\mathrm{v}}\left(\widehat{\mathit{\beta }}\right),\end{array}$$

where the first term was obtained numerically over $\mathit{\beta }=[\mathrm{0},\widehat{\mathit{\beta }}]$ and the remainder, $R_{\mathrm{v}}\left(\widehat{\mathit{\beta }}\right)$, is an approximation to the integral over $\mathit{\beta }=[\widehat{\mathit{\beta }},\mathrm{\infty }]$. $R_{\mathrm{v}}\left(\widehat{\mathit{\beta }}\right)$ is

$$\begin{array}{}\text{(48)}& R_{\mathrm{v}}\left(\widehat{\mathit{\beta }}\right)=x{R}_{\mathrm{\infty }}^{\mathrm{2}}/\left(\mathrm{2}{\widehat{\mathit{\beta }}}^{\mathrm{3}}{p}^{\mathrm{4}}\right).\end{array}$$

The remainder for I_a is independent of x:

$$\begin{array}{}\text{(49)}& R_{\mathrm{a}}\left(\widehat{\mathit{\beta }}\right)=R_{\mathrm{\infty }}^{\mathrm{2}}/\left(\mathrm{2}{\widehat{\mathit{\beta }}}^{\mathrm{2}}{p}^{\mathrm{3}}\right).\end{array}$$

This result also follows from Eq. (41) when z≫R_∞.

It was found that $\widehat{\mathit{\beta }}=\mathrm{200}\mathit{\pi }$ was sufficient to ensure six-figure accuracy of the integrals over the range of x considered below. I_v converged faster than I_a, reaching 99 % of the final value by β=2π for any x.

The Biot–Savart integrands in Eq. (43) are plotted in Figs. 2 and 3 for x close to the blade tip, in terms of axial distance z=pβ, where β is the vortex angle starting from zero at the rotor. The figures also show the small-β asymptotes in Eqs. (45) and (46) and the large-β remainders defined in Eqs. (48) and (49). If the tip vortex radius t remains at 1, Eq. (44) gives $I_{\mathrm{a}}=\mathrm{1}/p$ for any x, and the conventional momentum equation (Eq. 6) remains valid. We assume that $\mathrm{1}/p$ is the minimum value of I_a, and, as explained in the Appendix, we impose a≤a_∞ so that $\mathrm{1}/p\le I_{\mathrm{a}}\le \mathrm{2}/p$. For maximum power, the familiar derivation of the Betz–Joukowsky limit suggests $R_{\mathrm{\infty }}^{\mathrm{2}}\approx \mathrm{2}$, so we investigate R_∞ around that value. Note, however, the use of Eq. (6) to derive this limit means that it is applicable only to a wake that expands either very slowly, as explained above, or very rapidly to $t=\sqrt{\mathrm{2}}$, as $I_{\mathrm{a}}=\mathrm{1}/p$ for any constant t. We will show that generic wind turbine wakes at high λ expand at a rate that is intermediate between these extremes, which causes Eq. (6) to be inaccurate. There is no direct maximization of power output in the following analysis. Instead, the wake model is constrained as we now describe.

Figure 2 Integrand, i_v, for p=0.1, $R_{\mathrm{\infty }}^{\mathrm{2}}=\mathrm{1.597}$. ◦, x=0.9; □, $x=\mathrm{0.99},$; and ⋄, x=0.999 from Eq. (42). i_v increases with x. × is the integrand in Eq. (45). For clarity, only every second data point is plotted. The solid line shows the remainders from Eq. (48). The differences with varying x are within the thickness of the line.

Figure 3 Integrand, i_a, for p=0.1, $R_{\mathrm{\infty }}^{\mathrm{2}}=\mathrm{1.597}$. ◦, x=0.9; □, $x=\mathrm{0.99},$; and ⋄, x=0.999 from Eq. (43). i_a increases with x. × is the integrand in Eq. (46). For clarity, only every second data point is plotted. The solid line shows the remainders from Eq. (49).

Solving Eq. (37) for I_a and I_v requires p and the tip vortex trajectory. We used the very simple form:

$$\begin{array}{}\text{(50)}& t=R_{\mathrm{\infty }}-\left(R_{\mathrm{\infty }}-\mathrm{1}\right)\exp (-k\mathit{\beta }),\end{array}$$

which satisfies three necessary conditions: t=1 when β=0, t→R_∞ for large β, and t and its derivatives are continuous. The fourth condition is that k must satisfy the reduced version of Eq. (2):

$$\begin{array}{}\text{(51)}& \underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\left(I_{\mathrm{v}}^{\mathrm{2}}-I_{\mathrm{a}}^{\mathrm{2}}\right)x\mathrm{d}x=\mathrm{0}.\end{array}$$

This integral will be called the “expansion integral”. It uniquely fixes k for any choice of $R_{\mathrm{\infty }}^{\mathrm{2}}$ and p. I_a and I_v were obtained using the MATLAB function integral over $\mathit{\beta }=[\mathrm{0},\mathrm{200}\mathit{\pi }]$ to an absolute tolerance of 10⁻⁶. The remainders, Eqs. (48) and (49), were then added. The expansion integral and the mass flux integral described below were found by trapezoidal integration using the points shown in Fig. 4. The expansion integral is large for small k as v is (not obviously) maximized when there is very little vortex expansion near the rotor.

Figure 4 Axial induction, a, and radial velocity, v, for the conditions in Table 1. p=0.05: a, ⋄; v, ×. The corresponding results for p=0.10 are the dashed lines. Note that −a is shown to make it distinct from v. The solid lines are −a and v for no expansion, i.e. k=0, R_∞=2. The x axis is logarithmic.

The mass flux through the rotor, using Eq. (33) to remove NΓ, determines a_∞:

$$\begin{array}{}\text{(52)}& \mathrm{1}-a_{\mathrm{\infty }}p\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{I}_{\mathrm{a}}x\mathrm{d}x=\left(\mathrm{1}-a_{\mathrm{\infty }}\right)R_{\mathrm{\infty }}^{\mathrm{2}}.\end{array}$$

Equation (32) then yields λ. A number of possible methods were considered for solving the integral in Eq. (52). i_a(x,θ) can be written as

$$\begin{array}{}\text{(53)}& i_{\mathrm{a}}(x,\mathit{\theta })={\displaystyle \frac{\mathrm{d}}{\mathrm{d}x}}\left({\displaystyle \frac{x}{\mathrm{d}}}\right)-{\displaystyle \frac{p^{\mathrm{2}}{\mathit{\beta }}^{\mathrm{2}}}{d^{\mathrm{3}}}},\end{array}$$

which allows an analytic integration of i_a(x,θ)x in x. The resulting expression is complicated and probably requires numerical integration in θ and β to obtain the mass flux. Further, the integrand is singular at a point that varies with θ and β. The simpler alternative of numerical integration of I_ax was used.

To find the unique $R_{\mathrm{\infty }}^{\mathrm{2}}$, we impose the further condition that k must match the slope of the vortex surface at the rotor. Then k in Eq. (50) equals k_*, given by

$$\begin{array}{}\text{(54)}& {\displaystyle \frac{\mathrm{d}t}{\mathrm{d}z}}(\mathit{\beta }=\mathrm{0})={\displaystyle \frac{v(x=\mathrm{1})}{\mathrm{1}-a(x=\mathrm{1})}}=\left(R_{\mathrm{\infty }}-\mathrm{1}\right){\displaystyle \frac{k_{*}}{p}}.\end{array}$$

4.2 Results

The results in Table 1 were obtained using the MATLAB routine patternsearch to minimize the single objective function that combined the magnitude of the expansion integral and $\mathrm{|}k-k_{*}\mathrm{|}$. This, surprisingly, occurred at a constant value of $k_{*}/p$, implying that the vortex expansion to the far-wake radius happens over a fixed distance and the surface containing the vortices has the same shape, independent of p or λ.

Table 1 Results for the expanding Joukowsky wake with constant pitch.

C_T was calculated from Eq. (1) with w² ignored because λ is large:

$$\begin{array}{}\text{(55)}& C_{\mathrm{T}}\approx N\mathrm{\Gamma }\mathit{\lambda }/\mathit{\pi }\approx \mathrm{2}{a}_{\mathrm{\infty }}p\mathit{\lambda }\approx \mathrm{2}{a}_{\mathrm{\infty }}\left(\mathrm{1}-a_{\mathrm{\infty }}/\mathrm{2}\right),\end{array}$$

using Eqs. (32) and (33). We note that Eq. (1) makes the high-λ blade-element thrust constant across the rotor, whereas the familiar form involving the axial velocity equation in Eq. (6) requires a significant variation near the tip. From conservation of angular momentum, and finding the power as the product of torque and angular velocity,

$$\begin{array}{}\text{(56)}& C_{\mathrm{P}}\approx C_{\mathrm{T}}\left(\mathrm{1}-a_{\mathrm{\infty }}\right)R_{\mathrm{\infty }}^{\mathrm{2}},\end{array}$$

so the power extraction also decreases significantly near the tip. Equation (56) and the third component of Eq. (55) also hold for the conventional analysis that leads to the Betz–Joukowsky limit.

From Table 1, the biggest change from the familiar Betz–Joukowsky wake is the 20 % reduction in $R_{\mathrm{\infty }}^{\mathrm{2}}$, which occurs because $a>a_{\mathrm{\infty }}/\mathrm{2}$ for much of the rotor (Fig. 4). In other words, more of the expansion occurs upwind of the rotor. In contrast, the maximum C_P is reduced by only 6 % to 0.557 and the mass flux is increased by 2 %. The second to last column in Table 1 gives $C_{\mathrm{T}}^{\prime }$, the conventional thrust coefficient evaluated from Eq. (6):

$$\begin{array}{}\text{(57)}& C_{\mathrm{T}}^{\prime }=\mathrm{8}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}a(\mathrm{1}-a)x\mathrm{d}x.\end{array}$$

The positive values of $\mathrm{\Delta }{C}_{\mathrm{T}}=C_{\mathrm{T}}^{\prime }-C_{\mathrm{T}}$ in Table 1 indicate the conventional equation overestimates the thrust by around 5 %. From Eqs. (29) and (30), we get the integrated asymmetry of the disc pressure distribution at high λ as

$$\begin{array}{}\text{(58)}& \underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\left({\displaystyle \frac{P_{\mathrm{U}}}{\mathit{\rho }}}+{\displaystyle \frac{P_{\mathrm{D}}}{\mathit{\rho }}}\right)x\mathrm{d}x\approx \underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\left(a^{\mathrm{2}}-v^{\mathrm{2}}\right)x\mathrm{d}x.\end{array}$$

Using the data in Table 1, the integral of δ_P is 0.087, so the magnitude of P_D is generally significantly less than that of P_U. It was shown in the previous section that the pressure integrals are equal in magnitude in the minimally expanding wake when λ=0, but the analysis in this section shows divergence in the expanding Joukowsky wake at high λ.

The integrands i_v and i_a Figs. 2 and 3 are large in the vicinity of the rotor. Their size implies that the simple assumed shape of the tip vortex trajectory, Eq. (50), is reasonable and that adding a term or terms, say, to control the approach to the far wake would not change the analysis significantly. The low- and high-β asymptotic approximations to i_v and i_a are accurate at sufficiently large x for an appropriate range of β but are not sufficiently general to yield analytic approximations to the integrals. The final figure, Fig. 5, plots the ratio of a and a_∞. The large differences from a_∞=2a occur near the edge of the rotor.

Figure 5 Ratio of axial induction at the rotor, a, to value in the far wake, a_∞. p=0.05; ⋄. p=0.10; dashed line. Note that the x axis is logarithmic.

Figure 4 shows that a and v at the rotor for the cases in Table 1 are independent of p. a(0)=0.296, which is less than the Betz–Joukowsky value of $\mathrm{1}/\mathrm{3}$, and v(0)=0 as it must. The results for v and $a=\mathrm{1}/\mathrm{3}$ in an unexpanding wake ($k_{*}=\mathrm{0}$) are plotted as solid lines. Note that a=0 for x>1, but v is very large near the edge of the rotor, and it is clear that the expansion integral cannot be satisfied. The limit a≤a_∞ was applied near the blade tip for the expanding wake, where v has increased to be nearly equal but smaller than a. Outside the wake, v>a and takes until x=3 to fall to 0.03. Similarly shaped distributions of a and v for a Joukowsky wake are shown in Fig. 5 of van Kuik (2020), who also found that Eq. (2) was satisfied in his low-λ simulations.

In the final figure, Fig. 5, the ratio $a/a_{\mathrm{\infty }}$ documents the departure from the conventional relationship a_∞=2a near the rotor edge. Note that a_∞ is constant for a Joukowsky wake, so the ratio is non-zero for x>1.

5 Discussion

The pressure in the expanding flow ahead of a wind turbine contributes to the axial force on the rotor and a momentum deficit in the flow outside the rotor. Researchers have been aware of these two effects for many years, but the present analysis provides the first quantitative determination of them in Eqs. (14) and (15) derived by approximating the rotor as an actuator disc. The effects on the disc in integral and incremental form depend on a²−v², where v is the normalized radial velocity and a is the usual axial induction factor. Further, v²−a² can be used to quantify the external flow disturbed by the wind turbine and so may be useful to the study of multiple rotors in close proximity, as analyzed by, for example, Branlard and Meyer Forsting (2020). One way to do this is by defining I_E as

$$\begin{array}{}\text{(59)}& I_{\mathrm{E}}=\underset{x_{\mathrm{BS}}}{\overset{\mathrm{\infty }}{\int }}\left(v^{\mathrm{2}}-a^{\mathrm{2}}\right)x\mathrm{d}x=-\underset{\mathrm{0}}{\overset{x_{\mathrm{BS}}}{\int }}\left(v^{\mathrm{2}}-a^{\mathrm{2}}\right)x\mathrm{d}x,\end{array}$$

where x_BS is the radius of BS at any z≤0. The last integral is a consequence of Eq. (2) being valid for S_U lying anywhere in the upwind flow. I_E must be zero in the undisturbed upwind flow. It then increases to its maximum value at the rotor according to the present analysis. I_E then decreases in the wake to be zero in the far wake. In other words, the perturbation to the external flow is complete by the time the far wake is reached. Note that the second equality in Eq. (59) does not hold in the wake.

The impulse analysis of Limacher and Wood (2020) (LW) showed that the Kutta–Joukowsky (KJ) equations for rotor thrust, Eq. (1), and for the blade-element contributions to the thrust, Eq. (3), are exact in the presence of wake expansion, where “exact” means using no more assumptions or approximations than the eight listed in the Introduction. The KJ equations, containing only the circumferential velocity and tip speed ratio, are not equivalent to the conventional equation involving only the axial velocity, when the flow expands. This is the outcome of the analysis in Sect. 4, where an expanding Joukowsky wake comprising tip and hub vortices of constant pitch was analyzed. The conventional thrust equation is altered by around 5 %–10 %, depending on the trajectory of the tip vortices because the geometrical relation in Eq. (4) is modified by the expansion.

The first three sections of the paper used only the standard form of control volume (CV) analysis for axial momentum to determine the thrust of the rotor and the incremental thrust of the blade elements comprising the rotor. To clarify the effects of expansion, most analysis in this paper used CVs with downwind faces in the immediate vicinity of the rotor, as opposed to their common placement in the far wake. The rotor and the flow are assumed to be circumferentially uniform. We argued in the Introduction that the impulse analysis provides a simple and novel perspective on the role of the pressure. The thrust equations derived in Sect. 3 for the rotor, and in Sect. 4 for the local flow at any radius, contain the pressure acting on the downwind face of the actuator disc, which must be removed to make the equations suitable for actual blade analysis. Removal can be done accurately only for very low tip speed ratios where the expansion and its effects are small.

To the rotor thrust, the pressure along the bounding streamsurface adds a term containing the integral of a²−v² over the rotor. This integral is equal and opposite the integral for the flow outside the wake, so there is no net contribution to the thrust determined using the CVs shown in Fig. 1. Unsurprisingly, the corresponding term in the local thrust equation at any x also contains a²−v². It follows that the conventional local thrust equation implies a≈v, but a² is generally larger than v² over the rotor, but more precise estimates of v do not appear to be possible. a²>v² implies that the pressure adds to the rotor thrust and is associated with a momentum deficit in the external flow. The common derivation of the axial momentum equation, which leads to the Betz–Joukowsky limit, ignores the interaction of pressure and external momentum and then ignores the radial velocity in relating the pressure at the rear of the disc to the far wake. These errors cancel, so the main failing of the conventional equation is the breakdown of the relation a_∞=2a when expansion is significant. The previous section suggests this breakdown is due to the expanding tip vortices at high λ in the Joukowsky wake. At the rotor, the slope of the streamsurface containing the tip vortices is 53^∘ for maximum power extraction, Table 1, so their trajectory is intermediate between very slow and very rapid expansion, either of which would require a_∞=2a. This analysis used Eq. (32) for the pitch of the tip vortex, found by moving the CV outlet to the far wake and using LW's impulse equation for thrust. We note that van Kuik (2020) estimated the streamsurface angle at the rotor edge to be 46^∘, which is close to the present value. The effect of the expansion on a was constrained so that a≤a_∞ as an alternative to using a cut-off in the Biot–Savart integral. Figure 4 shows that a increased with radius to reach a_∞ in the streamtube bounded by the tip vortex, suggesting a very substantial effect of expansion. Qualitatively, this large value of a is in agreement with the wind tunnel measurements of Krogstad and Adaramola (2012), who found that a increased across the rotor to reach 0.8 at the tip at high λ.

Including v in the axial momentum equation effectively adds an extra unknown to the conservation equations that may render them useless unless another equation for u, v, or w could be derived. Further, high v may cause significant alterations to the lift and drag of the blade elements near the tip. To our knowledge, radial velocity effects on airfoil lift and drag have not been studied in the context of blade-element theory.

The role of the radial velocity and flow expansion is probably more complicated in rotors with a limited number of blades than the actuator discs considered here. Eriksen and Krogstad (2017) measured u, v, and w immediately behind the rotor of a model three-bladed turbine out to a radius 20 % larger than the blade tip radius. They used phase-locked averaging to obtain the flow field as seen by an observer rotating with the blades. Significant phase variations occurred in a and v, showing that the averages a² and v² over a blade revolution could be large even if the mean values of a and v are small. Nevertheless, the magnitude of both a and v was largest near the angular location of the blades, suggesting that the issues with radial deflection will occur in real turbines. We hope that these comments, and the present analysis, will inspire further measurements to be made far enough outside the wake to help clarify the role of flow expansion and the disturbances to the external flow.

6 Conclusion

This analysis started from the impulse-derived Kutta–Joukowsky equation for wind turbine thrust, which does not involve the axial velocity. The equation is valid for any amount of expansion in the upwind flow and the wake and any distribution of bound circulation on the rotor. We were able to

• demonstrate the conventional thrust equation containing the axial velocity can be correct only when the tip speed ratio is large;

• derive an exact expression for the effects of flow expansion on the conventional momentum equation – this involves the axial induction factor and the radial velocity;

• apply the conventional and impulse thrust equations in the far wake to give the pitch of the tip vortices in the Joukowsky wake in terms of the tip speed ratio and the far-wake induction;

• find a semi-analytic solution of the Biot–Savart law for the induced velocities at the rotor by assuming the tip vortex had constant pitch – the axial velocity near the rotor tip approached the far-wake value, but was prevented from exceeding it as an alternative to using the familiar cut-off in the Biot–Savart integrals, and the increase in the rotor value contradicts the familiar relation that the axial induction factor everywhere at the rotor is half that of the far wake; and

• derive in Sect. 5 the following results from the model of constant pitch, expanding tip vortices:

  1. the angle of the tip vortex surface to the wind direction was 53^∘ for maximum power production, independently of the tip speed ratio and vortex pitch;

  2. because it is neither very small nor very large, this expansion leads to an error of around 6 % in the conventional thrust equation, which would be accurate for both extreme expansions;

  3. the resulting wake expands less than the familiar Betz–Joukowsky wake – for two pitch values corresponding to tip speed ratios of 7 and 14, the far-wake area was 1.59 times the rotor area;

  4. we find the reduction in the rotor power and thrust due to expansion – the maximum power coefficient and corresponding thrust coefficient were 6 % less than the values given by the Betz–Joukowsky limit; and

  5. we quantify the influence of the expansion on the flow outside the rotor – for example, the radial velocity at three rotor radii is still 3 % of the wind speed when the rotor is producing maximum power, and the axial induction factor decays to zero more rapidly than the radial velocity as radius.