# Hovering flapping wings with dynamic twist

Florian Bouard, Thierry Jardin, Laurent David

► To cite this version:

Florian Bouard, Thierry Jardin, Laurent David. Hovering flapping wings with dynamic twist. Physics of Fluids, 2024, 36 (11), pp.111910. 10.1063/5.0234656. hal-04864113

HAL Id: hal-04864113

https://hal.science/hal-04864113v1

Submitted on 4 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Hovering flapping wings with dynamic twist

Florian Bouard,¹,² Thierry Jardin,¹ and Laurent David²

¹) ISAE-SUPAERO, Université de Toulouse, 10 Av. Edouard Belin, Toulouse cedex 4,
31055, France

²) Université de Poitiers, ISAE-ENSMA, CNRS UPR 3346, Institut PPrime,
11 Bld. M. et P. Curie, TSA 41123, Poitiers cedex 9, 86073,
France

(Dated: 13 December 2024)

The role of (active) dynamic wing twist on aerodynamic performance of three-dimensional hovering flapping flight is explored using numerical simulations. A variety of cases with different pitch angles and with (flexible wings) or without (rigid wings) dynamic twist are compared. The results show that changes in aerodynamic performance due to dynamic twist are comparable to those obtained without twist (rigid wing cases) by pitching the whole wing and that lift and lift-to-power ratio generally collapse onto a single curve when plotted as a function of the mid-stroke pitch angle at 2/3 wing radius. However, in some cases dynamic twist yields enhanced time-averaged efficiency. Using the force and power partitioning method (FPPM), it is shown that this enhancement results from the absence of vortical structures near the wing root lower surface and to the presence of an extended leading edge vortex on the wing upper surface, when compared to the most efficient rigid wing case. These differences in flow topology lead to enhanced lift during the early phase of the strokes without changes in power consumption.

<footnote>(*Electronic mail: florian.bouard@univ-poitiers.fr)</footnote>
# I. INTRODUCTION

The development of micro-electronics has allowed the continuous reduction in the size of unmanned aerial vehicles. At extremely small scales, typically of the order of the centimeter and below, flapping wings appear as an alternative to more conventional rotating and fixed wings. Accordingly, they have fostered considerable interest from researchers these past three decades. Over this period, most studies focused on rigid wings undergoing flapping kinematics both in forward and hovering flight¹⁻³. Pioneer works revealed the predominant role of vortices (induced by the flapping motion) in lift production as recently summarized by Liu et al.⁴, for example. Precisely, the formation of a strongly attached, stable leading-edge vortex (LEV)⁵⁻¹⁰ on the suction side of the wing has been recognized as the principal source of lift allowing flapping flight.

Despite obvious representation of wing bending¹¹ and twisting¹² in nature, flexible flapping wings have been much less studied. Passive deformation of the wing received more and more attention in recent years³ with studies aiming at understanding the effects of structural parameters on aerodynamic performance of flapping wings. One general conclusion is that flexibility has a major influence on aerodynamic performance. That is, the force production increases with the flexibility (e.g. with the increase in mass or frequency ratio) until a specific flexibility is reached¹³⁻¹⁵. Further increasing the flexibility appears detrimental to force production. For example, Zhao et al.¹³ showed that at high angles of attack, flexible wings generate more lift due to the orientation of the net force vector induced by the chordwise flexion of the wing. Dai et al.¹⁴ found that when flapping frequency is lower than 0.3 times the natural frequency of the wing, the efficiency can be significantly improved. The mechanisms at play rely on the natural oscillation of the wing working with the aerodynamic forces to bend the chord, which leads to larger pitch angle with respect to rigid wings. Shah et al.¹⁵ suggested that an optimum level of flexibility inhibits chaotic transition of the wake, otherwise encountered in rigid and highly flexible wing cases. They ultimately correlated the quasi-periodicity of the wake with enhanced efficiency of flexible foils, contrasting the chaotic nature of the flow, and lower efficiency, of rigid wings.

The deformation of flexible wings undergoing flapping kinematics can be decomposed into three different modes¹⁶: spanwise bending, spanwise twisting and chordwise cambering illustrated in Figure 1 a), b) and c), respectively, where the undeformed wing geometry is shown using black plain lines and the deformed wing geometries are depicted using dotted, grey lines.

The role of chordwise flexibility has been investigated by Heatcote et al.¹⁷, Vanella et al.¹⁸ and
Hovering flapping wings with dynamic twist

<figure><img src="image_2.png" /><figcaption>FIG. 1. Illustration of flapping flexible wing undergoing different deformation$^{16}$: a) spanwise bending, b) spanwise twisting and c) chordwise cambering. Plain lines depict rigid wings while dotted lines illustrate the flexible wing undergoing different deformation mechanisms.</figcaption></figure>

Eldredge et al.¹⁹, for example. Heatcote et al.¹⁷ considered an airfoil started from rest and under-
going pure plunging motion with no freestream velocity. They showed that depending on plunge
frequency and amplitude, an optimal flexibility was found to promote the ratio between thrust and
input-power. Vanella et al.¹⁸ and Eldredge et al.¹⁹ considered 2D flapping wings under hovering
conditions. The conclusion from both studies was that the flexibility significantly modifies the
interactions between the wing and its wake from the previous stroke. Specifically, the LEV from
previous stroke has almost no influence on the development of a new LEV which results in larger
lift production with respect to that obtained with rigid wings.

The role of spanwise flexibility was addressed by Heatcote et al.²⁰, Nakata and Liu²¹, Addo-
Akoto et al.²² and Diaz-Arriba et al.²³ for various flight conditions (forward and hovering flight).
Overall, it was demonstrated that moderate spanwise flexibility promotes force production (thrust
or lift) by tuning the phase difference between revolving (affected by spanwise bending) and pitch-
ing (affected by spanwise twisting) motions and/or promoting the stability of the leading edge
vortex (LEV).

On the other hand, active wing deformation has been poorly investigated. By considering both
coupled and uncoupled wing cambering and twisting for insect wings under hovering flight condi-
tions, Gang and Sun²⁴,²⁵ concluded that camber has a large impact on lift production while twist
has a negligible impact. Conversely, they showed that the power required to achieve flapping
flight is reduced with wing twisting. This conclusion is in line with the ones from Thielicke and
Stamhuis²⁶ who conducted an experimental study in the context of a slow-speed forward flapping
Hovering flapping wings with dynamic twist

flight. They showed that wing twisting, by lowering the effective angle of attack (with respect to a rigid wing), is most likely to increase the efficiency at the cost of lower aerodynamic forces production when compared to a rigid wing, emphasizing the trade-off between flapping flight efficiency and force production. More recently, Joshi and Bhattacharya²⁷ performed experimental characterization of an accelerating flat plate with active spanwise bending. The wing was held normal to the flow and bent into or away the incoming flow. While the into-flow bent wing experienced larger unsteady drag, the away-flow bent wing experienced significantly lower initial force peak. This suggests that wing active deformation might be used to reach specific aerodynamic forces at a given instant for maneuverability purposes, for example. Soto and Bhattacharya²⁸ conducted experiments on an actively twisted heaving flat plate in an incoming flow with a constant angle of attack of 15°. They concluded that actively twisting the wing such that the effective angle of attack is reduced results in lowering significantly the drag down to almost zero at the cost of lower lift production. On the contrary, twisting the wing such that the effective angle of attack is increased increases both lift and drag during the downstroke phase while the wing experiences negative lift during the upstroke phase. Dong et al.²⁹ performed numerical simulations on actively twisted flapping wings in forward flight. Wing twisting was achieved such that the effective angle of attack is constant along the span. By maintaining a proper spanwise distribution of the effective angle of attack, the LEV stability was enhanced resulting in larger time-averaged thrust and vertical force production, which contrasts previous studies²⁴⁻²⁶. Finally, Bouard et al.³⁰ numerically investigated the effect of prescribed bending deformation on 3D hovering flapping wings. The bending-induced phase-lag has been systematically varied to assess its influence and out-of-phase bending (i.e. wing bending is maximum at stroke reversal) was found to promote lift production. It has been shown that this performance enhancement is correlated with larger flapping amplitude of the morphing wing with respect to the rigid wing, which increases wing velocity (for a given flapping frequency) and thus lift production.

Despite these few studies, the role of active deformation on aerodynamic performance of flapping wings remains poorly addressed. A general conclusion is that such morphing may enhance efficiency at the cost of lower force production, for both hovering conditions²⁴,²⁵,³¹ and low-speed forward flight²⁶. Yet, there is no general consensus on the mechanisms that drive these changes in aerodynamic performance. Therefore, in this paper, we systematically evaluate the influence of dynamic wing twist over a wide range of pitch angles and show that aerodynamic performance scales, to leading order, with the time-averaged pitch angle at 2/3 span. Furthermore, we pro-
Hovering flapping wings with dynamic twist

vide a detailed analysis of the correlation between changes in aerodynamic performance and flow
topology. For that purpose, we adapt the force and moment partitioning method (FMPM) to the
analysis of force generating and power consuming structures (that we refer to as the force and
power partitioning method (FPPM)), which helps explain in an unambiguous way the increase in
efficiency without loss of force production in some specific cases.

## II. NUMERICAL METHODS

We consider a three-dimensional flapping wing under hovering flight conditions. The wing profile is a NACA0012 airfoil and the aspect ratio AR is $R/c = 4$, which was found to be the optimal lift aspect ratio for a revolving wing at low Reynolds number$^{32}$. The wing chord $c = 0.01\text{m}$ is constant along the span. The setup is similar to that in our previous studies$^{23,30,33}$, where rigid wings, actively bending wings and passively deforming wings have been considered. Here, we focus on actively twisting wings and compare their aerodynamic performance to that of rigid wings. In rigid wing cases, the flapping kinematics consist in combined, sinusoidal revolving and pitching motions. Their respective angular speeds are denoted by $\phi$ and $\dot{\alpha}$, as illustrated in Figure 2 (bottom right). The wing flaps along a horizontal plane, over a revolving amplitude $\phi_0 = 120^\circ$ during half a flapping period $T/2$ (with $T \approx 0.098\text{s}$). The same, minimum, pitch angle $\alpha_0$ is attained at mid-downstroke and mid-upstroke, as shown in Figure 2 (top right), rendering the motion symmetric. The pitch axis is positioned one quarter chord away from the leading edge, i.e. between the leading edge and the mid-chord line, which has been shown to be the optimal range for pitch axis location$^{34}$. In actively twisting wing cases, a sinusoidal dynamic twist motion is added to the revolving and pitching motions. Linear spanwise twist is considered, which has been observed in natural species such as mosquitoes$^{35}$. The twist axis is located on the pitch axis (one quarter chord away from the leading edge). Dynamic twist is here in phase with the pitching motion. In other words, the morphing wing position coincides with that of the rigid wing at stroke reversal and undergoes maximum twist at mid-downstroke and mid-upstroke. The twist angle $\epsilon$ is defined such that the wing tip pitch angle is $\alpha_{tip} = \alpha + \epsilon$. This way, $\alpha_{tip}$ can be both greater or lower than $\alpha$ depending on the sign of $\epsilon$. $\epsilon_0$ (amplitude of the dynamic twist motion) is bounded such that $\alpha_{tip}$ can not be lower than 0° or exceed 90°. Figure 2 (bottom right) illustrates a morphing wing flapping motion with $\alpha_0 = 45^\circ$ and $\epsilon_0 = -25^\circ$ (depicted in lightgrey). Note that, similarly to rigid wing cases, the kinematics of the morphing wing case is symmetric. We further stress that
Hovering flapping wings with dynamic twist

