1 Introduction↩︎

Actuator disk model is one of the most fundamental models in fluid dynamics for analysing turbine rotors and propellers [1]. First proposed by Rankine and later refined into its classical form by Froude in 1889 [2], it is now commonly referred to as Froude’s actuator disk theory. The model represents the complex loading of a rotor as a simple pressure jump across an infinitely thin, permeable disk of the same diameter, known as an actuator disk. Despite its simplicity, the theory remains a cornerstone of rotor fluid dynamics because it provides clear and intuitive insights into rotor performance. Froude’s actuator disk theory is typically one of the very first conceptual models introduced to students and researchers interested in the fluid mechanics of turbomachinery. It appears in almost every major wind energy textbook [3]–[8], and classical aerodynamics and hydrodynamic textbooks on propellers and rotors [9]–[12].

Despite its historical importance and value as a simple conceptual model, Froude’s theory has important limitations. First, it does not describe how flow variables evolve with streamwise distance and only relates conditions far upstream and far downstream to those on the actuator disk. Secondly, the model does not describe the turbulent wake behaviour since the predicted downstream velocity asymptotes to \(U_0 (1 - 2a)\), where \(U_0\) is the incoming velocity and \(a\) is the induction factor, and the wake width becomes constant after an initial expansion. In reality, following the initial velocity reduction due to pressure recovery, the wake gradually recovers through turbulent mixing, and wake expansion continues downstream. For these reasons, Froude’s theory is generally regarded as conceptually relevant only in the region immediately downstream of the disk, before turbulence dominates.

The limitations become especially severe for highly-loaded actuator disks (with induction factors \(a > 0.5\)), where the model even predicts negative wake velocities. Likewise, the well-known relation for the thrust coefficient \(C_T = 4a(1-a)\) breaks down in this regime and fails to provide meaningful results [5]. This breakdown is believed to result from strong shear that drives the wake into a fully turbulent state and promotes strong interactions with the ambient flow, effects absent from Froude’s formulation. To mitigate these limitations, empirical corrections are often introduced to adjust the \(C_T\)–\(a\) relationship [3], [13]. For a highly-loaded actuator disk (i.e., a disk with low porosity), where a significant portion of the fluid bypasses the disk rather than passing through it, the actuator disk behaves more like a solid plate. In this case, intense wake turbulence and flow separation generate a strong negative pressure behind the disk, contributing to pressure drag [14]. [15] accounted for this negative wake pressure by formulating governing equations for a control volume (CV) surrounding the disk that has non-zero pressure at the outlet, unlike Froude’s original formulation. The flow inside the CV is assumed to be inviscid and irrotational, and the outlet is placed immediately before the region where turbulent mixing becomes significant. However, because introducing a non-zero outlet pressure adds an extra unknown, additional information is needed to close the system of equations. To address this, they modelled the actuator disk using potential flow theory. As representing the disk as a distribution of potential sources results in an unrealistic velocity discontinuity at the disk plane, they rescaled the wake velocities to enforce mass conservation and validated their predictions against water-channel experiments of porous flat plates. More recently, [16] adopted a similar CV analysis but determined the negative pressure at the CV outlet by solving the two-dimensional pressure Poisson equation proposed by [17]. They assumed that the CV outlet coincides with the end of the near wake, whose length was estimated using the method of [18]. They also extended this framework to yawed rotors to predict a \(C_T\)–\(a\) relationship that shows good agreement with large-eddy simulation (LES) data for various cases. While these approaches provide improved predictions for highly-loaded disks, they still do not overcome the important limitation of Froude’s actuator disk theory; namely, its inability to incorporate effects of turbulent mixing induced by strong velocity shear that may occur in the immediate downstream vicinity of a highly-loaded disk.

Turbulent wake development is more realistically represented in far-wake analytical models, which are widely used in wind energy applications to characterise wake interactions in wind farms [19], [20]. These models are generally based on momentum theory, where pressure effects are neglected under the far-wake assumption, so that the thrust force balances the streamwise flux of momentum deficit [21]. Although recent wake models aim to use more realistic velocity profiles in the near wake [22]–[25], their neglect of pressure effects prevents them from accurately capturing the flow immediately behind the disk, where such effects remain significant. Consequently, analytical wake models typically require assumptions about near-wake conditions as inputs. Examples include specifying the onset of the far-wake region, the initial wake width, or the downstream location where the velocity deficit reaches its maximum. Moreover, commonly used far-wake models typically fail to provide reliable predictions for induction factors \(a > 0.25\) in the far-wake region [26]. For instance, the widely used top-hat model of [27] and the Gaussian model of [21] may predict a decrease in the maximum wake velocity deficit as the induction factor \(a\) increases beyond 0.25 (i.e., \(C_T > 0.75\)), which is not physically expected. Although more recent studies [26] have incorporated pressure effects into the far-wake evolution of highly-loaded disks, they remain limited to the far-wake region, where velocity monotonically increases with streamwise distance from the disk. This behaviour is not valid in the region immediately behind the disk.

This paper aims to bridge the gap between two traditionally separate areas of research, actuator disk theory and turbulent wakes. We show that it is possible to overcome the limitations outlined above by developing a new actuator disk theory that includes wake recovery due to turbulent entrainment. While still idealised, the proposed approach provides a more accurate and physically meaningful representation of actuator-disk flows.

2 Generalised actuator disk theory including wake development↩︎

We consider the “hybrid” CV around a disk as shown in Figure 1. Upwind of the disk, as in Froude’s theory, we assume there is no mass exchange between the CV and the surroundings, so the CV is the streamtube encompassing the rotor disk. This assumption implies that turbulent mixing does not play a significant role in shaping the flow distribution upstream of the disk. This notion is supported by previous studies showing that inviscid solutions often provide satisfactory predictions in this region [28]–[30]. Downstream of the rotor disk, however, the CV is not a streamtube but is assumed to coincide with the outer boundary of the wake, where turbulence drives entrainment of the ambient flow through the lateral surface of the CV, until the wake is fully recovered. For this CV, at each streamwise position \(x\), the velocity \(U(x)\) and pressure \(P(x)\) are assumed to be uniform, with the diameter of the CV cross-section denoted by \(\sigma(x)\). The disk is located at \(x=0\), and the diameter of the CV at \(x=0\) is equal to the disk diameter \(D\). Subscripts \(0\) and \(\infty\) denote the asymptotic far-upstream (\(x\!\to\!-\infty\)) and far-downstream (\(x\!\to\!\infty\)) values, respectively.

In Froude’s original formulation, the CV extends from far upstream to far downstream. Here, to capture the streamwise variations of the flow quantities, we retain the inlet at far upstream but place the CV outlet at an arbitrary streamwise location \(x\), where \(x\) can take on positive or negative values. The goal of our analysis is to determine \(P(x)\), \(U(x)\), \(\sigma(x)\), and thrust coefficient \(C_T\) for a given induction factor \(a\) (or for a known \(C_T\), to determine \(a\)).

2.1 Mass conservation↩︎

The conservation of mass for this CV is written as \[\label{eq:continuity95new} \dot{m}(x)=\dot{m}_0+\dot{m}_e(x),\tag{1}\] where \(\dot{m}_0=\frac{\pi}{4} \rho \sigma^2_0 U_0\) and \(\dot{m}(x)=\frac{\pi}{4} \rho \sigma^2(x) U(x)\) are the mass flow rates at the inlet and outlet, respectively, and \(\dot{m}_e(x)\) represents the mass flow rate of the entrained flow across the lateral surface. Also, \(\rho\) is the air density. If we define the entrainment velocity as \(U_e(x)\), the axial rate of change of \(\dot{m}_e(x)\) can be written as the product of perimeter times entrainment velocity and is given by \[\label{eq:entrainment95mass} \frac{\mathrm{d}\dot{m}_e}{\mathrm{d}x}=\rho\pi \,\sigma(x) \, U_e(x) \, \mathrm{H}(x),\tag{2}\] where \(\mathrm{H}(x)\) is the Heaviside function, which ensures that the entrainment driven by turbulence only takes place downwind of the actuator disk. Taking the derivative of equation 1 with respect to \(x\) and using equation 2 , we obtain \[\label{eq:sigma2u95new} \frac{\mathrm{d}\left(\sigma^2 U\right)}{\mathrm{d}x}=4\sigma U_e\mathrm{H}(x),\tag{3}\] which can be expanded and rearranged as \[\label{eq:dsigmadx} \frac{\mathrm{d}\sigma}{\mathrm{d}x} = \frac{1}{2U} \left( - \sigma\frac{\mathrm{d}U}{\mathrm{d}x} + 4U_e\mathrm{H}(x)\right).\tag{4}\]

