1. Introduction
In nature, birds, insects and fish utilize flapping wings or fins to generate thrust for locomotion. Related issues intrigue lots of researchers and considerable breakthroughs and advances have been achieved over the past few decades through experimental measurements, theoretical analysis and numerical modelling (Triantafyllou, Triantafyllou & Yue Reference Triantafyllou, Triantafyllou and Yue2000; Lauder Reference Lauder2015; Smits Reference Smits2019). These advances are useful not only to biologists for better understanding of the underlying biology of fish and aquatic mammals, but also to engineers for the design of efficient biomimetic underwater vehicles (Dai et al. Reference Dai, He, Zhang and Zhang2018; Smits Reference Smits2019; Lin, Wu & Zhang Reference Lin, Wu and Zhang2021).
Previous studies have focused on the scaling laws of swimming performances (e.g. thrust and power expenditure). Quinn et al. (Reference Quinn, Moored, Dewey and Smits2014) investigated a rigid tethered pitching foil near a solid boundary through experiments and numerical simulations at a chord-based Reynolds number (
$Re$) of approximately 4700. They have shown that the thrust 
$T$ can be scaled as 
$(\,f A_t)^2$ and the power input 
$P$ can be scaled as 
$(\,f A_t)^{2.7}$ if the foil is far from the ground, where 
$f$ and 
$A_t$ are the flapping frequency and the amplitude of the trailing edge, respectively. Floryan et al. (Reference Floryan, Van Buren, Rowley and Smits2017) studied rigid tethered pitching and heaving foils by experiments on a nominally two-dimensional flow with 
$Re = 4780$. They found that the foil performances depend on the Strouhal number and reduced frequency. Further, they proposed that the thrust scaling is 
$T \sim (\,f A_t)^2$ if the offset drag is negligible (Floryan, Van Buren & Smits Reference Floryan, Van Buren and Smits2018), which is consistent with Quinn et al. (Reference Quinn, Moored, Dewey and Smits2014). They also suggested that the power scaling is 
$P \sim f L( V^2 - V_h V_{\theta })$ (Floryan et al. Reference Floryan, Van Buren and Smits2018), where 
$L$ is the length of the propulsor, 
$V = f A_t$ is the characteristic speed of the transverse motion, and 
$V_h = f h_0$ and 
$V_{\theta } = f L \theta _0$ are the transverse velocity scales characteristic of the heaving and pitching motions, respectively. The thrust scaling 
$T \sim (\,fA_t)^2$ is also derived by Gazzola, Argentina & Mahadevan (Reference Gazzola, Argentina and Mahadevan2014) through theoretical analysis and they also found that the cruising Reynolds number 
$Re_c$ can be scaled as 
$Re_f^{4/3}$, where 
$Re_c = U L/\nu$ based on the cruising speed 
$U$, and 
$Re_f = 2 {\rm \pi}f A_t L/ \nu$ is the flapping Reynolds number (also termed as the swimming number) (Vandenberghe, Zhang & Childress Reference Vandenberghe, Zhang and Childress2004; Gazzola et al. Reference Gazzola, Argentina and Mahadevan2014). Here 
$\nu$ is the kinematic viscosity of the fluid. By numerical simulations of rigid unconstrained pitching foils, Lin et al. (Reference Lin, Wu and Zhang2021) indicated that 
$Re_c \sim Re_f^{5/3}$ based on a new thrust scaling 
$T \sim (\,f A_t)^{5/2}$. Nevertheless, this new thrust scaling has never been confirmed or validated yet.
However, the rigid foils are quite different from the real wings/fins of birds/fish. These flapping wings/fins undergo large active or passive deformation in nature (Wootton Reference Wootton1992; Lauder Reference Lauder2015). Evidence has shown that flexibility has significant effect on the performance of the propeller (Thiria & Godoy-Diana Reference Thiria and Godoy-Diana2010; Kang et al. Reference Kang, Aono, Cesnik and Shyy2011; Marais et al. Reference Marais, Thiria, Wesfreid and Godoy-Diana2012; Shoele & Zhu Reference Shoele and Zhu2012; Zhu, He & Zhang Reference Zhu, He and Zhang2014b; Floryan & Rowley Reference Floryan and Rowley2018; Peng, Huang & Lu Reference Peng, Huang and Lu2018a). For instance, the thrust for the tethered pitching flexible foil is found to be up to three times larger than that of the rigid foil (Marais et al. Reference Marais, Thiria, Wesfreid and Godoy-Diana2012) and for the self-propulsive heaving flexible plate, the propulsive efficiency is much higher under certain conditions compared with the rigid plate (Zhu et al. Reference Zhu, He and Zhang2014b). In addition, the self-propelled flexible flapping wing may take on a more aerodynamic shape to achieve high efficiency (Ramananarivo, Godoy-Diana & Thiria Reference Ramananarivo, Godoy-Diana and Thiria2011). Dewey et al. (Reference Dewey, Boschitsch, Moored, Stone and Smits2013) experimentally studied the performances of pitching flexible panels at 
$Re \approx 7200$ and found that the global maximum efficiency for the flexible panels is twice as much as that of the rigid panels.
Besides, the drag experienced by a flexible body is quite different from that of a rigid body. Alben, Shelley & Zhang (Reference Alben, Shelley and Zhang2002, Reference Alben, Shelley and Zhang2004) showed that a flexible body experiences a drag proportional to 
$U^{4/3}$ by experiments and theoretical analysis with 
$Re$ in the range 
$2000 - 40\,000$. Neither 
$U^{3/2}$ scaling of a flat plate aligned with the flow (Gazzola et al. Reference Gazzola, Argentina and Mahadevan2014), nor classical 
$U^2$ drag scaling of rigid bodies (Batchelor Reference Batchelor1967) is applicable to the flexible body. This means that a flexible body can achieve drag reduction through self-similar bending or shape deformations due to flexibility (Alben et al. Reference Alben, Shelley and Zhang2002). Later, Zhu (Reference Zhu2008) numerically studied the same problem as Alben et al. (Reference Alben, Shelley and Zhang2002) and found that the 
$4/3$ power law of drag also occurs at lower 
$Re$ (
$\approx 800$). The 
$U^{4/3}$ drag scaling of a flexible body has also been confirmed by other studies further (Gosselin, de Langre & Machado-Almeida Reference Gosselin, de Langre and Machado-Almeida2010; Luhar & Nepf Reference Luhar and Nepf2011).
It is noteworthy that, in the experimental studies mentioned above, the object is fixed in an oncoming flow and cannot propel themselves freely. Although Van Buren et al. (Reference Van Buren, Floryan, Wei and Smits2018) suggested that tethered studies can be used to make robust conclusions about swimming performance, the self-propelled study may be better since the dynamic response of the swimmer to the surrounding flow is considered. In the present study, the performances of self-propelled flexible plates are investigated numerically. The plate can move freely in the propulsion direction. Two locomotion styles are considered, namely, the pitching and heaving motions. Meanwhile, the bending stiffness, the flapping amplitude and frequency are variable and their corresponding effects are investigated. We mainly focus on the scaling laws for propulsive speed, thrust and input power of the plates.
The remainder of this paper is organized as follows. The physical problem and mathematical formulation are presented in § 2. The numerical method and validation are described in § 3. Detailed results are discussed in § 4 and concluding remarks are addressed in § 5.
2. Physical problem and mathematical formulation
The schematic diagram of pitching and heaving flexible plates that we considered are shown in figure 1. The leading edge of the plate is forced to pitch (figure 1a) and heave (figure 1b) in the lateral direction, respectively. The motions are described as
$$\begin{gather} \theta(t)=\theta_{0} \sin (2 {\rm \pi}f t), \end{gather}$$
$$\begin{gather}h(t)=h_{0} \sin (2 {\rm \pi}f t), \end{gather}$$
where 
$\theta (t)$ and 
$h(t)$ are the instantaneous pitching and heaving motions, respectively, 
$\theta _{0}$ and 
$h_{0}$ are the pitching and heaving amplitudes, respectively, and 
$f$ is the flapping frequency.
Figure 1. Schematic diagram of the (a) pitching and (b) heaving flexible plates. The red curved lines represent the flexible plates and the thick dashed lines represent the rigid plates. Here 
$\theta _0$ and 
$h_0$ are the active pitching and heaving amplitudes, respectively, and 
$\theta _p$ is the passive pitching angle associated with the deformation of the flexible plate.
Note that only the leading edge of the plate is restricted with the prescribed lateral motion while other parts of the plate can move freely due to fluid–structure interaction. It is noticed that for the heaving motion, there is no active pitching but the plate may pitch passively.
The incompressible Navier–Stokes equations are adopted to solve the fluid flow:
$$\begin{gather} \frac{\partial \boldsymbol{v}}{\partial t}+\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v}={-}\frac{1}{\rho} \boldsymbol{\nabla} p+\frac{\mu}{\rho} \nabla^{2} \boldsymbol{v}+\boldsymbol{f}_{b}, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v}=0, \end{gather}$$
where 
$\boldsymbol {v}$ is the velocity, 
$p$ is the pressure, 
$\rho$ is the density of the fluid, 
$\mu$ is the dynamic viscosity and 
$\boldsymbol {f}_{b}$ is the Eulerian momentum force acting on the surrounding fluid due to the immersed boundary. To enforce the no-slip boundary condition, a virtual force is applied to each Lagrange point of the immersed boundary (Peskin Reference Peskin2002; Mittal & Iaccarino Reference Mittal and Iaccarino2005). The virtual force is distributed to the surrounding neighbouring fluid nodes through a delta function, i.e. at a relevant fluid node, there is a non-zero 
$\boldsymbol {f}_{b}$. It is noticed that at most fluid nodes far away from the Lagrange points, 
$\boldsymbol {f}_{b}=0$.
The structural equation is employed to describe the deformation and motion of the plate (Zhu & Peskin Reference Zhu and Peskin2002; Connell & Yue Reference Connell and Yue2007; Hua, Zhu & Lu Reference Hua, Zhu and Lu2013):
\begin{equation} \rho_{l} \frac{\partial^{2} \boldsymbol{X}}{\partial t^{2}}-\frac{\partial}{\partial s} \left[E h\left(1-\left|\frac{\partial \boldsymbol{X}} {\partial s}\right|^{{-}1}\right) \frac{\partial \boldsymbol{X}}{\partial s}\right]+E I \frac{\partial^{4} \boldsymbol{X}}{\partial^{4} s}=\boldsymbol{F}_{s}, \end{equation}
where 
$s$ is the Lagrangian coordinate along the plate, 
$\rho _{l}$ is the structural linear mass density, 
$\boldsymbol {X}(s, t)=(X(s, t), Y(s, t))$ is the position vector of the plate, 
$\boldsymbol {F}_{s}$ is the Lagrangian force exerted on the plate by the surrounding fluid. Here 
$Eh$ and 
$EI$ denote the structural stretching rigidity and bending rigidity, respectively. At the leading edge of the plate, for the heaving motion, the clamped boundary condition is adopted, i.e.
\begin{equation} -E h\left(1-\left|\frac{\partial \boldsymbol{X}}{\partial s}\right|^{{-}1}\right) \frac{\partial X}{\partial s}+E I\frac{\partial^{3} X}{\partial^{3} s}=0,\quad Y(t)=h(t),\quad \frac{\partial \boldsymbol{X}}{\partial s}=(1,0), \end{equation}while for the pitching motion, the boundary condition is expressed as
\begin{equation} -E h\left(1-\left|\frac{\partial \boldsymbol{X}}{\partial s}\right|^{{-}1}\right) \frac{\partial X}{\partial s}+E I\frac{\partial^{3} X}{\partial^{3} s}=0,\quad Y(t)=0,\quad \frac{\partial \boldsymbol{X}}{\partial s}=(\cos\theta,\sin\theta). \end{equation}Note that the first boundary condition in (2.6a–c) and (2.7a–c) represents the horizontally unconstrained condition. At the free end, the boundary condition is expressed as
\begin{equation} -E h\left(1-\left|\frac{\partial \boldsymbol{X}}{\partial s}\right|^{{-}1}\right) \frac{\partial \boldsymbol{X}}{\partial s}+EI \frac{\partial^{3} \boldsymbol{X}}{\partial^{3} s}=0,\quad \frac{\partial^{2} \boldsymbol{X}}{\partial s^{2}}=0. \end{equation}
In addition, 
$\boldsymbol {X}(s, 0)=(s, y(0)), \partial \boldsymbol {X} / \partial t(s, 0)=(0,0)$ are the initial conditions for the plate.
In our study, the fluid density 
$\rho$, the dynamic viscosity 
$\mu$, the dimensional length of the plate 
$L$ and the Reynolds number (
$Re=200$) are fixed. To normalize the above equations, the characteristic quantities 
$\rho$, 
$L$ and 
$U_{r e f}$ are chosen where 
$U_{ref} = \mu Re/(\rho L)$. Therefore, the characteristic time is 
$T_{r e f}=L / U_{r e f}$. Based on dimensional analysis, the following dimensionless governing parameters are introduced: the heaving amplitude 
$h_0^*$, the pitching amplitude 
$\theta _0$, the flapping frequency 
$f ^*$, the Reynolds number 
$Re=\rho U_{r e f} L / \mu$, the mass ratio of the plate to the fluid 
$M^*=\rho _{l} / (\rho L)$, the stretching stiffness 
$S^*=E h / (\rho U_{r e f}^{2} L)$ and the bending stiffness 
$K^*=E I / (\rho U_{r e f}^{2} L^{3})$.
3. Numerical method and validation
The governing equations of the fluid–plate problem are solved numerically by an immersed boundary-lattice Boltzmann method for the fluid flow and a finite element method for the motion of the flexible plate. The body force term 
$\boldsymbol {f}_{b}$ in (2.3) represents an interaction force between the fluid and the immersed boundary to enforce the no-slip velocity boundary condition. Equation (2.5) for the plate is discretized by a finite element method, and deformations with a large displacement of the plate are handled by the corotational scheme (Doyle Reference Doyle2001). More details on numerical methods can be found in our previous papers (Hua et al. Reference Hua, Zhu and Lu2013; Huang, Wei & Lu Reference Huang, Wei and Lu2018; Zhang, Huang & Lu Reference Zhang, Huang and Lu2020).
To validate the numerical method, a single plate in isolated swimming (Zhu et al. Reference Zhu, He and Zhang2014a) was simulated with 
$Re=200$, 
$h_0^*=0.5$, 
$M^*=0.2$, 
$K^*=0.8$ and 
$S^*=1000$. In the simulations, the computational domain for fluid flow is chosen as 
$[-15,25] \times [-15,15]$ in the 
$x$ and 
$y$ directions, which is sufficiently large so that the blocking effects of the boundaries are eliminated. A constant pressure with 
$\boldsymbol {v}=0$ is imposed at all boundaries except for the outlet where 
$\partial \boldsymbol {v} / \partial x=0$ with constant pressure is imposed (Zou & He Reference Zou and He1997). Initially, the fluid velocity field is zero in the entire computational domain. In the 
$x$ and 
$y$ directions the mesh is uniform with spacing 
${\rm \Delta} x={\rm \Delta} y=0.01 L$, where 
$L$ is the dimensional length of the plate. The time step is 
${\rm \Delta} t=T_f/10\,000$ for the simulations of fluid flow and plate deformation, where 
$T_f=1/ f$ is the flapping period. Besides, a finite moving computational domain (Hua et al. Reference Hua, Zhu and Lu2013) is used in the 
$x$-direction to allow the plate to move for a sufficiently long time. As the plate travels one lattice in the 
$x$-direction, the computational domain is shifted, i.e. one layer is added at the inlet and another layer is removed at the outlet (Hua et al. Reference Hua, Zhu and Lu2013).
Figure 2(a) shows the streamwise velocity of the leading edge as a function of time. It is seen that the present result is consistent with that of Zhu et al. (Reference Zhu, He and Zhang2014a). The results of grid independence and time-step independence for the case of 
$Re_c=100$ are shown in figure 2(b). It is seen that 
${\rm \Delta} x / L=0.01$ and 
${\rm \Delta} t / T_f=0.0001$ are sufficient to achieve accurate results. The results of grid and time-step independence studies for higher 
$Re_c$ are presented in table 1. It is seen that 
${\rm \Delta} x / L=0.01$ and 
${\rm \Delta} t / T_f=0.0001$ are accurate enough for 
$Re_c = 350$, while 
${\rm \Delta} x = 0.0075$ and 
${\rm \Delta} t/T_f = 0.000075$ are better for 
$Re_c = 860$. Hence, in most of the present simulations (i.e. cases with 
$Re_c < 400$), 
${\rm \Delta} x / L=0.01$ and 
${\rm \Delta} t /T_f=0.0001$ were adopted. For the cases with larger 
$Re_c$ (
$> 400$), a finer mesh 
${\rm \Delta} x / L=0.0075$ and a smaller time step 
${\rm \Delta} t /T_f=0.000075$ were adopted.
Figure 2. (a) Validation for the case of a self-propelled heaving plate with the non-dimensional governing parameters: 
$Re=200$, 
$h_0^*=0.5$, 
$M^*=0.2$, 
$K^*=0.8$ and 
$S^*=1000$ (Zhu, He & Zhang Reference Zhu, He and Zhang2014a). (b) Grid independence and time-step independence studies for a self-propelled pitching plate with 
$K^* =1$, 
$f^* = 1.0$, 
$\theta _0 = 20^{\circ }$ and 
$Re_c = 100$. The streamwise velocity of the leading edge as a function of time is presented.
Table 1. Grid independence and time-step independence studies for two typical cases with larger 
$Re_c$. Case1: 
$Re_c = 350$ (
$K^*=1$, 
$f^*=1$, 
$h_0^* = 0.5$). Case2: 
$Re_c = 860$ (
$K^*=5$, 
$f^*=2$, 
$h_0^* = 0.5$). The cruising speed 
$U^*$ and input power 
$P^*$ are presented.
Besides, our numerical strategy used in this study has been validated and successfully applied to investigate many flow problems, such as the coupling performance of tandem flexible inverted flags in a uniform flow (Huang et al. Reference Huang, Wei and Lu2018), the effect of the trailing-edge shape on the self-propulsive performance of heaving flexible plates (Zhang et al. Reference Zhang, Huang and Lu2020) and the intermittent locomotion performance of a self-propelled flapping plate (Liu, Huang & Lu Reference Liu, Huang and Lu2020).
4. Results and discussion
In present simulations, three main parameters are fixed, namely, the Reynolds number 
$Re=200$, mass ratio 
$M^*=0.2$ and stretching stiffness 
$S^*=1000$, which are identical to those in the previous studies (Zhu et al. Reference Zhu, He and Zhang2014a; Peng et al. Reference Peng, Huang and Lu2018a; Liu et al. Reference Liu, Huang and Lu2020). Note that 
$S^*$ is large enough so that the plate is almost inextensible. While other key parameters are variable, i.e. the bending stiffness 
$K^* \in [0.5,10]$, the heaving amplitude 
$h_0^* \in [0.1,0.5]$, the pitching amplitude 
$\theta _0 \in [5,30]$ (deg.) and the flapping frequency 
$f^* \in [0.5,2]$. It is noted that for real fish, 
$K^* \approx O(1)$, e.g. 
$K^*$ of the tail fin of a goldfish (
$Carassius\ auratus$) is within the range of 
$2.5 - 23$ (Hua et al. Reference Hua, Zhu and Lu2013; Peng, Huang & Lu Reference Peng, Huang and Lu2018b). Hence, 
$K^* \in [0.5,10]$ is adopted here. In the following discussion, variables with a superscript ‘
$*$’ are dimensionless variables. The symbol ‘
$\sim$’ means that the two variables on the left and right-hand sides are proportional, and their dimensions are not necessarily consistent.
4.1. Propulsive performances and scaling laws
The mean propulsive speed 
$U$, thrust 
$T$ and input power 
$P$ are paid special attention. Specifically, 
$U$ is calculated by the time-averaged speed of the leading edge (
$s = 0$) within one cycle at the equilibrium state (Zhu et al. Reference Zhu, He and Zhang2014a; Peng et al. Reference Peng, Huang and Lu2018a), i.e.
\begin{equation} U=\frac{1}{T_f} \int_{t^{\prime}}^{t^{\prime}+T_f} u(t)\,\mathrm{d} t={-}\frac{1}{T_f} \int_{t^{\prime}}^{t^{\prime}+T_f}\left(\left.\frac{\partial X}{\partial t}\right|_{s=0}\right) \mathrm{d} t, \end{equation}
Here u(t) is the instantaneous horizontal speed of the leading edge of the platea and 
$P$ is defined as (Zhu et al. Reference Zhu, He and Zhang2014a)
\begin{equation} P=\frac{1}{T_f} \int_{t^{\prime}}^{t^{\prime}+T_f} P(t)\,\mathrm{d} t=\frac{1}{T_f} \int_{t^{\prime}}^{t^{\prime}+T_f}\left[\int_{0}^{1} \boldsymbol{F}_{r}(s, t) \cdot \frac{\partial \boldsymbol{X}(s, t)}{\partial t} \mathrm{d} s\right] \mathrm{d} t, \end{equation}
where 
$\boldsymbol {F}_{r}$ represents the force on the surrounding fluid by the plate.
The thrust and drag are defined based on force decomposition. The Lagrangian force 
$\boldsymbol {F}_{s}$, exerted on the plate by the surrounding fluid, can be decomposed into two parts: one is the normal force 
$\boldsymbol {F}^{n}$, in which the pressure component dominates; the other is the tangential force 
$\boldsymbol {F}^{\tau }$, which mainly comes from the viscous effects. These forces at the 
$i$th Lagrangian node are defined as follows (Peng et al. Reference Peng, Huang and Lu2018a; Liu et al. Reference Liu, Huang and Lu2020):
$$\begin{gather} \boldsymbol{F}_{s,i}=[{-}p \boldsymbol{I}+\boldsymbol{T}] \boldsymbol{\cdot} \boldsymbol{n}=\boldsymbol{F}_i^n+\boldsymbol{F}_i^{\tau}, \end{gather}$$
$$\begin{gather}\boldsymbol{F}_i^{n}=\left(\boldsymbol{F}_{s,i} \boldsymbol{\cdot}\boldsymbol{n}\right) \boldsymbol{n}=\left(F_{x,i}^{n},F_{y,i}^{n}\right), \end{gather}$$
$$\begin{gather}\boldsymbol{F}_i^{\tau}=\left(\boldsymbol{F}_{s,i} \boldsymbol{\cdot} \boldsymbol{\tau}\right) \boldsymbol{\tau}=\left(F_{x,i}^{\tau}, F_{y,i}^{\tau}\right), \end{gather}$$
where 
$\boldsymbol {I}$ is the unit tensor, 
$\boldsymbol {T}$ is the viscous stress tensor, 
$\boldsymbol {\tau }$ is the unit tangential vector toward the trailing edge, 
$\boldsymbol {n}$ is the unit normal vector and 
$[\cdot ]$ denotes the jump in a quantity across the immersed boundary. Figure 3(a) shows the schematic diagram for force decomposition. The total 
$x$-component of 
$\boldsymbol {F}^{n}$ and 
$\boldsymbol {F}^{\tau }$ are 
$F_{x}^{n}=\sum _{i} F_{x, i}^{n}$ and 
${F}_x^\tau =\sum _{i} F_{x, i}^{\tau }$, respectively. Figure 3(b) presents the dimensionless forces 
$F_{x}^{n*}$ and 
$F_{x}^{\tau *}$ as functions of time within one cycle. It is seen that, for the heaving plate, 
$F_x^n$ contributes much to a thrust since 
$F_x^n<0$, while 
$F_x^{\tau }$ contributes to a drag for 
$F_x^{\tau } >0$. However, for the pitching plates, 
$F_x^n$ may contribute to a drag and 
$F_x^{\tau }$ may contribute to a thrust at some time. Hence, the thrust 
$T$ and drag 
$D$ can be defined as (Bottom et al. Reference Bottom, Borazjani, Blevins and Lauder2016)
\begin{equation} \left.\begin{gathered} -T(t) =\tfrac{1}{2}\left(F_{x}^{n}-\left|F_{x}^{n}\right|\right)+ \tfrac{1}{2}\left(F_{x}^{\tau}-\left|F_{x}^{\tau}\right|\right), \\ D(t) =\tfrac{1}{2}\left(F_{x}^{n}+\left|F_{x}^{n}\right|\right)+ \tfrac{1}{2}\left(F_{x}^{\tau}+\left|F_{x}^{\tau}\right|\right). \end{gathered}\right\} \end{equation}
The time-averaged thrust and drag are defined as 
$T=({1}/{T_f}) \int _{t^{\prime }}^{t^{\prime }+T_f} T(t)\,\mathrm {d} t$ and 
$D=({1}/{T_f}) \int _{t^{\prime }}^{t^{\prime }+T_f} D(t)\,\mathrm {d} t$, respectively.
Figure 3. (a) Schematic diagram for force decomposition. Black curved line represents the plate. Red dot represents the 
$i$th Lagrangian node. Here 
$\boldsymbol {\tau }$ and 
$\boldsymbol {n}$ denote the local tangential and normal vectors, respectively, and 
$\theta$ is the angle between 
$\boldsymbol {\tau }$ and the positive direction of the 
$x$-axis. (b) The total dimensionless 
$x$-component of 
$\boldsymbol {F}^{n}$ and 
$\boldsymbol {F}^{\tau }$ (i.e. 
$F_{x}^{n *}=\sum _{i} F_{x, i}^{n *}$ and 
${F}_x^{\tau *}=\sum _{i} F_{x, i}^{\tau *}$) for the pitching (
$K^* =1, \theta _0 = 20^{\circ }, f^*=1.0$) and heaving (
$K^*=1, h_0^* = 0.5, f^*=1.0$) plates as functions of time within one flapping period. Solid and dashed lines represent 
$F_{x}^{n*}$ and 
$F_{x}^{\tau *}$, respectively. The forces are normalized by 
$F_{r e f}=\frac {1}{2} \rho U_{r e f}^{2} L$.
The propulsive performances of the plates with the bending stiffness 
$K^*=1$ are investigated first, as shown in figure 4. It is seen that 
$U$, 
$T$ and 
$P$ all increase with the increase of the flapping amplitude and frequency for both pitching and heaving plates. This is because that, generally, larger flapping amplitude and frequency mean that the plate may consume more energy, achieve larger thrust and be propelled faster. Besides, we have checked parameters of all cases in the present study and found that, for most of the cases, the flapping frequency 
$f$ is significantly below the first natural frequency 
$f_1$ of the system. For the other few cases, even when the frequency ratio 
$f_r={f}/{f_1}$ is larger than unity, it is close to unity. Hence, our result shows that when 
$f< f_1$, the speed, thrust and power increase monotonically with frequency. That is consistent with the results in the literature (Ramananarivo et al. Reference Ramananarivo, Godoy-Diana and Thiria2011).
Figure 4. (a) Dimensionless propulsive speed 
$U^*$, (b) dimensionless thrust 
$T^*$ and (c) dimensionless input power 
$P^*$ as functions of flapping amplitude for pitching plates with different flapping frequency. Panels (d–f) are the corresponding results for heaving plates. In all these cases, 
$K^* =1$. The velocity, force and power are normalized by 
$U_{ref} = \mu Re/(\rho L)$, 
$F_{ref}= (1/2) \rho U_{ref}^{2} L$ and 
$P_{ref}= \rho U_{ref}^{3} L$, respectively.
Actually, besides the effects of 
$\theta _0$ (
$h_0^*$) and 
$f^*$, the bending stiffness 
$K^*$ may also have a significant effect on the performances of the plate (Peng et al. Reference Peng, Huang and Lu2018a; Liu et al. Reference Liu, Huang and Lu2020). Hence, it is necessary to adopt new parameters to concisely describe the results with different 
$K^*$, 
$h_0^*$ or 
$\theta _0$ and 
$f^*$. Inspired by previous works (Vandenberghe et al. Reference Vandenberghe, Zhang and Childress2004; Gazzola et al. Reference Gazzola, Argentina and Mahadevan2014; Lin et al. Reference Lin, Wu and Zhang2021), the flapping Reynolds number 
$Re_f = 2{\rm \pi} f A_t L/\nu$ is introduced (i.e. the swimming number 
$Sw$ in the study of Gazzola et al. Reference Gazzola, Argentina and Mahadevan2014). Besides, the cruising Reynolds number 
$Re_c = U L/\nu$, based on the cruising speed 
$U$, can be defined.
The cruising Reynolds number 
$Re_c$ as a function of 
$Re_f$ is presented in figure 5. It is astonishing to see that, for both pitching and heaving plates, 
$Re_c$ values of all cases follow the same scaling law approximately, i.e. 
$Re_{c} \sim Re_{f}^{3/2}$. This can be interpreted as follows. By a simple dimensional analysis reasoning (see Appendix A.1) and the elongated-body theory (EBT) (Lighthill Reference Lighthill1960, Reference Lighthill1971) (see Appendix A.2), we get that the thrust scaling is 
$T \sim (\,f A_t)^2$. This scaling is also identical to that in previous studies, e.g. Gazzola et al. (Reference Gazzola, Argentina and Mahadevan2014) and Quinn et al. (Reference Quinn, Moored, Dewey and Smits2014). As for the drag scaling, ample evidence (including theoretical analysis, experiments and simulations) indicates that, for a flexible body in an oncoming flow, 
$D$ is proportional to 
$U^{4/3}$ (Alben et al. Reference Alben, Shelley and Zhang2002; Zhu Reference Zhu2008; Luhar & Nepf Reference Luhar and Nepf2011). In Appendix B, we further show that this drag scaling is still valid for the inclined plates, which are more like our flexible plates. More detailed discussions about the thrust and drag scaling are presented in § 4.2. Balancing the thrust and drag, one can get that 
$U^{4/3} \sim (\,fA_t)^2$, i.e. 
$U \sim (\,fA_t)^{3/2}$. From the definitions of 
$Re_c$ and 
$Re_f$, the formula can be further written as 
$Re_c \sim Re_f^{3/2}$. The result is highly consistent with our observation from figure 5.
Figure 5. The cruising Reynolds number 
$Re_c$ as a function of the flapping Reynolds number 
$Re_f$ for (a) pitching and (b) heaving plates. Each symbol denotes a case that we simulated. Each kind of symbol represents the cases with a specific 
$K$. The green, red, blue and black symbols represent the cases with 
$f=0.5, 1.0, 1.5$ and 
$2.0$, respectively.
We would like to check whether the present scaling law is supported by the biological observations. Figure 6 shows that our scaling law is consistent with the biological data from the supplementary information of Gazzola et al. (Reference Gazzola, Argentina and Mahadevan2014). Meanwhile, the present scaling law seems to be applicable to a wider range of 
$Re_f$ (up to 
${\sim }2.5 \times 10^5$).
Figure 6. Biological data and the scaling 
$Re_c \sim Re_f^{3/2}$. The data sets are obtained from figures S4(d), S5(d), S8(c) and S10(a,b) in the supplementary information of Gazzola et al. (Reference Gazzola, Argentina and Mahadevan2014).
Next, our scaling laws for the thrust 
$T$ and the drag 
$D$ are discussed. Figure 7 plots the dimensionless thrust 
$T^*$ vs 
$Re_f$ and the dimensionless drag 
$D^*$ vs 
$U^*$. It is seen that there also exist some scaling laws. For the pitching plates, 
$T^*$ (figure 7a) is proportional to 
$Re_f^{1.8}$ and 
$D^*$ (figure 7c) is proportional to 
$U^{*1.2}$. While for heaving plates, 
$T^*$ and 
$D^*$ (see figure 7b,d) show approximately 
$Re_f^{2.1}$ and 
$U^{*1.4}$ growth, respectively. It is surprising to find that the thrust and drag scaling are very close to 
$T^* \sim Re_f^{2}$ and 
$D^* \sim U^{*4/3}$, respectively, which we adopted above. Besides, due to 
$T = D$ at equilibrium state, one can get that 
$Re_f^{1.8} \sim U^{*1.2}$ for the pitching plates and 
$Re_f^{2.1} \sim U^{*1.4}$ for the heaving plates, which can both be further transformed to 
$Re_c \sim Re_f^{3/2}$. This confirms the scaling law in figure 5 from another view. It is also noticed that there is a small discrepancy between the numerical and theoretical scalings. This may be due to the nonlinearity of the system and the complex fluid–structure interaction, which are greatly simplified in the theoretical analyses. For example, in the theoretical analysis (see Appendix A), we assumed that the effect of stiffness 
$K^*$ is fully included in 
$A_t$.
Figure 7. (a) The dimensionless thrust 
$T^*$ as a function of 
$Re_f$ and (c) the dimensionless drag 
$D^*$ as a function of 
$U^*$ for pitching plates. Panels (b,d) are the corresponding results of the heaving plates. Symbols are identical to those in figure 5. The velocity and force are normalized by 
$U_{ref} = \mu Re/\rho L$ and 
$F_{ref}= (1/2) \rho U_{ref}^{2} L$, respectively.
It is seen that the data sets of 
$T$ and 
$D$ for the pitching plates overlap very well (figure 7a,c). Note that due to the flexibility of the plate, the direct result is the emergence of the passive pitching angle 
$\theta _p$ (see figure 1). Nevertheless, the locomotion style has not changed for the pitching plates (compared with the heaving plates, see below). Under the same locomotion style, the growing forms of 
$T$ and 
$D$ should be identical. Therefore, all pitching data sets collapse into a single curve.
On the other hand, for the flexible heaving plates, there are both heaving (active, amplitude 
$h_0$) and pitching (passive, amplitude 
$\theta _p$) motions. When the plate is almost rigid, the heaving motion is prominent. When the plate becomes more flexible, the passive pitching motion becomes more prominent, i.e. there may be a large 
$\theta _p$. Actually, here 
$K^*$ plays an important role in the thrust generation. For the rigid plate, the thrust is mainly generated by the leading-edge vortex (Vandenberghe et al. Reference Vandenberghe, Zhang and Childress2004), while for the flexible plate, the normal force contributes much to the thrust (Peng et al. Reference Peng, Huang and Lu2018a; Liu et al. Reference Liu, Huang and Lu2020). Due to the different flow mechanism, the flexibility influences the thrust scaling considerably (and drag scaling, since 
$T = D$). Therefore, the data sets do not overlap very well. It is also noticed that 
$T$ at smaller 
$K^*$ (i.e. 
$K^*=0.5$ and 1, the square and circle symbols in figure 7b) is higher than that at larger 
$K^*$ (i.e. 
$K^*=5$ and 10). This is because the thrust can be significantly enhanced by appropriately increasing the flexibility (Marais et al. Reference Marais, Thiria, Wesfreid and Godoy-Diana2012).
Finally, the power expended is considered, as presented in figure 8. Results show that the dimensionless input power 
$P^*$ can be scaled as 
$P^* \sim Re_f^{2.7}$ for both pitching and heaving plates, which is consistent with the results of Quinn et al. (Reference Quinn, Moored, Dewey and Smits2014). Besides, as described in the introduction, Floryan et al. (Reference Floryan, Van Buren and Smits2018) showed that the power scaling is 
$P \sim fL( V^2 - V_h V_{\theta })$. Note that 
$V_h V_{\theta } =0$ in the present study, hence 
$P \sim fL V^2 \sim fL(\,f A_t)^2$. If the flapping amplitude is large, i.e. 
$A_t \sim L$, the formula becomes 
$P \sim (\,f A_t)^3$. The same power scaling can also be obtained by dimensional analysis (see Appendix A.1). In the derivation, if 
$A_t$ is large, the effect of lateral velocity 
$f A_t$ is more significant and therefore more suitable to be chosen as the characteristic velocity. The derivation from the EBT of Lighthill (Reference Lighthill1971) also confirms our power scaling (see Appendix A.2). On the other hand, for small amplitude, i.e. 
$A_t / L \ll 1$, the model of Theodorsen (Reference Theodorsen1935) indicates 
$P \sim (\,f A_t)^2$. It is noticed that in our study, 
$A_t/L \in [0.09, 0.55]$. Therefore, the amplitude in our study is moderate, which is between the conditions 
$A_t / L \ll 1$ and 
$A_t \sim L$. Therefore, approximately the exponent should between two and three. In the experimental study of rigid tethered pitching foils by Quinn et al. (Reference Quinn, Moored, Dewey and Smits2014), their flapping amplitude range is 
$A_t/L \in [0.095, 0.16]$ and they have 
$P \sim (\,f A_t)^{2.7}$. In general, our 
$A_t$ range is a little bit larger than that of Quinn et al. (Reference Quinn, Moored, Dewey and Smits2014). The exponent should be closer to three. However, compared with the rigid foils, our flexible plates should be favourable for energy consumption (Zhu et al. Reference Zhu, He and Zhang2014b). Hence, the power scaling of a flexible foil with larger 
$A_t$ may be close to that of the rigid foil with smaller 
$A_t$. In other words, our power scaling 
$(\,f A_t)^{2.7}$ is reasonable.
Figure 8. The dimensionless time-average input power 
$P^*$ as a function of 
$Re_f$ for the (a) pitching plates and (b) heaving plates. Symbols as in figure 5. The powers are normalized by 
$P_{ref}= \rho U_{ref}^{3} L$.
4.2. Comparison with results in the literature
In the above, we proposed a new scaling law, i.e. 
$Re_c \sim Re_f^{3/2}$. It is also noted that other two scaling laws, namely, 
$Re_c \sim Re_f^{4/3}$ and 
$Re_c \sim Re_f^{5/3}$, are given by Gazzola et al. (Reference Gazzola, Argentina and Mahadevan2014) and Lin et al. (Reference Lin, Wu and Zhang2021), respectively. In the following, we would like to make a thorough comparison between these scaling laws.
Indeed, the three scaling laws may look close in the qualitative analysis of a log–log plot. However, a quantitative analysis of our data fittings indicates that the optimal 
$Re_c$-scaling is 
$Re_c \sim Re_f^{3/2}$ with 
$R^2 = 0.99$ for the pitching cases (figure 5a) and 
$R^2 = 0.97$ for the heaving cases (figure 5b).
Besides, note that all three scaling laws can be derived from the balance of thrust and drag. Hence, a more essential analysis should start from thrust and drag scaling. Gazzola et al. (Reference Gazzola, Argentina and Mahadevan2014) have derived a scaling law of thrust 
$T$, i.e. 
$T \sim (\,f A_t)^2$. The thrust scaling has been found and confirmed further in other studies (Quinn et al. Reference Quinn, Moored, Dewey and Smits2014; Floryan et al. Reference Floryan, Van Buren and Smits2018; Gibouin et al. Reference Gibouin, Raufaste, Bouret and Argentina2018; Van Buren et al. Reference Van Buren, Floryan, Wei and Smits2018) and our physical reasoning (see Appendix A). In addition, the quantitative results of 
$T$ in figures 7(a) and 7(b) are also consistent well with this thrust scaling.
Under the circumstances of 
$T \sim (\,f A_t)^2$, the swimming speed scale is sensitive to the assumed drag law. For the drag scaling, Gazzola et al. (Reference Gazzola, Argentina and Mahadevan2014) suggested that 
$D \sim U^{3/2}$, which is estimated from the viscous drag on a flat and rigid plate aligned with the flow. However, for a real fish, due to its flexibility, when it flaps, its body continuously deforms. The curved boundary layer and the corresponding drag should be quite different from those of the rigid flat plate. Therefore, the drag scaling of Gazzola et al. (Reference Gazzola, Argentina and Mahadevan2014) may be not universally applicable. The flexible curved boundary-layer effect has been fully considered in previous studies and found that a stationary curved flexible fibre or plate in an oncoming flow experiences a drag proportional to 
$U^{4/3}$ (Alben et al. Reference Alben, Shelley and Zhang2002; Zhu Reference Zhu2008; Gosselin et al. Reference Gosselin, de Langre and Machado-Almeida2010; Luhar & Nepf Reference Luhar and Nepf2011). Generally speaking, 
$U^{4/3}$ scaling denotes a drag reduction for flexible curved plates (Alben et al. Reference Alben, Shelley and Zhang2002). Besides, as we have shown in Appendix B, the 
$U^{4/3}$ scaling is still valid for inclined plates, whose shapes are also curved and streamlined due to deformation. Therefore, in the present study, the drag of the deformed curved flapping plates should be close to that of plates in Alben et al. (Reference Alben, Shelley and Zhang2002), rather than that of flat plates. Moreover, it can be observed directly from figures 7(c) and 7(d) that the 
$U^{4/3}$ scaling of drag can fit our data well. Hence, it is adopted to perform the analysis.
In summary, there are strong grounds to support our scalings 
$T \sim (\,f A_t)^2$ and 
$D \sim U^{4/3}$, and further support the scaling 
$Re_c \sim Re_f^{3/2}$.
In the study of Lin et al. (Reference Lin, Wu and Zhang2021), the drag scaling of Gazzola et al. (Reference Gazzola, Argentina and Mahadevan2014), i.e. 
$D \sim U^{3/2}$, was adopted. They then derived a new thrust scaling, i.e. 
$T \sim (\,f A_t)^{5/2}$. However, the new thrust scaling has never been confirmed or validated yet. One might conjecture that the thrust scaling is only applicable to their rigid self-propulsive flapping model.
It should be pointed out that 
$Re$ in the present study is close to that of larvae. Although our 
$Re$ is much lower than those of many adult swimming animals, our scaling may extend to surprisingly higher 
$Re$ (up to 
${\sim }2.5 \times 10^5$, see figure 6). While 
$Re$ in the study of Gazzola et al. (Reference Gazzola, Argentina and Mahadevan2014) went up past 
$10^8$, around 1000 higher than what is presented here. Indeed, 
$Re$ may impact the nature of drag. If 
$Re$ is high enough, biological data may switch to a different scaling.
5. Conclusions
The propulsive performances of self-propelled flapping flexible plates are investigated numerically. Both pitching and heaving motions are considered. Results indicate that the propulsive speed 
$U$, thrust 
$T$ and input power 
$P$ increase with the flapping amplitude and frequency in both locomotion styles. Further, some simple scaling laws are found. Specifically, for both pitching and heaving plates, the cruising Reynolds number 
$Re_c$ can be scaled as 
$Re_c \sim Re_f^{3/2}$, and the dimensionless input power 
$P^*$ can be scaled as 
$P^* \sim Re_f^{2.7}$. The scaling laws of 
$T$ and 
$P$ are also derived using dimensional analysis reasoning and the EBT (Lighthill Reference Lighthill1971). The derived laws are consistent with or close to our findings. For the thrust and drag calculated in our simulations, the scaling laws are very close to 
$(\,fA_t)^2$ and 
$U^{4/3}$, respectively. They are consistent with those in the literature. The 
$U^{4/3}$ drag law is found to be still applicable for inclined plates. Besides, the flexibility influences the thrust scaling considerably for the heaving plates, while it does not for the pitching plates. The 
$Re_c$ scaling can be simply derived by balancing the thrust and drag. Our present scaling law is also well supported by biological data. This study may deepen our understanding of the locomotion of aquatic animals and may be helpful for bionic design.
Acknowledgements
We appreciate all the anonymous reviewers for their constructive suggestions and comments.
Funding
This work was supported by the Natural Science Foundation of China (NSFC) grant nos. 11972342, 11772326 and 11621202.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Physical reasoning for the scaling laws
A.1. Dimensional analysis reasoning for the thrust and power scalings
Dimensional analysis reasoning is a broadly applicable technique for developing scaling laws. In isothermal fluid mechanics problems, there are only three basic dimensions, i.e. mass 
$\mathcal {M}$, length 
$\mathcal {L}$ and time 
$\mathcal {T}$. Therefore, only three repeating parameters are required, and the determinant of the dimensional matrix formed from these three parameters must be non-zero (Kundu, Cohen & Dowling Reference Kundu, Cohen and Dowling2012). The common choices are characteristic length, velocity and density. In the present flow system, the length of the plate 
$L$ and the fluid density 
$\rho$ are naturally chosen as characteristic length and density, respectively. For the characteristic velocity, there are two choices, i.e. the cruising speed 
$U$ and 
$f A_t$ based on the trailing-edge amplitude 
$A_t$. However, 
$U$ is determined by the characteristic lateral velocity of the tail motion 
$f A_t$, i.e. the most important velocity scale is 
$f A_t$ instead of 
$U$ (Van Buren et al. Reference Van Buren, Floryan, Wei and Smits2018; Smits Reference Smits2019). Hence, the characteristic velocity 
$f A_t$ is adopted here.
It is noted that there are four variable dimensionless input parameters in our study, i.e. the bending stiffness 
$K^*$, the heaving amplitude 
$h_0^*$, the pitching amplitude 
$\theta _0$ and the flapping frequency 
$f^*$ (see § 4). Actually, the effects of these four parameters are directly or indirectly included in the characteristic velocity 
$f A_t$. Specifically, the flapping frequency 
$f$ appears explicitly in 
$f A_t$; the amplitude of the trailing edge 
$A_t$ is mainly affected by the flapping amplitude (
$h_0^*$ or 
$\theta _0$) and modulated by the bending stiffness 
$K^*$. From this point of view, 
$f A_t$ is indeed a key characteristic quantity of the system.
Hence, we assume that the effects of the flapping amplitude (
$h_0^*$ or 
$\theta _0$) and the rigidity (
$K^*$) can be all quantified by 
$A_t$ as a whole. Through this assumption, we can get a much simpler form of scalings using only three parameters (i.e. 
$\rho$, 
$L$ and 
$f A_t$). Besides, we finally get the same thrust scaling as that of Gazzola et al. (Reference Gazzola, Argentina and Mahadevan2014). Also note that in the experimental study about flexible flapping plates of Floryan, Van Buren & Smits (Reference Floryan, Van Buren and Smits2019), when 
$f$ is small (i.e. the Strouhal number 
$St=2f A_t/U_{\infty }$ is small), the data of thrust and power are well collapsed (see their figures 8(b), 9(b), 10(b) and 11(b)). In other words, 
$A_t$ can capture the effects of flexibility and describe the results well when 
$f$ is small. Therefore, this assumption is reasonable to a certain extent. Indeed, 
$A_t$ cannot fully represent the effect of 
$K^*$. As we have mentioned in § 4.1, this simplified analysis may lead to a small discrepancy between the numerical and theoretical scalings.
According to the 
$\varPi$ theorem of Buckingham (Reference Buckingham1914), a dimensionless group of thrust is expressed as
\begin{equation} \varPi_1 = T \rho^a (\,f A_t)^b L^c. \end{equation}
Note that, for a two-dimensional problem, the dimensions of 
$T$, 
$\rho$, 
$f A_t$ and 
$L$ are 
$[ T ] = [\textrm {N}\ \textrm {m}^{-1}]$, 
$[ \rho ] = [\textrm {kg}\ \textrm {m}^{-3}]$, 
$[\,f A_t] = [\textrm {m}\ \textrm {s}^{-1}]$ and 
$[ L ] = [\textrm {m}]$, respectively. Since 
$\varPi _1$ is dimensionless, one can easily obtain that the exponents 
$a=-1$, 
$b=-2$ and 
$c=-1$. Therefore,
\begin{equation} \varPi_1 =\frac {T} {\rho (\,f A_t)^2 L}. \end{equation}
Holding 
$\varPi _1$ constant, one can get that 
$T \sim (\,f A_t)^2$, i.e. 
$T^* \sim Re_f^2$ when written in dimensionless form.
As for the input power 
$P$, the dimensionless group is
\begin{equation} \varPi_2 = P \rho^a (\,f A_t)^b L^c. \end{equation}
The dimension of 
$P$ is 
$[ P ] = [\textrm {kg}\ \textrm {m}\ \textrm {s}^{-3}]$ for a two-dimensional problem. Because 
$\varPi _2$ is dimensionless, it is easy to find that 
$a = -1$, 
$b = -3$ and 
$c = -1$. Hence,
\begin{equation} \varPi_2 =\frac {P} {\rho (\,fA_t)^3 L}. \end{equation}
Holding 
$\varPi _2$ constant, we have 
$P \sim (\,f A_t)^3$, which can be written in non-dimensional form as 
$P^* \sim Re_f^{3}$.
A.2. Elongated-body theory for the scaling laws
Elongated-body theory has been widely adopted in previous studies about two- and three-dimensional problems (Bottom et al. Reference Bottom, Borazjani, Blevins and Lauder2016; Floryan et al. Reference Floryan, Van Buren and Smits2018; Akbarzadeh & Borazjani Reference Akbarzadeh and Borazjani2019). For an undulating plate, the lateral motion of the whole body is expressed by the travelling-wave equation
\begin{equation} Y(z,t) = A(z)\sin(2{\rm \pi} ft - 2{\rm \pi} z / \lambda), \end{equation}
where 
$A(z)$ is the amplitude of the wave, 
$z$ is the axial direction from the leading edge (
$z=0$) to the trailing edge (
$z=L$), 
$f$ is the frequency and 
$\lambda$ is the wavelength. According to the EBT (Lighthill Reference Lighthill1960, Reference Lighthill1971), the thrust of an undulating plate is expressed as
\begin{equation} T=\left[m w\left(W-\frac{1}{2} w\right)\right]_{z=L}-\frac{\partial}{\partial t} \int_{0}^{L}\left(m w \frac{\partial Y}{\partial z}\right) \mathrm{d} z, \end{equation}
where 
$m$ is the virtual mass per unit length (Lighthill Reference Lighthill1971; Akbarzadeh & Borazjani Reference Akbarzadeh and Borazjani2019). Note that for a three-dimensional body, 
$m = \frac {1}{4} {\rm \pi}\beta \rho s^2$, where 
$\beta \approx 1$ is a non-dimensional parameter, 
$\rho$ is the fluid density and 
$s$ is the span (Lighthill Reference Lighthill1971). For two-dimensional problems, although 
$m$ may have a different expression, it is independent of kinematic parameters (i.e. 
$f$ and 
$A_t$), i.e. we can take 
$m$ as a constant. Hence, 
$m$ does not affect the scaling laws. Here 
$W = {\partial Y }/ {\partial t} = 2{\rm \pi} f A(z)\cos (2{\rm \pi} ft - 2{\rm \pi} z / \lambda )$ is the lateral velocity of the plate, and 
$w = {\textrm {D} Y}/{\textrm {D} t} = {\partial Y}/{\partial t}+ U ({\partial Y}/{\partial z})$ is the lateral velocity of fluid adjacent to the plate. Using (A5) and assuming the variation of 
$A(z)$ along the axial direction is very small (especially near the trailing edge (
$z=L$)), i.e. 
$A^{'}(z) \approx 0$, we can get that
\begin{equation} w = 2{\rm \pi} f A(z) \left(1 - \frac{U}{f \lambda} \right) \cos(2{\rm \pi} ft - 2{\rm \pi} z / \lambda) = W \left(1 - \frac{U}{f \lambda} \right). \end{equation}
Actually, just as indicated by Lighthill (Reference Lighthill1971), we have 
$w=W(V-U) / V$, where 
$V = f \lambda$ is the wave speed.
Note that the time-averaged value within one cycle of the integral item in (A6) is zero and 
$A(z) = A_t$ at 
$z=L$. Hence, (A6) becomes
\begin{equation} T \sim m \left(\,f A_t \right)^2 \left(1- \left(\frac{U}{f \lambda}\right)^2 \right). \end{equation} For the flapping plate in our study, there is no obvious undulation during swimming (see figure 1). In other words, 
$\lambda$ is large in our flapping cases and the item 
$(U/f \lambda )$ may be small enough to be neglected. Hence, (A8) is reduced to 
$T \sim (\,f A_t)^2$.
We further estimate the power consumed by the plate. The acceleration of fluid adjacent to the plate is 
$a ={\textrm {D} w}/{\textrm {D} t}$. Considering 
$A^{'}(z) \approx 0$, we have
\begin{equation} a ={-}4 {\rm \pi}^2 f^2 A(z) \left( 1 - \frac{U}{f \lambda} \right)^2 \sin(2{\rm \pi} ft - 2{\rm \pi} z / \lambda). \end{equation}
Similarly, 
$A(L) = A_t$. Neglecting the item 
$(U/f \lambda )$, we have 
$a \sim f^2 A_t$. The mass of fluid pushed by the plate is scaled as 
$\rho L A_t$. The plate does work at a rate equal to its lateral velocity 
$W \sim f A_t$. Hence, the power expenditure is scaled as 
$P \sim (\rho L A_t) \cdot (\,f^2 A_t) \cdot (\,f A_t)$, i.e. 
$P \sim (\,f A_t)^3$.
Appendix B. Drag scaling for the inclined plates
One may notice that the fibre in the study of Alben et al. (Reference Alben, Shelley and Zhang2002) is normal to the incoming flow, i.e. the inclination angle 
$\theta$ in figure 9(a) is 
$90^{\circ }$. For a fixed 
$\theta$, i.e. in Alben's model, increased speed is associated with reduced projected frontal area. People may have an illusion that in the self-propulsive cases, the increased speed is associated with increased frontal area, which is opposite to Alben's cases. However, note that in our study, 
$\theta$ (i.e. pitching amplitude) is variable. Hence, it is unreasonable to directly compare the relationship between speed and projected frontal area. Nevertheless, we would like to illustrate why Alben's drag law can be applied in our study.
Figure 9. (a) Schematic diagram of the inclined plate in oncoming flow, where 
$\theta$ is the inclination angle of the plate and 
$U$ is the oncoming flow velocity. (b–d) The dimensionless drag 
$D^*$ of inclined plates as functions of Reynolds number (
$Re_U = U L/ \nu \sim U$) with different 
$\theta$ and bending stiffness 
$K^*$: (b) 
$\theta =20^{\circ }$, 
$K^*=0.01$, (c) 
$\theta =30^{\circ }$, 
$K^*=0.01$ and (d) 
$\theta =40^{\circ }$, 
$K^*=0.005$.
First, we try to test whether the drag scaling of Alben et al. (Reference Alben, Shelley and Zhang2002) can be applied to the cases with different fixed inclination angle 
$\theta$ (see figure 9
$(a)$), since our flexible plates are more like inclined plates (see figure 1). Here many cases for different 
$\theta$ and bending rigidity 
$K^*$ are simulated. The computational domain, grid spacing and time step are identical to these in § 3. The oncoming flow velocity is U and the Reynolds number is 
$Re_U = U L/ \nu$. The horizontal force experienced by the plate is the drag 
$D$. Note that, under these circumstances, the plate is completely passive and 
$K^*$ should be relatively small to ensure that the bending deformation of the plate is relatively large (Alben et al. Reference Alben, Shelley and Zhang2002; Zhu Reference Zhu2008), so that the plate shape is more streamlined.
The dimensionless drag 
$D^*$ as a function of 
$Re_U$ for cases of different 
$\theta$ and 
$K^*$ are presented in figure 9(b–d). It is seen that the growing forms of drag are very close to that of Alben et al. (Reference Alben, Shelley and Zhang2002) (i.e. 
$D \sim U^{4/3}$). In other words, even when the plate is inclined, Alben's law is still valid. Therefore, if a flexible plate has a streamlined curved shape in an oncoming flow due to the fluid–structure interaction (FSI), 
$U^{4/3}$ drag law seems applicable. This may be regarded as a simple generalization of Alben's drag law.
Besides, we can change our views on the pitching and heaving cases in new inertial and non-inertial frames, respectively, to see the analogies between our cases and those of Alben et al. (Reference Alben, Shelley and Zhang2002).
For a pitching plate, because the cruising speed 
$U$ is almost a constant, we can view the locomotion in an inertial frame (see figure 10a) which is moving with 
$U$. In the new frame, at any time, the oncoming flow with 
$U$ is passing through an inclined flexible plate. It is noticed in the inertial frame, except the head of the plate experiences a rotational moment, the other part of the plate is free to move passively and form a streamlined shape due to the FSI. According to this viewpoint, Alben's law can be applied.
Figure 10. The pitching (a) and heaving (b) plates in an inertial frame. (c) The heaving plate in a non-inertial frame moving with the head of the plate.
For a heaving plate, suppose the cruising speed is 
$U$. In the inertial frame which is moving with 
$U$, the plate is heaving up and down (see figure 10b). We can set up a non-inertial frame moving with the head of the plate. Suppose at time 
$t$, the plate moves upward with a lateral component velocity 
$v$. In this new frame, we can see that the inclined oncoming flow with velocity magnitude 
$\sqrt {U^2+v^2}$ is passing through the flexible plate (see figure 10c). It is easy to imagine that at different times during a period, the oncoming flow from different angles is passing through the flexible plate. In the non-inertial frame, except the head of the plate experiences an inertia force, the other part of the plate is free to move passively and form a streamlined shape due to the FSI. According to this viewpoint, Alben's law should also be valid in this situation.