<figure><img src="image_3.png" /><figcaption>FIG. 2. Sketch of the computational domain (left), revolving and pitching speeds for both downstroke (white area) and upstroke (grey area) phases in a reference $\alpha_0 = 45^\circ$ case (top right) and kinematics of the flapping wing during the downstroke phase (bottom right). Grey and lightgrey colours in the bottom right plot are used for the rigid and morphing wings, respectively.</figcaption></figure>

wing morphing is here prescribed and does not result from a fluid-structure interaction problem.

In all cases, the Reynolds number based on the chord length $c$ and the mean velocity over a stroke $\bar{U}_{ref} = 1.221\text{m.s}^{-1}$ at the $r = 0.72R$ radial position is set to 840.

The flow is computed by directly solving the incompressible Navier-Stokes equations using StarCCM+ ver. 17.06. This solver uses a cell-centered finite volume method. An overset mesh technique (see Appendix A for further details) is employed where the wing is positioned within a small cylindrical domain, which undergoes the flapping motion within a larger, fixed background cylindrical domain. Figure 2 (left) displays the computational domain considered here. The wing is treated as a non-slip wall and flaps along an horizontal plane that is parallel to the top and bottom boundaries of the background cylindrical domain. The top boundary is modeled as a stagnation inlet while the bottom boundary is modeled as a pressure outlet. Lateral boundaries are treated as slip walls. Table I gives an overview of the boundary conditions used in the present study. Second-order schemes are employed for both spatial and temporal discretizations. A co-located variable arrangement and a Rhie-and-Chow-type pressure-velocity coupling combined with a SIMPLE-type algorithm are used.
<table><caption>TABLE I. Boundary conditions for the numerical domain considered in the present study.</caption><thead><tr><th>Boundary</th><th>Prescribed condition</th></tr></thead><tbody><tr><td>Wing</td><td>Non-slip wall</td></tr><tr><td>Top boundary</td><td>Stagnation inlet</td></tr><tr><td>Bottom boundary</td><td>Pressure outlet</td></tr><tr><td>Lateral boundaries</td><td>Slip walls</td></tr></tbody></table>

Low resolution cases are first considered to conduct the parametric study. These cases exhibit typical cell dimensions of $\Delta x_s = 0.02c$ and $\Delta x_v = 0.04c$ at the wing surface and in the small moving domain, respectively and a time-step of $\Delta t = T/250$. After identification of configurations of interest (i.e. parametric study), higher fidelity numerical simulations ($\Delta x_s = 0.005c$ and $\Delta t = T/1000$) are used to analyze and correlate aerodynamic performance to flow topology. Previous works on a similar setup$^{23,30,33}$ (i.e. similar geometry and solver) have shown that an increase in spatial and temporal resolutions from $\Delta x_s = 0.02c$ and $\Delta t = T/250$ do not yield significant changes in period-averaged lift coefficient $\bar{C}_L$. Furthermore, it has been shown that initial transients have sufficiently decayed after three flapping cycles. Nonetheless, a convergence study is performed here for the active twisting wing case. Similar conclusions can be drawn from the results, reported in Appendix B (Figure 20), where it is shown that, once initial transients have passed, decreasing $\Delta x_s$ from 0.02c to 0.01c yields changes in period-averaged lift $\bar{C}_L$ below 2.2%, and decreasing $\Delta x_s$ from 0.01c to 0.005c yields changes below 1.1% (similarly, changes in $\bar{C}_L$ as $\Delta t$ decreases from $T/250$ to $T/500$ and from $T/500$ to $T/1000$ are below 2% and 0.4% respectively).

Finally, it is worth mentioning that the solver has been previously validated on a variety of low Reynolds number flows such as that past axisymmetric bluff bodies$^{36}$, revolving wings$^{37}$ and perching airfoils$^{38}$.

In what follows, results are analyzed in terms of both mean and instantaneous lift coefficient $C_L$, power coefficient $C_P$, lift-to-power coefficient ratio $C_L/C_P$ and power associated to revolving motion coefficient $C_{P,rev}$. These coefficients are defined as follows:

$$C_L = 2L/(\rho \bar{S} \bar{U}_{ref}^2), \qquad (1)$$

$$C_P = -2P/(\rho \bar{S} \bar{U}_{ref}^3), \qquad (2)$$

$$C_{P,rev} = -2P_{rev}/(\rho \bar{S} \bar{U}_{ref}^3). \qquad (3)$$
<figure><img src="image_4.png" /><figcaption>FIG. 3. (a) Computational domain and definitions used to compute auxiliary potential fields (not at scale) and (b) example of the $\phi_{e_y}$ field in a cross-section of the wing.</figcaption></figure>

Where $S$ is the wing surface area, $\rho$ is the fluid density, and $L$, $P$ and $P_{rev}$ are the wing lift, global power and power associated to revolving motion respectively. The global power is computed taking into account the three revolving, pitching and twisting (for the morphing wing cases) torques. Mean values are obtained by time-averaging instantaneous values over the 3rd flapping period, following conclusions from Appendix B.

### III. FORCE AND POWER PARTITIONING METHOD

As previously discussed, vortices generated by the flapping motion of a wing play a crucial role in aerodynamic force production. The analysis of aerodynamic performance hence benefits from correlating aerodynamic loads and vortical structures. To achieve this goal, we employ a force-moment partitioning method (FMPM). This method is based on the work of Quartapelle and Napolitano$^{39}$, who demonstrated that by projecting the Navier-Stokes equations onto the gradient of a suitably constructed auxiliary potential field, the pressure-induced forces and moments acting on an immersed body can be expressed in terms of velocity-field gradients.

More specifically, we adopt the formulation derived by Menon and colleagues$^{40-43}$ where FMPM was applied to various 2D pitching and oscillating airfoil configurations. Note that this method has also been applied to 3D cases, for example by Zhang et al.$^{44}$ on data from a 3D flapping wing simulation, and by Zhu and Breuer$^{45}$ on PIV data from 3D pitching swept wings experiments. Here, in addition to force partitioning, we propose to partition the power with the aim of highlighting power-consuming structures in the flow past flapping wings.

The method is based on the resolution of auxiliary potentials $\phi_i$, $\psi_k$ with $i$ and $k$ the direction
Hovering flapping wings with dynamic twist

of the force and power to be partitioned, respectively. The auxiliary potentials for the force and
power are defined as:

$$
\begin{align}
\nabla^2 \phi_i &= 0 \text{ in } V_f, \quad \text{and} \nonumber \\
\hat{n} \cdot \nabla \phi_i &= \begin{cases} n_i, & \text{on the profile } B \\ 0, & \text{on the domain boundaries } \Sigma \end{cases} \tag{4}
\end{align}
$$

with $i = e_x, e_y, e_z$,

and:

$$
\begin{align}
\nabla^2 \psi_k &= 0 \text{ in } V_f, \quad \text{and} \notag \\
\hat{n} \cdot \nabla \psi_k &= \begin{cases} [(\vec{X} - \vec{X}_C) \times \Omega_k \hat{n}] \cdot \hat{e}_k, & \text{on } B \\ 0, & \text{on } \Sigma \end{cases} \tag{5}
\end{align}
$$

with $k = e_x, e_y, e_z$,

where $\vec{X}$ is the position vector of a point on the surface of the immersed body and $\vec{X}_C$ is the point about which the moment is computed. $\Omega_k$ is the angular velocity in the direction $k$. The velocity of the immersed body $\vec{U}_B$ is shown in Figure 3(a). Figure 3(b) displays an example of the $\phi_{e_y}$ field computed in a section of the fluid volume $V_f$.

As derived by Menon and Mittal⁴⁰, the instantaneous pressure-induced force on a wing can be
expressed as:

$$
\begin{equation}
\begin{split}
F_{p,i} = -\rho \underbrace{\iint_B \vec{n} \cdot \left( \frac{d\vec{U}_B}{dt} \right) \phi_i}_{I} dS - \rho \underbrace{\iiint_{V_f} 2Q\phi_i dV}_{II} \\
\qquad + \mu \underbrace{\iiint_{V_f} \vec{\nabla}^2 \vec{u} \cdot \vec{\nabla} \phi_i dV}_{III} - \rho \underbrace{\iint_{\Sigma} \vec{n} \cdot \left( \frac{d\vec{u}}{dt} \right) \phi_i}_{IV} dS
\end{split}
\tag{6}
\end{equation}
$$

, i = e_x, e_y, e_z.

Similarly, the instantaneous pressure-induced power can be expressed as:

$$
\begin{equation}
\begin{split}
P_{p,k} = -\rho \underbrace{\iint_B \vec{n} \cdot \left( \frac{d\vec{U}_B}{dt} \right) \psi_k}_{I} dS - \rho \underbrace{\iiint_{V_f} 2Q\psi_k dV}_{II} \\
\qquad + \mu \underbrace{\iiint_{V_f} \vec{\nabla}^2 \vec{u} \cdot \vec{\nabla} \psi_k dV}_{III} - \rho \underbrace{\iint_{\Sigma} \vec{n} \cdot \left( \frac{d\vec{u}}{dt} \right) \psi_k}_{IV} dS
\end{split}
\tag{7}
\end{equation}
$$

, k = e_x, e_y, e_z.
Hovering flapping wings with dynamic twist

