2.1. Wing model and kinematics

A rectangular wing model is commonly used to study flapping wing kinematics (Shyy et al. Reference Shyy, Aono, Chimakurthi, Trizila, Kang, Cesnik and Liu2010; Diaz-Arriba et al. Reference Diaz-Arriba, Jardin, Gourdain, Pons and David2021; Gehrke & Mulleners Reference Gehrke and Mulleners2021): even though insect wings have more complex shapes, a rectangular planform is a relevant simplified model allowing the accurate reproduction of the main physical phenomena at stake (namely the leading and trailing edge vortex dynamics). Thus, the wing model considered in this study is rectangular with a NACA0012 profile, where $c=0.01$ m is the chord length and $b=0.04$ m is the wingspan.

Two angles are used to describe the kinematics: the stroke, or revolution, angle $\phi$ and the pitch, or rotation, angle $\psi$. The elevation, or flap, angle is kept at zero as in actual insect motions; it has a low amplitude compared with the other two angles (Shyy et al. Reference Shyy, Aono, Chimakurthi, Trizila, Kang, Cesnik and Liu2010). Figure 1 shows the angles’ definition. The pitch axis $z_{w}$ is attached to the wing, at quarter-chord, and stays in the $xz$ plane throughout the motion.

Figure 1. Definition of stroke and pitch angles, respectively $\phi$ and $\psi$, used to describe the wing motion.

The Reynolds number $Re$ describes the ratio between inertial and viscous forces. It is defined as

(2.1)\begin{equation} Re=\frac{\rho\bar{U}_{\phi} c}{\mu}, \end{equation}

where $\bar {U}_{\phi }=4bf \phi _0$ is the mean velocity at the tip of the wing during a stroke, induced only by the revolving motion of amplitude $2\phi _0$. In all cases throughout the paper, the Reynolds number is fixed to 1000, so that the flapping frequency $f$ can be deduced as

(2.2)\begin{equation} f=\frac{\mu Re}{4\rho cb\phi_{0}}. \end{equation}

This work is restricted to periodic motions of frequency $f$, where kinematics $(\psi, \phi )$ are parametrised by a few Fourier coefficients (typically the $f$ and $3f$ components), such as

(2.3)$$\begin{gather} \psi(t)= a_1 \cos(2{\rm \pi} ft) + b_1 \sin (2{\rm \pi} ft) + b_3 \sin(6{\rm \pi} f t), \end{gather}$$

(2.4)$$\begin{gather}\phi(t)=\phi_0 \sin(2{\rm \pi} ft-\delta), \end{gather}$$

where $\delta /2{\rm \pi} f$ corresponds to the time delay between the revolution and rotation motions. The stroke motion amplitude is thus defined as $2 \phi _0$. Note that such a parametrisation restricts the possible kinematic candidates for optimisation, yet the present definition allows for satisfying approximations of real insect flapping kinematics, further described in § 2.2. Additionally, the deep active learning framework proposed in § 3 can be extended to any other kinematic parametrisation if needed.

This study focuses on the hovering regime, where the motion objective is to produce a positive time-average lift, denoted $\bar {L}$. To achieve this lift performance, the system requires a rotation (respectively revolution) power $P_{rot}=T_{rot} \dot {\psi }$ (respectively $P_{rev}=T_{rev} \dot {\phi }$), $T_{rot}$ (respectively $T_{rev}$) being the rotation (respectively revolution) torque applied to the wing root. As a consequence, the efficiency $\eta$ can be defined as the ratio between the mean lift coefficient ($\overline {C_L}$) generated and the mean total power coefficient ($\overline {C_P}$), determined as

(2.5)$$\begin{gather} \overline{C_L}= \frac{\bar{L}}{\frac{1}{2}\rho b c \bar{U}^2}, \end{gather}$$

(2.6)$$\begin{gather}\overline{C_P}= \frac{\overline{P_{rot}}+\overline{P_{rev}}}{\frac{1}{2}\rho b c \bar{U}^3}, \end{gather}$$

(2.7)$$\begin{gather}\eta= \frac{\overline{C_L}}{\overline{C_P}}. \end{gather}$$

It is usual to use $\bar {U}_{\phi }$ as the normalisation velocity (Shyy et al. Reference Shyy, Aono, Chimakurthi, Trizila, Kang, Cesnik and Liu2010; Gehrke & Mulleners Reference Gehrke and Mulleners2021). Yet, the stroke amplitude is a degree of freedom in this study, and for kinematics with low amplitude, the stroke-induced velocity does not dominate the pitch-induced velocity. If the latter is not taken into account, a pure pitch motion would have an infinite lift coefficient, and the reference speed would not be adequate since $\bar {U}_{\phi }=0$ in that case. Even though the lower bound of the study is $\phi _0=10^\circ$, for such kinematics, the relevant speed governing the unsteady aerodynamics is no longer $\bar {U}_{\phi }$. As a consequence, another normalisation velocity $\bar {U}$ is thus introduced, different from that used for the Reynolds number calculation: it is defined as the time-averaged leading-edge velocity at the tip of the wing, induced by both the stroke and pitch motions:

(2.8) \begin{align} \bar{U}=\overline{\lVert \boldsymbol{u_{\phi}}+\boldsymbol{u_{\psi}} \rVert}=\frac{1}{T} \int_{0}^{T} \begin{bmatrix}\Big(\underbrace{b \dot{\phi}(t)}_{\boldsymbol{u_{\phi}}}+ \frac{c}{2} \dot{\psi}(t) \cos(\psi(t))\Big)^2 + \left(\frac{c}{2} \dot{\psi}(t) \sin(\psi(t)) \right)^2\end{bmatrix}^{1/2}\, \mathrm{d}t. \end{align}