2.2 Choice of entrainment velocity \(U_e\)↩︎

2.2.1 Wake-shear driven entrainment↩︎

In the absence of ambient turbulence, turbulent entrainment is primarily driven by wake shear. In this case, the entrainment velocity is modelled as \[\label{eq:entrainment95wake95shear} \frac{U_e^w(x)}{U_0} = E_1 \frac{\Delta U(x)}{U_0},\tag{5}\] where \(\Delta U(x) = U_0 - U(x)\), the superscript \(w\) stands for wake-shear driven entrainment, and \(E_1\) is the entrainment coefficient with a typical value of 0.1 [31]–[33]. Entrainment is expected to be primarily wake driven in the near wake where the wake shear is strong. However, further downstream it is increasingly influenced by background turbulence as the wake recovers in the far wake. This mechanism is described in the next section.

2.2.2 Background-turbulence driven entrainment↩︎

Regarding the influence of background turbulence on wake development, [34] showed that the incoming turbulence intensity \(I\) becomes the dominant factor governing entrainment in the far wake, while being only weakly dependent on the inflow integral length scale. Therefore, we assume that the entrainment velocity driven by background turbulence is proportional to the incoming turbulence intensity. We note that the turbulence intensity can be defined based on either the total turbulent kinetic energy or its streamwise component. While the former may be more physically grounded, the latter is easier to measure and is typically used in most current wind energy research [20]. We therefore write

\[\frac{U_e^{b}}{U_0} = E_2 I, \label{eq:entrainment-atmospheric}\tag{6}\] where \(E_2\) is the background-driven entrainment coefficient, and superscript \(b\) stands for background turbulence. From the analysis in section 3.2, we find that \(E_2\) should be in the range of \(0.25\)–\(0.4\), with \(E_2 = 0.3\) used in this work as a representative value.

2.2.3 Total entrainment driven by wake shear and background turbulence↩︎

We propose that the total entrainment velocity is proportional to the generalised mean of the normalized velocity deficit and the incoming turbulence intensity: \[\frac{U_e(x)}{U_0} = \left[ \left( \frac{U_e^w}{U_0}\right)^n + \left( \frac{U_e^{b}}{U_0} \right)^n \right]^{1/n} = \left[ \left( E_1 \frac{\Delta U(x)}{U_0} \right)^n + \left( E_2 \, I \right)^n \right]^{1/n}, \label{eq:entrainment95velocity}\tag{7}\] where \(n\) is a positive integer. Here, we assume \(n=4\), noting that larger values do not lead to a considerable difference. The advantage of using equation 7 is that the entrainment velocity closely approaches the larger of the two contributing factors. Consequently, under typical conditions, the near-wake entrainment velocity is dominated by the wake shear, which is generally larger in that region, while in the far wake the entrainment velocity is dominated by background turbulence.

2.3 Momentum balance↩︎

