2. Problem description and methodology

A two-dimensional tear-shaped rigid foil with a semicircular leading edge of diameter $D$ and chord length $L$ ($D/L = 0.2$) is positioned in a quiescent fluid. The foil is free to move horizontally from right to left, and the foil's swimming direction is the negative $x$-axis (figure 1a). In the self-propelled pitching mode, the foil rotates around the centre of its semicircular leading edge, exhibiting a time-dependent pitching angle $\theta (t)$ that varies during the oscillation (figure 1b). In the self-propelled heaving mode, the foil oscillates transversely to the swimming direction ($y$), where its vertical movement ($y$-direction) is controlled by $h(t)$ (figure 1c).

Figure 1. (a) Sketch of the computational domain; geometry and kinematics of the (b) pitching foil and (c) heaving foil; definition of (d) basic motion (BM), (e) perturbation motion (PM) and (f) accumulated motion (AM). Here, $A_{\_max}$ denotes the maximal tailbeat amplitude.

We have established three modes: the basic mode (BM) characterised by a sinusoidal wave (figure 1d); the perturbation mode (PM) composed of high-frequency, low-amplitude oscillations (figure 1e); and the accumulated mode (AM), which combines the sinusoidal wave (BM) with perturbations (PM) to describe the rhythmic tailbeats of the foil (figure 1f). The definitions of these modes are as follows:

(2.1) \begin{equation} \left.\begin{array}{ll@{}} A_{b}(t) = A_{b\_{max}} \sin(2 {\rm \pi} f_{b} t) & \text{basic mode, BM},\\ A_{p}(t) = A_{p\_{max}} \sin(2 {\rm \pi} f_{p} t) & \text{perturbation mode, PM},\\ A_{a}(t) = A_{b}(t) + A_{p}(t) & \text{accumulated mode, AM}, \end{array}\right\} \end{equation}

where $A_b(t)$ and $A_p(t)$ denote the time-dependent taibeat amplitude in BM and PM, $A_{b\_{max}}$ and $A_{p\_{max}}$ describe the maximal taibeat amplitude in BM and PM, $f_{b}$ and $f_{p}$ refer to the oscillating frequency in BM and PM, respectively, and $t$ is the time. In the rest of the paper, subscripts $b$ and $a$ refer to BM and AM, respectively. We describe the foil's motions as follows:

(2.2) \begin{equation} \left.\begin{array}{ll@{}} A_{b}(t) = \theta_b(t) = \theta_{b\_{max}} \sin(2 {\rm \pi} f_{b} t) & \text{BM, Pitching mode}, \\ A_{p}(t) = \theta_p(t) = \theta_{p\_{max}} \sin(2 {\rm \pi} f_{p} t) & \text{PM, Pitching mode}, \\ A_{b}(t) = h_b(t) = h_{b\_{max}} \sin(2 {\rm \pi} f_{b} t) & \text{BM, Heaving mode}, \\ A_{p}(t) = h_p(t) = h_{p\_{max}} \sin(2 {\rm \pi} f_{p} t) & \text{PM, Heaving mode}, \end{array}\right\} \end{equation}

where $\theta _{b\_{max}}$ and $\theta _{p\_{max}}$ refer to the maximal pitching angle in BM and PM, respectively. Since the foil thickness is $D/L = 0.2$, we have $\sin (\theta _{b\_{max}}) = A_{b\_{max}} / (L - D/2) = A_{b\_{max}} / 0.9L$ and $\sin (\theta _{p\_{max}}) = A_{p\_{max}} / (L - D/2) = A_{p\_{max}} / 0.9L$ for the pitching mode, respectively. However, $h_{b\_{max}} = A_{b\_{max}}$ and $h_{p\_{max}} = A_{p\_{max}}$ are the maximal heaving distance for the heaving mode in BM and PM, respectively. Without loss of generality, both the self-propelled pitching and heaving modes employ $A_{b\_{max}} = 0.15L$ and $f_b = 1$ Hz in our simulations.

We have introduced non-dimensionalised flapping frequency and amplitude, denoted as

(2.3a,b)\begin{equation} \tilde{f} = f_{p}/f_{b},\quad \tilde{A} = A_{p\_{max}}/A_{b\_{max}}, \end{equation}

to control the superimposed perturbations. In the present study, we considered $\tilde {f} = 4- 10$ with an interval of $\Delta \tilde {f} = 1$ and $\tilde {A} = 0.01$, 0.02, 0.04, 0.05, 0.06, 0.08 and 0.10. We further defined the swimming number (Gazzola et al. Reference Gazzola, Argentina and Mahadevan2014) ratio as

(2.4)\begin{equation} \widetilde{Sw} = \frac{Sw_p}{Sw_b} = \frac{(2{\rm \pi} L f_p A_{p\_{max}})/\nu}{(2{\rm \pi} L f_b A_{b\_{max}})/\nu} = \tilde{f}\tilde{A}, \end{equation}

where $\nu$ is the fluid kinematic viscosity. To describe the swimming efficiency of the self-propelled foil, we employed the energy consumption coefficient $C_E$, defined as $C_E = CoT/\rho \nu ^2$, where $CoT = \bar {P} /\bar {u}$ is the mechanical cost of transport (Bale et al. Reference Bale, Hao, Bhalla and Patankar2014; Zhang, Zhang & Huang Reference Zhang, Zhang and Huang2022), given by the ratio of time-averaged total power ($\bar {P}$) to time-averaged steady swimming speed ($\bar {u}$), and $\rho$ is the fluid (foil) density.