Note that $\boldsymbol {u_{\phi }}$ is contained in the stroke plane, while $\boldsymbol {u_{\psi }}$ is perpendicular to the wing chord, thus the normalisation velocity $\bar {U}$ is not parallel to the stroke plane.

2.2. Kinematics generation

Real hovering kinematics from the hawkmoth, fruit fly and honeybee, described by Liu & Aono (Reference Liu and Aono2009), are used as a basis to generate the dataset used for neural network training. Figure 2 shows how their pitching motions are approximated using two or three Fourier coefficients. Using at least three Fourier coefficients is necessary to capture the reference pitching kinematics with reasonable accuracy, justifying the chosen pitching law in (2.3). Their stroke motion is close to sinusoidal, and their corresponding amplitude $2\phi _0$ is 120$^\circ$ for the hawkmoth, 90$^\circ$ for the fruit fly and 140$^\circ$ for the honeybee. The phase shift between revolving and pitching motions is 90$^\circ$ for all three cases.

Figure 2. Reference pitch $\psi$-kinematics (solid line) for (a) the fruit fly, (b) the honeybee and (c) the hawkmoth, as well as their approximation using ($a_1,b_1$) coefficients (dotted line), ($b_1,b_3$) (dash-dotted line) and ($a_1,b_1,b_3$) coefficients (dashed line) (see (2.3)).

To provide a large variety of data for the neural network training, 294 motions are initially created, of which 144 motions have a purely sinusoidal pitch profile (i.e. $a_1=b_3=0$) and 150 motions are composed using the three Fourier coefficients. For each motion, a uniform law between the bounds shown in table 1 is used to randomly generate each parameter ($a_1, b_1, b_3, \phi _0$ and $\delta$). The bounds were chosen to contain the parameters of the three bio-inspired reference kinematics for $a_1$, $b_1$, $b_3$ and $\delta$. For $\phi _0$, the higher bound is 60$^\circ$. We believe this bound does not affect significantly the findings in this study. For amplitudes higher than 120$^\circ$, the lift production mechanisms are similar to those described in § 4.2, but would need a mesh refined on a larger area to be simulated.

Table 1. Parameter bounds for the five parameters of the kinematics described by (2.3) and (2.4), used to generate randomly the 294 motions composing the initial dataset.

2.3. Computational set-up

DNS of the generated flapping motions is performed, and will later be used as the data for neural network training. The CFDs solver used to compute the flow around the flapping wing is StarCCM+. The three-dimensional Navier–Stokes equations are directly solved using a cell-centred finite volume method. Second-order numerical schemes are employed for spatial and temporal discretisations. The computational domain is a circular cylinder, 40$c$ high and 30$c$ wide, shown in figure 3. The wing moves parallel to the top and bottom ends of the far field. It is enclosed within a 2$c$-wide moving overset mesh. The displacement of the moving mesh is contained in a fixed control volume, in which the mesh is refined to resolve the wake with sufficient accuracy. The wing is modelled by a non-slip wall, the top of the domain is a stagnation inlet boundary condition, the bottom end is a pressure outlet and the lateral boundary conditions are treated as a slip wall.

Figure 3. Computational fluid domain with details on the refinement zones and the boundary conditions employed.

The spatial resolution chosen is $\Delta x=0.02c$ on the wing surface and $\Delta x=0.04c$ in the control volume, which yields a total number of cells of 1.9 million. The time resolution chosen is 250 time steps per period, where the period $T=1/f$ corresponds to an upstroke followed by a downstroke motion. The numerical set-up is similar to that used by Diaz-Arriba et al. (Reference Diaz-Arriba, Jardin, Gourdain, Pons and David2021), where it was shown that the results were converged with respect to both spatial and temporal resolutions. Four periods are simulated, and only the latest two are used in the dataset to avoid any transient effects in the data. Additional information on the decay of initial transients and convergence with respect to spatial and temporal resolutions is provided in the Appendix.

For flow analysis on optimal cases predicted by neural networks (in § 4), a finer mesh will be used. The spatial discretisation for these simulations is $\Delta x=0.005c$ on the wing surface and $\Delta x=0.01c$ in the control volume, leading to a total number of 36 million cells. The time discretisation is identical to the coarser mesh, i.e. set to $\Delta t=T/250$.

3.1. Unsteady aerodynamic model using deep learning

Lift, rotation torque and revolution torque are needed to compute $\eta$ and $\overline {C_L}$, and assess a motion performance. Thus, three neural networks are trained to predict those three outputs. The neural networks are feed-forward neural networks (FNNs), coded and trained using Pytorch (Paszke et al. Reference Paszke2019). The FNN employed here is an MLP, composed of successive layers of neurons. The output layer size of each neural network is 1, meaning each neural network will be trained to predict a scalar output at an instant $t$.

3.1.1. Neural network input

Multiple architectures exist to predict temporal data, such as recurrent neural network (RNN), LSTM or echo-state network (ESN) among many others. While being effective for speech recognition and other tasks involving temporal sequences, their use in computational physics is still sparse, mostly because only a short time history is needed to predict the further time steps. A typical example is CFD, where time-marching is achieved using time integration with very few previous time samples. As a consequence, a simpler MLP is chosen here where the input contains a few samples at previous time steps. To make a prediction for an output at time $t$, each neural network takes as input a vector $\boldsymbol {X}=[\boldsymbol {x}(t), \boldsymbol {x}(t-\tau \,\mathrm {d}t),\ldots, \boldsymbol {x}(t-(m-1)\tau \,\mathrm {d}t)]$, with $\boldsymbol {x}$ a state vector, $m$ and $\tau$ two integers chosen in the following, and $\mathrm {d}t=1/250f$. Such a method has been shown sufficient to replicate CFD results from high-order numerical codes, as for instance by Alguacil et al. (Reference Alguacil, Bauerheim, Jacob and Moreau2021, Reference Alguacil, Bauerheim, Jacob and Moreau2022) for acoustic waves propagation.