The momentum equation for the CV shown in figure 1 is given by \[\label{eq:momentum} \rho\frac{\pi}{4}\sigma^{2}(x)\, U(x) \, [U_{0}-U(x)] =\rho\frac{\pi}{8} C_{T}U_{0}^{2} D^{2} \mathrm{H}(x)+\frac{\pi}{4}\sigma^{2}(x) P(x) -F_{P_s}(x),\tag{8}\] where the momentum deficit flux term on the left-hand-side is written using the mass conservation (equation 1 ), i.e. \(\dot{m}_e=\frac{\pi}{4} \rho (\sigma^2 U- \sigma^2_0 U_0)\) and its associated momentum flux \(\dot{m}_e U_0\). In equation, 8 the term \(F_{P_s}(x)\) is the axial component of the force exerted by the pressure on the lateral surface of the streamtube, given by \[\label{eq:F95p95s} F_{P_s}(x)= \frac{\pi}{4}\int_{-\infty}^{x}P(x') \, \frac{\,\mathrm{d}\sigma^2(x')}{\,\mathrm{d}x'} \,\mathrm{d}x',\tag{9}\] where \(x'\) is a dummy variable. Equation 9 is obtained by integrating the pressure force exerted on an infinitesimal lateral area, shown in figure 2 (a), from far upstream to \(x\). In Froude’s actuator theory, the lateral pressure term is neglected because the CV’s outlet is assumed to be sufficiently far downstream, such that the lateral pressure force contribution downstream cancels that upstream, resulting in a net effect of zero. This zero net contribution in the far wake can be understood by considering a spherical control volume of very large radius, \(R_s\!\to\! \infty\), surrounding the disk, as shown in figure 2 (b). For this CV, the lateral pressure forces are all internal [1], and the pressure at the control surface is atmospheric everywhere. Consequently, equation 8 reduces to a simple balance between the thrust force and the net momentum flux across the control surface. However, it is important to note that for a finite \(x\) (either positive or negative) considered here, the net contribution of the lateral pressure term is non-zero and must be retained in the momentum equation.

Taking the derivative of equation 8 with respect to \(x\) and rearranging terms leads to \[\label{Eq:mom95new} U\frac{\mathrm{d}U}{\mathrm{d}x}+ \frac{1}{\rho}\frac{\mathrm{d}P}{\mathrm{d}x}=\left(\frac{2U}{\sigma}\frac{\mathrm{d}\sigma}{\mathrm{d}x}+\frac{\mathrm{d}U}{\mathrm{d}x}\right)(U_0-U)-\frac{1}{2} C_{T}U_{0}^{2} \mathrm{\delta}(x),\tag{10}\] where \(\delta(x)\) is the Dirac delta function. Using equation 4 to express \(\,\mathrm{d}\sigma/\,\mathrm{d}x\) allows equation 10 to be written as \[\label{eq:mom95new2} U(x)\frac{\mathrm{d}U}{\mathrm{d}x}+ \frac{1}{\rho}\frac{\mathrm{d}P}{\mathrm{d}x}=-\frac{1}{2} C_{T}U_{0}^{2} \mathrm{\delta}(x) + \frac{4U_e(x)[U_0-U(x)]\, \mathrm{H}(x)}{\sigma(x)}.\tag{11}\] Integrating from \(-\infty\) to \(x\) results in the following Bernoulli-type equation \[\label{eq:Bernoulli95modified2} \frac{1}{\rho}\,P(x)+\frac{1}{2}U^2(x)=\frac{1}{2}U_0^2\left(1- C_T \mathrm{H}(x) \right)+ 4 \int_0^{\mathrm{max}(x,0)} \frac{U_e(x')\,[U_0-U(x')]}{\sigma(x')} \, \mathrm{d}x',\tag{12}\] which represents the variation of the mechanical energy with \(x\). It can also be regarded as a generalised form of the Bernoulli equation, with two additional effects incorporated: the extraction of energy by the disk (i.e., \(-C_T\mathrm{H}(x)\) term on the right-hand side) and the entrainment of energy by turbulence (the last term on the right-hand side). When \(C_T=U_e = 0\), this equation reduces to a Bernoulli equation. It is also important to note that the derivation of this generalised form of Bernoulli equation from mass and momentum equations is possible when the momentum equation is written in its complete form including the lateral pressure term \(F_{P_s}\).

Finally, equation 11 can be rearranged as \[\label{eq:dudx} \frac{\mathrm{d}U}{\mathrm{d}x} = \frac{1}{U} \left( -\frac{1}{\rho} \frac{\mathrm{d}P}{\mathrm{d}x} - \frac{1}{2} C_{T} U_0^2\, \delta(x) + \frac{4U_e(U_0 - U)\, \mathrm{H}(x)}{\sigma} \right).\tag{13}\] We then have a system of two ordinary differential equations comprising equation 13 for \(\mathrm{d}U/\mathrm{d}x\) and equation 4 for \(\mathrm{d}\sigma/\mathrm{d}x\). However, there are three dependent variables, \(U(x)\), \(\sigma(x)\), and \(P(x)\), with \(x\) as the independent variable. The required third equation (for \(\mathrm{d}P/\mathrm{d}x\)) is obtained in the next section.

2.4 Pressure variations for lightly-loaded actuator disks↩︎

For a lightly-loaded actuator disk (i.e. small values of the induction factor, \(a\)), the axisymmetric pressure Poisson equation can be simplified to the Laplace equation \(\nabla^2P=0\) as a leading-order approximation (see Appendix 5 for details). In the half-space upstream of the rotor disk the Dirichlet boundary condition on the disk \(P(x=0^-)=P^+_D\) is used, while on the downstream half-space one uses \(P(x=0^+)=P^-_D\), where \(P_D^{+}\) and \(P_D^{-}\) are the values of pressure immediately upstream and downstream of the actuator disk, respectively. Along the centreline, the solution to the Laplace equation is given by \[\label{eq:PPE95solution95semiinf95f61g610} P(x)=P_D^{+}\left(1+\frac{x}{\sqrt{x^2+R^2}}\right)\mathrm{H}(-x)+P_D^{-}\left(1-\frac{x}{\sqrt{x^2+R^2}}\right)\mathrm{H}(x),\tag{14}\] where \(R=D/2\) is the rotor radius. We note that our approach here differs from the one in [17] and [16] who used a Cartesian 2D formulation, as we use the axisymmetric formulation of the pressure Poisson equation that is appropriate in three dimensions. The interested reader is referred to Appendix 5 for additional information about the derivation, underlying assumptions and simplifications leading to equation 14 . From equation 12 with \(U_D=U_0(1-a)\), \(P_D^{+}\) and \(P_D^{-}\) are respectively given by \[\label{eq:P95D} P_D^+ = \frac{1}{2}\rho U_0^2a(2-a), \quad P_D^- = \frac{1}{2}\rho U_0^2\left(a(2-a) - C_T\right).\tag{15}\] To provide support for the proposed \(P(x)\) expression in equation 14 , we compare the predicted results with those from numerical simulations. The interested reader is referred to [35] and [36] for more information on the numerical setup, which consists of a LES with laminar uniform inflow towards an actuator disk, and negligible effects of turbulence in the near-rotor region. Results are shown in figure 3 for three different values of \(a\). The model (dashed lines) requires specification of \(C_T\) for the downstream portion of the flow. As will be shown later in §2.6, for the case \(U_e=0\), the classical Froude relation \(C_T=4a(1-a)\) holds. Since for \(x/D<3\) the effects of turbulence were negligible in the simulations, we use \(C_T=4a(1-a)\) to specify \(C_T\) in this comparison with idealised simulation results. Results shown in figure 3 confirm that equation 14 yields satisfactory predictions of the pressure distribution, with accuracy improving as \(a\) decreases. It is worthwhile mentioning that, according to equation 15 , the pressure drop across the actuator disk is not necessarily symmetric. In fact, for typical values of \(C_T\), the magnitude of the immediate upstream pressure increase, \(|P_D^+|\), is greater than the immediate downstream pressure drop, \(|P_D^-|\). This subtle fact is also important for determining the pressure variations (see Appendix 5 for details).

Taking the derivative of equation 14 and inserting values of \(P_D^{+}\) and \(P_D^{-}\) from equation 15 gives \[\label{eq:dp47dx} \frac{\,\mathrm{d}P}{\,\mathrm{d}x} = \frac{1}{2} \rho U_0^2 \left( \frac{R^2}{(x^2 + R^2)^{3/2}} \left[ a(a - 2) \operatorname{sgn}(x) + C_T\,\mathrm{H}(x) \right] - C_T \, \delta(x) \right),\tag{16}\] where \(\operatorname{sgn}(x)\) is the sign function, defined as \(\operatorname{sgn}(x) = -1\) if \(x<0\), and \(\operatorname{sgn}(x) = 1\) if \(x>0\). We now have the three relations needed to obtain the streamwise evolution of velocity \(U(x)\), pressure \(P(x)\), and CV diameter \(\sigma(x)\) for given values of the induction factor \(a\) and thrust coefficient \(C_T\).

2.5 Solving for \(U(x)\), \(\sigma(x)\), and \(P(x)\)↩︎

Here, we determine \(\sigma(x)\), \(U(x)\), and \(P(x)\) for given \(C_T=C_T(a)\). Inserting \(\,\mathrm{d}P/\,\mathrm{d}x\) from 16 in equation 13 and using 4 , we obtain a system of equations that can be solved to find \(U(x)\) and \(\sigma(x)\) \[\tag{17} \begin{align} \frac{\mathrm{d}\sigma}{\mathrm{d}x} &= \frac{1}{2U} \left( - \sigma\frac{\mathrm{d}U}{\mathrm{d}x} + 4U_e\mathrm{H}(x) \right), \tag{18} \\[6pt] \frac{\mathrm{d}U}{\mathrm{d}x} &= \frac{U_0^2}{2U} \left( \frac{R^2}{(x^2 + R^2)^{3/2}} \left( a(2-a)\, \operatorname{sgn}(x) - C_T \, \mathrm{H}(x) \right) + \frac{8}{\sigma}\frac{U_e}{U_0}\frac{U_0 - U}{U_0}\, \mathrm{H}(x) \right). \tag{19} \end{align}\] For the upwind region (\(x<0\)), equation 17 is readily solved, yielding \[\tag{20} \begin{align} \frac{U(x<0)}{U_0} &= \sqrt{a(a-2)\left(\frac{x}{\sqrt{x^2+R^2}}+1\right)+1}, \tag{21} \\[6pt] \frac{\sigma(x<0)}{D} &= \sqrt{\frac{U_0(1-a)}{U(x)}} = \sqrt{1-a} \,\left[ a(a-2)\left(\frac{x}{\sqrt{x^2+R^2}}+1\right)+1 \right]^{-1/4}. \tag{22} \end{align}\]

For the downwind region (\(x>0\)), no analytical solution was found due to the presence of \(U_e(x)\). Instead, a basic forward marching numerical scheme (or more efficient Runge–Kutta methods) can be used, starting at \(x=0\). The initial values \(U(0)=U_0(1-a)\) and \(\sigma(0)=D\) can be used to solve equation 17 numerically and compute \(U(x)\) and \(\sigma(x)\). Pressure \(P(x)\) is directly obtained from equation 14 .

2.6 Determining the thrust coefficient \(C_T\)↩︎

So far, we have not determined a general relationship between \(a\) and \(C_T\). Next, we apply the proposed generalised actuator disk model with turbulent entrainment to develop a relationship for \(C_T\) and establish a framework in which all three flow quantities (\(U,\sigma\), and \(P\)) and \(C_T\) are determined for a given value of \(a\). The obtained relations can also be inverted to obtain \(a\) for any prescribed \(C_T\), since it is the latter quantity that is generally known a-priori for a given rotor. Section 2.3 discussed that the net contribution of the lateral pressure term in the momentum equation must asymptote to zero as \(x\!\to\!\infty\). By enforcing this condition, a new relation for \(C_T\) can be obtained. The total force from the side pressure contribution is given by \(F_{P_s}(x\to \infty)\) from equation 9 and implies: \[\label{eq:F95p95s2} \frac{\pi}{4}\int_{-\infty}^{0} P(x) \, \frac{\,\mathrm{d}\sigma^2(x)}{\,\mathrm{d}x} \,\mathrm{d}x + \frac{\pi}{4}\int_{0}^{\infty}P(x) \, \frac{\,\mathrm{d}\sigma^2(x)}{\,\mathrm{d}x} \,\mathrm{d}x = 0.\tag{23}\] Denoting the first and second integrals as \(F_{P_s}^+\) and \(F_{P_s}^-\), respectively, we can evaluate the upwind lateral pressure force as \[\label{eq:f95ps43} F_{P_s}^+=\frac{\pi}{4}\int_{-\infty}^{0} P(x) \, \frac{\,\mathrm{d}\sigma^2(x)}{\,\mathrm{d}x} \,\mathrm{d}x = \frac{\pi}{8}\rho D^2U_0^2a^2,\tag{24}\] where we have used equations 14 and 15 for \(P(x)\) and 22 for \(\sigma(x)\). From equations 23 (\(F_{P_s}^- = - F_{P_s}^+\)) and 24 , we conclude that \[\label{eq:f95ps4361f95ps-} \int_0^{\infty}P(x) \frac{\,\mathrm{d}\sigma^2(x)}{\,\mathrm{d}x}\,\mathrm{d}x= -\frac{1}{2}\rho D^2U_0^2a^2,\tag{25}\] which can be expanded using equations 14 and 15 to \[\int_0^{\infty}\left(a(2-a) - C_T\right)\left(1-\frac{x}{\sqrt{x^2+R^2}}\right)\frac{\,\mathrm{d}\sigma^2(x)}{\,\mathrm{d}x}\,\mathrm{d}x=-D^2a^2.\] Solving for \(C_T\) gives \[\label{eq:CT95new} C_T=2a+\left(\frac{1}{Y}-1\right)a^2,\tag{26}\] where \[\label{eq:Y} Y=\frac{1}{D^2}\int_0^{\infty}\left(1-\frac{x}{\sqrt{x^2+R^2}}\right)\frac{\,\mathrm{d}\sigma^2(x)}{\,\mathrm{d}x}\,\mathrm{d}x.\tag{27}\] To solve equation 26 with \(Y\) given by equation 27 , \(\sigma(x)\) is required. However, from equation 17 , determining \(\sigma(x)\) itself requires knowledge of \(C_T\). Consequently, equations 17 and 26 must be solved iteratively. We begin by guessing a value for \(C_T\) (e.g., using \(C_T = 4a(1-a)\) or other models such as [13], [15], [16]) and use it in equation 17 to compute \(\sigma(x)\). The resulting \(\sigma(x)\) is then substituted into equations 27 and 26 to obtain an updated \(C_T\). This updated value is fed back into equation 17 , and the process is repeated until convergence is achieved.

If instead of \(a\), we are provided with \(C_T\), we begin by guessing a value of \(a\) and use it in equation 17 to compute \(\sigma(x)\). The resulting \(\sigma(x)\) is then substituted into equation 27 to find the value of \(Y\), which is inserted in equation 26 to find a new value for \(a\), and the process is repeated until convergence is achieved. Implementations of these iterative procedures can be found in the computational notebooks associated with figures 6 or 8 (see the figure captions for details).

2.7 Reduction to Classical Froude Theory for \(U_e=0\)↩︎

Solving equation 17 for \(U_e=0\) at \(x>0\), we obtain \[\label{eq:downind95solution95ue610} \frac{U(x>0)}{U_0}= \sqrt{(1-a)^2+\big[(2a-a^2)-C_T\big]\frac{x}{\sqrt{x^2+R^2}}},\tag{28}\] and \(\sigma(x)=D\sqrt{U_0(1-a)/U(x)}\). For the upwind region (\(x<0\)), the solution is the same as equation 20 . From equation 28 , the asymptotic velocity at \(x\to\infty\) is \(U_{\infty}=U_0\sqrt{1-C_T}\), or \(U_{\infty}=U_0(1-2a)\) if \(C_T=4a(1-a)\). The well-known relation of \(C_T=4a(1-a)\) can be also derived from our new generalised theory for \(U_e=0\). Inserting equation 28 into equation 27 yields \(Y = -1 + \tfrac{2(1-a)}{\,a(2-a)-C_T\,}\left(\sqrt{1-C_T} - (1-a)\right)\), which can then be substituted into equation 26 , whose only non-trivial solution becomes \(C_T = 4a(1-a)\).

3 Results and discussions↩︎

3.1 Model predictions and trends↩︎

The influence of entrainment velocity on model predictions is first illustrated in figure 4, which presents results for various values of the entrainment coefficients \(E_1\), with \(E_2 = 3E_1\) and \(I = 5\%=0.05\). Results are presented for two different values of induction factor \(a\) of \(0.25\) and \(0.45\). When \(E_1 =E_2= 0\), there is no turbulent entrainment, and the far-wake velocity approaches \(U_0(1 - 2a)\), as shown in figure 4. As expected, with increasing the level of turbulent entrainment, the wake velocity recovers more rapidly for both cases shown. However, the effect of entrainment on \(\sigma\) in the near wake differs between the two induction factors. At \(x < 2.5D\), turbulent entrainment slightly increases \(\sigma\) for the lower induction factor (\(a=0.25\)), whereas for the higher induction factor (\(a = 0.45\)) turbulent entrainment significantly reduces the near-wake width compared to the case without turbulence. This behaviour can be interpreted using equation 18 . According to this relation, \(\sigma\) in the near-wake region may increase due to two mechanisms: (i) turbulent entrainment (second term on the right-hand side, \(4U_e\mathrm{H}(x)\)), and (ii) flow deceleration occurring behind the disk (first term on the right-hand side, \(-\sigma \,\,\mathrm{d}U/\,\mathrm{d}x\)). However, as shown in figure 4, flow deceleration behind the disk is less in the presence of turbulent entrainment. Consequently, it has opposing effects on \(\sigma\), and the net outcome, whether \(\sigma\) increases or decreases, depends on which one dominates.

Figure 5 shows how the entrainment velocity varies with the streamwise distance downwind of the actuator disk for a typical case with \(a = 0.3\) and \(I = 8\%\), using values of \(E_1 = 0.1\) and \(E_2 = 0.3\). According to equation 7 , the total entrainment velocity \(U_e\) consists of two components: (i) the wake-shear driven entrainment velocity \(U_e^w\) (equation 5 ), and (ii) the background-turbulence-driven entrainment velocity \(U_e^b\) (equation 6 ). The figure shows that, as expected, the turbulent entrainment and wake recovery are mainly driven by wake shear in the near wake, while background turbulence becomes the dominant mechanism of wake recovery in the far-wake region. The transition between these two regimes occurs at approximately \(5D\) for this case, which lies within the typically reported range of transition from the near wake to the far-wake region [18].

Model predictions for five different values of the induction factor \(a\) are shown in figure 6. For all cases, the entrainment coefficients are set to \(E_1 = 0.1\) and \(E_2 = 0.3\), and the incoming turbulence intensity is \(I = 5\%\). Figure 6a shows the streamwise variation of velocity. In all cases, as expected, the velocity upstream decreases from \(U_0\) to \(U_0(1-a)\) at the disk. Immediately downstream, it drops further, mainly due to pressure recovery (\(\,\mathrm{d}p/\,\mathrm{d}x > 0\) resulting in \(\,\mathrm{d}U/\,\mathrm{d}x<0\) according to equation 13 ), until reaching a minimum value (i.e., maximum deficit) approximately one rotor diameter downstream. Our model predicts the maximum deficit to be smaller than \(2aU_0\) predicted by Froude’s theory, due to the inclusion of turbulent entrainment already beginning at \(x=0\). Further downstream, when pressure is mostly recovered, flow entrainment becomes the dominant mechanism, resulting in \(\,\mathrm{d}U/\,\mathrm{d}x > 0\) according to equation 13 . Moreover, as expected, higher values of \(a\) produce greater velocity deficits.

As shown in figure 6b, the cross-sectional width \(\sigma\) of the CV generally increases with \(x\) in both the upwind and downwind regions. According to equation 4 , velocity reduction (\(\,\mathrm{d}U/\,\mathrm{d}x < 0\)) contributes to cross-sectional expansion both upstream and downstream, with downstream growth further enhanced by flow entrainment. Figure  6b also shows that \(\sigma\) increases with \(a\) downstream, indicating a wider wake for actuator disks with higher loading. In far-wake engineering models (for instance, [21]), this impact of disk loading on the wake width is often incorporated through an semi-empirical initial wake width that depends on the thrust coefficient \(C_T\).

An interesting feature is observed for \(a = 0.5\) and \(a = 0.6\), where the wake undergoes a rapid initial expansion followed by a contraction and then gradually re-expands due to turbulent entrainment. Such non-monotonic behaviour has been extensively reported in the literature ([6], [37], [38], among others) and is often attributed to the formation of complex turbulent wake structures, including vortex rings, characteristic of the turbulent state with \(a>0.4\). It is interesting that our simpler model reproduces such features. Nonetheless, we should bear in mind that the developed pressure relation, equation 14 , is an approximate solution of the pressure Poisson equation for lightly-loaded disks. Therefore, model predictions for large values of \(a\) should be interpreted with caution. Finally, figure 6c shows the pressure distribution, where the asymmetric pressure drop on the disk plane, discussed earlier in section 2.4, is evident, particularly for large values of \(a\). The figure suggests that pressure is recovered downstream within the first two rotor diameters. This behaviour may however differ from some findings in the literature on highly-loaded disks [26], [38], where pressure persists longer. As mentioned before, it is likely because the simplified pressure relation used here is a first-order approximation for lightly-loaded disks.

3.2 Far-wake solution↩︎

Here, we examine model predictions for the asymptotic far-wake case, where \(x \to \infty\). In this case, both the normal and lateral pressure force terms in the momentum equation 8 vanish, and \(U\to U_0\). From the momentum equation 8 for large values of \(x\), we obtain \[\mathrm{d} \left(\sigma^{2} U (U_{0}-U) \right)\approx \mathrm{d} \left(\sigma^{2} (U_{0}-U)\right)\approx 0.\] This means that \[\label{eq:sigma95asymp95intermediate} U_0\mathrm{d}\sigma^{2}\approx \textrm{d}\left(\sigma^2 U\right).\tag{29}\] Inserting equation 29 in the conservation of mass expression in 3 leads to \[\label{eq:sigma95far95wake} \frac{\,\mathrm{d}\sigma}{\,\mathrm{d}x}= 2\frac{U_e}{U_0}\quad \mathrm{if}\; x\to\infty.\tag{30}\]

3.2.1 Entrainment driven by wake shear↩︎

First, we assume the ambient turbulence is negligible, and the main driving force for turbulent entrainment is the wake shear. In this case, \(E_2=0\), so \(U_e=E_1(U_0-U)\). From the simplified form of the momentum equation 8 for \(x\to \infty\) (i.e., no pressure effects and \(U\to U_0\)), we substitute \(U_0-U\) to have \[\label{eq:U95e:case1} U_e=E_1(U_0-U)=\left(\frac{1}{2}E_1C_TU_0D^2\right)\sigma^{-2}.\tag{31}\] Inserting equation 31 in equation 30 and integrating leads to \[\label{eq:sigma95inf95case1} \frac{\sigma}{D}=\left(3E_1C_T\right)^{1/3}\left(\frac{x-x_0}{D}\right)^{1/3},\tag{32}\] where \(x_0<0\) is the virtual origin defined as the location where \(\sigma(x_0)=0\). There are two interesting points about equation 32 . First, it shows that the far wake grows with downstream distance as \(x^{1/3}\) in the absence of ambient turbulence. This is consistent with the extensive literature of three-dimensional wakes subject to laminar free-stream flows (see [39] and references therein). Moreover, this shows that the spreading rate of the wake depends on \(C_T\). This supports George’s theory [40], [41] that although canonical free shear flows may reach universal asymptotic states in terms of self-similarity, their spreading rate is not universal, and it depends on upstream (i.e., initial) conditions. This is clear in figure 7 (a), which shows variations of \(\sigma^3\) with respect to \(x\) for different values of \(a\), where \(E_1=0.1\) and \(E_2=0\). To verify equation 32 , we rearrange it assuming that \(x_0\) is chosen such that \(\sigma(x=0)=D\). Under this assumption, equation 32 simplifies to \((\sigma/D)^3 - 1 \propto x\), which further reduces to \((\sigma/D)^3 \propto x\) when \(\sigma \gg D\), as shown in figure 7 (a).

3.2.2 Entrainment driven by both wake shear and background turbulence↩︎

For the case where ambient turbulence is not negligible, equation 7 simplifies to \(U_e/U_0 = E_2 I\) as \(x \to \infty\). Substituting this relation into equation 30 and integrating yields \[\label{eq:sigma95inf95case2} \frac{\sigma}{D} = 2k_\infty \left(\frac{x - x_0}{D}\right),\tag{33}\] where \(k_{\infty} = \frac{1}{2} \,\mathrm{d}\sigma_{\infty} / \,\mathrm{d}x = E_2 I\) is the expansion rate of the far-wake radius. The parameter \(k_{\infty}\) is an empirical constant widely used in far-wake engineering models of wind turbines [42]. [43] suggested \(k_{\infty} = 0.4 I\), which corresponds to \(E_2 = 0.4\). [44] reported that \(k_{\infty} = 0.02\) provides results in closer agreement with operational data from offshore wind farms. For offshore turbulence intensity at turbine hub height in the range of \(5\%\)–\(7.5\%\), this value of \(k_{\infty}\) translates to \(E_2 = 0.25\)–\(0.4\). Here, we use \(E_2=0.3\). Variations of \(\sigma/D-1\) with \(x/D\) are shown in figure 7 (b) for different values of incoming turbulence intensity \(I\). In this figure, \(E_1 = 0.1\), \(E_2 = 0.3\), and \(a = 0.3\) are used. Similar to the approach used in figure 7 (a) and to verify equation 33 , here we plot \(\sigma/D - 1\), which is basically the normalised increase of the wake width with respect to its initial value. The results illustrate that linear wake growth begins earlier for higher values of incoming turbulence intensity.

For the same parameter set, the corresponding values of \(k(x) = \frac{1}{2} \,\mathrm{d}\sigma / \,\mathrm{d}x\) are shown in figure 7 (c) as solid lines, with the asymptotic value \(k_{\infty}\) included as a dashed line for reference. This figure confirms that for small values of incoming turbulence intensity \(I\), the wake expansion rate \(k\) asymptotes to \(k_{\infty}\) only at large downwind distances. The wake initially expands at a higher rate, and for example for \(I = 4\%\) the wake expansion rate has not yet reduced to its asymptotic linear-growth value even after \(20D\). This variable wake recovery rate is consistent with observations in the literature [45], [46].

3.3 Relation between thrust coefficient and induction factor↩︎

Next, we discuss the predictions of thrust coefficient \(C_T\) and power coefficient \(C_P = C_T (1-a)\) based on the developed actuator disk model and the iterative process explained in Section 2.6. Figure 8 shows variations of \(C_T\) and \(C_P\) with the induction factor \(a\) for different values of \(E_1 = (0.02, 0.05, 0.1)\) and \(E_2 = (0, 0.3)\). For comparison, the figure also shows Froude’s relation (\(C_T=4a(1-a)\)), LES data of [38], and the models of [15] and [16]. The figure shows that in the windmill state (\(a < 0.4\)), model predictions are not very sensitive to the value of entrainment coefficients and provide results similar to Froude’s theory, albeit with slightly smaller values. On the other hand, in the turbulent state (\(a > 0.4\)), results are highly sensitive to the value of \(E_1\), and the inclusion of entrainment considerably improves predictions of actuator disk models in comparison to the unphysical predictions of Froude’s theory in this region, where higher values of \(E_1\) lead to higher values of \(C_T\). A value of \(E_1 = 0.05\) seems to provide results in satisfactory agreement with the model of [16], validated against LES data, while a value of \(E_1 = 0.1\) provides results in very good agreement with the model of [15], validated against laboratory experiments. It is worth noting that the iterative method described in Section 2.6 does not converge for cases with high values of \(a\) (e.g., \(a > 0.6\)) and small values of \(E_1\) (e.g., \(E_1 = 0.02\)). This may be due to the fact that, in highly loaded cases, turbulent mixing and entrainment play such an important role in shaping the wake structure downstream of the disk that minimising their effect does not lead to meaningful predictions. Furthermore, as mentioned earlier, the simplified pressure relation used in the model may not provide accurate estimates for highly-loaded disks. Nevertheless, figure 8 clearly illustrates how including the effect of turbulent entrainment improves predictions of \(C_T\) and \(C_P\) significantly beyond the classical Froude theory’s limiting value of \(a \sim 0.4\).

In reality, the entrainment coefficient \(E_1\) may not be strictly constant and may have some dependencies on the actuator disk characteristics, inflow properties and streamwise location (i.e., \(E_1 = E_1(x/D,a,I)\)). However, the key physical insight provided by figure 8 is that a universal relation between \(a\) and \(C_T\) may not exist; rather, their exact relationship depends on the interaction of the actuator disk flow with the ambient flow and the turbulent entrainment, particularly immediately downstream of the disk. This is especially important for highly-loaded disks, where the interaction between the wake and the surrounding flow is strong and cannot be neglected. This dependence may also explain why observations of \(C_T = C_T(a)\) performed in different studies typically show good agreement in the windmill state, while the data become widely scattered in the turbulent state, as shown in figure 8 and reported in previous works [13].

Finally, the figure indicates that the value of \(E_2\) has minimal impact on the variation of \(C_T\) and \(C_P\). This is expected, as \(E_2\) mainly governs entrainment in the far-wake region where background turbulence effects dominate, whereas the near-wake region, where wake-shear driven entrainment mainly occurs, has a more direct influence on the thrust force and power output. Background-turbulence driven entrainment is expected to have some possible impact on \(C_T\) only in the asymptotic case of \(a \to 0\), since in this case the entrainment induced by ambient turbulence can dominate even in the region immediately behind the turbine.

3.4 Optimal power and revised Betz limit↩︎

Here, we examine predictions of the new theory for the maximum power coefficient \(C_{P,\max}\) and its corresponding optimal induction factor \(a_{\text{opt}}\). Figures 9 (a) and (b) show zoomed-in portions of the \(C_T\)–\(a\) and \(C_P\)–\(a\) curves near \(a_{\text{opt}}\), while figures 9 (c) and (d) illustrate the variations of \(C_{P,\max}\) and \(a_{\text{opt}}\) as functions of \(E_1\). The results indicate that as the turbulence parameter \(E_1\) increases, \(a_{\text{opt}}\) exceeds its classical value of \(1/3\), and \(C_{P,\max}\) slightly surpasses the Betz limit of \(16/27 \approx 0.593\).

Since \(C_P = C_T(1-a)\), it is useful to analyse how \(C_T\) varies with \(E_1\) for \(a\) close to \(a_{\text{opt}}\). By definition, \(\tfrac{1}{2}\rho U_0^2 C_T = P_D^+ - P_D^-\), where the upstream pressure \(P_D^+\) is unaffected by downstream turbulent mixing (see equation 15 ). Thus, an increase in \(|P_D^-|\) leads to a larger \(C_T\). From equations 14 and 25 , we conclude that \(|P_D^-|\), and therefore \(C_T\), increases when turbulent mixing reduces the initial near-wake expansion, which occurs for \(a > 1/3\). An example of this reduction by turbulent mixing for a high value of \(a\) is shown in figure 4 (b). A larger \(C_T\) for \(a > 1/3\) in turn yields a larger \(a_{\text{opt}}\) and a slightly higher \(C_{P,\max}\), which was also reported in [16]. Exceeding the Betz limit has been reported in other contexts, such as flow confinement [47], whereas here the mechanism arises from turbulent mixing that lessens near-wake expansion and thereby increases the pressure drop behind the disk. A higher value of \(C_{P}\) in the presence of free-stream turbulence has been also reported in previous experimental works [48].

It is worth noting that for \(a < 1/3\), the effect of turbulent mixing is the opposite, as seen in figures  9 (a) and (b). In this case, mixing slightly increases the near-wake expansion (see figure 4 (a)), which reduces \(|P_D^-|\) according to equation 25 and consequently lowers both \(C_T\) and \(C_P\). As already shown, however, the influence of turbulence mixing on the rotor performance for \(a < 1/3\) is much weaker than for \(a > 1/3\). In reality, the impact of wake turbulent mixing on rotor performance for small induction factors (\(a < 1/3\)) may be even smaller than shown here if entrainment does not begin immediately behind the disk; for instance, if shedding tip vortices act as a barrier between the wake and the ambient flow [49].

4 Summary and conclusions↩︎

In this work, we develop a generalised actuator disk theory that extends the classical formulation by: (i) describing how flow properties vary with the streamwise coordinate \(x\), rather than limiting the analysis to far-upstream and far-downstream conditions; (ii) incorporating the effects of turbulent mixing in the wake and its contribution to wake recovery; and (iii) providing physically consistent predictions of the thrust coefficient \(C_T\) for a given induction factor \(a\) (or vice-versa) for highly-loaded disks.

The mass and momentum equations are formulated for a hybrid CV encompassing the streamtube and the wake, upstream and downstream of the actuator disk, respectively. In the downwind region, wake recovery is modelled as the result of turbulent entrainment across the lateral surface of the CV with an entrainment velocity that depends on both the wake shear and the ambient turbulence. Moreover, it is shown that, at finite streamwise positions, pressure forces acting on the lateral surface of the CV contribute a non-negligible term to the axial momentum balance and must be included for a self-consistent formulation.

By integrating mass and momentum balances between far upstream and an arbitrary downstream position \(x\), we derive a generalised Bernoulli-type equation that accounts for both energy extraction by the disk and energy injection into the wake via turbulent mixing. An additional pressure equation is obtained by solving a simplified, axisymmetric form of the pressure Poisson equation, yielding accurate predictions for lightly-loaded disks when compared with LES data. Together, the system consists of three governing relations: two differential equations (obtained from mass and momentum equations) and one algebraic equation (obtained from an approximate solution of the pressure Poisson equation), which can be solved to obtain the streamwise evolution of velocity \(U(x)\), pressure \(P(x)\), and CV diameter \(\sigma(x)\) for given values of the induction factor \(a\) and thrust coefficient \(C_T\). The upwind region can be solved analytically, while the downwind region must be solved numerically using a forward marching scheme. The model predicts a physically consistent wake evolution, with an initial flow deceleration due to pressure recovery followed by re-acceleration from turbulent entrainment, along with a gradual wake expansion. Also, the asymptotic far-wake behaviour is examined. If entrainment is driven solely by wake shear, the wake width expands as \(x^{1/3}\). When background turbulence is also included, the wake expansion becomes linear with \(x\) (i.e., a constant wake expansion rate). However, for small values of the incoming turbulent intensity, the wake expands linearly only at sufficiently large downstream distances, and at intermediate downwind distances, the wake experiences a variable expansion rate.

The new theory provides a new framework to determine the \(C_T\)–\(a\) relationship. This relation is achieved by enforcing the condition that the axial component of the net pressure force on the lateral control surface asymptotically vanishes as \(x \to \infty\). A simple iterative method then yields \(U\), \(P\), \(\sigma\), and \(C_T\) as functions of \(a\). The results show that for small values of \(a\) or \(C_T\), the dependence between \(C_T\) and \(a\) is largely unaffected by wake turbulent mixing, whereas for larger values of \(a\) or \(C_T\), it becomes highly sensitive to the level of turbulent entrainment in the near-rotor downstream wake region. The maximum power coefficient can slightly exceeds the Betz limit, and we analyse the causes for this effect. Notably, the model reduces to the Froude classical actuator disk relationship between \(a\) and \(C_T\) when turbulence entrainment is neglected.

Future extensions of this framework may incorporate more realistic wake profiles to capture the transition from a top-hat shape in the near wake to Gaussian-like shapes in the far wake observed in experiments and high-fidelity simulations. Additionally, retaining non-linear terms in the pressure Poisson equation could improve predictions for highly-loaded actuator disks.

Acknowledgements↩︎

MB acknowledges the support of the EPSRC Impact Acceleration Account at Durham University. CM and DG acknowledge support from the National Science Foundation and the Department of Energy (via NSF grant CBET-2401013).

Declaration of Interests.↩︎

The authors report no conflict of interest.

5 Pressure Poisson equation for actuator disk flows↩︎

Let us consider an axisymmetric, steady, incompressible flow without swirl around an actuator disk of radius \(R = D/2\), subjected to an incoming flow with velocity \(U_0\). Note that in this Appendix, unlike the rest of the paper, our variables (e.g., \(P\) and \(U\)) depend on both streamwise \(x\) and radial \(r\) coordinates. The flow is described using cylindrical coordinates \((x, r)\), and \(U\) and \(V\) represent the time-averaged streamwise and radial velocity components, respectively. Taking the divergence of the momentum equation (including the actuator disk force and an eddy-viscosity term for turbulence) and combining with the continuity equation, yields the mean pressure Poisson equation which can written as

\[\label{eq:final95poisson32equation} \frac{1}{\rho}\nabla^2P=- \frac{1}{\rho}\Delta P_D \frac{\partial \delta(x)}{\partial x}H(R-r) + f_1(x,r)+f_2(x,r),\tag{34}\] where \[\nabla^2=\frac{\partial^2}{\partial x^2} + \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial}{\partial r}\right),\] and \(\Delta P_D=0.5\rho C_T U_0^2\) is the total pressure drop due to forcing at \(x=0\). The non-homogenous terms \(f_1(x,r)\) and \(f_2(x,r)\) are respectively defined as \[\label{eq:f40x44r41} f_1(x,r)=-\left(\frac{\partial U}{\partial x}\right)^2 - \frac{1}{r}\frac{\partial (rU)}{\partial r} \frac{\partial V}{\partial x},\tag{35}\] and \[\label{eq:g40x44r41} f_2(x,r)= \frac{\partial \nu_t}{\partial x} \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial U}{\partial r} \right)+\frac{\partial \nu_t}{\partial r}\frac{\partial}{\partial r}\left(\frac{1}{r}\frac{\partial \left(rV\right)}{\partial r}\right).\tag{36}\]