All swimming speeds observed in our simulations are normalised by the foil length $L$. The non-dimensionalised time-averaged swimming speed is defined as

(2.5)\begin{equation} \tilde{u} = \bar{u}_a/\bar{u}_b, \end{equation}

to evaluate the effect of perturbations on the swimming speed generated by self-propelled foils. Furthermore, we defined

(2.6a,b)\begin{equation} \tilde{P} = \bar{P}_a/\bar{P}_b,\quad \tilde{C}_E = C_{E\_a}/C_{E\_b}, \end{equation}

to assess whether AM incurs greater energy cost and enhances swimming efficiency, respectively. In (2.5) and (2.6), $\bar {u}_a$ ($\bar {u}_b$), $\bar {P}_a$ ($\bar {P}_b$) and $C_{E\_{a}}$ ($C_{E\_{b}}$) denote the time-averaged swimming speed, time-averaged power cost and swimming efficiency generated by the self-propelled flapping foil undergoing the AM (BM) when the swimming is steady, respectively. A larger $\tilde {u}$ and $\tilde {P}$ indicate AM enjoys a faster swimming speed and greater energy consumption than BM, respectively. For AM, $\tilde {C}_E < 1$ and $\tilde {C}_E > 1$ correspond to improvements and decreases in swimming efficiency compared with the non-perturbed conditions (BM), respectively.

The unsteady flow around the self-propelled foil was simulated using the immersed boundary (IB) method (Peskin Reference Peskin2002). The IB method (also referred to as ‘Cartesian grid methods’) employs a non-body conformal Cartesian grid, where the Cartesian grids are generated regardless of the body surface. As a result, the body surface cuts through this Cartesian grid. Because the grid does not conform to the body surface, imposing the boundary conditions will require modifying the governing equations in the vicinity of the solid boundary. With the precise modifications, the governing equations can be discretized without resorting to coordinate transformation or complex discretization operators. More details about the IB method and relevant mathematical formulation can be found from Mittal & Iaccarino (Reference Mittal and Iaccarino2005) and Griffith & Patankar (Reference Griffith and Patankar2020). We employed the open-source CFD software IBAMR to conduct the numerical investigations (IBAMR 2014). This software has been widely employed to gain insights into the hydrodynamics of fish-like swimming (Bhalla et al. Reference Bhalla, Bale, Griffith and Patankar2013a; Bhalla, Griffith & Patankar Reference Bhalla, Griffith and Patankar2013b; Tytell et al. Reference Tytell, Leftwich, Hsu, Griffith, Cohen, Smits, Hamlet and Fauci2016; Hoover et al. Reference Hoover, Cortez, Tytell and Fauci2018; Zhang et al. Reference Zhang, Pan, Chao and Zhang2018; Yang & Wu Reference Yang and Wu2022; Chao, Bhalla & Li Reference Chao, Bhalla and Li2023). In our research, we adopted a computational domain in the form of a rectangular box measuring $80 \times 20 L$, with periodic boundary conditions applied along the axial direction and no-slip boundary conditions imposed in the lateral direction (figure 1a). To ensure statistically steady swimming, all cases were simulated for twenty normalised swimming periods $t^{*} = tf_b$.

The grid and time step size convergence study is conducted for both the self-propelled pitching and heaving mode with the most strenuous parameter values of $(\,\tilde {f}, \tilde {A}) = (10, 0.10)$. For the grid convergence study, three different grids with uniform mesh spacings of $\Delta x/L = \Delta y/L = 0.008$ (M1), 0.005 (M2) and 0.002 (M3) are considered. The time step size is fixed at $\Delta t \times f_p = 0.001$. Figures 2(a) and 2(c) present the instantaneous swimming velocity $u_{pitching}$ (self-propelled pitching mode) and $u_{heaving}$ (self-propelled heaving mode) for the three grids, respectively. Table 1 lists the time-averaged steady swimming speeds derived from the self-propelled pitching ($\bar {u}_{pitching}$) and heaving ($\bar {u}_{heaving}$) mode. There is a relatively small difference between M2 and M3, i.e. $\Delta \bar {u}_{pitching} = 0.712\,\%$ and $\Delta \bar {u}_{heaving} = 1.270\,\%$. For the remainder of the simulations, mesh M2 is selected. Using M2, the time step size convergence study is performed by selecting three values of $\Delta t \times f_p$: $\Delta t \times f_p = 0.01$ ($\Delta t1$), 0.001 ($\Delta t2$) and 0.0005 ($\Delta t3$). The difference between the time-dependent swimming speed using $\Delta t2$ and $\Delta t3$ is negligible (figure 2b,d), where $\Delta \bar {u}_{pitching}$ and $\Delta \bar {u}_{heaving}$ is less than $1.0\,\%$ (table 2). Grid M2 and time step size $\Delta t2$ are selected based on accuracy and computational resources for the remainder of the simulations. More convergence studies related to the numerical method can be found in our previous study (Chao et al. Reference Chao, Bhalla and Li2023).

Figure 2. (a) Grid and (b) time step size convergence study of self-propelled pitching motion; (c) grid and (d) time step size convergence study of self-propelled heaving motion, respectively. $t^{*} = tf_b$. M1, $\Delta x/L = \Delta y/L = 0.008$; M2, $\Delta x/L = \Delta y/L = 0.005$; M3, $\Delta x/L = \Delta y/L = 0.002$. $\Delta t1$, $\Delta t \times f_p = 0.01$; $\Delta t2$, $\Delta t \times f_p = 0.001$; $\Delta t3$, $\Delta t \times f_p = 0.0005$.