The state vector used is $\boldsymbol {x}=[\psi,\dot {\psi }, \ddot {\psi },\dot {\phi }]$ to provide information about the pitch and stroke kinematics to the neural network. Its components were chosen by testing different combinations. It appears that adding additional derivatives of the pitch angle improves the prediction performance, while for the revolving angle, only the first derivative is sufficient, as it follows a simpler sinusoidal law. Here, $m$ and $\tau$ are determined by performing a partial gridsearch. Increasing $m$ causes an increase in neural network size, as the input is larger and more neurons are needed to process it. It affects the training duration and prediction times once trained, thus $m$ should be chosen as low as possible. Figure 4 shows the $m$ and $\tau$ effect on the neural network performance, evaluated by the mean $R2$ score for lift prediction over a test set of 45 motions. For each motion in the test dataset, the $R2$ score is evaluated on both flapping periods (i.e. from $t=0$ to $t=2T$) and defined as

(3.1)\begin{equation} R2 = 1 - \frac{\displaystyle \sum_{t=0}^{2T} (y(t) - \hat{y}(t) )^2}{\displaystyle \sum_{t=0}^{2T} (y(t) - \bar{y} )^2}, \end{equation}

with $\bar {y}$ the average value of the output $y$ and $\hat {y}$ the prediction from the neural network.

Figure 4. Optimal $m$ and $\tau$ research, evaluated by the mean $R2$ score on the test set for lift prediction. The final choice $m=3$ and $\tau =9$ is indicated by dashed lines.

The chosen values are $m=3$ and $\tau =9$. As a consequence, for a prediction at the time step $t$, data are needed from $18\,\mathrm {d}t=0.072T$ earlier. This is in line with the results of Bayiz & Cheng (Reference Bayiz and Cheng2021): they demonstrated that aerodynamic states depended on the wing motion in less than half a cycle, thanks to cross-correlation analysis of the latent states of their model and kinematic variables. Here, $0.072T$ is slightly lower than the delay they found for maximal correlation, often approximately $0.1\unicode{x2013}0.2T$. The input vector is therefore $\boldsymbol {X}=[\boldsymbol {x}(t), \boldsymbol {x}(t- 9 \,\mathrm {d}t), \boldsymbol {x}(t-18\,\mathrm {d}t)]$. It is sampled from the motion data and then fed to the trained networks that produce the three needed outputs.

As detailed in § 2, a dataset of 297 simulated motions is created, containing the 294 randomly generated kinematics as well as the three reference ones from the fruit fly, the hawkmoth and the honeybee. This dataset of motions is divided in three parts. First, 70 % of the motions are used to train the neural network by updating the weights and biases in each neuron. Then, 15 % for validation to check the neural network performance during training on a set of data on which the weights are not updated. This is useful to prevent the overfitting (i.e. the memorisation rather than learning) of the training data. The last 15 % are used to test the final network on unseen data and compare the trained models. The three reference kinematics (figure 2) are placed in the testing dataset. Each motion, having been assigned to a dataset, is then split in input vectors $\boldsymbol {X}$ of length $4m=12$.

3.1.2. Prediction performance

The architecture chosen is identical for the three feed-forward neural networks: three hidden layers, with respectively 512, 1024 and 512 neurons, which correspond to a total of approximately 1 million tunable parameters per network. Each neuron activation function is a ReLU. The optimiser used to train the weights is Adam, with a learning rate of $10^{-4}$. The loss is the mean square error (MSE). Each neural network is trained for 50 epochs, with a batch size of 512.

Figure 5 shows the neural network prediction on the three outputs for the three reference motions. The data from the CFD simulation are also displayed. The neural network prediction generally fits the numerical data accurately.

Figure 5. Neural network predictions on the reference motions: (a–c) hawkmoth, (d–f) honeybee and (g–i) fruit fly. Predictions are performed over a period [$0, T$] for (a,d,g) the lift, (b,e,h) the revolution torque and (c, f,i) the rotation torque. The solid curve is the simulation results and the dashed line is the neural network prediction.

The best predicted output is the revolution torque with a mean $R2$ score over the test dataset of $\overline {R2}=0.985$, then the lift with $\overline {R2}=0.983$ and third the rotation torque with $\overline {R2}=0.948$. Overall, the neural networks show good performance in predicting unsteady aerodynamics for the three reference cases and over a wide range of motions as the $\overline {R2}$ scores show. These predictions of the unsteady aerodynamics take the order of 1 s to run on an 8-core CPU. This makes the neural networks suitable to be used as a surrogate model for an optimisation of the kinematics, as each evaluation of performance will be inexpensive compared with a CFD simulation (approximately 200 core-hours needed for each simulation in the dataset).

3.2. Multi-objective optimisation and active learning

3.2.2. Active learning algorithm

Neural networks are known to lose accuracy when generalising outside the training dataset. Since the kinematics have been generated randomly, it is unlikely they correspond to situations close to an optimum. As a consequence, the predicted Pareto front can be inaccurate because of the poor network generalisation capabilities. To evaluate this issue, CFD validations are performed on seven kinematics lying on the predicted Pareto front. The efficiency and lift coefficient obtained by CFD ($\triangledown$) are compared with the network predictions ($\circ$) in figure 7 (step 1, top left). Results show that the global trend of the Pareto front is predicted correctly by the networks, but large local errors can be also observed, especially for cases with high lift coefficients (the differences between the CFD and network predictions are highlighted by the dotted lines in figure 7). In fact, no training data exist in this zone, the maximum lift coefficient of the training dataset is ${\overline {C_L}} \approx 0.9$, whereas the neural network predicts cases with coefficients up to ${\overline {C_L}} \approx 1.5$.