We note that the auxiliary potential fields $\phi_i$ and $\psi_k$ only depends on the instantaneous shape, position and velocity of the wing during the flapping motion. We further note that these potentials are computed at each time-step.

The first term (I) in equations 6 and 7 is related to added-mass effects. The second term (II) is the vortex-related term, as evidenced by its relation to the Q-criterion. Recall that Q is the second-invariant of the velocity-gradient tensor, defined as :

$$Q = \frac{1}{2} (\|\Omega\|^2 - \|S\|^2), \qquad (8)$$

where $\Omega$ and $S$ are the anti-symmetric and symmetric parts of the velocity-gradient tensor respectively. $Q > 0$ identifies regions where rotation dominates over shear, referred to as vortices following the definition from Hunt et al.⁴⁶. The third term (III) in equations 6 and 7 is the viscous diffusion term. Finally, the fourth term (IV) is related to flow acceleration at the outer boundary of the computational domain and is negligible in our case (because of the large domain and the hovering flight condition).

The global force coefficient $C_F$ can be decomposed into two contributions, namely the pressure-induced $C_{F,p}$ and shear-induced $C_{F,s}$ contributions:

$$C_F = C_{F,p} + C_{F,s}. \qquad (9)$$

Using FPM, the global force coefficient $C_F$ can be expressed as a function of the contributions described in equation 6:

$$C_F = C_{F,\kappa} + C_{F,\omega} + C_{F,\sigma} + C_{F,s}, \qquad (10)$$

with $C_{F,\kappa}$ the added-mass-induced force coefficient, $C_{F,\omega}$ the vortex-induced force coefficient and $C_{F,\sigma}$ the viscous-induced force coefficient. Note that $C_{F,s}$ is negligible in most cases and is thus not shown in plots displaying quantities computed with the FPM. $C_{F,\kappa}$, $C_{F,\omega}$ and $C_{F,\sigma}$ are then the force coefficients related to term I, II and III presented in equation 6, respectively.

## IV. RESULTS

In this section, we report our results obtained from the parametric study of the influence of dynamic twist on aerodynamic performance of 3D flapping wings under hovering flight conditions.
Hovering flapping wings with dynamic twist

<figure><img src="image_5.png" /><figcaption>FIG. 4. $\bar{C}_L$ as a function of $\bar{C}_L/\bar{C}_P$ for different pitch ($\alpha_0$) and twist ($\varepsilon_0$) angles. Data obtained for positive and negative $\varepsilon_0$ are shown in (a) and (b), respectively. White squares show rigid wing performance. Light-grey, dark grey and black squares show performance obtained with $\varepsilon_0 = \pm 15^\circ, 30^\circ$ and $45^\circ$, respectively.</figcaption></figure>

A. Overall performance

Figure 4 shows the mean lift coefficient $\bar{C}_L$ as a function of the mean lift coefficient to mean power coefficient ratio $\bar{C}_L/\bar{C}_P$ obtained for different pitch ($\alpha_0$) and twist ($\varepsilon_0$) angles (similar plots showing $\bar{C}_L$ as a function of $\bar{C}_P$ are provided in Appendix D). White squares depict cases without twist (i.e. rigid wing) and with $\alpha_0$ varying from 0° to 90° (5° step between 0° to 75°). Recall that $\alpha_0 = 0°$ corresponds to the wing being horizontal at mid-downstroke (and mid-upstroke). Similarly, $\alpha_0 = 90°$ corresponds to the wing being vertical during the whole simulation (no pitch motion). Light grey, dark grey and black squares depict twisted wing cases with $\varepsilon_0 = \pm 15°, 30°$ and $45°$, respectively. Cases with positive $\varepsilon_0$ are shown in Figure 4 (a) while cases with negative $\varepsilon_0$ are shown in Figure 4 (b) to ease the reading. Plain black lines are used to connect cases without twist and dashed lines are used to connect cases with similar $\alpha_0$.

First, for cases without twist, it is observed that $\bar{C}_L$ increases with $\alpha_0$ from $\alpha_0 = 0°$ to $45°$ and then drops as $\alpha_0$ further increases to $90°$. On the other hand, $\bar{C}_L/\bar{C}_P$ increases with $\alpha_0$ from $\alpha_0 = 0°$ to $15°$ and then drops as $\alpha_0$ further increases. That is, $\bar{C}_L/\bar{C}_P$ reaches its maximum (around 1.15 at $15°$) well before that of $\bar{C}_L$ (around 1.2 at $45°$). The maximum lift pitch angle is consistent with previous data from the literature$^{2,33}$. Although the effective angle-of-attack of the wing is lower than $45°$ due to the downwash induced as a reaction to vertical force production during previous strokes, stall (in the time-averaged sense) is delayed to larger pitch angles than those conventionally observed on 2D wings. According to previous studies$^{6,8}$, the attachment of
Hovering flapping wings with dynamic twist

the leading edge vortex (LEV) which forms on the suction side of the wing is promoted by strong
rotational effects. This LEV, comparable to that observed on delta wings, supports a low pressure
region that contributes to generating a strong, normal aerodynamic force. The strength of the
LEV increases with α₀, hence the normal force. However, C̄_L decreases at some point, here after
α₀ = 45°, because the relative contribution of the normal force to the vertical force decreases with
α₀.

Second, it is shown that a positive twist (Figure 4 (a)), when applied in-phase with the pitching motion, does not benefit aerodynamic performance because none of the twisted wing cases outperforms rigid wing cases in either C̄_L or C̄_L/C̄_P. In other words, positive wing twist does not help in reaching the portion of the (C̄_L, C̄_L/C̄_P) space located on the right and above the plain lines that connect α₀ = 15° to 45° (indicated by the shaded region). Nonetheless, positive twist can significantly affect C̄_L and C̄_L/C̄_P and can hence be used to adapt aerodynamic performance during flight. For example, it can be seen that twisting the α₀ = 30° wing with ε₀ = 15° yields comparable changes in aerodynamic performance to those obtained by pitching the whole wing to 45°. The overall increase in lift and reduction in lift-to-power coefficient ratio observed as ε₀ increases corroborates previous observations²⁴,²⁵,³¹ and can here be correlated with a general increase in angle of attack along the span and hence an increase in LEV strength, which results in a larger normal force on the wing. However, as previously mentioned, the increase in lift is limited since, concomitantly, the projection of the normal force on the vertical axis decreases with ε₀.

Conversely, it is shown that a negative twist (Figure 4 (b)), again when applied in-phase with
the pitching motion, helps increase aerodynamic performance when compared to rigid wing cases.
Maximum C̄_L/C̄_P is here obtained for α₀ = 45° and ε₀ = -45° and reaches 1.22 i.e. 5.07% larger
than that obtained for α₀ = 15° and ε₀ = 0°. In addition to the α₀ = 45° and ε₀ = -45° case, two
other twisted wing cases lie within the shaded region in Figure 4 (b) (i.e. improving performance
with respect to rigid wing cases, ε₀ = 0°): the α₀ = 45° and ε₀ = -30° case, and the α₀ = 30°
and ε₀ = -15° case. More generally, and similarly to that observed for positive twist, negative
twist can significantly affect aerodynamic performance and can thus help in adjusting lift during
flight. Here again, it is observed, for example, that twisting the α₀ = 45° wing with ε₀ = -15°
yields comparable changes in aerodynamic performance to those obtained by pitching the whole
wing to 35°. Similarly to what has been discussed for rigid wings and wings with positive twist,
changes in lift and lift-to-power coefficient ratio with decreasing, negative twist angle ε₀ can be
correlated with changes in LEV strength, hence normal force, and the increasing contribution of
Hovering flapping wings with dynamic twist

<figure><img src="image_6.png" /><figcaption>FIG. 5. Aerodynamic performance as a function of mid-stroke pitch angle at the $r = 2/3R$ radial position $\alpha_{2/3}$ for positive $\varepsilon_0$ cases. $\bar{C}_L$ and $\bar{C}_L/\bar{C}_P$ as a function of $\alpha_{2/3}$ are shown in (a) and (b), respectively. Symbols are similar to those used in Figure 4.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 6. Same plots as in Figure 5 but for negative $\varepsilon_0$ cases.</figcaption></figure>

the projection of the normal force on the vertical axis.

These trends will be further discussed in the next sections following two perspectives. First, it will be shown that data are driven, to leading order, by a characteristic angle of attack along the span, providing further evidence for the correlation between wing twist and normal force production. Second, relevant cases will be discussed in detail using the force and power partitioning method to further highlight the correlation between vortical structures and force production.

## B. Scaling

Figures 5 (a) and (b) present $\bar{C}_L$ and $\bar{C}_L/\bar{C}_P$ as functions of the mid-stroke pitch angle at the $r = 2/3R$ radial position, $\alpha_{2/3}$ (for rigid wings, $\alpha_{2/3} = \alpha_0$), for positive $\varepsilon_0$ values, respectively.
Hovering flapping wings with dynamic twist

Note that this radial position approximately corresponds to that where $\int_0^r Ldr = \int_r^R Ldr$, as verified at mid-downtroke for the reference, rigid wing case with $\alpha_0 = 45°$ (not shown here for the sake of conciseness). The same colors and symbols as in previous plots are used for consistency with the previous subsection, i.e. rigid wing performance is indicated using white squares, while light grey, dark grey, and black squares depict positive values $\epsilon_0 = 15°, 30°$ and $45°$, respectively.

It is observed from Figure 5 (a) that data generally collapse onto a single curve, with some exceptions, typically for $\alpha_0 = 15°$ and $\epsilon_0 = 45°$, $\alpha_0 = 30°$ and $\epsilon_0 = 30°$, and $\alpha_0 = 30°$ and $\epsilon_0 = 45°$. In these cases, the pitch angle passes from below to above $45°$ along the span. That is, an increase in pitch angle will promote lift coefficient in inboard sections while being detrimental in outboard sections (see Figure 4). The evolution of wing lift with pitch angle is thus different to that obtained for cases where the pitch angle is either always below or always above $45°$ along the span. A similar observation can be made from Figure 5 (b), where cases with similar $\alpha_{2/3}$ generally yield similar $\overline{C}_L/\overline{C}_P$. Again, this is true for most cases except for $\alpha_0 = 0°$ and $\epsilon_0 = 30°$, and $\alpha_0 = 0°$ and $\epsilon_0 = 45°$, for example. In these cases, the pitch angle passes from below to above $15°$ along the span. That is, an increase in pitch angle promotes efficiency in inboard sections while being detrimental in outboard sections (see Figure 4). The evolution of wing efficiency with pitch angle is thus different to that obtained for cases where the pitch angle is either always below or always above $15°$ along the span. Finally, it is again visible from those plots that wings with positive twist applied in-phase with the pitching motion never outperform rigid wings, both in $\overline{C}_L$ and $\overline{C}_L/\overline{C}_P$.