Table 1. Grid convergence study with $\Delta t \times f_p = 0.001$.

Table 2. Time-step size convergence study with $\Delta x/L = \Delta y/L = 0.005$.

4. Scaling laws

4.1. Scaling for swimming speeds

To better understand how perturbations ($\,\tilde {f}$ and $\tilde {A}$) affect the non-dimensionalised swimming speed ($\tilde {u}$) of a self-propelled pitching/heaving foil, we derived scaling laws by balancing time-averaged inertial forces $\bar {F}_I = ({1}/{nt^*}) \int _{\alpha }^{\alpha +nt^*} F_I \,\textrm {d}t^*$ and viscous drag $\bar {F}_D = ({1}/{nt^*}) \int _{\alpha }^{\alpha +nt^*} F_D \,\textrm {d}t^*$ when the time-averaged swimming speed is steady, where $\alpha$ denotes arbitrary normalised time. For the per unit depth scales, we have $F_I \sim \rho V_{lateral}^2 L$ and $F_D \sim \mu ({u}/{\delta })L$, where $\mu = \rho \nu$ is the fluid viscosity, $\delta \sim L/\sqrt {Re}$ is the boundary layer thickness (Gazzola et al. Reference Gazzola, Argentina and Mahadevan2014), $Re = u L/ \nu$ is the Reynolds number, and $V_{lateral}$ denotes the foil's lateral velocity and can be derived from the foil's kinematic equations:

(4.1) \begin{align} \left.\begin{array}{ll@{}} V_{lateral} = 2{\rm \pi} f_b A_{b\_max} \cos(2 {\rm \pi} f_b t) \times \cos[\theta_b(t)] & \text{BM, Pitching mode},\\ V_{lateral} = 2{\rm \pi} [f_b A_{b\_max}\cos(2 {\rm \pi} f_b t) + f_p A_{p\_max}\cos(2 {\rm \pi} f_p t) ] & \\ \qquad\qquad \times \cos[\theta_a(t)] & \text{AM, Pitching mode},\\ V_{lateral} = 2{\rm \pi} f_b A_{b\_max}\cos(2 {\rm \pi} f_b t) & \text{BM, Heaving mode}, \\ V_{lateral} = 2{\rm \pi} [f_b A_{b\_max}\cos(2 {\rm \pi} f_b t) + f_p A_{p\_max}\cos(2 {\rm \pi} f_p t) ] & \text{AM, Heaving mode}. \end{array}\right\} \end{align}

Considering a small $\theta _{b\_max}$ and even smaller $\theta _{p\_max}$ in pitching mode, we get $\cos [\theta _b(t)] \approx \cos [\theta _a(t)] \approx 1$. Equation (4.1) further suggests

(4.2) \begin{equation} \left.\begin{array}{ll@{}} \displaystyle\dfrac{1}{nt^*} \int_{\alpha}^{\alpha+nt^*} V_{lateral}^2 \,\text{d}t^* \sim f_b^2 A_{b\_max}^2 & \text{BM,}\\ \displaystyle\dfrac{1}{nt^*} \int_{\alpha}^{\alpha+nt^*} V_{lateral}^2 \,\text{d}t^* \sim f_b^2 A_{b\_max}^2 + f_p^2 A_{p\_max}^2 & \text{AM.} \end{array}\right\} \end{equation}

By balancing $\bar {F}_I$ and $\bar {F}_D$, we obtain

(4.3) \begin{equation} \left.\begin{array}{ll@{}} \rho (f_b A_{b\_{max}})^2 L \sim \rho \nu \bar{u}_b \left(\dfrac{\bar{u}_b L}{\nu} \right)^{1/2}\nu & \text{BM,}\\ \rho [(f_b A_{b\_{max}})^2 + (f_p A_{p\_{max}})^2 ] L \sim \rho \nu \bar{u}_a\left(\dfrac{\bar{u}_a L}{\nu} \right)^{1/2}\nu & \text{AM,} \end{array}\right\} \end{equation}

resulting in

(4.4) \begin{equation} \left.\begin{array}{ll@{}} \bar{u}_b \sim (f_b A_{b\_{max}})^{4/3} L^{1/3} \nu^{{-}1/3}\sim \dfrac{(f_b A_{b\_{max}})^{4/3} L^{4/3}}{\nu^{4/3}} L^{{-}1} \nu & \text{BM,}\\ \bar{u}_a \sim [(f_b A_{b\_{max}})^2 + (f_p A_{p\_{max}})^2]^{2/3} L^{1/3} \nu^{{-}1/3} \\ \quad \ \sim \left[\dfrac{ (f_b A_{b\_{max}})^{2} L^{2}}{\nu^{2}} + \dfrac{ (f_p A_{p\_{max}})^{2} L^{2}}{\nu^{2}}\right]^{2/3} L^{{-}1} \nu & \text{AM.} \end{array}\right\} \end{equation}

Recalling $Sw_p = {2{\rm \pi} L f_p A_{p\_{max}}}/{\nu }$ and $Sw_b = {2{\rm \pi} L f_b A_{b\_{max}}}/{\nu }$ (2.4), we naturally have