To improve the accuracy of neural networks, and thus of the multi-objective optimisation, an active learning approach is used. Lye et al. (Reference Lye, Mishra, Ray and Chandrashekar2021) described an iterative process to improve partial differential equations optimisation using surrogate modelling. They tested it on an example of airfoil shape optimisation. One issue affecting the performance of deep learning surrogate models is that the training dataset is generated before the optimisation step, and the data points may be very different from the Pareto optimal points. The iterative surrogate model optimisation (ISMO) algorithm described by Lye et al. (Reference Lye, Mishra, Ray and Chandrashekar2021) allows for an iterative construction of the dataset with data points from successive optimisation steps. This process was adapted here to the gradient-free multi-objective optimisation algorithm NSGA II, described in § 3.2.1.

In this study, an initial training dataset of 208 motions ($70\,\%$ of the total dataset of 297 motions) was generated randomly to train the neural networks. The predictions from these first trained networks were detailed in § 3.1.2, and the validation of the Pareto front obtained with the NSGA II algorithm in figure 7 (step 1, top left). Now, $N=7$ motions from the predicted Pareto front are arbitrarily selected, and numerical simulations are performed to check the accuracy of the surrogate model on kinematics unseen during the training. These new data are then added to the training dataset for another step of neural network training and multi-objective optimisation. This active learning approach is illustrated in figure 6.

Figure 6. Complete framework used to discover optimal kinematics, involving the training of neural networks through an active learning approach to provide accurate predictions of the Pareto front by the NSGA II algorithm.

Two iterations of the active learning algorithm were needed to obtain an accurate prediction of the Pareto front. The corresponding Pareto fronts are displayed in figure 7. The predicted Pareto front is compared with the DNS of seven selected individuals. Step 1 corresponds to the surrogate model trained on the initial dataset. As already mentioned, it shows a good accuracy until ${\overline {C_L}} \approx 0.8$. For higher ${\overline {C_L}}$, the predicted Pareto front and DNS are further apart. Moreover, looking at the CFD simulations of the last two selected motions, which are linked to their initial predicted performance by a dotted line, the actual highest-${\overline {C_L}}$ motion is not actually the predicted one. The predicted Pareto front is not accurate in the high-${\overline {C_L}}$ area, leading to a wrong hierarchy between motions.

Figure 7. In the left column, Pareto fronts from the three steps of the active learning approach. The Pareto front predicted using the deep learning surrogate model ($\circ$) is compared with DNS ($\triangledown$) on $N=7$ individuals selected on the predicted Pareto front. The difference between network and DNS results is highlighted by the dotted lines. In the right column, the evolution of the prediction error of the model on the testing dataset in the ($\eta$, $\overline {C_L}$) plane. Dots ($\bullet$) are the DNS prediction and the lines point at the predicted position by the model.

Figure 8 displays the DNS of the Pareto fronts and the initial dataset ($+$). The reference kinematics (hawkmoth, fruit fly and honeybee) used to create the dataset are shown ($\star$). They are located near the Pareto fronts in the high efficiency area. The initial dataset was randomly generated, as described in § 2.2, and has no samples with $\overline {C_L} > 0.9$. This explains why the neural network surrogate model is struggling with accuracy at Step 1. The model behaviour on its testing set in figure 7 (step 1, top right) already hinted at a lack of precision towards the high lift. On this plot, the dots ($\bullet$) are the DNS and a line links it to the predicted position. For ${\overline {C_L}} >0.6$, errors are already clearly visible. Yet, it was able to predict high lift motions in an area where no data were available: even if errors are made, DNS shows that some optimal candidates have indeed a very high lift coefficient, up to ${\overline {C_L}} \approx 1.4$ when validated by DNS, while the network was predicting for this specific case ${\overline {C_L}} \approx 1.2$. As a result, the highest lift obtained by DNS at Step 1 is showing an increase of $55\,\%$ in $\overline {C_L}$ from the initial training dataset. At this point, it is thus remarkable that the deep learning model from Step 1 was able to discover new high lift kinematics, and predict their performance with a prediction error of $15\,\%$ on the lift coefficient. Similarly, on the other end of the Pareto front, the most efficient motions were improved by $12\,\%$, yet more training samples lie in this zone, which results in less accuracy problems for the neural networks.

Figure 8. Comparison between the DNS of the seven selected motions from each Pareto front and the initial dataset ($+$). The reference motions used to generate the dataset are highlighted ($\star$).

Following the active learning approach, these $N=7$ new motions are added to the dataset for Step 2. The predicted Pareto front is now more accurate throughout the whole range (figure 7b). There is a small over-estimation of the efficiency and lift coefficients for most of the front. The best motions for each objective were not improved. Figure 8 shows that, except around ${\overline {C_L}} \simeq 0.75$, the motions in Step 2 dominate those from Step 1. Intermediate motions were thus improved, in addition to the accuracy of the model. A third step leads to the most accurate predicted Pareto front, shown in figure 7(c): the CFD simulations match the predictions throughout the whole Pareto front. There is also a local improvement compared to the previous steps: in the area ${\overline {C_L}} \in [0.55,0.9]$, an increase in efficiency is observed compared with similar-${\overline {C_L}}$ motions from the first two steps, as visible in figure 8.

Overall, looking at the right column of figure 7, the accuracy over the whole testing set was improved by the two steps of active learning. The mean error (i.e. the mean size of the lines connecting DNS to neural net predictions on the plot) went from 0.081 for Step 1 to 0.056 for Step 3. The training dataset has been densified near the Pareto front by adding twice seven new motions, but it also benefited for the model performance over the whole range of motions.