Similar to Figure 5, Figure 6 (a) and (b) show $\overline{C}_L$ and $\overline{C}_L/\overline{C}_P$ as functions of the mid-stroke pitch angle at the $r = 2/3R$ radial position $\alpha_{2/3}$ for negative $\epsilon_0$ values, respectively. Again, most data are found to collapse onto a single curve, with some exceptions. In Figure 6 (a) these include, for example, cases with $\alpha_0 = 90°$ and $\epsilon_0 = -45°$, and $\alpha_0 = 75°$ and $\epsilon_0 = -45°$, where the pitch angle along the span passes from above to below $45°$. In Figure 6 (b), there is no clear outliers. Only small data scattering is observed near maximum $\overline{C}_L/\overline{C}_P$ (around $\alpha_{2/3} = 15°$) where cases with negative twist are found to outperform rigid wing cases.

The physical reasons behind these changes in aerodynamic performance with $\alpha_0$ and $\epsilon_0$ will be discussed in the next subsections.
<figure><img src="image_8.png" /><figcaption>FIG. 7. Surface-integral and FPM-computed $C_L$ as a function of non-dimensional time $t/T$ for the rigid wing case with $\alpha_0 = 45^\circ$. Added-mass $C_{L,\kappa}$, vortex-induced $C_{L,\omega}$ and viscous $C_{L,\sigma}$ contributions to lift from FPM are also shown.</figcaption></figure>

### C. Analysis at maximum $\bar{C}_L$

We first focus on the flapping wing case for which maximum $\bar{C}_L$ is obtained, i.e. rigid wing with mid-stroke pitch angle $\alpha_0 = 45^\circ$. The analysis relies on the time history of lift and power coefficients, correlated with flow topology using the force partioning method (FPM) described in section III.

#### 1. Lift and lift-producing structures

Figure 7 shows the instantaneous lift coefficient $C_L$ as a function of the non-dimensional time $t/T$ obtained during the fourth downstroke phase. Both $C_L$ computed by integrating pressure and shear at the wing surface and from FPM are shown, using plus symbols and plain line respectively.

First, it is verified that the lift computed using FPM is similar to that obtained through wing surface integration. We thus leverage the use of FPM to decompose lift into different contributions associated with various terms in equation (6). Accordingly, added-mass, vortex and viscous-induced lift are displayed in Figure 7 using dotted, dashed and dashdotted lines respectively. Recall that the shear stress contribution is not displayed because it is small compared to other contributions. The lift induced by added-mass, denoted as $C_{L,\kappa}$, exhibits a sinusoidal behavior in line with the sinusoidal kinematics of the wing. $C_{L,\kappa}$ is approximately zero at the start of the motion and then increases until reaching its maximum value at $t/T = 3.12$. Following this peak, $C_{L,\kappa}$ decreases to negative values, with the minimum observed at $t/T = 3.40$. Eventually, $C_{L,\kappa}$ increases again
Hovering flapping wings with dynamic twist

until reaching approximately zero at the end of the downstroke phase. The lift induced by viscous
diffusion, denoted as $C_{L,\sigma}$, is of comparable magnitude to that of $C_{L,\kappa}$. It is also null at the start
of the motion and then increases to reach its maximum around midstroke before decreasing back
to zero at the end of the stroke. Conversely, the vortex-induced lift, denoted as $C_{L,\omega}$, has a much
larger magnitude than that of the added-mass and viscous diffusion terms and largely dominates
lift production. It exhibits a nearly sinusoidal shape, in phase with the viscous diffusion lift. That
is, it globally increases from 0 at $t/T = 3.00$ to its maximum at $t/T = 3.28$, and decreases back to
zero at $t/T = 3.50$. The evolution during the first half-stroke is however not strictly monotonous,
with a small kink observed around $t/T = 3.12$.

Because $C_{L,\omega}$ is found to dominate lift production, we correlate the evolution of lift with flow
topology by plotting a time-sequence of $Q$-criterion isosurfaces $\frac{Qc^2}{\overline{U}_{ref}^2} = 10$ colored by the
amplitude of the vortex lift term of equation (6) ($-2\rho Q\phi_{ey}$) in Figure 8. That is, we show vortex
structures in the flow and how they locally contribute to lift.

At $t/T = 3.00$, the flow is characterized by small-scale structures shed in the wake during
the previous upstroke phase. As will be shown thereafter, these are the traces of the leading edge
vortex (LEV) and tip vortex (TV) that formed during previous instants. Their white color indicates
that their contribution to lift, at this instant, is negligible. This is in line with the zero lift force
observed at $t/T = 3.00$ on Figure 7. From this instant, the wing starts to revolve about its root
from left to right and thus interacts with wake structures. By $t/T = 3.05$, the detrimental effect of
this interaction becomes apparent as wake structures sweeping on the lower surface of the wing
support negative values of $-2\rho Q\phi_{ey}$ (blue contours). Moving to $t/T = 3.10$, this detrimental effect
intensifies due to wing orientation which contributes to increasing $\phi_{ey}$ on the lower surface (the y
projection of the normal to the surface increases in magnitude). That is, relatively intense vortices
(with large $Q$ values) interact with regions of relatively large positive $\phi_{ey}$. On the other hand, a
new LEV starts forming on the wing's upper surface, contributing to lift production (as shown
by red contours). Similarly to the mechanism described above, the LEV is a region of large $Q$
values interacting with regions of large negative $\phi_{ey}$. The LEV is the main contributor to lift and
is responsible for the increase in lift observed on Figure 7. However, it is slightly counterbalanced
by wing-wake interactions such that, although the LEV continues to grow between $t/T = 3.10$
and $t/T = 3.15$, a kink in the lift curve is observed near $t/T = 3.12$. At $t/T = 3.15$ both red
and blue contours are observed on the upper and lower surface of the wing, respectively. At
$t/T = 3.20$, the newly formed LEV and a tip vortex (TV) are clearly visible, with the TV seemingly
<figure><img src="image_9.png" /><figcaption>FIG. 8. Q-criterion isosurfaces $Qc^2/\bar{U}_{ref}^2 = 10$ colored by the contribution to the vortex-induced lift ($-2\rho Q\phi_{ey}$) for the rigid wing case with $\alpha_0 = 45^\circ$. The wing surface is shown in grey and the wing outline is shown using dotted lines on the snapshot at $t/T = 3.05$.</figcaption></figure>
Hovering flapping wings with dynamic twist

having minimal influence on lift production. At mid-downstroke ($t/T = 3.25$), it is observed
that both structures connect. The connection, or merging, of the LEV and TV is known to be
associated with the interaction of the LEV core flow (directed outboard) with the tip vortex flow
(directed inboard)$^{47}$, and the subsequent bursting of these structures into small scale structures
near the tip. These small structures are observed to advect in the wake of the wing. This instant
coincides with maximum lift production. Interestingly, the unbursted, conical part of the LEV
inboard supports larger $-2\rho Q\phi_{e_y}$ values than its bursted part near the tip. This provides evidence
that LEV burst (due to its interaction with the TV) is detrimental to lift production (i.e. with respect
to the coherent, unbursted region of the LEV). By $t/T = 3.30$, as the wing begins decelerating,
the LEV is found to burst over a wider spanwise region. In other words, deceleration contributes
to LEV burst and to the extent of the bursted region inboard. The LEV continues to grow due
to the revolving motion but the conical portion that contributes the most to lift is reduced. In
addition, the wing starts to pitch up (back to 90° at $t/T = 3.50$) and hence the y projection of
the normal to the surface decreases in magnitude, reducing the magnitude of $\phi_{e_y}$. The lift thus
starts to decrease (see Figure 7). At $t/T = 3.35$, the pitch angle is large and the wing continues
to decelerate, which contributes to further decreasing lift. Moving to $t/T = 3.40$, LEV burst is
observed over a significant portion of the wing. Because of the large pitch angle, $\phi_{e_y}$ is low and
the contribution of the LEV to lift drastically decreases, in line with the evolution of $C_L$. Finally,
at $t/T = 3.45$, as the wing approaches stroke reversal, the LEV has a very weak contribution to
lift, as evidenced by its white contours.

In further elucidating the relationship between the structure of the LEV and global lift produc-
tion, Figure 9 shows the pressure lift coefficient $C_{L,p}$ as a function of the non-dimensional span-
wise location $r/R$ and non-dimensional time $t/T$ for the rigid wing case ($\epsilon_0=0^\circ$) with $\alpha_0 = 45^\circ$.
Positive and negative contours are depicted by plain and dashed lines, respectively.

During early stages, $C_{L,p}$ increases with $t/T$ for all spanwise sections. Near the tip, this in-
crease rapidly saturates and $C_{L,p}$ then exhibits an unsteady behavior. The first evidence of a local
maximum in the $C_{L,p}$ map as $t/T$ increases, which marks the onset of unsteady flow and sectional
lift oscillations, is observed at $r/R = 0.975$ (indicated with the white dot). The spanwise portion of
the wing that supports an unsteady flow then extends inboard with time. For example, local max-
ima in $C_{L,p}$ are observed at $r/R = 0.87, 0.79$ and $0.715$ at $t/T = 3.16, 3.21$ and $3.252$, respectively.
Those locations are indicated using white stars. This is in line with the occurrence and the extent of
LEV burst observed previously on Figure 8. For spanwise locations lower than $r/R = 0.70$, there is
Hovering flapping wings with dynamic twist

<figure><img src="image_10.png" /><figcaption>FIG. 9. Pressure-induced lift coefficient $C_{L,p}$ as a function of non-dimensional spanwise location $r/R$ and non-dimensional time $t/T$ for the rigid wing case with $\alpha_0 = 45^\circ$. Plain and dashed lines display positive and negative $C_{L,p}$ isolines respectively. White dot and stars indicate first maxima in $C_{L,p}$ obtained in four cross sections $r/R = 0.975, 0.87, 0.79$ and $0.715$.</figcaption></figure>

no evidence of lift oscillations. In this inboard region, the flow reaches a quasi-steady state, which is in line with the quasi-steady, conical LEV shape observed on Figure 8. With this representation, it is visible that, at a given instant, the LEV has maximum contribution to lift immediately before bursting (as $r/R$ increases) and that bursting is detrimental to lift production, again corroborating previous results from Figure 8.