(4.5) \begin{equation} \left.\begin{array}{ll@{}} \bar{u}_b \sim Sw_{b}^{4/3} L^{{-}1} \nu & \text{BM}, \\ \bar{u}_a \sim (Sw_{b}^{2} + Sw_{p}^{2})^{2/3} L^{{-}1} \nu & \text{AM}. \end{array}\right\} \end{equation}

Equation (4.5) further reveals $Re_{b} \sim Sw_{b}^{4/3}$ for BM, conforming with the works of Gazzola et al. (Reference Gazzola, Argentina and Mahadevan2014), while $Re_{a} \sim (Sw_{b}^{2} + Sw_{p}^{2})^{2/3}$ for AM in the viscous flow.

Based on (2.5) and (4.5), we naturally obtain

(4.6)\begin{equation} \tilde{u} \sim (1 + \widetilde{Sw}^{2})^{2/3}, \end{equation}

where $\tilde {u}$ is always larger than 1 (4.6, figures 3a and 4a), suggesting that perturbations can constantly improve the time-averaged swimming speeds of the self-propelled foil.

4.2. Scaling for power costs

As mentioned in § 3.2, the phase lag between the foil motion and the resulting fluid dynamics response plays an important role in affecting the power cost of a self-propelled flapping foil. Therefore, we first introduce the definitions of the phase lag in the self-propelled pitching and heaving mode before the derivation of the power cost scaling. Figure 5(a) reveals that the flow-induced torque $M_f(t)$ provided by BM remains at the same frequency with foil pitching acceleration $\ddot {\theta }(t)$, while the perturbation dominates the $M_f(t)$ profile in AM. Particularly, a phase lag $\phi$ exists between $M_f(t)$ and $\ddot {\theta }(t)$. Similarly, we can also observe a phase lag between flow-induced lift $L_f(t)$ and foil heaving acceleration $\ddot {h}(t)$ in the self-propelled heaving mode (figure 5d). As the state of a vibrating curve can be modelled as a point in a Hilbert space, we use the Hilbert transform to calculate $\phi _{b}$ and $\phi _{a}$ (Williamson & Govardhan Reference Williamson and Govardhan2004) referred to as the BM-based and AM-based phase lag, respectively. Figure 5(b,e) reveals that $\sin (\phi _{a})/\sin (\phi _{b})$ decreases with an increase of $\tilde {f}$ and/or $\tilde {A}$. This tendency is attributed to the high-frequency components introduced by superimposed perturbations, which induce rapid transitions in flapping velocity and thus accelerate the response of the fluid dynamics to the flapping motion. Particularly, the minimal $\sin (\phi _{a})/\sin (\phi _{b})$ values at the specific $\tilde {f}$ agree well with the data-driven $\tilde {P} = 1$ line in self-propelled heaving mode (figure 5e), suggesting this phase lag is a critical parameter in affecting the power cost. Naturally, $\sin (\phi _{a})/\sin (\phi _{b})$ is a function of $\tilde {f}$ and $\tilde {A}$, and these functions corresponding to self-propelled pitching and heaving mode are fitted as

(4.7) \begin{equation} \left.\begin{array}{ll@{}} \dfrac{\sin(\phi_{a})}{\sin(\phi_{b})} \sim {\rm e}^{{-}2\tilde{f}^2\tilde{A}/3} & \text{Pitching mode}, \\ \dfrac{\sin(\phi_{a})}{\sin(\phi_{b})} \sim {\rm e}^{-\tilde{f}^2\tilde{A}} & \text{Heaving mode}. \end{array}\right\} \end{equation}

Figure 5(c,f) demonstrates a good match between the numerical results and the fitted equation. However, the fitted equations also indicate that the relationship between phase difference and superimposed perturbation is highly complex.

Figure 5. (a) Normalised foil pitching acceleration $\ddot {\theta }(t)$ and flow-induced torque $M_f(t)$ in BM and AM, where $(\,\tilde {f}, \tilde {A}) = (5, 0.05)$. (b) The $\sin (\phi _{a})/\sin (\phi _{b})$ contour in self-propelled pitching mode. (c) Fitting for $\sin (\phi _{a})/\sin (\phi _{b})$ in self-propelled pitching mode, where $\sin (\phi _{a})/\sin (\phi _{b}) \sim \textrm {e}^{-2\tilde {f}^2\tilde {A}/3}$. (d) Normalised foil heaving acceleration $\ddot {h}(t)$ and flow-induced lift $L_f(t)$ in BM and AM, where $(\,\tilde {f}, \tilde {A}) = (5, 0.05)$. (e) The $\sin (\phi _{a})/\sin (\phi _{b})$ contour in self-propelled heaving mode, where the dashed line refers to the $\tilde {P} = 1$ line in figure 4(b) and the green symbol denotes the location of the minimal $\sin (\phi _{a})/\sin (\phi _{b})$ at the specific $\tilde {f}$. (f) Fitting for $\sin (\phi _{a})/\sin (\phi _{b})$ in self-propelled heaving mode, where $\sin (\phi _{a})/\sin (\phi _{b}) \sim \textrm {e}^{-\tilde {f}^2\tilde {A}}$.

We now go back to the power cost scaling. For the pitching foil, the time-averaged power cost is considered as the product of the net torque $M(t^*) = M_{i}(t^*) - M_{f}(t^*)$ and $\dot {\theta }(t^*)$, i.e.

(4.8)\begin{equation} \bar{P} = \frac{1}{nt^*} \int_{\alpha}^{\alpha+nt^*} M(t^*) \dot{\theta}(t^*) \,\text{d}t^* = \frac{1}{nt^*} \int_{\alpha}^{\alpha+nt^*} [M_i(t^*) - M_f(t^*)]\dot{\theta}(t^*) \,\text{d}t^*, \end{equation}