The active deep learning process needed 314 CFD simulations of around 2 hours (running on 4 CPU nodes Intel Skylake 6126, 24 cores, 96Go RAM). As mentioned in § 3.2.1, the genetic algorithm converges after 1800 evaluations of the cost function. Using the surrogate model, each evaluation takes the order of a second. Thus, the total cost for the active learning process and optimisation is approximately 630 h. If we were to use CFD for the 1800 evaluations of the optimisation process, it would take approximately 3600 h. Our active deep learning model has thus allowed us to perform a two-objective optimisation for 17 % of the CFD cost.

Figure 9 allows to compare the evolution of three parameters and the differences of their distributions from Step 1 to Step 3: stroke amplitude $2\phi _0$ (a), the phase shift $\delta$ (b) and the stroke-average angle of attack $\bar {\alpha }$ (c). The latter is defined as the average angle during a stroke between the wing and the stroke plane $xz$. On the Step 1 Pareto front, there is a harsh change in stroke amplitude $2\phi _0$ around $\overline {C_L}=0.5$ from 120$^\circ$ to 20$^\circ$, whereas on the Step 3 front, there is a gradual decrease, passing by motions with an intermediate amplitude. For the phase shift, the decrease is monotonous when increasing the lift coefficient along the Step 3 front, while for Step 1 around $\overline {C_L}=0.8$, there is a shift in trend. Similarly, the high values of phase shift around $\overline {C_L}=1.2$ are attained more gradually in Step 3. The same remarks can be made for the average angle of attack. Overall, the evolution of those three quantities is smoother and monotonous on Step 3. The difference in trends in the Step 1 Pareto front is linked to its lack of precision, mixing the hierarchy between motions. The active learning process thus allowed a smoother transition between high-efficiency and high-lift motions to be found.

Figure 9. Colour-coded Pareto fronts highlighting the evolution of three parameters $\phi _0$, $\delta$ and $\bar {\alpha }$, for both the Step 1 and Step 3 Pareto fronts. Two regimes are identified (dashed lines) in the Step 3 Pareto front, one at high efficiency (1) and the other at high lift coefficient (2).

4.1. Pareto front analysis

To better understand what parameters drive these optimal kinematics, the evolutions of the stroke amplitude, phase shift and average angle of attack on the Step 3 Pareto front are observed in figure 9.

First, the stroke amplitude colouring in figure 9(a) highlights the difference between the motions at each end of the Pareto front. For high-efficiency motions, until $\overline {C_L} \simeq 0.6$ (regime 1), the stroke amplitude takes values between 90 and $120^\circ$. Then, there is a drop in amplitude to ${\simeq }50^\circ$ and a gradual decrease reaching 20$^\circ$ from $\overline {C_L} = 0.9$ onwards (regime 2). This stroke amplitude difference between high-lift and high-efficiency motion hints at two distinct lift production mechanisms. The high-lift kinematics having a shorter stroke amplitude compared with the high-efficiency ones is in line with the results from Altshuler et al. (Reference Altshuler, Dickson, Vance, Roberts and Dickinson2005). They concluded that at constant wingtip velocity $U_\phi$, decreasing stroke amplitude promoted lift while decreasing efficiency. The two regimes are highlighted by the dashed lines: (1) corresponds to high-efficiency kinematics and (2) corresponds to high-lift motions. They can be characterised by their slope in the (${\overline {C_L}}, \eta$) plane. For regime 1 (i.e. from $\overline {C_L}=0.32$ to $0.61$) the efficiency decreases rapidly, with a slope $\Delta \eta /\Delta \overline {C_L}=-2.12$. This indicates that to gain a certain amount in lift coefficient, it costs roughly twice the amount in efficiency. For regime 2 (i.e. from $\overline {C_L}=0.6$ to $1.4$), the decline in efficiency is slower, with a typical slope $\Delta \eta /\Delta \overline {C_L}=-0.74$. A gain in lift coefficient is not as costly in efficiency. Motions from those two regimes will be discussed in § 4.2 for regime 1 and § 4.3 for regime 2. It is also to be noted that though the Reynolds number based on the revolution speed $\overline {U_\phi }$ is 1000 for all motions, the Reynolds number based on the force-normalisation velocity $\bar {U}$ described in (2.8) varies between 400 (for high-efficiency cases) and 1300 (high-lift cases) along the Pareto front.

Second, the phase-shift colour-coded Pareto front (figure 9b) reveals a gradual evolution throughout the front. For efficient motions until $\overline {C_L}\simeq 0.65$, the phase shift is between 90$^\circ$ and 110$^\circ$. That is, the wing angle passes through vertical slightly before stroke reversal. This kinematic property is commonly referred to as ‘advanced rotation’ (e.g. Sane & Dickinson Reference Sane and Dickinson2001). Then for ${\overline {C_L}} \in [0.65,1.2]$, $\delta$ varies between 110$^\circ$ and 135$^\circ$. The highest lift motions (${\overline {C_L}} >1.2$) have higher values of phase shift, from 150$^\circ$ to 170$^\circ$. Those high values of phase shift may seem counter-intuitive at first, and imply that the trailing edge of the wing is ahead of the leading edge during a large portion of the stroke. This relatively uncommon flapping regime will be further studied in § 4.3.

Third, the evolution of the stroke-average angle of attack $\bar {\alpha }$ on figure 9(c) shows an overall smooth increase from values of 37$^\circ$ for high-$\eta$ motions to 81$^\circ$ for high-lift motions. This angle is influenced by the pitch kinematics and the phase-shift evolution. While values of 37$^\circ$ for high-efficiency motions are in line with previous results from Sane & Dickinson (Reference Sane and Dickinson2001), Gehrke & Mulleners (Reference Gehrke and Mulleners2021) and Diaz-Arriba et al. (Reference Diaz-Arriba, Jardin, Gourdain, Pons and David2021), for example, values of 81$^\circ$ appear very large. This, again, suggests different mechanisms of force production between high-efficiency and high-lift motions. We note that the effective angle of attack of the wing is in fact much lower than $\bar {\alpha }$ due to the downwash induced during previous strokes (as a result of lift production).