## 2. Power and power-consuming structures

Figure 10 displays the time-history of total and revolving power coefficients, $C_p$ and $C_{P,rev}$ respectively, obtained by integrating pressure and viscous stresses at the wing surface. It is shown that the revolving power accounts for most of the total power and that pitching power is comparatively negligible. Therefore, we focus our attention on the revolving power and further plot the time-history of $C_{P,rev}$ computed from the power partitioning method (PPM). Again, good agreement is observed between both $C_{P,rev}$ (obtained from wing surface integration and PPM) and we thus leverage the use of PPM to decompose revolving power into different contributions associated with various terms in equation 7. Specifically, revolving power arising from added-mass, vortex and viscous terms are displayed in Figure 10 using dotted, dashed and dashdotted lines
<figure><img src="image_11.png" /><figcaption>FIG. 10. Surface-integral and PPM-computed $C_{P_{rev}}$ as a function of non-dimensional time $t/T$ for the rigid wing case with $\alpha_0 = 45^\circ$. Surface-integral $C_P$ as well as added-mass $C_{P_{rev},\kappa}$, vortex-induced $C_{P_{rev},\omega}$ and viscous $C_{P_{rev},\sigma}$ contributions to revolving power from PPM are also shown.</figcaption></figure>

respectively. The time evolution of the different contributions mirror those observed for the lift. That is, the added-mass-induced power $C_{P_{rev},\kappa}$ exhibits a sinusoidal shape in phase with the pitching motion. Similarly, the viscous-induced power $C_{P_{rev},\sigma}$ also manifests a sinusoidal shape, albeit with a phase shift relative to $C_{P_{rev},\kappa}$. Again, $C_{P_{rev},\kappa}$ and $C_{P_{rev},\sigma}$ have comparable, relatively weak magnitudes. Conversely, the vortex-induced power $C_{P_{rev},\omega}$ dominates power consumption and dictates the overall trend in the time-evolution of $C_{P_{rev}}$ for most of the downstroke phase.

Again, because $C_{P_{rev},\omega}$ is found to dominate power consumption, we correlate the evolution of power with flow topology by plotting a time-sequence of $Q$-criterion isosurfaces $Qc^2/\overline{U}_{ref}^2 = 10$ colored by $-2\rho Q\psi_{ey}$ in Figure 11. That is, we show vortex structures in the flow and how they locally contribute to power consumption. We do not show the whole sequence displayed in Figure 8 but rather focus on 4 relevant instants.

At $t/T = 3.00$, the vortical structures have negligible impact on power consumption, as indicated by their white contours. At this instant, the revolving speed is null and hence $\psi_{ey}$ is weak (but not null because of the pitch velocity). At $t/T = 3.10$, the development of the new LEV on the suction side of the wing is identified as a power consuming structure (red contours), while strong interactions between the wing and its own wake are characterized by power recovery (blue contours). At $t/T = 3.20$, the LEV continues to grow both in strength and size, resulting in the interaction between large $Q$ values and relatively large positive $\psi_{ey}$ values. It is observed that $-2\rho Q\psi_{ey}$ slightly drops near the very tip of the wing but that there is no clear drop with LEV burst, contrasting previous observations from Figure 8. Here, $\psi_{ey}$ increases with the radial loca-
Hovering flapping wings with dynamic twist

<figure><img src="image_12.png" /><figcaption>FIG. 11. Q-criterion isosurfaces $Qc^2/\bar{U}_{ref}^2 = 10$ colored by the contribution to the vortex-induced revolving power ($-2\rho Q\psi_{ey}$) for the rigid wing case with $\alpha_0 = 45^\circ$.</figcaption></figure>

tion $r$ (see Figure 22) which favors power consuming structures outboard. Accordingly, power consumption continues to increase, as previously discussed from Figure 10. By $t/T = 3.30$, the burst has expanded inboard. It is again clearly visible that the bursted portion of the LEV has a significant contribution to power consumption. Because it was shown that, on the other hand, it is not conducive to high lift production, this suggests that LEV burst is detrimental to aerodynamic efficiency.

### D. Analysis at maximum $\bar{C}_L/\bar{C}_P$

We now focus on the flapping motion for which maximum $\bar{C}_L/\bar{C}_P$ is obtained, i.e. twisted wing with mid-stroke pitch angle $\alpha_0 = 45^\circ$ and twist angle $\epsilon_0 = -45^\circ$. This configuration is also compared with the maximum $\bar{C}_L/\bar{C}_P$ obtained for a rigid wing, i.e. wing with mid-stroke pitch angle $\alpha_0 = 15^\circ$. The analysis relies on the time history of lift and power coefficients, correlated with flow topology using the force partitioning method described in section III and already used in the previous subsection.
Hovering flapping wings with dynamic twist

<figure><img src="image_13.png" /><figcaption>FIG. 12. Surface-integral and FPM-computed $C_L$ as a function of non-dimensional time $t/T$ for the rigid wing case with $\alpha_0 = 15^\circ$ (a) and the twisted wing case with $\alpha_0 = 45^\circ$ and $\varepsilon_0 = -45^\circ$ (b). Added-mass $C_{L,\kappa}$, vortex-induced $C_{L,\omega}$ and viscous $C_{L,\sigma}$ contributions to lift from FPM are also shown.</figcaption></figure>

# 1. Unsteady lift and power

Figures 12 (a) and 12 (b) show the time history of $C_L$ computed from surface integration and FPM, as well as its added-mass $C_{L,\kappa}$, vortex $C_{L,\omega}$ and viscous diffusion $C_{L,\sigma}$ contributions obtained from FPM for the rigid $\alpha_0 = 15^\circ$ and twisted $\alpha_0 = 45^\circ$, $\varepsilon_0 = -45^\circ$ wing cases, respectively.

We first note that both cases exhibit very similar $C_{L,\kappa}$ and $C_{L,\sigma}$ contributions, and that differences in $C_L$ mainly arises from differences in vortex-induced lift $C_{L,\omega}$.

During initial stages, from $t/T = 2.00$ to $t/T = 2.05$, $C_L$ increases in both cases due to added-mass effects ($C_{L,\kappa}$, as further discussed thereafter). Then, $C_L$ decreases until $t/T = 2.12$. At this instant, the rigid wing experiences negative lift, while the twisted wing case exhibits positive lift. This time interval (from $t/T = 2.05$ to $t/T = 2.12$, i.e immediately following stroke reversal) typically involves interactions between the wing and the wake from the previous stroke. Hence, differences in $C_L$ between rigid and twisted wing cases at those instants may result from differences in wing-wake interactions. This point will be discussed on the basis of the evolution of $C_{L,\omega}$ and of the flow topology presented in the next subsection. After $t/T = 2.12$, $C_L$ increases again for both cases. This increase can again be associated with (i) the development of a LEV on the upper surface of the wing as previously discussed for the $\alpha_0 = 45^\circ$ case and/or (ii) a reduced negative effect of the interaction of the wing with its wake on the lower surface. The increase is stronger in the rigid wing case such that the corresponding lift becomes larger than that of the twisted wing case near $t/T = 2.22$. $C_L$ being almost exclusively driven by $C_{L,\omega}$ at this time, this indicates
Hovering flapping wings with dynamic twist

a larger LEV on the suction side of the rigid wing with respect to that of the twisted wing case. Maxima in $C_L$ are reached around $t/T = 2.35$ for both cases, with the rigid wing experiencing 14% more lift than the twisted wing case. Finally, $C_L$ rapidly drops to zero from the time of maximum lift to the end of the stroke ($t/T = 2.50$).

As already mentioned, the increase in $C_L$ at the beginning of the stroke is driven by added-mass effects for both rigid and morphing wing cases. Accordingly, $C_{L,\kappa}$ increases from $t/T = 2.00$ to $t/T = 2.08$. This increase is due to the combination of the pitch-down motion about the quarter-chord and revolving acceleration. Specifically, revolving acceleration is maximum at $t/T = 0$ and the corresponding normal force is thus maximum, but its projection on the vertical axis is around zero. As $t/T$ increases, the revolving acceleration decreases but the projection of the corresponding added-mass force on the vertical axis increases, leading to the increase in $C_{L,\kappa}$. The added-mass induced lift then decreases and becomes negative slightly after $t/T = 2.25$, where the revolving acceleration becomes negative. It reaches minimum value at $t/T \approx 2.42$ and increases again until reaching zero (no lift production) at the end of the stroke, where the revolving deceleration is maximum but the projection of the corresponding normal force on the vertical axis is zero.

$C_{L,\omega}$ is also zero at $t/T = 2.00$ because of the orientation of the wing (i.e. projection of the vortex-induced normal force on the vertical axis is zero). The vortex-induced lift then drops for both cases, from $t/T = 2.00$ to $t/T \approx 2.12$. This drop is related to the wing experiencing interactions with its own wake, resulting in negative vortex-lift production. The negative influence of wake interactions on lift production have been observed on Figure 8. On the other hand, a new LEV starts forming on the upper side of the wing. However, this new LEV is not strong enough yet to counterbalance the negative impact of wing-wake interactions. Interestingly, the actively twisted wing experiences larger $C_{L,\omega}$, yet still negative. This indicates a different flow topology between both configurations which will be detailed later. After $t/T \approx 2.12$, both $C_{L,\omega}$ increase, indicating the development of a strong LEV on the wing's suction side. Moreover, as observed for the $\alpha_0 = 45^\circ$ case, after this time, the wing has passed over most of the vortical structures formed during the previous stroke. In other words, the influence of wing-wake interactions decreases as the wing revolves. The maximum in $C_{L,\omega}$ is reached around $t/T = 2.38$ for both configurations. Vortex-induced lift then drops as the wing approaches stroke-reversal because the leading edge vortex structure loses coherency and the wing pitches up such that the projection of the corresponding normal force on the vertical axis decreases.
Hovering flapping wings with dynamic twist

<figure><img src="image_14.png" /><figcaption>FIG. 13. Surface-integral $C_L$ and $C_p$ (a) and $C_L/C_p$ (b) as a function of non-dimensional time $t/T$ for the rigid wing case with $\alpha_0 = 15^\circ$ and the twisted wing case with $\alpha_0 = 45^\circ$, $\epsilon_0 = -45^\circ$.</figcaption></figure>

