2. A quasi-steady model for flapping flight dynamics

2.1. Model of forward flight dynamics

We introduce a quasi-steady dynamical model for the driven heaving of a rigid thin plate which is free to move in the direction perpendicular to heaving. To derive such a model, we adopt the form of the QSM that was developed for the passive free flight of a thin wing (Li et al. Reference Li, Goodwill, Wang and Ristroph2022; Pomerenk & Ristroph Reference Pomerenk and Ristroph2025) and which takes the form of the Newton–Euler equations with forces and torques due to gravity, buoyancy and fluid dynamical effects such as lift, drag, and added mass. The particular case of horizontal motion in a plane and in response to vertical heaving greatly simplifies the aerodynamic considerations by eliminating several degrees of freedom and many of the participating effects. The rotational dynamics need not be evolved because the plate is fixed to lie horizontally in the laboratory frame as in figure 1(a): $\theta =\omega =0$ , where $\theta$ is the plate orientation angle relative to the horizontal and $\omega = \dot {\theta }$ is the angular velocity. This condition simplifies the translational dynamics by eliminating $\omega$ -dependent reference frame terms. The vertical motion is also imposed, rather than evolved, and the position is $y(t) = a\sin (2\pi \textit{ft})$ and velocity $v_y(t) = \dot {y} = 2\pi a f \cos (2\pi \textit{ft})$ for flapping frequency $f$ and amplitude $a$ . The horizontal flight motion is therefore the only output, leading to a two-variable dynamical system for $(x,v_x)$ that is non-autonomous due to the time-dependent input $v_y(t)$ :

(2.1) \begin{align} \begin{aligned} &\dot {x} = v_{x},\\ &\dot {v}_{x}= \frac {\rho \ell }{2(m+m_{11})}\sqrt {v_x^2+v_y^2} \left [C_{\!L}(\alpha ,\textit{Re})v_y - C_{\!D}(\alpha ,\textit{Re})v_x\right ]\!. \end{aligned} \end{align}

Here, $x$ and $y$ are laboratory frame coordinates of the plate centre or any other arbitrary reference point on the body. The plate has chord length $\ell$ and mass $m$ per unit length along the span. The fluid has density $\rho$ , and the added mass coefficient $m_{11}$ (also per unit span) is associated with edgewise motions of the plate (Andersen et al. Reference Andersen, Pesavento and Wang2005b ).

Figure 1. Aerodynamic force coefficients. (a) Definitions of dynamical quantities. Shown here during the downstroke, the plate has horizontal and vertical velocity components $v_x$ and $v_y$ that determine the instantaneous attack angle $\alpha$ and Reynolds number $\textit{Re}$ . The total aerodynamic force is resolved into lift and drag. Panels (b) and (c) show colour maps of lift and drag coefficients $C_{\!L}(\alpha ,\textit{Re})$ and $C_{\!D}(\alpha ,\textit{Re})$ . Panels (d) and (e) show transects of lift and drag coefficients at fixed $\textit{Re}=10^3$ . Dashed magenta curves and the accompanying equations represent the mathematical forms applicable to the attached flow regime at small $\alpha$ and the separated flow regime at high $\alpha$ . The grey band shows the transitional region where stall occurs. Dashed blue lines and accompanying labels highlight some parameter values. Open black circles indicate experimental measurements from Li et al. (Reference Li, Goodwill, Wang and Ristroph2022). ( f) The lift-to-drag ratio $C_{\!L}/C_{\!D}$ with respect to $\alpha$ at logarithmically spaced $\textit{Re}$ . The graph of $\cot \alpha$ (black dashed curve) and the marked intersections relate to an equilibrium analysis.

Notably, the dominant aerodynamic effect that remains is the force associated with translation, which is decomposed into lift and drag components and expressed in the usual high- $\textit{Re}$ form with force coefficients $C_{\!L}$ and $C_{\!D}$ . Figure 1(a) shows a force diagram at an arbitrary instant. The force coefficients are in general expected to vary with both the attack angle $\alpha =\arctan (v_y/v_x)$ and Reynolds number $\textit{Re}= \rho \ell \sqrt {v_x^2+v_y^2} / \mu$ , where $\mu$ is the fluid viscosity coefficient. These quantities vary in time and, due to the involvement of $v_x$ , are outcomes of the solution. The mathematical forms for the coefficients $C_{\!L}(\alpha ,\textit{Re})$ and $C_{\!D}(\alpha ,\textit{Re})$ dictate the behaviour of the system and thus whether such a model reproduces the key dynamical features of flapping locomotion. If at all successful, one expects such a model to reproduce forces averaged over times longer than a wing stroke period $1/f$ and hence aspects of the flight dynamics over similarly long time scales. As a quasi-steady treatment, many effects are treated only approximately through the specified force coefficients, and many are neglected altogether. The key assumption that the instantaneous force during unsteady motions is taken as the time-averaged force for the corresponding steady condition (i.e. same speed and attack angle) is not rigorously justified a priori and hence must be assessed by comparing model outputs with known results. Intrinsically unsteady phenomena such as fluctuating forces from vortex shedding are not included (Wang Reference Wang2000).

2.2. Dimensionless form of the model and its outputs

