4.1. Lift dynamics

First, we focus on the lift dynamics; from which, we will construct the dynamics of drag/thrust and the separation point. Second, we follow a feedback linearization approach (Sastry Reference Sastry1999). In this approach, given a nonlinear dynamical system (2.1), the following question is addressed: does there exist (nonlinear) transformations $\boldsymbol {z}=\boldsymbol {T}(\boldsymbol {x})$ and $\boldsymbol {v}=\boldsymbol {\beta } (\boldsymbol {x},\boldsymbol {u})$ such that the dynamics in the transformed domain is linear: $\dot {\boldsymbol {z}}=\boldsymbol {A}\boldsymbol {z}+\boldsymbol {B}\boldsymbol {v}$ (Nijmeijer & Van der Schaft Reference Nijmeijer and Van der Schaft1990)? That is, while the lift (and its dynamics) are essentially nonlinear with respect to the angle of attack (in the stall regime), does there exist a nonlinear map $v=\beta (\boldsymbol {x},\alpha )$ such that the lift dynamics is linear in $v$? Interestingly, although the seminal results of Von Kármán & Sears (Reference Von Kármán and Sears1938) were developed approximately 40 years before the concept of feedback linearization (Brockett Reference Brockett1978; Jakubczyk & Respondek Reference Jakubczyk and Respondek1980), they can shed some light onto this question. Recall that, in the earlier efforts of Wagner (Reference Wagner1925) and Theodorsen (Reference Theodorsen1935), the authors wrote the main governing equation (integral equation of wake circulation) in terms of the normal speed of the airfoil

(4.1)\begin{equation} \frac{1}{2{\rm \pi} b} \int_b^\infty \gamma_w(z,t)\sqrt{\frac{z+b}{z-b}}\textrm{d}z=U\sin\alpha(t), \end{equation}

where $b$ is the half-chord length, $U$ is the forward speed of the airfoil and $\gamma _w$ is the strength of the wake vortex sheet per unit span. Equation (4.1) is an infinite-dimensional linear dynamical system – represented by an integral equation in the wake circulation $\gamma _w$ – whose input is the airfoil normal speed $U\sin \alpha (t)$. This formulation represents an impasse against generalization to unconventional lift mechanisms; it will always result in the classical lift coefficient $2{\rm \pi} \sin \alpha$ with some transient (dynamics). In contrast, Von Kármán & Sears (Reference Von Kármán and Sears1938) arrived at a fundamentally different representation of the same dynamical system

(4.2)\begin{equation} \int_b^\infty \gamma_w(z,t)\sqrt{\frac{z+b}{z-b}}\textrm{d}z=\varGamma_0(t), \end{equation}

where $\varGamma _0$ is the quasi-steady circulation

(4.3)\begin{equation} \varGamma_0(t)=U b C_{L,s}(\alpha(t)) + b^2 k_{\dot{\alpha}}\dot{\alpha}(t), \end{equation}

and $k_{\dot {\alpha }}$ is a coefficient for the rotational circulation, which depends on the hinge location (Dickinson, Lehmann & Sane Reference Dickinson, Lehmann and Sane1999). It is well known that the lift dynamics resulting from the two dynamical systems (4.1)–(4.2) are identical. However, the first formulation is nonlinear in its input $\alpha$ while the latter is linear in its input $\varGamma _0$. That is, the nonlinear transformation (4.3) provides the sought feedback linearization; and $\varGamma _0$ represents the new input $v$. More importantly, the first formulation captures only the conventional lift mechanism (e.g. $2{\rm \pi} \sin \alpha$) while the latter allows for unconventional (possibly nonlinear) lift mechanisms. To illustrate this point, recall that in (4.3) for the quasi-steady circulation, $C_{L,s}$ is the steady lift coefficient, which is a function of the angle of attack. This functional dependence $C_{L,s}(\alpha )$ is arbitrary, allowing treatment of unconventional lift mechanisms (e.g. a leading-edge vortex).

The main assumption in the present modelling effort is that the lift dynamics is linear in the quasi-steady circulation. That is, the circulatory lift can be written as

(4.4)\begin{equation} \left.\begin{gathered} \dot{\boldsymbol{x}}_c(t) = [\boldsymbol{A}]_{n\times n} \boldsymbol{x}_c(t) + [\boldsymbol{B}]_{n\times 1}\varGamma_0(t) \\ L_C(t) = \rho U \left\{ [\boldsymbol{C}]_{1\times n}\boldsymbol{x}_c(t) + [D]_{1\times 1}\varGamma_0(t)\right\}, \end{gathered}\right\} \end{equation}

where $\boldsymbol {x}_c\in \mathbb {R}^n$ represents the internal aerodynamic states used to realize such a linear dynamical system. The above assumption may not be a severe one given that no limitation is set on the nonlinear function $C_{L,s}(\alpha )$, the order $n$ of the system (4.4) or its nature (i.e. stable/unstable); its eigenvalues may even change with the operating condition (i.e. with the mean angle of attack). In potential flow, such a system can be given by any suitable finite-dimensional approximation of Wagner or Theodorsen response functions – and $D={1}/{2}$ is the high-frequency gain (Taha & Rezaei Reference Taha and Rezaei2020).

Since the focus here is to reconstruct the dynamics of the output, the nature (physical meaning) of the states $\boldsymbol {x}_c$ is not important. Note that one can always use a canonical similarity transformation $\boldsymbol {x}_{c,1}=\boldsymbol {Tx}_{c,2}$, where $\boldsymbol {T}$ is invertible, to transform between a given set of states $\boldsymbol {x}_{c,1}$ to another $\boldsymbol {x}_{c,2}$; in fact, we have infinitely many of these sets. While the internal aerodynamic states may not have an obvious physical meaning in many finite-dimensional approximations of the unsteady lift dynamics, few efforts presented realizations of the dynamical system (4.4) with states of physical relevance; Peters (Reference Peters2008) presented his model in terms of the Fourier coefficients/modes of the wake-induced normal velocity on the airfoil (inflow). In fact, the internal states of any realization of (4.4) can always be thought of as weighted modes of the wake vorticity.

Given (i) the steady $C_{L,s}$-$\alpha$ curve corresponding to some unconventional lift mechanism, and (ii) the airfoil motion represented by $U$ and $\alpha (t)$, this approach provides the unsteady lift dynamics that captures such an unconventional lift mechanism. Figure 3 shows a schematic for this procedure. In fact, this approach was successfully adopted to model the unsteady nonlinear lift dynamics in insect flight (Taha et al. Reference Taha, Hajj and Beran2014). The model captures the leading-edge vortex (LEV) nonlinear contribution to lift in an unsteady fashion, relying on Wang's empirical formula for the steady lift due to a stabilized LEV: $C_{L,s}=A\sin 2\alpha$ (Wang, Birch & Dickinson Reference Wang, Birch and Dickinson2004), where $A$ is a constant. The unsteady model achieves such a challenging task with a very cheap computational cost and compact form (state space form).

Figure 3. A schematic diagram for the circulatory lift dynamics and its linear dependence on the quasi-steady circulation.

The non-circulatory lift $L_{NC}$ can be easily modelled linearly in the normal acceleration $a_{1/2}$ of the mid-chord point

(4.5)\begin{equation} L_{NC}(t) = m_v a_{1/2}(t) \cos\alpha(t), \end{equation}

where $m_v$ is the virtual (added) mass, which is given by $m_v={\rm \pi} \rho b^2$ in a potential-flow setting, and the airfoil normal acceleration is given by

(4.6)\begin{equation} a_{1/2}(t)=\frac{\textrm{d}}{\textrm{d}t}\left(U\sin\alpha(t)+\dot{h}(t)\cos\alpha(t)-ab\dot{\alpha}(t)\right), \end{equation}

where the hinge location is at a distance $ab$ behind the mid-chord point. The relation (4.5) would make the output (lift) have a direct dependence on the input (acceleration); i.e. $\boldsymbol {\varPsi }$ would have dependence on the input $\boldsymbol {u}$ as well as the states $\boldsymbol {x}$. Furthermore, from a different perspective, our recent efforts indicate that viscosity induces a phase lag between the normal acceleration and the non-circulatory lift force (Taha & Rezaei Reference Taha and Rezaei2020). Therefore, we introduce a first-order lag between the two

(4.7)\begin{equation} \left.\begin{gathered} \tau_v\dot{x}_v(t) +x_v(t) = a_{1/2}(t) \\ L_{NC}(t) = m_v x_v(t) \cos\alpha(t), \end{gathered}\right\} \end{equation}