where inertial torque $M_i(t^*) = J\ddot {\theta }(t^*)$ with the mass moment of inertia of the foil $J$ and flow-induced torque $M_{f}(t^*) = - c_{1}\ddot {\epsilon }(t^*) - c_{2}u \dot {\epsilon } (t^*)$ with two positive coefficients $c_1 \sim L^{4}$ and $c_2 \sim L^{3}$. Here, $\ddot {\epsilon }(t^*) = \ddot {\theta }(t^* + \phi )$ and $u\dot {\epsilon }(t^*) = u\dot {\theta }(t^* + \phi )$ are sourced from the Euler angular acceleration and Coriolis acceleration of the fluid around the foil (Alam & Muhammad Reference Alam and Muhammad2020), respectively. This phase lag $\phi$, driven by diffusion-convection of the vorticity layer generated on the foil surface, corresponding to the flow retardation effect (Granger & Païdoussis Reference Granger and Païdoussis1996; Chao et al. Reference Chao, Jia, Wang, Liberzon, Ravi, Couzin and Li2024), is crucial in transferring power between the body and fluid (Han et al. Reference Han, Huang, Qin, Wang and Zhao2023).

Note $\int _{\alpha }^{\alpha +nt^*} M_i(t^*)\dot {\theta }(t^*) \,\textrm {d}t^* = \int _{\alpha }^{\alpha +nt^*} J\ddot {\theta }(t^*)\dot {\theta }(t^*) \,\textrm {d}t^* = 0$, and the power cost is actually the integral of the product of $M_f (t^*)$ and $\dot {\theta }(t^*)$. We first calculate the flow-induced torque, as

(4.9)\begin{align} \left.\begin{array}{ll@{}} M_{b}(t^*) ={-} (4{\rm \pi} ^2) J f_{b}^2 \theta_{b\_max} \sin(2{\rm \pi} f_{b}t^*) \\ \qquad\qquad -\, (4{\rm \pi} ^2)c_1f_{b}^2 \theta_{b\_max} \sin(2{\rm \pi} f_{b}t^* + \phi_{b}) \\ \qquad\qquad +\, (2{\rm \pi} )c_{2}u_{b} f_{b}\theta_{b\_max}\cos(2{\rm \pi} f_{b}t^* + \phi_{b}) & \text{BM,} \\ M_{a}(t^*) ={-} (4{\rm \pi} ^2) J[f_{b}^2 \theta_{b\_max} \sin(2{\rm \pi} f_{b}t^*) + f_{p}^2 \theta_{b\_max} \sin(2{\rm \pi} f_{p}t^*)]\\ \qquad\qquad -\, (4{\rm \pi} ^2)c_1[f_{b}^2 \theta_{b\_max} \sin(2{\rm \pi} f_{b}t^* + \phi_{a}) \!+ f_{p}^2 \theta_{p\_max} \sin(2{\rm \pi} f_{p}t^* + \phi_{a}) ]\\ \qquad\qquad +\, (2{\rm \pi} )c_{2}u_{a}[f_{b} \theta_{b\_max} \cos(2{\rm \pi} f_{b}t^* + \phi_{a}) \!+ f_{p} \theta_{p\_max} \cos(2{\rm \pi} f_{p}t^* \!+ \phi_{a})], & \text{AM,} \end{array}\right\} \end{align}

where $\phi _b$ and $\phi _a$ refer to the phase lag in BM and AM, respectively. The $\dot {\theta }(t)$ can be directly obtained from the kinematics, as

(4.10)\begin{equation} \left.\begin{array}{ll@{}}\dot{\theta}_{b}(t) = 2{\rm \pi} f_{b} \theta_{b\_max} \cos(2{\rm \pi} f_{b}t) & \text{BM},\\ \dot{\theta}_{a}(t) = 2{\rm \pi} [f_{b} \theta_{b\_max} \cos(2{\rm \pi} f_{b}t) + f_{p} \theta_{p\_max}\cos(2{\rm \pi} f_{p}t )] & \text{AM}. \end{array}\right\} \end{equation}

Therefore, time-averaged power costs read as