4.2. High-efficiency kinematics

The most efficient flapping motion is studied here, as a typical example of kinematics of regime 1. Its parameters are $a_1=0.001, b_1=-1.34, b_3=-0.22, \phi _0=55^\circ$ and $\delta =94^\circ$. This motion is characterised by a high stroke amplitude of 110$^\circ$ and an overall low angle of attack. The time evolution of the motion's parameters are displayed in figure 10, along with the lift and power coefficient evolution over a flapping period, obtained by DNS and the neural network model. The neural network error on the prediction is 2.8 % for $\eta$ and 1.8 % for $\overline {C_L}$. It is interesting to note that this motion highly resembles the three reference kinematics observed in nature (figure 2). The pitch profile is trapezoidal: the pitch angle is almost constant at 65$^\circ$ during most of the stroke. This corresponds to an angle of attack of 35$^\circ$, which was also observed in real hovering kinematics for small insects by Ellington (Reference Ellington1984).

Figure 10. Time evolution of the (a) kinematics parameters, and (b) lift and power coefficients for the most efficient motion, obtained from DNS (black) and neural networks (grey). Key instants are located by blue dots.

The lift coefficient evolution (figure 10) also has an almost trapezoidal shape during a single stroke ($t\in [0,T/2]$): a fast increase at the start of the stroke, then little variations in lift during the almost constant pitch angle phase, followed by a fast decrease at the end of the stroke. The lift seems mainly controlled by $\dot {\phi }$ and $\psi$, which suggests that quasi-steady effects are dominant. Note that for this high-efficiency motion, the lift and the power coefficients evolution are in phase. During the part of the stroke when lift is produced, the power coefficient is half the value of the lift coefficient. This leads to the efficiency value $\eta = {\overline {C_L}} / \overline {C_P}=1.79$.

To further analyse the dependency of lift to $\dot {\phi }$ and $\psi$, results are compared with predictions from the quasi-steady model by Lee et al. (Reference Lee, Lua, Lim and Yeo2016), adapted to use our normalisation velocity $\bar {U}$ defined in (2.8). The model includes effects from pure revolving motions, combined revolving and rotating motions, as well as from added mass. The total quasi-steady lift coefficient reads

(4.1)\begin{equation} C_L = C_{L,rev} + (C_{L,rot 1}+C_{L, rot 2}) \cos{\alpha} + C_{L, am} \cos{\alpha}, \end{equation}

where $C_{L,rev}$ is the quasi-steady lift force of a wing under pure revolving motion and is expressed as

(4.2)\begin{equation} C_{L,rev}=(f_{AR,tr}) (f_{Ro,tr}) \frac{b^2 C_{L,tr}}{3\bar{U}^2} \dot{\phi}^2, \end{equation}

where $C_{L,tr}$ is the quasi-steady lift coefficient resulting from a translation motion and is equal to $1.966 - 3.94 Re^{-0.429} \sin ({2\alpha })$. Here, $f_{AR,tr}$ and $f_{Ro,tr}$ are correction factors associated with aspect ratio and Rossby number effects (see Lee et al. (Reference Lee, Lua, Lim and Yeo2016) for further details).

Additionally, $C_{L,rot 1}$ and $C_{L,rot 2}$ are the quasi-steady lift forces of a wing undergoing combined rotating and revolving motion. They are derived as

(4.3)\begin{equation} C_{L,rot 1}=(f_{\alpha}) (f_{r}) C_{rot} \frac{ c b |\dot{\phi}| \dot{\alpha} }{\bar{U}^2}|\dot{\phi}| \dot{\alpha} \end{equation}

and

(4.4)\begin{equation} C_{L,rot 2}= 0.778 \frac{c^2 |\dot{\alpha}| \dot{\alpha} }{\bar{U}^2}, \end{equation}

where $f_{\alpha }$ and $f_{r}$ are correction factors linked to the rotational axis position and the instantaneous angle of attack. Here, $C_{rot}$ is equal to $0.842 - 0.507Re^{-0.1577}$.

Finally, the added mass coefficient is derived as

(4.5)\begin{equation} C_{L,am}= (f_{a}) (f_{AR,a}) \frac{\rm \pi}{4 \bar{U}^2} \left(c b \dot{\phi} \sin{\alpha}+\frac{c^2}{2}\dot{\alpha} \right), \end{equation}

with $f_{a}$ and $f_{AR,a}$ correction factors empirically determined.

Figure 11 compares the lift predicted by the quasi-steady model with that obtained from the present simulations during the first and third flapping cycles. The lift coefficients obtained with the numerical simulation, during the first flapping cycle, and the quasi-steady model are found to match reasonably well. In particular, the lift generally increases up to $C_L\approx 0.9$ slightly before $T/4$ and then decreases before increasing again up to $C_L\approx 0.8$ around $3T/4$. The quasi-steady model suggests that the first increase is principally driven by quasi-steady effects from the revolving motion, and hence by kinematic parameters $\dot {\phi }$ and $\psi$. Added mass and rotational effects tend to shift the maximum lift to earlier times, when compared with that obtained from the sole contribution of the revolving motion (which occurs at $T/4$). However, the occurrence of the second peak in lift appears to be principally driven by rotational (pitch-up) effects. Overall, we note that quasi-steady effects from the revolving motion contribute positively to lift production throughout the whole stroke while rotational/added mass effects contribute negatively/positively during the first half of the stroke and positively/negatively during the second half, respectively. Finally, by comparing the lift obtained during the first and third cycles, it is shown that wake effects (present during the third cycle but virtually absent during the first cycle) contribute to reducing lift. These effects can be classified into two contributions: one associated with the overall downwash induced by lift production during previous strokes and another one associated with local effects of wing–vortex interactions. The overall downwash tends to reduce the effective angle of attack of the wing and hence can be viewed as a general offset of the lift curve towards lower values. Wing–vortex interactions, however, may have very local effects, especially at stroke reversal where the lift is here found to reach negative values.