where $\tau _v$ is the time constant of such a dynamics. That is, the normal acceleration $a_{1/2}$ does not impact the non-circulatory lift $L_{NC}$ instantaneously (i.e. through an algebraic relation: (4.5)). Instead, it affects $L_{NC}$ after passing through a first-order lag whose internal state is $x_v$.

One must stress that the above modelling approach is not meant to provide an accurate quantitative assessment of the unsteady nonlinear lift dynamics in the stall regime. Rather, this modelling approach is conceptual/qualitative. For example, conceptually, the non-circulatory lift is indeed proportional to the airfoil acceleration, and the proportionality constant ($m_v$) is not stipulated here; it can take any arbitrary value. In fact, its value is irrelevant for the intended symbolic analysis; it would rather serve as a bookkeeping parameter. That is, a resulting term proportional to $m_v$ would imply that it is due to added/virtual mass effects; another term proportional to $k_{\dot {\alpha }}$ would point to rotational contributions. Similarly, the steady $C_{L,s}$-$\alpha$ curve can assume any shape; a term proportional to ${\textrm {d}C_{L,s}}/{\textrm {d}\alpha }$ would point to the role of the lift curve slope, and so on. Moreover, the dynamics between $\varGamma _0$ and $L_C$ can take an arbitrary order $n$ or nature (arbitrary $\boldsymbol {A}$, $\boldsymbol {B}$, $\boldsymbol {C}$, $D$). From this conceptual perspective, the above model is credible for the ensuing symbolic analysis; the specific values of the model parameters ($n$, $\boldsymbol {A}$, $\boldsymbol {B}$, $\boldsymbol {C}$, $D$, $C_{L,s}(\alpha )$, $k_{\dot {\alpha }}$, $m_v$, $\tau _v$) are irrelevant.

Recall that one can always use a canonical similarity transformation to transform a given linear system into a special form, we assume without loss of generality (to simplify symbolic computations) that the linear model ($\boldsymbol {A}$, $\boldsymbol {B}$, $\boldsymbol {C}$, $D$) is in the observable canonical form (Ogata & Yang Reference Ogata and Yang1970)

(4.8)\begin{equation} \left.\begin{gathered} \boldsymbol{A}=\left[\begin{array}{cccc} 0 & \cdots & 0 & -a_0 \\ 1 & \cdot & 0 & -a_1 \\ \vdots & \ddots & \vdots & \vdots\\ 0 & \cdots & 1 & -a_{n-1} \end{array} \right], \quad \boldsymbol{B}=\left(\begin{array}{c} b_0-b_n a_0 \\ b_1-b_n a_1 \\ \vdots \\ b_{n-1}-b_n a_{n-1} \end{array} \right),\\ \boldsymbol{C}=\left[ 0 \quad \dots \quad 0 \quad 1\right], \quad D=b_n, \end{gathered}\right\} \end{equation}

where the coefficients $a$ and $b$ are the coefficients of the denominator and numerator of the corresponding transfer function

(4.9)\begin{equation} \frac{L_C}{\rho U \varGamma_0}(s)=\frac{b_n s^n + b_{n-1}s^{n-1} + \cdots + b_1 s+b_0}{s^n + a_{n-1}s^{n-1}+\cdots+a_1 s + a_0}. \end{equation}

From a dynamical-system perspective, the ratio ${b_0}/{a_0}$ is the steady-state (DC) gain of the system: the ratio of the steady-state output (lift) to the input $\rho U \varGamma _0$. Clearly, this is unity by construction. Therefore, in the ensuing symbolic computations, $b_0$ is taken equal to $a_0$. Also, $D=b_n$ is the high-frequency gain $k_{hf}$ of the system.

4.2. Drag/thrust dynamics

Continuing with the conceptual modelling approach, and ignoring the skin friction drag, the pressure drag is typically given by

(4.10)\begin{equation} \mathcal{D}(t)=L(t)\tan\alpha(t) - F_S(t), \end{equation}

where $L$ is the total lift force ($L=L_C+L_{NC}$) and $F_S$ is the leading-edge suction force. Indeed, without leading-edge suction, the resultant force $N$ is normal to the wing (Schlichting & Truckenbrodt Reference Schlichting and Truckenbrodt1979); the lift and drag would be given as $N\cos \alpha$ and $N\sin \alpha$, respectively; i.e. $\mathcal {D}=L\tan \alpha$.

The dynamics of the first drag term ($L\tan \alpha$) is automatically linked to the lift dynamics. Interestingly, the dynamics of the second term (the suction force $F_S$) can also be connected with the circulatory lift dynamics. Using potential-flow aerodynamics, Garrick (Reference Garrick1937) wrote the suction force as

(4.11)\begin{equation} F_S(t)=2{\rm \pi}\rho b \left[v_{3/4}(t) C(k)-\frac{b}{2}\dot{\alpha}(t)\right]^2, \end{equation}

where $v_{3/4}$ is the airfoil normal speed at the three-quarter chord and $C(k)$ is Theodorsen's lift frequency response function (Theodorsen Reference Theodorsen1935), which is given in terms of the reduced frequency $k$ as

(4.12)\begin{equation} C(k)=\frac{H_1^{(2)}(k)}{H_1^{(2)}(k)+iH_0^{(2)}(k)}, \end{equation}

where $H_n^{(m)}$ is the Hankel function of the $m\textrm {{th}}$ kind of order $n$. Note that the multiplication $v_{3/4}(t) C(k)$ in (4.11) is traditionally interpreted as the output of the transfer function whose frequency response is $C(k)$ due to the input signal $v_{3/4}$; i.e. it is the unsteady version of $v_{3/4}$.

To generalize Garrick's relation (4.11), we replace $v_{3/4}$ with the quasi-steady circulation $\varGamma _0/(2{\rm \pi} b)$, as done with the lift dynamics above. Moreover, its unsteady version $v_{3/4} C(k)$ will be replaced by the unsteady version of $\varGamma _0$; i.e. the output of the dynamical system (4.4): $\boldsymbol {Cx}_c+D\varGamma _0$. As such, the suction force is modelled in this conceptual approach as

(4.13)\begin{equation} F_S(t)=k_S\rho b \left[\frac{\boldsymbol{Cx}_c(t)+D\varGamma_0(t)}{2{\rm \pi} b}-\frac{b}{2}\dot{\alpha}(t)\right]^2, \end{equation}

where $k_S$ is a coefficient of suction: bookkeeping parameter for suction; its potential-flow value is $k_S=2{\rm \pi}$.

4.3. Dynamics of the separation point

Unsteady separation is a complex phenomenon whose dynamics is quite difficult to model. In fact, a mere criterion for unsteady separation is controversial – the common criterion of vanishing shear is not accurate (Sears Reference Sears1956). However, as stated above, the intent is not to develop a quantitative model for accurate prediction, but rather a conceptual model. A fairly straightforward conceptual model of the dynamics of the separation point is that of Goman & Khrabrov (Reference Goman and Khrabrov1994). Simply, the unsteady separation point $x_s$ follows the steady one $x_0$ (which is a function of the angle of attack) via a simple lag $\tau _1$ and a delay $\tau _2$

(4.14)\begin{equation} \tau_1\dot{x}_s(t) +x_s(t) = x_0\left(\alpha(t)-\tau_2\dot{\alpha}(t)\right). \end{equation}

In this conceptual model, the specific values of $\tau _1$, $\tau _2$ are irrelevant; the functional dependence $x_0(\alpha )$ is arbitrary. Moreover, to relate the dynamics of separation to the lift dynamics, we adopt the Kirchhoff steady model, which relates a nonlinear lift coefficient to separation

(4.15)\begin{equation} C_{L,s}(\alpha)=2{\rm \pi}\sin\alpha\left( \frac{1+\sqrt{x_0(\alpha)}}{2}\right)^2. \end{equation}

That is, given a steady $C_L$-$\alpha$ curve, $C_{L,s}(\alpha )$, the corresponding variation, $x_0(\alpha )$, of the separation point with the angle of attack can be determined from the Kirchhoff model (4.15).

4.4. The final dynamical ROM

Applying the above conceptual modelling efforts to a pitching–plunging wing, any dependence on the angle of attack $\alpha$ would be replaced by the effective angle of attack: $\alpha _{{eff}}=\alpha +\arctan ({\dot {h}}/{U})$, where $h$ is positive downward. As such, the quasi-steady circulation is written as

(4.16)\begin{equation} \varGamma_0(\alpha,\dot{\alpha},\dot{h})=U b C_{L,s}\left(\alpha(t)+\arctan\frac{\dot{h}(t)}{U}\right) + b^2 k_{\dot{\alpha}}\dot{\alpha}(t). \end{equation}