5.1 Neglecting non-homogenous terms in pressure Poisson equation↩︎

We here show that the terms \(f_1\) and \(f_2\) are small compared to the expected impact of the actuator disk force which must be balanced by the pressure Laplacian. The term \(f_1\) arises from the advective nonlinearity of the RANS equations (see equation 35 ), while the term \(f_2\) originates from spatial variations in turbulent viscosity (see equation 36 ). Therefore, removing \(f_1\) and \(f_2\) terms is equivalent to solving the linearised RANS equations with constant turbulent viscosity. To estimate the order of magnitude of each term, we introduce characteristic velocity and length scales: \(\mathcal{U}_0\) represents the magnitude of the inflow velocity, while \(\mathcal{U}_d\) characterises the velocity deficit. The cross-stream and streamwise length scales are denoted by \(\mathcal{L}_r\) and \(\mathcal{L}_x\), respectively. The following equation expresses the order of magnitude of each term in the pressure Poisson equation: \[\label{eq:PPE95dimension} \underbrace{\frac{1}{\rho}\frac{\partial^2 P}{\partial x^2}}_{\text{\large \frac{\mathcal{U}_0\mathcal{U}_d}{\mathcal{L}_x^2}}} + \underbrace{\frac{1}{\rho}\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial P}{\partial r}\right)}_{\text{\large \frac{\mathcal{U}_0\mathcal{U}_d}{\mathcal{L}_r^2}}} = - \frac{1}{\rho}\Delta P_D \frac{\partial \delta(x)}{\partial x}H(R-r) \, + \underbrace{f_1(x, r)}_{\text{\large \frac{\mathcal{U}_d^2}{\mathcal{L}_x^2}}} + \underbrace{f_2(x, r)}_{\text{\large \frac{\mathcal{U}_d^2}{\mathcal{L}_x \mathcal{L}_r}}}\tag{37}\] To estimate the pressure variation scale, we assume based on Bernoulli’s equation that \(\delta P/\rho\propto U_0 \, \delta U \sim \mathcal{U}_0\mathcal{U}_d\) when \(\mathcal{U}_d \ll \mathcal{U}_0\). For the turbulent viscosity \(\nu_t\) in \(f_2(x,r)\), we assume \(\mathcal{O}(\nu_t) = \mathcal{U}_d \mathcal{L}_r\) [50]. From Equation 37 , it follows that for lightly loaded disks (i.e., small values of \(a\)) where \(\mathcal{U}_d/\mathcal{U}_0 \to 0\), the terms \(f_1\) and \(f_2\) become negligible compared to the terms in the Laplacian, as a leading-order approximation. In this case, the two Laplacian terms on the left-hand side must balance each other away from the disk (\(\mathcal{L}_x \propto \mathcal{L}_r\)) and each still much larger than \(f_1\) and \(f_2\) near the disk. Therefore, while we present the complete solution of the pressure Poisson equation in the following for the sake of completeness, we omit the effects of \(f_1\) and \(f_2\) for simplicity in the final result.

