# Thrust force is tuned by the rigidity distribution in insect-inspired flapping wings

Roméo Antier, Benjamin Thiria, Ramiro Godoy-Diana

► To cite this version:

Roméo Antier, Benjamin Thiria, Ramiro Godoy-Diana. Thrust force is tuned by the rigidity distribution in insect-inspired flapping wings. Journal of Fluids and Structures, 2023, 124, pp.104043. 10.1016/j.jfluidstructs.2023.104043. hal-04290440

HAL Id: hal-04290440

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

Submitted on 16 Nov 2023

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 CC BY 4.0 - Attribution - International License
# Thrust force is tuned by the rigidity distribution in insect-inspired flapping wings

Roméo Antier, Benjamin Thiria, and Ramiro Godoy-Diana

Laboratoire de Physique et Mécanique des Milieux Hétérogènes (PMMH),

CNRS UMR 7636, ESPCI Paris-Université PSL,

Sorbonne Université, Université Paris Cité, F-75005 Paris, France

We study the aerodynamics of a flapping flexible wing with a two-vein pattern that mimics the elastic response of insect wings in a simplified manner. The experiments reveal a non-monotonic variation of the thrust force produced by the wings when the angle between the two veins is varied. An optimal configuration is consistently reached when the two veins are spaced at an angle of about 20 degrees. This value is in the range of what has been measured in the literature for several insect species. The deformation of the wings is monitored during the experiment using video recordings, which allows to pinpoint the physical mechanism behind the non-monotonic behaviour of the force curve and the optimal distribution of the vein network in terms of propulsive force.

## I. INTRODUCTION

Quoting Wootton [1]: "In considering insect wings, whether for comparative illustration or aerodynamic analysis, some simplifications are inevitable. Two in particular are common: to regard the wing as essentially flat, and as effectively rigid. Neither is true, and the latter can be seriously misleading." And the same can be said for most flapping wings and fins, where the structural deformation that accompanies the back and forth motion is a fundamental element of the dynamical balance [2, 3]. In particular, the periodic stroke reversals of flapping wings and their associated cycle of acceleration and deceleration give rise to a rich variety of vortex structures that are crucial players in the unsteady aerodynamic mechanisms inherent to flapping flight—see e.g. [4–7], for a review. Another noteworthy point is that these mechanisms are tuned with the deformation dynamics of the wings [8, 9], where specific features such as the passive wing pitch reversal observed in some insects [10] or the active camber control used by bats [11] are determinant in the cycle of aerodynamic force production. Several works [8, 12, 13] have addressed the problem of wing deformation (see e.g. [14] for a review) and models have usually decomposed the main deformation modes as a combination of spanwise and chordwise bending [15–17]. In the case of insects, a network of veins confers their wings an anisotropic rigidity [3, 18, 19], which governs the passive responses of the wings to aerodynamic, inertial and occasional impact forces [20, 21]. The structural function of veins is not straightforward, and it coexists with their other roles as transmission conduits of air and hemolymph, and as sensory elements [19]. However, a few main features are recurrent, such as the veins being the thickest closest to the wing root, tapering towards the tip and trailing edge of the wings [22, 23]. Another observation is that most insect wings present a zone near the leading edge stiffened by thick veins and relief (see e.g. [20]), with thinner veins elsewhere that will let the wing membrane deform more easily during the flapping motion. A secondary stiffened axis oriented obliquely at an angle with the leading edge is also present in many insect

wings. Interestingly, the angle between the leading edge and this oblique stiffened area is narrowly-distributed around 15° to 30° (see e.g. for dipteran wings [18, 24]); the stiffness of this zone is provided by the combined effect of veins and corrugation [24]. The details of the local venation are essential here, and it has been suggested that extant vein patterns have resulted from evolutionarily convergent vein fusions [25].

The goal of this paper is to examine the effect of different patterns of wing rigidity on the flapping wing aerodynamics by using a minimal model that mimics the elastic response of insect wings in a simplified manner. We use a two-vein pattern (see Fig. 1) with one main vein along the leading edge of the wing, and a secondary vein also attached at the root of the wing but that extends obliquely at a specific angle with respect to the leading edge. The angle between the two veins is the main experimental parameter explored. The experiments reveal a non-monotonic variation of the average thrust force produced by the flapping wings with a local optimum when the two veins are spaced at an angle of about 20°, which is in the range of the typical angles observed in insects. An explanation of the physical mechanisms involved is proposed using observations of the instantaneous kinematics of the flapping wing deformation.

## II. EXPERIMENTAL SETUP AND METHODS

### A. Wings and flapping system