Let us summarize the assumptions behind the developed model:

1. (i) The dynamics (due to vortex shedding) from the quasi-steady circulation $\varGamma _0$ (the input) to the circulatory lift $L_C$ (the output) is linear (4.4), but arbitrary.

2. (ii) The quasi-steady circulation $\varGamma _0$ is given by (4.16) for a pitching–plunging wing, where $C_{L,s}(\alpha )$ is the steady lift curve, which depends arbitrarily (nonlinearly) on the angle of attack $\alpha$.

3. (iii) The non-circulatory lift $L_{NC}$ depends linearly on the wing's normal acceleration at the half-chord point with a simple lag (4.7).

4. (iv) The dynamics of the skin friction drag is neglected; hence, the unsteady drag $\mathcal {D}$ is constructed from the unsteady lift $L$ and suction force $F_S$ (4.10).

5. (v) The unsteady suction force can be determined from the circulatory lift dynamics (4.13).

6. (vi) The Goman–Khrabrov model (4.14) of unsteady separation is adopted; and the steady location $x_0(\alpha )$ of the separation point is determined from the arbitrary steady lift curve $C_{L,s}(\alpha )$ according to the Kirchhoff model (4.15).

Combining all these elements, we write the following ROM of the unsteady nonlinear aerodynamics of a pitching–plunging wing:

(4.17)\begin{align} \frac{\textrm{d}}{\textrm{d}t}\left(\begin{array}{c}\boldsymbol{x}_c \\ x_v \\ x_s \\ \alpha \\ \dot{\alpha} \\ \dot{h}\end{array}\right) &= \left(\begin{array}{c} \boldsymbol{A} \boldsymbol{x}_c+\boldsymbol{B}\varGamma_0(\alpha,\dot{\alpha},\dot{h}) \\ \dfrac{-1}{\tau_v} \left[ x_v+\left(U\cos\alpha-\dot{h}\sin\alpha\right)\dot{\alpha}\right]\\ \dfrac{-1}{\tau_1} \left[ x_s-x_0\left(\alpha+\arctan\dfrac{\dot{h}}{U}-\tau_2\dot{\alpha}\right)\right] \\ \dot{\alpha} \\ 0 \\ 0 \end{array}\right) + \left(\begin{array}{c} \boldsymbol{0}_{n\times 1} \\ \dfrac{ab}{\tau_v} \\ 0 \\ 0\\ 1 \\ 0 \end{array}\right) u_\alpha \nonumber\\ &\quad + \left(\begin{array}{c} \boldsymbol{0}_{n\times 1} \\ \dfrac{\cos\alpha}{\tau_v} \\ 0 \\ 0 \\ 0 \\ 1 \end{array}\right) u_h, \end{align}

which can be written in an abstract form as

(4.18)\begin{equation} \dot{\boldsymbol{x}}(t) = \boldsymbol{f}(\boldsymbol{x}(t)) + \boldsymbol{g}_\alpha(\boldsymbol{x}(t)) u_\alpha(t) + \boldsymbol{g}_h(\boldsymbol{x}(t)) u_h(t), \end{equation}

where $\boldsymbol {x}=[\boldsymbol {x}_c,x_v,x_s,\alpha ,\dot {\alpha },\dot {h}]$ is the state vector of dimension $n+5$, and the control inputs $u_\alpha$, $u_h$ are the pitching and plunging accelerations $\ddot \alpha$, $\ddot {h}$, respectively. The system (4.17) is exactly in the form (2.1), which is amenable to geometric-control analysis. The output equation is written as

(4.19)\begin{equation} \left.\begin{gathered} L(t)=\varPsi_L(\boldsymbol{x}(t))=\rho U \left(\boldsymbol{C}\boldsymbol{x}_c(t)+D\varGamma_0(\alpha,\dot{\alpha},\dot{h})\right)+m_v x_v(t) \cos\alpha(t),\\ \mathcal{D}(t)=\varPsi_D(\boldsymbol{x}(t))=L(t)\tan\alpha(t)-k_S\rho b \left(\frac{\boldsymbol{C}\boldsymbol{x}_c(t)+D\varGamma_0(\alpha,\dot{\alpha},\dot{h})}{2{\rm \pi} b}-\frac{b}{2}\dot{\alpha}(t)\right)^2,\\ x_s(t)=\varPsi_s(\boldsymbol{x}(t))= x_s(t). \end{gathered}\right\} \end{equation}

5.1. A special technique for averaging of high-amplitude periodically forced nonlinear systems

Consider a nonlinear system subjected to a high-frequency, high-amplitude, periodic forcing in the form

(5.4)\begin{equation} \dot{\boldsymbol{x}}(t) = \boldsymbol{f}(\boldsymbol{x}(t)) + \frac{1}{\epsilon} \boldsymbol{G}\left(\boldsymbol{x}(t),\frac{t}{\epsilon}\right). \end{equation}

Theorem A.1 If $\boldsymbol {G}$ is a $T$-periodic, zero-mean vector field and both $\boldsymbol {f}$, $\boldsymbol {G}$ are continuously differentiable, then the averaged dynamical system corresponding to the system (A3a,b) is written as

(5.5)\begin{equation} \dot{\bar{\boldsymbol{x}}}(t)= \bar{\boldsymbol{F}}(\bar{\boldsymbol{x}}(t)), \end{equation}

where $\bar {\boldsymbol {F}} (\bar {\boldsymbol {x}}(t))=({1}/{T}) \int _0^T \boldsymbol {F}(\boldsymbol {x}(t),\tau )\,\textrm {d}\tau$, and $\boldsymbol {F}$ is the pullback of $\boldsymbol {f}$ along the flow $\varPhi _t^{\boldsymbol {G}}$ of the time-varying vector field $\boldsymbol {G}$.

Using the chronological calculus formulation of Agrachev & Gamkrelidze (Reference Agrachev and Gamkrelidze1978), Bullo (Reference Bullo2002) showed that, for a time-invariant $\boldsymbol {f}$ and time-varying $\boldsymbol {G}$, the pullback vector field $\mathcal {\boldsymbol {F}}(\boldsymbol {x}(t),t)$ can be written explicitly in terms of iterated Lie brackets of $\boldsymbol {f}$ and $\boldsymbol {G}$ as

(5.6)\begin{equation} \mathcal{\boldsymbol{F}}(\boldsymbol{x}(t),t)= \boldsymbol{f}(\boldsymbol{x}(t))+ \sum_{k=1}^{\infty} \int_{0}^{t}\cdots\int_{0}^{s_{k-1}} \left(ad_{\boldsymbol{G}(\boldsymbol{x}(t),~s_k)}\dots ad_{\boldsymbol{G}(\boldsymbol{x}(t),s_1)} \boldsymbol{f}(\boldsymbol{x})\right) \textrm{d}s_k\ldots \textrm{d}s_1, \end{equation}

where $ad_{\boldsymbol {G}}\boldsymbol {f}=[\boldsymbol {G}, \boldsymbol {f}]$.

5.2. Application to the unsteady fluid dynamics system

In this section, we apply the above averaging theorem to the developed ROM (4.17), equivalently (4.18) or (5.3). In particular, we use the series (A5) to obtain the pullback vector field $\boldsymbol {F}$ as

(5.7)\begin{equation} \boldsymbol{F}(\boldsymbol{x},t)=\boldsymbol{f}(\boldsymbol{x}) + \boldsymbol{\mathcal{G}}(\boldsymbol{x},t), \end{equation}

where

(5.8)\begin{align} \boldsymbol{\mathcal{G}}(\boldsymbol{x},t) &= \sum_{j=\alpha,h} [\boldsymbol{g}_j,\boldsymbol{f}]\int_0^t u_j(s_1)\,\textrm{d}s_1\notag\\ &\quad + \sum_{j=\alpha,h} \sum_{\ell=\alpha,h} [\boldsymbol{g}_j,[\boldsymbol{g}_\ell,\boldsymbol{f}]]\int_0^t\int_0^{s_1} u_j(s_1)u_\ell(s_2)\,\textrm{d}s_2 \,\textrm{d}s_1+\dots \end{align}

As such, the averaged dynamics is written as