Finally, $C_{L,\sigma}$ slowly increases from $t/T = 2.00$ to 2.37. The viscous diffusion induced lift is correlated with the formation of the LEV, that is with $C_{L,\omega}$ between $t/T = 2.20$ and 2.50. Differences between both $C_{L,\sigma}$ from rigid and twisted wing cases are however weak and do not influence the overall differences in $C_L$. Again, $C_{L,\sigma}$ is around zero at the beginning and end of the stroke because of the corresponding normal force projection on the vertical axis.

Overall, the lower lift observed near $t/T = 2.12$ in the rigid wing case slightly compensates for the larger lift observed near $t/T = 2.35$, such that differences in time-averaged lift between both rigid and twisted wing cases remain relatively weak, as previously observed on Figure 4. On the contrary, differences in time-averaged lift-to-power ratio are substantial and we thus now analyze both unsteady lift and power coefficients in Figure 13 (a) and unsteady lift-to-power coefficients ratio in Figure 13 (b).

From Figure 13 (a), it is observed that $C_p$ values from both rigid and twisted wing cases are similar from $t/T = 2.00$ to $t/T = 2.25$. During this interval, they both smoothly increase from approximately 0 to 0.5. Conversely, as previously shown, $C_L$ from the twisted wing case exceeds that of the rigid wing case, which leads to the significantly larger $C_L/C_p$ ratio displayed in Figure 13 (b). After $t/T = 2.25$, $C_p$ from both cases increase more strongly until reaching a maximum near $t/T = 2.38$ and then drop to approximately zero at the end of the stroke. The increase is however more pronounced in the rigid wing case such that the corresponding $C_p$ is significantly larger than that of the twisted wing case during this time interval. On the other hand, similar differences are observed on $C_L$, such that $C_L/C_p$ from both cases are relatively similar, see Figure 13 (b).
Hovering flapping wings with dynamic twist

These observations demonstrate that the time interval between $t/T = 2.05$ and $t/T = 2.20$ is crucial in explaining the larger $\bar{C}_L/\bar{C}_P$ obtained for the twisted wing case compared to the rigid wing case. Specifically, while $C_P$ from both configurations exhibit similar trends and values during this interval, the larger $C_L$ from the twisted wing, compared to the rigid one, is responsible for the larger $\bar{C}_L/\bar{C}_P$. As previously discussed, differences in lift between both cases arise from differences in vortex-induced lift. We thus focus on the analysis of vortex-induced lift in the next subsection.

## 2. Lift producing structures

The left column of Figure 14 shows a time sequence of $Q$-criterion isosurfaces $Qc^2/\bar{U}_{ref}^2 = 10$ colored by vortex-induced lift ($-2\rho Q\phi_{ey}$), obtained for the rigid wing case with $\alpha_0 = 15^\circ$.

At $t/T = 2.00$, the flow is characterized by the LEV and TV from the previous stroke. A root vortex (RV) is also visible near the wing root. Their white contour indicates that they have negligible contribution to lift at this instant. As already discussed for the rigid wing case with $\alpha_0 = 45^\circ$, this is due to the orientation of the wing and the resulting weak $\phi_{ey}$ field. At $t/T = 2.10$, the LEV from the previous stroke clearly appears as a source of negative lift (blue contours) as it sweeps on the lower surface of the wing. Similarly, the TV and the RV also act as sources of negative lift near the wing tip and root, respectively. On the other hand, the development of a new LEV, extending from midspan to the wing tip (i.e., from $r/R = 1/2$ to $r/R = 1$) on the upper surface of the wing, counteracts those negative contributions. This is further highlighted on Figure 15 which shows the corresponding pressure-induced lift coefficient $C_{L,p}$ as a function of non-dimensional spanwise location $r/R$ and non-dimensional time $t/T$. At this instant ($t/T = 2.10$), negative $C_{L,p}$ is observed inboard, highlighting the interaction between the wing and the LEV from the previous stroke. On the contrary, positive $C_{L,p}$ is observed outboard, highlighting the role of the new LEV in counteracting this negative contribution. Negative and positive contributions from inboard and outboard regions cancel out such that $C_L$ is equal to 0, see Figure 13. At $t/T = 2.15$, as the wing continues to revolve, the new LEV grows and the orientation of the wing favors large $\phi_{ey}$ values, hence interacting with relatively large $Q$ values. Accordingly, $C_{L,p}$ increases in the outboard region, see Figure 15. The negative $C_{L,p}$ region observed inboard is still significant, again highlighting strong interactions between the wing and the LEV from the previous stroke. The following snapshots on Figure 14, at $t/T = 2.175$ and $t/T = 2.20$, show similar flow topology.
<figure><img src="image_15.png" /><figcaption>FIG. 14. Q-criterion isosurfaces $Qc^2/\bar{U}_{ref}^2 = 10$ colored by the contribution to the vortex-induced lift $-2\rho Q\phi_{ey}$ for the rigid wing case with $\alpha_0 = 15^\circ$ (left column) and the twisted wing case with $\alpha_0 = 45^\circ$, $\epsilon_0 = -45^\circ$ (right column).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>FIG. 15. Pressure-induced lift coefficient $C_{L,p}$ as a function of non-dimensional spanwise location $r/R$ and non-dimensional time $t/T$ for the rigid wing case with $\alpha_0 = 15^\circ$. Plain and dashed lines display positive and negative isolines respectively.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIG. 16. Pressure-induced lift coefficient $C_{L,p}$ as a function of non-dimensional spanwise location $r/R$ and non-dimensional time $t/T$ for the twisted $\alpha_0 = 45^\circ, \epsilon_0 = -45^\circ$ case. Plain and dashed lines display positive and negative isolines respectively.</figcaption></figure>

Specifically, the LEV continues to expand inboard and strengthen outboard, resulting in larger $C_{L,\omega}$, and consequently, larger $C_{L,p}$ (Figure 15) and $C_L$ (Figure 13). The wing has passed over the structures from the previous stroke and hence the impact of wing wake interactions on lift has decayed.

The right column of Figure 14 shows similar snapshots of Q-criterion isosurfaces colored by
Hovering flapping wings with dynamic twist

$-2\rho Q\phi_{e_y}$, obtained for the twisted wing case with $\alpha_0 = 45^\circ$ and $\epsilon_0 = -45^\circ$.

At $t/T = 2.00$, the flow is relatively similar to that observed for the rigid wing case, with the presence of a LEV and a TV that formed during the previous stroke. However, there is no clear evidence of a RV. At $t/T = 2.10$, as the wing revolves and pitches down, the LEV from the previous stroke acts as a source of negative lift and a new LEV forms on the upper side of the wing, acting as a source of positive lift. This is again very similar to what has been observed for the rigid wing case. However, slight differences exist between both cases. In addition to the absence of a RV, the LEV in the twisted wing case extends slightly more inboard than that in the rigid wing case. This is visible from the close-up view at $t/T = 2.15$ and 2.175. As previously addressed in section IV A, the LEV supports a low pressure region on the wing upper surface which contributes to generating a strong, normal aerodynamic force pointing upward (i.e. generating lift). That is, the extended leading-edge vortex in the twisted wing configuration lead to enhanced, instantaneous lift with respect to the rigid wing case. Similarly, the presence of the RV in the rigid wing case is associated with a low pressure region near the wing lower surface which induces a normal aerodynamic force pointing downward. This results in the negative contribution to lift observed near the wing root in Figure 15 for the rigid wing case. As a consequence, for the twisted wing configuration, there is no negative contribution from RV and more positive contribution from LEV to lift in the inboard region, which leads to a reduced negative $C_{L,p}$ region in Figure 16. This reduced negative lift region explains the difference in $C_L$ observed between both cases on Figure 13 (a), in the time interval $t/T \in [2.05; 2.20]$, and hence the difference in $C_L/C_P$ observed on Figure 13 (b) and the difference in $\bar{C}_L/\bar{C}_P$ observed on Figure 4.

We stress that this reasoning also relies on the observation that power coefficients from both cases are similar during this time interval, see Figure 13 (a). In the next subsection, we thus briefly analyze power consumption to provide physical explanations for this similarity.

### 3. Power consuming structures

Figures 17 (a) and 17 (b) show the time history of $C_{P_{rev}}$ computed from surface integration and PPM, as well as added-mass $C_{P_{rev},\kappa}$, vortex $C_{P_{rev},\omega}$ and viscous diffusion $C_{P_{rev},\sigma}$ contributions obtained from PPM for the rigid wing case with $\alpha_0 = 15^\circ$ and twisted wing case with $\alpha_0 = 45^\circ$ and $\epsilon_0 = -45^\circ$, respectively.

It is observed that both cases exhibit very similar $C_{P_{rev},\kappa}$ and $C_{P_{rev},\sigma}$ contributions during the
Hovering flapping wings with dynamic twist

<figure><img src="image_18.png" /><figcaption>FIG. 17. Surface-integral and PPM-computed $C_{P,rev}$ as a function of non-dimensional time $t/T$ for the rigid wing case with $\alpha_0 = 15^\circ$ (a) and twisted wing case with $\alpha_0 = 45^\circ$, $\epsilon_0 = -45^\circ$ (b). Added-mass $C_{P,rev,\kappa}$, vortex-induced $C_{P,rev,\omega}$ and viscous $C_{P,rev,\sigma}$ contributions to revolving power are also shown.</figcaption></figure>

whole downstroke phase. Moreover, vortex-induced power are also very similar for both cases during the majority of the stroke. Specifically, the only difference is observed on the magnitude of the peak around $t/T \approx 2.38$ which accounts for the difference in $C_{P,rev}$ observed at this instant. Therefore, it is observed that, during the time interval $t/T \in [2.05; 2.20]$, differences in flow topology observed in the previous subsection do not significantly impact power consumption.

The left column on Figure 18 shows a time-sequence of $Q$-criterion isosurfaces $Qc^2/\bar{U}_{ref}^2 = 10$ colored by the contribution of vortex structures to the vortex-induced revolving power ($-2\rho Q\psi_{ey}$) for the rigid wing case with $\alpha_0 = 15^\circ$. The right column on Figure 18 displays the same quantity but for the twisted wing case with $\alpha_0 = 45^\circ$ and $\epsilon_0 = -45^\circ$. It can be seen that the RV in the rigid wing case has a negligible contribution to power $C_{P,rev,\omega}$, as indicated by its white contour. This is due to the amplitude of $\psi_{ey}$, which is low in the inboard region (see Appendix C). For the same reason, it can also be seen that the extended LEV in the twisted wing case does not significantly contribute to power consumption. Hence, during the time interval $t/T \in [2.05; 2.20]$, differences in flow topology do not significantly affect power consumption, while they do affect lift production.