A flapping system with two wings is used for the experiments. The wings are composed of a 3D-printed skeleton with two veins and a thin membrane. The model wing is a Zimmerman planform (see Fig.1 C) [26, 27], that represents a simplified hummingbird or insect wing shape. The wings are shaped like two quarter ellipses. A quarter ellipse with half major axis the span and half minor axis one quarter of the chord. The leading edge of the wing follows the curve of this first quarter ellipse. The second quarter ellipse is connected to the first one
<figure><img src="image_2.png" /><figcaption>FIG. 1: (A) and (B) Two examples of the rigidity distribution in insect wings (Figures from [1]). Supporting areas (stippled), deformable areas (unstippled) and flexion lines (dashed) in (A) *Syrphus ribesii* (Diptera); (B) *Vespula germanica* (Hymenoptera). m.f.l., median flexion line; cl.f., claval furrow; tr.f.l., transverse flexion line. Scale lines = 5 mm. (C) Model wing used in the present work. (D) and (E) Frontal and side views, respectively, of the system mounted on the force sensor. In (D) several snapshots are superposed to illustrate the flapping wing motion.</figcaption></figure>

along the span. The half major axis is also the length of the wingspan, and the half minor axis is three quarters of the chord. The length of the wing is $R = 75$ mm, and the root chord measures $c = 25$ mm, the mean chord thus being $c_m = 19.6$ mm. The aspect ratio of the wing is $AR = R/c_m = 3.82$. The wing shape is cut on a polyethylene terephthalate (PET) film of thickness 30 µm. The cutting is done with a laser cutting machine (CO2 with infrared ray, Epilog Laser, type Helix). The Young's modulus of the membrane is 4 GPa.

The membrane is supported by two reinforcements disposed respectively on the leading edge and along a direction making an angle $\beta$ with the leading edge (see Fig.1 C). The reinforcements play the role of veins. The material of the veins is polylactic acid (PLA) with a Young's modulus of $4 \pm 0.3$ GPa. In order to have a symmetric deformation, a vein skeleton is glued on each side of the membrane using a Teroson SB2444 rubber adhesive. This adhesive is very elastic and allows the membrane to slide between the veins without detaching. The veins are 1 mm wide. The thicknesses on each side of the membrane are of 480 µm for the leading edge and 240 µm for the radial vein. The angle between the two veins $\beta$, varies between 10° and 90° in 5° steps. We also made a wing consisting of a membrane and a single vein at the leading edge. This wing is referenced by the case $\beta = 0°$. The leading edge vein is extended by 3.5 mm toward the center of the ellipse so that it can be connected

to the flapping mechanism. The length of the radial vein, $R_v$, evolves as a function of the angle $\beta$ through the relationship:

$$R_v = \frac{1}{\cos(\beta)} \left( \frac{1}{R^2} + \frac{\tan^2 \beta}{c^2} \right)^{-1/2}$$

Because $R_v$ diminishes with increasing $\beta$, the total mass of the wing also diminishes slightly. The difference in mass between the heaviest wing ($\beta = 5°$) and the lightest wing ($\beta = 90°$) is 16%.

Two wings are mounted on a flapping system based on the DelFly design [28] obtained by dismantling a commercially-available flapping-wing bird toy to keep only the motor and crank mechanism. The system, powered externally, allows to generate a sinusoidal planar flapping motion with an amplitude of 32° for frequencies $f$ ranging from 5 to 20 Hz. Front view and side view photos of the system are shown in Fig. 1 (D) and (E), respectively.

## B. Wings scaling analysis

In order to asses how far, dynamically, our rudimentary model wings are from the case of an insect wing, it is convenient to examine a few dimensionless quantities. The main elements of the flapping wing problem are the aerodynamic force, the elastic bending rigidity of the wing, and its inertia [32–34]. We can use as a basic model a flexible plate (dimensions: mean chord $c_m$, span $R$ and thickness $h$, density $\rho_s$ —i.e. surface mass density $\mu_s = \rho_s h—$, and elastic modulus $E$) that will bend under the action of its own inertia and of the aerodynamic forces. For a flapping motion in hovering characterised by an angular amplitude $\psi_0$ and frequency $\omega = 2\pi f$, before considering the wing deformation we can already define [2]: (i) the Reynolds number, written in terms of a reference flapping velocity $U_{ref} = 2\psi_0 R f$, the density $\rho$ and dynamic viscosity $\eta$ of air, and using the mean chord $c_m$ as reference length scale $L_{ref}$:

$$Re = \frac{\rho U_{ref} L_{ref}}{\eta} = \frac{\rho 2 \psi_0 R f c_m}{\eta}, \quad (1)$$

which governs the aerodynamic regime by setting the importance of fluid inertial versus viscous effects; and (ii), the reduced frequency