(5.9) \begin{equation} \frac{\textrm{d}}{\textrm{d}t}\left(\begin{array}{c}\bar{\boldsymbol{x}}_c \\ \bar{x}_v \\ \bar{x}_s \\ \bar{\alpha} \\ \bar{\dot{\alpha}} \\ \bar{\dot{h}}\end{array}\right) = \left(\begin{array}{c} \boldsymbol{A} \bar{\boldsymbol{x}}_c+\boldsymbol{B}\varGamma_0(\bar{\alpha},\bar{\dot{\alpha}},\bar{\dot{h}})\\ \dfrac{-1}{\tau_v} \left[ \bar{x}_v+\left(U\cos\bar{\alpha}-\bar{\dot{h}}\sin\bar{\alpha}\right) \bar{\dot{\alpha}}\right]\\[8pt] \dfrac{-1}{\tau_1} \left[ \bar{x}_s-x_0\left(\bar{\alpha}+\arctan \dfrac{\bar{\dot{h}}}{U}-\tau_2\bar{\dot{\alpha}}\right)\right] \\ \bar{\dot{\alpha}} \\ 0 \\ 0 \end{array}\right) + \left(\begin{array}{c} \bar{\boldsymbol{\mathcal{G}}}_{1\to n}(\bar{\alpha}) \\ 0 \\ \bar{\mathcal{G}}_{n+2}(\bar{\alpha}) \\ 0\\ 0 \\ 0 \end{array}\right), \end{equation}

where $\bar {\boldsymbol {\mathcal {G}}}_{1\to n}$ are the first $n$ components of the cycle average of the vector field $\boldsymbol {\mathcal {G}}$, and $\bar {\boldsymbol {\mathcal {G}}}_{n+2}$ is the $(n+2)\textrm {{th}}$-entry. It is found that these components are the only non-zero entries of $\bar {\boldsymbol {\mathcal {G}}}$ and that they depend only on the angle of attack.

Having obtained the averaged dynamics (5.9) of the unsteady fluid dynamical system (4.17), we need to study its equilibrium – the averaging theorem above relates properties of the original nonlinear time-periodic (NLTP) system to equilibrium of the averaged dynamics. That is, we set the left-hand side of (5.9) to zero and solve for $\boldsymbol {x}^*$ that satisfies the equation; i.e. we have

(5.10)\begin{equation} \boldsymbol{f}(\boldsymbol{x}^*) + \bar{\boldsymbol{\mathcal{G}}}(\boldsymbol{x}^*)=0. \end{equation}

Noting that the last three components of $\boldsymbol {\mathcal {G}}$ are zeros, it implies that equilibria of the last two states ($\bar {\dot {\alpha }}$, $\bar {\dot {h}}$) are automatically satisfied; and equilibrium of the $\bar {\alpha }$-dynamics stipulates that $\dot {\alpha }^*=0$, which is intuitively expected. This result, in addition to $\boldsymbol {\mathcal {G}}_{n+1}=0$ implies that equilibrium of $\bar {x}_v$ is attained at $x_v^*=0$. It remains to solve for $\boldsymbol {x}_c^*$ and $x_s^*$. Taking $\dot {h}^*=0$, leaving $\alpha ^*$ arbitrary (to study the effect of the mean angle of attack) and realizing that the first $n$ and the $n+2$ components of $\boldsymbol {\mathcal {G}}$ depend only on $\alpha$, we obtain

(5.11)\begin{equation} \left.\begin{gathered} \boldsymbol{x}_c^* ={-}[\boldsymbol{A}]^{{-}1}\left( \bar{\boldsymbol{\mathcal{G}}}_{1\to n}(\alpha^*)+UbC_{L,s}(\alpha^*)\boldsymbol{B}\right) \\ x_s^* = x_0(\alpha^*)+\tau_1 \mathcal{G}_{n+2}(\alpha^*). \end{gathered}\right\} \end{equation}

Having solved for the equilibrium point $\boldsymbol {x}^*$ of the averaged dynamics (5.9), we perform averaging of the output variables (lift, drag and separation point), which are nonlinear functions of the states $\boldsymbol {x}$. Due to these nonlinearities, the average of an output function $\varPsi _m(\boldsymbol {x})$, $m\in \{L,D,s\}$, is not simply $\varPsi _m(\boldsymbol {x}^*)$. As a first-order approximation, each state $x_i$, $i\in \{1,\ldots ,n+2\}$, can be written as $x_i(t)=x_i^*+A_{x_i}\cos (\omega t+\phi _i)$. Hence, we expand each output function $\varPsi _m$ around $\boldsymbol {x}^*$ in a Taylor series expansion and average over the cycle to obtain

(5.12)\begin{equation} \bar{\varPsi}_m = \varPsi_m(\boldsymbol{x}^*)+\sum_{i=1}^{n+5}\sum_{q=1}^{n+5}\frac{\partial^2\varPsi_m}{\partial x_i \partial x_q} \frac{A_{x_i}A_{x_q}}{2}\cos(\phi_i -\phi_q) + \dots\end{equation}

5.3. Averaged lift

Applying the averaging analysis above to the lift output variable $\varPsi _L$ and normalizing by $\rho U^2 b$, we obtain

(5.13)\begin{equation} \bar{C}_L = C_{L,s}(\alpha^*) + \tfrac{1}{4}\left(H^2 k^2+k_{hf}A_\alpha^2\right) C_{L,s}''(\alpha^*) +O(\epsilon^4;m_v,A_\alpha), \end{equation}

where the infinite series is truncated after $O(\epsilon ^2)$. The result (5.13) is quite revealing. Recall that $\bar {C}_L$ represents the mean of the unsteady lift coefficient due to the pitching–plunging inputs (5.1a,b) and $C_{L,s}(\alpha ^*)$ is the (mean) steady lift coefficient at the mean angle of attack $\alpha ^*$. From a linear system perspective, they may be equal – oscillating the airfoil around a mean angle of attack results in a periodic lift whose mean corresponds to the steady lift at the mean angle of attack. As such, the difference between the two would represent lift enhancement/deficiency. The present analysis implies that, to the leading order, the lift enhancement due to pitching and/or plunging is proportional to the curvature $C_{L,s}''\equiv {\textrm {d}^2 C_{L,s}}/{\textrm {d} \alpha ^2}$ of the steady lift curve. Aerodynamicists typically concern themselves with the slope $C_{L,s}'$ of such a curve; perhaps, this is the first time the role of its curvature has been pointed out. Moreover, the analysis revealed a higher-order pitching contribution to lift enhancement, which is due to added-mass effects: the last term $O(\epsilon ^4;m_v,A_\alpha )$ is proportional to $A_\alpha$ and $m_v$.

The above result is quite interesting in the sense that it identifies the role of the curvature $C_{L,s}''$ of the steady lift curve in controlling lift enhancement due to high-frequency, low-amplitude pitching–plunging. Since $C_{L,s}''$ is almost zero in the small-$\alpha$ linear regime, this result points to potential lift enhancement/deficiency in the stall and post-stall regimes: lift deficiency near the peak(s) in the stall regime; and lift enhancement near the trough in the post-stall regime, as shown in figure 4. In fact, there have been numerous experimental and computational studies on the lift dynamics due to pitching and/or plunging in the stall regime (e.g. Carr Reference Carr1988; Lee & Gerontakos Reference Lee and Gerontakos2004; Andro & Jacquin Reference Andro and Jacquin2009; Rival & Tropea Reference Rival and Tropea2010; Cleaver et al. Reference Cleaver, Wang, Gursul and Visbal2011, Reference Cleaver, Wang and Gursul2012, Reference Cleaver, Wang and Gursul2013; Baik et al. Reference Baik, Bernal, Granlund and Ol2012; Choi et al. Reference Choi, Colonius and Williams2015; Chiereghin et al. Reference Chiereghin, Cleaver and Gursul2019; Gupta & Ansell Reference Gupta and Ansell2019). However, none of these efforts revealed this role of $C_{L,s}''$.

Figure 4. The curvature $C_{L,s}''$ of the steady lift curve controls lift enhancement due to high-frequency, low-amplitude pitching–plunging. Note that the angles of attack at the peaks and troughs of the $C_{L,s}$-$\alpha$ curve change with $R$. The curve reported here is for $R=500$ K for which the angles of attack at the peaks and troughs are $15^\circ$ and $20^\circ$, respectively. For a smaller $R$ (in the range of 10 000–60 000), they would be smaller (approximately$10^\circ$ and $15^\circ$, respectively).

