Power synchronisations determine the hovering flight efficiency of passively pitching flapping wings

Abstract

The power exchange between fluid and structure plays a significant role in the force production and flight efficiency of flapping wings in insects and artificial flyers. This work numerically investigates the performance of flapping wings by using a high-fidelity fluid–structure interaction solver. Simulations are conducted by varying the hinge flexibility (measured by the Cauchy number, $Ch$, i.e. the ratio between aerodynamic and torsional elastic forces) and the wing shape (quantified by the radius of the first moment of area, $\bar {r}_1$). Results show that the lift production is optimal at $0.05 < Ch \leq 0.2$ and larger $\bar {r}_1$ where the minimum angle of attack is around $45^\circ$ at midstroke. The power economy is maximised for wings with lower $\bar {r}_1$ near $Ch=0.2$. Power analysis indicates that the optimal performance measured by the power economy is obtained for those cases where two important power synchronisations occur: anti-synchronisation of the pitching elastic power and the pitching aerodynamic and inertial powers and nearly in-phase synchronisation of the flapping aerodynamic power and the total input power of the system. While analysis of the kinematics for the different wing shapes and hinge stiffnesses reveals that the optimal performance occurs when the timing of pitch and stroke reversals are matched, thus supporting the effective transfer of input power from flapping to passive pitching and into the fluid. These results suggest that careful optimisation between wing shapes and hinge properties can allow insects and robots to exploit the passive dynamics to improve flight performance.

1. Introduction

The aerodynamics of flapping wings inspired by insects and birds has been studied extensively due to its importance in unsteady aerodynamic fundamental and micro-aerial vehicle (MAV) designs (Wang 2005; Shyy et al. 2010; Deng et al. 2013; Shyy et al. 2016; Eldredge & Jones 2019). One of the typical features of insect wings is flapping forward and backward, implementing pitch reversal to maintain a positive angle of attack throughout both half-strokes to maximise the lift production (Dickinson, Lehmann & Sane 1999). The lift production can be enhanced by employing a variety of mechanisms, as can be found in Shyy et al. (2008, 2010).

A few mechanisms of force production for biolocomotion work through a variety of synchronisations (see e.g. Somps & Luttges 1985; Kang & Shyy 2013; Van Buren, Floryan & Smits 2019; Wang et al. 2020; De & Sarkar 2021; Cai et al. 2022; Liu et al. 2022a). For example, the wake capture mechanism works by synchronising the translational reversal of the wing and the vortex wake created during the previous stroke, so that the effective flow velocity increases, resulting in an additional aerodynamic force peak (Dickinson et al. 1999; Birch & Dickinson 2003; Shyy et al. 2010; Bomphrey et al. 2017). Rapid pitching rotation mechanism synchronises the pitch and flapping stroke (translational) motions to form an advanced rotation, gaining additional lift in a way similar to the Magnus effect (Dickinson et al. 1999; Sun & Tang 2002; Shyy et al. 2010). For propulsive-type flapping foils, the motion is well synchronised to avoid uncontrolled flow separation and vortex shedding, and thus, to save the energy expended (Triantafyllou et al. 2002). Dragonflies often synchronise their forewings and hindwings, maintaining a specific phase relationship between the wings to enhance the lift generation (Somps & Luttges 1985). Detailed studies show that the anti-phase wing motion generates uniform forces with nearly minimal power, which is commonly observed in steady hovering. On the other hand, the in-phase motion generates higher lift, providing an additional force to accelerate, which is normally observed during takeoffs (Wang & Russell 2007; Hu & Deng 2014). In addition, the lift enhancement is obtained via the downwash and the leading-edge–vortex interactions (Hu & Deng 2014). The synchronisation between the foil motion and the vortex wake generated by a leading body can enhance the thrust of the foil to save energy, and to provide sufficient lift for the foil to maintain its location in the wake (see e.g. Beal 2003; Taguchi & Liao 2011; Tian et al. 2011; Stewart et al. 2016).

As discussed above, the synchronisations between motions (pitching and heaving), multiple wings (fore and hind wings), motion–self-generated vortex and motion–leading-body-generated vortex have been extensively studied. However, the synchronisation due to the elasticity of wing material has not been well explored. Actually, wing flexibility has a significant influence on force production and flight efficiency (see e.g. Shyy et al. 2010; Dai, Luo & Doyle 2012; Tian et al. 2013). Increasing evidence gained through direct kinematics measurements (Bergou, Xu & Wang 2007; Gehrke et al. 2022; Mathai et al. 2022), numerical (Bergou et al. 2007; Dai et al. 2012; Tian et al. 2013; Chen et al. 2016; Kolomenskiy et al. 2019; Cai et al. 2022) and analytical modelling (Ennos 1988; Whitney & Wood 2010) have shown that the wing flexibility and passive deformation could significantly enhance the force production and flight performance of flapping wings. Specifically, the flexible wing with medium flexibility enjoys a higher lift force (see e.g. Tian et al. 2013; Shahzad et al. 2018a,b; Gehrke et al. 2022; Mathai et al. 2022). Mathai et al. (2022) and Gehrke et al. (2022) particularly reported that camber is enhanced by the hydrodynamic forces during the flapping motion. Kang & Shyy (2013) found that the optimal lift is obtained when the wing deformation synchronises with the prescribed translational motion.

One of the important consequences of wing flexibility is passive pitching, which has received significant interest over the last decade (Bergou et al. 2007; Tian et al. 2013; Liu et al. 2022a). In MAV designs, purely passive pitching could reduce mechanical complexity and the system mass required in miniature robotic systems as it removes the need to actuate the wing in the pitching direction (Whitney & Wood 2010). To facilitate this application, passively pitching flapping wings have been implemented in several experimental models and MAV designs (Lentink, Jongerius & Bradshaw 2009; Keennon, Klingebiel & Won 2012; Ma et al. 2013; Farrell Helbling & Wood 2018). Passive pitching during a flapping stroke is a consequence of flexibility and mediated by the interplay of elastic restoring, wing inertial and aerodynamic forces, associated with energy transfer between the fluid and structural kinetic and elastic energies (Peng & Milano 2013; Tian et al. 2013; Ishihara & Horie 2016; Kolomenskiy et al. 2019; Cai et al. 2022), and thus, it is governed by many factors, including the architectural and material properties of the wings. The surface density, flexural stiffness as well as the vein distribution of real insect wings vary in the chordwise and spanwise directions (Combes & Daniel 2001; Shahzad et al. 2018b).