$$k = \frac{\omega c_m}{U_{ref}} = \frac{\pi}{\psi_0 AR}, \quad (2)$$

which in the present hovering case does not depend explicitly on the physical frequency because the reference velocity is the flapping velocity that is itself proportional to the frequency.

Now, to estimate the effects of aerodynamic loading and wing inertia measured against the elastic response
<table><caption>TABLE I: Wing morphological and material properties, kinematic parameters, and dimensionless numbers for a few insect species and for the model wings. Data from [2, 15, 29–31].</caption><thead><tr><th>Parameter</th><th>units</th><th>Hawkmoth<br>(Manduca sexta)</th><th>Hoverfly<br>(Eristalis tenax)</th><th>Bumblebee<br>(Bombus terrestris)</th><th>European honey bee<br>(Apis mellifera)</th><th>Model wings<br>of the present study</th></tr></thead><tbody><tr><td>R</td><td>(m)</td><td>0.049</td><td>0.009</td><td>0.016</td><td>0.0097</td><td>0.075</td></tr><tr><td>S</td><td>(m<sup>2</sup>)</td><td>8.91×10<sup>-4</sup></td><td>3.7×10<sup>-5</sup></td><td>1.1×10<sup>-4</sup></td><td>4.2×10<sup>-5</sup></td><td>0.0015</td></tr><tr><td>c<sub>m</sub></td><td>(m)</td><td>0.018</td><td>0.004</td><td>0.007</td><td>0.004</td><td>0.02</td></tr><tr><td>m<sub>w</sub></td><td>(kg)</td><td>4.7×10<sup>-5</sup></td><td>6.0×10<sup>-7</sup></td><td>1.25×10<sup>-6</sup></td><td>5.0×10<sup>-7</sup></td><td>1.5×10<sup>-4</sup></td></tr><tr><td>μ<sub>s</sub></td><td>(kg m<sup>-2</sup>)</td><td>0.0527</td><td>0.0162</td><td>0.0118</td><td>0.0119</td><td>0.1018</td></tr><tr><td><sup>†</sup> EI<sub>beam</sub></td><td>(N m<sup>2</sup>)</td><td>8.0×10<sup>-6</sup></td><td>7.7×10<sup>-6</sup></td><td>7.7×10<sup>-6</sup></td><td>1.82×10<sup>-6</sup></td><td>[0.77 - 47.5]×10<sup>-6</sup></td></tr><tr><td><sup>†</sup> EI&#773;</td><td>(N m)</td><td>1.63×10<sup>-4</sup></td><td>8.56×10<sup>-4</sup></td><td>4.81×10<sup>-4</sup></td><td>1.88×10<sup>-4</sup></td><td>[1 - 63]×10<sup>-5</sup></td></tr><tr><td>ψ<sub>0</sub></td><td>(°; rad)</td><td>57; 0.99</td><td>51; 0.90</td><td>60; 1.05</td><td>65; 1.13</td><td>16; 0.28</td></tr><tr><td>f</td><td>(Hz)</td><td>25</td><td>210</td><td>150</td><td>250</td><td>[12 - 20]</td></tr><tr><td>Re</td><td>-</td><td>2933</td><td>934</td><td>2204</td><td>1577</td><td>[653 - 1089]</td></tr><tr><td>C<sub>Y</sub></td><td>-</td><td>0.131</td><td>0.001</td><td>0.009</td><td>0.008</td><td>[0.002 - 0.311]</td></tr><tr><td>N<sub>ei</sub></td><td>-</td><td>2.34</td><td>0.02</td><td>0.11</td><td>0.14</td><td>[0.15 - 24.85]</td></tr><tr><td>k</td><td>-</td><td>1.17</td><td>1.58</td><td>1.24</td><td>1.23</td><td>2.95</td></tr><tr><td>A</td><td>-</td><td>2.69</td><td>3.10</td><td>2.18</td><td>2.41</td><td>3.82</td></tr></tbody></table>

of the wing, we can use, respectively, a Cauchy number [35, 36]:

$$C_Y = \frac{\frac{1}{2} \rho U_{\text{ref}}^2 L_{\text{ref}}^3}{EI} = \frac{2 \rho R^2 \psi_0^2 f^2 c_m^3}{EI}, \quad (3)$$

which characterizes the deformation of the wing under the effect of the fluid flow, and the elasto-inertial number [33]:

$$\mathcal{N}_{ei} = \frac{\mu_s a_{\text{ref}} L_{\text{ref}}^3}{EI} = \frac{4\pi^2 \mu_s R \psi_0 f^2 c_m^3}{EI}, \quad (4)$$