Focusing on the relevant studies corresponding to small-amplitude, high-frequency oscillations around the peaks and troughs of the steady lift curve, we find them corroborating the current result. In fact, as early as the efforts of McCroskey et al. (Reference McCroskey, McAlister, Carr and Pucci1982), one can observe this behaviour. For example, figure 6 of Carr (Reference Carr1988) shows a 4 % lift deficiency due to pitching around the peak. All other cases were either performed at large amplitudes or low frequencies. More recently, the efforts of Rival & Tropea (Reference Rival and Tropea2010) studying plunging of the SD 7003 airfoil at $R=60\,000$ around $\alpha ^*=0^{\circ },5^{\circ },8^{\circ },10^\circ$ also yielded similar results. Inspecting figure 4(a) in their effort, we find little-to-no change between the mean lift coefficient $\bar {C}_L$ and the steady lift coefficient $C_{L,s}(\alpha ^*)$ at the mean angle of attack in the cases of plunging in the linear regime ($\alpha ^*=0,5$); and a lift decrease in the cases of plunging near the peak in the stall regime: 26 % in the $\alpha ^*=8^\circ$ case, and 13 % in the $\alpha ^*=10^\circ$ case. Note that the peak of the steady $C_{L,s}$-$\alpha$ curve is at $10^\circ$ in this case. In general, the angles of attack at the peaks and troughs of the $C_{L,s}$-$\alpha$ curve change with $R$; they increase when $R$ is increased. Also, the results of Cleaver et al. (Reference Cleaver, Wang, Gursul and Visbal2011) are quite clear in this regard; they performed a plunging experiment of NACA 0012 at $R=10\,000$ around $\alpha ^*=15^\circ$, which corresponds to the trough of the $\bar {C}_L$-$\alpha$ curve at this value of $R$ (Chiereghin et al. Reference Chiereghin, Cleaver and Gursul2019). They observed significant lift enhancement in these conditions, which increases with the plunging amplitude, reaching up to 310 % (see figure 11 in their effort). The same group performed a more comprehensive study (Chiereghin et al. Reference Chiereghin, Cleaver and Gursul2019) investigating enhancement/deficiency in the mean and amplitude of the unsteady lift on a plunging NACA 0012 at $R=20\,000$ with different amplitudes, frequencies and mean angles of attack. Investigating figure 7 in their effort and confining ourselves to small amplitudes, one finds: (i) almost no net change in the mean lift coefficient when plunging in the linear regime at $\alpha ^*=0$, $5^\circ$; (ii) a small deficiency in the mean lift coefficient when plunging at $\alpha ^*=9^\circ$, which is the peak of the $\bar {C}_L$-$\alpha$ curve; and (iii) a significant enhancement in the mean lift coefficient when plunging at $\alpha ^*=15^\circ$, which is the trough of the $\bar {C}_L$-$\alpha$ curve. Indeed, these experimental results corroborate the current theoretical finding on the role of the curvature $C_{L,s}''$ of the steady lift curve in controlling lift enhancement, which will be validated in more detail below using relatively high-fidelity computational simulations (URANS).

Finally, it should be noted that the efforts of Cleaver et al. (Reference Cleaver, Wang and Gursul2012, Reference Cleaver, Wang and Gursul2013), although concerned with lift enhancement/deficiency, do not apply here since their observed bifurcations in the mean lift coefficient were at quite high Strouhal numbers (frequencies). Note that increasing the frequency (even with keeping the plunging amplitude fixed) leads to an increase in the effective angle of attack. Therefore, even though the current analysis and results should be valid for high frequencies, there is a limitation on the frequency in the plunging case; for a given plunging amplitude, the maximum frequency should be set to maintain the effective angle of attack within the intended range of $C_{L,s}''$. For example, Cleaver et al. (Reference Cleaver, Wang and Gursul2012) interestingly observed a bifurcation in the mean lift coefficient leading to a significant lift enhancement when plunging at $\alpha ^*=0$. However, this behaviour is not attained until a quite high Strouhal number of 2 (i.e. $k=2{\rm \pi}$) is reached at $H=0.2$, which corresponds to an effective angle of attack in the range $\pm \arctan kH=\pm 51.5^\circ$. Clearly, this is outside the scope of our analysis.

5.4. Averaged drag/thrust

Applying the averaging analysis above to the drag output variable $\varPsi _D$ due to a plunging input, and normalizing by $\rho U^2 b$, we obtain

