2.1. Introduction to modelling framework

Figure 1 shows the control volumes (CVs) that will be used for the unsteady AD analysis; figure 1(a) (CV1) will be used for the momentum balance, figure 1(b) (CV2) for the energy balance and figure 1(c) (CV3) for the mass balance. Figure 1(c) also illustrates the definition of the flow acceleration parameters $\alpha _2$ and $\alpha _4$, such that the mean velocity at the aerofoil is given by $\alpha _2U_\infty$ and at the exit face by $\alpha _4U_\infty$, where $U_\infty$ is the free-stream velocity. These are analogous to the induction factors of conventional AD theory. Based on these we define a ‘global’, ‘foil’ and an ‘exit’ reduced frequency ($k_g$, $k_f$ and $k_e$) given by

(2.1)\begin{equation} k_g = \frac{\omega b}{U_\infty} = \alpha_2 k_f = \alpha_4 k_e, \end{equation}

where $\omega$ is the foil oscillation frequency in rad s$^{-1}$ and $b$ is the half-chord. Following AD theory convention, we assume the flow to be inviscid and the mean pressure to be fully recovered ($\bar {p}=p_\infty$) at the exit boundaries.

Figure 1. Control volumes used; (a) CV1, with mass-permeable side boundaries far from aerofoil, (b) CV2, with side boundaries far from aerofoil following the mean streamlines, (c) CV3, with side boundaries encompassing the AD following the mean streamlines.

2.2. Small-amplitude linear aerofoil theory: the Garrick function

The full form of the Garrick function will not be presented here; readers are referred to the original paper. For the purpose of this paper we present only an expression for the wake circulation distribution. For an aerofoil oscillating in a combination of pitch and heave in a fluid of density $\rho$, the wake circulation at downstream location $x$ and time $t$ is given by

(2.2)\begin{equation} \gamma(x,t) = A_0(t) \cos\left(\frac{k_ex}{b}\right) + B_0(t) \sin\left(\frac{k_ex}{b}\right), \end{equation}

where the exit reduced frequency $k_e$ is used by assuming a location far downstream of the aerofoil. The time-dependent variables $A_0(t)$ and $B_0(t)$ are given by

(2.3)$$\begin{gather} A_0 = 4[\zeta_1\sin(\omega t) - \zeta_2\cos(\omega t)], \end{gather}$$

(2.4)$$\begin{gather}B_0 = 4[\zeta_1\cos(\omega t) + \zeta_2\sin(\omega t)], \end{gather}$$

where $\zeta _1$ and $\zeta _2$ are functions of the aerofoil kinematics, and can be found in the original paper (also provided in a supplementary data sheet; see the Data availability statement below). They are also functions of the foil frequency $k_f$ and the local velocity $\alpha _2U_\infty$. In Garrick's original paper $k_f$ and $k_e$ were taken as equal to $k_g$ (that is, $\alpha _2=\alpha _4=1$). In the following section we will introduce the flow acceleration parameters through a cycle-averaged unsteady AD theory framework, which is subsequently coupled to the Garrick function to obtain the values of $k_f$, $k_e$, $\alpha _2$ and $\alpha _4$.

2.3. Cycle-averaged unsteady actuator disc theory

The principles of steady-flow AD theory are well known and can be found in textbooks such as Hansen (Reference Hansen2015). In this paper we will conduct an unsteady CV analysis assuming potential flow, and applying cycle averaging to predict the effect of unsteady flow on the mean thrust and propulsive efficiency. The unsteadiness is assumed to be generated entirely by aerofoil and wake vorticity, meaning that the unsteady components of bulk flow acceleration terms $\alpha _2$ and $\alpha _4$ are assumed negligible. Based on the results of Yu et al. (Reference Yu, Tavernier, Ferreira, van Kuik and Schepers2019), who used an unsteady AD model to estimate bulk flow oscillations in the wake of a wind turbine, this assumption is most likely to hold at high reduced frequencies. We retain the assumptions of linear aerofoil theory that the wake is planar and moving at the local mean velocity.

Similarly to Young et al. (Reference Young, Tian, Liu, Lai, Nadim and Lucey2020) but omitting the viscous terms, we begin with the integral equations for momentum and energy balance on a fluid CV with a control surface (CS)