Figure 11. Quasi-steady lift model (black line) from Lee et al. (Reference Lee, Lua, Lim and Yeo2016) and its different components (translational, rotational and added mass effects), compared with the DNS results during the the first (blue line) and third (red line) upstrokes for the most efficient motion.

The occurrence of wing–vortex interactions can be investigated by looking at the flowfield in the vicinity of the wing. To this aim, isosurfaces of non-dimensional $Q$-criterion, defined as $Qc^2/\overline {U_{\phi }}^2$, are displayed in figure 12 at different instants throughout the motion for half a period (i.e. $t/T$ from $0$ to $0.5$, corresponding to instants (a)–( f) in figure 10a). Two isosurfaces corresponding to values of 1 (grey) and 10 (blue) are shown to visualise the evolution of the vortex structures. The video corresponding to these snapshots is available in the supplementary movies available at https://doi.org/10.1017/jfm.2023.832. At the start of the stroke in panel (a), the wing is at its leftmost position and pitches clockwise. A trailing edge vortex (TEV) is visible in the vicinity of the trailing edge as a result of the rotation of the wing between strokes. The leading-edge vortex (LEV) generated during the previous stroke is also visible in the vicinity of the wing. In panel (b), the wing revolves from left to right. The LEV from the previous stroke sweeps on the lower surface, which contributes to negative lift by acting as a low-pressure region beneath the wing. In panels (c) to (e), a new leading-edge vortex is formed. It is cone-shaped and attached to the wing, as the velocity at the tip of the wing is higher than that near the wing root. It reaches a quasi-steady state, which corresponds to the observations made on the lift in figure 10. It is known that the leading-edge vortex is key to delay stall on low aspect ratio revolving wings and contributes to increasing the lift production throughout the stroke, by generating a low pressure area on the upper surface of the wing (Ellington et al. Reference Ellington, Van Den Berg, Willmott and Thomas1996; Dickinson et al. Reference Dickinson, Lehmann and Sane1999; Eldredge & Jones Reference Eldredge and Jones2019). A tip vortex (TV) is also visible. In panels (c) and (d), the TV and LEV appear as two distinct structures. In panel (e), the wing starts its counter-clockwise pitching motion to reverse for the next stroke. The wing angle passes through vertical before the wing reaches its rightmost position. This (slightly advanced) rotation mechanism increases the angle of attack and the leading-edge vortex, initially cone-shaped in panel (d), expands chord-wise and covers most of the wing while still being attached. In panel (e), the TV and LEV merge at the tip of the wing. Then, in panel ( f), the leading-edge vortex merges with the 3-D structures of the tip vortex as the wing fully rotates for the next stroke. The decrease in revolution speed $\dot {\phi }$ of the wing weakens the spanwise pressure gradient stabilising the leading-edge vortex, as described by Medina & Jones (Reference Medina and Jones2016). The axial flow along the leading edge decreases as the pitch movement dominates in this reversal part of the flapping. The leading edge vortex loses its quasi-steady state, and this is linked to the fast and sudden drop in lift after the second peak.

Figure 12. Non-dimensional $Q$-criterion isosurfaces obtained at different instants $t/T$, from $0$ to $0.5$, for the most efficient motion. The grey isosurfaces correspond to a value of 1 and the blue to a value of 10.

4.3. High-lift kinematics

The flapping motion generating the highest lift is investigated in this section, as a typical example of motions in regime 2. The parameters are $a_1=0.186$, $b_1=-0.586$, $b_3=-0.3$, $\phi _0=10^\circ$ and $\delta =167^\circ$. The corresponding $\psi$ and $\phi$ time evolutions are displayed in figure 13, along with the lift and power coefficient evolution over a flapping period, obtained by DNS and the neural network model. The neural network error on the prediction is 1.9 % for $\eta$ and 0.9 % for $\overline {C_L}$. The motion is very different from the most efficient one and the three reference kinematics: it is a low revolution amplitude motion, with a complex pitching motion during the stroke. Due to the large phase shift between pitching and revolving motions, the pitch angle $\psi$ changes sign very early in the stroke, i.e. soon after $t/T = 0.25$, meaning that the wing is oriented backwards during a large portion of the stroke. This is in contrast with previously studied motions (regime 1), where $\psi$ changed sign around stroke reversal, i.e. at $t/T \approx 0.5$. Moreover, the pitching motion during a stroke can be divided in two phases: a main pitch-up motion ($t/T \in [0.2-0.4]$), and smaller pitching variations at the beginning and end of a stroke.

Figure 13. Time evolution of the (a) kinematics parameters, and (b) lift and power coefficients for highest-lift-generating motion, obtained from DNS (black) and neural networks (grey). Key instants are located by blue dots.

The lift coefficients produced in the upstroke (i.e. $t/T \in [0,T/2]$) and downstroke (i.e. $t/T \in [T/2,T]$) are different, due to the pitching motion not being fully symmetrical. The mean lift coefficient in the upstroke is 1.67 and is 1.19 in the downstroke. In both strokes, there are two peaks of lift production, the first one around half of the stroke and a second peak, smaller and flatter, towards the end of the stroke. During the upstroke, the first peak reaches $C_L=7.5$, while the second one reaches $C_L=0.96$. For the downstroke, the peaks reach maximum values of $C_L=4.75$ and $2$, respectively. The maximum lift value may change, but the overall lift evolution is similar between both strokes, thus a single stroke is analysed to understand the underlying flow characteristics leading to this performance.