(5.14)\begin{equation} \bar{C}_D = C_{D,s}(\alpha^*) + \frac{H^2 k^2}{4}\underbrace{\left[ C_{L,s}''\tan\alpha-\frac{k_S}{2{\rm \pi}^2}\left(C_{L,s}C_{L,s}''+k_{hf}^2 C_{L,s}'^2\right)\right]_{\alpha=\alpha^*}}_{-\chi_T} +O(\epsilon^3), \end{equation}

where the infinite series is truncated after $O(\epsilon ^2)$. The result (5.14) is also quite revealing in the sense that any difference between $\bar {C}_D$ and $C_{D,s}(\alpha ^*)$ would represent either a drag increase or decrease (i.e. thrust production). The analytical nature of the current study presents itself here. For example, the last term ($-({k_S k_{hf}^2}/{2{\rm \pi} ^2})C_{L,s}'^2$) between the square brackets is negative indicating thrust production, proportional to the lift curve slope $C_{L,s}'$ (i.e. active in the linear regime), and proportional to the bookkeeping parameter $k_S$ (i.e. due to suction). Therefore, it represents the classical thrust-production mechanism in the linear regime due to plunging. In fact, when using potential-flow characteristics: $k_S=2{\rm \pi}$, $C_{L,s}'=2{\rm \pi}$, and recalling that $k_{hf}$ is the high-frequency gain of the lift transfer function (i.e. Theodorsen's $C(k)$), this term reduces to Garrick's potential-flow result on the mean thrust coefficient $C_T$ due to plunging (Garrick Reference Garrick1937; Taha Reference Taha2018)

(5.15)\begin{equation} C_T={\rm \pi} H^2 k^2 |C(k)|^2. \end{equation}

Therefore, the current result represents a high-frequency generalization of Garrick's classical result to regimes with an arbitrary lift dynamics (not necessarily governed by Theodorsen) and unconventional lift mechanisms (not necessarily the classical $2{\rm \pi} \sin \alpha$). Hence, more contributions are revealed in our result in (5.14) besides Garrick's. In particular, the first two terms are quite interesting as they point to a nonlinear thrust generation mechanism in the stall or post-stall regime where the effect of $C_{L,s}''$ is prominent. The former indicates an increased drag inherited from lift enhancement in regions of positive curvature and vice versa (i.e. an increase in the resultant normal aerodynamic force, similar to the findings of Choi et al. Reference Choi, Colonius and Williams2015). The latter is due to suction and is opposing to the former: it would provide thrust production in regions of positive curvature. Since these two contributions are competing, we performed the following study to assess the role of the curvature $C_{L,s}''$ of the lift curve in drag increase or thrust production.

We performed a URANS simulation (detailed below) of a NACA 0012 at $R=500$ K at different static angles of attack to obtain the mean lift and drag coefficients $C_{L,s}$, $C_{D,s}$, as shown in figure 5(a). Indeed, in the absence of leading edge suction, the resultant force (ignoring skin friction) is normal to the wing surface (Schlichting & Truckenbrodt Reference Schlichting and Truckenbrodt1979), leading to the geometric relation $C_D=C_L\tan \alpha$ between the lift and drag coefficients $C_L$, $C_D$, which is the case for wings making use of a stabilized LEV such as delta wings (Polhamus Reference Polhamus1966) and insects (Ellington et al. Reference Ellington, Van Den Berg, Willmott and Thomas1996; Dickinson et al. Reference Dickinson, Lehmann and Sane1999). Note that it has been usually argued that LEV and leading-edge suction do not coexist (Polhamus Reference Polhamus1966; Ramesh et al. Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014). In fact, the recent concept of the leading-edge suction parameter (LESP) is mainly based on this hypothesis, which was postulated by Ramesh et al. (Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014) as a criterion for the formation of a LEV: an airfoil can sustain leading-edge suction up to a certain limit (depending on the geometry and $R$), and if the LESP exceeds this maximum value, leading-edge separation will occur (i.e. a LEV will form). Therefore, a plausible definition of the suction force coefficient $C_S$ in terms of the computed lift and drag coefficients is

(5.16)\begin{equation} C_S=C_L\tan\alpha-\left(C_D-C_{D,f}\right), \end{equation}

where $C_{D,f}$ is the skin friction drag coefficient. Using the URANS computational results of $C_{L,s}$ and $C_{D,s}$, we use (5.16) to determine the suction force coefficient $C_S$, as shown in figure 5(a).

Figure 5. URANS simulations of NACA 0012 at $R=500$ K at different static angles of attack along with the variation of the thrust control parameter $\chi _T$ due to low-amplitude, high-frequency plunging with the mean angle of attack.

Figure 5(a) shows a large difference between $C_L\tan \alpha$ and $C_D$ at small angles of attack, indicating a large suction force over this range. However, this difference attains a maximum at $\alpha =15^\circ$, supporting the LESP claim: there is a maximum suction force that can be sustained by an airfoil at a given $R$ (Ramesh et al. Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014). Beyond this angle of attack, figure 5(a) shows that the suction force diminishes as the angle of attack increases in accordance with previous findings in literature (Usherwood & Ellington Reference Usherwood and Ellington2002; Dickson & Dickinson Reference Dickson and Dickinson2004; Yan et al. Reference Yan, Taha and Hajj2014). Moreover, at larger angles of attack (above $20^\circ$), indeed the behaviour of $C_L\tan \alpha$ is quite close to $C_D$ with an almost constant shift ($C_{D,f}$), which indicates vanishing of the suction force. These results are quite plausible and intuitive.

Focusing on the change in the mean drag due to plunging, we normalize this change by the effect of kinematics (plunging amplitude) as it appears in (5.14): $H^2 k^2$. In particular, we denote the term between square brackets by $-\chi _T$, and call $\chi _T$ the thrust control parameter: if it is positive, it indicates drag reduction or thrust production; and if it is negative, it indicates drag increase. It simply represents the thrust-production capability due to plunging; if it is larger at condition $A$ than condition $B$, it implies that more thrust would be produced at $A$ than $B$ when oscillating by the same effective-angle-of-attack amplitude: $A_{\alpha _{{eff}}}=\arctan H k$. In fact, the actual thrust coefficient $C_T$ would be given from $\chi _T$ through multiplying by $\tan ^2 A_{\alpha _{{eff}}}/4$.

Normalizing (4.13) by $\rho U^2 b$, we obtain the following relation between $C_S$ and $C_L$:

(5.17)\begin{equation} C_S=\frac{k_S}{4{\rm \pi}^2}C_L^2, \end{equation}

which can be used to estimate the suction parameter $k_S$ from the computed $C_L$ and $C_D$, as shown in figure 5(a). Interestingly, the suction parameter increases beyond its potential-flow value $2{\rm \pi}$, although the suction force itself is less than its potential-flow counterpart $2{\rm \pi} \sin ^2\alpha$. In fact, figure 5(a) implies that the parameter $k_S$ actually varies with $\alpha$; i.e. $k_S=k_S(\alpha )$. Taking this variation into account in our averaging analysis, we obtain more contributions to the thrust control parameter $\chi _T$ besides (5.14), due to the $k_S$-derivatives

(5.18)\begin{equation} \Delta \chi_T = \frac{1}{{\rm \pi}^2}C_{L,s}(\alpha^*)\left[ k_{hf} k_S' C_L'+k_S'' C_L/4 \right]|_{\alpha=\alpha^*} +O(\epsilon^3). \end{equation}

Equations (5.14), (5.18) imply that $\chi _T$ depends on the mean angle of attack $\alpha ^*$; that is, the ability to generate thrust from low-amplitude, high-frequency plunging depends on the mean angle of attack.

Figure 5(b) shows a comparison among several estimates of the $\chi _T$-variation with the angle of attack: (i) the total (combined) $\chi _T$ from (5.14), (5.18), which considers $k_S=k_S(\alpha )$; (ii) the definition (5.14), ignoring the functional variation of $k_S$ with $\alpha$ (i.e. no $k_S$-derivatives) in the averaging analysis; (iii) a no-suction solution: $k_S=0$; and (iv) a full-suction solution considering a constant $k_S=2{\rm \pi}$. The difference between (i) and (ii) is clearly the effects of $k_S$-derivatives; the difference between (ii) and (iv) represents the departure of $k_S$ from the ideal value $2{\rm \pi}$ (mainly in the stall and post-stall regimes); the difference between (ii) and (iii) represents the contribution of suction force; and finally the no-suction solution represents the effect of the curvature $C_{L,s}''$: it is simply $-C_{L,s}''\tan \alpha ^*$. This comparison reveals some interesting points. First, at small angles of attack, $\chi _T>0$, mainly due to suction effects; this is the classical thrust-production mechanism in the linear regime (Garrick Reference Garrick1937). Second, as the angle of attack increases towards stall, there is a strong effect of suction as implied from the small difference between the full-suction solution and the (5.14) definition. However, this strong effect of suction is negative, leading to an increased drag ($\chi _T<0$) rather than thrust production. This behaviour is due to the weak slope $C_{L,s}'$ and strong negative curvature $C_{L,s}''$ of the lift curve near stall, leading to a positive coefficient of $k_S$ in (5.14); i.e. a negative $\chi _T$. In other words, while the suction force is the main thrust generator in the linear regime, it has a negative effect in the stall regime; the no-suction solution is actually thrust producing in this regime (due to the negative curvature $C_{L,s}''$), as shown in figure 5(b). Having said that, it must be pointed out that the effect of $k_S$-derivatives is significant in this range because the suction force increases at a large rate near stall before dropping. Therefore, the net result still ensures a thrust-producing capability: $\chi _T>0$. Computational validation below confirms these findings and shows a thrust-producing capability near stall: at angles of attack ($14^\circ$) quite close to the stall angle ($15^\circ$). Third, right after stall (the suction force stalls slightly after the lift), all the factors contribute in the same direction of a drag increase. The suction force collapses as shown in figure 5(a), which is also observed in the small difference between the no-suction solution and the (5.14) in figure 5(b). Moreover, this rapid decrease of the suction force results in a negative contribution of the $k_S$-derivatives: the total $\chi _T$ is considerably below the (5.14)-definition and the no-suction solution over this range, as shown in figure 5(b). As a result, the thrust-producing capability is annihilated; a huge rise in the mean drag coefficient occurs due to plunging in this range, which is validated below by computational simulations: URANS simulations show that the largest mean drag due to plunging with $k=0.5$, $A_{\alpha _{{eff}}}=5^\circ$ over the range $0\leqslant \alpha ^*\leqslant 20^\circ$ occurs at $\alpha ^*=16,17^\circ$, as will be shown in § 6 below. Finally, figure 5(b) interestingly shows a recovery of $\chi _T$ in the post-stall regime after the negative peak. This fact is also validated below; the URANS simulations indeed show a thrust production at $\alpha ^*=18^\circ$ (see figure 14). The comparison in figure 5(b) implies that this recovery is due to the $k_S$-derivatives; (5.18) suggests that this behaviour is due to the negative slopes of $k_S$ and $C_L$ (the first term).

The effect of pitching on the drag dynamics can be studied in a similar way

(5.19)\begin{align} \bar{C}_D &= C_{D,s}(\alpha^*) + \frac{A_\alpha^2}{4} \left[ 2\frac{C_{L,s}\tan\alpha+k_{hf} C_{L,s}'}{\cos^2\alpha}+k_{hf}C_{L,s}''\tan\alpha+\right.\nonumber\\ &\quad - \frac{k_S}{2{\rm \pi}^2}\left(k_{hf}C_{L,s}C_{L,s}''+k_{hf}^2 C_{L,s}'^2+k^2\left[k_{hf}k_{\dot{\alpha}}(k_{hf}k_{\dot{\alpha}}-2{\rm \pi} )+{\rm \pi}^2\right]\right)\nonumber\\ &\quad - \left. \frac{C_{L,s}}{{\rm \pi}^2} \left(k_{hf} k_S' C_{L,s}'+k_S'' C_{L,s}/4\right) \right]_{\alpha=\alpha^*} +O(\epsilon^3). \end{align}

Equation (5.19) for the pitching effects on drag/thrust production shows similar terms to plunging effects in (5.14). However, unlike plunging, the effect of reduced frequency appears explicitly in the thrust control parameter $\chi _T$ mainly through rotational contributions (proportional to $k_{\dot {\alpha }}$). Therefore, there are non-zero kinematic effects on $\chi _T$ due to pitching, as shown in figure 6: the larger the reduced frequency $k$ and the more backward the pitching axis (smaller $k_{\dot {\alpha }}$), the larger the thrust-producing capability (i.e. $\chi _T$). However, these effects of kinematics are not observed in the post-stall regime where the suction force is diminished. Note that (5.19) implies that the $k$-effects are associated with suction (i.e. proportional to $k_S$). Another notable difference between plunging and pitching effects on thrust production is the well-known fact that, unlike plunging, pure pitching is thrust producing only beyond a certain threshold of frequency (Garrick Reference Garrick1937): figure 6 shows a drag rise over small $\alpha ^*$ when pitching at $k=1$ vs a slight thrust production when pitching at $k=2$. The current results also indicate that the pitching thrust-production capability increases with the mean angle of attack up to stall. That is, pitching is significantly more efficient in thrust production near stall than at small angles of attack. Finally, it should be reported that the coupled pitching–plunging (low-amplitude, high-frequency) effects are of higher orders in $\epsilon$ than the individual effects.

Figure 6. Variation of the thrust control parameter $\chi _T$ due to low-amplitude, high-frequency pitching with the mean angle of attack.

5.5. Averaged separation

Applying the averaging analysis above to the output variable $\varPsi _s$ describing the chord-normalized location $x_s$ of the separation point, we obtain

(5.20)\begin{equation} \bar{x}_s = x_0(\alpha^*) + \frac{k^2}{4}x_0''(\alpha^*)\left[ H^2 +A_\alpha^2\hat{\tau}_2^2\right] +O(\epsilon^4), \end{equation}

where $x_0(\alpha )$ is the steady separation point and $\hat \tau _2={U\tau _2}/{b}$ with $\tau _2$ being the Goman–Khrabrov delay (Goman & Khrabrov Reference Goman and Khrabrov1994), as defined in (4.14). The Kirchhoff model, which relates the steady lift coefficient $C_{L,s}$ to the separation point $x_0$, as given in (4.15), can be used to estimate the variation of the separation point with the angle of attack from a given steady lift curve

(5.21)\begin{equation} x_0(\alpha)=\left(2\sqrt{\frac{C_{L,s}(\alpha)}{2{\rm \pi}\sin\alpha}} - 1 \right)^2. \end{equation}

If the flow is fully attached (i.e. $x_0=1$), then Kirchhoff's model results in the classical potential-flow lift $C_{L,s}=2{\rm \pi} \sin \alpha$; the larger the deviation of $C_{L,s}$ from this classical relation, the smaller $x_0$ (i.e. earlier separation).

Using the URANS $C_{L,s}$ computations, we use the Kirchhoff model (5.21) to estimate the variation of separation point $x_0$ with the angle of attack, as shown in figure 7. The figure also shows the averaged unsteady separation point $\bar {x}_s$ due to plunging with $A_{{eff}}=5^\circ$. The current analysis implies that plunging promotes separation in the pre-stall regime, which is intuitively expected, but provides a significant delay of separation if performed in the post-stall regime near the trough of the $C_{L,s}$-$\alpha$ curve.

Figure 7. Variation of the steady separation point $x_0$ with the mean angle of attack along with the mean unsteady separation point due to plunging with $A_{{eff}}=5^\circ$.

5.6. Added-mass effects

The above analysis did not reveal any added-mass effects up to $O(\epsilon ^2)$ because of the additional lag $\tau _v$; only higher-order effects will be captured. However, if such a lag is relaxed ($\tau _v=0$), the analysis will reveal some important added-mass effects. For example, pure pitching is found to provide lift enhancement only when the pitch axis is ahead of the mid-chord ($a<0$), while a combined pitching–plunging motion can result in lift enhancement independent of the hinge location

(5.22)\begin{equation} \varDelta_{m_v} \bar{C}_L = \frac{\rm \pi}{4}\mu_v k^2 A_\alpha [ H \sin2\alpha^*\sin\phi- 2 a A_\alpha \sin\alpha^*], \end{equation}

where $\varDelta _{m_v}$ denotes an added-mass contribution, $\mu _v$ is the mass ratio ($\mu _v={m_v}/{{\rm \pi} \rho b^2}$: added-mass normalized by its potential-flow value) and $\phi$ is the phase difference between the two harmonic motions, as defined in (5.1a,b). It should be noted that this lift enhancement mechanism due to added mass may not be well known; it is actually a non-intuitive mechanism in the sense that it is purely inertial, so one might think that purely harmonic inertial loads may not have net/averaged effects. However, the current geometric nonlinearities in $\alpha$ promoted such a mechanism. Definitely, it cannot be captured by a purely linear analysis.

Similarly, pure pitching is found to produce thrust when the pitch axis is ahead of the mid-chord ($a<0$), while a combined pitching–plunging motion can be thrust producing at any hinge location if the relative phase is set appropriately

(5.23)\begin{equation} \varDelta_{m_v} \bar{C}_D ={-}\frac{\rm \pi}{2}\mu_v k^2 A_\alpha [ H \cos^2\alpha^*\sin\phi - a A_\alpha \cos\alpha^*]. \end{equation}

Unlike the lift enhancement mechanism due to added mass, this thrust-production mechanism is well known (Garrick Reference Garrick1937). In fact, the second term reduces to Garrick's potential-flow result on the added-mass contribution when taking the limit to small angles of attack ($\cos \alpha ^*\to 1$) and considering potential-flow characteristics ($\mu _v=1$). As shown in figure 8, the added-mass contributions (when pitching at the quarter-chord) may provide significant thrust particularly at high frequencies (it varies quadratically with $k$). Finally, it should be reported that no net added-mass effects due to pure plunging are observed.

Figure 8. Added-mass effects on the thrust control parameter $\chi _T$ due to pitching at the quarter-chord at different reduced frequencies.

6.1. Computational set-up

All simulations are performed using the ANSYS FLUENT commercial package, while meshes are generated using the ICEM CFD software. To handle the airfoil plunging motion, we split the mesh domain (with a C-grid topology) into two different zones: an inner moving circular zone and a static outer one. The inner zone moves with the airfoil as a rigid body, while the outer deforms (using the Dynamic Mesh feature in FLUENT) to accommodate this motion. The inlet semicircle has a 15 chord radius, measured from the quarter chord, while the downstream extends to 30 chords. The inner circular zone radius is 7 chords. The meshes of the two zones are generated using Multi Blocking in ICEM CFD to obtain a two-dimensional quad structured mesh. A grid independence study is performed to reach a final mesh of 102 000 elements with 420 over a sharp trailing edged airfoil, 200 over each of the semi-circle inlet, outlet and upper and lower edges. A boundary layer mesh treatment is implemented with a first layer height of $h=2.15\,\textrm {e}\text {-}{5}$ chord, which is half of the value estimated from Schlichting et al. (Reference Schlichting, Gersten, Krause, Oertel and Mayes1960); this choice ensured $y^+<1$ for proper turbulence modelling.

The simulations are performed using the pressure-based solver. Turbulence is modelled using the Menter shear stress transport $k\text {-}\omega$ SST two-equation complete model (Menter Reference Menter1994). The model possesses a good capability in predicting separation, which is important for the current study. Moreover, it can handle high pressure gradients in near-wall viscous layers (Wilcox Reference Wilcox1998; Pope Reference Pope2000; Wilcox Reference Wilcox2008). The upstream semi-circle and upper and lower edges are set as velocity inlets; and the downstream edge is set as a pressure outlet. Finally, a no-slip wall is set for the airfoil surface. A sliding mesh interface is used between the moving and deforming zone. The airfoil, inner zone and interface edge motion has a rigid body user defined function written to implement the motion. The smoothing dynamic mesh method is used to deform the mesh as springs; and re-meshing is not used because the motion amplitude is relatively small compared with the domain size. The pressure–velocity coupled scheme is adopted. The second-order upwind scheme is selected for spatial discretization and a second-order dual time stepping implicit scheme is adopted for temporal integration with a time step adjusted to ensure 430 steps per motion cycle; and 20 iterations per time step were adopted to ensure convergence in the inlet–outlet net-mass residuals. Finally, simulations are performed until periodicity of the lift coefficient is observed.

6.2. Lift enhancement

The above theoretical analysis, (5.13), revealed the role of the curvature $C_{L,s}''$ of the steady $C_{L,s}$-$\alpha$ curve in controlling enhancement (or deficiency) in the mean lift due to low-amplitude, high-frequency oscillations. To validate this finding, we first constructed the steady $C_{L,s}$-$\alpha$ curve using the URANS computational set-up discussed above, as shown in figure 9(a). Below $\alpha =16^\circ$, no vortex shedding is observed at this Reynolds number. For $\alpha >16^\circ$, the mean values of the steady lift curve are considered for analysis, after smoothing by cubic splines. According to figure 9(a), the curvature is most negative around the peak at $\alpha =15^\circ$ and is most positive around the trough at $\alpha =20^\circ$. So, these two points may be good candidates to test the hypothesis by performing plunging oscillations around these two mean angles of attack with a high frequency ($k=O(1)$) and low amplitude ($H=O(\epsilon )$). There is also another limitation on the plunging amplitude so that the effective angle of attack $\alpha _{{eff}}=\alpha ^* \pm \arctan H k$ stays within the designated range (of signed curvature). Therefore, we performed plunging simulations around $\alpha ^*=0$, 15, $20^\circ$ with $k=0.5$ and $A_{\alpha _{{eff}}}=\arctan H k=5^\circ$, as shown in figure 9. The figure shows 40.6 % enhancement in the mean lift coefficient when plunging around $\alpha ^*=20^\circ$ (i.e. $C_{L,s}''>0)$; 5 % decrease in the mean lift coefficient when plunging around $\alpha ^*=15^\circ$ (i.e. $C_{L,s}''<0$); and almost no change in the mean lift coefficient when plunging around $\alpha ^*=0^\circ$ (i.e. $C_{L,s}''=0$), which is in accordance with the theoretical result from the combined geometric-control-averaging analysis. Indeed, geometric-control theory provides heuristic tools for discovery of nonlinear force-generation and stabilization mechanisms.

Figure 9. URANS steady $C_{L,s}$-$\alpha$ curve for NACA 0012 at $R=500$ K along with plunging simulations around $\alpha ^*=0$, 15, $20^\circ$ with $k=0.5$ and $A_{\alpha _{{eff}}}=\arctan HK=5^\circ$. (a) $C_{L,s}$-$\alpha ^*$ curve. (b) Lift enhancement/deficiency.

It may be interesting to investigate why this happens. Why does plunging around $\alpha ^*=20^\circ$ provide such a significant lift enhancement while plunging around $\alpha ^*=15^\circ$ with the same amplitude and frequency results in a lift deficiency instead? To investigate this point, we show in figure 10 the time history of the unsteady lift coefficient $C_L$ along with the quasi-steady one $C_{L,s}(\alpha _{{eff}})$ for the two cases, and also show snapshots of the vorticity contours over the cycle for the two cases in figures 11, 13. The following list provides a chronological description of the airfoil motion during the cycle:

1. (i) (a) Start at the mid-stroke (i.e. maximum $\dot {h}$ with $\alpha _\textrm {{eff}}=\alpha ^*-A_{\alpha _{{eff}}}$).

2. (ii) (a)–(c) Upward deceleration: until the airfoil reaches the top dead centre after (c) at $t/T=0.25$.

3. (iii) (c)–(e) Downward acceleration: until the airfoil reaches the mid-stroke point again after (e) at $t/T=0.5$.

4. (iv) (e)–(g) Downward deceleration: until the airfoil reaches the bottom dead centre after (g) at $t/T=0.75$.

5. (v) (g)–(i) Upward acceleration: until the airfoil reaches the mid-stroke point after (i) at $t/T=1.0$.

Figure 10. Variations of the unsteady lift coefficients along with the quasi-steady ones over the plunging cycles around $\alpha ^*=15$, $20^\circ$ with $k=0.5$ and $A_{\alpha _{{eff}}}=5^\circ$. (a) $\alpha^{\ast}=20^\circ$ and (b) $\alpha^{\ast}=15^\circ$.

Figure 11. Vorticity contours over the plunging cycle around $\alpha ^*=20^\circ$ with $k=0.5$ and $A_{\alpha _{{eff}}}=5^\circ$.

Figures 11, 13 show that, as the airfoil decelerates upward between (a–c), the flow reattaches. During this phase, there is a trailing-edge vortex (TEV) in the two cases, although significantly stronger in the $20^\circ$ case. This TEV explains why the unsteady lift is less than its quasi-steady counterpart during this phase, as shown in figure 10; this difference is significantly larger in the $20^\circ$ case. It is the classical Wagner wake deficiency effect (Wagner Reference Wagner1925). As the airfoil moves upward, the TEV gets fainter and fainter, consequently, the unsteady lift gets closer to the quasi-steady behaviour. It increases at a faster rate in the $20^\circ$ case. In this case, at some early moment in the downstroke – between (d,e) at $t/T\simeq 0.4$ – a LEV forms, as shown in figure 11. It then grows and moves relatively slowly; it remains over the airfoil for a considerable part of the cycle (approximately a third of the cycle). Therefore, this attached LEV explains the lift enhancement in the $20^\circ$ case: the $C_L$ reaches almost double the quasi-steady value specifically during the period ($t/T=0.4\text {--}0.7$) where the LEV is attached to the airfoil. On the other hand, this strong LEV is not observed in the $15^\circ$ case, as shown in figure 13, which answers the posed question. It should be noted that at these large angles of attack ($15^\circ$ or $20^\circ$), whether static or dynamic, there is a continuous flow of small-scale shear layers separating from the leading edge. However, the current URANS simulation captures only the large coherent structures.

We then pose another question: Why does a strong LEV form in the $20^\circ$ case and not in the $15^\circ$ case? Since the question here is about the formation of a LEV as a large coherent structure (not leading-edge separation in general), the LESP concept (Ramesh et al. Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014) may be a convenient tool in this case. While the LESP concept was introduced for low Reynolds number applications, there are some recent efforts to extend it to high Reynolds number flows (Narsipur et al. Reference Narsipur, Hosangadi, Gopalarathnam and Edwards2016). Unlike the low Reynolds number case, the critical LESP has some kinematic dependence at high Reynolds numbers. In particular, it increases with $k$ until it reaches a maximum value around $k=0.5$. It is found that the critical LESP for NACA 0012 at a relatively high $k$ (${\simeq }0.5$) ranges from 0.32 at $R=30\,000$ up to 0.5 at $R=3\,000\,000$ (Narsipur et al. Reference Narsipur, Hosangadi, Gopalarathnam and Edwards2016). Performing a simple interpolation between the two values, we postulate a 0.35 critical LESP for the current case of $R=500\,000$. Moreover, for pure plunging, one can show that the LESP may be approximated as (Ramesh Reference Ramesh2020)

(6.1)\begin{equation} \mbox{LESP}(t)=\sin\alpha^* + \frac{\dot{h}(t)}{U}C(k), \end{equation}

where $C(k)$ is Theodorsen's lift transfer function. Figure 12 shows the LESP variation over the plunging cycle in the two cases of $\alpha ^*=15$, $20^\circ$ with $k=0.5$ and $A_{\alpha _{{eff}}}=5^\circ$ along with the LESP critical value (0.35). The figure clearly shows that the LESP in the $20^\circ$ case exceeds the critical value around $t/T\simeq 0.31$, indicating the formation of a LEV around this instant, which matches the vorticity snapshots shown in figure 11. On the other hand, the LESP (equivalently the unsteady effective angle of attack) in the $15^\circ$ is not strong enough to initiate a LEV.

Figure 12. LESP variation over the plunging cycle in the two cases of $\alpha ^*=15$, $20^\circ$ with $k=0.5$ and $A_{\alpha _{{eff}}}=5^\circ$ along with the critical LESP for NACA 0012 at $R=500$ K.

Figure 13. Vorticity contours over the plunging cycle around $\alpha ^*=15^\circ$ with $k=0.5$ and $A_{\alpha _{{eff}}}=5^\circ$.

6.3. Drag reduction/thrust generation

The above theoretical analysis, (5.14), (5.18) and figure 5(b), revealed an interesting behaviour of the thrust control parameter $\chi _T$ with the mean angle of attack $\alpha ^*$. It shows a positive $\chi _T$ (i.e. an ability to generate thrust from plunging oscillations) up to the stall angle of attack beyond which there is a significant rise in the mean drag coefficient due to the collapse of the suction force. Moreover, the thrust-producing capability is recovered when increasing the mean angle of attack beyond the stall regime. To validate these findings, we performed plunging simulations at several mean angles of attack ($\alpha ^*=0^\circ , 12^\circ , 14^\circ , 15^\circ , 16^\circ , 17^\circ , 18^\circ$) with $k=0.5$. For each simulation, we compute the mean drag coefficient $\bar {C}_D$ and use it to estimate the thrust control parameter $\chi _T$ as

(6.2)\begin{equation} \chi_T(\alpha^*)=\frac{\bar{C}_D-C_{D,s}(\alpha^*)}{H^2k^2/4}, \end{equation}

where $C_{D,s}(\alpha ^*)$ is the mean drag coefficient at a static angle of attack equal to $\alpha ^*$. Although we do not seek a quantitative validation, the comparison between theoretical predictions (5.14), (5.18) and the URANS computations, shown in figure 14, is encouraging – to capture the drag dynamics in the critical stall regime with such an accuracy using the present simple modelling and analysis is more than satisfactory.

Figure 14. Validation of the thrust control parameter $\chi _T$ due to plunging against URANS computational simulations of a plunging NACA 0012 with $k=0.5$ at $R=500$ K.