which characterizes the deformation of the wing under the effect of its own inertia. N<sub>ei</sub> is written in terms of a reference acceleration a<sub>ref</sub> = Rψ<sub>0</sub>ω<sup>2</sup>. EI in Eqs. 3 and 4 is an average plate bending rigidity.

These dimensionless numbers can be used to give an indicative picture of the model wings in comparison with insect (or other) wings in a global parameter space, as shown in Table I, and comfort the idea of using the present experiment to examine the effect of the rigidity distribution of the wing on its aerodynamic performance.

### C. Force sensor

The system is mounted on a Schunk FT-Nano 17 6-axis force sensor as shown in Fig. 1 D, such that the average

propulsive force produced by the flapping wings points towards the x-direction. In a right-handed cartesian reference frame, the weight of the device is thus directed towards the negative z-direction. In what follows, we focus on the forward component of the force F<sub>x</sub>, which in the present setup is the most relevant concerning the aerodynamic force production because the forward component is perpendicular to the stroke plane (as in the merry-go-round setup of [33, 37]). The reciprocal motion of the wings and the symmetry of the setup determine that the F<sub>y</sub> and F<sub>z</sub> components of the force as well as the M<sub>x</sub> and M<sub>z</sub> components of the moment average to zero over each flapping period. The M<sub>y</sub> component has non-zero mean, but what can be learned from its dynamics in the tethered frame of the present experiment is redundant from what is obtained from the analysis of the forward force F<sub>x</sub>.

A typical time series of the measured force signal is shown in Fig. 2 (A). The signal is noisy because the forces produced by the wings were close to the limit of the measurement range of our sensor. Nonetheless, the periodicity driven by the flapping motion is clearly visible, as highlighted in the figure by the running average also plotted. These time-resolved measurements were robustly repeatable. The Fourier transform of the signal—shown in Fig. 2 (B)—has its largest peak at twice the flapping frequency. This is expected because two peaks of force are produced over one period when the wing instantaneous velocity is highest during the upstroke and

<footnote>†the values EI of average plate bending rigidity were obtained as EI = EI<sub>beam</sub>/R, where EI<sub>beam</sub> is the chord-wise flexural stiffness. EI<sub>beam</sub> values for insects [15] come from an indirect measurement of an equivalent beam performed by applying a point force to bend the wing and using the measured force F and wing displacement δ to calculate an overall flexural stiffness EI<sub>beam</sub> = Fl<sup>3</sup>/3δ, l being the effective beam length.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 2: Thrust force measurements. (A) Time series $F_x(t)$ and (B) its frequency content obtained by FFT for a typical case ($\beta = 20^\circ$ and $f = 17.5$ Hz). In (A) a running mean of the signal is shown (solid blue line) as well as its average value (dashed line). (C) Average force $\bar{F}_x$ as a function of the radial vein angle $\beta$ for several frequencies ($f = 12.5, 14.5, 16.5, 17.5$ and $18.5$ Hz, from darkest to lightest color respectively). Each point is the mean value of several runs with identical parameters and the corresponding standard deviations are represented as error bars. (D) Thrust coefficient $C_f$ for the same experimental data computed using Eq. 5; in the inset, $C_f$ is normalized by the maximum value $C_{f\max}$ of each frequency series.</figcaption></figure>

the downstroke. The average force, marked as a dashed horizontal line in Fig. 2 (A), is the main output used as a performance probe as the $(\beta, f)$-parameter space is explored. The results are summarised in Fig. 2 (C), where this time-averaged force $\bar{F}_x$ is plotted as a function of the radial vein angle $\beta$ for several frequencies.

### D. Kinematics tracking

In addition to the force measurements, the motion of the wings was tracked using a Phantom Miro M120 high-speed camera recording 1920 × 1200 pixel² images at

800Hz. Fig. 3 shows time series of a side view (A) and a front view (B) to give a qualitative picture of the deformation of the wing during the flapping cycle. In order to quantify the wing deformation, four points of interest were tracked using ImageJ [38]: two at the leading edge (at the root and at the tip), and two at the trailing edge (at the point where the radial vein ends, and at end of the largest chord section).
<figure><img src="image_4.png" /><figcaption>FIG. 3: (A) Side and (B) front views of a flapping sequence of a wing with $\beta = 40^{\circ}$ at 17.5 Hz.</figcaption></figure>

### III. RESULTS AND DISCUSSION

#### A. Thrust coefficient

As mentioned above, the time average of the force $\bar{F}_x$ presented in Fig. 2 (C) is the main performance indicator of the system as a function of the parameter space constituted by the radial vein angle $\beta$ and the flapping frequency $f$. A first step in the analysis is naturally to find a dimensionless representation of the data that is plotted in physical units in Fig. 2. To do so, we define a thrust coefficient