To simplify the complicated variation, the torsional flexibility can be lumped together and be modelled by an elastic spring (see e.g. Bergou et al. 2007; Ishihara et al. 2009; Zhang, Liu & Lu 2010; Lei & Li 2020). To understand the passive pitching mechanism of flapping wings, various experimental and numerical studies have been conducted (see e.g. Ennos 1988; Bergou et al. 2007; Ishihara et al. 2009; Eldredge, Toomey & Medina 2010; Spagnolie et al. 2010; Whitney & Wood 2010; Zhang et al. 2010; Ishihara, Horie & Niho 2014; Beatus & Cohen 2015; Chen et al. 2016; Wang, Goosen & van Keulen 2017; Bluman, Sridhar & Kang 2018; Kolomenskiy et al. 2019; Wu, Nowak & Breuer 2019; Lei & Li 2020; Mazharmanesh et al. 2021. Specifically, Chen et al. (2016) performed experiments on an insect-scale passively pitching robotic flapper and compared the results with a quasi-steady dynamic model and a computational fluid dynamic solver incorporating fluid–structure interaction (FSI). They showed that the wing kinematics and flapping efficiency depend on the hinge stiffness and found that stiffer wing hinges achieve favourable pitching kinematics leading to larger mean lift forces. Lei & Li (2020) numerically investigated the effects of different flapping trajectories on the wing's passive pitching dynamics for a fruit fly wing, finding that the optimal lift and lift-to-power ratio are achieved with medium flexural stiffness (i.e. with Cauchy number of approximately 0.3). A special note is given to Bergou et al. (2007), who numerically analysed the aerodynamic pitching power expenditures in four different insect species with passive pitching kinematics. They calculated the pitching power about the torsion axis due to aerodynamic and wing inertial forces and found that the net pitching rotational power is negative, suggesting the feasibility of passive wing pitching without the additional rotational power input from the muscles. The time and rate of the elastic energy released during the supination of a flexible wing can significantly affect its performance. For example, there is a delayed effective pitching motion (related to the active pitching component) for a lower-mass-ratio wing compared with a higher-mass-ratio one, resulting in a higher power economy (Tian et al. 2013). Therefore, the synchronisation between fluid-mechanic, structural kinetic and elastic energies is very important in determining the aerodynamics and efficiency of flexible flapping wings and needs to be carefully studied.

Apart from the flexibility, the wing planform (i.e. shape and aspect ratio) can also significantly affect the passive pitching and thus energy/power synchronisation and flight efficiency (Stanford et al. 2012). In the MAV designs, it is technically easier to change the wing planform compared with implementing changes in the wing kinematics (Ansari, Knowles & Zbikowski 2008a,b). Therefore, several studies have been conducted to seek the optimal wing shape for various sizes and Reynolds number scales (Luo & Sun 2005; Ansari et al. 2008b; Young et al. 2009; Stanford et al. 2012; Li & Dong 2016; Shahzad et al. 2016, 2018a; Bhat et al. 2019b; Wang & Tian 2020; Ji et al. 2022; Wang, Tian & Liu 2022). Specifically, Ansari et al. (2008b) numerically studied the effects of various synthetic wing planforms with fully prescribed kinematics on the aerodynamic performance. They found that wing planforms having larger area outboard produce the highest lift. The recent work of Bhat & Thompson (2022) showed that increasing the leading-edge curvature can further enhance the lift. Shahzad et al. (2016) numerically investigated the effects of wing shape on the aerodynamic performance in hover. The wing shapes were defined by the radius of the first moment of the wing area $\bar {r}_{1}$ and aspect ratio ($AR$) while the power economy $PE$ was defined as the ratio of the mean lift coefficient to the mean aerodynamic power coefficient. They found that maximum $PE$ was achieved at $AR=2.96$. Although the maximum lift was observed at high $\bar {r}_{1}$ (i.e. large area outboard of the wing) and high-$AR$ wings, they recommended low $\bar {r}_{1}$ (i.e. large area inboard of the wing) and high $AR$ to maximise $PE$ for a given lift at the Reynolds numbers of insects. Similar conclusions were obtained for flexible flapping wings (Shahzad et al. 2018a,b). Despite the above studies, it is unclear how the planform would affect power synchronisation.

This work numerically investigates the performance of flexible flapping wings, with a focus on power synchronisation. The wings undergo prescribed flapping (stroke) motion and passive pitching motion, of which the latter is determined by the torsional flexibility. Following previous studies (Spagnolie et al. 2010; Zhang et al. 2010; Lei & Li 2020), the torsional flexibility is lumped together and modelled by an elastic torsional spring at the wing root, with the wing itself being modelled as a rigid plate. The FSI system is solved by an in-house FSI solver based on an immersed-boundary–lattice Boltzmann method. Two parameters are particularly considered. The first parameter is the Cauchy number ($Ch$), which is defined by the ratio of aerodynamic forces acting on the wing and the elastic torsional spring force, and is varied from $0.05$ to $0.6$ covering low-, medium- and high-flexible cases. The other parameter is the radius of the first moment of the wing area normalised with the wingspan ($\bar {r}_1$), which is used to describe the wing shape and is varied from 0.39 to 0.63. The forces and the power economy are studied, and power synchronisation is discussed. Compared with our previous work on flexible flapping wings (see e.g. Tian et al. 2013; Shahzad et al. 2018a,b), this work particularly focuses on how power synchronisations determine the hovering flight efficiency of passively pitching flapping wings. Here, hovering flight is selected as it is considered a vital flight profile in artificial and natural flyers alike.

This paper is organised as follows. The physical and mathematical models are described in § 2 with more details of the derivation provided in Appendix A. The numerical method is given in § 3 with validation presented in Appendix B. The numerical results are presented and discussed in § 4, and concluding remarks are provided in § 5.

2. Model description

As the flight speed of insects is low (generally less than 20 m s$^{-1}$) and the flapping frequency ranges from $10$ to 1000 Hz, the Mach number is much lower than $0.3$ and the flapping period is much larger than the time of sound propagation over the characteristic length (Taylor, Nudds & Thomas 2003; Wang 2005; Landau & Lifshitz 2013). Therefore, the flow around the flapping wing is considered incompressible and is governed by

(2.1)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{u}}= 0, \end{gather}$$

(2.2)$$\begin{gather}\rho\frac{\partial{\boldsymbol{u}}}{\partial{t}}+\rho\boldsymbol{u} \boldsymbol{\cdot}{\boldsymbol{\nabla}{\boldsymbol{u}}} ={-}\boldsymbol{\nabla}{p}+\mu \nabla^2{\boldsymbol{u}}+\boldsymbol{f}, \end{gather}$$

where $\boldsymbol {u}$ is the fluid velocity, $\rho$ is the fluid density, $p$ is the pressure, $\mu$ is the dynamic viscosity and $\boldsymbol {f}$ is the body force density. A wing of mean chord $c$ and span $L$ is modelled to undergo two-angle flapping kinematics, where the two angles of concern are the flapping (stroke) position angle $\phi$ and the pitching angle $\theta$, as shown in figure 1. The deviation angle outside the stroke plane is taken as zero as a simplified case of normal hovering (Ellington & Lighthill 1984). In such a case, the straight leading edge of the wing remains parallel to the horizontal plane throughout the flapping motion. The flapping axis is situated along the $z$-axis, passing through the wing root, and the pitching axis is aligned with the leading edge of the wing. Such a wing model is used because Ansari et al. (2008b) showed that wings with straight leading edges produce higher lift and that the pitching axis located within 0–0.25$c$ from the leading edge provides an optimised compound performance for the wings similar to those in the current study. Thus, the leading edge is used as the pitching axis for simplicity.

Figure 1. The schematic shows the two-angle flapping kinematics of an insect wing, with the stroke angle ($\phi$) and the pitching angle ($\theta$). Here, $\dot {\phi }$ and $\dot {\theta }$ are the stroke and pitching angular velocities, respectively; $c$ is the mean wing chord and $L$ is the wingspan; $x^{\prime }$, $y^{\prime }$ and $z^{\prime }$ are the coordinate axes in the laboratory coordinate system, whereas $x$, $y$ and $z$ are the axes in the body-fixed coordinate system.

The wing flapping motion is prescribed by the stroke angle according to

(2.3)\begin{equation} \phi={-}\frac{\phi_{A}}{2}\cos (2 {\rm \pi}f t), \end{equation}

where $\phi _{A}$ is the peak-to-peak flapping amplitude, and $f$ is the flapping frequency. The stroke amplitude is maintained to be $\phi _A=140^{\circ }$, which corresponds to that of fruit flies (Altshuler et al. 2005; Berman & Wang 2007). The flapping motion profile used here is similar to that used in several previous studies (Li, Dong & Zhao 2018; Lei & Li 2020) as it nominally matches the kinematics of real insects (Li et al. 2018).

Initially, the wing surface is vertically oriented. The passive pitching dynamics is modelled using a torsional spring at the wing hinge, which provides the restoring torque to the wing and returns it to its neutral position (vertical orientation). The passive pitching angle $\theta$ is determined by solving the dynamical equation for the pitching motion of the wing, as described in Kolomenskiy et al. (2019),

(2.4)\begin{equation} \underbrace{-J_{yy}\ddot{\theta}+J_{zy} \ddot{\phi} \cos \theta + \tfrac{1}{2} J_{yy} \dot{\phi}^2 \sin 2\theta}_{T_{iner}^{pitch}}+T_{aero}^{pitch} \underbrace{-K_{s}(\theta - \theta_0)}_{T_{spri}^{pitch}} -C\dot{\theta} =T_{input}^{pitch},\end{equation}

where $J_{yy}$ is the moment of inertia around the pitching axis $y$, $J_{zy}$ is the moment of inertia around the pitching axis when the wing is rotated around the $z$-axis (i.e. flapping axis), $\dot {\theta }$ and $\ddot {\theta }$ are the pitching angular velocity and acceleration, respectively, $\dot {\phi }$ and $\ddot {\phi }$ are the flapping angular velocity and acceleration, respectively, $C$ is the damping coefficient of the spring and set as zero here, $K_s$ is the torsional stiffness of the spring, $\theta _0$ is the rest pitching angle, $T_{iner}^{pitch}$ is the pitching inertial torque, $T_{aero}^{pitch}$ is the aerodynamic pitching moment on the wing, $T_{spri}^{pitch}$ is the torque of the spring and $T_{input}^{pitch}$ is the input torque to actuate the wing pitch (note, $T_{input}^{pitch}=0$ for purely passive pitching). Please note that the angles are defined in the body-fixed coordinate system. Then, the pitching inertial power ($P_{iner}^{pitch}$), the pitching aerodynamic power ($P_{aero}^{pitch}$), the pitching elastic power ($P_{spri}^{pitch}$) and the total input power for pitching ($P_{input}^{pitch}$) can be calculated, respectively, as

(2.5a–d)\begin{align} P_{iner}^{pitch} = T_{iner}^{pitch} \dot{\theta},\quad P_{aero}^{pitch} = T_{aero}^{pitch} \dot{\theta},\quad P_{spri}^{pitch} = T_{spri}^{pitch} \dot{\theta},\quad \mathrm{and} \quad P_{input}^{pitch} = T_{input}^{pitch} \dot{\theta}. \end{align}

