Coupled unsteady actuator disc and linear theory of an oscillating foil propulsor

Abstract

Linear unsteady aerofoil theory, while successfully used for the prediction of unsteady aerofoil lift for many decades, has yet to be proven adequate for predicting the propulsive performance of oscillating aerofoils. In this paper we test the hypothesis that the central shortcoming of linear small-amplitude models, such as the Garrick function, is the failure to account for the flow acceleration caused by aerofoil thrust. A new analytical model is developed by coupling the Garrick function to a cycle-averaged actuator disc model, in a manner analogous to the blade-element momentum theory for wind turbines and propellers. This amounts to assuming the Garrick function to be locally valid and, in combination with a global control volume analysis, enables the prediction of flow acceleration at the aerofoil. The new model is demonstrated to substantially improve the agreement with large-eddy simulations of an aerofoil in combined heave and pitch motion.

1. Introduction

The Theodorsen function has been successfully used over the last century for the prediction of unsteady harmonic aerofoil lift in applications requiring analytical solutions, low computational cost or fast computations. An extension to the Theodorsen function was derived by Garrick (1937) to also include the propulsive thrust of a foil oscillating in heave and/or pitch. The function is derived based on the same underlying assumptions as those of Theodorsen: potential flow, the aerofoil represented by a flat plate, small-amplitude motion and the wake assumed to be co-planar with the aerofoil and moving with the free-stream velocity. However, the Garrick function has been demonstrated to severely over-predict the propulsive efficiencies of oscillating foils relative to experiments and simulations (e.g. Fernandez-Feria 2016; Faure et al. 2022), leading to the supposition that the inviscid small-amplitude assumptions are inappropriate for the propulsive foil problem.

In this paper we hypothesise that the shortcomings of small-amplitude linear aerofoil theory are largely explained by the neglecting of the axial flow acceleration that results from the aerofoil thrust – a well-known phenomenon that can be modelled using actuator disc (AD) theory. Coupling the Garrick function to an AD model is analogous to blade-element momentum (BEM) theory used for e.g. wind turbines. This amounts to assuming the Garrick function to be locally valid, and coupled to a global control volume analysis through an AD (representing the frontal area swept by the oscillating foil) via the axial induction factor. We develop a cycle-averaged unsteady AD model and demonstrate, by comparing with large-eddy simulations (LES), that both steady and cycle-averaged AD coupling give substantial improvements to the Garrick function. The results suggest that linear small-amplitude models are able to make reliable predictions for foil propulsion trends, as long as the local flow acceleration is taken into account. The outcome of this study is a simple and fully analytical model for oscillating foil propulsion.

2. Theoretical model

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 (2015). 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. (2019), 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. (2020) 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 (2001), 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 (2001), 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 2001, 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 2001, 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 2015) 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 1937), 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. Numerical procedure and results

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 2015; Zurman-Nasution, Ganapathisubramani & Weymouth 2020). 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 2019). 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.

3.2. Results

Figure 2(a–c) compares results for the global propulsive efficiency ($\eta _g$, (2.29)) and the foil power and thrust coefficients ($C_{Pf}= \overline{W}_f/\rho U_\infty ^3b$ and $C_{Tf}= F_x/\rho U_\infty ^2b$) normalised by $k_g^2 a^2=St^2{\rm \pi} ^2/4$, where $a=h_0/2b$ is the non-dimensional heave amplitude. Figure 2(d–g) shows the acceleration parameters ($\alpha _2$ and $\alpha _4$), local foil efficiency ($\eta _l$) and added mass parameter ($\eta _{am}$), all for the same set of cases and plotted against the global reduced frequency $k_g$. The original Garrick predictions (black lines) have significant errors in both global efficiency (figure 2a) and power (figure 2b) relative to the LES (circles), while the AD-coupled models substantially improve the agreement of both. The AD coupling also improves the thrust prediction (figure 2c), although it is marginal compared with its effect on the power.

Figure 2. (a) Global propulsive efficiency. (b) Power coefficient with visualisations of LES results (aerofoil in green and isosurfaces of $\lambda _2$-criterion coloured by spanwise vorticity). (c) Thrust coefficient. (d) Acceleration parameter at the foil. (e) Acceleration parameter at the exit face. (f) Local foil efficiency. (g) Added mass parameter.