$$C_f = \frac{\bar{F}_x}{\frac{1}{2} \rho u_{wing}^2 S} \qquad (5)$$

where $u_{wing} = 2\pi f A$ is the characteristic flapping speed defined by the frequency $f$ and amplitude $A$ of the flapping motion. We define a nominal amplitude $A = R \sin(32^{\circ})$ based on the wing length and the flapping angular amplitude. The reference surface $S$ is the area of the two wings. The thrust coefficient is presented in Fig. 2 (D). Two main observations can be made: on the one hand, the measurements for each frequency constitute a non-monotonic curve. As the angle $\beta$ of the radial vein is increased from zero, a clear maximum occurs

around $\beta \approx 20^{\circ}$, followed by a minimum at $\beta \approx 35^{\circ}$. Further increasing $\beta$ makes the propulsive force grow again until it reaches a similar value to that of the first maximum observed at $\beta \approx 20^{\circ}$. The second observation is that when the frequency is increased, the performance curve is shifted to higher values whilst keeping a fairly similar shape. This can be seen clearly by normalising the curve corresponding to each flapping frequency by its maximum value $C_{fmax}$, as shown in the inset of Fig. 2 (D). That the thrust force increases with the flapping frequency is of course an expected result, which has been reported for similar systems in the literature [39, 40].

In what follows we analyse these results in light of the wing deformation observations.

#### B. Wing deformation kinematics

An overview of the flapping motion captured from side and top views is presented in Fig. 3. The main point of this visualisation is to examine the typical behaviour of the wing and to identify the basic elements of its deformation dynamics. The phase lag between the leading edge and trailing edge has been used in the literature [37] to explain the performance increase of flexible wings
<figure><img src="image_5.png" /><figcaption>FIG. 4: (A) and (B) Schematic representations of the flexible wing moving at speed $u_{wing}$. In (A), a section in the $xz$ plane is pictured, the thick blue arrow represents the total aerodynamic force which points more in the forward direction than the wing is bent. The trailing angle $\theta$ and the projected surface $S_x$ are also shown. (B) shows a three-dimensional sketch. (C) and (D) Measured values of the projected surface $S_x$ (C) and the trailing angle (D) as a function of $\beta$ for $f = 14.5$ and $17.5$ Hz.</figcaption></figure>

with respect to rigid wings. Considering that the average main component of the aerodynamic force is perpendicular to the stroke plane, the advantage of the flexible wing comes from the redirection of the force in the useful direction. A representation of this idea for a section of a flexible wing is shown in Fig 4 (A), which defines the two parameters that will be used in the following: the projection of the wing surface on the stroke plane $yz$, defined as $S_x$ —see also Fig 4 (B)—, and the trailing angle $\theta$, both measured at the instant of maximum flapping velocity (i.e. when the wing passes the horizontal position). Note that the trailing angle $\theta$ is different from the usual aircraft trailing edge angle, defined as the angle between the tangents of the upper and lower airfoil at trailing edge, indicating trailing edge sharpness. Now, $S_x$ can be used as a measure of the aforementioned force redirection. If the wing is considered as a homogeneous plate bending under its first mode of deformation, this picture is sufficient to explain the basic physical mechanism driving the performance of a flexible wing [8, 33]. The radial vein complicates the picture because each section of the wing behaves now differently —see Fig. 5 (B) and (C). In particular, the phase lag of the trailing edge becomes different depending on the span-wise position. Nonetheless, we can still examine how $S_x$ changes with

$\beta$. This is shown in Fig. 4 (B) for the cases of $f = 14.5$ and $17.5$ Hz. As $\beta$ increases from zero, the radial vein starts preventing part of the wing to bend and $S_x$ diminishes. This trend saturates at $\beta \approx 40^{\circ}$ until $\approx 60^{\circ}$, after what $S_x$ increases again to larger values. Recalling the force measurements of Fig. 2, where the thrust minima are observed around $\beta \approx 40^{\circ}$, reinforces the idea of a larger projected surface $S_x$ being a necessary feature for increasing thrust production.

### C. Projected wing surface and trailing angle

To go further, it is useful to describe the wing in terms of its two sections: the first section is the area comprised between the leading edge and the radial vein and the second one that between the radial vein and the trailing edge. We define these, respectively, as the *inter-vein area* and the *trailing area*. Because the wing is built with portions of ellipses, the areas of these two sections can be expressed analytically, as a function of $\beta$, which gives the curves shown in Fig. 5 (A). Panels (B), (C), and (D) in Fig. 5 show three examples for different angles of the radial vein, when the wing passes the horizontal position, with the perimeters of the two areas highlighted. Since
<figure><img src="image_6.png" /><figcaption>FIG. 5: (A) Areas of the inter-vein and trailing sections of the wing as a function of $\beta$; (B), (C), and (D) frontal snapshots of the wing for different values of $\beta$ with the wing sections highlighted.</figcaption></figure>