The dimensional torsional spring stiffness is defined using a non-dimensional Cauchy number ($Ch$)

(2.6)\begin{equation} Ch=\frac{\rho \phi_A^{2}\, f^{2} c^{3} L^{2}}{K_s}, \end{equation}

which provides a relative measure of the aerodynamic forces acting on the wing and the elastic torsional spring force at the wing hinge. Note that the mean wing chord $c$ and wing length $L$ are maintained as constant across all wing planforms in this study. Thus, the parameters $Ch$ and $\bar {r}_1$, by definition, are independent of each other in this study. For all wing planforms, $Ch$ is systematically varied from $0.05$ to $0.6$ by changing the stiffness $K_s$ of the torsional spring. This range of $Ch$ includes values for realistic passive pitching flapping-wing kinematics (Ishihara et al. 2009; Lei & Li 2020) to well beyond those tested in previous studies, covering low-, medium- and high-flexible cases. The other non-dimensional parameters include the Reynolds number $Re$ and mass ratio $M$, which are respectively given by

(2.7a,b)\begin{equation} Re=\frac{\rho U L}{\mu}=300,\quad M=\frac{m_s}{\rho L}=\frac{m_s}{\rho c}\frac{c}{L}=\frac{M_c}{AR} = 0.338, \end{equation}

where $U=2 f \phi _A L$ is the mean wingtip velocity, $m_s$ is the surface density of the wing material, $M_c$ is the mass ratio based on the mean chord and $AR$ is the aspect ratio of the wing. The value of $M$ dictates the relative effects of the inertial force vs the aerodynamic force (Tian et al. 2013). Here, $M_c$ is maintained to be $1.0$ as this value is close to the mass ratios found in a variety of insect species, such as fruit flies and honeybees (Lei & Li 2020).

To quantify the flapping inertial and aerodynamic power expenditures, the torque used to actuate the wing flapping motion in (2.3) is given (see the derivation in Appendix A) as

(2.8)\begin{equation} \underbrace{(J_{z z}+J_{y y} \sin ^{2} \theta ) \ddot{\phi} + J_{y y} \sin (2 \theta) \dot{\theta} \dot{\phi} - J_{y z}(\ddot{\theta} \cos \theta-\dot{\theta}^{2} \sin \theta)}_{T_{iner}^{flap}} - T_{aero}^{flap} = T_{input}^{flap},\end{equation}

where $J_{zz}$ is the moment of inertia around the $z\hbox{-}$axis, $J_{yz}$ is the moment of inertia around the $z\hbox{-}$axis when the wing is rotated around the pitching axis $y$ (note, $J_{yz}=J_{zy}$), $T_{iner}^{flap}$ is the flapping inertial torque, $T_{aero}^{flap}$ is the flapping aerodynamic torque and $T_{input}^{flap}$ is the total input torque to actuate the wing flapping motion in (2.3). Then, the flapping inertial power ($P_{iner}^{flap}$), the flapping aerodynamic power ($P_{aero}^{flap}$) and the total input power for flapping ($P_{input}^{flap}$) can be calculated, respectively, as

(2.9a–c)\begin{equation} P_{iner}^{flap} = T_{iner}^{flap} \dot{\phi}, \quad P_{aero}^{flap} ={-}T_{aero}^{flap} \dot{\phi}, \quad \mathrm{and} \quad P_{input}^{flap} = T_{input}^{flap} \dot{\phi}. \end{equation}

The combination of (2.4), (2.5a–d), (2.8) and (2.9a–c) can provide a measure of the total power ($=P_{input}^{flap} + P_{input}^{pitch}$) required for flapping (stroke) and pitching as well as the contributions of each parameter to the total power. Here, the total power consumption for the passively pitching wing is only the total input power for flapping $P_{input}^{flap}$, which has no net power consumption for pitching ($P_{input}^{pitch}=0$). The positive power indicates the work done by the wing on the fluid. For example, a positive flapping aerodynamic power means the work done by the wing on the fluid to overcome the drag.

The drag, lift and aerodynamic power coefficients are defined, respectively, as

(2.10a–c)\begin{equation} C_{D}=\frac{2 F_{x}}{\rho U^{2} A}\frac{\dot{\phi}}{|\dot{\phi|}},\quad C_{L}=\frac{2 F_{z}}{\rho U^{2} A},\quad \mathrm{and} \quad C_{P}=\frac{-2 P_{aero}}{\rho U^{3} A}, \end{equation}

where $F_x$ and $F_z$ are the forces acting on the wing by the ambient fluid in the $x$ and $z$ directions, respectively, $P_{aero}$ is the aerodynamic power and $A$ is the surface area of the wing ($A=cL$ for a rectangular plate). Note that $C_D$ is defined in such a way as to measure the force coefficient along $x$ in the opposite direction to the flapping motion. Here, $P_{aero}$ consists of the flapping aerodynamic power and the pitching aerodynamic power ($P_{aero}= P_{aero}^{flap} + P_{aero}^{{pitch }}$). An example is provided in Appendix B.1 to demonstrate the calculation of the aerodynamic powers. The cycle-averaged lift, drag and power coefficients are respectively denoted as $\bar {C}_L$, $\bar {C}_D$ and $\bar {C}_P$. The power economy $PE = \bar {C}_L/\bar {C}_P$ is introduced to discuss the hovering efficacy, which measures the lift production per unit aerodynamic power.

Eight different wing shapes are considered to systematically explore the interaction between wing-hinge properties and wing planform, as shown in figure 2. Each wing shape is generated using the beta function ($\beta (p, q)$) by varying the radius of the first moment of area $\bar {r}_{1}$. The equations for the wing-shape generation are adopted from Ellington (1984a)

(2.11a,b)$$\begin{gather} \bar{r}_k^k=\int_{0}^{1} \bar{c} (\bar{r}^k) \,{\rm d} \bar{r}, \quad \bar{r}_2=0.929\left(\bar{r}_{1}\right)^{0.732}, \end{gather}$$

(2.12a,b)$$\begin{gather}p=\bar{r}_{1}\left(\frac{\bar{r}_{1}\left(1-\bar{r}_{1}\right)}{\bar{r}_{2}^{2}-\bar{r}_{1}^{2}}-1\right),\quad q=\left(1-\bar{r}_{1}\right)\left(\frac{\bar{r}_{1}\left(1-\bar{r}_{1}\right)}{\bar{r}_{2}^{2}-\bar{r}_{1}^{2}}-1\right), \end{gather}$$

(2.13a,b)$$\begin{gather}\beta(p, q)=\int_{0}^{1} \bar{r}^{p-1}(1-\bar{r})^{q-1} \,{\rm d} \bar{r}, \quad \bar{c}=\frac{\bar{r}^{p-1}(1-\bar{r})^{q-1}}{\beta(p, q)}, \end{gather}$$

where $\bar {c}$ is the wing chord normalised by $c$ and $\bar {r}$ is the spanwise distance from the wing root normalised by $L$; $\bar {r}_k$ is the non-dimensional radius of the $k$th moment of the wing area. The values of $\bar {r}_1$ chosen for wing-shape generation are varied between $0.39$ and $0.63$, including the range of 0.43 and 0.6 for most insect wings (Ellington 1984a; Bhat et al. 2019a; Wang & Tian 2020). The wing shapes are modelled with a constant planform area $A$, aspect ratio $AR$, wingspan $L$ and the mean chord length $c$ ($c = L/AR$) for all cases to easily cross-compare performance metrics between them. As $\bar {r}_1$ increases, the planform area of the wing shifts further outboard from the wing root. The value $AR=2.96$ is chosen because previous studies on rigid wing shapes with prescribed flapping kinematics (Shahzad et al. 2016) have shown that maximum power economy is achieved at this aspect ratio independent of $Re$ and $\bar {r}_1$. The values of $J_{yy}$, $J_{zz}$ and $J_{yz}$ for different wings are provided in Appendix B.2.

Figure 2. Eight different wing planforms have been obtained by varying $\bar {r}_{1}$: (a) $\bar {r}_1=0.39$; (b) $\bar {r}_1=0.43$;(c) $\bar {r}_1=0.48$; (d) $\bar {r}_1=0.53$; (e) $\bar {r}_1=0.565$; ( f) $\bar {r}_1=0.58$; (g) $\bar {r}_1=0.605$; (h) $\bar {r}_1=0.63$.

The far-field boundary condition is applied at the external boundaries of the computational domain and the no-slip boundary condition is applied on the wing surface. Zero velocity and constant pressure are specified as the initial condition of the flow. The initial stroke position of the wing is at $-\phi _A/2$ with zero velocity, and the initial pitching angle is the rest pitching angle. Note that we have conducted simulations to verify that the results are independent of the initial condition of the wing (see Appendix B.3).

3. Numerical method