The steady AD ($\eta _{am}=\eta _l=1$, dashed lines) and the cycle-averaged AD (solid coloured lines) give similar trends in foil performance, the latter agreeing better with the LES, especially at high frequencies. The similarity between steady and cycle-averaged AD predictions is due to the non-ideal energy input by the foil ($\eta _l<1$) being largely balanced by the wake energy exiting the CV ($\eta _{am} > 1$). Figure 2(d,e) indicates the significance of local flow accelerations; the velocity at the foil is up to approximately $4$ times $U_\infty$, and velocity at the CV exit even higher. The steady AD under-predicts $\alpha _2$ and over-predicts $\alpha _4$ compared with the cycle-averaged AD.

Figure 2(f) shows that $\eta _l$ approaches $\approx$ 0.58 to 0.6 with increasing $k_g$, suggesting that the local efficiency in the AD model is close but not equal to the unmodified Garrick efficiency (black line in figure 2a). Figure 2(g) shows that $\eta _{am}$ increases with $k_g$, and the added mass energy becomes larger than the mean-flow energy (that is, $\eta _{am}>2$) at $k_g>4$ to 4.5. There are small increases in $\eta _l$ and $\eta _{am}$ with increasing amplitude $h_0$.

The remaining discrepancies in foil performance between the AD models and LES are likely due to the factors not accounted for in the former, such as viscous effects and bulk flow oscillations. The LES results also capture strong vortex instabilities in the near wake when the Strouhal number is above the optimum range of $0.2 < St < 0.5$ (Zurman-Nasution et al. 2020) as shown in the flow field insets in figure 2(b), which may partly explain the discrepancies at higher frequencies. Furthermore, in all LES cases the aerofoil motion was found to result in the generation of leading-edge vortices (LEVs), which increased in strength with the oscillation amplitude. In a recent comparative study of low-order models for flapping foil propulsion by Faure et al. (2022), models implementing dynamic stall or LEV corrections are shown to give improved agreement with experiments and high-order simulations over inviscid models. Considering their results, it is likely that the remaining discrepancies in figures 2(a–c) are largely due to the absence of stall effects in the Garrick-AD models. The form drag induced by trailing-edge vortex rollup is also not included, which may affect the prediction of both lift and thrust (Ayancik et al. 2019).

Despite these discrepancies, the coupled Garrick-AD theory provides a fully analytical solution for inviscid foil propulsors that correctly represents efficiency trends, and is a substantial improvement on the original Garrick theory. The trends of the coupled models are similar to those of the numerical inviscid panel method evaluated by Faure et al. (2022): figure 5 in their paper shows the panel method giving an error in $\eta _g$ of approximately 50 % at $k_g=3$ (the Garrick theory error is about 200 %) relative to high-order simulations of heaving foils. The present method has similar accuracy but at a fraction of the computational time, without requiring a numerical panel solver. More generally, the AD coupling method also opens the possibility for further analytical modelling of unsteady foils, and for using small-amplitude inviscid models to evaluate finite-amplitude viscous problems to first-order accuracy.

4. Conclusions

We have developed a method for coupling a linear unsteady aerofoil theory (Garrick 1937) and an inviscid AD theory, analogous to the BEM theory for wind turbines and propellers, to improve analytical prediction of the propulsive performance of oscillating foils. By cycle averaging the integral forms of the inviscid momentum and energy conservation equations for three different CVs, we have derived concise analytical expressions linking the mean foil thrust and power to the local flow acceleration at the foil. The cycle-averaged AD model deviates from steady AD theory through only two additional parameters, $\eta _l$ and $\eta _{am}$. The former accounts for the ‘non-ideal’ energy input by the foil, and the latter for the added mass energy in the wake, both of which are obtained from the Garrick theory. Both the steady and cycle-averaged AD models coupled to the Garrick theory were shown to substantially improve agreements with LES, although some discrepancies remained, especially in the thrust prediction. It is likely that these discrepancies are largely due to the effects of LEV formation, which is not accounted for in the present model. The results demonstrate the applicability of small-amplitude inviscid unsteady aerofoil theory to finite-amplitude foil propulsion problems, as long as the local flow acceleration at the foil is taken into account.