(2.5)$$\begin{gather} f_i = \int_{CS} pn_i \,{\rm d} A + \int_{CS} \rho u_i u_j n_j \,{\rm d} A + \frac{\partial}{\partial t}\int_{CV} \rho u_i \,{\rm d} V \end{gather}$$

(2.6)$$\begin{gather}W = \int_{CS} \left(p+\rho\frac{u_i u_i}{2} \right)u_j n_j \,{\rm d} A + \frac{\partial}{\partial t}\int_{CV} \rho\frac{u_i u_i}{2} {\rm d} V , \end{gather}$$

with vector quantities given in tensor form. Here, $f_i$ is the force acting on the fluid and $W$ is the power input to the fluid by the aerofoil. Because we assume potential flow conditions, we can use the unsteady Bernoulli equation (neglecting gravity) given by

(2.7)\begin{equation} \frac{p}{\rho} + \frac{u_iu_i}{2} + \frac{\partial \varPhi}{\partial t} = \chi(t)= \left[ \frac{p}{\rho} + \frac{u_iu_i}{2} + \frac{\partial \varPhi}{\partial t} \right]_{ref}. \end{equation}

We introduce the time-dependent parameter $\chi (t)$ which denotes the reference value for the Bernoulli equation taken from a single point in the flow field, usually far from the aerofoil. The term $\partial \varPhi /\partial t$ is the time derivative of the velocity potential, and $\rho \partial \varPhi /\partial t$ is the added mass pressure of potential flows (see e.g. Katz & Plotkin (Reference Katz and Plotkin2001), chapter 13.7). Substituting for the pressure $p$ in (2.5) and (2.6)

(2.8)$$\begin{gather} f_i = \int_{CS} \rho\left( \chi(t)-\frac{u_ju_j}{2} - \frac{\partial\varPhi}{\partial t} \right)n_i \,{\rm d} A + \int_{CS} \rho u_i u_j n_j \,{\rm d} A + \frac{\partial}{\partial t}\int_{CV} \rho u_i \,{\rm d} V \end{gather}$$

(2.9)$$\begin{gather}W = \int_{CS} \rho\left( \chi(t)-\frac{u_iu_i}{2} - \frac{\partial\varPhi}{\partial t} +\frac{u_i u_i}{2} \right)u_j n_j \,{\rm d} A + \frac{\partial}{\partial t}\int_{CV} \rho\frac{u_i u_i}{2} {\rm d} V. \end{gather}$$

It is immediately clear that the velocity terms in the first integral of (2.9) cancel. We now divide the variables into mean and fluctuating components, such as $u=U+u'$ where the capital letter implies the mean value and a dash implies the fluctuating component. Alternatively, time averages are also denoted with overbars. Cycle averaging the momentum and energy balance equations, we simplify to

(2.10)$$\begin{gather} F_i = \rho \int_{CS} \left[\left(\overline{\chi} - \frac{U_jU_j + \overline{u'_ju'_j}}{2} \right)n_i + (U_i U_j +\overline{u'_i u'_j})n_j \right] {\rm d} A \end{gather}$$

(2.11)$$\begin{gather}\overline{W} = \rho \int_{CS} \left( \overline{\chi}U_j +\overline{\chi' u'_j} - \overline{\frac{\partial \varPhi}{\partial t}u'_j} \right) n_j \,{\rm d} A. \end{gather}$$

We can now solve (2.10) and (2.11) for CV1 and CV2, respectively, noting that only the exit faces of CV1 and CV2 will be affected by the unsteady terms.

2.4. Fluctuating velocity terms

In order to derive analytical expressions for the fluctuating velocity terms we make the following assumptions: the wake is planar, it travels at the local free-stream velocity, the vortex circulation is given by the Garrick function (2.2) and the exit face is far enough downstream of the aerofoil so that the wake can be approximated as extending to positive and negative infinity along the horizontal axis. Based on these assumptions, the fluctuating components of the velocity at a point $(x,y)$ on the exit face, induced by the wake vortex circulation along the horizontal axis (with vortices located at $x'$), are given by the Biot–Savart law (see Katz & Plotkin (Reference Katz and Plotkin2001), chapter 2) as

(2.12)$$\begin{gather} u'_x ={-}\int_{-\infty}^{\infty} \frac{A_0 \cos\left(\dfrac{k_ex'}{b}\right) + B_0 \sin\left(\dfrac{k_ex'}{b}\right)}{2{\rm \pi}} \frac{y}{[(x-x')^2 + y^2]} {{\rm d}\kern0.7pt x}' \end{gather}$$