The flow over a flapping wing is solved by an immersed boundary-lattice Boltzmann method (IB–LBM) FSI solver based on the three-dimensional nineteen discrete velocity (D3Q19) lattice Boltzmann method (LBM) with a multi-relaxation-time (MRT) model for modelling the fluid dynamics and the feedback immersed-boundary method (IBM) for handling the boundary conditions at the fluid-solid interface. In the LBM, the macroscopic dynamics of the fluid is the result of the statistical behaviour of the fluid particles, which is described by the distribution function $g_i(\boldsymbol {x},t)$ according to (Lallemand & Luo 2000; Luo et al. 2011)

(3.1)\begin{equation} g_i\left(\boldsymbol{x}+\boldsymbol{e}_i{\rm \Delta} t, t+{\rm \Delta} t\right)-g_i \left(\boldsymbol{x},t\right)=\varOmega_i\left(\boldsymbol{x},t\right)+ {\rm \Delta} t G_i, \end{equation}

where $g_i(\boldsymbol {x},t)$ is the distribution function for particles with velocity $\boldsymbol {e}_i$ at position $\boldsymbol {x}$ and time $t$, ${\rm \Delta} t$ is the time increment, $\varOmega _i(\boldsymbol {x},t)$ is the collision operator and $G_i$ is the forcing term accounting for the body force $\boldsymbol {f}$. The D3Q19 model (D'Humieres et al. 2002) is used on a cube lattice. Compared with the single-relaxation-time collision model, the MRT model has been proven to be numerically more stable (Lallemand & Luo 2000; Xu et al. 2018). Therefore, the MRT collision model is adopted here and is given by (Lallemand & Luo 2000)

(3.2)\begin{equation} \varOmega_i={-}(\boldsymbol{M}^{{-}1}\boldsymbol{S}\boldsymbol{M})_{ij} [g_i(\boldsymbol{x},t)-g_i^{eq}(\boldsymbol{x},t)],\end{equation}

where $\boldsymbol {M}$ is a $19\times 19$ transform matrix, and $\boldsymbol {S}=diag(\tau _0,\tau _1,\ldots,\tau _{18})^{-1}$ is a non-negative diagonal relaxation matrix. The determination of $\boldsymbol {S}$ in a three-dimensional model can be found in D'Humieres et al. (2002). The equilibrium distribution function $g_i^{eq}$ is defined as

(3.3)\begin{equation} g^{eq}_i=\rho\omega_i\left[1+\frac{\boldsymbol{e_i} \boldsymbol{\cdot}\boldsymbol{u}}{c^2_s}+ \frac{\boldsymbol{uu}:(\boldsymbol{e}_i\boldsymbol{e}_i-c^2_s\boldsymbol{I})}{2c^4_s}\right], \end{equation}

where $c_s={\rm \Delta} x/(\sqrt {3}{\rm \Delta} t)$ is the speed of sound, ${\rm \Delta} x$ is the lattice spacing, $\boldsymbol {I}$ is the unit tensor and the weighting factors $\omega _i$ are given by $\omega _0=1/3$, $\omega _{1-6}=1/18$ and $\omega _{7-18}=1/36$. The velocity $\boldsymbol {u}$, mass density $\rho$ and pressure $p$ can be obtained according to

(3.4a–c)\begin{equation} \rho=\sum_{i}g_i,\quad p=\rho c_s^2, \quad \mathrm{and} \quad \boldsymbol{u}=\left(\sum_{i}\textbf{e}_i g_i+\frac{1}{2}\, \boldsymbol{f}{\rm \Delta} t\right)\Bigg/\rho, \end{equation}

respectively. The force scheme proposed in Guo, Zheng & Shi (2002) is adopted to determine $G_i$

(3.5)$$\begin{gather} G_i=[\boldsymbol{M}^{{-}1}(\boldsymbol{I}-\boldsymbol{S}/2)\boldsymbol{M}]_{ij}F_i, \end{gather}$$

(3.6)$$\begin{gather}F_i=\left(1-\frac{1}{2\tau}\right)\omega_i\left[\frac{\boldsymbol{e}_i-\boldsymbol{u}}{c^2_s}+ \frac{\boldsymbol{e}\boldsymbol{\cdot}\boldsymbol{u}}{c^4_s} \boldsymbol{e}_i\right]\boldsymbol{\cdot}\boldsymbol{f}, \end{gather}$$

where $\tau$ is the non-dimensional relaxation time.

The feedback IBM (e.g. Goldstein, Handler & Sirovich 1993; Huang & Tian 2019; Huang et al. 2021b, 2022) is adopted to handle the no-slip boundary conditions on the flapping wing. In this method, the body force $\boldsymbol {f}$ is added in the Navier–Stokes equations to mimic a boundary condition according to

(3.7)$$\begin{gather} \boldsymbol{f}(\boldsymbol{x},t)={-}\int \boldsymbol{F}_{ib}(s,t) \delta (\boldsymbol{x}-\boldsymbol{X}(s,t))\,{\rm d} A, \end{gather}$$

(3.8)$$\begin{gather}\boldsymbol{F}_{ib}(s,t)=\kappa \rho(\boldsymbol{x},t) (\boldsymbol{U}_{ib}(s,t)-\boldsymbol{U}(s,t)), \end{gather}$$

(3.9)$$\begin{gather}\boldsymbol{U}_{ib}(s,t)=\int \boldsymbol{u}(x,t)\delta (\boldsymbol{x}-\boldsymbol{X}(s,t))\,{\rm d}\kern0.7pt\boldsymbol{x}, \end{gather}$$

where $\boldsymbol {F}_{ib}(s,t)$ is the Lagrangian force density, $dA$ is the element surface area of the immersed boundary, $\delta (\boldsymbol {x}-\boldsymbol {X}(s,t))$ is Dirac's delta function, $\boldsymbol {x}=(x^{\prime },y^{\prime },z^{\prime })$ is the coordinate of the fluid lattice nodes, $\boldsymbol {X}$ is the coordinate of the structure (i.e. the flapping wing here), $\kappa$ is the feedback coefficient and is set to $\kappa =2$ m s$^{-1}$ in the LBM simulations (Huang et al. 2021b, 2022), $\boldsymbol {U}_{ib}(s,t)$ is the immersed-boundary velocity obtained by an interpolation at the immersed boundary and $\boldsymbol {U}(s,t)$ is the velocity of the wing. The 4-point discrete delta function $\delta _h(\boldsymbol {x}-\boldsymbol {X}(s,t))$ is used to approximate the Dirac delta function (Peskin 2002)

(3.10)$$\begin{gather} \delta_h(\boldsymbol{x}-\boldsymbol{X}(s,t))= \frac{1}{{\rm \Delta} x^{\prime} {\rm \Delta} y^{\prime} {\rm \Delta} z^{\prime}}\zeta \left(\frac{x^{\prime}-X(s,t)}{{\rm \Delta} x^{\prime}}\right) \zeta\left(\frac{y^{\prime}-Y(s,t)}{{\rm \Delta} y^{\prime}}\right) \zeta\left(\frac{z^{\prime}-Z(s,t)}{{\rm \Delta} z^{\prime}}\right), \end{gather}$$

(3.11)$$\begin{gather}\zeta({r})=\left\{\begin{array}{ll} \dfrac{1}{8}\left(3-2|r|+\sqrt{1+4|r|-4r^2}\right), & 0\leq |r|\leq 1, \\ \dfrac{1}{8}\left(5-2|r|-\sqrt{-7+12|r|-4r^2}\right), & 1\leq |r|\leq 2, \\ 0, & |r|>2. \end{array}\right. \end{gather}$$

The numerical method used here has been extensively validated (e.g. Xu et al. 2018; Ma et al. 2020; Huang 2021; Huang et al. 2021a,b, 2022; Liu, Tian & Feng 2022b) for external and internal flows. The velocity error, the spurious flow penetration and their consequences in external and internal flows of the numerical method have been discussed in Huang et al. (2022). Here, the simulations of a rectangular flapping wing and a fruit fly flapping wing are used to further validate the computations for three-dimensional flapping-wing cases with grid and time-step-independence tests (see Appendix B). To reduce the computational effort, a multi-block LBM (Yu, Mei & Shyy 2002; Xu et al. 2018; Liu et al. 2022b) has been implemented to provide a high-resolution Cartesian grid near the solid body, with relatively lower resolutions in the far field. The computational domain has a size of $30c \times 30c \times 30c$. Four blocks of grid are used, with a minimum grid size of $0.04c$ around the wing and the maximum grid size of $0.32c$ in the far field, resulting in a total grid number of $7.25\times 10^6$. The determination of the finest fluid grid size is based on the two validation cases shown in Appendix B. The grid size of the wing is half of the finest fluid grid size (i.e. $0.02c$) which is required by the IBM employed here. The far-field boundary conditions along six sides of the computational domain are set to be the Dirichlet boundary conditions for the velocity and pressure, which produce almost identical results as the Neumann boundary conditions due to the large computational domain used. The in-house solver is parallelised by a hybrid open multi-processing (OpenMP) and open message-passing interface (OpenMPI). Six stroke cycles are simulated to ensure that all force histories (e.g. $C_D$ and $C_L$) and kinematics are periodic, and the time-averaged values are calculated over the last cycle. Validation has been conducted to show that the cycle-to-cycle variations in $C_L$ and $C_D$ are negligible, as shown in Appendix B.4.

4. Results and discussion

4.1. Flight performance on the plane of $\bar {r}_1$ and $Ch$

The flight performance of the wing measured by cycle-averaged lift coefficient $\bar {C}_L$ and power economy $PE$ on the plane of $\bar {r}_1$ and $Ch$ is shown in figure 3, where several interesting results are obtained. Firstly, $\bar {C}_L$ increases with $\bar {r}_1$ when $Ch<0.15$ where the passive pitching is small. This behaviour is similar to the discussion for rigid wings (Ellington 1984b; Shahzad et al. 2016; Bhat et al. 2019b), which is understandable as the area of the wing is shifted towards the wingtip when $\bar {r}_1$ increases (see also figure 2) generating increasing aerodynamic forces. Secondly, when $Ch$ is large (e.g. ${\ge }0.2$), $\bar {C}_L$ first increases and then decreases with $\bar {r}_1$. The decrease of $\bar {C}_L$ is caused by the large passive pitching and the inertia force during stroke reversal, which have a significant negative effect on lift generation. This observation is consistent with the results obtained with wings modelled by flexible plates of low mass ratio ($M_c\le 1$) and aspect ratio ($AR<3.0$) (Shahzad et al. 2018a,b), and the compliant membrane wings (Gehrke et al. 2022; Mathai et al. 2022). Such complex behaviour is because flexibility becomes important at large $Ch$, and the FSI system undergoes complex power synchronisation, which will be further discussed in § 4.2. Thirdly, when $Ch$ increases from $0.05$ to $0.6$, $\bar {C}_L$ first increases and then decreases, with the peak located in the band of $0.05\le Ch \le 0.2$. This observation is consistent with the studies with wings modelled by flexible plates (see e.g. Dai et al. 2012; Tian et al. 2013; Shahzad et al. 2018b). Overall, a high $\bar {r}_1$ (${\geq }0.565$) and a low $Ch$ (${\approx }0.1$) is useful in maximising $\bar {C}_L$ with an optimal value of ${\approx }0.67$. Finally, the behaviours of $PE$ in relation to $Ch$ are similar to $\bar {C}_L$, but its maximum value (${\approx }0.71$) is obtained with $\bar {r}_1\approx 0.45$ and $Ch\approx 0.2$. This $\bar {r}_1$ value is very close to that of hawkmoth wings, which has been reported to be $\bar {r}_1\approx 0.47$ (Willmott & Ellington 1997). However, in relation to $\bar {r}_1$, $PE$ shows negligible changes for $Ch<0.15$. This will be further investigated in § 4.3.

Figure 3. Flight performance on the plane of $\bar {r}_1$ and $Ch$: (a) mean lift coefficient $\bar {C}_L$, and (b) power economy $PE$.

We should point out that, even though the above observations are obtained at $Re=300$, $M_c=1.0$ and a specific wing shape, they can be generalised to other conditions (e.g. different $Re$ and $M_c$), with the optimal values and locations dependent on these conditions. This statement is made based on the results from previous studies. Specifically, the optimal $\bar {C}_L$ is larger, but the optimal $PE$ is smaller for higher $AR$, as reported by Shahzad et al. (2018a), who considered flexible wings at $Re=400$. Moreover, the optimal locations of $\bar {C}_L$ and $PE$ are shifted to larger $\bar {r}_1$ for higher $AR$ (Shahzad et al. 2018a). In addition, insects may fly close to, but not exactly in, the optimal region. For example, the wings of fruit flies are found to be of $\bar {r}_1 \approx 0.53$ (Meng, Liu & Sun 2017) and $Ch \approx 0.27$ (Lei & Li 2020). According to figure 3, $\bar {C}_L=0.43$ might be just sufficient to balance their weight and $PE \approx 0.58$ might be the near-optimal power economy. Overall, the current analysis highlights the strong influence of the coupled effects of wing shape and pitching spring stiffness on wing performance, which will be helpful for designing the flapping-wing flyers and further explain the mechanisms of insect flight.

To summarise this section, two distinct regions of maximum lift production and efficiency are identified. It is found that a high $\bar {r}_1$ (${\geq }0.565$) and a low $Ch$ (${\approx }0.1$) is useful in maximising $\bar {C}_L$, and the power economy is maximised for wings with lower $\bar {r}_1$ (${\approx }0.45$) near $Ch=0.2$.

4.2. Effects of the torsional spring stiffness

To discuss the mechanisms of high lift generation and flight performance, we present the results obtained by varying $Ch$ for three values of $\bar {r}_1$ (0.39, 0.48 and 0.58). The forces, power economy, passive pitching and power components are particularly analysed.

We first discuss $\bar {C}_L$, $\bar {C}_D$ and $PE$, as shown in figure 4. The trends of $\bar {C}_L$ and $PE$ are similar when $Ch$ increases from 0.05 to 0.6, as introduced in § 4.1. The peak values of $\bar {C}_L$ are approximately 0.36, 0.52 and 0.65 for $\bar {r}_1=0.39$, 0.48 and 0.58, respectively (see figure 4a). For all three $\bar {r}_1$ values, the highest value of $\bar {C}_L$ is observed between $Ch=0.1$ and $0.15$, suggesting an optimum range of $Ch$ to achieve the highest possible $\bar {C}_L$. The exact value of the optimum $Ch$ is found to vary within this small range with a change in the wing shape. For all the wings, $\bar {C}_D$ is observed to decrease monotonically with $Ch$, see figure 4(b). This observation is consistent with that of the two-dimensional hovering wing studied, for example, by Yin & Luo (2010). The decrease in $\bar {C}_D$ for a higher $Ch$ is caused by the larger passive pitching motion resulting in a smaller angle of attack. On the other hand, $PE$ increases with $Ch$, reaching its maximum value at $Ch\approx 0.15$ before decreasing with $Ch$, similar to the trend observed in $\bar {C}_L$ (see figure 4c). From $Ch=0.05$ to $0.15$, $PE$ is increased due to the increase of lift and the decrease of aerodynamic power ($\bar {C}_P$) caused by the softening of the spring. For $Ch>0.15$, $PE$ is found to decrease as the rate of the decrease in $\bar {C}_L$ is higher than that in $\bar {C}_P$. The significant and quick decrease of $\bar {C}_L$ is due to the negative lift caused by the delayed pitch for soft springs and the reduced size of the leading-edge vortex, which will be demonstrated in the analysis of wing kinematics and vortex structures.

Figure 4. Variations in (a) $\bar {C}_L$, (b) $\bar {C}_D$ and (c) $PE$ are shown as functions of $Ch$ for three different wing shapes of $\bar {r}_1=0.39$, 0.48 and 0.58.

We then present the flapping and pitching power expenditures from various sources, as shown in figure 5, from which several important observations are obtained. Firstly, the flapping aerodynamic power dominates power consumption with peak values roughly 3 times those of the flapping inertial power. This is consistent with the definition of the mass ratio, $M=M_c/AR\approx 1/3$, which indicates that the scale of the inertial forces is approximately 1/3rd of that of the aerodynamic forces. The exception is the case of $Ch=0.6$, where the spring is very soft, resulting in a much smaller aerodynamic force (thus much smaller aerodynamic power) and larger pitching inertial power compared with that of lower $Ch$ cases. Secondly, the flapping aerodynamic power almost always remains positive due to a low mass ratio considered and the continuous work required to overcome the drag. The exception is the case of a very stiff wing (i.e. $Ch=0.05$) with a negative aerodynamic power just before the end of each half-stroke, similar to that observed with a two-dimensional flapping wing in Tian et al. (2013), and is caused by the deceleration of the stroke motion. According to two-dimensional simulations (Yin & Luo 2010; Tian et al. 2013), the cases with high mass ratios could also generate negative aerodynamic power due to a delayed rotation. Normally, the flapping aerodynamic power reaches its maximum value at around midstroke (e.g. $tf=5.25$ and $tf=5.75$), where the drag and the stroke velocity are maximal. Thirdly, for pitching motion, the elastic power of the spring balances the pitching–inertial and pitching–aerodynamic powers and, thus, no net power consumption is observed. Fourthly, the flapping aerodynamic power decreases monotonically with $Ch$ similar to the trend observed in $C_D$ (figure 4b), while the flapping and pitching inertial powers increase monotonically with $Ch$ because the wing of higher $Ch$ travels over a larger $\theta$ range, resulting in higher $\dot \theta$ and $\ddot \theta$. As a result, the pitching elastic power increases correspondingly to balance the increased pitching inertial and aerodynamic powers.

Figure 5. Time traces of power coefficients ($CP$ values) at $\bar {r}_1=0.48$ for six different spring stiffness values indicated by $Ch=0.05$, $Ch=0.1$, $Ch=0.15$, $Ch=0.2$, $Ch=0.3$ and $Ch=0.6$. $CP_{iner}^{flap}$: flapping inertial power coefficient; $CP_{aero}^{flap}$: flapping aerodynamic power coefficient; $CP_{iner}^{pitch}$: pitching inertial power coefficient; $CP_{aero}^{pitch}$: pitching aerodynamic power coefficient; $CP_{spri}^{pitch}$: pitching elastic power coefficient; $CP_{total}$: total power coefficient. $CP=power/(0.5 \rho U^3 A)$, where $A$ is the surface area of the wing. Here, $tf$ is non-dimensional time, normalised by flapping frequency $f$.

Finally, for $Ch\in [0.1$, $0.2]$, a nearly in-phase synchronisation (with a phase difference ${\rm \Delta} \phi _{flap} \in [18^\circ, 30^\circ ]$ as measured in figure 6a) is observed between the flapping aerodynamic power and the total system input power. Simultaneously, an anti-phase synchronisation (with a mean phase difference ${\rm \Delta} \phi _{pitch} \approx 180^\circ$ as measured in figure 6a) between the pitching elastic power and the pitching aerodynamic and inertial powers is observed. The nearly in-phase synchronisation between the flapping aerodynamic power and the total system input power ensures that the total system input power is mainly used to drive the flapping motion without being wasted in other power components, while the anti-phase synchronisation between the pitching powers ensures that the pitching aerodynamic and inertial powers are fully stored in the elastic pitching power and vice versa. The phase difference between flapping aerodynamic power and the total system input power shows an increasing trend as the increase of $Ch$. The nearly in-phase synchronisation between the flapping aerodynamic power and the total system input power is also reflected by the timing of pitch rotation, which will be further analysed below. As shown in figure 6(b), the maximal absolute value of the peak pitching elastic power occurs at $Ch\in [0.1$, $0.2]$, indicating an efficient storage and release process of the elastic energy in this region. This is consistent with the work by Cai et al. (2022) who reported that the elastic storage minimises the high energetic cost of flapping wings. For $Ch=0.05$, the wing is stiff, and while the synchronisation between flapping and pitching is similar to that required for high lift, the wing's minimal angle of attack is much higher than $45^{\circ }$ where the optimum lift is produced. In addition, the drag is high for $Ch=0.05$ resulting in a low power economy. Wings with $Ch \in [0.1, 0.2]$ produce an optimal lift and a relatively small drag as the wing experiences the minimal angle of attack ${\approx }45^{\circ }$, which will be analysed and discussed in the following parts of this section.

Figure 6. (a) Phase difference (${\rm \Delta} \phi$) of power coefficients, as a function of $Ch$. Here, ${\rm \Delta} \phi _{flap}$: phase difference between the flapping aerodynamic power and the total system input power; ${\rm \Delta} \phi _{pitch}$: mean phase difference between the pitching elastic power and the pitching aerodynamic and inertial powers. (b) The absolute value of the peak pitching elastic power.

To further study the reasons for high $\bar {C}_L$ and $PE$ with various $Ch$, the time histories of $C_L$, $C_D$, the pitching angle $\theta$ and the corresponding angle of attack $\alpha$ are shown in figure 7; the minimal angle of attack $\alpha _{min}$ and the pitch rotation delay $\delta _{\alpha }$ are shown in figure 8. It should be noted that the pitching angle $\theta$ about the pitching axis is measured with reference to the $z$-axis. Hence, the angle of attack $\alpha$ is defined as

(4.1)\begin{equation} \alpha=\begin{cases} 90^{{\circ}}-\theta, & \text{if $\dot{\phi}<0$},\\ 90^{{\circ}}+\theta, & \text{if $\dot{\phi}\geq 0$}. \end{cases} \end{equation}

Various quasi-steady models correlate $C_L$ and $\alpha$ as $C_L=f[\sin (2\alpha )]$ (Dickinson et al. 1999; Sane & Dickinson 2001, 2002; Wang, Goosen & van Keulen 2016). Thus, $C_L$ is expected to be maximum when $\alpha \approx 45^\circ$, which is slightly higher than that for a translated flat-plate wing (Taira & Colonius 2009). On the other hand, $C_D$ is expected to increase with $\alpha$. As the softer springs at higher $Ch$ result in lower $\alpha$ values, the corresponding $C_D$ is also lower, irrespective of $\bar {r}_1$. For a higher $\bar {r}_1$, a larger area is available outboard to support a larger leading-edge vortex (LEV), which results in higher suction magnitudes beneath the LEV. Thus, values of both $C_L$ and $C_D$ are slightly more enhanced in magnitude compared with those at lower $\bar {r}_1$.

Figure 7. Time traces of $\theta$, $\alpha$, ${C}_{L}$ and $C_D$ are shown for passively pitching wings with four different $Ch$ values. The left column shows the results for the wing of $\bar {r}_1=0.39$ and the right column shows the results for the wing of $\bar {r}_1=0.58$. Data for fruit fly kinematics are obtained from Fry, Sayaman & Dickinson (2005). Here, $tf$ is the non-dimensional time. Note that the results are shown for the 6th flapping cycle at which point the forces and kinematics are periodic. Markers are shown at every 10th time step in the calculations. (a) Time traces of $\theta$ for $\bar {r}_1=0.39$. (b) Time traces of $\theta$ for $\bar {r}_1=0.58$. (c) Time traces of $\alpha$ for $\bar {r}_1=0.39$. (d) Time traces of $\alpha$ for $\bar {r}_1=0.58$. (e) Time traces of $C_L$ for $\bar {r}_1=0.39$. ( f) Time traces of $C_L$ for $\bar {r}_1=0.58$. (g) Time traces of $C_D$ for $\bar {r}_1=0.39$. (h) Time traces of $C_D$ for $\bar {r}_1=0.58$.

Figure 8. Contours of $\alpha _{min}$ reached in a flapping stroke are shown in (a) on the map of $\bar {r}_1$ and $Ch$, where the white colour represents $\alpha _{ref}\approx 45^{\circ }$. The contours of pitch-rotation delay ($\delta _{\alpha }$) as a percentage of the flapping period are shown in (b) on the map of $\bar {r}_1$ and $Ch$.

The instantaneous values of $C_L$ are also affected by the synchronisation in the flapping and pitching kinematics. The delay in pitch reversal $\delta _{\alpha }$ is calculated as the phase difference between the flapping stroke reversal and the time when $\alpha$ reaches the neutral position of $90^{\circ }$. From figures 7 and 8, a few interesting results are obtained. Pitch reversal is more delayed with increasing $Ch$, as seen in figures 7(a) and 7(b). The delayed pitch causes the instantaneous $\alpha$ after the flapping stroke reversal to be greater than $90^{\circ }$, which results in negative $C_L$ due to the wing being pushed down by the fluid on the wing pressure side. Additionally, the rapid pitch-down motion observed prior to $t/T=5.2$ results in high circulation around the wing that contributes to a significantly high downward force. Thus, $Ch=0.6$ shows negative values of $C_L$ during a part of a stroke, which is responsible for a lower $\bar {C}_L$ of a complete cycle, see figures 7(e) and 7( f). Close to $t/T=5.2$, the pitch-down motion decelerates, followed by the pitch-up motion, which contributes to a positive lift even at low values of $\alpha$ close to $t/T=5.35$. Irrespective of the wing shape, the pitching amplitude of the wing during the midstroke is observed to increase with $Ch$; consequently, the angle of attack $\alpha$ reached near the midstroke is lower, see figures 7(c) and 7(d). The case $Ch=0.15$ experiences the highest instantaneous $C_L$, as the wing nearly maintains $\alpha \approx 45^{\circ }$ (which is also very close to $\alpha _{min}$, see figures 7(c) and 7(e)), which has been shown in previous studies to be the most effective $\alpha$ for lift generation (e.g. Dickinson 1994; Dickinson et al. 1999). The wing of $\bar {r}_1=0.58$ approaches significantly lower values of $\alpha _{min}$. Combined with a delayed pitch, this behaviour results in large negative instantaneous values of $C_L$ as well as an overall lower $\bar {C}_L$. Typically, fruit fly wings have $\bar {r}_1\approx 0.53$ (Meng et al. 2017) and $Ch\approx 0.27$ (Lei & Li 2020). The time traces of $\theta$ and $\alpha$ for a hovering fruit fly from Fry et al. (2005) are compared in figures 7(b) and 7(d) along with a wing with $\bar {r}_1=0.58$. The variations in $\theta$ and $\alpha$ for a fruit fly wing appear to follow the time traces between the cases $Ch=0.15$ and 0.3, which are not in, but close to, the optimal lift region. For all wing shapes, $C_D$ decreases monotonically with $Ch$, caused by the reduced angle of attack for softer springs, see figures 7(g) and 7(h).

The minimal angle of attack $\alpha _{min}$ is found to decrease continuously with increasing $Ch$, see figure 8(a). This is expected since a softer spring will allow the wing to deform more in pitch from its neutral position. Wings with $Ch$ in the range of $[0.1, 0.2]$ produce a maximum $C_L$ (figure 4a) and experience $\alpha _{min}\approx 45^{\circ }$. The value of $Ch$ required to achieve $\alpha _{min}=45^{\circ }$ reduces with $\bar {r}_1$. This is because, for a higher $\bar {r}_1$, the larger area available outboard to support the LEV results in increased suction pressure, whose resultant force makes the wing tilt more than that for a lower $\bar {r}_1$. Consequently, the region of $\alpha _{min}=45^{\circ }$ in figure 8(a) is observed to decrease linearly on the map of $\bar {r}_1$ and $Ch$. This also explains why the region of high $\bar {C}_L$ in figure 3(a) is shifted to lower $Ch$ at high $\bar {r}_1$ values.

The delay in pitch reversal $\delta _{\alpha }$ is seen to increase with $Ch$ for all the examined wing shapes, see figure 8(b). For stiffer springs with $Ch<0.15$, the pitch rotation is advanced with respect to the flapping stroke (i.e. $\delta _{\alpha }<0$), whereas for softer springs, the pitch rotation is delayed (i.e. $\delta _{\alpha }>0$). As has been previously reported by Dickinson et al. (1999) and Sane & Dickinson (2001), the advanced pitch can be useful in achieving higher lift and the delayed pitch may be detrimental. On the other hand, the advanced pitch also results in a higher $C_D$. This confirms the in-phase synchronisation between the flapping aerodynamic power and the total system input power. Similarly, the anti-phase synchronisation between the elastic pitching power and the pitching aerodynamic and inertial powers is reflected in $\delta _{\alpha }\approx 0^{\circ }$ for optimal $PE$ cases (see the $Ch=0.15$ case in figure 5). The region of $\delta _{\alpha }=0$ represents the exact synchronisation between the flapping and pitching motions. This is the region where the highest power economy can be observed for a given wing shape or stiffness, which explains the variations in $PE$ shown in figure 3(b). Note that the cases with high flexibility, i.e. $Ch>0.45$ and $\bar {r}_1>0.6$ have the wing rotated up to and beyond $\alpha _{min}=0^{\circ }$, causing the lift and drag to be highly unstable during a stroke. Overall, $Ch \in [0.1,0.2]$ may be ideal for simultaneously obtaining $\alpha _{min}=45^{\circ }$ to maximise $\bar {C}_L$ and $\delta _{\alpha }\leq 0$ to maximise $PE$.

To summarise this section, the lift production is optimal at $0.05 < Ch \leq 0.2$ and larger $\bar {r}_1$, where the minimum angle of attack is around $45^\circ$ at the midstroke, while the best configuration for the maximum $PE$ is achieved at $\bar {r}_1=0.45$ and $Ch=0.2$, where the flapping aerodynamic power is approximately three times in magnitude of the flapping inertial power.

The optimal performance measured by the power economy is obtained for those cases where two important power synchronisations occur: anti-synchronisation of the pitching elastic power and the pitching aerodynamic and inertial powers and nearly in-phase synchronisation of the flapping aerodynamic power and the total input power of the system.

4.3. Wing-shape effects

The effects of wing shape on the aerodynamic performance are explored by varying $\bar {r}_1$. The variation in the cycle-averaged lift coefficient ($\bar {C}_L$) as a function of $\bar {r}_1$ using three different torsional springs corresponding to $Ch = 0.05$, $0.15$ and $0.45$, is shown in figure 9(a). For the case of the stiffest spring at the wing hinge (i.e. $Ch=0.05$), $\bar {C}_L$ increases monotonically with $\bar {r}_1$. This is in concurrence with previous observations on actively pitching rigid wing shapes with prescribed kinematics. It has been shown that the larger outboard area of a high-$\bar {r}_1$ wing supports the larger LEV that is responsible for creating higher suction on the wing surface, resulting in higher lift (Shahzad et al. 2018a; Bhat et al. 2019b). However, this observation does not hold true for passively pitching wings with softer torsional springs, see figure 9(a). For wings with softer springs (i.e. $Ch=0.15$ and $Ch=0.45$), $\bar {C}_L$ increases with $\bar {r}_1$ only until a certain $\bar {r}_1$, beyond which it starts decreasing. The value of $\bar {r}_1$ giving a maximum $\bar {C}_L$ changes with $Ch$: $\bar {C}_{Lmax}=0.573$ at $\bar {r}_1=0.565$ for $Ch=0.15$ and $\bar {C}_{Lmax}=0.187$ at $\bar {r}_1=0.48$ for $Ch=0.45$.

Figure 9. Variations in (a) $\bar {C}_{L}$, (b) $\bar {C}_{D}$ and (c) $PE=\bar {C}_{L}/\bar {C}_{P}$ with $\bar {r}_1$ (wing shape) are shown for passively pitching wings with three different $Ch$ values.

Similar trends are also observed in the cycle-averaged drag coefficient ($\bar {C}_D$), as shown in figure 9(b). For the stiffest hinge (i.e. $Ch=0.05$), $\bar {C}_D$ appears to increase monotonically with $\bar {r}_1$ but for less stiff hinges such as $Ch=0.15$ and $0.45$, $\bar {C}_D$ increases with $\bar {r}_1$ to a maximum after which it decreases with $\bar {r}_1$. Previous studies (e.g. Shahzad et al. 2018a) have shown that, for actively pitching rigid wings, $PE$ decreases with an increase in $\bar {r}_1$ due to an increased $\bar {C}_D$. However, for passively pitching wings, the observations are different, as can be seen in figure 9(c). For $Ch=0.05$, the increase in the outboard area with $\bar {r}_1$, together with a reducing angle of attack $\alpha$ approaching $45^{\circ }$, results in a greater increase in $\bar {C}_L$ than in $\bar {C}_D$, causing $PE$ to increase with $\bar {r}_1$ for $\bar {r}_1>0.43$. However, for $Ch=0.45$, both $\bar {C}_L$ and $\bar {C}_D$ are lower, resulting in the highest $PE$ ($0.466$) at $\bar {r}_1= 0.43$. For the medium-soft spring case ($Ch=0.15$) the balance between spring stiffness and aerodynamic contribution results in the highest $PE=0.662$, which remains virtually constant for $0.43\leq \bar {r}_1 \leq 0.565$.

Overall, the variation in $PE$ appears to directly depend on the pitch-rotation delay caused by the power synchronisations, as discussed earlier. For $Ch=0.05$, the pitch rotation is advanced compared with the flapping motion, as can be seen in figure 8(b). At higher $\bar {r}_1$, the rotation is less advanced, which causes an increase in $PE$ seen here in figure 9(c). At $Ch=0.15$, the delay of $\delta _{\alpha }\approx 2.5\,\%$ is mostly unaffected by $\bar {r}_1$. Thus, the value of $PE$ also remains mostly unchanged. However, at a higher $Ch$, the pitch rotation is more delayed. It first decreases with $\bar {r}_1$, followed by an increase. Hence, the resulting $PE$ also shows a variation as shown for the case of $Ch=0.45$ in figure 9(c).

It can be seen that the choice of the wing shape to achieve the best possible hover performance, both in terms of $\bar {C}_L$ and $PE$, is determined by the spring stiffness used for the passively pitching wings. The reasons behind the trends in $\bar {C}_L$ with various $\bar {r}_1$ are further analysed by comparing their time histories of the pitching angle $\theta$ and the corresponding angle of attack $\alpha$, as shown in figure 10.

Figure 10. Time traces of the pitching angle $\theta$, angle of attack $\alpha$ and lift coefficient ${C}_{L}$ are shown for passively pitching wings with three different $\bar {r}_1$ values. The left column shows the results for $Ch=0.05$ and the right column shows the results for $Ch=0.45$. Here, $tf$ is non-dimensional time, normalised by flapping frequency $f$. Note that the results are shown for the 6th cycle at which point the forces and kinematics are periodic. Markers are shown at every 10th time step in the calculations. (a) Time traces of $\theta$ for $Ch=0.05$. (b) Time traces of $\theta$ for $Ch=0.45$. (c) Time traces of $\alpha$ for $Ch=0.05$. (d) Time traces of $\alpha$ for $Ch=0.45$. (e) Time traces of $C_L$ for $Ch=0.05$. ( f) Time traces of $C_L$ for $Ch=0.45$.

The pitching amplitude is observed to increase with $\bar {r}_1$, as can be seen in figures 10(a) and 10(b). The corresponding variations in the angle of attack $\alpha$ are shown in figures 10(c) and 10(d). As $\alpha$ approaches $45^{\circ }$, $C_L$ is also observed to increase, in accordance with the quasi-steady models (Dickinson et al. 1999; Sane & Dickinson 2001, 2002), as shown in figures 10(e) and 10( f). On the other hand, for a softer spring, i.e. $Ch=0.45$, the elastic energy from the spring is low and the pitching amplitudes are high, resulting in $\alpha <45^{\circ }$ during the midstroke (e.g. $tf=5.25$). Hence, the corresponding values of $C_L$ at midstroke are lower than those for $Ch=0.05$. Please note that for most cases, $\alpha$ reaches its minimum value before the midstroke and reduces thereafter even though the stroke velocity is still increasing. This is caused by the inertial force of the wing during stroke reversal (see e.g. Dai et al. 2012; Tian et al. 2013). Due to the large pitching motions experienced by the wing with a softer spring, there are also significant changes in the phase of pitch reversal, identified by the phase at which the wing crosses the $\alpha =90^{\circ }$ threshold. Instantaneous $C_L$ can assume significantly negative values near the midstroke for $\bar {r}_1=0.63$ and $Ch=0.45$ (see figure 10f) which contributes to the low $\bar {C}_L$ produced by the wing. The negative $C_L$ is due to the high negative pressure that develops near the trailing edge of the wing pressure side during stroke reversal (see supplementary movie 1 available at https://doi.org/10.1017/jfm.2023.821 for an animation of the time history of the lift coefficient and the corresponding vortex structures). The animation for $\bar {r}_1=0.63$ and $Ch=0.05$ is also provided for comparison (see supplementary movie 2).

The wing-shape effects can be better understood by observing the flow structures around the wing, as shown in figure 11. Here, the flow over wings of various $\bar {r}_1$ for $Ch=0.05$ and $Ch=0.45$ is shown at the midstroke, i.e. at $tf=5.25$. The LEV, the dominant flow feature shown to be responsible for the high lift produced by flapping wings, is highlighted in all cases. Since the value of $\alpha$ at midstroke is higher with a stiffer spring ($Ch=0.05$) and is closer to 45$^{\circ }$, the size of the LEV over the wings is larger than the LEV that forms over the wing with a softer spring ($Ch=0.45$). The larger LEV creates a higher suction on the wing surface beneath the vortex. Moreover, the LEV is smaller near the root and increases in size in the spanwise direction. Hence, a larger outboard area of the wing with $\bar {r}_1=0.63$ supports a larger LEV, which results in a stronger suction. This explains the increase in $C_L$ with $\bar {r}_1$ for the wing with a stiff hinge (figure 9(a) for $Ch=0.05$) or fully prescribed kinematics (Combes & Daniel 2001; Shahzad et al. 2018a). For a softer spring with $Ch=0.45$, the suction pressure beneath the LEV increases with $\bar {r}_1$. However, beyond $\bar {r}_1=0.48$, the value of $\alpha$ at the midstroke also reduces significantly, causing the LEV to be smaller and weaker. In addition, the wing of $\bar {r}_1=0.63$ and $Ch=0.45$ rapidly pitches up (i.e. the trailing edge swings up) before $tf=5.2$, leading to a strong vortex at the bottom surface near the trailing edge (see figure 11f), and a strong deceleration occurs thereafter until $tf=5.25$ causing a strong inertia force downwards. Therefore, even with a larger outboard area and a positive angle of attack, the $\bar {r}_1=0.63$ wing produces smaller $\bar {C}_L$, as seen in figure 11( f). Overall, the variation in $C_L$ with $\bar {r}_1$ is found to be dependent largely on the torsional spring stiffness.

Figure 11. Vortices are visualised by semi-transparent iso-surfaces of the normalised Q-criterion ($Q^{\ast }=100$) and the wing suction sides are coloured by the normalised pressure ($p^{\ast }=p/(\rho U^2)$) contours at the midstroke $tf=5.25$. (a–c) Represent the wings of $\bar {r}_1=0.39$, 0.48 and 0.63, respectively, with $Ch=0.05$.(d–f) Represent the same wings with $Ch=0.45$.

Power expenditures from various sources of four wing shapes $\bar {r}_1=0.39$, $0.48$, $0.58$ and $0.63$ are further analysed in figure 12. For all four wing shapes, the flapping inertial power traces two nearly sinusoidal curves with the positive part being larger than the negative part, see figure 12. The positive flapping inertial power is used to accelerate the wing at the beginning of each stroke cycle (i.e. from $tf=5$ to $tf=5.25$), and the negative inertial power indicates the work done by the fluid on the wing to decelerate it after the midstroke (i.e. from $tf=5.25$ to $tf=5.5$). The pitching inertial power increases with $\bar {r}_1$, due to the significantly increased inertial contribution from flapping (i.e. $J_{zy}$). As a result, the pitching elastic power increases with $\bar {r}_1$ to balance the increased pitching inertial and pitching aerodynamic power.

Figure 12. Time traces of power coefficients ($CP$ values) at $Ch=0.15$ for four different wing shapes of ${\bar {r}_1=0.39}$, ${\bar {r}_1=0.48}$, $\bar {r}_1=0.58$ and $\bar {r}_1=0.63$. Here, $CP_{iner}^{flap}$: flapping inertial power coefficient; $CP_{aero}^{flap}$: flapping aerodynamic power coefficient; $CP_{iner}^{pitch}$: pitching inertial power coefficient; $CP_{aero}^{pitch}$: pitching aerodynamic power coefficient; $CP_{spri}^{pitch}$: pitching elastic power coefficient; $CP_{total}$: total power coefficient. Power coefficient $CP=power/(0.5 \rho U^3 A)$, where $A$ is the surface area of the wing. Here, $tf$ is non-dimensional time, normalised by flapping frequency $f$.

As $\bar {r}_1$ increases from $0.39$ to $0.58$, the peak value of the flapping aerodynamic power increases by $42.50\,\%$ (i.e. from $1.20$ to $1.71$), followed by a minor decrease ($9.36\,\%$) as $\bar {r}_1$ increases to $0.63$. On the other hand, the flapping inertial power increases monotonically and significantly with $\bar {r}_1$ as the moment of inertia around the flapping axis (i.e. $J_{zz}$) is doubled ($J_{zz}=5.64$ and $11.38$ for $\bar {r}_1=0.39$ and $0.63$, respectively). At $Ch=0.15$, as shown in figure 9, $\bar {C}_D$ and $PE$ are roughly the same for different $\bar {r}_1$, with a $10.15\,\%$ increase of $\bar {C}_L$ for $\bar {r}_1=0.48$ and $\bar {r}_1=0.58$. The slight increase of $\bar {C}_L$ is at the cost of a significant increase in the flapping inertial power ($50.52\,\%$ in the peak value).

5. Conclusions

A numerical investigation of the aerodynamic performance of passively pitching flapping wings for different combinations of wing shapes (defined by the radius of the first moment of area $\bar {r}_1$) and the pitching hinge stiffness (defined by Cauchy number $Ch$) is conducted using an IB–LBM FSI model. The two parameters ($\bar {r}_1$ and $Ch$) of the wing are found to have strong coupled effects on aerodynamic performance in hover. For passively pitching flapping wings, the variations in the power economy ($PE$) with $\bar {r}_1$ are found to be dependent on the hinge stiffness and significantly different from that reported for wings with prescribed kinematics. The variations in the spring stiffness are observed to result in variations in the phase synchronisations between the elastic, aerodynamic, inertial and total power. The optimum wing performance occurs when there is in-phase synchronisation between the flapping aerodynamic and total system input power, simultaneously with the anti-phase synchronisation between the pitching elastic and the pitching aerodynamic and inertial powers where the absolute value of the peak pitching elastic power occurs. Across all wing shapes, the optimum value of $Ch$ to achieve the highest possible cycle-averaged lift coefficient ($\bar {C}_L$) occurs in the range of $0.05 < Ch \leq 0.2$, where the wing's minimum angle of attack approaches $45^{\circ }$. The combination of $\bar {r}_1$ and $Ch$ is also found to affect the phase between the pitch reversal and flapping motions that could either positively or detrimentally affect performance and $PE$. For cases where the power economy was highest, there is also the synchronisation between the pitching and flapping kinematics of the wings resulting in them being in phase during stroke reversal. The best configuration for the maximum $PE$ is achieved at $\bar {r}_1=0.45$ and $Ch=0.2$, where the flapping aerodynamic power is approximately three times in magnitude of the flapping inertial power. Finally, comparing the wing shape and kinematics of insects with the results presented here suggests that passive pitching may allow insects like fruit flies to support their weight while maximising power economy. This study provides insight into passive pitching hovering performance that may be used to optimise the design of flapping-wing flyers.