(4.11)\begin{align} \left.\begin{array}{ll@{}}\bar{P}_{b}(t) \sim{-}c_3 L f_{b}^3 \theta_{b\_max}^2\sin(\phi_{b}) + \bar{u}_b f_{b}^2\theta_{b\_max}^2\cos(\phi_{b})\\ \qquad \ \sim \left[{-}c_3 \dfrac{ L^3 f_{b}^3 \theta_{b\_max}^3 }{\nu^{3}}\theta_{b\_max} ^{{-}1} \nu \sin(\phi_{b}) + \bar{u}_b \dfrac{ L^2 f_{b}^2 \theta_{b\_max}^2 }{\nu^{2}}\cos(\phi_{b})\right]L^{{-}2}\nu^{2}\\ \qquad \ \sim [{-}c_3 Sw_{b}^3 \theta_{b\_max}^{{-}1} \nu \sin(\phi_{b}) + \bar{u}_b Sw_{b}^2 \cos(\phi_{b})]L^{{-}2}\nu^{2} &\text{BM},\\ \bar{P}_{a}(t) \sim{-}c_3 L (f_{b}^3\theta_{b\_max}^2 \!+ f_{p}^3\theta_{p\_max}^2)\sin(\phi_{a}) \!+ \bar{u}_a (f_{b}^2\theta_{b\_max}^2 \!+ f_{p}^2\theta_{p\_max}^2)\cos(\phi_{a})\\ \qquad \ \sim \left[{-}c_3 \left( \dfrac{ L^3 f_{b}^3 \theta_{b\_max}^3 }{\nu^{3}}\theta_{b\_max}^{{-}1} + \dfrac{ L^3 f_{p}^3 \theta_{p\_max}^3 }{\nu^{3}}\theta_{p\_max}^{{-}1}\right) \nu \sin(\phi_{a}) \right]L^{{-}2}\nu^{2}\\ \qquad \quad +\, \bar{u}_a \left[ \left( \dfrac{ L^2 f_{b}^2 \theta_{b\_max}^2 }{\nu^{2}} + \dfrac{ L^2 f_{p}^2 \theta_{p\_max}^2 }{\nu^{2}}\right) \cos(\phi_{a}) \right]L^{{-}2}\nu^{2}\\ \qquad \ \sim [{-}c_3 ( Sw_{b}^3 \theta_{b\_max}^{{-}1} + Sw_{p}^3 \theta_{p\_max} ^{{-}1} ) \nu \sin(\phi_{a})]L^{{-}2}\nu^{2}\\ \qquad \quad + \,\bar{u}_a [( Sw_{b}^2 + Sw_{p}^2)\cos(\phi_{a})]L^{{-}2}\nu^{2} & \text{AM}, \end{array}\right\} \end{align}

with the positive coefficient $c_3$. Then, we have the $\tilde {P}$ scaling for the pitching foil as

(4.12)\begin{align} \tilde{P} & \sim \frac{-c_3 ( Sw_{b}^3 \theta_{b\_max}^{{-}1} + Sw_{p}^3 \theta_{p\_max}^{{-}1} ) \nu \sin(\phi_{a}) + \bar{u}_a ( Sw_{b}^2 + Sw_{p}^2 )\cos(\phi_{a})}{-c_3 Sw_{b}^3 \theta_{b\_max}^{{-}1} \nu \sin(\phi_{b}) + \bar{u}_b Sw_{b}^2 \cos(\phi_{b})}\nonumber\\ & \sim \frac{-c_3 (Sw_{b}\theta_{b\_max}^{{-}1} + \widetilde{Sw}^2 Sw_{p} \theta_{p\_max}^{{-}1} ) \nu \sin(\phi_{a}) + \bar{u}_a (1 + \widetilde{Sw}^2 )\cos(\phi_{a})}{-c_3 Sw_{b} \theta_{b\_max}^{{-}1} \nu \sin(\phi_{b}) + \bar{u}_b \cos(\phi_{b})}\nonumber\\ & \sim \frac{-c_3 Sw_{b} (\tilde{\theta} + \widetilde{Sw}^3 ) \nu \sin(\phi_{a}) + \theta_{p\_max} \bar{u}_a (1 + \widetilde{Sw}^2 )\cos(\phi_{a})}{-c_3 Sw_{b} \tilde{\theta} \nu \sin(\phi_{b}) + \theta_{p\_max} \bar{u}_b \cos(\phi_{b})}\nonumber\\ & \sim \frac{-c_3 Sw_{b} (\tilde{\theta} + \widetilde{Sw}^3 ) \nu \sin(\phi_{a}) + c_4 \theta_{p\_max} \bar{u}_b (1 + \widetilde{Sw}^2 )^{5/3}\cos(\phi_{a})}{-c_3 Sw_{b} \tilde{\theta} \nu \sin(\phi_{b}) + \theta_{p\_max} \bar{u}_b \cos(\phi_{b})}, \end{align}

where $c_4$ is a positive coefficient. For the low-amplitude pitching, we have $\tilde {\theta } = \theta _{p\_max} / \theta _{b\_max} = A_{p\_max}/ A_{b\_max} = \tilde {A}$ and $\theta _{p\_max} \approx 0$, and therefore

(4.13)\begin{equation} \tilde{P} \sim \frac{-c_3 Sw_{b} (\tilde{\theta} + \widetilde{Sw}^3 ) \nu \sin(\phi_{a})}{-c_3 Sw_{b} \tilde{\theta} \nu \sin(\phi_{b}) } \sim (1 + \tilde{A}^{{-}1}\widetilde{Sw}^3) \cdot \frac{\sin(\phi_{a})}{\sin(\phi_{b})}. \end{equation}

Similar to the pitching motion, we have the time-averaged power cost of the heaving foil as

(4.14)\begin{equation} \bar{P} = \frac{1}{nt^*} \int_{\alpha}^{\alpha+nt^*} [L_i(t^*) - L_f(t^*)]\dot{h}(t^*) \,\text{d}t^*, \end{equation}

where $L_i(t^*)$ and $L_{f}(t^*)$ refer to the inertial and flow-induced lift, respectively, and $L_i(t^*) \sim \ddot {h}(t^*)$. Different from the condition in self-propelled pitching mode, we have $L_{f}(t^*) = -c_{5}\ddot {\epsilon }(t^*)$ with a coefficient $c_5$. Therefore, we have

(4.15)\begin{equation} \tilde{P} \sim \frac{ Sw_{b} (\tilde{h} + \widetilde{Sw}^3 ) \nu \sin(\phi_{a})}{ Sw_{b} \tilde{h} \nu \sin(\phi_{b})} \sim (1 + \tilde{A}^{{-}1}\widetilde{Sw}^3) \cdot \frac{\sin(\phi_{a})}{\sin(\phi_{b})}, \end{equation}