Therefore, the inboard region of the twisted wing case with $\alpha_0 = 45^\circ$, $\epsilon_0 = -45^\circ$ is more efficient than that of the rigid wing case and this ultimately explains the enhanced efficiency of the twisted wing case observed on Figure 4.
<figure><img src="image_19.png" /><figcaption>FIG. 18. Q-criterion isosurfaces $Qc^2/\overline{U}_{ref}^2 = 10$ colored by the contribution to the vortex-induced revolving power ($-2\rho Q\psi_{ey}$) for the rigid wing case with $\alpha_0 = 15^\circ$ (left column) and twisted wing case with $\alpha_0 = 45^\circ$, $\epsilon_0 = -45^\circ$ (right column).</figcaption></figure>

V. CONCLUSION

Flapping wings have fostered the interest of both biologists and aerodynamicists these past three decades. Works were mainly focused on rigid wings for the purpose of understanding lifting mechanisms allowing insect and bird flight. Researches have revealed the key role of vortex dynamics generated by the flapping motion of the wing on lift production. Further drawing inspir-
Hovering flapping wings with dynamic twist

ration from natural flyers, flexibility was then introduced as an additional parameter in flapping
wing kinematics. Wing flexibility induces deformation of the wing which can be decomposed into
three mechanisms: spanwise bending, twisting and chordwise cambering. While passive defor-
mation of the wing are now relatively well documented, active wing deformation remains poorly
addressed. Yet, active wing deformation allows decoupling deformation mechanisms, which can
help understand their independant role on aerodynamic performance. Furthermore, it can be used
to enhance aerodynamic performance and actively control instantaneous forces for enhanced ma-
neuverability. The present study focused on the effect of active dynamic in-phase twisting on the
aerodynamic performance of a three-dimensional wing undergoing a flapping motion (revolving
and pitching) under hovering flight conditions. The twist amplitude was systematically varied for
different pitch angles.

Overall, it was shown that active twist yields comparable changes to pitching the whole wing.
That is, a change in wing twist generally led to time-averaged lift and lift-to-power ratio obtained
without twist but at a different pitch angle, clarifying previous works where it was observed that
twist generally leads to an increase in efficiency at the cost of lower force production²⁴,²⁵,³¹. Ac-
cordingly, most data obtained within this study were found to collapse onto a single curve when
plotted as a function of the mid-stroke pitch angle at 2/3 wing radius. The scaling of the forces
with this characteristic pitch angle can be directly related to LEV strength and LEV induced nor-
mal force.

Although data were found to scale, to leading order, with the mid-stroke pitch angle at 2/3 wing
radius, additional physical mechanisms that tend to enhance efficiency were identified. Specifi-
cally, positive wing twist (i.e. the wingtip angle of attack is larger than that of the wingroot) does
not lead to enhanced lift coefficient or efficiency when compared to most lift producing and most
efficient rigid wing motions. However, some cases with negative wing twist (i.e. the wingtip
angle of attack is lower than that of the wingroot) lead to enhanced efficiency with respect to
the most efficient rigid motions. The physical mechanisms behind this enhanced efficiency have
been identified and discussed using the force-power partitioning method (FPPM, an adaptation
of FMPM). Overall, it was shown that enhancement in efficiency arises from (i) the absence of
vortical structures near the wing root lower surface and (ii) the presence of an extended leading
edge vortex, when compared to the most efficient rigid wing case. In particular, these differences
in flow topology lead to enhanced lift during the early phase of the strokes without changes in
power consumption, which ultimately results in enhanced efficiency in the active twisting wing
Hovering flapping wings with dynamic twist

configuration.

As a conclusion, it was found that dynamic wing twist can be used as an alternative to wing
pitch to adapt both instantaneous and time-averaged lift and lift-to-power ratio. Furthermore, while
morphing wings do not outperform rigid wings in terms of mean lift production, it was shown that
they do provide enhanced efficiency. Future studies should explore to what extent those results are
relevant for forward flight conditions, for lower and higher Reynolds numbers and for wings with
larger aspect ratio where the LEV is known to be (partly) unstable.

ACKNOWLEDGMENTS

This work was supported by the French government program “Investissements d’Avenir” (EUR INTREE, Reference No. ANR-18-EURE-0010). The authors are grateful to GENCI-IDRIS, GENCI-CINES (Grant A0122A07178) and CALMIP (Grant 2022- p1425) for providing HPC resources.

DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author
upon reasonable request.

**Appendix A: Overset mesh technique**

Similarly to Diaz-Arriba et al.³³, the inner cylindrical domain, in which the wing is embedded,
is designed such that most of the vorticity remains inside of it. In our case, this is made possible
by the fact that the LEV remains attached to the wing along most of the wing span and during
the whole flapping motion. This way, interpolations between the inner, moving, domain and outer,
fixed, domain have only a limited influence on flow dynamics. Nonetheless, some vortex structures
do pass through the interpolation region that connects the solutions of the inner and outer meshed
regions. Figure 19 shows vorticity contours obtained in the r = 0.72R section at different instants
for the rigid α₀ = 45° case. Black line depicts the position of the interpolation region. It can be
observed that vortex structures that pass through this line do not suffer from significant numerical
dissipation due to the interpolation scheme. Note that the latter relies on a distance weighted
Hovering flapping wings with dynamic twist

<figure><img src="image_20.png" /><figcaption>FIG. 19. Spanwise vorticity contours in the $r = 0.72R$ spanwise cross-section at different instants $t/T$ for the rigid $\alpha_0 = 45^\circ$ case. The black line depicts the position of the interpolation region between the inner moving and background fixed meshes.</figcaption></figure>

approach with a typical ratio of two between inner and outer domain cell sizes.

## Appendix B: Meshgrid and timestep convergence

Figure 20 presents the period-averaged lift coefficient $\bar{C}_L$ as a function of the number $Nb$ of the flapping period for different spatial (a) and temporal (b) resolutions. Results are obtained for the twisted wing configuration with $\alpha_0 = 45^\circ$ and $\epsilon_0 = -45^\circ$. Horizontal lines represent $\pm 3\%$ of the $\bar{C}_L$ value obtained over the $10^{th}$ flapping period with resolutions $\Delta x_s = 0.005c$ and $\Delta t = T/1000$ (largest spatial resolution for convection tests on spatial resolution) in (a) and with resolutions $\Delta x_s = 0.02c$ and $\Delta t = T/1000$ (largest temporal resolution for convection tests on temporal resolution) in (b). Similarly to rigid wing configurations, it is shown that $\bar{C}_L$ is not severely affected by an increase in spatial and temporal resolutions with respect to the case with $\Delta x_s = 0.02c$ and $\Delta t = T/250$. Furthermore, it is observed that initial transients have sufficiently decayed after two flapping cycles.

## Appendix C: Auxiliary potentials fields

Figures 21 and 22 show the three-dimensional auxiliary potentials fields $\phi_{ey}$ (used for lift computation) and $\psi_{ey}$ (used for power associated with revolving motion), respectively. The distribution of $\phi_{ey}$ highlights the fact that vortices in the vicinity of the wing will strongly affect lift (see $\phi_{ey} = -1.5e-3$ isosurface). Note that two vortices with the same topology and strength but located one at wing root and the other at wing tip will have the same influence on lift production. On the contrary, distribution of $\psi_{ey}$ highlights the fact that vortical structures located outboard will
Hovering flapping wings with dynamic twist

<figure><img src="image_21.png" /></figure>

FIG. 20. Period-averaged lift coefficient $\bar{C}_L$ as a function of the number $Nb$ of the flapping period for different spatial and temporal resolutions. Data are obtained for the twisted wing configuration with $\alpha_0 = 45^\circ$ and $\epsilon_0 = -45^\circ$. (a) shows results obtained for different spatial resolutions $\Delta x_s = 0.02c, 0.01c$ and $0.005c$ for $\Delta t = T/1000$. (b) presents $\bar{C}_L$ for different temporal resolutions $\Delta t = T/1000, T/500$ and $T/250$ for $\Delta x_s = 0.02c$. Horizontal lines represent $\pm 3\%$ of the $\bar{C}_L$ value obtained over the $10^{th}$ period for the largest spatial (a) and temporal (b) resolutions.

<figure><img src="image_22.png" /></figure>

FIG. 21. 3D $\phi_{ey}$ field obtained for the rigid wing case with $\alpha_0 = 45^\circ$. 6 isosurfaces are used to illustrate the distribution.
Hovering flapping wings with dynamic twist

<figure><img src="image_23.png" /><figcaption>FIG. 22. 3D $\psi_{ey}$ field obtained for the rigid wing case with $\alpha_0 = 45^\circ$. 6 isosurfaces are used to illustrate the distribution.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>FIG. 23. $\bar{C}_L$ as a function of $\bar{C}_P$ for different pitch ($\alpha_0$) and twist ($\varepsilon_0$) angles. Data obtained for positive and negative $\varepsilon_0$ are shown in (a) and (b), respectively. White squares show rigid wing performance. Lightgrey, dark grey and black squares show performance obtained with $\varepsilon_0 = \pm 15^\circ, 30^\circ$ and $45^\circ$, respectively.</figcaption></figure>

more strongly affect revolving power than those located inboard.
**Appendix D: $\bar{C}_L$ as a function of $\bar{C}_P$**

Figure 23 shows the mean lift coefficient $\bar{C}_L$ as a function of the mean power coefficient $\bar{C}_P$ obtained for different pitch ($\alpha_0$) and twist ($\varepsilon_0$) angles. As observed on Figure 4 (a), positive twist does not benefit aerodynamic performance. In other words, positive twist does not help in reaching the large $\bar{C}_L$ and/or low $\bar{C}_P$ region indicated by the shaded region. Conversely, negative twist (Figure 23 (b)) does help in reaching the shaded region. Specifically, cases with $\alpha_0 = 45°$ and $\varepsilon_0 = -45°$, $\alpha_0 = 45°$ and $\varepsilon_0 = -30°$, $\alpha_0 = 30°$ and $\varepsilon_0 = -15°$ and $\alpha_0 = 30°$ and $\varepsilon_0 = -30°$ are found to lie in the shaded region. Yet, it is worth mentioning that none of the twisted wing cases exhibit $\bar{C}_P$ lower than that of the rigid wing case with $\alpha_0 = 0°$.

**REFERENCES**