the snapshots are frontal views of the wing, the addition of the highlighted areas is actually $S_x$.

For small angles between the veins, typically $\beta < 20^\circ$, the trailing area is larger than the inter-vein area, this implies that $S_x$ is dominated by the deformation of the trailing area —see Fig. 5 (B). The larger this free surface, the larger its deformation. On the contrary, for large $\beta$, typically $\beta > 50^\circ$, the trailing area is very small and hardly deforms at all. Its influence is then small in the generation of aerodynamic forces. The surface between the veins is the largest and its swelling is at the origin of the redistribution of the aerodynamic forces —see Fig. 5 (D) for the limit case of $\beta = 90^\circ$ where the trailing area has vanished and the whole wing is the inter-vein area.

Summarising, the changes in the projected area as a function of the radial vein angle give us a first physical insight to explain the non-monotonic behaviour of the propulsive force observed in Fig 2. For lower values of $\beta$, the force production is dominated by the trailing area, whereas for higher values, typically above $\beta \approx 40^\circ$, it is the inter-vein area that contributes the most.

Now, considering firstly the lower angles, say $\beta \lesssim 40^\circ$, the force measurements have a maximum value at $\beta \approx 20^\circ$. This means that the thrust force does not solely depend on the surface (or projected surface) of the trailing area, which is largest at $\beta = 0^\circ$, but also on the way the wing is bent relatively to the incident wind. Examining side views of the flapping wing (see Fig. 6) brings evidence of the reason for the suboptimal performance of the wings with very low values of $\beta$. Because the radial vein also imposes the maximum length of the trailing

area that can bend, a small angle $\beta$ means that the trailing surface is long enough to bend on itself as shown in Fig. 6, thus losing the aerodynamic benefit of flexibility. In practice, the excessive bending of the trailing area determines that its orientation is suboptimal during large portions of the flapping cycle.

This can be examined quantitatively by tracking the trailing angle—see Fig. 4 (A)—at the instant of maximum flapping velocity, as presented in Fig. 4 (D) as a function of $\beta$ for two different frequencies. We use the dimensionless representation $\theta/\phi$ introduced by [37], where the angle of the incoming wind $\phi$ is in the present case equal to $90^\circ$ since there is no incident velocity on the static wing because the system is fixed in the lab reference frame. The measurements of the trailing angle in figure Fig. 4 (D) bring a clear explanation underlying the maximum of aerodynamic force measured for the wing with the radial vein at $\beta \approx 20^\circ$: as in [37], this optimum coincides with the best alignment of the trailing angle $\theta$ and the angle $\phi$ of the local wind seen by the translating wing. As $\beta$ increases, the trailing angle goes to zero because it is measured at the tip, which is part of the trailing area that deforms less and less and behaves as in a rigid wing when $\beta$ tends to $90^\circ$. For these larger values of $\beta$ the influence of the trailing area concerning thrust production diminishes, so the measurement of the trailing edge as in Fig. 4 (D) becomes irrelevant. As $\beta$ becomes larger, the main part of the force production is ensured by the inter-vein area, which represents most of the total wing surface.

Coming back to Fig. 4 (C), the increase of the projected area $S_x$ for $\beta > 50^\circ$ is driven by the swelling of the inter-vein area, which can be seen in Fig. 5 (D) for the limit case of $\beta = 90^\circ$. We can hypothesise that the physical mechanism enhancing thrust at these large values of the radial vein angle $\beta$ should be similar to the case described for the trailing area. However, a quantitative picture would need the tracking of the whole trailing edge and not just of a single point as we have done to produce Fig. 4 (D), which is out of reach of the experiments reported here.

## IV. CONCLUDING REMARKS

We have examined the role of non-homogeneous wing stiffness in the aerodynamic force production by flapping wings, using a simple model with a two-vein skeleton. The main observation is that the thrust produced by the wings varies non-monotonically with changes in the angle $\beta$ between the two veins that constitute the skeleton: the leading edge vein and the radial vein. A local optimum of the aerodynamic performance is observed for $\beta \approx 20^\circ$, which is compatible with the typical angles observed in several insect wings [18]. The radial vein in the model used here is of course a crude simplification of the complex patterns found in real insect wings, but it serves the purpose of separating the wing surface in
<figure><img src="image_7.png" /><figcaption>FIG. 6: Side view snapshots showing the trailing area bending for wings beating at 17.5Hz. (A) $\beta = 0°$. (B) $\beta = 10°$. (C) $\beta = 15°$. (D) $\beta = 20°$. The wing with $\beta = 20°$ does not bend on itself at no instant in the flapping cycle.</figcaption></figure>