where $\tilde {A} = \tilde {h} = h_{p\_max} / h_{b\_max}$.

Overall, the scaling for power costs is

(4.16)\begin{equation} \tilde{P} \sim (1 + \tilde{A}^{{-}1}\widetilde{Sw}^3) \cdot \frac{\sin(\phi_{a})}{\sin(\phi_{b})}. \end{equation}

Equation (4.16) suggests that the power costs increase in AM is not only dependent on the input perturbations ($\,\tilde {f}$ and $\tilde {A}$), but also affected by the phase lag between the foil motion and the resulting fluid dynamic response. Even with the same input perturbations, the different flapping kinematics may generate different phase lags (figure 5b,e), therefore resulting in different power costs (figures 3b,4b).

4.3. Scaling for swimming efficiencies

Based on (4.6) and (4.16), the non-dimensionalised swimming efficiency $\tilde {C}_E$ can be scaled as

(4.17)\begin{equation} \tilde{C}_E = \tilde{P} \cdot \tilde{u}^{{-}1} \sim \frac{ 1 + \tilde{A}^{{-}1}\widetilde{Sw}^3}{(1+\widetilde{Sw}^{2})^{2/3}} \cdot \frac{\sin{(\phi_a)}}{\sin{(\phi_b)}}. \end{equation}

Equation (4.17) suggests that the efficiency improvement sourced from adding perturbations is affected by the non-dimensionalised amplitude, non-dimensionalised swimming number and phase lag. Referring to the research conducted by Floryan, Van Buren & Smits (Reference Floryan, Van Buren and Smits2018), the authors noted that the efficiency of a flapping foil is governed by both the flapping amplitude and the Strouhal number, a pattern akin to our (4.17). In particular, we demonstrated that the efficiency of a self-propelled foil can be enhanced by reducing the phase lag between the foil motion and the resulting fluid dynamics response.

To assess the compatibility of the scaling laws mentioned above with the simulated data of the self-propelled flapping foil, we created plots of $\tilde {u}$, $\tilde {P}$ and $\tilde {C}_E$ against the parameters associated with perturbations: $(1+\widetilde {Sw}^{2})^{2/3}$, $1 + \tilde {A}^{-1}Sw^{3}$ and $(1 + \tilde{A}^{-1}Sw^{3})/(1+\widetilde {Sw}^{2})^{2/3}$. These are presented in figure 6(a–c) for the self-propelled pitching mode and in figure 6(d–f) for the self-propelled heaving mode. The numerical data align remarkably well on a single linear curve, suggesting that the current scaling laws can accurately predict the hydrodynamic performance of self-propelled foils undergoing flapping with superimposed perturbations.

Figure 6. Scaling for the (a) $\tilde {u}$, (b) $\tilde {P}$ and (c) $\tilde {C}_E$ in the self-propelled pitching mode and (d) $\tilde {u}$, (e) $\tilde {P}$ and (f) $\tilde {C}_E$ in the self-propelled heaving mode, respectively.

5. Wake structures

Several distinguished wake structures generated by the self-propelled pitching foil are identified: reverse Kármán vortex (RKV) wake, symmetric wake (SW) I, symmetric wake (SW) II-A, symmetric wake (SW) II-B, asymmetric wake (AsW) I, asymmetric wake (AsW) II and 2P wake (Williamson & Roshko Reference Williamson and Roshko1988). The representative flow structures and their presence in the $\tilde {f} - \tilde {A}$ domain are shown in figure 7. The classical RKV wake is observed when the pitching foil employs BM, where the vortices from a side of the foil take the opposite side of the wake (figure 7b). When $\tilde {f}$ is an odd number ($\,\tilde {f} = 5,\ 7,\ 9$), the symmetric wake, including SW I, SW II-A and SW II-B, emerges behind the foil depending on both $\tilde {f}$ and $\tilde {A}$. The SW I forming as the RKV wake is observed at $\tilde {f} = 5$ and $(\,\tilde {f}, \tilde {A}) = (7, 0.01)$ (figure 7c). The emergence of the SW I is found at $(\,\tilde {f}, \tilde {A}) = (7, 0.02)$ and ($\,\tilde {f} = 9$, $\tilde {A} \leq 0.05$), where the wake involves the primary vortex (PV) and secondary vortex (SV) (figure 7d). The PV consists of a stronger main vortex pair occupied the position behind the foil's symmetric line, and weaker negative/positive vorticity located at the upper and bottom sides of the main vortex pair, respectively, while the SV forms the RKV wake. With a little increase in $\tilde {f}$ and /or $\tilde {A}$, the SW II-A transmutes to SW II-B (figure 7e), where the SV transits from the reverse Kármán vortex wake to the Kármán vortex (KV) wake. The asymmetric wake (AsW I and II) and 2P wake are observed when $\tilde {f}$ is an even number ($\,\tilde {f} = 4,\ 6,\ 8,\ 10$). The AsW I dominates the majority of the $\tilde {f} - \tilde {A}$ domain, manifesting as a downward deflection of the RKV wake (figure 7f). When $(\,\tilde {f}, \tilde {A}) = (8, 0.10)$, the vortices close to the foil's trailing edge fail to rapidly coalesce into a stronger vortex core, giving rise to the PV resembling 2P wake. Conversely, vortices farther from the foil's trailing edge evolve and experience viscous effects, forming a deflected RKV wake (figure 7g). The 2P wake is observed at $(\,\tilde {f}, \tilde {A}) = (10, 0.10)$. It is noted that the deflection of wake structure is influenced by whether $\tilde {f}$ is an odd or even number. When $\tilde {f}$ is an odd number, PM and BM are in phase, meaning that perturbations amplify the maximum flapping amplitude at $t/T_b = 0.25$. However, the flapping amplitude under AM is equal to that under BM when $\tilde {f}$ is an even number and $t/T_b = 0.25$, while the emergency of the maximum flapping amplitude is advanced (figure 8).