Comparison across different parameter values is facilitated by non-dimensionalisation of (2.1). Following earlier works, we choose characteristic scales for length $\ell$ , time $1/f$ , and thus speed $\ell f$ (Alben & Shelley Reference Alben and Shelley2005; Spagnolie et al. Reference Spagnolie, Moret, Shelley and Zhang2010). Let $v_x' = v_x/(\ell f)$ , $x'=x/\ell$ , likewise for $v'_y$ and $y'$ , and $t'=\textit{tf}$ . Then the imposed driving takes the form $y' = A\sin (2\pi t')$ and $v_y' = \dot {y'}$ for the dimensionless amplitude $A=a/\ell$ . The dynamics transform as

(2.2) \begin{align} \begin{aligned} \dot {x'} &= v_x' ,\\ \dot {v}_x' &= \frac {1}{2M}\sqrt {(v_x')^2 + (v_y')^2}\big [C_{\!L}(\alpha ,\textit{Re})v_y'-C_{\!D}(\alpha ,\textit{Re})v_x'\big ], \end{aligned} \end{align}

where $M=(m+m_{11})/\rho \ell ^2$ is the dimensionless mass that normalises the combined solid and added masses (per unit span) by a relevant mass of fluid.

The coefficients $C_{\!L}$ and $C_{\!D}$ are themselves dimensionless and expressed in terms of dimensionless quantities. The angle of attack $\alpha$ assumes the same form whether dimensional or dimensionless. The dynamical or instantaneous Reynolds number transforms as $\textit{Re}= \rho \ell \sqrt {v_x^2+v_y^2} / \mu = \rho \ell ^2 f \sqrt {v_x^{\prime 2}+v_y^{\prime 2}} / \mu = (\textit{Re}_{\!f} / A) \sqrt {v_x^{\prime 2}+v_y^{\prime 2}}$ , where the flapping Reynolds number is $\textit{Re}_{\!f} = \rho \ell a f / \mu$ . The system is therefore characterised by three dimensionless input parameters: $M=(m+m_{11})/\rho \ell ^2$ , $A=a/\ell$ , and $\textit{Re}_{\!f} = \rho \ell a f / \mu$ . The outputs are the dynamics $(x',v_x') = (x/\ell ,v_x/\ell f)$ , the results of which will be analysed to extract summary quantities.

2.3. Specification of aerodynamic force coefficients

Lift and drag coefficients for thin rectangular plates have been reported for certain ranges of attack angles and isolated Reynolds numbers (Tomotika & Aoi Reference Tomotika and Aoi1953; Sane & Dickinson Reference Sane and Dickinson2001; Andersen et al. Reference Andersen, Pesavento and Wang2005a ; Ehrenstein & Eloy Reference Ehrenstein and Eloy2013; Bhati et al. Reference Bhati, Sawanni, Kulkarni and Bhardwaj2018; Li et al. Reference Li, Goodwill, Wang and Ristroph2022). Based on these partial characterisations, we propose mathematical forms for $C_{\!L}(\alpha ,\textit{Re})$ and $C_{\!D}(\alpha ,\textit{Re})$ that capture the known flow phenomenology and which involve constant parameters whose values can be informed by reported aerodynamic data. These forms are:

(2.3) \begin{align} &C_{L,D}(\alpha ,\textit{Re}) = f(\alpha ,\textit{Re}) C^{A}_{L,D}(\alpha ,\textit{Re}) + \left [1-f(\alpha ,\textit{Re})\right ] C^{S}_{L,D}(\alpha ,\textit{Re}) \nonumber\\ & \quad \quad \textrm {where} \quad f(\alpha ,\textit{Re}) = \frac {1}{2}\left [1-\tanh \frac {\alpha -\alpha _S(\textit{Re})}{\delta _S(\textit{Re})}\right ] \quad \textrm {with} \quad \alpha _S = 15^\circ, \quad \delta _S = 5^\circ. \nonumber\\ &\textrm {Attached flow coefficients:} \nonumber\\ & \quad C^A_L(\alpha ,\textit{Re}) = C^I_L(\textit{Re}) \sin \alpha \quad \textrm {with} \quad C^I_L = 5 \nonumber\\ & \quad C^A_D(\alpha ,\textit{Re}) = C^0_D(\textit{Re}) + C^I_D(\textit{Re}) \sin ^2 \alpha \quad \textrm {with} \quad C^I_D = 5 \nonumber\\ & \quad \quad \quad \textrm {where} \quad C^0_D(\textit{Re}) = C^{0,\textit{FT}}_D + \frac {C^{0,\textit{SF}}_D}{\sqrt {\textit{Re}}} \quad \textrm {with} \quad C^{0,\textit{FT}}_D = 0.1, \quad C^{0,\textit{SF}}_D = 1.3. \nonumber\\ &\textrm {Separated flow coefficients:} \nonumber\\ & \quad C^S_L(\textit{Re}) = C^{\pi /4}_L(\textit{Re}) \sin 2\alpha \quad \textrm {with} \quad C^{\pi /4}_L = 1 \nonumber\\ & \quad C^S_D(\textit{Re}) = C^{\pi /2}_D(\textit{Re}) \sin ^2 \alpha \nonumber\\ &\quad \quad \textrm {where} \quad C^{\pi /2}_D(\textit{Re}) = C^{\pi /2,\textit{HR}}_D+\frac {C^{\pi /2,\textit{LR}}_D}{\textit{Re}} \quad \textrm {with} \quad C^{\pi /2,\textit{HR}}_D=1.8, \quad C^{\pi /2,\textit{LR}}_D=20. \end{align}

These equations represent surfaces $C_{\!L}(\alpha ,\textit{Re})$ and $C_{\!D}(\alpha ,\textit{Re})$ that are plotted in figures 1(b) and 1(c). The attack angle $\alpha \in [0,\pi /2]=[0,90^\circ ]$ is shown over the conventional aerodynamic range, and the coefficients are readily extended over all possible orientations $\alpha \in [0,2\pi )$ by consideration of the dual symmetries of the plate (Li et al. Reference Li, Goodwill, Wang and Ristroph2022). The displayed range of the Reynolds number $\textit{Re} \in [10,10^5]$ reflects the intermediate regime over which the model can be expected to apply. The lift map is strictly uniform with respect to $\textit{Re}$ , and the drag map is largely so except for notable variations at lower values. The forces vary with $\alpha$ in ways that are more clearly seen in the transects of figures 1(d) and 1(e) at a representative $\textit{Re}=10^3$ , and these panels graphically interpret the terms and parameters in the model equations. Comparison with the experimental measurements (open circles) of Li et al. (Reference Li, Goodwill, Wang and Ristroph2022) shows good agreement in the forms of the curves while not faithfully reproducing all features of the data. Note that these experiments reflect only pressure-based or normal forces and hence the measured $C_{\!D}$ at low angles $\alpha \lt \alpha _S$ (faded data points) are underestimates that neglect skin friction as a force tangential to the plate (Li et al. Reference Li, Goodwill, Wang and Ristroph2022). Additional checks against experimental and computational studies at $\textit{Re}=10^2$ (Wang et al. Reference Wang, Birch and Dickinson2004) and $\textit{Re}=10^5$ (Ananda, Sukumar & Selig Reference Ananda, Sukumar and Selig2015) indicate good agreement for low angles and generally at the level of trends in the curves and rough magnitudes.

The equations (2.3) above are structured to reflect the program we have used for specifying the coefficient formulas. First, we assume two distinct aerodynamic regimes which are respectively associated with attached flow for low $\alpha$ and fully separated flow at high $\alpha$ . The function $f(\alpha ,\textit{Re})$ serves the role of selecting the regime. It is the sigmoid-shaped logistic curve that takes on values near 1 for low $\alpha$ , drops near $\alpha _S$ , and approaches 0 for high $\alpha$ . The critical angle $\alpha _S$ marks the stall transition, and the factor $\delta _S$ specifies the angular width or range over which the transition occurs. Both may in principle depend on $\textit{Re}$ , but available data at $\textit{Re}=10^3$ and $10^5$ show stall to occur over the range $\alpha = 10^\circ$ to $20^\circ$ , hence our choices for the constant values $\alpha _S = 15^\circ$ and $\delta _S = 5^\circ$ (Ananda et al. Reference Ananda, Sukumar and Selig2015; Li et al. Reference Li, Goodwill, Wang and Ristroph2022). The transitional region is marked by grey stripes in figures 1(d) and 1(e).

Next, we set about specifying the coefficients $C^{A}_{L,D}(\alpha ,\textit{Re})$ for the attached flow regime, focusing first on the dependence on $\alpha$ then $\textit{Re}$ . The lift coefficient $C^A_L \propto \sin \alpha$ follows the form of Kutta–Joukowski theory for thin wings (Anderson Reference Anderson2011) and is consistent with measurements at various $\textit{Re}$ showing approximately linear increase in lift with attack angle for low $\alpha$ (Ananda et al. Reference Ananda, Sukumar and Selig2015; Li et al. Reference Li, Goodwill, Wang and Ristroph2022). The prefactor $C^I_L(\textit{Re})$ is the slope or inclination of the lift curve at small $\alpha$ , and it may in principle depend on $\textit{Re}$ . Its value is $2\pi$ in potential flow theory (Anderson Reference Anderson2011; Gutierrez-Castillo et al. Reference Gutierrez-Castillo, Aguilar-Cabello, Alcalde-Morales, Parras and del Pino2021), and experimental studies across widely ranging $\textit{Re} =10^2$ to $10^5$ report somewhat lower values between 4 and 6 with no systematic dependence on $\textit{Re}$ (Ananda et al. Reference Ananda, Sukumar and Selig2015; Gutierrez-Castillo et al. Reference Gutierrez-Castillo, Aguilar-Cabello, Alcalde-Morales, Parras and del Pino2021; Li et al. Reference Li, Goodwill, Wang and Ristroph2022). Hence we choose the constant $C^I_L(\textit{Re})=5$ as a representative value across all $\textit{Re}$ . The drag model $C^{A}_{D}(\alpha ,\textit{Re})$ has two terms, the first of which represents the fluid resistance present for edgewise motions at $\alpha =0$ . It involves the constant term that represents the (small) pressure drag associated with the finite thickness of the plate (Anderson & Bowden Reference Anderson and Bowden2005; Andersen et al. Reference Andersen, Pesavento and Wang2005b ). It is expected to be of the same order as the thickness-to-chord ratio, and the value $C^{0,\textit{FT}}_D = 0.1$ is chosen as typical of wing plates considered in previous studies (Taylor et al. Reference Taylor, Nudds and Thomas2003; Vandenberghe et al. Reference Vandenberghe, Zhang and Childress2004; Andersen et al. Reference Andersen, Pesavento and Wang2005a ; Alben & Shelley Reference Alben and Shelley2005; Li et al. Reference Li, Goodwill, Wang and Ristroph2022). The second term captures the effect of skin friction, and the prefactor value $C^{0,\textit{SF}}_D = 1.3$ and scaling as $1/\sqrt {\textit{Re}}$ are the classical results of the Blasius boundary-layer theory calculation (Tomotika & Aoi Reference Tomotika and Aoi1953; Bhati et al. Reference Bhati, Sawanni, Kulkarni and Bhardwaj2018). The second term in $C_{\!D}^A(\alpha ,\textit{Re})$ represents lift-induced drag, which arises due to the change in effective attack angle due to lift generation (Anderson & Bowden Reference Anderson and Bowden2005; Anderson Reference Anderson2011). Its form as $\sin ^2 \alpha$ reflects the quadratic dependence on the lift coefficient according to Prandtl’s lifting-line theory (Anderson & Bowden Reference Anderson and Bowden2005; Gutierrez-Castillo et al. Reference Gutierrez-Castillo, Aguilar-Cabello, Alcalde-Morales, Parras and del Pino2021), which also predicts the constant prefactor to depend on the span-to-chord aspect ratio. The chosen value $C_{\!D}^I=5$ was determined empirically in Li et al. (Reference Li, Goodwill, Wang and Ristroph2022) for plates with aspect ratios of 5 to 10.

We carry out similar procedures for the coefficients $C^{S}_{L,D}(\alpha ,\textit{Re})$ in the separated flow regime applicable to high angles beyond the stall transition. Considering lift, the dependence on $\alpha$ of $C^{S}_{L}\propto \sin 2\alpha$ has been adopted in previous modelling work and reflects the observation that lift tends to be maximal near $\alpha = \pi /4 = 45^\circ$ and fall off toward lower and higher values, as expected by symmetry (Andersen et al. Reference Andersen, Pesavento and Wang2005b ; Li et al. Reference Li, Goodwill, Wang and Ristroph2022). The maximal value of the lift coefficient has been reported to be 0.7 to 1.3 in studies across widely ranging $\textit{Re}$ , motivating our choice of the constant $C_{\!L}^{\pi /4}=1$ (Andersen et al. Reference Andersen, Pesavento and Wang2005b ; Ananda et al. Reference Ananda, Sukumar and Selig2015; Li et al. Reference Li, Goodwill, Wang and Ristroph2022). Considering drag, the dependence $C^{S}_{D} \propto \sin ^2 \alpha$ has similarly been used previously (Andersen et al. Reference Andersen, Pesavento and Wang2005b ; Li et al. Reference Li, Goodwill, Wang and Ristroph2022) and models the observed monotonic increase up to $\alpha = \pi /2 = 90^\circ$ . The maximal value of $C^{\pi /2}_D(\textit{Re})$ is the drag coefficient for a plate in normal flow, which for bluff bodies is constant at higher $\textit{Re}$ and scales as $1/\textit{Re}$ at low $\textit{Re}$ (Jones & Knudsen Reference Jones and Knudsen1961; Bhati et al. Reference Bhati, Sawanni, Kulkarni and Bhardwaj2018). The associated values $C^{\pi /2,\textit{HR}}=1.8$ and $C^{\pi /2,\textit{LR}}=20$ are chosen as fits to measurements across $\textit{Re} \in (10,10^8)$ (Bhati et al. Reference Bhati, Sawanni, Kulkarni and Bhardwaj2018).

In summary, the lift and drag coefficients $C_{L,D}(\alpha ,\textit{Re})$ have functional dependencies on attack angle and Reynolds number that are informed by aerodynamic theory. They involve 9 constants whose values are informed by previous measurements: $\alpha _S = 15^\circ$ , $\delta _S = 5^\circ$ , $C^I_L = 5$ , $C^{0,\textit{FT}}_D = 0.1$ , $C_{\!D}^{0,\textit{SF}} = 1.3$ , $C^I_D = 5$ , $C^{\pi /4}_L = 1$ , $C^{\pi /2,\textit{HR}}_D = 1.8$ , $C^{\pi /2,\textit{LR}}_D = 20$ . The first two define the stall transition, the following four define the attached flow regime for small $\alpha$ , and the last three define the separated flow regime for large $\alpha$ .

3. Characterisation of flight dynamical behaviours

Here, we characterise the dynamical behaviours resulting from numerical solutions of the proposed QSM and thereby demonstrate that it closely reproduces the key characteristics of flapping-based propulsion. To this end, equations (2.2) are integrated numerically in time via MATLAB’s built-in solver ode15s that is optimised for systems of stiff ordinary differential equations. Its algorithm is based on numerical differentiation formulas up to fifth order, and variable time stepping ensures a specified degree of accuracy. Convergence tests determine the values of the relative $\textrm {RelTol}=10^{-10}$ and absolute $\textrm {AbsTol}=10^{-9}$ tolerances used for all results presented here. The initial condition for the translational speed is $v_x'(0)=0.01$ , which represents a small perturbation whose purpose is to break the left–right symmetry. We choose $M=A=1$ for simplicity here, and the effects of varying these parameters will be systematically assessed in the next section.

Figure 2. Model dynamics reproduce characteristic flapping flight behaviours. (a) Dimensionless position $x'$ versus $t'$ for the indicated values of $\textit{Re}_{\!f}$ . (b) Dimensionless horizontal speed against time. Inset: log–linear plot for early times. (c) Terminal flight speed versus flapping frequency, shown in dimensionless forms as the horizontal Reynolds number $\textit{Re}_{U}$ against $\textit{Re}_{\!f}$ . Low $\textit{Re}_{\!f}$ leads to a stationary state, while higher $\textit{Re}_{\!f}$ yields forward flight following a linear relation $\textit{Re}_{U} \propto \textit{Re}_{\!f}$ . (d) Magnified view for low $\textit{Re}_{\!f}$ showing hysteresis and bistability at the transition to forward flight. Stability of the stationary state is lost at the critical value $\textit{Re}_{\!f}^*=25$ . (e) Strouhal number $\textit{St}$ plotted against $\textit{Re}_{\!f}$ , with the value $\textit{St}^* = 0.2$ characterising $\textit{Re}_{\!f}\gg \textit{Re}_{\!f}^*$ . ( f) Dimensionless exponential time scale $\tau '$ plotted against $\textit{Re}_{\!f}$ , where the sign of $\tau '$ indicates decay $(-)$ or growth $(+)$ of speed and hence stability or instability of the stationary state. The change of sign indicates the state transition at $\textit{Re}_{\!f}^*=25$ , and the limiting value $\tau '^*=2.5$ characterises $\textit{Re}_{\!f}\gg \textit{Re}_{\!f}^*$ .

Figure 2(a) displays curves of the dimensionless displacement $x' = x/\ell$ over time for logarithmically increasing $\textit{Re}_{\!f}=10,10^2,\dots ,10^5$ , where each curve is associated with a different fixed $\textit{Re}_{\!f}$ that increases from left to right. The initial positions are spaced apart for ease of visibility. For small $\textit{Re}_{\!f}=10$ (left vertical line, cyan), the plate remains stationary for all time. In contrast, higher $\textit{Re}_{\!f} \geqslant 100$ (right curves, darker cyan and magenta) induces the plate to take off into forward flight. Figure 2(b) presents these same simulation results in terms of dimensionless speeds $v_x'$ over time for the same choices of $\textit{Re}_{\!f}$ . Evidently, flapping of a plate induces one of two possible terminal or long-time behaviours determined by the value of the flapping Reynolds number. For low $\textit{Re}_{\!f}$ , the plate assumes a stationary state in which it flaps in place with no forward progression. For moderate to high $\textit{Re}_{\!f}$ , it reaches a forward flight state characterised by constant stroke-averaged translational speed with relatively small fluctuations within each stroke. The inset of figure 2(b) is a rescaled log–linear plot of the early-time data, from which it is evident that the speed either decays (negatively sloped cyan curve, $\textit{Re}_{\!f} = 10$ ) or grows (positively sloped curves, $\textit{Re}_{\!f} = 10^2$ to $10^5$ ) exponentially from its initial value.

Further details regarding the terminal states and the transition are provided in figure 2(c) and the magnified view of panel (d). Here the mean speed $U = \overline {v}_x = \overline {v}_x' \ell f$ during the terminal state is used to compute the forward flight Reynolds number $\textit{Re}_{U}=\rho U\ell /\mu$ based on the horizontal motion. Panel (c) reports the output motion $\textit{Re}_U$ versus the input flapping Reynolds number $\textit{Re}_{\!f}$ , revealing a linear relationship for sufficiently large $\textit{Re}_{\!f}$ . The more complex behaviour at low $\textit{Re}_{\!f}$ is clarified by the zoomed-in view in (d), which reveals a hysteresis loop that brackets a range of bistability in which the stationary and forward flight states coexist. This structure is revealed by a procedure in which we incrementally increase $\textit{Re}_{\!f}$ in a stepwise fashion, where each step is initialised with the previous step’s flight speed and $\textit{Re}_{\!f}$ is held constant for sufficiently long that the forward dynamics equilibrate and $\textit{Re}_U$ can be extracted. The results of this upward sweep are shown as the magenta curve, and it is followed by a downward sweep shown in cyan. These data identify a lower value $\textit{Re}_{\!f} = 18$ below which only the stationary state is observed and a higher value of $25$ above which only the forward flight state is seen. As summarised in table 1, these results can be compared with the hysteretic dynamics reported in previous experiments to occur around $\textit{Re}_{\!f} = 20$ to 55 (Vandenberghe et al. Reference Vandenberghe, Zhang and Childress2004), as well as the more complex transition to locomotion seen in numerical simulations around $\textit{Re}_{\!f} = 15$ to $50$ (Alben & Shelley Reference Alben and Shelley2005).

Table 1. Comparison of predictions from the current study for the quantities $\textit{Re}_{\!f}^*$ and $\textit{St}^*$ against previously reported values from experiments and direct numerical simulations.

These same data are recast in figure 2(e) in terms of the Strouhal number $\textit{St} = 2af/U = 2\textit{Re}_{\!f}/\textit{Re}_U$ attained across varying $\textit{Re}_{\!f}$ . Beyond the onset to locomotion, $\textit{St}$ rapidly drops to approximately 0.3 and thereafter changes little as it asymptotically approaches a constant value near $\textit{St}^* = 0.2$ for large $\textit{Re}_{\!f}$ , where $\textit{St}^* = \textit{St}(\textit{Re}_{\!f} \gg \textit{Re}_{\!f}^*)$ . As summarised in table 1, these results can be compared with the range $\textit{St} = 0.15$ to 0.4 that has been reported to apply to the cruising flight and steady swimming of many diverse animals and which has been associated with efficient flapping locomotion (Taylor et al. Reference Taylor, Nudds and Thomas2003; Nudds et al. Reference Nudds, Taylor and Thomas2004; Alben & Shelley Reference Alben and Shelley2005).

Lastly, we characterise the time scale associated with the observed exponential decay or growth of the flight speed at early times. Building on the analysis shown in the inset of figure 2(b), we fit lines to the log–linear data of $v_x'(t')$ across values of $\textit{Re}_{\!f}$ and associate the inverse of the slope with the decay/growth time scale $\tau '$ . Negative values $\tau '\lt 0$ are associated with decay and hence stability of the stationary state, and positive $\tau '\gt 0$ with growth and thus instability that drives the body into forward flight. The extracted data $\tau '(\textit{Re}_{\!f})$ are shown as the curve in figure 2( f). These results identify the critical value $\textit{Re}_{\!f}^*=25$ (vertical dotted line) below which $\tau '\lt 0$ and above which $\tau '\gt 0$ . This value matches to the upper bound for the hysteretic region as detailed in figure 2(d), and this correspondence is understood by noting that both results pertain to the loss of stability of the stationary state. These findings indicate that the sign change of $\tau '$ can be used generally to identify the transitional value $\textit{Re}_{\!f}^*$ that is shown here to be 25 for $M=A=1$ and whose characterisation across the parameter space is taken up in the next section. Further, for large $\textit{Re}_{\!f} \gg \textit{Re}_{\!f}^*$ the time scale asymptotically approaches a constant value $\tau '^*=\tau '(\textit{Re}_{\!f} \gg \textit{Re}_{\!f}^*)$ , and its parametric dependence will be assessed in what follows. We are not aware of previous studies characterising this time scale, and hence the results presented here are predictions that may be compared with the results of future experiments and simulations.

4. Characterisations across physical and kinematic parameters

Having confirmed the model’s ability to reproduce the key phenomena for the intermediate values $M=A=1$ , we next conduct a sweep throughout the full space of the physical and kinematic input parameters $M$ , $A$ and $\textit{Re}_{\!f}$ . The output quantities of interest include: the terminal Strouhal number $\textit{St}$ attained, the characteristic growth or decay time scale $\tau '$ , and the critical bifurcation parameter $\textit{Re}_{\!f}^*$ that marks the transition from the stationary state to forward flight. It will also be of interest to consider the high- $\textit{Re}_{\!f}$ asymptotic values of the former two parameters: $\textit{St}^* = \textit{St}(\textit{Re}_{\!f} \gg \textit{Re}_{\!f}^*)$ and $\tau '^* = \tau '(\textit{Re}_{\!f} \gg \textit{Re}_{\!f}^*)$ . The results are compiled in figure 3. The explored values $M\in \lbrace 0.1,1,10\rbrace$ , $A\in [0.1,10]$ and $\textit{Re}_{\!f}\in [10,10^5]$ are chosen to span the wide-ranging parameters applicable to flapping-based swimming and flight of animals (Pennycuick Reference Pennycuick1996; Taylor et al. Reference Taylor, Nudds and Thomas2003; Vandenberghe et al. Reference Vandenberghe, Zhang and Childress2004; Fish et al. Reference Fish, Schreiber, Moored, Liu, Dong and Bart-Smith2016; Kang et al. Reference Kang, Cranford, Sridhar, Kodali, Landrum and Slegers2018).

Displayed in the panels of figure 3(a) are terminal values of $\textit{St}$ , where each panel scans the space of $(\textit{Re}_{\!f},A)$ and the three panels correspond to $M = 0.1$ , 1 and 10. For lower values of $\textit{Re}_{\!f}$ below the transition, the stationary state prevails and $\textit{St} \rightarrow \infty$ (dark blue) across all values of $M$ and $A$ . For higher values of $\textit{Re}_{\!f}$ beyond the transition, the Strouhal number abruptly reaches finite values and saturates to the selected value of approximately $\textit{St}^* = 0.2$ (yellow-orange). In this regime, $\textit{St}$ shows only minor variations with $M$ and $A$ . These results seem to corroborate the view that forward flapping flight dynamics are universal in the sense that this key metric is conserved across the vast majority of the parameter space (Taylor et al. Reference Taylor, Nudds and Thomas2003; Nudds et al. Reference Nudds, Taylor and Thomas2004; Alben & Shelley Reference Alben and Shelley2005; Ramananarivo et al. Reference Ramananarivo, Fang, Oza, Zhang and Ristroph2016).

Figure 3. Tableau characterising the dependencies of the key output quantities on the input parameters $M$ , $A$ , and $\textit{Re}_{\!f}$ . Each pixel in a given map corresponds to a run of the simulation at the indicated parameter values. Panels (a) and (b) show the Strouhal number $\textit{St}$ and characteristic growth/decay time $\tau '$ . The former is nearly constant with value $\textit{St}^*=0.2$ for large $\textit{Re}_{\!f}$ and across all $M$ and $A$ . (c) The critical value $\textit{Re}_{\!f}^*$ , defined as the point at which $\tau '$ changes sign. This is a nearly constant field in $M$ and $A$ with $\textit{Re}_{\!f}^*=25$ . (d) The rescaled growth/decay time scale $A\tau '^*/M$ for fixed $\textit{Re}_{\!f}=10^5\gg \textit{Re}_{\!f}^*$ displays a nearly uniform map, implying that $\tau '^* \approx 2.5 M/A$ .

Displayed in the panels of figure 3(b) are the dimensionless growth/decay times $\tau '$ . Across all $M$ and $A$ , it is seen that $\tau '\lt 0$ (blue tones) for sufficiently low $\textit{Re}_{\!f} \lt \textit{Re}_{\!f}^*$ and $\tau '\gt 0$ (red tones) for higher $\textit{Re}_{\!f} \gt \textit{Re}_{\!f}^*$ . Thus, the existence of the transition from the stationary state to locomotion is maintained for all conditions, and the transitional value $\textit{Re}_{\!f}^*$ itself appears similarly preserved. To assess this quantitatively and systematically across the space, we extract $\textit{Re}_{\!f}^*(M,A)$ via the change of sign in $\tau '$ and compile the results in the map of figure 3(c). The panel is notably uniform in colour, indicating that the value of the critical flapping Reynolds number is highly conserved near $\textit{Re}_{\!f}^* = 25$ . Additional investigations indicate that the value $\textit{Re}_{\!f}=18$ marking the lower boundary of the hysteresis window is also conserved across $M$ and $A$ .

Further, comparisons within each panel of figure 3(b) show that $|\tau '|$ generally decreases (lighter tones) with $A$ , and comparisons across the panels show that $|\tau '|$ increases (darker tones) with $M$ . This suggests a relationship of the form $\tau '(M,A) \sim M/A$ , and indeed rescaling as $A\tau '/M$ proves to largely remove the variations, as shown by the map of figure 3(d). Here, we fix $\textit{Re}_{\!f} = 10^5$ as representative of $\textit{Re}_{\!f} \gg \textit{Re}_{\!f}^*$ , and it is seen that $A\tau '^*/M$ falls within a narrow range of 1 to 3 as $M$ and $A$ are each varied over two orders of magnitude.

In summary, these parametric sweeps reveal three quantities that are highly conserved: (i) the selected Strouhal number $\textit{St}^* = 0.2$ holds for sufficiently high $\textit{Re}_{\!f}\gt \textit{Re}_{\!f}^*$ and across all $M$ and $A$ ; (ii) the rescaled dynamical time scale $A\tau '^*/M = 2\pm 1$ similarly holds; and (iii) the critical Reynolds number $\textit{Re}_{\!f}^*=25$ holds for all $M$ and $A$ . The significance of these results is that they support the view that some key dynamical features are common to flapping-based propulsion across wide ranges of physical conditions, parameter values, and scales. Such universality has been recognised previously with respect to the Strouhal number and to some degree the critical Reynolds number, and future studies may now investigate the newly identified time scale.

5. Sensitivity analysis of model parameters

We next conduct a sensitivity analysis of the constant parameters appearing in the aerodynamic model in order to determine their influence on the dynamical outputs. Specifically, the full slate of parameters include the 9 constants $C_{\!L}^I$ , $C_{\!D}^{0,\textit{FT}}$ , $C_{\!D}^{0,\textit{SF}}$ , $C_{\!D}^I$ , $\alpha _S$ , $\delta _S$ , $C_{\!L}^{\pi /4}$ , $C_{\!D}^{\pi /2,\textit{HR}}$ , and $C_{\!D}^{\pi /2,\textit{LR}}$ , the first 4 of which pertain to the attached flow regime, the next 2 to the stall transition and the last 3 to the separated flow regime. The dynamical outputs of interest are $\textit{Re}_{\!f}^*$ , $\tau '^*$ , and $\textit{St}^*$ . For the purpose of this analysis, we consider fixed values $M=A=1$ , which should be viewed as generally representative given that the outputs are highly conserved. Similarly, we consider fixed $\textit{Re}_{\!f}=10^5$ as representative of conditions $\textit{Re}_{\!f} \gg \textit{Re}_{\!f}^*$ where the asymptotic values $\textit{St}^*$ and $\tau '^*$ are attained. We define a dimensionless sensitivity score to each output quantity $Q$ with respect to each input parameter $p$

(5.1) \begin{equation} S_{Q,p} = \left | \frac {\partial Q}{\partial p}\frac {p_0}{Q_0}\right | \approx \left |\frac {\Delta{\kern-1pt}Q/Q_0}{\Delta p / p_0}\right |\!. \end{equation}

Here, $p_0$ is the nominal parameter value and $Q_0 = Q(p_0)$ is the associated nominal value of the output quantity extracted from the numerical solution. The difference approximation to the derivative allows the score to be estimated by running computations for a perturbed parameter value $p = p_0 + \Delta p$ and extracting the change $\Delta{\kern-1pt}Q = Q - Q_0$ in the output quantity. We implement a two-sided approximation by averaging the results of perturbations of size $\pm 10\,\%$ in each parameter, i.e. $\Delta p = \pm 0.1 p_0$ .

Figure 4. Sensitivity of the dynamical outputs to the constant parameters in the aerodynamic model. The relevant quantities include the critical $\textit{Re}_{\!f}^*$ (green), characteristic growth/decay time scale $\tau '^*$ (cyan), and terminal Strouhal number $\textit{St}^*$ (purple). Here, the dimensionless mass and amplitude are fixed with $M=A=1$ , and $\textit{Re}_{\!f}=10^5\gg \textit{Re}_{\!f}^*$ for the results pertaining to $\textit{St}^*$ and $\tau '^*$ . Sensitivity scores less than $1\,\%$ are not displayed and instead marked with dots.

The results of this analysis are displayed as the bar chart of figure 4 showing the sensitivities of $\textit{Re}_{\!f}^*$ (green), $\tau '^*$ (cyan), and $\textit{St}^*$ (purple) with respect to the 9 aerodynamic parameters. Here, scores less than $0.1=10\,\%$ can be considered low and hence the output is robust to changes in the parameter, and scores above $1=100\,\%$ are high and signify a strong dependency. Moderate sensitivity occurs for intermediate values, and very low scores below $0.01=1\,\%$ are not displayed and instead marked with dots. Evidently, the quantities $\textit{Re}_{\!f}^*$ and $\tau '^*$ display sensitivity only to the coefficients $C_{\!L}^{\pi /4}$ , $C_{\!D}^{\pi /2,\textit{HR}}$ , and $C_{\!D}^{\pi /2,\textit{LR}}$ , which are the subset that govern the separated flow regime. This suggests that the aerodynamics associated with attached flow and stall does not play a role in the stationary-to-locomotion transition. In contrast, the Strouhal number $\textit{St}^*$ is broadly sensitive to nearly all parameters at levels that are at least moderate. Thus, it seems that the mechanisms associated with attached flow, stall and separated flow are all important determinants of $\textit{St}^*$ . These observations are explained by the equilibrium and stability analyses provided in the next section.

It should be noted that nearly all 9 parameters affect at least one output quantity at a moderate level of $10\,\%$ or greater. The exception is $\delta _S$ , which defines the angular breadth of the stall transition and displays uniformly low sensitivities. Even this parameter is necessary to include in the model in order to mathematically specify the cross-over between the flow regimes, but apparently the outputs are robust to its numerical value. Taken together, these results indicate that the model is minimal in the sense of containing no unnecessary parameters.

6. Interpretations via equilibrium and stability analyses

The tractability of the quasi-steady formulation allows for analytical insights into key features of the dynamics, and we first seek interpretations for the critical Reynolds number $\textit{Re}_{\!f}^*$ for take-off from rest and the characteristic dimensional time scale $\tau$ to reach the terminal locomotion state. Our theory is based on the linear stability of the stationary state, for which the attack angle is high and we may safely assume the separated flow forms of the lift and drag coefficients. It proves adequate to focus on one stroke in the flapping cycle and assume steady vertical motion $v_y \cong \overline {v}_y = 4af$ which here is chosen to be the mean speed. Perturbation of the flight speed introduces the small parameter $\varepsilon = v_x/v_y \ll 1$ whose dynamics we seek through equation (2.1). The terms involve $v = \sqrt {v^2_x+v^2_y} = v_y\sqrt {\varepsilon ^2+1} \cong v_y$ and $\textit{Re} = \rho \ell v/\mu \cong \rho \ell v_y/\mu \cong 4 \rho \ell a f/\mu = 4 \textit{Re}_{\!f}$ , where the symbol ` $\cong$ ’ is used wherever the steady motion approximation or linearisation in $\varepsilon$ is invoked. Inputting the attack angle $\alpha \cong \pi /2-\varepsilon$ yields the coefficients of lift $C^S_L(\textit{Re},\alpha ) = C^{\pi /4}_{L}\sin 2\alpha \cong C^{\pi /4}_{L} \sin (\pi -2\varepsilon ) \cong 2 C^{\pi /4}_{L} \varepsilon$ and drag $C^S_D(\textit{Re},\alpha ) = (C^{\pi /2,\textit{HR}}_D + C^{\pi /2,\textit{LR}}_D / \textit{Re}) \sin ^2 \alpha \cong (C^{\pi /2,\textit{HR}}_D + C^{\pi /2,\textit{LR}}_D / 4\textit{Re}_{\!f}) \sin ^2 (\pi /2-\varepsilon ) \cong C^{\pi /2,\textit{HR}}_D + C^{\pi /2,\textit{LR}}_D / 4\textit{Re}_{\!f}$ . Inserting all these forms into equation (2.1) and simplifying yields an expression for the linearised dynamics:

(6.1) \begin{align} \begin{aligned} \dot {\varepsilon } = \tau ^{-1} \varepsilon \,\Rightarrow \, \varepsilon (t) = \varepsilon _0 e^{t/\tau } \textrm {,}\, \tau ^{-1} = \frac {1}{2} \frac {\rho \ell v_y}{m+m_{11}}\left [2 C^{\pi /4}_{L} - \left (C^{\pi /2,\textit{HR}}_D + \frac {C^{\pi /2,\textit{LR}}_D}{4\textit{Re}_{\!f}}\right )\right ]\!. \end{aligned} \end{align}

Hence, the perturbation is predicted to decrease or increase according to the sign of $\tau$ , which is the associated exponential decay or growth time scale. Notably, $\tau \lt 0$ for sufficiently small $\textit{Re}_{\!f}$ , indicating that the stationary state is stable, whereas $\tau \gt 0$ for larger $\textit{Re}_{\!f}$ and the wing takes off into forward flight. The critical value $\textit{Re}^*_{\!f}$ marking the transition is determined by $\tau ^{-1}(\textit{Re}^*_{\!f})=0$ :

(6.2) \begin{align} \begin{aligned} \textit{Re}^*_{\!f} = \frac {C^{\pi /2,\textit{LR}}_D}{4(2 C^{\pi /4}_{L} - C^{\pi /2,\textit{HR}}_D)} = 25, \end{aligned} \end{align}

where the computed result reflects the nominal set of model parameters. This matches the critical value identified in the numerical results presented in §§ 3 and 4.

These results interpret the take-off condition for forward flight as requiring sufficiently high Reynolds number such that the horizontal component of lift exceeds drag at high attack angles, thereby generating a forward directed force vector that destabilises the stationary state. The expression (6.2) also explains why the critical Reynolds number is almost entirely independent of $M$ and $A$ per the results of figure 3(c), and why it is sensitive to the three aerodynamic coefficients $C^{\pi /4}_{L}$ , $C^{\pi /2,\textit{HR}}_D$ , and $C^{\pi /2,\textit{LR}}_D$ of the separated flow regime and insensitive to the other model parameters per the results of figure 4. Similarly, this analysis interprets the time scale $\tau$ via equation (6.1) as arising from the inertial resistance to motion in response to the aerodynamic forces. The dimensionless form and its limit for $\textit{Re}_{\!f} \gg \textit{Re}^*_{\!f}$ are

(6.3) \begin{align} \begin{aligned} (\tau ')^{-1} = 2\frac {A}{M}\left [2 C^{\pi /4}_{L} - \left (C^{\pi /2,\textit{HR}}_D + \frac {C^{\pi /2,\textit{LR}}_D}{4\textit{Re}_{\!f}}\right )\right ] \rightarrow 2\frac {A}{M}\left (2 C^{\pi /4}_{L} - C^{\pi /2,\textit{HR}}_D\right ) = \frac {A}{2.5M}. \end{aligned} \end{align}

This explains why the rescaled quantity $A\tau '^*/M$ is nearly conserved as shown in figure 3(d), and the predicted value of 2.5 corresponds well with the range of 1 to 3 seen in the numerical results. The above expression also explains the sensitivity of $\tau '^*$ to the separated flow parameters and insensitivity to all others per figure 4.

A similar attempt to account for the selected Strouhal number views the forward flight state as arising from a balance of horizontal forces in which the forward-inclined lift vector balances the rearward-pointing drag. We again assume constant $v_y \cong 4af$ , but $v_x$ cannot be assumed to be small. Equilibrium requires that $C_{\!L}(\alpha ,\textit{Re})\sin \alpha =C_{\!D}(\alpha ,\textit{Re})\cos \alpha$ and hence $C_{\!L}(\alpha ,\textit{Re})/C_{\!D}(\alpha ,\textit{Re})=\cot \alpha$ . This condition represents an implicit relation for the unknown $v_x$ , which appears in the attack angle $\tan \alpha = v_y/v_x \cong 4af/v_x = 2\textit{St}$ and Reynolds number $\textit{Re} = \rho \ell \sqrt {v^2_x + v^2_y} / \mu$ . The equilibrium solutions are those $\alpha$ marking the intersections shown in figure 1( f). Further analytical progress is hindered because the dynamically selected values $\alpha = 18^\circ$ to $27^\circ$ fall in or near the stall range, necessitating consideration of the full forms of $C_{\!L}$ and $C_{\!D}$ including the stall parameters and those of the attached and separated flow regimes. Nonetheless, the relatively narrow range of $\alpha$ across all $\textit{Re} \gt \textit{Re}_{\!f}^*$ implies a similarly conserved Strouhal number $\textit{St}^*\cong (\tan \alpha )/2 = 0.16$ to $0.21$ , which aligns with the findings of figures 2(e) and 3(a). Further, that the associated attack angles span the flow regimes explains why $\textit{St}^*$ is sensitive to the full array of model parameters per figure 4. We also note that these complexities similarly hinder any analytical accounting for the value $\textit{Re}_{\!f}=18$ marking the lower boundary of the hysteresis window and which would presumably be determined by loss of equilibrium for the locomoting state.