5.2 Solving the Poisson equation in an infinite domain↩︎

First, we solve \[\label{qnyrvbha} \frac{1}{\rho}\nabla^2P=- \frac{1}{\rho}\Delta P_D \frac{\partial \delta(x)}{\partial x}H(R-r)\tag{38}\] for an infinite domain \(-\infty<x<\infty\) using the Greens function method. For our model we seek the solution \(P(x, r)\) evaluated at the centreline \(r = 0\), i.e., \(P(x, 0)\). The domain of the variables is \[x \in (-\infty, \infty), \quad r \in [0, \infty),\] with the following boundary conditions: \[\begin{align} &\lim_{x \to \pm\infty} P(x, r) = 0, \\ &\lim_{r \to \infty} P(x, r) = 0, \\ &\left. \frac{\partial P}{\partial r} \right|_{r = 0} = 0 \quad \text{(due to axisymmetry)}. \end{align}\] In Cartesian coordinates, the general solution of \((1/\rho)\nabla^2 p= q(x,y,z)\) on \(\mathbb{R}^3\), subject to zero Dirichlet conditions at infinity is: \[\label{infsol} \frac{1}{\rho}p(x,y,z)=\int_{x'=-\infty}^\infty\int_{y'=-\infty}^\infty\int_{z'=-\infty}^\infty \frac{-q(x',y',z')}{4\pi\sqrt{(x-x')^2+(y-y')^2+(z-z')^2}}\, \,\mathrm{d}z' \,\mathrm{d}y' \,\mathrm{d}x'.\tag{39}\] In polar coordinates, with \(P(x,r,\theta)=p(x,r\cos\theta,r\sin\theta)\), \(Q(x,r,\theta)=q(x,r\cos\theta,r\sin\theta)\), this amounts to \[\frac{1}{\rho}P(x,r,\theta)=\int_{x'=-\infty}^\infty\int_{r'=0}^\infty\int_{\theta'=0}^{2\pi} \frac{-\,r'Q(x',r',\theta')}{4\pi\sqrt{(x-x')^2+(r')^2-2rr'\cos(\theta-\theta')+r^2}}\, \,\mathrm{d}\theta' \,\mathrm{d}r' \,\mathrm{d}x'.\] If \(Q\) is axisymmetric, integrate over \(\theta'\) and set \(a=(x-x')^2+(r')^2+r^2,\;b=2rr',\;\phi=\theta'-\theta\) to get \[\begin{align} \int_{\theta'=0}^{2\pi}\frac{\,\mathrm{d}\theta'}{\sqrt{a-b\cos(\theta-\theta')}}=&\int_{\phi=-\theta}^{2\pi-\theta}\frac{\,\mathrm{d}\phi}{\sqrt{a-b\cos\phi}}=\int_{\phi=0}^{2\pi}\frac{\,\mathrm{d}\phi}{\sqrt{a-b\cos\phi}}=\int_{\psi=0}^{\pi}\frac{2\,\mathrm{d}\psi}{\sqrt{a+b-2b\cos^2\psi}}\\ =&\,\frac{4K(\sigma)}{\sqrt{a+b}}\,,\quad \text{where } \sigma=\sqrt{2b/(a+b)}. \end{align}\] Here \(K(\sigma)\) is the complete elliptic integral of the first kind. For computational purposes, the identity \[\label{kagm} K(\sigma)=\frac{\pi}{2 M(1,\sqrt{1-\sigma^2})}\tag{40}\] is useful, where \(M(m,n)\) is the arithmetic-geometric mean (which is computed iteratively, with quadratic convergence). This yields the axisymmetric solution \[\label{axiinfsol} \frac{1}{\rho}P(x,r)=\int_{x'=-\infty}^\infty\int_{r'=0}^\infty \frac{-\,r'Q(x',r')}{2\sqrt{(x-x')^2+(r+r')^2}\;\,M\!\left(1,\frac{\sqrt{(x-x')^2+(r-r')^2}}{\sqrt{(x-x')^2+(r+r')^2}}\right)}\, \,\mathrm{d}r' \,\mathrm{d}x'.\tag{41}\] On the axis of symmetry, 41 reduces to \[\label{axiinfcl} \frac{1}{\rho}P(x,0)=\int_{x'=-\infty}^\infty\int_{r'=0}^\infty \frac{-\,r'Q(x',r')}{2\sqrt{(x-x')^2+(r')^2}}\, \,\mathrm{d}r' \,\mathrm{d}x',\tag{42}\] which can also be obtained directly from 39 by setting \(y=z=0\) and then transforming to polar coordinates. Note that \(Q(x,r)\) must decay to zero as \(r\rightarrow\infty\) in order for \(P(x,0)\) to be finite.

By linear superposition, each summand in the function \(Q(x,r)\) can be integrated separately to give its contribution to the solution. For instance, taking \(Q(x,r)= -(1/\rho) \Delta P_D\, \delta'(x)\, H(R-r)\), where \(\delta'\) denotes the formal derivative of the delta function, gives \[\begin{align} \label{eq:infi95sol95f61g61095App} P(x,0)=&\int_{x'=-\infty}^\infty \Delta P_D\, \delta'(x')\int_{r'=0}^R \frac{r'}{2\sqrt{(x-x')^2+(r')^2}}\, \,\mathrm{d}r' \,\mathrm{d}x'\nonumber\\ =&\,\tfrac{1}{2}\,\Delta P_D\left(\frac{x}{\sqrt{x^2+R^2}}\,-\,\frac{x}{|x|}\right). \end{align}\tag{43}\] An important limitation of the infinite-domain solution is that for a first-order approximated solution where we assume \(f_1 = f_2 = 0\), equation 43 predicts a symmetric pressure discontinuity across the rotor disk, i.e., \(|P_D^-| = |P_D^+| = 0.5\Delta P_D\) at \(x = 0\), which is incorrect. As discussed in section 2.4, the pressure drop across the rotor is not symmetrical, and generally \(|P_D^-|<|P_D^+|\). This discrepancy indicates that the non-homogenous terms in the Poisson equation must be retained for this solution to capture the asymmetric pressure jump across the disk. In this regard, the semi-infinite domain solution (i.e., solving for \(x>0\) and \(x<0\) separately) may be more suitable, as it allows specification of pressure right upstream or downstream of the disk from equation 15 , offering a more realistic modelling framework even without considering the non-homogenous terms of \(f_1\) and \(f_2\).

5.3 Solving Laplace equation in two semi-infinite domains↩︎

We here consider the domain \(x > 0\), i.e. starting immediately after the disk, where its forcing has no contribution and can be eliminated. We again set \(f_1 =f_2 =0\)) but consider that the most important effect of \(f_1+f_2\) on the pressure distribution is approximated by setting the boundary condition to the pressure’s known value at the centerline immediately behind the disk, i.e., \(P_{D}^-\) as specified by equation 15 . This value is imposed as a Dirichlet boundary condition at \(x = 0\) and \(r<R\) for the Laplace equation. We thus solve: \[\label{eq:PPE95smiinf} \nabla^2P = 0,\tag{44}\] with boundary conditions: \[\begin{align} P(x = 0, r) &= P_D^- \, H(R - r), \\ \lim_{x \to \infty} P(x, r) &= 0, \\ \lim_{r \to \infty} P(x, r) &= 0, \\ \left. \frac{\partial P}{\partial r} \right|_{r = 0} &= 0 \quad \text{(axisymmetry)}. \end{align}\] where \(P_D^-<0\) is a constant. Again, we are interested in the solution along the axis of symmetry, i.e., \(P(x>0, r = 0)\). In Cartesian coordinates, the general solution of \(\nabla^2 p= q(x,y,z)\) on the half-space \(x\geq 0\), subject to \(p(0,y,z)=s(y,z)\) and zero Dirichlet conditions at infinity is: \[\begin{align} \frac{1}{\rho}p(x>0,y,z)=&\iiint \frac{q(x',y',z')}{4\pi}\left\{\frac{1}{\sqrt{(x+x')^2+(y-y')^2+(z-z')^2}}-\frac{1}{\sqrt{(x-x')^2+(y-y')^2+(z-z')^2}}\right\} \,\mathrm{d}x' \,\mathrm{d}y' \,\mathrm{d}z'\nonumber\\ &+\frac{1}{\rho}\iint\frac{xs(y',z')}{2\pi}\left\{x^2+(y-y')^2+(z-z')^2\right\}^{-3/2}\,\mathrm{d}y' \,\mathrm{d}z'.\label{hspoisscart} \end{align}\tag{45}\] In terms of polar coordinates, with \(S(r,\theta)=s(r\cos\theta,r\sin\theta)\) and the assumption of axisymmetry, this solution amounts to \[\begin{align} \frac{1}{\rho}P(x>0,r)=&\int_{x'=0}^\infty\int_{r'=0}^\infty \frac{r'Q(x',r')\,K\!\left(\frac{2\sqrt{r'r}}{\sqrt{(x+x')^2+(r+r')^2}}\right)}{\pi\sqrt{(x+x')^2+(r+r')^2}}-\frac{r'Q(x',r')\,K\!\left(\frac{2\sqrt{r'r}}{\sqrt{(x-x')^2+(r+r')^2}}\right)}{\pi\sqrt{(x-x')^2+(r+r')^2}}\, \,\mathrm{d}r' \,\mathrm{d}x'\nonumber\\[4pt] &+\frac{1}{\rho}\int_{r'=0}^\infty\frac{2x\,r'S(r')\,E\!\left(\frac{2\sqrt{r'r}}{\sqrt{x^2+(r+r')^2}}\right)}{\pi\sqrt{x^2+(r+r')^2}\,\{x^2+(r-r')^2\}}\, \,\mathrm{d}r', \end{align}\] where \(E\) is the complete elliptic integral of the second kind. Here, if we assume the sum of non-homogenous terms in equation 44 is zero (i.e. \(Q(x,r)=f_1(x,r)+f_2(x,r)=0\)), the solution simplifies to \[\label{eq:P40x44r41} P(x>0,r)=\int_{r'=0}^\infty\frac{2x\,r'S(r')\,E\!\left(\frac{2\sqrt{r'r}}{\sqrt{x^2+(r+r')^2}}\right)}{\pi\sqrt{x^2+(r+r')^2}\,\{x^2+(r-r')^2\}}\, \,\mathrm{d}r'.\tag{46}\] Figure 10 shows the variation of \(P(x,r)\), computed from equation 46 , at several streamwise locations downstream of the disk. Immediately behind the disk, where the pressure is relatively high, the radial pressure profiles exhibit a top-hat shape. Further downstream, the pressure drops sharply, and the profiles become increasingly smooth. For simplicity, here we use the pressure along the axis of symmetry as a representative measure of how pressure evolves with streamwise position. This provides the necessary closure relation for the actuator disk theory developed in the main part of the paper. On the axis of symmetry, equation 46 simplifies to \[\label{eq:PPE95solution95semiinf95f61g61095x620} P(x>0,0)=\int_{r'=0}^R\frac{P_D^{-}x\,r'}{\{x^2+(r')^2\}^{3/2}} \,\mathrm{d}r'=P_D^{-}\left(1-\frac{x}{\sqrt{x^2+R^2}}\right).\tag{47}\]

As an alternate method of solving the Laplace equation, the same result can be also obtained using the Hankel transform method. For the \(x>0\) domain, this leads to \(P(x,r) = P_D^{+} \, R \int_0^\infty J_1(\rho R) \,J_0(\rho r) \, e^{-\rho x} \, \,\mathrm{d}\rho\), where \(J_0\) and \(J_1\) are zero and first order Bessel functions. Evaluating at \(r=0\) and performing the integral of the exponentials-weighted \(J_1\) function yields the same result as in equation 47 .

Equation 47 can also be applied to the upwind region by substituting \(-x\) for \(x\), and using \(P_D^+\) as the corresponding boundary condition. Thus, the solution for \(P(x, 0)\) over the entire domain \(-\infty < x < \infty\), \(x\neq 0\), assuming \(f_1=f_2=0\), is given by: \[\label{eq:PPE95solution95semiinf95f61g61095app} P(x,0)=P_D^{+}\left(1+\frac{x}{\sqrt{x^2+R^2}}\right)\mathrm{H}(-x)+P_D^{-}\left(1-\frac{x}{\sqrt{x^2+R^2}}\right)\mathrm{H}(x),\tag{48}\]

We have replaced the actual effects of \(f_1+f_2\) with specification of the distinct values of \(P_D^+\) and \(P_D^-\) that are obtained instead from an “integral” condition (integrated Bernoulli equation) of the problem. Finally, it is worth noting that the simplified solution presented here does not capture the complexities of the pressure distribution around the disk or outside the control volume, and should be regarded only as an approximation of pressure variation along the disk axis.

References↩︎