Similarly to the previous subsection, the quasi-steady model from Lee et al. (Reference Lee, Lua, Lim and Yeo2016) is used to estimate the main contributions to the lift peaks. Figure 14 compares the lift predicted by the quasi-steady model with that obtained from the present simulations during the first and third flapping cycles. First, it is observed that the quasi-steady model predicts a lift peak slightly before $t=T/4$. This peak occurs at a similar instant to that obtained with numerical simulations but reaches lower levels. It appears to be dominated by both added mass effects and rotational terms. The difference in peak amplitude between simulations and the quasi-steady model suggests that a part of the lift is produced by unsteady, vortex-induced effects. This point will be further addressed later in this subsection. Second, the second small peak is well predicted by the quasi-steady model and seems to be almost exclusively driven by added mass effects. Finally, it is shown that the amplitude of the lift peak is larger during the third flapping cycle than during the first flapping cycle. That is, wake effects appear to promote unsteady effects responsible for lift production near mid-stroke. However, the lift is slightly reduced before $t\approx T/6$ and very slightly increased after $t\approx T/3$. It is striking that the overall downwash induced by lift production during previous strokes does not have a detrimental influence on lift, conversely to what has been observed for the most efficient case analysed in the previous subsection. This is in line with the weak contribution of revolving motion-induced quasi-steady forces to the overall lift.

Figure 14. Quasi-steady lift model (black line) from Lee et al. (Reference Lee, Lua, Lim and Yeo2016) and its different components (translational, rotational and added mass effects), compared with the DNS results during the the first (blue line) and third (red line) upstrokes for the most efficient motion.

Figure 15 displays non-dimensional Q-criterion isosurfaces throughout one stroke, for values of $Qc^2/\overline {U_{\phi }}^2=60$ (grey) and $200$ (blue). The video corresponding to these snapshots is available in the supplementary movies. Because of the low revolution amplitude and the complex 3-D vortical structures, spanwise vorticity obtained in the $r=0.75b$ spanwise section is also displayed at various times in figure 16 to ease the visualisation. In the snapshot in panel (a), the wing is at its leftmost position and it initiates revolution from left to right. Three vortices are visible, generated in the previous strokes: one leading-edge vortex (LEV1) and two counter-rotating trailing edge vortices (TEV1 and TEV2). The TEV positions throughout the motion are highlighted using coloured dashed lines in figure 15. The line is pink if the vortex is created by the main pitch motion (panels c–e), and it is green if created by variations of smaller amplitude between strokes. In panels (a) and (b), the wing pitches clockwise and generates negative lift. LEV1 strongly interacts with the leading edge and promotes the generation of counter-rotating vorticity which appears to shed in panel (c). From panels (c)–(e), the wing performs a counter-clockwise pitching motion while revolving from left to right. The angle of attack increases from 25$^\circ$ to 107$^\circ$. This pitch-up motion creates the main lift peak. The combination of the rotation and revolution implies a high trailing edge speed where the flow separates and rolls up into a trailing edge vortex TEV3. The fast generation of this intense TEV3 is clearly visible from panels (c) and (d). A leading-edge vortex LEV2 is also generated during this phase, as highlighted in panel (d). The fast formation of both leading and trailing edge vortices from panels (c)–(e) generates a large time variation in the first moment of vorticity which, following the impulse theory (Wu Reference Wu1981), is directly related to force production. Interestingly, this mechanism resembles that recently observed for mosquitoes (Liu et al. Reference Liu, Du and Sun2020), which flap their wings with a relatively small amplitude, i.e. typically 45$^\circ$. In panels (c) and (e), the leading edge vortex develops close to the wing but remains small due to the limited revolving amplitude. Lift generation is here very different to that observed in the high-efficiency case where the LEV on its own, by acting as a low-pressure region on the upper surface of the wing, is the main contributor to lift generation. In other words, the ‘delayed stall’ mechanism is not relevant for this low stroke amplitude motion because the leading-edge vortex does not have enough time to form (Wang et al. Reference Wang, Birch and Dickinson2004; Zhu & Sun Reference Zhu and Sun2017; Liu et al. Reference Liu, Du and Sun2020). In panel (e), LEV2 is wrapped by vorticity braids around the mid-span that result from the development of 3-D spanwise instabilities. The latter may be triggered by the relative proximity of 3-D instabilities that have already developed in LEV1 and TEV2.

Figure 15. $Q$-criterion isosurfaces obtained at different instants $t/T$ for the highest-lift-generating motion. The grey isosurfaces correspond to a value of 60, and the blue to a value of 200. Red arrows indicate the direction of the pitching motion of the wing, and dashed lines track the evolution of particular vortices.

Figure 16. Spanwise vorticity contours obtained in the $r=0.75b$ spanwise cross-section at different instants $t/T$ for the highest-lift-generating motion.

Between panels (e) and ( f), the revolving speed decreases to zero and the wing pitches clockwise. TEV3 and LEV2 can be viewed as a vortex dipole that generates a fluid jet oriented towards the trailing edge. Together with the trailing edge speed induced by the pitching motion, this jet leads to the formation of TEV4 as the flow separates past the trailing edge. This mechanism was also observed for mosquitoes (Bomphrey et al. Reference Bomphrey, Nakata and Phillips2017). Again, the generation of TEV4 contributes to the time rate of change of the first moment of vorticity, and hence aerodynamic force. This mechanism, which is favoured by the interaction with a TEV3 and LEV2 induced jet, is found to be correlated with the occurrence of the second small lift peak in figure 14, although it was suggested that the latter is driven by added mass effects. At the end of the stroke, LEV2 seems to burst into smaller scale structures with vorticity braids developing along its full span.

Then the next stroke begins, involving similar phenomena. We note from panel (g), and similar to panel (b), that the contra-rotating vortex dipole formed by TEV2 and TEV4 rapidly advects downstream while LEV2 interacts with the leading edge of the wing.