two areas that have a dynamic equivalence to what is observed in nature. What we have called the inter-vein area functions in a similar manner to the part of the wing close to the leading edge in insects, which is rather stiff, while the trailing area deforms much more during the flapping cycle. Coupling thrust force measurements with visualisation of the wing deformation lets us explain the physical origin of the non-monotonic behaviour of the aerodynamic force production: increasing the radial vein angle makes the wing change from a regime dominated by the trailing area, with its associated strong deformations, to another regime where the inter-vein area pilots everything. It is the former case with lower angles that allows us to come back to the case of insects mentioned above: the observed optimum angle $\beta \approx 20°$ constitutes a trade-off between using the aerodynamic benefit of deformation that redirects the average force to have a stronger thrust component, and avoiding an excessive folding of the flexible wing that diminishes its effective surface. A word of caution should be said about the limitations of the present artificial flapping wings to represent the far more complex cases of real insect wings. A first point concerns the simple up-and-down movement of the wings

used here, which does not involve any of the mechanisms insects use to control wing kinematics through their thoracic muscles and hinge joints [19]. One major feature that is thus missing is the wing rotation that accompanies flapping. In addition, the two-vein wing design does not tightly control the wing shape. This leads, firstly, to camber profiles that are not representative of real wings (see, e.g., [13] for the case of hoverflies), and, secondly, to an exaggerated lack of constrain to bending of the trailing area. These issues should be considered in the design of future insect-inspired flapping robots. Ongoing work is concerned with the study of three-dimensional wing deformation dynamics with simultaneous measurement of aerodynamic forces.

## ACKNOWLEDGMENTS

This work was supported by the Agence Nationale de la Recherche and the ASTRID program through the project ANR-19-ASTR-002 (NANOFLY).

[1] R. J. Wootton, “Support and deformability in insect wings,” Journal of Zoology **193**, 447–468 (1981).

[2] W. Shyy, H. Aono, C. K. Kang, and H. Liu, *An introduction to flapping wing aerodynamics*, Vol. 37 (Cambridge University Press, 2013).

[3] T. L. Hedrick, S. A. Combes, and L. A. Miller, “Recent developments in the study of insect flight,” Canadian Journal of Zoology (2014), 10.1139/cjz-2013-0196.

[4] S. P. Sane, “The aerodynamics of insect flight,” The Journal of experimental biology **206**, 4191–4208 (2003).

[5] Z. J. Wang, “Dissecting insect flight,” Annual Review of Fluid Mechanics **37**, 183–210 (2005).

[6] M. Sun, “Insect flight dynamics: Stability and control,” Reviews of Modern Physics **86**, 615–646 (2014).

[7] D. D. Chin and D. Lentink, “Flapping wing aerodynamics: From insects to vertebrates,” Journal of Experimental Biology **219**, 920–932 (2016).
[8] L. Zhao, Q. Huang, X. Deng, and S. P. Sane, "Aerodynamic effects of flexibility in flapping wings," Journal of the Royal Society Interface 7, 485-497 (2010).

[9] X. Cai, Y. Xue, D. Kolomenskiy, R. Xu, and H. Liu, "Elastic storage enables robustness of flapping wing dynamics," Bioinspiration & Biomimetics 17, 045003 (2022).

[10] A. J. Bergou, S. Xu, and Z. J. Wang, "Passive wing pitch reversal in insect flight," Journal of Fluid Mechanics 591, 321-337 (2007).

[11] F. T. Muijres, L. C. Johansson, R. Barfield, M. Wolf, G. R. Spedding, and A. Hedenstrom, "Leading-edge vortex improves lift in slow-flying bats," Science (New York, NY) 319, 1250 (2008).

[12] M. Hamamoto, Y. Ohta, K. Hara, and T. Hisada, "Application of fluid-structure interaction analysis to flapping flight of insects with deformable wings," Advanced Robotics 21, 1-21 (2007).

[13] S. M. Walker, A. L. R. Thomas, and G. K. Taylor, "Deformable wing kinematics in free-flying hoverflies," Journal of The Royal Society Interface 7, 131-142 (2009).

[14] R. J. Bomphrey and R. Godoy-Diana, "Insect and insect-inspired aerodynamics: unsteadiness, structural mechanics and flight control," Current Opinion in Insect Science 30, 26-32 (2018).

[15] S. A. Combes and T. L. Daniel, "Flexural stiffness in insect wings I. Scaling and the influence of wing venation," Journal of Experimental Biology 206, 2979 (2003).