¹M. H. Dickinson, F.-O. Lehmann, and S. P. Sane, “Wing rotation and the aerodynamic basis of insect flight,” *Science* **284**, 1954–1960 (1999).

²S. P. Sane and M. H. Dickinson, “The control of flight force by a flapping wing: lift and drag production,” *Journal of Experimental Biology* **204**, 2607–2626 (2001).

³W. Shyy, H. Aono, S. K. Chimakurthi, P. Trizila, C.-K. Kang, C. E. S. Cesnik, and H. Liu, “Recent progress in flapping wing aerodynamics and aeroelasticity,” *Progress in Aerospace Sciences* (2010), 10.1016/j.paerosci.2010.01.001.

⁴H. Liu, S. Wang, and T. Liu, “Vortices and forces in biological flight: Insects, birds, and bats,” *Annual Review of Fluid Mechanics* **56**, 147–170 (2024).

⁵C. P. Ellington, C. van den Berg, A. P. Willmott, and A. L. R. Thomas, “Leading-edge vortices in insect flight,” *Nature* **384**, 626–630 (1996).

⁶D. Lentink and M. H. Dickinson, “Rotational accelerations stabilize leading edge vortices on revolving fly wings,” *Journal of Experimental Biology* **212**, 2705–2719 (2009).

⁷D. Garmann and M. Visbal, “Dynamics of revolving wings for various aspect ratios,” *Journal of Fluid Mechanics* **748**, 932–956 (2014).

⁸T. Jardin, “Coriolis effect and the attachment of the leading edge vortex,” *Journal of Fluid Mechanics* **820**, 312–340 (2017).

⁹J. D. Eldredge and A. R. Jones, “Leading-edge vortices: Mechanics and modeling,” *Annual Review of Fluid Mechanics* **51**, 75–104 (2019).
<footnote n="10">¹⁰L. Chen, B. Cheng, and J. Wu, “Vorticity dynamics and stability of the leading-edge vortex on revolving wings,” Physics of Fluids **35**, 091301 (2023).</footnote>

<footnote n="11">¹¹K. N. Lucas, N. Johnson, W. T. Beaulieu, E. Cathcart, G. Tirrell, S. P. Colin, B. J. Gemmell, J. O. Dabiri, and J. H. Costello, “Bending rules for animal propulsion,” Nature Communications **5** (2014), 10138–101429.</footnote>

<footnote n="12">¹²L. Zheng, T. L. Hedrick, and R. Mittal, “Time-varying wing-twist improves aerodynamic efficiency of forward flight in butterflies,” PLOS ONE **8**, 1–10 (2013).</footnote>

<footnote n="13">¹³L. Zhao, Q. Huang, X. Deng, and S. Sane, “The effect of chord-wise flexibility on the aerodynamic force generation of flapping wings: Experimental studies,” in *2009 IEEE International Conference on Robotics and Automation* (2009) pp. 4207–4212.</footnote>

<footnote n="14">¹⁴H. Dai, H. Luo, and J. F. Doyle, “Dynamic pitching of an elastic rectangular wing in hovering motion,” Journal of Fluid Mechanics **693**, 473–499 (2012).</footnote>

<footnote n="15">¹⁵C. L. Shah, D. Majumdar, C. Bose, and S. Sarkar, “Controlling the chaotic wake of a flapping foil by tuning its chordwise flexibility,” Journal of Fluids and Structures **127**, 1041–1054 (2024).</footnote>

<footnote n="16">¹⁶R. J. Wootton, “Support and deformability in insect wings,” Journal of Zoology **193**, 447–468 (1981).</footnote>

<footnote n="17">¹⁷S. Heathcote, D. Martin, and I. Gursul, “Flexible flapping airfoil propulsion at zero freestream velocity,” AIAA Journal **42**, 2196–2204 (2004).</footnote>

<footnote n="18">¹⁸M. Vanella, T. Fitzgerald, S. Preidikman, E. Balaras, and B. Balachandran, “Influence of flexibility on the aerodynamic performance of a hovering wing,” Journal of Experimental Biology **212**, 95–105 (2009).</footnote>

<footnote n="19">¹⁹J. D. Eldredge, J. Toomey, and A. Medina, “On the roles of chord-wise flexibility in a flapping wing with hovering kinematics,” Journal of Fluid Mechanics **659**, 94–115 (2010).</footnote>

<footnote n="20">²⁰S. Heathcote, Z. Wang, and I. Gursul, “Effect of spanwise flexibility on flapping wing propulsion,” Journal of Fluids and Structures **24**, 183–199 (2008).</footnote>

<footnote n="21">²¹T. Nakata and H. Liu, “Aerodynamic performance of a hovering hawkmoth with flexible wings: a computational approach,” Proceedings of the Royal Society B: Biological Sciences **279**, 722–731 (2012).</footnote>

<footnote n="22">²²R. Addo-Akoto, J.-S. Han, and J.-H. Han, “Roles of wing flexibility and kinematics in flapping wing aerodynamics,” Journal of Fluids and Structures **104**, 10331–10354 (2021).</footnote>

<footnote n="23">²³D. Diaz-Arriba, T. Jardin, N. Gourdain, F. Pons, and L. David, “Experiments and numerical simulations on hovering three-dimensional flexible flapping wings,” Bioinspiration and Biomimetics</footnote>
17, 065006 (2022).

²⁴D. Gang and M. Sun, “Effects of unsteady deformation of flapping wing on its aerodynamic forces,” *Applied Mathematics and Mechanics (English Edition)* **29**, 731–743 (2008).

²⁵D. Gang and M. Sun, “Effects of wing deformation on aerodynamic forces in hovering hover-flies,” *The Journal of experimental biology* **213**, 2273–83 (2010).

²⁶W. Thielicke and E. J. Stamhuis, “The effects of wing twist in slow-speed flapping flight of birds: trading brute force against efficiency,” *Bioinspiration and Biomimetics* **13**, 056015 (2018).

²⁷K. Joshi and S. Bhattacharya, “The unsteady force response of an accelerating flat plate with controlled spanwise bending,” *Journal of Fluid Mechanics* **933**, A56 (2022).

²⁸C. Soto and S. Bhattacharya, “The effect of dynamic twisting on the flow field and the unsteady forces of a heaving flat plate,” *Bioinspiration and Biomimetics* **18**, 026010 (2023).

²⁹Y. Dong, B. Song, W. Yang, and D. Xue, “A numerical study on the aerodynamic effects of dynamic twisting on forward flight flapping wings,” *Bioinspiration and Biomimetics* **19**, 026013 (2024).

³⁰F. Bouard, T. Jardin, and L. David, “Aerodynamics of flapping wings with passive and active deformation,” *Journal of Fluids and Structures* **128**, 104139 (2024).

³¹Y.-J. Han, H.-H. Yang, and J.-H. Han, “Twist-Coupled Flapping Mechanism for Bird-Type Flapping-Wing Air Vehicles,” *Journal of Mechanisms and Robotics* **15**, 051017 (2023).

³²T. Jardin and T. Colonius, “On the lift-optimal aspect ratio of a revolving wing at low reynolds number,” *Journal of The Royal Society Interface* **15**, 20170933 (2018).

³³D. Diaz-Arriba, T. Jardin, N. Gourdain, F. Pons, and L. David, “Numerical investigation of three-dimensional asymmetric hovering flapping flight,” *Physics of Fluids* **33**, 111907 (2021).

³⁴Q. Wang, J. F. L. Goosen, and F. van Keulen, “Optimal pitching axis location of flapping wings for efficient hovering flight,” *Bioinspiration and Biomimetics* **12**, 056001 (2017).

³⁵L.-G. Liu, G. Du, and M. Sun, “Aerodynamic-force production mechanisms in hovering mosquitoes,” *Journal of Fluid Mechanics* **898**, A19 (2020).

³⁶Y. Bury and T. Jardin, “Transitions to chaos in the wake of an axisymmetric bluff body,” *Physics Letters A* **376**, 3219–3222 (2012).

³⁷T. Jardin and L. David, “Root cutout effects on the aerodynamics of a low-aspect-ratio revolving wing,” *AIAA Journal* **55**, 2717–2726 (2017).

³⁸T. Jardin and N. Doué, “Influence of pitch rate on freely translating perching airfoils,” *Journal of Fluid Mechanics* **873**, 49–71 (2019).
<footnote n="39"><sup>39</sup>L. Quartapelle and M. Napolitano, “Force and moment in incompressible flows,” AIAA Journal **21**, 911–913 (1983).</footnote>

<footnote n="40"><sup>40</sup>K. Menon and R. Mittal, “On the initiation and sustenance of flow-induced vibration of cylinders: insights from force partitioning,” Journal of Fluid Mechanics **907**, A37 (2021).</footnote>

<footnote n="41"><sup>41</sup>K. Menon and R. Mittal, “Quantitative analysis of the kinematics and induced aerodynamic loading of individual vortices in vortex-dominated flows: A computation and data-driven approach,” Journal of Computational Physics **443**, 110515 (2021).</footnote>

<footnote n="42"><sup>42</sup>K. Menon and R. Mittal, “Significance of the strain-dominated region around a vortex on induced aerodynamic loads,” Journal of Fluid Mechanics **918**, R3 (2021).</footnote>

<footnote n="43"><sup>43</sup>K. Menon, S. Kumar, and R. Mittal, “Contribution of spanwise and cross-span vortices to the lift generation of low-aspect-ratio wings: Insights from force partitioning,” Phys. Rev. Fluids **7**, 114102 (2022).</footnote>

<footnote n="44"><sup>44</sup>C. Zhang, T. L. Hedrick, and R. Mittal, “Centripetal acceleration reaction: An effective and robust mechanism for flapping flight in insects,” PLOS ONE **10**, 1–16 (2015).</footnote>

<footnote n="45"><sup>45</sup>Y. Zhu and K. Breuer, “Flow-induced oscillations of pitching swept wings: stability boundary, vortex dynamics and force partitioning,” Journal of Fluid Mechanics **977**, A1 (2023).</footnote>

<footnote n="46"><sup>46</sup>J. Hunt, A. Wray, and P. Moin, “Eddies, streams, and convergence zones in turbulent flows,” Center for turbulence research report CTR-S88, 193–208 (1988).</footnote>

<footnote n="47"><sup>47</sup>A. Medina and A. R. Jones, “Leading-edge vortex burst on a low-aspect-ratio rotating flat plate,” Physical Review Fluids **1**(4), 044501 (2016).</footnote>