Figure 7. (a) Wake structure map in the self-propelled pitching mode. Typical instantaneous vorticity structures: (b) reverse Kármán vortex (BM); (c) symmetric wake (SW) I ($\,\tilde {f} = 5$, $\tilde {A} = 0.05$); (d) symmetric wake (SW) II-A ($\,\tilde {f} = 9$, $\tilde {A} = 0.04$); (e) symmetric wake (SW) II-B ($\,\tilde {f} = 7$, $\tilde {A} = 0.08$); (f) asymmetric wake (AsW) I ($\,\tilde {f} = 4$, $\tilde {A} = 0.05$); (g) asymmetric wake (AsW) II ($\,\tilde {f} = 8$, $\tilde {A} = 0.10$); (h) 2P wake ($\,\tilde {f} = 10$, $\tilde {A} = 0.10$).

Figure 8. (a) Varying of foil's pitching amplitude of BM and PM; (b) varying of foil's pitching amplitude of BM and AM.

Compared with the wake structures generated by the self-propelled pitching foil, the relationship between flow patterns provided by the self-propelled heaving foil and superimposed perturbations ($\,\tilde {f}$, $\tilde {A}$) is intricate (figure 9a). When the heaving foil undergoes the BM, it generates a symmetric wake composed of the PV, SV and tertiary vortex (TV): the PV takes shape as a 2P wake, the SV incorporates a main KV wake positioned behind the foil's symmetry line, with weaker positive/negative vorticity located on the upper and lower sides of the main KV wake, respectively, and the TV manifests as the RKV wake (figure 9b). As the flow patterns generated by the self-propelled heaving foil are complicated, we considered the symmetry of wake structures only from the PV. Only PV and SV are considered when the wake structures produced by the AM are identified. When the self-propelled heaving foil generates the symmetric wake (SW I, II and III), the PV constantly forms as the 2P wake, while the SV can be the deflected RKV wake (SW I, figure 9c), the main RKV wake positioned behind the foil's symmetry line (SW II, figure 9d), and the chaotic wake (SW III, figure 9e). The AsW describes the bifurcation flow deflected upwards or downwards (figure 9f). The CW (chaotic wake) is obscure, and we cannot identify meaningful $m$P$+ n$S vortex structures (Williamson & Roshko Reference Williamson and Roshko1988) from this flow field. Furthermore, The influence of the odd or even nature of $\tilde {f}$ on the flow structures is only significant when $\tilde {A} = 0.01$.

Figure 9. (a) Wake structure map in the self-heaving pitching mode. Typical instantaneous vorticity structures: (b) symmetric wake (BM); (c) symmetric wake (SW) I ($\,\tilde {f} = 8$, $\tilde {A} = 0.01$); (d) symmetric wake (SW) II ($\,\tilde {f} = 7$, $\tilde {A} = 0.01$); (e) symmetric wake (SW) III ($\,\tilde {f} = 9$, $\tilde {A} = 0.06$); (f) asymmetric wake (AsW) ($\,\tilde {f} = 5$, $\tilde {A} = 0.04$); (g) chaotic wake (CW) II ($\,\tilde {f} = 9$, $\tilde {A} = 0.10$).

In previous studies, researchers reported a smooth increase in propulsive force with an increase in flapping frequency and/or amplitude, accompanied by a smooth transition of the flow structures provided by the tethered flapping foil (KV $\rightarrow$ RKV $\rightarrow$ deflected RKV wake) (Godoy-Diana, Aider & Wesfreid Reference Godoy-Diana, Aider and Wesfreid2008; Andersen et al. Reference Andersen, Bohr, Schnipper and Walther2017; Chao, Alam & Ji Reference Chao, Alam and Ji2021). However, we found that the swimming speed produced by the self-propelled pitching/heaving foil with perturbations smoothly increases with the increment of $\tilde {f}$ and $\tilde {A}$, while the transition of the wake structures is not smooth. This observation raises a question: Could swimmers’ wake structures be reliable for predicting swimming performance? In the studies on the tethered flapping foil, some researchers suggested that the wake structures can be the footprints of a flapping foil and be used to predict the force generation (Zhang Reference Zhang2017), while others suggested that the swimmers’ wake structures are not reliable indicators of swimming performance (Floryan, Van Buren & Smits Reference Floryan, Van Buren and Smits2020). Due to our results indicating that the speed generation of self-propelled flapping foil largely aligns with the predictions of linear theory (Garrick Reference Garrick1936), we argue that the mode of the vortex wakes behind a self-propelled flapping foil with perturbations seems to be a remnant of speed generation, rather than a driver (Mackowski & Williamson Reference Mackowski and Williamson2015). Wake structures generated by the self-propelled flapping foil are shown in Appendix B (figures 12 and 13) and supplementary movie 1 available at https://doi.org/10.1017/jfm.2024.262.