[16] S. A. Combes and T. L. Daniel, "Flexural stiffness in insect wings II. Spatial distribution and dynamic wing bending," Journal of Experimental Biology 206, 2989-2997 (2003).

[17] S. Heathcote, Z. Wang, and I. Gursul, "Effect of spanwise flexibility on flapping wing propulsion," Journal of Fluids and Structures 24, 183-199 (2008).

[18] A. R. Ennos, "Comparative functional morphology of the wings of diptera," Zoological journal of the Linnean Society 96, 27-47 (1989).

[19] R. J. Wootton, "Functional morphology of insect wings," Annual review of Entomology 37, 113-140 (1992).

[20] R. Wootton, "The geometry and mechanics of insect wing deformations in flight: A modelling approach," Insects 11, 446 (2020).

[21] H. V. Phan and H. C. Park, "Mechanisms of collision recovery in flying beetles and flapping-wing robots," Science 370, 1214-1219 (2020).

[22] C. J. C. Rees, "Form and function in corrugated insect wings," Nature 256, 200-203 (1975).

[23] S. R. Jongerius and D. Lentink, "Structural analysis of a dragonfly wing," Experimental Mechanics 50, 1323-1334 (2010).

[24] H. Tanaka, J. P. Whitney, and R. J. Wood, "Effect of Flexural and Torsional Wing Flexibility on Lift Generation in Hoverfly Flight," Integrative and Comparative Biology 51, 142-150 (2011).

[25] G. Cousin and O. Béthoux, "The steps towards an inconspicuous vein fusion documented in Stenosmylinae forewings (Neuroptera: Osmylidae)," Organisms Diversity & Evolution (2016).

[26] P. Wu, P. Ifju, and B. Stanford, "Flapping Wing Structural Deformation and Thrust Correlation Study with Flexible Membrane Wings," AIAA journal. (2012), 10.2514/1.J050310.

[27] C. K. Kang, H. Aono, C. E. S. Cesnik, and W. Shyy, "Effects of flexibility on the aerodynamic performance of flapping wings," Journal of Fluid Mechanics 689, 32-74 (2011).

[28] G. C. H. E. de Croon, M. Perçin, B. D. W. Remes, R. Ruifsink, and C. De Wagter, "Mechanical design and materials," in *The DelFly: Design, Aerodynamics, and Artificial Intelligence of a Flapping Wing Robot* (Springer Netherlands, Dordrecht, 2016) pp. 31-39.

[29] A. Magnan, *La locomotion chez les animaux: Vol. 1: Le vol des insectes* (Hermann & Cie, Éditeurs Paris, 1934).

[30] C. H. Greenewalt, *Dimensional relationships for flying animals* (Smithsonian miscellaneous collections, 1962).

[31] R. Dudley and C. P. Ellington, "Mechanics of forward flight in bumblebees: I. kinematics and morphology," Journal of Experimental Biology 148, 19-52 (1990).

[32] T. L. Daniel and S. A. Combes, "Flexible wings and fins: Bending by inertial or fluid-dynamic forces?" Integrative and Comparative Biology 42, 1044 (2002).

[33] B. Thiria and R. Godoy-Diana, "How wing compliance drives the efficiency of self-propelled flapping flyers," Physical Review E 82, 015303 (2010).

[34] K. S. Wu, J. Nowak, and K. S. Breuer, "Scaling of the performance of insect-inspired passive-pitching flapping wings," Journal of The Royal Society Interface 16, 20190609 (2019).

[35] E. de Langre, "Effects of wind on plants," 40, 141-168 (2008).

[36] D. Ishihara, Y. Yamashita, T. Horie, S. Yoshida, and T. Niho, "Passive maintenance of high angle of attack and its lift generation during flapping translation in crane fly wing," Journal of Experimental Biology 212, 3882-3891 (2009).

[37] S. Ramananarivo, R. Godoy-Diana, and B. Thiria, "Rather than resonance, flapping wing flyers may play on aerodynamics to improve performance," Proceedings of the National Academy of Sciences 108, 5964-5969 (2011).

[38] C. A. Schneider, W. S. Rasband, and K. W. Eliceiri, "NIH Image to ImageJ: 25 years of image analysis," Nature methods 9, 671-675 (2012).

[39] K. Mazaheri and A. Ebrahimi, "Experimental investigation of the effect of chordwise flexibility on the aerodynamics of flapping wings in hovering flight," Journal of Fluids and Structures 26, 544-558 (2010).

[40] T. A. Nguyen, H. V. Phan, T. K. L. Au, and H. C. Park, "Experimental study on thrust and power of flapping-wing system based on rack-pinion mechanism," Bioinspiration & Biomimetics 11, 046001 (2016).