(2.13)$$\begin{gather}u'_y = \int_{-\infty}^{\infty} \frac{A_0\cos\left(\dfrac{k_ex'}{b}\right) + B_0 \sin\left(\dfrac{k_ex'}{b}\right)}{2{\rm \pi}} \frac{x-x'}{[(x-x')^2 + y^2]} {{\rm d}\kern0.7pt x}'. \end{gather}$$

The evaluation of these integrals for $k_e > 0$ and $b > 0$ gives

(2.14)$$\begin{gather} u'_x ={-}\frac{\exp\left({-\dfrac{k_e}{b}|y|}\right)}{2}\frac{|y|}{y}\left[A_0\cos\left(\frac{k_ex}{b}\right) +B_0\sin\left(\frac{k_ex}{b}\right)\right] \end{gather}$$

(2.15)$$\begin{gather}u'_y = \frac{\exp\left({-\dfrac{k_e}{b}|y|}\right)}{2}\left[ A_0\sin\left(\frac{k_ex}{b}\right) - B_0\cos\left(\frac{k_ex}{b}\right)\right]. \end{gather}$$

From (2.10) we know that for the $x$-component of momentum at the exit face we will require expressions for $\overline {u'_xu'_x}$ and $\overline {u'_yu'_y}$ only. To evaluate these, for simplicity we say $x=0$ at the exit face, such that

(2.16)$$\begin{gather} \overline{u'^2_x} = \frac{\exp\left({-2\dfrac{k_e}{b}|y|}\right)}{4}\overline{A_0^2} \end{gather}$$

(2.17)$$\begin{gather}\overline{u'^2_y} = \frac{\exp\left({-2\dfrac{k_e}{b}|y|}\right)}{4}\overline{B_0^2}. \end{gather}$$

Based on the definitions of $A_0$ and $B_0$ in (2.3) and (2.4), we see that $\overline {A_0^2} = \overline {B_0^2}$, and thus that $\overline {u'_xu'_x} = \overline {u'_yu'_y}$ over the CV exit face. Note that it can be demonstrated that $[u_x']_{y=0_\pm } = \mp \gamma /2$ (Katz & Plotkin Reference Katz and Plotkin2001, chapter 3), which suggests the equality $\overline {u'_xu'_x} = \overline {u'_yu'_y}$ holds also at $y=0$. Since $n_y=0$ for the exit face, the $\overline {u'_iu'_j}$ term in (2.10) is $\overline {u'_xu'_x}$ for $i=x$. Thus the cycle-averaged fluctuating terms in (2.10) cancel on the exit face.

2.5. The Bernoulli equation reference term

We introduced $\chi (t)$ to represent the reference parameter used in the unsteady Bernoulli equation (2.7). This term can be taken as the total far-field pressure $\chi = p_\infty /\rho + U_\infty ^2/2$ at all CV boundaries, except for at the exit boundary of CV3 (figure 1c). Here, the energy discontinuity created by the AD means that the reference point must be taken downstream of the AD and between the two mean-flow streamlines that define the CV boundary. However, by choosing the reference point on the internal face of the CV boundary, indicated by ‘ref-i’ in figure 1(c), we can simplify further. We define an additional reference point, marked ‘ref-e’ in figure 1(c), at the same position but on the external face of the CV boundary. The assumption of fully developed flow at the exit face of the CV suggests that the pressure at the two reference points must be equal. Applying (2.7) to obtain the pressure at ref-e, noting that $\chi =p_\infty /\rho +U_\infty ^2/2$ outside the streamtube, we get

(2.18)\begin{equation} \frac{p_{ref\text{-}i}}{\rho} = \frac{p_{ref\text{-}e}}{\rho} = \frac{p_\infty}{\rho} + \frac{U_\infty^2}{2} - \left[ \frac{(U_\infty+u'_x)^2+u'^2_y}{2} + \frac{\partial \varPhi}{\partial t} \right]_{ref\text{-}e}. \end{equation}

Using (2.18) as the pressure term in (2.7), we obtain $\chi$ downstream of the AD in CV3 as

(2.19)\begin{align} \chi(t) &= \frac{p_\infty}{\rho} + \frac{U_\infty^2}{2} - \left[\frac{(U_\infty+u'_x)^2+u'^2_y}{2} + \frac{\partial \varPhi}{\partial t} \right]_{ref\text{-}e} \nonumber\\ &\quad +\left[ \frac{(\alpha_4U_\infty+u'_x)^2+u'^2_y}{2} + \frac{\partial \varPhi}{\partial t} \right]_{ref\text{-}i}. \end{align}

Because we assume that the vortex-induced flow is the only source of unsteadiness, and this has no discontinuity across the CV boundary, the expressions in brackets in (2.19) are equal except for the $\alpha _4$ terms. Thus (2.19) can be simplified to

(2.20)\begin{equation} \chi(t) = \frac{p_\infty}{\rho} +\frac{\alpha_4^2U_\infty^2}{2} + [u_x']_{ref}U_\infty(\alpha_4-1), \end{equation}

where $[u_x']_{ref}$ indicates the fluctuating axial velocity at the reference location.

2.6. Added mass energy term

Since the time derivative of the potential field is needed, only the unsteady component of the flow potential will be considered, which is assumed fully determined by the wake vorticity. Only the added mass on the exit face is needed. The potential of a free vortex is given by $\varPhi =\gamma \theta /2{\rm \pi}$ (Katz & Plotkin Reference Katz and Plotkin2001, chapter 3), where $\theta$ is the angle between the point of interest and the horizontal axis intersecting the vortex core. The potential field induced by the wake circulation (2.2) distributed along $x'$, at a point $(x,y)$ on the exit face, is then given by

(2.21)\begin{equation} \varPhi = \frac{1}{2{\rm \pi}}\int_{-\infty}^{\infty} \left[A_0 \cos\left(\frac{k_ex'}{b}\right) + B_0 \sin\left(\frac{k_ex'}{b}\right)\right] \tan^{{-}1} \left(\frac{y}{x-x'} \right) {{\rm d}\kern0.7pt x}'. \end{equation}

Evaluating the integral at $x=0$, the potential is

(2.22) \begin{equation} \varPhi(x=0,y) ={-} b\frac{|y|}{y}\frac{\left[1-\exp\left({-\dfrac{k_e}{b}|y|}\right)\right]}{2k_e} B_0. \end{equation}

Taking the time derivative and noting that $\partial B_0/\partial t = -\omega A_0$, we obtain

(2.23) \begin{equation} \frac{\partial\varPhi}{\partial t}(x=0,y) = \frac{b\omega }{2k_e} \frac{|y|}{y}A_0 \left[1-\exp\left({-\dfrac{k_e}{b}|y|}\right)\right].\end{equation}

We then evaluate $\overline {\partial \varPhi /\partial (t u_x')}$ for the energy balance in (2.11) using the expression for $u'_x$ from (2.14)

(2.24) \begin{equation} \overline{\frac{\partial\varPhi}{\partial t} u_x'}(x=0,y) ={-} \frac{b\omega }{4k_e} \overline{A_0^2} \left[1-\exp\left({-\frac{k_e}{b}|y|}\right)\right]\exp\left({-\frac{k_e}{b}|y|}\right). \end{equation}

2.7. Momentum balance

To evaluate the cycle-average momentum balance (2.10) we use CV1 (figure 1a). The CV is rectangular and all boundaries (shown by dashed black lines) are mass permeable. The upper and lower boundaries are assumed to be far from the aerofoil and wake, such that unsteady-flow effects are negligible everywhere except over the central part of the exit face. We demonstrated in § 2.4 that the fluctuating components of velocity in (2.10) cancel on the exit face. This removes all unsteady terms from (2.10), reducing to the momentum balance for steady AD theory. The evaluation procedure is well known (see e.g. Hansen Reference Hansen2015) and we can obtain

(2.25)\begin{equation} \frac{F_x}{\tfrac{1}{2}\rho U_\infty^2\mathcal{A}} = C_{Tg} = 2\alpha_2 (\alpha_4 - 1). \end{equation}

Here, we define the global thrust coefficient $C_{Tg}$ using the AD area $\mathcal {A}$, given by the frontal area swept by the oscillating aerofoil (see figure 1c).

2.8. Energy balance

We use CV2 (figure 1b) to evaluate (2.11). The upper and lower boundaries follow the mean-flow streamlines, such that the exit area can be found (through mean-flow mass conservation) to be $\mathcal {A}_{2} - \mathcal {A}({\alpha _2}/{\alpha _4})(\alpha _4 - 1)$, where $\mathcal {A}_2$ is the inlet area. Again, the upper and lower boundaries are assumed far enough from the aerofoil so that unsteady effects are negligible everywhere except at the central part of the exit face. There are two unsteady terms in (2.11), one related to $\chi '$ and the other to the added mass. At the exit face, as shown in § 2.5, $\chi '=[u_x']_{ref}U_\infty (\alpha _4-1)$, which is constant in $y$. From (2.14) we see that $u'_x$ is anti-symmetric in $y$. Thus the integral of $\overline {\chi ' u'_x}$ over the exit face is zero. This leaves the added mass term as the only unsteady-flow contribution to the energy balance. Evaluating (2.11) for each CV boundary, recalling the expressions for $\chi$ inside and outside the wake streamtube from § 2.5, we obtain

(2.26)\begin{align} \frac{\overline{W}}{\rho} &={-}\mathcal{A}_{2} U_\infty \left[\frac{p_\infty}{\rho} + \frac{U_\infty^2}{2} \right] + \left(\mathcal{A}_{2} - \mathcal{A}\frac{\alpha_2}{\alpha_4}(\alpha_4 - 1) - \mathcal{A}\frac{\alpha_2}{\alpha_4} \right) U_\infty\left[ \frac{p_\infty}{\rho} + \frac{U_\infty^2}{2} \right] \nonumber\\ &\quad +\mathcal{A}\frac{\alpha_2}{\alpha_4}\alpha_4U_\infty\left[ \frac{p_\infty}{\rho} + \frac{\alpha_4^2U_\infty^2}{2}\right] - \int_{-\infty}^{\infty} \overline{\frac{\partial \varPhi}{\partial t}u'_x} {{\rm d}y}. \end{align}

We can integrate the added mass term over ${\pm }\infty$ without loss of generality since there are no unsteady wake effects at the upper and lower CV boundaries. Cancelling terms and normalising to obtain the global power coefficient $C_{Pg}$ gives

(2.27)\begin{equation} \frac{\overline{W}}{\tfrac{1}{2}\rho U_\infty^3\mathcal{A}} = C_{Pg} = \alpha_2 ( \alpha_4^2 - 1) - \frac{1}{\tfrac{1}{2}\rho U_\infty^3\mathcal{A}} \int_{-\infty}^{\infty} \rho\overline{\frac{\partial \varPhi}{\partial t}u'_x} {{\rm d}y}. \end{equation}

Evaluating the added mass integral from (2.24), noting that it is symmetric in $y$

(2.28) \begin{align} 2\int_{0}^{\infty} \rho\overline{\frac{\partial \varPhi}{\partial t}u'_x} {{\rm d}y} ={-}\rho\frac{b\omega}{2 k_e} \overline{A_0^2} \int_{0}^{\infty} \left[1 - \exp\left({-\frac{k_e}{b}y}\right)\right] \exp\left({-\frac{k_e}{b}y}\right) {{\rm d}y} ={-} \rho\frac{\omega b^2}{4k_e^2}\overline{A_0^2}. \end{align}

In steady AD theory the energy input $\overline{W}$ to the CV is the energy required for generating the thrust of an ideal disc propulsor, that is $\overline{W} = \alpha _2U_\infty F_x$. However, for the present non-ideal case the total oscillation energy of the aerofoil $\overline{W}_f$ must be considered, which is equal to the thrust energy plus the energy required to generate the wake (Garrick Reference Garrick1937), i.e. $\overline{W} = \overline{W}_f \geq \alpha _2U_\infty F_x$. Here, $\overline {W}_f$ is obtained from the chordwise integration of lift force times the vertical aerofoil velocity, and is evaluated analytically by Garrick. Note that for $\alpha _2=\alpha _4=1$, (2.28) is equivalent to the Garrick wake energy. To account for $\overline{W}_f$ we define the ‘local’ efficiency $\eta _l$ in relation to the ‘global’ efficiency $\eta _g$

(2.29)\begin{equation} \eta_l = \frac{\alpha_2 U_\infty F_x}{\overline{W}_f} = \alpha_2\eta_g. \end{equation}

Incorporating the expressions from (2.28)–(2.29) into (2.27), and noting that $k_e = k_g/\alpha _4$, we obtain

(2.30)\begin{equation} \frac{F_x}{\tfrac{1}{2}\rho U_\infty^2\mathcal{A}} = C_{Tg} = \eta_l [\alpha_4^2\eta_{am} - 1], \end{equation}

where we have introduced the parameter $\eta _{am}$ to represent the added mass term

(2.31)\begin{equation} \eta_{am} = 1 + \frac{1}{2\alpha_2k_g} \frac{\overline{A_0^2} b }{U_\infty^2\mathcal{A}}. \end{equation}

Equations (2.25) and (2.30) thus define the cycle-averaged inviscid AD equations for this system. As mentioned, (2.25) is equivalent to the steady formulation, but (2.30) differs through the $\eta _l$ and $\eta _{am}$ terms. The standard steady AD theory result for an ‘ideal’ propulsor (that is when $\overline{W}=\alpha _2U_\infty F_x$) is recovered by setting $\eta _l=\eta _{am}=1$.

3.1. Numerical procedure

Equations (2.25), (2.29), (2.30) and (2.31), and the Garrick function expressions for thrust $F_x$ and power $\overline{W}_f$, are solved iteratively using the MATLAB function fsolve until convergence of all variables, adjusting the Garrick function for the local flow acceleration $\alpha _2$ and reduced frequency $k_f$. The cycle-averaged AD model is only evaluated for $C_{Tg}>0$. To validate the new model we evaluate the case of an aerofoil flapping in combined heave and pitch, and compare the results with LES.

We use an immersed-boundary implicit-LES solver called the boundary data immersion method (BDIM) to solve the three-dimensional incompressible Navier–Stokes equations. The solver has been validated in several previous studies of flapping foils with chord-based Reynolds numbers (Re) up to $Re=50\,000$ (Maertens & Weymouth Reference Maertens and Weymouth2015; Zurman-Nasution, Ganapathisubramani & Weymouth Reference Zurman-Nasution, Ganapathisubramani and Weymouth2020). The foil is a NACA0016, and the kinematics consist of a combination of heave $H(t)$ and pitch $\theta (t)$ motions, with functional form

(3.1)$$\begin{gather} H(t) = h_0 \sin(\omega t), \end{gather}$$

(3.2)$$\begin{gather}\theta(t) = \alpha_{0} \sin(\omega t + \psi), \end{gather}$$

(3.3)$$\begin{gather}\alpha_{0} = \sin^{{-}1}\left(\frac{0.7h_0}{1.5b}\right), \end{gather}$$

where $h_0$ is heave amplitude and $\alpha _{0}$ is pitch amplitude in radians, and the heave–pitch phase difference is $\psi =90^\circ$. The pitch axis is located at the quarter-chord from the leading edge. For a given heave amplitude in the range of $0.4b< h_0< b$, we vary the Strouhal number $St=\omega h_0/({\rm \pi} U_\infty$) in the range $0.15< St<0.80$ to vary the frequency $\omega$. The AD area is found from the heave amplitude as $\mathcal {A} = 2h_0$.

The simulations were performed at $Re=10\,000$, which was deemed sufficiently high since the thrust coefficient of a flapping foil is almost invariant at $Re>10\,000$ (Senturk & Smits Reference Senturk and Smits2019). The domain extends horizontally from $-6b$ to $24b$ and vertically from $-8b$ to $8b$. To ensure domain height independence, two representative cases were re-run at double domain height; the resulting propulsive efficiencies changed by less than 1 %. The foil has a spanwise width of $2b$ and periodic boundary conditions applied to both sides. The foil and its near wake are simulated within a sub-domain using a uniform Cartesian grid with a resolution of $\Delta Y = 2b/128$ to reach (on average) $y^+=y_n u_{\tau }/\nu \approx 5$, $\Delta X=2\Delta Y$ and $\Delta Z=2\Delta Y$. Here, $y_n$ is the wall-normal distance, $u_{\tau }$ is friction velocity, $\nu$ is kinematic viscosity and $X, Y, Z$ are global coordinates for horizontal, vertical and spanwise directions, respectively.