2.1. A nonlinear model for turbine rotation rate and power extraction

We build upon the linear modelling approach of Wei & Dabiri (Reference Wei and Dabiri2022), who describe the time-varying dynamics of a turbine using the swing equation (i.e. Newton's second law for rotation)

(2.1)\begin{equation} J\,\frac{{\rm d}\omega}{{\rm d}t} = \tau_{aero} - \tau_{gen}, \end{equation}

where $J$ is the moment of inertia of the turbine system about its axis of rotation, and $\omega$ is the rotation rate of the turbine. Deviations from equilibrium between the aerodynamic and generator torques $\tau _{aero}$ and $\tau _{gen}$ produce changes in the turbine rotation rate. Wei & Dabiri (Reference Wei and Dabiri2022) model the generator torque using the ordinary differential equation for the torque from a permanent-magnet motor,

(2.2)\begin{equation} \tau_{gen} = K_2\,\frac{{\rm d}\omega}{{\rm d}t} + K_1 \omega + K_0, \end{equation}

where $K_2$ is the moment of inertia of the generator about its axis of rotation, $K_1$ is the generator constant, and $K_0$ is an empirical offset. These parameters can be established empirically for a given generator over a range of resistive loads.

While Wei & Dabiri (Reference Wei and Dabiri2022) utilized a linearized model for the aerodynamic torque, in this work we include explicitly the nonlinear relationship between the turbine coefficient of power,

(2.3)\begin{equation} C_p = \frac{\mathcal{P}}{\frac{1}{2}\rho {\rm \pi}R^2 u_\infty^3}, \end{equation}

and the tip-speed ratio,

(2.4)\begin{equation} \lambda = \frac{R\omega}{u_\infty}, \end{equation}

where $\mathcal {P}$ is the power extracted by the turbine from the flow, $\rho$ is the fluid density, $R$ is the radius of the turbine, and $u_\infty$ is the free-stream flow velocity. Any given turbine has a power curve defined as $C_p = C_{p,0}(\lambda )$, which has a local maximum at a power-maximizing tip-speed ratio $\lambda _{opt}$. Note that in this work, we will use the subscript $0$ to refer to steady-flow quantities. Since the turbine power is determined by the product of the torque on the turbine and its rotation rate, the torque on the turbine can be written in terms of the power curve as

(2.5)\begin{equation} \tau = \frac{1}{2}\,\rho {\rm \pi}R^2 u_\infty^3\,\frac{C_{p,0}(\lambda)}{\omega}. \end{equation}

For a given surge-velocity profile $U(t)$, the effective free-stream velocity is $u_\infty = u_1 - U(t)$, where we define $u_1$ as the far-field wind speed relative to a ground-fixed frame. Thus we may write the aerodynamic torque as

(2.6)\begin{equation} \tau_{aero} = \tau_{aero}(\omega,U,t) = \frac{1}{2}\,\rho {\rm \pi}R^2 \frac{\left(u_1-U(t)\right)^3}{\omega}\,C_{p,0}\left(\frac{R\omega}{u_1-U(t)}\right). \end{equation}

Substituting (2.2) and (2.6) into (2.1) results in a nonlinear ordinary differential equation for the turbine rotation rate that is first-order in time and depends on the surge velocity as an input forcing parameter:

(2.7)\begin{equation} \frac{{\rm d}\omega}{{\rm d}t} = \frac{1}{J+K_2}\left[\frac{1}{2}\,\rho {\rm \pi}R^2 \frac{\left(u_1-U(t)\right)^3}{\omega}\,C_{p,0}\left(\frac{R\omega}{u_1-U(t)}\right) - K_1 \omega - K_0\right]. \end{equation}

This model is a nonlinear and non-autonomous ordinary differential equation, which precludes straightforward mathematical characterization, but it can be integrated forward in time from an initial condition $\omega (t=0)$ until it reaches a period-averaged equilibrium. The model can therefore yield numerical predictions of the time-varying and time-averaged rotation rate, torque and power of a turbine under surge motions or dynamic inflow conditions. This stands in contrast to the linearized model developed by Wei & Dabiri (Reference Wei and Dabiri2022), which can be written as a transfer function for convenient analysis but is unable to capture changes in time-averaged quantities.

As a limiting case, we may consider a quasi-steady solution to the model in which ${{\rm d}\omega }/{{\rm d}t}=0$ and $C_{p,0}$ is constant as a function of time. The time-averaged power is defined as $\bar {\mathcal {P}} = \overline {\tau _{gen}\omega } = K_1 \overline {\omega ^2} + K_0 \bar {\omega }$, which suggests that for a sinusoidal surge-velocity waveform with amplitude $u^* = \hat {U}/u_1$, the time-averaged power is

(2.8)\begin{equation} \bar{\mathcal{P}} = \frac{f}{2{\rm \pi}} \int_0^{{2{\rm \pi}}/{f}} \frac{1}{2}\,\rho {\rm \pi}R^2 C_{p,0} \left(u_1 - U(t)\right)^3 {\rm d}t = \mathcal{P}_0\left(1+\frac{3}{2}\,{u^*}^2\right). \end{equation}

This result is identical to that derived by Wen et al. (Reference Wen, Tian, Dong, Peng and Zhang2017) and Johlas et al. (Reference Johlas, Martínez-Tossas, Churchfield, Lackner and Schmidt2021) for a surging turbine, and is equivalent to that of a stationary turbine with constant $C_p$ in an oscillating inflow.

2.2. Modelling the relationship between turbine dynamics and upstream flow conditions

The 1-D momentum theory derived by Betz (Reference Betz1920) can be used to infer flow properties upstream of the turbine rotor plane from the power extracted by the turbine. This theoretical framework employs conservation relations over a control volume composed of an axisymmetric streamtube surrounding the turbine, which is modelled as an actuator disc. An axial induction factor, representing the decrease in velocity from far upstream of the turbine to the upstream face of the actuator disc (i.e. location 2 in figure 1), is defined as

(2.9)\begin{equation} a = \frac{u_1 - u_2}{u_1}. \end{equation}

The induction factor is the single free parameter needed to compute the coefficient of power within this framework, which is given by

(2.10)\begin{equation} C_p = 4a(1-a)^2. \end{equation}

This yields a theoretical maximum for the efficiency of a wind-energy system, $C_{p,Betz} = 16/27$, which is attained at $a=1/3$. For the surging-turbine system, then, we have the similar relation

(2.11)\begin{equation} \frac{\tau_{gen}\omega}{\frac{1}{2}\rho {\rm \pi}R^2 u_1^3} = 4a(1-a)^2. \end{equation}

Figure 1. Schematic of the parameters and control volumes referenced in § 2. The actuator disc is located instantaneously at $x_2(t)$ and moves with velocity $U(t)$ relative to the inertial frame defined by the $x$- and $r$-axes. Circled numbers denote streamwise interrogation locations (1–4).

One-dimensional momentum theory assumes that the flow is inviscid, incompressible, irrotational and steady. Upstream of the turbine rotor plane, the first three assumptions are reasonable to make. If we further assume that the flow upstream of the turbine can be modelled in a quasi-steady manner, then we can use the rotation rate given by (2.7) to compute an instantaneous coefficient of power $C_p(\omega,t)$, and by inverting (2.11), we can estimate the time-varying induction factor $a(t)$. Since the induction factor is physically constrained as $a\in [0,1]$, two solutions are possible for $C_p < C_{p,Betz}$. The upper solution for $a$ represents a heavily loaded turbine, and for $a\gtrsim 0.37$, the theoretical framework breaks down (Wilson & Lissaman Reference Wilson and Lissaman1974). In a majority of cases, including the experiments presented in this work, the turbine is not heavily loaded. Thus the lower solution for $a$ is assumed to represent the system.

Equations (2.7) and (2.11) therefore connect the time-varying dynamics of a turbine under dynamic axial flow conditions to the flow properties just upstream of the rotor disc. To propagate these predictions further upstream, an induction model for the turbine is needed. As mentioned previously, a common modelling approach involves representing the wake of the turbine as a cylindrical vortex sheet and performing Biot–Savart integration to compute the induced velocity from this wake model at any point in the domain (Johnson Reference Johnson1980; Branlard & Gaunaa Reference Branlard and Gaunaa2015). This is known as vortex-sheet or vortex-cylinder theory (VCT). Evaluating this integral upstream of the turbine along its rotational axis yields a model for the induced velocity along the upstream centreline of the turbine (Medici et al. Reference Medici, Ivanell, Dahlberg and Alfredsson2011):

(2.12)\begin{equation} \frac{u(r=0,x)}{u_\infty} = 1 - a\left[1+\frac{x-x_2}{R}\left(1+\left(\frac{x-x_2}{R}\right)^2\right)^{{-}1/2}\right], \end{equation}

where $x$ is the streamwise coordinate along the axis of the turbine (originating at the turbine and positive downstream), and $x_2$ is the instantaneous location of the turbine, as shown in figure 1. For a surging turbine, the induction effect in the second term should scale with the effective free-stream velocity $u_1-U(t)$, since in the limiting case where the turbine is moving downstream at velocity $U(t) = u_1$, the turbine should have no effect on the flow, and the flow everywhere upstream of the turbine should be equal to $u_1$. Involving the time-varying induction factor $a(t)$ provided by 1-D momentum theory, we thus obtain the flow velocity at any upstream location $x< x_2$ and at any time as

(2.13)\begin{align} u(r=0,x,t) = u_1 - a(t)\,(u_1-U(t))\left[1+\frac{x-x_2(t)}{R}\left(1+\left(\frac{x-x_2(t)}{R}\right)^2\right)^{{-}1/2}\right]. \end{align}

To complete the description of flow properties upstream of the surging turbine, the pressure along the centreline may be modelled by substituting the model for $u(r=0,x,t)$ into the steady Bernoulli equation,

(2.14)\begin{equation} p(r=0,x) = p_1 + \frac{1}{2}\,\rho (u_1^2 - {u(r=0,x,t)}^2), \end{equation}

where $\rho$ is the density of the fluid, and $p_1$ is the ambient pressure in the free stream. If the velocity potential of the induction model is known, then the unsteady Bernoulli equation may be applied instead. For the purposes of this work, however, we maintain quasi-steady assumptions for the flow physics in the induction zone, in keeping with the quasi-steady aerodynamics parametrized by the turbine model and 1-D momentum theory. As mentioned previously, the effect of this unsteady potential will be considered in more detail in § 5.2. While the expressions presented here have been confined to the centreline, Branlard & Gaunaa (Reference Branlard and Gaunaa2015) provide full relations for $u(r,x)$ that can be employed in place of (2.12) to allow this modelling framework to cover the entire upstream induction zone.

Alternatively, we can model the effect of the surging turbine on the flow using a porous-disc representation in potential flow, which does not rely on parametrizations of the turbine wake geometry and vorticity. This is inspired by the work of Taylor (Reference Taylor1944) and Koo & James (Reference Koo and James1973), which has recently been extended by Steiros & Hultmark (Reference Steiros and Hultmark2018) for flat plates of arbitrary porosity, and by Bempedelis & Steiros (Reference Bempedelis and Steiros2022) for wind turbines at arbitrary loading conditions. The porous-disc approach, which has not been investigated widely in the wind-energy literature, is presented to demonstrate the robustness of the overall modelling framework put forth in this work to the choice of induction model. Its velocity potential is also more readily accessible than that of VCT, which will be advantageous when we consider the effects of unsteady fluid mechanics in § 5.2. Finally, the model provides a convenient generalization to rotors with arbitrary radial distributions of streamwise velocity, which, though not explored in detail in this study, could be exploited to integrate this model with blade element momentum (BEM) computations that generate radially varying induction-factor profiles. We will refer to the model throughout this paper as the porous-disc theory (PDT).

2.3. A porous-disc induction model for a surging turbine

We first consider a circular porous disc with radius $R$ located at streamwise coordinate $\xi =0$, represented as a distribution of sources with velocity potential $\phi (r,\xi =0) = C\sqrt {R^2-r^2}$ for $r< R$, and $\phi (r,\xi =0)=0$ for $r>R$, where $C>0$ is an arbitrary constant that represents the aggregate strength of the source distribution. Using this distribution as a boundary condition at $\xi =0$, we may solve the Laplace equation $\nabla ^2\phi =0$ in cylindrical coordinates to obtain the velocity potential of a porous disc (cf. Lamb Reference Lamb1916; Tranter Reference Tranter1968):

(2.15) \begin{equation} \phi\left(r,\xi\right) ={-}\sqrt{\frac{\rm \pi}{2}}\,C R^{3/2} \int_0^\infty s^{-1/2}\,{\rm J}_{3/2}(Rs)\,{\rm J}_0(rs)\,{\rm e}^{{-}s \xi} \,{\rm d}s,\quad \xi>0. \end{equation}

Here, ${\rm J}_\nu (z)$ is a Bessel function of the first kind, and $s$ is a dummy integration variable. The choice $C=({2}/{{\rm \pi} })V$ gives the velocity potential of a solid disc moving at axial velocity $V$ in a quiescent fluid (Lamb Reference Lamb1916, § 102.4). More generally, the velocity $V$ represents the velocity of the disc relative to that of the fluid in the far field. For a porous disc, we may define a representative source term $a$, corresponding directly to the induction factor defined in (2.9), such that $C=({2}/{{\rm \pi} })a V$. The choice of $a$ dictates the porosity of the disc: $a=0$ represents a fully permeable disc, $a=1$ yields a fully solid disc, and intermediate values ($0< a<1$) reduce the source strength from the solid disc solution so that a non-zero mass flux through the disc is established. Evaluating (2.15) along the centreline yields

(2.16)\begin{equation} \phi\left(r=0,\xi\right) ={-}a V\,\frac{2}{\rm \pi}\left[R - \xi\arctan\left(\frac{R}{\xi}\right)\right],\quad \xi>0. \end{equation}

This solution is valid only for $\xi >0$. To describe the other half of the domain as well, one might follow the ansatz of Taylor (Reference Taylor1944) and use the even extension of $\phi$ to represent $\xi <0$. The velocity discontinuity across the disc that this extension creates could then be removed using the base suction correction of Steiros & Hultmark (Reference Steiros and Hultmark2018). However, since in this work we are concerned only with the upstream region, we leave these derivations for future consideration.

We differentiate the velocity potential with respect to $\xi$ to obtain the streamwise velocity along the centreline:

(2.17)\begin{equation} u\left(r=0,\xi\right) ={-}a V\, \frac{2}{\rm \pi}\left[\frac{R\xi}{{\xi}^2+R^2}-\arctan\left(\frac{R}{\xi}\right)\right],\quad \xi > 0. \end{equation}

This relation emphasizes the effect of the porosity parameter (or equivalently, the induction factor) on the behaviour of the model. For $a=0$, the flow is everywhere unaffected by the motion of the disc. For $a=1$, $u(r=0,\xi \rightarrow 0^+)=V$, which satisfies the surface boundary condition for a moving solid disc.

We now apply this expression in an inertial frame containing a uniform flow with free-stream velocity $u_1$, in which the disc translates at velocity $U(t)$ relative to the frame. In this frame, we define the downstream-oriented axial coordinate $x$ and the instantaneous position of the disc $x_2(t)$ as shown in figure 1, such that $\xi =x_2 - x$ and $U(t) = {{\rm d}\kern0.06em x_2}/{{\rm d}t}$. The velocity of the disc relative to the far-field flow velocity is thus $V = U(t) - u_1$. Applying these definitions to (2.17), we arrive at the following model for the centreline velocity in the upstream induction zone ($x< x_2$):

(2.18)\begin{align} u(r=0,x,t) = u_1 - a(t) \left(u_1-U(t)\right) \frac{2}{\rm \pi}\left[\frac{R(x-x_2(t))}{\left(x-x_2(t)\right)^2+R^2} -\arctan\left(\frac{R}{x-x_2(t)}\right)\right]. \end{align}

We reiterate that the model cannot be used to predict the velocity downstream of the porous disc, given the constraint $\xi >0$ on the velocity potential. Additionally, though the solution is technically valid if the turbine moves downstream faster than the wind speed, i.e. $u_1-U(t) < 0$, we expect that the model will cease to be valid in this case because the rotor will interact with its own wake.

As with the expression obtained from VCT, (2.15) can in principle be integrated at any point upstream of the porous disc so that the velocity and pressure can be computed throughout the entire induction zone. We note that in both models, the streamwise velocity $u_2$ on the upstream face of the rotor plane is predicted to be constant over $r$, i.e.

(2.19)\begin{equation} u(x=x_2,r)=u_1-a(t)\,(u_1-U(t)). \end{equation}

This is generally a poor approximation for real turbines, whose streamwise velocities typically increase towards the free-stream value with increasing radial distance from the hub (cf. Medici et al. Reference Medici, Ivanell, Dahlberg and Alfredsson2011; Troldborg & Meyer Forsting Reference Troldborg and Meyer Forsting2017). An additional advantage of the porous-disc modelling approach is that arbitrary radial velocity profiles can be modelled by changing the source strength distribution in (2.15). For example, if the turbine blade geometry is known, then radial variations in the induction factor can be computed from BEM theory as a function $a(r)$ and integrated to obtain a modified velocity potential. A correction factor for the effects of non-uniform velocities at the disc face can then be derived by defining a rotor-averaged induction factor,

(2.20)\begin{equation} {\langle a \rangle} = \frac{2}{R^2}\int_0^R a(r)\,r \,{\rm d}r = \left[\frac{2}{R^2}\int_0^R \frac{a(r)}{a(r=0)}\,r \,{\rm d}r\right] a(r=0) \equiv \kappa\,a(r=0). \end{equation}

The scaling constant $\kappa$, which maps the centreline induction factor to the equivalent rotor-averaged induction factor, can be computed from a known velocity-deficit radial profile, such as the self-similar solution of Troldborg & Meyer Forsting (Reference Troldborg and Meyer Forsting2017) or a BEM computation. For a top-hat velocity-deficit profile, $\kappa =1$. If the radial induction-factor distribution does not change much throughout the surge cycle (i.e. the velocity-deficit profile remains self-similar), then $\kappa$ can be assumed to be independent of the surge kinematics for a given loading condition, and thus can be treated as a constant for that particular wind speed and mean tip-speed ratio.

The distinction between the centreline and rotor-averaged induction factors is particularly important for a quantitative comparison between the two-dimensional (2-D) axisymmetric induction theories developed in this section and the 1-D axial momentum theory of Betz. Because 1-D momentum theory by definition does not account for radial differences in streamwise velocity, the induction factor estimated using (2.11) must represent the rotor-averaged induction factor. By the same logic, the induction factor in the VCT and PDT expressions for the streamwise velocity is the centreline induction factor, and is what will be measured by a point anemometer placed along the rotational axis of the turbine. The parameter $\kappa$ thus serves as a bridge between the 1-D and 2-D analyses employed in this modelling framework.

For the remainder of this paper, we will focus on measurements taken along the rotational axis of the turbine ($r=0$), therefore the flow quantities $u(x)$, $p(x)$ and $a$ will implicitly refer to these centreline quantities.

The modelling approach presented in this work may be summarized as follows.

1. (i) Identify the turbine power curve $C_{p,0}(\lambda )$, and generator and inertia constants $K_0$, $K_1$, $K_2$ and $J$, from steady-flow measurements and manufacturer specifications.

2. (ii) Integrate (2.7) in time with a given surge-velocity forcing $U(t)$ to obtain the turbine rotation rate $\omega (t)$.

3. (iii) From $\omega (t)$, use (2.2) and (2.11) to predict the coefficient of power $C_p(t)$ and the rotor-averaged induction factor $\langle a(t) \rangle$.

4. (iv) Calculate the centreline induction factor $a(t)$ using an analytical or empirical correction factor $\kappa$.

5. (v) Include $a(t)$ in an induction model (e.g. (2.13) for VCT, or (2.18) for PDT) to obtain the velocity field upstream of the turbine.

6. (vi) Use the steady Bernoulli equation to obtain the pressure field from the velocity field.

We reiterate that this modelling framework invokes a quasi-steady assumption for the aerodynamics of the turbine, and the time dependence of the model comes from an unsteady treatment of the turbine rotation rate as a function of the aerodynamic and generator torques. The practical benefit of this approach is that time-resolved predictions of the turbine dynamics and upstream flow properties can be obtained solely on the basis of steady-flow measurements; the model does not depend on empirical calibrations from unsteady test cases. From a fluid mechanics perspective, the approach provides an instructive disambiguation between the rotational dynamics of the turbine and the actuator-disc aerodynamics associated with the rotor and its motions. Furthermore, this analytical foundation allows the effects of unsteady flow physics to be characterized more directly.

2.4. On the problem of a stationary turbine in an oscillating inflow

The preceding analysis has focused on the case of a periodically surging turbine in a uniform inflow. From the work of Wen et al. (Reference Wen, Tian, Dong, Peng and Zhang2017), El Makdah et al. (Reference El Makdah, Ruzzante, Zhang and Rival2019), Johlas et al. (Reference Johlas, Martínez-Tossas, Churchfield, Lackner and Schmidt2021) and others, it is apparent that the case of a stationary turbine in an axial gust can be made equivalent to the surging-turbine case by shifting from a ground-fixed to a turbine-fixed frame of reference. The time-averaged power is not affected by this transformation. To determine the effect of non-inertial-frame accelerations, we consider the force on a body oscillating with velocity $W_i(t)$ in an oscillating inflow $U_i(t)$, which is given by Brennen (Reference Brennen1982) as

(2.21)\begin{equation} F_i ={-}M_{ij}\,\frac{{\rm d}W_j}{{\rm d}t} + \left(M_{ij} + \rho V_D \delta_{ij}\right) \frac{{\rm d}U_j}{{\rm d}t};\quad j = 1,2,3. \end{equation}

Here, $M_{ij}$ is the added-mass tensor of the body, $V_D$ is the volume of the body, $\delta _{ij}$ is the Kronecker delta operator, and the flow is assumed to be inviscid. This expression thus quantifies the influence of added-mass effects and an unsteady buoyancy force, which comes from the oscillating pressure gradient that drives the oscillating inflow (Granlund et al. Reference Granlund, Monnier, Ol and Williams2014). We assume that neither the added-mass tensor nor the volume of the body changes as a function of time for a porous disc, and that $U_i(t)$ and $W_i(t)$ are periodic. Since $U_i(t)$ and $W_i(t)$ are periodic, the time averages over a single period of ${{\rm d}U_j}/{{\rm d}t}$ and ${{\rm d}W_j}/{{\rm d}t}$ are both zero. It thus follows that the time-averaged force on the body due to these two types of unsteady contributions is also zero, and therefore neither of these unsteady effects creates a theoretical difference between the time-averaged performance of an oscillating turbine and a stationary turbine in an oscillating inflow. We also note that the volume of a 2-D actuator disc is effectively zero, thus the unsteady buoyancy force that differentiates the two scenarios should be negligible. These considerations do not exclude the possibility of additional differences between these two cases, for example due to viscous time-history effects, vortex dynamics, or other flow physics not captured by the potential-flow assumption.

3.1. Experimental apparatus

To characterize the range of conditions over which the ideal flow model holds, velocity and pressure measurements were conducted in a $2.88\times 2.88\,{\rm m}^2$ open-circuit fan-array wind tunnel at the Caltech Center for Autonomous Systems and Technologies (CAST). A three-bladed, fixed-pitch horizontal-axis wind turbine (Primus Wind Power AIR Silent X) with rotor diameter $D = 1.17$ m was mounted on a traverse that translated along 2 m long rails (NSK NH-series) and was actuated by a magnetic piston-type linear actuator (LinMot PS10-70x320U). A diagram of this apparatus is given in figure 2. The hub height of the turbine was 1.97 m above the floor of the facility, and the farthest downstream position of the turbine (defined as $x=0$) was 3.09 m downstream of the fan array. The electrical load on the turbine was provided by 10, 20 and 40 $\Omega$ resistors (TE Connectivity TE1000-series). A rotary torque transducer (FUTEK TRS300) and rotary encoder (US Digital EM2) were used to measure the power produced by the turbine. The estimated blockage of the swept area of the turbine and all support structures, relative to the surface area of the fan array, was 14 %. Further details regarding the dimensions and capabilities of the apparatus may be found in Wei & Dabiri (Reference Wei and Dabiri2022).

Figure 2. Schematic of the experimental apparatus, including the fan-array wind tunnel (left) and surging turbine (centre-right). The turbine is illustrated at its maximum upstream position relative to the origin ($x = -0.6$ m). The inset (top right) shows the two types of surge-velocity waveforms used in these experiments.

A constant-temperature hot-wire anemometry system (Dantec MiniCTA 54T42) and differential pressure transducer (MKS Baratron 398-series with Type 270B signal conditioner) were used to measure flow properties at two locations along the turbine centreline, one upstream of the turbine at $x_u = -0.840 D$, and one downstream at $x_d = 0.810 D$. The hot-wire probe was placed approximately on the centreline, while the input line of the pressure transducer was located 3.8 cm to the side. The transducer's reference line was placed in a shielded area outside the flow of the wind tunnel. Data were collected at sampling rate 20 kHz and were low-pass filtered using a sixth-order Butterworth filter with cutoff frequency 100 Hz. The hot-wire anemometer was calibrated in the wind tunnel against a Pitot probe using the same pressure transducer. Because the facility was exposed to the atmosphere, the temperature and relative humidity were recorded during all experiments to estimate the air density and correct the hot-wire calibration for temperature changes.

3.2. Experimental procedure

Experiments were conducted over two nights in March 2022, in which the free-stream velocities in the wind tunnel were $u_1 = 7.79\pm 0.10$ and $7.96\pm 0.11\,{\rm m\,s}^{-1}$, corresponding to an average diameter-based Reynolds number $Re_D = 6.27\times 10^5$. The hot-wire anemometer was calibrated at the beginning and end of each set of experiments. The turbine was operated at three tip-speed ratios, $\lambda _0 = 6.48\pm 0.25$, $7.84\pm 0.28$ and $8.77\pm 0.27$, with corresponding coefficients of power $C_{p,0} = 0.298\pm 0.013\approx C_{p,max}$, $C_{p,0}=0.248\pm 0.012$ and $C_{p,0}=0.165\pm 0.010$. The generator constants for these cases were obtained from the data of Wei & Dabiri (Reference Wei and Dabiri2022), as were the turbine and generator moments of inertia ($J=0.0266\pm 0.0008\,{\rm kg}\,{\rm m}^2$, $K_2 = 6.96\times 10^{-4}\,{\rm kg}\,{\rm m}^2$). The operating parameters of the turbine are summarized in table 1. The turbine was actuated in sinusoidal and trapezoidal motions (see inset of figure 2) with amplitude $A = 0.3$ m ($0.257D$) and periods between $T = 1$ s and $T=6$ s, corresponding to non-dimensional surge-velocity amplitudes between $u^* \equiv fA/u_1 = 0.039$ and $0.242$. Data were phase-averaged over 100 motion periods. The amplitudes and phases of each quantity of interest were computed from a fast Fourier transform of the phase-averaged signal. Upstream and downstream flow measurements were collected in separate tests. Additionally, a series of quasi-steady flow measurements was obtained for each tip-speed ratio by placing the turbine at six equally spaced streamwise locations between $x/D = -0.514$ and $0$, and recording measurements over 120 s. To correct against differences in the ambient conditions between measurement sessions, and facilitate more direct comparisons, quasi-steady measurements taken at $x/D = 0$ on both sessions were used to scale the measured velocities and pressures from one session to match those from the other session.

Table 1. Performance characteristics and model constants for the three loading conditions investigated in this study.

3.3. Determination of model parameters from steady-flow measurements

3.3.1. Power-curve parametrization and model integration

To compare the analytical modelling framework with the experimental data, it was first necessary to parametrize the steady-flow power curve of the turbine, $C_{p,0}(\lambda )$. The data in the measured power curve of the turbine, shown in figure 3, were fitted to a type of exponential function used frequently for wind turbine modelling:

(3.1)\begin{equation} C_p = \left(\frac{c_1}{\lambda+c_2}-c_3\right)\exp{\left(\frac{-c_4}{\lambda+c_2}\right)}. \end{equation}

Several coefficients from the general model form given by Heier (Reference Heier2014) were omitted to simplify the model; the four remaining fitted coefficients were $c_1=16.784$, $c_2=-1.510$, $c_3=1.702$ and $c_4=8.764$. This parametrization was employed (as opposed to e.g. polynomial fits) to ensure that the slope and concavity at the extremes of the power curve would be captured reliably, since it will be shown in § 4.1 that the performance of the modelling framework is sensitive to these factors.

Figure 3. Steady power curve for the turbine used in these experiments, measured over a range of resistive loads and wind speeds. Some of these data points are reproduced from Wei & Dabiri (Reference Wei and Dabiri2022). The result of the exponential fit given by (3.1) is shown as a solid red line.

To obtain time-resolved predictions of the turbine rotation rate, torque and power, (2.7) was integrated numerically over ten surge periods using a fourth-order Runge–Kutta scheme. Time steps were kept no larger than $0.001T$ to maintain numerical stability and accuracy. The steady-flow turbine rotation rate $\omega _0$ was used as the initial condition, and convergence was established typically within a few forcing periods. The model predictions for the amplitude, phase and time average of each quantity were computed from the final period in the simulation.

3.3.2. Steady-flow induction data

To estimate the steady induction-factor values of the turbine at the three tip-speed ratios tested in this study, quasi-steady measurements of the streamwise velocity $u(x)$ (described in § 3.2) were plotted as a function of streamwise distance $x/D$, and two-parameter fits for the wind speed $u_1$ and centreline induction factor $a_0$ were applied to these data. Dimensional data for the highest tip-speed ratio tested are shown in figure 4 with fits using the VCT (2.13) and PDT (2.18) models. This test case demonstrates the slight differences between the modelling frameworks: for the same data, the PDT model predicts a stronger induction effect (as $x/D \rightarrow 0$) and a slightly lower free-stream velocity ($x/D \leq -2$) than the VCT model. The fit results for all three tip-speed ratios are shown in figures 5(a) (VCT) and 5(b) (PDT), where the measured and modelled flow velocities are normalized by the wind speeds identified by each two-parameter fit. These steady-flow tests demonstrate that within the range of streamwise distances tested, the agreement of both models with the trends observed in the data is reasonably good.

Figure 4. Streamwise velocity measurements at six different streamwise distances upstream of the turbine for the highest tip-speed ratio tested, compared with two-parameter fits based on the VCT (dotted line) and PDT (dash-dotted line) induction models. The PDT model shows a slightly more aggressive drop in streamwise velocity close to the turbine and a lower predicted free-stream velocity than the VCT model.

Figure 5. Streamwise velocity measurements at six different streamwise distances upstream of the turbine and three tip-speed ratios, compared with fits based on (a) the VCT and (b) the PDT induction models. The velocity data are normalized by the free-stream velocities obtained from the two-parameter fits.

The centreline induction-factor values shown in figure 5 were compared with the rotor-averaged values estimated using (2.11) and a $C_p$ calculated from the average of the measured turbine torque and rotation rate over all six streamwise locations. The ratio between these estimates of $\langle a_0 \rangle$ and the fitted values of $a_0$ gave an empirical estimate for the correction factor $\kappa$ for each tip-speed ratio. The estimated values of $a_0$ and $\kappa$ are given in table 1 for the two induction model frameworks. For comparison, an analytical correction factor can be calculated by integrating the self-similar solution of Troldborg & Meyer Forsting (Reference Troldborg and Meyer Forsting2017), which yields $\kappa =0.649$. This is larger than the values of $\kappa$ estimated empirically. However, the ratio of the turbine hub diameter to the turbine diameter for the simulations used to calibrate the self-similar solution was 2.4 %, whereas for the turbine used in this study, the ratio was around 12 %. It is therefore expected that the flow would decelerate more strongly along the centreline of this particular turbine, due to the increased blockage effect from the larger hub and nacelle, thus lowering the ratio between the rotor-averaged and centreline induction factors.

4.1. Power generation

The predictions of the nonlinear ordinary differential equation given in (2.7) are compared with rotation rate, torque and power data collected from the surging-turbine experimental apparatus. For this purpose, we reuse the experimental torque and power measurements of Wei & Dabiri (Reference Wei and Dabiri2022), and plot the results of the present modelling framework against these data. This dataset contained measurements at three lower tip-speed ratios, in addition to those investigated in the present study, and measurements were conducted over a wider range of surge-velocity amplitudes and frequencies.

Figure 6 shows the measured torque amplitude and phase from these data, plotted against surge frequency. The torque amplitude was scaled by the surge-velocity amplitude $u^*=fA$ in the manner of a transfer-function gain, and this was non-dimensionalized by the steady reference torque $\tau _{gen,0}$ and the free-stream velocity. The frequency was non-dimensionalized by a characteristic frequency $f_c$ derived from the linear model of Wei & Dabiri (Reference Wei and Dabiri2022). The trends in the data are relatively well captured by the model predictions, and the nonlinear model shows some improvement over the linear model at the lowest tip-speed ratios tested (cf. Wei & Dabiri Reference Wei and Dabiri2022, figures 7a and 9a). This suggests that the nonlinear model derived in this work is an effective generalization of the linearization developed and validated in the preceding study.

Figure 6. Generator torque (a) amplitude and (b) phase, for a series of sinusoidal surge-velocity waveforms, plotted against normalized frequency and compared with model predictions (dashed lines). Data are reproduced from Wei & Dabiri (Reference Wei and Dabiri2022).

The benefit of the nonlinear model is evident when predictions for time-averaged quantities of interest are required. According to the preceding linearized model, the time-averaged rotation rate and power of the turbine will not deviate from their corresponding steady-flow quantities. Departures from this ansatz are clearly evident in figure 7, as the mean rotation rate and power decrease as a function of increasing surge-velocity amplitude for low tip-speed ratios, and increase for high tip-speed ratios. The predictions of the nonlinear model, however, are able to follow these trends. The largest enhancements in time-averaged power over the steady case are predicted at the highest tip-speed ratios, while the greatest decreases in time-averaged power occur at the lowest tip-speed ratios. The nonlinear model still overpredicts these power decreases; it is likely that these additional losses are a result of the inception of flow separation on the turbine blades. The turbine stalled and stopped spinning when forced to operate at or below $\lambda _0\approx 5$, so the decrease in power as a function of decreasing tip-speed ratio was in reality much sharper than that suggested by the fit to the turbine power curve in figure 3. The lack of a parametrization for these dynamics in the current modelling framework is thus probably responsible for the lack of quantitative agreement at lower tip-speed ratios.

Figure 7. Time-averaged (a) rotation rate and (b) power, plotted against surge-velocity amplitude. Data are reproduced from Wei & Dabiri (Reference Wei and Dabiri2022); model predictions derived from (2.7) are plotted as dashed lines and coloured by tip-speed ratio. Circles represent sinusoidal surge-velocity waveforms, and diamonds and other markers represent trapezoidal waveforms. Error bars are plotted on every sixth point for the sake of clarity.

The discussion of flow separation and the extent to which stall is captured in the turbine power curve emphasizes the key point that the critical nonlinearity in (2.7) is the functional form of the power curve itself. In other words, the present quasi-steady modelling framework hinges on a reliable characterization of the turbine in steady flow conditions. This is a considerable advantage of the modelling approach, since it precludes the need for unsteady calibration and can thus be applied directly to the design of turbines in unsteady flow conditions when only steady-flow data are available. It also implies, however, that particular attention must be paid to the parametrization of the steady-flow power curve of the turbine. This dependence and its implications will be discussed in § 5.1.

4.2. Upstream flow properties

We now investigate the extension of the nonlinear model for the turbine dynamics to the flow properties in the upstream induction region of the turbine. The unsteady and quasi-steady data from three selected experimental cases, all measured at $x=x_u$, are shown in figure 8. The streamwise velocity signals showed a phase lead and increased amplitude with respect to the quasi-steady measurements. Also shown in this figure are the VCT and PDT model predictions, which align well with the shape of the phase-averaged data and anticipate the increased amplitudes and phase leads as well. For all of these cases, both models show good agreement with the velocity and pressure data, as well as with each other.

Figure 8. Phase-averaged velocity and pressure profiles for (a,b) a sinusoidal surge-velocity waveform with $\lambda _0=6.48$, (c,d) a trapezoidal waveform with $\lambda _0=6.48$, and (e, f) a trapezoidal waveform with $\lambda _0=8.77$. All surge-velocity waveforms had $u^* = 0.242$. The solid red lines represent unsteady measurements, the blue squares represent quasi-steady (QS) measurements, the dotted lines show the VCT model predictions, and the dash-dotted lines show the PDT model predictions.

The differences between the induction predictions of the models is reflected in the induction factors estimated from the power data. The time-averaged induction factors estimated using the VCT and PDT models and correction factors are shown in figures 9(a) and 9(b), respectively. The induction factors all increase with increasing surge-velocity amplitude, but at different rates depending on tip-speed ratio. The trends are consistent between the two models; the main difference is that the estimated induction factors from the PDT model are slightly higher, in accordance with the model's sharper induction profile (noted previously in figure 4).

Figure 9. Time-averaged induction-factor estimates from (a) the VCT model and (b) the PDT model. For these and all following figures, circles indicate sinusoidal waveform data, and diamonds denote trapezoidal waveform data. Points are estimated from measured power data, while lines are estimated from the nonlinear dynamical model of the turbine.

The increases in time-averaged induction factor with surge-velocity amplitude are tied to increases in the thrust force exerted by the turbine on the incoming flow. According to 1-D momentum theory, the thrust coefficient of the turbine is given as

(4.1)\begin{equation} C_t \equiv \frac{F_t}{\frac{1}{2}\rho {\rm \pi}R^2 u_\infty^2}= 4{\langle a \rangle}\left(1-{\langle a \rangle}\right), \end{equation}

which, for the relatively low values of $ {\langle a \rangle }$ considered in this study, increases with increasing $ {\langle a \rangle }$. The induction-factor estimates in figure 9 therefore suggest that the thrust force on the turbine increases with increasing surge-velocity amplitude.

This aligns with a simplified analysis of the thrust force of a turbine in an oscillating inflow. If $C_t$ is held constant and the dimensional thrust is integrated over a sinusoidal surge-velocity waveform, as done for $C_p$ in (2.8), then the time-averaged thrust enhancement is

(4.2)\begin{equation} \frac{\overline{F_t}}{F_{t,0}} = 1 + \frac{1}{2}\,{u^*}^2, \end{equation}

which is also an increasing function with $u^*$. Although direct thrust measurements were not possible in this study due to the large inertial forces associated with the turbine motions, the induction-factor estimates and constant-$C_t$ analysis both suggest that unsteady surge motions increase the time-averaged thrust loading on the turbine rotor. In accordance with the trends observed in the time-averaged power data in figure 7(b), the time-averaged thrust may decrease with increasing surge-velocity amplitudes for tip-speed ratios lower than those investigated here.

The case studies in figure 8 and induction-factor estimates in figure 9 suggest that the VCT and PDT models produce very similar results. For the sake of clarity, figures 10–12 will show only PDT predictions. The VCT predictions are qualitatively similar, and are provided for completeness in Appendix A.

Figure 10. (a) Amplitude and (b) phase of the estimated induction factors using the PDT model, plotted against surge-velocity amplitude. Model predictions are given as dash-dotted lines.

Figure 11. (a) Amplitude and (b) phase of the measured flow velocity at $x=x_u$, plotted against surge-velocity amplitude. PDT model predictions are given as dash-dotted lines. Error bars are plotted on every fourth point for clarity.

Figure 12. (a) Amplitude and (b) phase of the measured relative pressure at $x=x_u$, plotted against surge-velocity amplitude. PDT model predictions are given as dash-dotted lines. Error bars are plotted on every fourth point.

In figure 10, the amplitude and phase (relative to the surge-velocity waveform) of the estimated induction factors are presented as a function of surge-velocity amplitude. The induction-factor amplitudes and phase offsets decrease in magnitude with increasing tip-speed ratio; the highest amplitudes observed in the measured data do not exceed 12 % of the time-averaged values. Good agreement between the estimates from the measured data and the estimates from the nonlinear dynamical model is observed.

Good agreement between the model predictions and measured data is also observed in the streamwise velocity data, shown in terms of amplitude and phase (relative to the quasi-steady measurements) in figure 11. As noted previously, the surge motions slightly increase the amplitude of the velocity oscillations above the quasi-steady case, and a slight phase lead accrues with increasing surge-velocity amplitude. There is some quantitative mismatch between the phase measurements and model predictions, including a plateau in the phase data at the highest surge-velocity amplitude that is not reflected by the model. Given the large experimental uncertainties in the streamwise velocity measurements, this is not surprising. However, the slopes of the model prediction lines still align with the slopes of the measured data (except for the highest surge-velocity amplitude tested), which differ slightly across tip-speed ratios. This suggests that the modelling framework is able to capture some of the more subtle differences in flow properties in the upstream induction region as the tip-speed ratio of the turbine changes.

Similar agreement in trend is observed in the pressure data, shown in amplitude and phase (relative to the quasi-steady measurements) in figure 12. The slight increase in pressure amplitude as a function of surge-velocity amplitude is generally reflected in the model predictions, and as in the streamwise velocity phase data, the model predictions of the pressure phase follow slopes similar to those evident in the data. Again, some quantitative differences are apparent in the plots, but given the measurement uncertainties and the relatively small magnitude of the signals being quantified and predicted, the qualitative agreement between the model predictions and measured data suggests that the modelling framework is parametrizing the salient dynamics of the system.

Finally, the time-averaged difference between the upstream and downstream pressure measurements at $x=x_u$ and $x=x_d$, respectively, is shown as a function of surge-velocity amplitude in figure 13 for both the VCT and PDT models. The corresponding model predictions are calculated by evaluating the model at $x=x_u$ and $x=-x_d$, and assuming that the pressure exhibits odd symmetry about the turbine rotor plane so that $p(x=-x_d) = -p(x=x_d)$, as is the case for 1-D momentum theory. For a real turbine with a low-pressure wake region, this is likely not a tenable assumption. However, the model predictions still follow the trends in the data within measurement uncertainty. The measured pressure difference can be understood as an analogue to the thrust force on the turbine, and as implied previously by the time-averaged induction-factor estimates, it increases with increasing surge-velocity amplitude. As well as further demonstrating the predictive capabilities of the modelling framework, these data suggest that the model could also be extended to serve as an initial condition for wake models of the surging turbine, which will depend on the thrust force and pressure drop across the rotor.

Figure 13. Pressure difference across the turbine, $\Delta p = p(x=x_u)-p(x=x_d)$, normalized by the quasi-steady pressure difference and plotted against surge-velocity amplitude. Model predictions are given as dotted lines for the VCT model in (a), and as dash-dotted lines for the PDT model in (b). Error bars are shown on every fourth point.

In summary, both induction model frameworks are able to reproduce the trends observed in the centreline flow measurements recorded upstream of the turbine, and the modelling approach may have some bearing on the near-wake region downstream of the turbine as well. This observed agreement implies that the proposed modelling framework is capturing the dominant dynamics of the surging-turbine system, in spite of incorporating only empirical data from the turbine power curve and the induction-factor correction parameter $\kappa$ – both of which can be obtained from steady-flow theories, simulations or experiments. While the nonlinear dynamical model accounts for unsteady rotation accelerations of the turbine, the aerodynamics of the turbine itself is treated in a quasi-steady manner, as a function of the power curve and the induction models. The fact that this quasi-steady first principles model still manages to align with data collected from a real surging turbine might suggest that unsteady flow physics is not strongly present in this system. However, the uncertainties in these experiments preclude us from ruling out this possibility, and we will thus present a theoretical analysis of the effects of unsteady flows in § 5.2.

5.1. Time-averaged power enhancements

The results presented in § 4.1 suggested that the amount of enhancement in the time-averaged power extraction of a surging turbine relative to the steady case depends on the characteristics of the power curve of the turbine. To explore this connection further, a parametric survey using simulations of the turbine model from (2.7) with simplified power curves was instantiated. For four known equilibrium points along the power curve, corresponding to four of the tip-speed ratios tested in experiments, local quadratic power curves were constructed with the form

(5.1)\begin{equation} C_p(\lambda) = C_p(\lambda_0) + \left.\frac{{\rm d}C_p}{{\rm d}\lambda}\right\vert_{\lambda_0} (\lambda-\lambda_0) + \frac{1}{2}\,\frac{{\rm d}^2C_p}{{\rm d}\lambda^2}\,(\lambda-\lambda_0)^2, \end{equation}

where the local slope ${{\rm d}C_p}/{{\rm d}\lambda }\vert _{\lambda _0}$ was taken from the fitted power curve in (3.1), and the concavity ${{\rm d}^2C_p}/{{\rm d}\lambda ^2}$ was varied between $-0.1$ and 0.02. The four tip-speed ratios used in this exploration were $\lambda _0=5.47$, 6.33, 7.67 and 8.69, and the local slope of the second tip-speed ratio was approximated as zero. This parametrization created a set of power curves for each equilibrium point with identical local slope but different concavities, and the effect of changing local slope could be ascertained by comparing the simulation results across tip-speed ratios. A fixed surge-velocity amplitude of $u^*=0.24$ was used for the simulations; other relevant parameters were $u_1 = 8\,{\rm m\,s}^{-1}$, $\rho =1.19\,{\rm kg}\,{\rm m}^{-3}$ and $\Delta t = 0.001$ s. As before, a fourth-order Runge–Kutta scheme was used to integrate the model forward in time over ten surge-oscillation periods. A strict convergence metric required the difference between successive periods to decrease monotonically; test cases that failed this criterion were deemed unstable and were not plotted.

The results of these simulations are shown in figure 14, where the relative power enhancement $\overline {C_p}/C_{p,0} - 1$ is plotted as a function of power curve concavity for the four selected tip-speed ratios. Filled circles show the concavity of the full power curve from (3.1), corresponding approximately to the power gains and losses observed in figure 7(b) at $u^*\approx 0.24$. At zero slope and zero concavity, the system reduces to the quasi-steady prediction derived by Wen et al. (Reference Wen, Tian, Dong, Peng and Zhang2017) and Johlas et al. (Reference Johlas, Martínez-Tossas, Churchfield, Lackner and Schmidt2021) (i.e. (2.8), shown as a red $\times$ in the figure). Where the concavity is zero, differences are still evident across local slopes: the two highest tip-speed ratios have negative local slopes and larger power enhancements relative to the constant-$C_p$ case. Conversely, the lowest tip-speed ratio (which has a positive local slope) exhibits the lowest power enhancement. As concavity decreases from zero, these power gains decrease and eventually become losses, finally becoming unstable below some critical concavity. Positive concavities, by contrast, show increasing power enhancements.

Figure 14. Fractional enhancements in the time-averaged coefficient of power, plotted against the concavity of the power curve (approximated as a quadratic function with fixed local slope and variable concavity), for four tip-speed ratios. A sinusoidal surge-velocity waveform with $u^*=0.24$ was used for these simulations. The red $\times$ shows the constant-$C_p$ solution given in (2.8), and the solid circles denote the approximate concavity of the actual power curve at each reference tip-speed ratio.

The influence of power curve concavity can be understood geometrically. The dynamics of the turbine is constrained in the present model to its power curve, and thus periodic forcings represent periodic excursions along the power curve centred at some equilibrium point $\lambda _0$. First, let us assume that the local slope at $\lambda _0$ is zero. If the curve is concave down in a neighbourhood about $\lambda _0$, then the value of $C_p$ at $\lambda _0\pm \epsilon$ will be lower than $C_p(\lambda _0)$. Therefore, the integrated value over a periodic excursion away from $\lambda _0$ will be lower than the equilibrium value at $\lambda _0$. The opposite is true when the power curve is concave up; since $C_p(\lambda _0)$ is now a local minimum, non-zero perturbations away from $\lambda _0$ will lead to an increased time-averaged $C_p$. These arguments also hold qualitatively for a non-zero local slope at $\lambda _0$, though the local slope does have an influence on the time-averaged value across a periodic perturbation.

This analysis underscores the point that the geometry of the steady power curve and the equilibrium operating point of the turbine dictate the time-averaged power enhancements or losses that a surging turbine will experience relative to the steady-flow case. In general, it will be favourable to operate in regions of a power curve that have minimal concavity. In terms of turbine design, this implies that turbines whose power curves exhibit a relatively flat maximal region will benefit the most from time-averaged power enhancements relative to the steady case. Flattening the curve in this manner may be achieved through traditional turbine design methods; alternatively, the pitch of the turbine blades could be varied within a single oscillation period to produce the same topological effect on the $C_p$ manifold. Finally, while concave-up power curves are not found typically in current wind-energy systems, this analysis does suggest that such designs (if physically possible) would lead to even greater power enhancements in unsteady flow conditions.

5.2. The role of unsteady flow physics

Until this point, we have neglected the contribution of unsteady fluid mechanics to the system in question. We now investigate these effects analytically. The analysis in this subsection is purely theoretical, and is included to complete the conceptual picture of wind-energy systems in dynamic-inflow conditions that has been presented in this work. Thus quantitative predictions and comparisons with the experimental results shown previously are not pursued.

5.2.1. Specifying velocity potentials for an unsteady extension to 1-D momentum theory

To characterize the influence of streamwise unsteadiness on the theoretical efficiency of wind-energy systems, we extend the analytical framework of Dabiri (Reference Dabiri2020) using the potential-flow modelling approach proposed in § 2.3. This unsteady extension to the 1-D momentum theory of Betz uses the same control volume shown in figure 1, but makes use of two additional unsteady terms:

(5.2)\begin{equation} \frac{{\rm d}}{{\rm d}t}\left[KE\right] = \frac{{\rm d}}{{\rm d}t}\left[-\frac{1}{2}\,\rho \iint_{A_2} \phi \hat{n}\boldsymbol{\cdot} \boldsymbol{\nabla}\phi \,{\rm d}A\right], \end{equation}

which represents the unsteady power associated with changes in the streamwise kinetic energy of the actuator disc, and

(5.3)\begin{equation} \varPhi_t = \frac{\partial \phi_2}{\partial t} - \frac{\partial \phi_3}{\partial t}, \end{equation}

which is the difference in the unsteady potential across the face of the actuator disc. These quantities are accounted for in energy conservation relations to obtain expressions for the time-averaged unsteady power. Note that we use $A_i$ to refer to the cross-sectional area at streamwise location $i$; this is not to be confused with the surge-motion amplitude $A$ used elsewhere in this work.

By specifying a velocity potential for the moving actuator disc, we aim to close these two terms and parametrize the instantaneous and time-averaged coefficients of power as a function of the induction factor (as in the steady Betz analysis) and the surge kinematics of the actuator disc. For the sake of simplicity, we will assume that the induction factor $a$ varies in a quasi-steady manner, and will not consider time derivatives of $a$. Additionally, we consider only the velocity potential associated with the surge velocity $U(t)$, and assume that the velocity potential connected to the free-stream velocity $u_1$ has no unsteady contribution. We begin by completing the analysis for a moving porous disc, and then explore the effects of more general classes of velocity potentials. For the sake of brevity, we omit intermediate steps in the derivations here; several of these details are provided in Appendix B.

From (2.15), the velocity potential for a moving porous disc located instantaneously at $x=x_2$ can be written as

(5.4) \begin{equation} \phi(r,x) = a(u_1-U)\,\sqrt{\frac{2}{\rm \pi}}\,R^{3/2} \int_0^\infty s^{-1/2}\,{\rm J}_{3/2}(Rs)\,{\rm J}_0(rs)\,{\rm e}^{s (x-x_2)} \,{\rm d}s,\quad x< x_2, \end{equation}

and from this expression, the kinetic energy of the disc moving at velocity $U$ can be derived as

(5.5)\begin{equation} KE_{disc} = \frac{4}{3}\,\rho a^2 U^2 R^3 \end{equation}

(cf. Lamb Reference Lamb1916, Art. 102, Eq. (20)). The time derivative is therefore

(5.6)\begin{equation} \frac{{\rm d}}{{\rm d}t}[KE_{disc}] = \frac{8}{3}\,\rho R^3 a^2 U\,\frac{{\rm d}U}{{\rm d}t}. \end{equation}

The unsteady-potential term $\varPhi _t$ can be computed by applying the chain rule to (5.4), which gives

(5.7)\begin{equation} \frac{\partial \phi_2}{\partial t} ={-}\frac{2}{\rm \pi}\,a\,\sqrt{R^2-r^2}\,\frac{{\rm d}U}{{\rm d}t} + aU^2. \end{equation}

Using the odd extension to model the region immediately downstream of the disc, we obtain an expression for ${\partial \phi _3}/{\partial t}$ that is identical except that the first term is positive. Averaging this over the area of the disc gives

(5.8)\begin{equation} \langle \varPhi_t \rangle \equiv \varPhi_{t,disc} ={-}\frac{8R}{3{\rm \pi}}\,a\,\frac{{\rm d}U}{{\rm d}t}. \end{equation}

These results may now be introduced into the framework of Dabiri (Reference Dabiri2020), with some additional considerations that reduce the parameter space of the theory. The unsteady power associated with the motion of the disc leads to a difference between the available power upstream and downstream of the disc (cf. Dabiri Reference Dabiri2020, Eq. (11)):

(5.9)\begin{equation} \frac{1}{2}\,\rho A_2(u_2^3-u_3^3) = \frac{{\rm d}}{{\rm d}t}[KE]. \end{equation}

Defining additional induction factors $b = 1-u_3/u_1$ and $c = 1-u_4/u_1$, this expression can be written as

(5.10)\begin{equation} (1-a)^3-(1-b)^3=\frac{2}{\rho A_2 u_1^3}\,\frac{{\rm d}}{{\rm d}t}[KE], \end{equation}

which we can solve for $b$ in terms of $a$:

(5.11)\begin{equation} b = 1-\left[(1-a)^3-\frac{2}{\rho A_2 u_1^3}\,\frac{{\rm d}}{{\rm d}t}[KE]\right]^{1/3}. \end{equation}

A non-dimensional form of the momentum equation for this problem can be written by means of the unsteady Bernoulli equation as

(5.12)\begin{equation} -2c(1-a) = (1-c)^2 + (1-a)^2 - (1-b)^2 - 1 + \frac{2\varPhi_t}{u_1^2}, \end{equation}

which can be solved for the remaining induction factor $c$ as

(5.13)\begin{equation} c = a \pm \sqrt{2a - 1 + (1-b)^2 - \frac{2\varPhi_t}{u_1^2}}. \end{equation}

In this expression, the larger root is taken to satisfy the physical requirement that the flow must slow down between locations 2 and 4, which implies $c > a$. Having written the additional induction factors $b$ and $c$ in terms of the original induction factor $a$ and the known unsteady contributions $({{\rm d}}/{{\rm d}t})[KE]$ and $\varPhi _t$, we can find the instantaneous coefficient of power using the relation

(5.14) \begin{equation} C_p = \frac{1}{2}(4c-4c^2+c^3) + \frac{1}{2}\,(2-c)[(1-b)^2-(1-a)^2] - (2-c)\left(\frac{\varPhi_t}{u_1^2}\right). \end{equation}

This expression is limited by the physical constraints $b\in [0,1]$, $c\in [0,1]$ and $C_p \geq 0$, which enforce that the velocities at locations 3 and 4 cannot be negative or exceed the free-stream velocity, and that the power extracted from the actuator disc cannot be negative. Note that this analysis does not account for quasi-steady changes in $C_p$ due to the normalization by the effective inflow velocity $u_1 - U(t)$; therefore, changes in $\overline {C_p}$ predicted by this theoretical framework will appear in practice as adjustments to the time-averaged power predictions of the quasi-steady modelling approach outlined in § 2.1.

5.2.2. Phase-plane analysis of a surging porous disc

While Dabiri (Reference Dabiri2020) assumed that the parameters $a$, $b$, $c$ and $\varPhi _t$ were independent, we now have a theoretical framework that depends only on the induction factor $a$ and the surge kinematics of the actuator disc, $U$ and ${{\rm d}U}/{{\rm d}t}$. For a given induction factor, we may therefore use the system of equations described above to construct a phase portrait for $C_p$ in terms of the surge kinematics. Contours of $C_p$ in the $U$–${{\rm d}U}/{{\rm d}t}$ phase plane are shown for four values of $a$ (0.21, 0.27, $1/3$ and 0.40) in figure 15. The orange contour marks the Betz efficiency ($C_p/C_{p,Betz}=1$), and grey regions represent locations in the phase plane where one of the physical constraints on $C_p$ is violated.

Figure 15. Contours of $C_p/C_{p,Betz}$ calculated from (5.14) for a moving porous disc, plotted in the surge kinematics phase plane as a function of surge velocity and acceleration for induction factors (a) 0.21, (b) 0.27, (c) $1/3$, and (d) 0.40. The orange line indicates $C_p=C_{p,Betz}$, and grey regions denote dynamics that violate one or more assumptions of the theoretical framework. To non-dimensionalize the surge acceleration, $T$ is taken to be 1 s.

If $a$ is constant, then a surge waveform will appear as a periodic loop in the phase plane that must be centred on the origin so that the turbine has no net displacement. A sinusoidal waveform will follow an elliptical trajectory about the origin, while a trapezoidal waveform will appear as a rectangular trajectory. The time-averaged coefficient of power is calculated by evaluating the line integral of $C_p$ along one closed cycle of this trajectory. The depiction of $C_p$ in the phase plane thus allows the effects of the unsteady velocity potential and power terms from the porous disc to be evaluated by topological reasoning.

To illustrate this line of argument, we consider the case where $a=1/3$ (figure 15c). A trajectory that is centred on the origin will experience instantaneous values of $C_p$ in excess of the steady-flow Betz limit for positive surge accelerations, while $C_p$ will decrease below $C_{p,Betz}$ for negative surge accelerations. However, the slope of the contours of $C_p$ is greater below the Betz-limit contour than above it. Thus for a non-zero surge trajectory, the lower $C_p$ values sampled below the zero-acceleration axis will outweigh the higher $C_p$ values sampled above the axis, and the time-averaged coefficient of power will be lower than the steady coefficient of power. This mathematical effect is qualitatively similar to the concavity-based arguments described in § 5.1.

For values of $a$ that are above or below the steady-flow optimal value $1/3$, the same arguments hold: the topology of $C_p$ in the phase plane implies that non-trivial surge motions will yield a decrease in time-averaged efficiency as a result of unsteady effects. Furthermore, if $a$ is allowed to vary in a quasi-steady manner, the surge trajectory will exist in a three-dimensional phase space spanned by $U$, ${{\rm d}U}/{{\rm d}t}$ and $a$. Ascertaining the precise differences in time-averaged power extraction may be done numerically for a given profile of $a(t)$, e.g. one obtained from the amplitude and phase data shown in figure 10. Still, as the region in the phase plane where efficiencies above the steady Betz limit occur moves away from the origin for $a\neq 1/3$, it is apparent that oscillations in $a$ will further decrease the time-averaged efficiency of the system.

5.2.3. General velocity potentials and the effects of fore–aft asymmetry

We have shown that the unsteady contribution to the time-averaged efficiency of a wind-energy system modelled as a surging porous disc is negative. However, efficiency enhancements may be attained with a different choice of velocity potential. Consider a general velocity potential for a moving body with translation velocity $\boldsymbol {U}$, angular velocity $\boldsymbol {\varOmega }$, and circulation $\varGamma$, located instantaneously at a point $\boldsymbol {x}=\boldsymbol {x_0}$ in an otherwise quiescent fluid:

(5.15)\begin{equation} \phi(\boldsymbol{x})=\boldsymbol{U}\boldsymbol{\cdot}\boldsymbol{\varXi}+\boldsymbol{\varOmega}\boldsymbol{\cdot}\boldsymbol{\varTheta}+\varGamma\left(\frac{\theta}{2{\rm \pi}}+\psi\right) \end{equation}

(Batchelor Reference Batchelor2000, Eq. (6.4.10)). Here, boldfaced variables refer to vectors in 3-space. The functions $\boldsymbol {\varXi }$, $\boldsymbol {\varTheta }$ and $\psi$ are geometric descriptors that depend only on location relative to the body, $\boldsymbol {x}-\boldsymbol {x_0}$. For bodies producing zero circulation, $\varGamma = 0$. We also ignore the effect of rotation, and set $\boldsymbol {\varOmega }=0$. (The effect of rotation is explained in Appendix B.) The translational kinetic energy associated with such a moving body scales as $U^2$, thus the unsteady power term will always scale as

(5.16)\begin{equation} \frac{{\rm d}}{{\rm d}t}[KE]\sim U\,\frac{{\rm d}U}{{\rm d}t}, \end{equation}

irrespective of the body geometry parametrized by the function $\boldsymbol {\varXi }$.

The time derivative of the velocity potential, on the other hand, is affected by body geometry. This can be written as

(5.17)\begin{equation} \frac{\partial \phi}{\partial t} = \frac{{\rm d}\boldsymbol{U}}{{\rm d}t}\boldsymbol{\cdot} \boldsymbol{\varXi} - \boldsymbol{U}\boldsymbol{\cdot}\boldsymbol{u} \end{equation}

(Batchelor Reference Batchelor2000, Eq. (6.4.22)), where $\boldsymbol {u}$ is the local flow velocity at $\boldsymbol {x}$. To compute $\varPhi _t$ from the unsteady theoretical framework, we note that the local flow velocity in the quiescent fluid frame at streamwise locations 2 and 3 scales with the velocity of the body, $U$. Additionally, since the velocity on either side of the translating body must decrease with increasing distance from the body, the geometric function $\boldsymbol {\varXi }$ must exhibit odd symmetry about the body plane.

For a translating symmetric body in potential flow, $\boldsymbol {u}$ on either side of the body will be identical. Therefore, when taking the difference of the unsteady-potential terms upstream and downstream of the body, the $\boldsymbol {U}\boldsymbol {\cdot }\boldsymbol {u}$ term will cancel. The odd symmetry of the geometric function $\boldsymbol {\varXi }$ will retain the acceleration-dependent term, thus for a symmetric body, we have

(5.18)\begin{equation} \varPhi_t \sim R\,\frac{{\rm d}U}{{\rm d}t}. \end{equation}

For a body that exhibits fore–aft asymmetry about its flow-normal centre plane, however, $\boldsymbol {u}$ may differ across the upstream and downstream sides of the body. In this case, the quadratic velocity term may not cancel, suggesting that

(5.19)\begin{equation} \varPhi_t \sim R\,\frac{{\rm d}U}{{\rm d}t} \pm U^2. \end{equation}

The additional dependence on $U^2$ will change the topology of $C_p$ in the phase plane and affect the time-averaged coefficients of power that are possible in the theoretical framework.

To illustrate these effects, we conduct a phase-plane analysis by assuming the existence of a velocity potential for a moving asymmetric body that yields

(5.20)\begin{equation} \varPhi_t ={-}\frac{8}{3{\rm \pi}}\,a\left( R\,\frac{{\rm d}U}{{\rm d}t} + U^2\right). \end{equation}

We keep the same scaling coefficients as for the surging porous disc for the sake of comparison, and use the same expression for the kinetic energy of the body, (5.6). Contours of $C_p$ for the same four values of $a$ shown previously are given for this representative model in figure 16. We again first consider the case where $a=1/3$ (figure 16c). Here, we observe that the addition of the $U^2$ term in the expression for $\varPhi _t$ has curved the contours of $C_p$ such that the region where $C_p>C_{p,Betz}$ now extends below the zero-acceleration axis. This implies that a trajectory centred on the origin can have a time-averaged efficiency that exceeds the steady Betz limit. This may also be possible if we allow for small oscillations in $a$, depending on the surge kinematics applied and their corresponding trajectories in the phase space.

Figure 16. Contours of $C_p/C_{p,Betz}$ for a moving porous asymmetric body (parametrized by (5.20)), plotted in the surge-kinematics phase plane as a function of surge velocity and acceleration for induction factors (a) 0.21, (b) 0.27, (c) $1/3$, and (d) 0.40. The orange line indicates $C_p=C_{p,Betz}$, and grey regions denote dynamics that violate one or more assumptions of the theoretical framework.

It is important to note that if the sign of the $U^2$ term in (5.20) were reversed, then the concavity of the contours would also be reversed, leaving the region where $C_p>C_{p,Betz}$ on the interior of a concave-up parabola with respect to the surge velocity. In this scenario, a periodic trajectory will spend less time in the efficiency-enhancing region relative to the efficiency-depleting region, and the time-averaged efficiency will therefore drop below the steady Betz limit. Thus fore–aft asymmetry in the velocity potential is a necessary but not sufficient condition for enhancements in the time-averaged efficiency of a periodically translating actuator body.

This conceptual exercise demonstrates that for symmetric bodies in potential flow, the unsteady contribution to time-averaged efficiency will be negative. As the modelling efforts and experiments presented previously seem to suggest that a moving porous-disc model captures the dominant dynamics of a periodically surging horizontal-axis wind turbine, it is likely that unsteady effects would act against the time-averaged power enhancements described in § 5.1. This hypothesis appears to be consistent with the time-averaged power data shown in figure 7(b), which at high surge-velocity amplitudes tend to be overpredicted by the nonlinear turbine model (which assumes quasi-steady aerodynamics). However, wind-energy systems need not be symmetric. Streamwise asymmetries across the rotor plane could be introduced either mechanically, through the geometry of the turbine, or dynamically, through intracycle blade-pitch or generator-load control. These non-traditional design and control paradigms could create beneficial asymmetries in the equivalent velocity potential of the system that could be leveraged to achieve higher time-averaged efficiencies than steady-flow or quasi-steady analyses would suggest. In the case of dynamic induction control schemes, the influence of unsteady induction-factor variations ${{\rm d}a}/{{\rm d}t}$ may no longer be negligible, and the analysis presented here would need to be expanded to include these variations in a four-dimensional phase space. Still, the present analysis can serve as a helpful theoretical framework for characterizing trends in unsteady contributions to the efficiency of wind-energy systems.