# Design and Optimisation of a Vibrating Wing Insect-Size Air Vehicle with Lumped Parameter Models and Compliant Links

Marguerite de La Bigne, Mathieu Colin, Eric Cattan, Sofiane Ghenna, Marie Zwingelstein, Sébastien Grondel, Olivier Thomas

► To cite this version:

Marguerite de La Bigne, Mathieu Colin, Eric Cattan, Sofiane Ghenna, Marie Zwingelstein, et al. Design and Optimisation of a Vibrating Wing Insect-Size Air Vehicle with Lumped Parameter Models and Compliant Links. Journal of Bionic Engineering, 2025, 10.1007/s42235-025-00761-1. hal-05310578

HAL Id: hal-05310578

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

Submitted on 13 Oct 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
# Design and Optimisation of a Vibrating Wing Insect-Size Air Vehicle with Lumped Parameter Models and Compliant Links

Marguerite de La Bigne<fnref n="1" /> ID · Mathieu Colin<fnref n="1" /> · Éric Cattan<fnref n="2" /> · Sofiane Ghenna<fnref n="2" /> · Marie Zwingelstein<fnref n="2" /> · Sébastien Grondel<fnref n="2" /> · Olivier Thomas<fnref n="1" />

Received: 24 January 2025 / Revised: 3 July 2025 / Accepted: 10 August 2025
© The Author(s) 2025

## Abstract

This article presents the design of a microfabricated bio-inspired flapping-wing Nnano Aaerial Vvehicle (NAV), driven by an electromagnetic system. Our approach is based on artificial wings composed of rigid bodies connected by compliant links, which optimise aerodynamic forces through replicating the complex wing kinematics of insects. The originality of this article lies in a new design methodology based on a triple equivalence between a 3D model, a multibody model, and a mass/spring model (0D) which reduces the number of parameters in the problem. This approach facilitates NAV optimisation by using only the mass/spring model, thereby simplifying the design process while maintaining high accuracy. Two wing geometries are studied and optimised in this article to produce large-amplitude wing motions (approximately 40°), and enabling flapping and twisting motion in quadrature. The results are validated thanks to experimental measurements for the large amplitude and through finite element simulations for the combined motion, confirming the effectiveness of this strategy for a NAV weighing less than 40 mg with a wingspan of under 3 cm.

**Keywords** Resonant motion · Micro/nano electromechanical system (M/NEMS) · Modal combination · Design of structural vibration · Optimisation · Compliant links

## 1 Introduction

Insects are characterised by their agility in flight, including hovering, backward and forward motions, as well as sudden accelerations. Moreover, their diversity provides abundant sources of bio-inspiration to address scientific issues, particularly in understanding their flight mechanisms. It has been demonstrated that lift force results both from a wide range of unconventional aerodynamic mechanisms [1] and from specific kinematics induced by the thorax and the flexible characteristics of the wings [2]. Since the kinematics can be approximated by a combination of twisting and flapping motions [3–5], some research projects use this principle to build unmanned air vehicles (UAVs) that mimic the kinematics of an insect wing [6–8].

In this article, we focus on the smallest category of UAVs, those weighing less than 10 g, usually refered as nano-air vehicles (NAVs) [9–14]. In the literature, there are two main concepts of those vibrating wings NAVs. The first is based on rigid wings linked with hinges to a thorax and excited through a four-bar system to amplify their motion [15]. Two methods of manufacture were reported, proposing

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

Olivier Thomas
olivier.thomas@ensam.eu

Marguerite de La Bigne
marguerite.de_la_bigne@ensam.eu

Mathieu Colin
mathieu.r.colin@gmail.com

Éric Cattan
Eric.cattan@uphf.fr

Sofiane Ghenna
Sofiane.Ghenna@uphf.fr

Marie Zwingelstein
marie.zwingelstein-colin@uphf.fr

Sébastien Grondel
sebastien.grondel@uphf.fr

<footnote n="1">¹ Arts et Métiers Institute of Technology, LISPEN, HESAM Université, 59000 Lille, France</footnote>

<footnote n="2">² Univ. Polytechnique Hauts-de-France, CNRS, Univ. Lille, UMR 8520-IEMN-Institut d'Electronique de Microélectronique et de Nanotechnologie, 59313 Valenciennes, France</footnote>
two families of 60 mg, 3 cm wingspan NAV prototypes. The first, called “smart composite microstructures”, is based on laser micromachining of composite materials and polymers, which are then stacked and cured together [10]. The second, called monolithic fabrication, involves the laser cut of a carbon/polyimide composite which is folded and locked to form the desired three-dimensional shape [9]. The second concept, considered in the present article, relies on a fully flexible monolithic thorax/wing structure, made in clean room using microfabrication technologies [16, 17], giving a smaller and lighter structure of 22 mg and 2.2 cm wingspan [14]. It also takes advantage of a resonant motion to increase the vibration amplitude, based on a combination of flapping and twisting resonances to create an insect bio-inspired wing kinematics [13]. Comparing those two NAV concepts, the use of microfabrication enables very light and small prototypes to be produced as well as massive production at low cost. However, the design of a fully flexible structure relies on a large number of parameters, leading most of the time to a complex and time consuming optimisation process [14].

In this article, a new design methodology for microfabricated resonant NAVs is proposed. It is based on designing the fully flexible structure of the NAV wings such that the deformation is localised in small size regions, in order to create a global deformed shape of the wings in the form of rigid bodies (naturally called rib throughout the paper) articulated with compliant links. This new methodology, leading to new architectures for the NAVs under concern, has three main advantages, described hereafter. The first is that the optimisation of the complete structure can be achieved using simple mass/spring vibrating systems, denoted here by “zero-dimensional” (0D) models, with very few degrees of freedom, thanks to the triple equivalence illustrated in Fig. 1. If the real system (model 0) is composed of a thick rib with a compliant hinge made of a thin elastic part, it is reasonable to consider the rib as a rigid body and the compliant link as a perfect hinge with a torsional spring (model 1). Then, since this system has only one (rotation) degree of freedom, its dynamics is perfectly equivalent to the one of a mass/spring system (model 2). The main idea of this article, originally proposed in the PhD of the second author [18], is to validate and apply this methodology to model the dynamics of NAVs using simple mass/spring systems with very few degrees of freedom. Compliant mechanisms are widely used in many applications thanks to advantages such as monolithic design, increased motion precision, no errors

from assembly, simplified manufacturing process, reduced weight and less maintenance [19, 20], especially in robotics [21] and also in the present subject of NAVs [9, 10]. However, most applications of compliant mechanisms are confined to static analyses, with very few addressing their dynamical behaviour in the design process. An interesting review on works dealing with dynamical compliant mechanisms is proposed in [22], noticing that the resonance frequencies of a mechanism may have importance in some applications and that they can be included in the design process using several strategies, from full finite-element numerical computations to analytical and semi-analytical approaches. However, to the knowledge of the authors, the strategy proposed in the present article, based on the triple equivalence of Fig. 1, is original, especially for the design of NAVs. The second advantage of our architecture is that the thinness of the bodies is localised at the compliant links, enabling a very high amplitude of vibration without damaging the structure [23]. This point will be addressed in Sect. 2.2. Then, the third advantage is that the aerodynamic force produced by a rigid wing attached to a compliant link appears slightly greater than the one observed for an equivalent flexible wing. This is quantitatively shown in Sect. 2.3 for a simple bending motion, and assumed to be still valid for a more realistic combined bending/twisting motion of a wing.

The main goal of this paper is to demonstrate the benefits of these new NAV architectures, that include flexible links and rigid bodies, in term of design methodology and performance. In parallel, these new NAV architectures are also the pretext for exploring and validating equivalences between fully flexible structures, multibody models and, at the end, 0D models. The paper is divided into four main parts. The first one (Sect. 2) introduces the basic architectures of our NAV, the vibration modes combination used to generate lift and proves the benefit, in term of lift force, of using hinged rigid bodies for the wings. Then, Sect. 3 specifically addresses the equivalence between a flexible structure, a multibody model and a 0D model, by considering several architectures of increasing complexity. Finally, based on the equivalence, 0D models are used in Sect. 4 to optimise the design of several NAV architectures, to obtain a final optimised configuration described in Sect. 5.

<figure><img src="image_3.png" /><figcaption>**Fig. 1** Equivalent models for NAV design</figcaption></figure>
# 2 Basic Concepts of our NAVs

## 2.1 Basic Architecture and Working Actuation Principle

The NAV basic structure, shown in Fig. 2, is composed of two wings with their roots connected to the free end of a central beam, at which a magnet drives the system in vibrations, thanks to a coil attached to a thorax, the two latter being not shown in Fig. 2. Each wing is composed of two ribs, assumed rigid, the leading edge and the trailing edge, connected to the free end of the central beam by compliant links, made of two short and very thin parts. A very thin membrane is deposited between the leading edge and the trailing edge to create the lifting surface coupled with the fluid and thus the thrust force of the NAV. The central beam can be considered as the motor of the NAV. Fixed at its origin, it is designed to have enough flexibility to allow the electromagnetic force on the magnet to deform the beam, inducing a resonant motion of the two wings.

The NAVs under study are manufactured in polymer by dedicated microfabrication techniques, as described in [16, 17]. The fabrication principle is based on the photosensitive characteristic of the SU-8 epoxy-based resin that is used, which polymerises when exposed to Ultra-Violet (UV) radiation and becomes insoluble to liquid developers. The NAVs are manufactured on 3 inch diameter silicon wafers, on which eight NAVs can be placed. At the start of manufacturing, a sacrificial layer in aluminium is applied to the wafer to prevent the prototypes from sticking to the silicium. Then, successive layers of SU-8 between 15 and 200 µm are deposited with a precise geometry defined by several masks, to form the structure of the NAV. After removing all non-isolated SU-8 parts with a liquid developer, Parylene C is pulverised onto the wafer to form the NAV's wings, resulting in thin membranes of 400 nm. A mask is also used to keep Parylene C only in the two triangle shaped areas of the wing membrane, between the leading and the trailing edges.

Many studies of insect flight [1, 2, 24–26] have shown that it is related to complex aerodynamic phenomena linked to the particular kinematics of the wings. Basically, it is usually described as the combination of four successive motions: the downstroke, the supination, the upstroke, and the pronation (see [3, 13, 27]). The main characteristics of these motions can be approximated by the simple combination of a flapping motion in phase quadrature with a twisting motion, i.e. one being maximal when the other is minimal. Based on this result and as proposed in [13, 14, 28], the NAV design principle of this work focuses on optimising a flexible structure in a manner that the flapping and the twisting motions are both induced by its resonant dynamic behaviour. More precisely, the flapping motion is created by a resonant mode where the leading and trailing edges of the wing are in phase (Fig. 3a), whereas the twisting motion corresponds to a resonant mode where the two edges are in phase opposition (Fig. 3b). To obtain such a kinematics, the elastic characteristics of the wings (in our case, we adjust only the geometry since the material, the SU-8 polymer, is the same for all the structural parts of the NAV, except for the Parylene C membranes, assumed to have a negligible structural effect) must be optimised in order to place the natural frequencies of the flapping and twisting modes close to each other. Then, driving the NAV with a sine signal at a particular frequency in the middle of the two resonances creates the desired quadrature motion, as shown in [13] and in Fig. 28.

## 2.2 Large Amplitude Motion

Insects are capable of achieving very large amplitude angles with their wings, this property being essential for generating lift. As a preliminary, we consider here a simple “toy” prototype, shown in Fig. 4a–c, with its wing restricted to a leading edge only, without trailing edge and membrane. The leading edge is connected to a flexible central beam

<figure><img src="image_4.png" /><figcaption>Fig. 2 a Basic geometry and b photograph of a NAV prototype under study</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 4** a–c “T” prototype with only a leading edge; “thick” means “thickness”. d Sketch of the compliant link in large amplitude deformation</figcaption></figure>

with a thin compliant link, to exhibit significant angular displacements.

Due to the large amplitude vibrations in the compliant link, it is crucial to avoid any plasticity and possible cracking. This is ensured by limiting the maximum stress $\sigma$ in the compliant link such that it remains under the yield limit $\sigma_e$ of the material. Following the results of Appendix A, the compliant link is mainly submitted to a bending moment $M_L$ at its tip (see Fig. 4d). The equilibrium of any cross section of the beam shows that the bending moment is uniform over its length: $M(y) = M_L$, with y the axial coordinate. If the compliant link is considered as an elastic beam of length $L_l$, with a rectangular cross section of thickness $h_l$ and width $b_l$, with an Euler–Bernoulli kinematics, its constitutive law is $M = EI d\theta/dy$, with $\theta$ the cross-section rotation, E the Young’s modulus of the material and $I = b_l h_l^3 / 12$ the second moment of area of the cross section. This relation is valid even for extreme bending with strong geometrical

nonlinearities effects (see e.g. [29, 30]). By integration of $M(y) = M_L$, the rotation of the tip’s cross section is:

$$ \theta_L = \frac{M_L L_l}{EI} \qquad (1) $$

Notice that the compliant link must be thin (smaller than $\approx 1/10$) such that the beam theory is valid. The maximum stress appears on the upper and lower surfaces of the beam and is equal to $\sigma = h_l M(y) / (2I) = h_l M_L / (2I)$. Assuming that $\sigma < \sigma_e$ with $\sigma_e$ the yield limit of the material, the maximal rotation angle of the tip’s cross-section is:

$$ \theta_m = \frac{2L_l \sigma_e}{h_l E} \qquad (2) $$
<table><caption>**Table 1** Dimensions of the “T” prototype of Fig. 4</caption><thead><tr><th></th><th></th><th>Leading edge</th><th>Compliant link</th><th>Central beam</th></tr></thead><tbody><tr><td>Length</td><td>mm</td><td>L<sub>1</sub></td><td>L<sub>l1</sub></td><td>L<sub>3</sub></td></tr><tr><td>Width</td><td>mm</td><td>b<sub>1</sub></td><td>b<sub>l1</sub></td><td>b<sub>3</sub></td></tr><tr><td>Thickness</td><td>μm</td><td>h<sub>1</sub></td><td>h<sub>l1</sub></td><td>h<sub>3</sub></td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 5** Snapshots of the large amplitude motion of a “T” prototype of Fig. 4a–c</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6** Comparison of NAV’s wings made with **a** flexible ribs or **b** rigid bodies</figcaption></figure>

According to the SU-8 data-sheet, $\sigma_e/E \approx 1/33$ [23]. With the dimensions of the prototype presented in Table 1, Eq. (2) leads to a theoretical value of the maximum angle of 38.4°.

To verify these results, the “T” prototype of Fig. 4 was manufactured by microfabrication technology [14, 16, 17]. It has only a leading edge and no membrane. It is fixed by a clamp at the origin of its central beam and powered by an external coil (not shown). A high-speed camera (a Phantom Micro LC311) enables taking some snapshots of the motion, as shown in Fig. 5. An angle $\theta_L$ of approximately 40° was measured without damage, thus validating the order of magnitude of the theoretical calculations and showing the ability of our compliant link to enable large amplitude motion.

## 2.3 Aerodynamic Force

Unlike NAVs, insects can actively control the deformation of their wings with their muscles [5], allowing them to generate significant lift. However, for artificial wings, studies have shown that the curvature of the wing has little influence on lift generation [31] and that a rigid wing generates more lift than a flexible one [32–35]. Since this article specifically addresses rigid wings with compliant links, this section is dedicated to estimating and comparing the aerodynamic force created by rigid and flexible wings.

The two simple models of Fig. 6 are considered. For the first one, the wing is considered flexible and modeled as a cantilever beam oscillating on its first vibration mode shape. For the second model, the wing is considered rigid and connected to the ground with a hinge and a torsional spring. The transverse displacement field of both models is denoted by $w(x, t)$, with $x$ the axial coordinate and $t$ the time. The amplitude of the motion at the tip is considered identical for comparison purpose and denoted by $w_L(t) = w(L, t)$, with $L$ the length of the two wings.

The aerodynamic force vector is computed using the following formula [36]:

$$F = \frac{\rho_f c_0}{2} \iint_S ||v||(v \cdot n)n \, dS, \quad (3)$$

where $v = v(M, t)$ is the velocity vector of any point $M$ of the wing surface $S$, $n = n(M, t)$ is the unitary normal vector to this surface, is the dot product of two vectors, $|| \circ ||$ is the norm of vector $\circ$, $\rho_f$ is the fluid density and $c_0$ is an aerodynamic dimensionless coefficient. This formula assumes that the local aerodynamic force has three properties: (i) it is proportional to the square of the velocity of the fluid with respect to the wing; (ii) it is applied on the local
surface orthogonal to the fluid flow, thanks to the scalar product; (iii) its direction is orthogonal to the wing surface.

We consider only small displacements, neglecting geometrical nonlinearities, such that the displacement field direction is $\mathbf{y}$, orthogonal to the direction $\mathbf{x}$ of the wings at rest (undeformed). Assuming a modal motion for $w(x, t)$, one obtains:

$$w(x, t) = \Phi(x)w_L(t), \Rightarrow v(x, t) = \dot{w}(x, t)\mathbf{y} = \Phi(x)\dot{w}_L(t)\mathbf{y}, \quad (4)$$

where $\Phi(x)$ is the shape of the displacement field, normalised such that $\Phi(L) = 1$. For the two cases of Fig. 6, one has (see e.g. [37] for the normalised first mode shape of a cantilever beam):

$$\text{Flexible wing: } \Phi(x) = \frac{1}{2} [\cos \beta_1 \bar{x} - \cosh \beta_1 \bar{x} + \frac{\sin \beta_1 - \sinh \beta_1}{\cos \beta_1 + \cosh \beta_1} (\sin \beta_1 \bar{x} - \sinh \beta_1 \bar{x})], \quad (5a)$$

$$\text{Rigid wing: } \Phi(x) = \bar{x}, \quad (5b)$$

with $\bar{x} = x/L$ and $\beta_1 = 1.8751$.

If the position vector of any point of the wing's deformed middle line is $\mathbf{p} = x\mathbf{x} + wy$, the tangent vector is $\mathbf{t} = \partial\mathbf{p}/\partial x$, such that the unitary normal vector is:

$$\mathbf{n}(x, t) = \frac{-w'(x, t)\mathbf{x} + \mathbf{y}}{\sqrt{1 + w'^2(x, t)}} \approx \mathbf{y}, \quad (6)$$

where $w'(x, t) = \partial w/\partial x \ll 1$ since we consider a small amplitude motion and no geometrical nonlinear effects. The axial component (along $x_0$) of the force is thus neglected since it would be directly related to geometrical nonlinearities. We then restrict the discussion to the transverse (and main) component of the force, that reads:

<figure><img src="image_8.png" /><figcaption>Fig. 7 The two architectures of the NAVs under study: a the trailing edge (TE) is connected to the leading edge (LE) and b both edges are connected to the central beam (CB)</figcaption></figure>

$$F \approx F_y \mathbf{y} = \frac{\rho_f c_0 b L}{2} \dot{w}_L^2 \mathbf{y} \int_0^1 \Phi^2(\bar{x}) \, d\bar{x}, \quad (7)$$

with b the width of the wing (the dimension along direction $z$).

In the two considered cases of Fig. 6, only the integral term in the above equation changes. Its value is exactly 1/4 in the flexible wing case and 1/3 for the rigid wing (a consequence of Eqs. 5a, b).

As a conclusion, we recover here the announced result: at equal tip amplitude value, a rigid wing generates more transverse aerodynamic force in comparison to a fully flexible wing. With the two configurations of Fig. 6, the ratio is, for a small amplitude motion:

$$\frac{F_{\text{rw}}}{F_{\text{fw}}} = \frac{4}{3} \approx 1.33, \quad (8)$$

with subscripts “rw” and “fw” meaning respectively “rigid wing” and “flexible wing”. Even if this result is here demonstrated considering a very simple (bending/flapping) motion, we assume that it is also valid for a more complex and realistic motion of an insect wing, that involves an additional twisting motion in combination with the flapping (see Sect. 2.1). Doing the same aerodynamic force evaluation on this kind of more complex motion is left out of the scope of the present article and addressed elsewhere [36].

## 3 Equivalence Between 0D Models, Multibody Models and a Compliant Links Structure

In this section, we address the triple equivalence illustrated in Fig. 1 in the case of more complex models related to the architectures of our NAVs, shown in Fig. 7. Two slightly different NAV's architectures are under study. They both include a compliant central beam, with a magnet glued at

<figure><img src="image_9.png" /></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 8 Overview of the four models studied in Sect. 3</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9 a Real system of a single discontinuous beam (model A0), b multibody model with a torsional spring at the rigid body end (model A1₁) and c multibody model with the torsional spring at the middle of the compliant link's length (model A1₂)</figcaption></figure>

its tip, to realise the electromagnetic actuation of the NAV. For architecture 1 (Fig. 7a, see the detailed view), the trailing edge is connected to the leading edge: it is directly inspired by the geometry of the NAVs which are already taking off [9, 10, 38]. For the second architecture (Fig. 7b), both the trailing and the leading edges are attached to the central beam. The compliant central beam can be modelled by a mass/spring system (see Sect. 5.1). Then, considering only one wing and the equivalence illustrated in Fig. 1, the two edges with their compliant link can be modelled by two additional mass/spring systems, in series for architecture 1 and in parallel for architecture 2.

Before using those 3 degrees of freedom 0D models in Sect. 4 to optimise the dynamics of the two NAVs architectures, the equivalences of Figs. 1 and 7 are addressed quantitatively in this section. It leads to considering four families of models, by increasing complexity, as summarised in Fig. 8.

## 3.1 A Single Articulated Beam

This section addresses quantitatively the equivalence of Fig. 1 for a simple rib/hinge system, as sketched in Fig. 9. The real system (Fig. 9a, named model A0) is composed of a single discontinuous beam of rectangular cross section,

composed of two segments. The first, the compliant link, has a length $L_l$, a width $b_l$ and a thickness $h_l$. The second, the rib, has a length $L$, a width $b$ and a thickness $h$. It is thicker and longer, $h \gg h_l$ and $L \gg L_l$, such that it can be considered as a rigid body. It also enables to assume the deformation of the compliant link as quasi-static, neglecting its inertia. The equivalent multibody models (Fig. 9b,c) are composed of a rigid beam, corresponding to the rib, and a hinge/torsional spring assembly, to replace the compliant link. Two slightly different cases for the multibody model will be studied: model A1₁, where the torsional spring is at the origin of the rigid body (Fig. 9b) and model A1₂, where it is on the middle of the compliant link's length (Fig. 9c).

Let us first consider the multibody models A1₁ and A1₂. We denote by $\Gamma$ the stiffness of the torsional spring (in N/m/rad) and by $m, L$ and $J = mL^2/3$ the mass, length and moment of inertia of the rib, respectively. The natural frequency of model A1₁ is

$$ \omega_{rb} = \sqrt{\frac{\Gamma}{J}} = \sqrt{\frac{3\Gamma}{mL^2}} \qquad (9) $$

If now the length of the compliant link $L_l$ is taken into account (Fig. 9c), neglecting its inertia and using the
Huygens–Steiner theorem, the new moment of inertia $J'$ is therefore given by,

$$J' = \frac{m(4L^2 + 6LL_l + 3L_l^2)}{12}, \qquad (10)$$

leading to the following expression of the natural frequency of model A12:

$$\omega'_{rb} = \omega_{rb} \sqrt{\frac{1}{1 + \frac{3L_l}{2L} + \frac{3L_l^2}{4L^2}}.} \qquad (11)$$

For the value of $\Gamma$, the compliant link is considered as quasi-static and in pure bending such that (see [39] and Appendix A):

$$\Gamma = \frac{EI_l}{L_l}, \qquad (12)$$

where $I_l = b_l h_l^3 / 12$ is the compliant link's second moment of area and $E$ is the material's Young's modulus. Introducing Eq. (12) in Eq. (9) leads to a natural frequency for the multibody model of Fig. 9b,

$$\omega_{rb} = \frac{1}{2} \sqrt{\frac{E}{\rho}} \sqrt{\frac{b_l h_l^3}{L_l b h L^3}} \qquad (13)$$

where $\rho$ is the density of the material.

The equivalence of the models A0, A11 and A12 is based on the two assumptions of rigid body and quasi-static deformations, respectively for the rib and the compliant link. To evaluate the validity domain of those assumptions, we model the real system (model A0) as an Euler–Bernoulli beam with variable section. Introducing the time $t$ and defining $x$ as the position along the beam, the transverse displacement $w(x, t)$ satisfies the following partial differential equation, $\forall x \in [0, L + L_l], t \ge 0$ [40]:

$$\rho S \ddot{w} + (EIw'')'' = p(x), \qquad (14)$$

where $p(x)$ is the load per unit length, $\circ = \partial \circ / \partial t$, $\circ' = \partial \circ / \partial x$ are the time and space derivatives of $\circ(x, t)$, whereas $S(x) = b_b(x)h_b(x)$ is the area of the section and $I(x) = b_b(x)h_b(x)^3 / 12$ is the second moment of area. The problem depends on 8 parameters: 6 geometric ones ($L, b, h, L_l, b_l, h_l$) and 2 material ones ($E, \rho$). Equation (14) is rewritten in a dimensionless form to identify the independent parameters and simplify the analysis. The following dimensionless quantities, denoted by overbars, are introduced:

$$\begin{gathered} \bar{x} = \frac{x}{L_l}, \quad \bar{b} = \frac{b}{b_l}, \quad \bar{h} = \frac{h}{h_l}, \quad \bar{L} = \frac{L}{L_l}, \quad \bar{w} = \frac{w}{L_l}, \\ \bar{t} = \omega_0 t, \quad \bar{\omega} = \frac{\omega}{\omega_0}, \quad \omega_0 = \frac{h}{L_l^2} \sqrt{\frac{E}{12\rho}}. \end{gathered} \qquad (15)$$

Introducing Eq. (15) into Eq. (14) and enforcing $p = 0$, one shows that the eigenmodes $(\bar{\omega}_i, \Phi_i(\bar{x}))$, $i \in \mathbb{N}^*$, are solutions of the following dimensionless equation:

$$\bar{b}\bar{h}\ddot{\bar{w}} + (\bar{b}\bar{h}^3\bar{w}'')'' = 0 \Rightarrow (\bar{b}\bar{h}^3\Phi_i'')'' - \bar{\omega}_i^2\bar{b}\bar{h}\Phi_i = 0. \qquad (16)$$

These equations imply that the modes $(\omega_i, \Phi_i(x))$ depend on three independent geometric parameters only: $(\bar{h}, \bar{b}, \bar{L})$, defined in Eq. (16).

To evaluate the range of validity of equivalence of Fig. 9 between models A0, A11 and A12, Eq. (16) is solved with a finite element analysis (FEA, with the Abaqus software and quadratic Timoshenko B32 beam finite-elements) for a beam that respects $L/h > 20$ (to fullfill the beam model assumptions), in order to compute the first natural frequency $\bar{\omega}_1$ of model A0. Choosing beam elements is based on their less than 1% error in the first mode frequency compared to a simulation with 3D elements. Then, $\bar{\omega}_1$ is compared to the natural frequency of the two multibody models A11 and A12 thanks to Eqs. (9) and (11). Figure 10a–c shows the comparison as a function of parameters $(\bar{h}, \bar{b}, \bar{L})$. The red zone in each figure denotes a tolerance margin, representing a range of ±5% of the value of $\bar{\omega}_1$. In addition, the real system first mode shape $\Phi_1(\bar{x})$ is plotted in Fig. 10d–f, for several values of the three ratios $(\bar{h}, \bar{b}, \bar{L})$.

As a preliminary comment, since for a beam with homogeneous rectangular cross section, the natural frequencies do not depend on the width $b$ of the beam, one may have thought that it would have also been the case for a beam of constant height $h$ and varying width $b(x)$. However, in the latter case, the natural frequencies depend on $b(x)$, as shown by Eq. (16), by Fig. 10a (with $\bar{h} = 1$, $\bar{\omega}_1 \neq 1.875^2$ since $\bar{b} = 1.25 \neq 1$) and also in [41].

Figure 10b, e illustrate that the influence of the ratio $\bar{b}$ on models A11 and A12 is minimal. Similarly, Fig. 10c shows that for model A12, the ratio $\bar{L}$ has no impact, even though Fig. 10f displays a rigid body deformation for $\bar{L} = 2$. It is ratio $\bar{h}$ that is the most critical for both models, as shown in Fig. 10a, d: when $\bar{h}$ is low, the assumption of a rigid body for the rib is no longer valid due to its deformation.

As a main conclusion, to remain under 5% error and have a reasonable equivalence between the multibody model and the real system, for model A11 (where the torsional spring is at the origin of the rigid body), it is necessary to keep $\bar{h} = h/h_l > 7.5$ and $\bar{L} = L/L_l > 15$ (see Fig. 10a, c). For model A12 (where the torsional spring is at the middle of the compliant link's length) $\bar{h} = h/h_l > 3$ must be verified.
<figure><img src="image_12.png" /></figure>

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

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

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

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

<figure><img src="image_17.png" /><figcaption>**Fig. 10** (First row): First natural dimensionless frequency $\bar{\omega}_1$ of the two segment beam of Eq. (16) as a function of the dimensionless parameters $\bar{L} = L/L_l$, $\bar{h} = h/h_l$ and $\bar{b} = b/b_l$: **a** $\mathbf{a}\bar{\omega}_1 = f(\bar{h})$ with $\bar{L} = 16.4$ and $\bar{b} = 1.25$; **b** $\mathbf{b}\bar{\omega}_1 = f(\bar{b})$ with $\bar{L} = 16.4$ and $\bar{h} = 12.5$; **c** $\mathbf{c}\bar{\omega}_1 = f(\bar{L})$ with $\bar{b} = 1.25$ and $\bar{h} = 12.5$. The blue line is for model A11, the purple line for model A12 and the red line is for FEA. The</figcaption></figure>

red shaded area illustrates a variation of ±5% of the FEA solution. (second row): First mode shape $\Phi_1(x)$ for two extreme selected values of $(\bar{L}, \bar{h}, \bar{b})$ corresponding to the above plot: **d** $\bar{L} = 16.4, \bar{b} = 1.25, \bar{h} \in \{2, 8\}$; **e** $\bar{L} = 16.4, \bar{h} = 12.5, \bar{b} \in \{2, 20\}$; **f** $\bar{b} = 1.25, \bar{h} = 12.5, \bar{L} \in \{2, 16\}$. The dots shows the limit between the compliant link and the rib

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

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

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

<figure><img src="image_21.png" /><figcaption>**Fig. 11** Comparison between FEA and experiments for the first two mode shapes of the “T” prototype of Fig. 4. **a**, **b** operational deflection shapes at the first two resonance frequencies; **b** profile of the first mode</figcaption></figure>

shape $\Phi_1(x)$ on the dashed-dotted red line shown in (**a**); **c** idem for the second mode shape $\Phi_2(x)$

So, using these assumptions, it is possible to model the NAV wing by a multibody model.

The above numerical results are experimentally verified in Fig. 11, showing the experimental first two natural modes of the “T” prototype of Fig. 4, obtained with a scanning

laser vibrometer (Polytec PSV-400 with a close-up unit composed of a lens enabling scanning small objects).

They are obtained as the frequency and deformed shape of the structure at the maxima of the frequency response functions, shown in Fig. 12. Figures 11b, c and 12 show a
<figure><img src="image_22.png" /><figcaption>**Fig. 12** Comparison between FEA and experimental measurement of the FRF measured at the magnet (in blue) and at the tip of the leading edge (in orange). Solid lines correspond to the FEA simulation, dotted lines to experimental measurements, and dashed lines represent an FEA simulation with corrected average thicknesses (from 14 to 13µm for the compliant joints, and from 98 to 95µm for the central beam). The right ordinate axis corresponds to the FEA amplitude, and the left one to the experimental measurements. The damping used for the simulation is the one measured experimentally</figcaption></figure>

comparison with numerical simulation obtained by FEA. A good agreement is obtained, with discrepancies explained by a slight variability of the thickness of the members of the order of 3%, verified by careful profilometry measurements. Since the structure deformation is mainly bending a small change of the thickness results in a huge change in the stiffness and thus in the natural frequencies. Figure 12 also shows the FEA frequency response function with a corrected (though uniform) thickness, showing an excellent agreement between theory and experiments.

## 3.2 Two Articulated Beams

Before addressing the equivalence between a complete NAV and a three degrees of freedom multibody model in Sects. 3.3 and 3.4, the simple case of two articulated beams is considered here as a preliminary. As shown in Fig. 13, the system under study is composed of a central beam and a rib, inclined with an angle $\beta$, playing the role of either the leading ($\beta = 0$) or trailing ($\beta \neq 0$) edge. This approach is

useful to verify the equivalence between the 0D model and the real system by using the multibody model, with a geometry close to that of a complete NAV.

The central beam (in blue) has a mass $M_3$, a center of mass $G_3$ in its middle and a length $L_3$. It is hinged at point $O_3$ with respect to the ground and rotates around the axis ($O_3, y_0$) by an angle $\gamma(t) = (\hat{z}_0, \hat{z}_3)$, in parallel with a torsional spring of stiffness $\Gamma_3$. The rib (in purple) has a mass $M_2$, a center of mass $G_2$ in its middle and a length $L_2$. It forms a fixed angle $\beta = (\hat{y}_\beta, \hat{y}_0)$ along the axis ($O_2, x_3$) (this angle is the same as that in Fig. 2). It is hinged at point $O_2$ and rotates with respect to the central beam around the axis ($O_2, z_\beta$) by an angle $\alpha(t) = (\hat{y}_2, \hat{y}_\beta)$, in parallel with a torsional spring of stiffness $\Gamma_2$. One has $O_3G_3 = L_3/2 z_3, O_3O_2 = L_3 z_3 + c_3 y_0$ and $O_2G_2 = L_2/2 y_2$. Applying Lagrange's equations and linearising them around the equilibrium position ($\gamma = \alpha = 0$) (see Appendix B.2 for details), one obtains:

$$ \left[ \left( \frac{M_3}{3} + M_2 \right) L_3^2 + M_2 \left( \frac{L_2^2 \sin^2(\beta)}{3} - L_2 L_3 \sin(\beta) \right) \frac{L_2 M_2 \left( \frac{-2L_2 \sin(\beta) + 3L_3}{6} \right)}{L_2 M_2 \left( \frac{-2L_2 \sin(\beta) + 3L_3}{6} \right)} \right] \begin{bmatrix} \ddot{\gamma} \\ \dot{\alpha} \end{bmatrix} + \begin{bmatrix} \Gamma_3 & 0 \\ 0 & \Gamma_2 \end{bmatrix} \begin{bmatrix} \gamma \\ \alpha \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}. \quad (17) $$

The equivalent 0D model (model B2) shown in Fig. 13 is composed of a mass $m_3$ and a spring $k_3$ representing the central beam, while the rib is modelled by a mass $m_2$ and a

spring $k_2$. The displacements with respect to the equilibrium position of masses $m_3$ and $m_2$ are denoted by $u_3(t)$ and $u_2(t)$ respectively, whereas $u_\alpha(t) = u_2(t) - u_3(t)$ is the
<figure><img src="image_23.png" /><figcaption>Fig. 13 Real system with compliant links (model B0), multibody model (model B11) and their equivalent 0D model (model B2), for a NAV with only one trailing edge</figcaption></figure>

relative displacement of the masses. Applying Lagrange's equations leads to:

$$ \begin{bmatrix} m_3 + m_2 & m_2 \\ m_2 & m_2 \end{bmatrix} \begin{bmatrix} \ddot{u}_3 \\ \ddot{u}_\alpha \end{bmatrix} + \begin{bmatrix} k_3 & 0 \\ 0 & k_2 \end{bmatrix} \begin{bmatrix} u_3 \\ u_\alpha \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}. \quad (18) $$

Both equations of motions (17) and (18) are of the generic form:

$$ M\ddot{q} + Kq = 0, \quad (19) $$

with M and K the inertia and stiffness matrices and where the degrees of freedom's vector are respectively $q = [\gamma \ \alpha]^T$ and $[u_3 \ u_\alpha]^T$ with $\circ^T$ the transpose of vector $\circ$.

The equivalence between the multibody model and the 0D model is obtained if the equations of motion (17) and (18) are identical. The first step is to find a change of variables between the two sets of degrees of freedom $(\gamma, \alpha)$ and $(u_3, u_\alpha)$ and then to identify the members of the inertia and stiffness matrices. This can be made either by comparing Eqs. (17) and (18) or on the associated expressions of the kinetic and potential energies $T = \dot{q}^T M \dot{q}$ and $V = q^T K q$.

The angles $(\gamma, \alpha)$ are associated to equivalent displacements with:

$$ u_3(t) = \lambda_3 \gamma(t), \quad u_\alpha(t) = \lambda_2 \alpha(t), \quad (20) $$

where $(\lambda_3, \lambda_2)$ are two unknown lengths. Introducing Eqs. (20) into (18), multiplying the result by $[\lambda_3 \ \lambda_2]^T$ and identifying the components of the inertia matrix M leads to:

$$ \left(\frac{M_3}{3} + M_2\right) L_3^2 + M_2 \left(\frac{L_2^2 \sin^2(\beta)}{3} - L_2 L_3 \sin(\beta)\right) = (m_3 + m_2)\lambda_3^2, \quad (21) $$

$$ M_2 \frac{L_2^2}{3} = m_2 \lambda_2^2, \quad L_2 M_2 \frac{(-2L_2 \sin(\beta) + 3L_3)}{6} = m_2 \lambda_2 \lambda_3. $$

We then have an under-constrained system with three equations for four unknowns, namely $(M_2, M_3, \lambda_2, \lambda_3)$. By choosing $\lambda_3 = L_3$, one obtains:

$$ \begin{align} \lambda_2 &= \frac{2L_2L_3}{-2L_2 \sin(\beta) + 3L_3}, \quad m_3 = \frac{M_3}{3} + \frac{M_2}{4}, \tag{22} \\ m_2 &= \frac{M_2(-2L_2 \sin(\beta) + 3L_3)^2}{12L_3^2}. \nonumber \end{align} $$

Since $\lambda_2$ is a function of $\beta$, there is a $\beta$ value for which $\lambda_2 > L_2$. In this case, another equivalence can be used. Based on Eq. (21), be choosing $\lambda_2 = L_2$, the result is:

$$ \begin{align} \lambda_3 &= \frac{3L_3}{2} - L_2 \sin \beta, \quad m_3 = \frac{L_3^2 \left(M_2 + \frac{4M_3}{3}\right)}{\left(2L_2 \sin(\beta_2) - 3L_3\right)^2}, \tag{23} \\ m_2 &= \frac{M_2}{3}. \nonumber \end{align} $$

Then, the equivalent stiffnesses $(k_1, k_2)$ result directly from comparison of the stiffness matrices, such that:

$$ k_3 = \frac{\Gamma_3}{\lambda_3^2} \quad \text{and} \quad k_2 = \frac{\Gamma_2}{\lambda_2^2}. \quad (24) $$
<table><caption>**Table 2** Dimensions of the simple case with two articulated beams of Fig. 13</caption><thead><tr><th></th><th></th><th colspan="2">Trailing edge</th><th colspan="2">Compliant link (Γ<sub>2</sub>)</th><th colspan="2">Central beam</th><th colspan="2">Compliant link (Γ<sub>3</sub>)</th></tr></thead><tbody><tr><td>Length</td><td>mm</td><td>L<sub>2</sub></td><td>10</td><td>L<sub>l2</sub></td><td>0.250</td><td>L<sub>3</sub></td><td>8.51</td><td>L<sub>l3</sub></td><td>0.3</td></tr><tr><td>Width</td><td>mm</td><td>b<sub>2</sub></td><td>1.05</td><td>b<sub>l2</sub></td><td>1.05</td><td>b<sub>3</sub></td><td>1.62</td><td>b<sub>l3</sub></td><td>1.62</td></tr><tr><td>Thickness</td><td>μm</td><td>h<sub>2</sub></td><td>185</td><td>h<sub>l2</sub></td><td>25</td><td>h<sub>3</sub></td><td>145</td><td>h<sub>l3</sub></td><td>20</td></tr></tbody></table>

**Fig. 14** Natural frequencies of the two modes of models B of Fig. 13 as a function of angle $\beta$, for the FEA (model B0) and the two multibody/0D equivalent models (models B1<sub>1</sub> and B1<sub>2</sub>)

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

As a conclusion, using Eqs. (20), (22) and (24), the dynamics of the two-beam system, model B1<sub>1</sub> of Fig. 13, can be *exactly replaced* by model B2.

To analyse the accuracy of the equivalence, we now consider a real system (model B0 of Fig. 13) composed of two ribs, made of thick and long beams, linked by compliant links, made of thin and short beams. The values of their lengths, thickness and width are given in Table 2. The comparison is made on the first two eigenmodes (natural frequencies and mode shapes). For the real system (model B0), they are computed with a FEA (with the Abaqus software and quadratic Timoshenko B32 finite-elements) and for the two other models, the two eigenmodes ($\omega_i$, $\Phi_i$), $i = 1, 2$ of Eq. (19) are obtained by solving the eigenvalue problem:

$$ (K - \omega^2 M)\Phi = 0. \qquad (25) $$

We consider two slightly different cases for the multibody model, in the same way than in Sect. 3.1. The first case is model B1<sub>1</sub> sketched in Fig. 13. The second one, called model B1<sub>2</sub>, not sketched for the sake of conciseness, is model B1<sub>1</sub> corrected with the center of the hinges slightly shifted, in the same way as what is done between models A1<sub>1</sub> and A1<sub>2</sub> in Fig. 9. The equivalence with model B2 is still exact but the closed form expressions are more complex and gathered in Appendix B.3.

Figure 14 reports the evolution of the two natural frequencies ($f_i = \omega_i/(2\pi)$, in Hz) as functions of the angle $\beta$, for the real system B0 and the two multibody systems

B1<sub>1</sub> et B1<sub>2</sub>. One must recall that these two latter cases are identical, as explained above, to the 0D model B2. One can observe that, for the first eigenmode, the multibody/0D models are very close to the real system since $f_1$ lies inside a 5% tolerance margin (blue shaded area). For the second mode, a larger discrepancy is observed, such that $\beta$ must be smaller than 50° to guarant an error smaller than 5% for model B1<sub>1</sub> (orange shaded area). This discrepancy is largely corrected by model B1<sub>2</sub>.

The next step is to verify the equivalence on the mode shapes. From the FEA, we easily get the displacement field of the structure projected onto a fixed direction, chosen here as the transverse direction $x_0$. Considering that $\lambda_3 = L_3$, one has:

$$ O_3C_3 \cdot x_0 = \lambda_3 \sin \gamma \approx \lambda_3 \gamma = u_3, \qquad (26) $$

such that $u_3(t)$ is exactly the displacement of point $C_3$ (the tip of the central beam) projected onto $x_0$, as shown in Fig. 13. Then, defining point $A_2$ attached to the rib such that $O_2A_2 = \lambda_2 y_2$ (see Fig. 13), its displacement projected onto $x_0$ is:

$$ \begin{aligned} \bar{u}_2 &= O_3A_2 \cdot x_0, \\ &= \lambda_3 \sin \gamma + \lambda_2(\cos \gamma \sin \alpha - \sin \beta \cos \alpha \sin \gamma) \approx \lambda_3 \gamma \\ &+ \lambda_2(\alpha - \gamma \sin \beta) = u_{\alpha} + u_3(1 - \lambda_2 \sin \beta / \lambda_3), \end{aligned} \qquad (27) \\ \Rightarrow u_{\alpha} \approx \bar{u}_2 - u_3 \left(1 - \sin \beta \frac{\lambda_2}{\lambda_3}\right), \end{aligned} $$
where Eq. (20) have been used. Consequently, one has to correct the displacement $\tilde{u}_2$, stemming from the FEA, to recover the multibody model displacement $u_\alpha$.

The mode shapes $\Phi_i = [u_3 u_\alpha]^T$, $i = 1, 2$, are obtained as following. The mode shapes of model B0 are computed with FEA, as for the natural frequencies, and $(u_3, u_\alpha)$ are obtained by taking the displacement of the first two mode shapes of the structure at points $C_3$ and $A_2$ along $x_0$ and using Eqs. (26), (27). Then, for multibody model $B1_1$ and $B1_2$, the mode shapes are computed with Eq. (25), with $M$ and $K$ defined in Eq. (18) with for $\beta < 25^\circ$ the equivalence of Eq. (22) and for $\beta > 25^\circ$ the one of Eq. (23). Since the mode shapes are known up to a multiplicative constant, only the ratio $r = u_\alpha/u_3$ of the two components counts.

<figure><img src="image_25.png" /><figcaption>**Fig. 15** Mode shape of the two modes of models B in Fig. 13 as a function of angle $\beta$, for the FEA (model B0) and the multibody/0D equivalent model (models $B1_1$). Relative deformation $r = u_\alpha/u_3$ for the</figcaption></figure>

Figure 15a, b shows $r$ as a function of $\beta$, for both models B0 and $B1_1$.

One observes a behaviour analogous to what was observed for the natural frequencies in Fig. 14: the equivalence is excellent for small values of $\beta$, with an error that increases with $\beta$, larger for mode 2 than for mode 1.

Figure 15c–f also shows the mode shapes for two values of $\beta$ ($\beta = 20^\circ$ and $\beta = 70^\circ$). One recovers the excellent equivalence between model B0 and the multibody model, since the deformed shape in the ribs areas is almost undistinguishable from straight lines. This shows that the fully flexible model B0, thanks to its thin compliant links, deforms as if it was a multibody system.

In conclusion, it is possible to dimension a NAV wings with a 0D model since it generates an error less than 5% on

<figure><img src="image_26.png" /><figcaption>first mode (a) and the second mode (b). First mode shape for $c\beta = 20^\circ$ and $e\beta = 70^\circ$; second mode shape for $d\beta = 20^\circ$ and $f\beta = 70^\circ$</figcaption></figure>
the frequencies and mode shapes. The following sections extend these results to the more complicated cases of the NAV architectures shown in Fig. 7.

## 3.3 Equivalence for NAV Architecture 1

In this section, we address NAV architecture 1 (Fig. 7a), where the trailing edge is attached to the leading edge, which is itself fixed to the central beam free end. This architecture is inspired by other NAV of the literature [9, 10]. The equivalent multi-body system, model C1, is shown in Fig. 16.

The central beam (in blue) is identical to the one in Sect. 3.2. The leading edge (in orange) has a mass $M_1$, a length $L_1$ and a center of mass $G_1$ in its middle. It is hinged at point $O_1$ and rotates with respect to the central beam around the axis ($O_1, z_3$) by an angle $\theta(t) = (\hat{y}_1, \hat{y}_3)$, in parallel with a torsional spring of stiffness $\Gamma_1$. The trailing edge (in purple) has a mass $M_2$, a center of mass $G_2$ in its middle and a length $L_2$. It forms a constant angle $\beta = (\hat{y}_\beta, \hat{y}_1)$ along the axis ($O_2, x_1$) (this angle is the same as that in Fig. 2). It is hinged at point $O_2$, assumed merged with $O_1$, and rotates with respect to the leading edge around the axis ($O_2, z_\beta$) by an angle $\alpha(t) = (\hat{y}_2, \hat{y}_\beta)$, in parallel with a torsional spring of stiffness $\Gamma_2$. One has $O_3G_3 = L_3/2 z_3, O_3O_1 = L_3z_3 + c_3y_0, O_1G_1 = L_1/2 y_1$ and $O_2G_2 = L_2/2 y_2$.

Applying Lagrange's equations and linearising them around the equilibrium position ($\gamma = \theta = \alpha = 0$) (see Appendix B.2 for details), one obtains equations of motions of the form Eq. (19) with:

$$ K = \begin{bmatrix} \Gamma_3 & 0 & 0 \\ 0 & \Gamma_1 & 0 \\ 0 & 0 & \Gamma_2 \end{bmatrix} \qquad (28b) $$

and $q = [\gamma, \theta, \alpha]^T$ the degrees of freedom vector.

The equivalent 0D model of this architecture, model C2, is presented in Fig. 16. The central beam, leading edge, and trailing edge are represented by the blue, orange, and purple mass / spring systems, with masses ($m_3, m_1, m_2$) and stiffnesses ($k_3, k_1, k_2$), respectively. Their displacements with respect to the base are represented by ($u_3, u_1, u_2$), respectively. The relative displacements $u_\theta = u_1 - u_3$ and $u_\alpha = u_2 - u_1$ are also defined. The equations of motion are of the form (19) with:

$$ M = \begin{bmatrix} \frac{m_3 + m_1 + m_2}{m_1 + m_2 + M} & \frac{m_1 + m_2 + M}{m_1 + m_2} & \frac{m_2}{m_2} \\ \frac{m_3 + m_1 + m_2}{m_1 + m_2 + M} & \frac{m_1 + m_2 + M}{m_1 + m_2} & \frac{m_2}{m_2} \end{bmatrix}, \quad K = \begin{bmatrix} k_3 & 0 & 0 \\ 0 & k_1 & 0 \\ 0 & 0 & k_2 \end{bmatrix}, \qquad (29) $$

with $[u_3, u_\theta, u_\alpha]^T$ the degrees of freedom vector. The term $M$ is added here for the equivalence property that is addressed in the following and we have $M = 0$ for the equations of motion of model C2.

Using the same approach than in Sect. 3.2, three length parameters ($\lambda_3, \lambda_1, \lambda_2$) are defined to find the equivalence between models C1 and C2, with ($u_3, u_\theta, u_\alpha$) = ($\lambda_3\gamma, \lambda_1\theta, \lambda_2\alpha$). Choosing $\lambda_3 = L_3$, one obtains:

$$ M = \begin{bmatrix} L_3^2 M_1 + \frac{L_3^2 M_3}{3} + M_2 \left( \frac{L_2^2 \sin^2(\beta)}{3} - L_2 L_3 \sin(\beta) + L_3^2 \right) & \frac{L_1 L_3 M_1}{2} + \frac{L_2 M_2 (-L_2 \sin(2\beta) + 3L_3 \cos(\beta))}{\frac{L_1 L_3 M_1}{2} + \frac{L_2 M_2 (-2L_2 \sin(\beta) + 3L_3)}{\frac{L_1 L_3 M_1}{2} + \frac{L_2 M_2 \cos(\beta)}{3}} + \frac{L_2 M_2 (-2L_2 \sin(\beta) + 3L_3)}{\frac{L_1 L_3 M_1}{2} + \frac{L_2 M_2 \cos(\beta)}{3}} - \frac{L_2 M_2 (-2L_2 \sin(\beta) + 3L_3)}{\frac{L_1 L_3 M_1}{2} + \frac{L_2 M_2 \cos(\beta)}{3}}) & \frac{L_2 M_2 (-2L_2 \sin(\beta) + 3L_3)}{\frac{L_1 L_3 M_1}{2} + \frac{L_2 M_2 \cos(\beta)}{3}} \\ \frac{L_1 L_3 M_1}{2} + \frac{L_2 M_2 (-2L_2 \sin(\beta) + 3L_3)}{\frac{L_1 L_3 M_1}{2} + \frac{L_2 M_2 \cos(\beta)}{3}} & \frac{L_1 L_3 M_1}{2} + \frac{L_2 M_2 \cos(\beta)}{\frac{L_1 L_3 M_1}{2} + \frac{L_2 M_2 \cos(\beta)}{3}} & \frac{L_2 M_2 (-2L_2 \sin(\beta) + 3L_3)}{\frac{L_1 L_3 M_1}{2} + \frac{L_2 M_2 \cos(\beta)}{3}} \end{bmatrix} \qquad (28a) $$

<figure><img src="image_27.png" /><figcaption>Fig. 16 a Multibody model (model C1) for NAV architecture 1 and b equivalent 0D model (model C2)</figcaption></figure>
$$
\begin{align}
\lambda_1 &= \frac{2L_2L_3 \cos(\beta)}{-2L_2 \sin(\beta) + 3L_3}, & m_1 &= \frac{L_1^2 M_1 (2L_2 \sin(\beta) - 3L_3)^2}{12L_2^2 L_3^2 \cos^2(\beta)}, \nonumber \\
\lambda_2 &= \frac{2L_2L_3}{-2L_2 \sin(\beta) + 3L_3}, & m_2 &= \frac{M_2 (2L_2 \sin(\beta) - 3L_3)^2}{12L_2^2}, \tag{30} \\
m_3 &= -\frac{L_1^2 M_1 (2L_2 \sin(\beta) - 3L_3)^2}{12L_2^2 L_3^2 \cos^2(\beta)} + M_1 + \frac{M_2}{4} + \frac{M_3}{3}, & & \nonumber
\end{align}
$$

and

$$
k_1 = \frac{\Gamma_1}{\lambda_1^2}, \quad k_2 = \frac{\Gamma_2^2}{\lambda_2}, \quad k_3 = \frac{\Gamma_3}{\lambda_3^2}. \qquad (31)
$$

In the present case, on the contrary to model B1 of Sect. 3.2,
the equivalence is not perfect and we have to add a term *M*
in the mass matrix of Eq. (29), of value:

$$
M = L_1 M_1 \left( \frac{L_1 \tan(\beta)}{3} - \frac{L_1 L_3}{2L_2 \cos(\beta)} + \frac{L_3}{2} \right). \quad (32)
$$

$$
M = \begin{bmatrix}
\left(\frac{M_3}{3} + M_1 + M_2\right) L_3^2 + M_2 \left(\frac{L_2^2 \sin^2(\beta)}{3} - L_2 L_3 \sin(\beta)\right) & M_1 L_3 \frac{L_1}{2} \\
M_1 L_3 \frac{L_1}{2} & M_1 \frac{L_1^2}{3} \\
L_2 M_2 \frac{(-2L_2 \sin(\beta) + 3L_3)}{6} & 0
\end{bmatrix}
$$

with the same stiffness matrix **K** and degrees of freedom
vector **q** than for architecture (1) (Eq. 28b).

Since the trailing edge is now connected to the central
beam, the equivalent 0D model, called model D2 and illus-
trated in Fig. 17, is composed of two mass/spring systems in
parallel, with the relative displacements u<sub>θ</sub> = u<sub>1</sub> - u<sub>3</sub> and
u<sub>α</sub> = u<sub>2</sub> - u<sub>3</sub>. Contrary to architecture 1, the equivalence
between multibody model D1<sub>1</sub> and 0D model D2 is now

This non-perfect equivalence makes this architecture diffi-
cult to optimise. To solve this problem, a second architec-
ture is introduced in the following section.

**3.4 Equivalence for NAV Architecture 2**

Architecture 2 is illustrated in Fig. 7b, for which both the
leading and trailing edges are connected in parallel to the
free end of the central beam. Its multibody model, model
D1₁, is shown in Fig. 17. The main difference with the pre-
vious case (model C1₁ of Fig. 16) is that the trailing edge
geometry and motion are defined with respect to the central
beam. Identically to model B1₁, it forms a constant angle
β = (ỹβ, ỹ0) along the axis (O₂, x₃) and rotates around the
axis (O₂, zβ) by an angle α(t) = (ỹ₂, ỹβ). The correspond-
ing mass matrix is (see Appendix B.2 for details):

$$
\begin{bmatrix}
L_2 M_2 \frac{(-2L_2 \sin(\beta) + 3L_3)}{6} & 0 \\
0 & M_2 \frac{L_2^2}{3}
\end{bmatrix}, \qquad (33)
$$

exact for architecture 2. With the notations of Sect. 3.3, the
stiffness matrix of model D2 is the one of Eq. (29) and the
mass matrix is:

$$
M = \begin{bmatrix} m_3 + m_1 + m_2 & m_1 & m_2 \\ m_1 & m_1 & 0 \\ m_2 & 0 & m_2 \end{bmatrix}, \qquad (34)
$$

<figure><img src="image_28.png" /><figcaption>**Fig. 17** Real system with compliant links (model D0), multibody model (model D1₁) and their equivalent 0D model (model D2), for NAV architecture 2</figcaption></figure>
<table><caption>**Table 3** Dimensions of a prototype of architecture 2, Fig. 17, for an angle of $\beta = 35^{\circ}$</caption><thead><tr><th></th><th>Leading edge</th><th>Compliant link (Γ<sub>1</sub>)</th><th>Trailing edge</th><th>Compliant link (Γ<sub>2</sub>)</th><th>Central beam</th><th>Compliant link (Γ<sub>3</sub>)</th></tr></thead><tbody><tr><td>Length</td><td>mm</td><td>10</td><td>8</td><td>0.249</td><td>8.51</td><td>0.3</td></tr><tr><td>Width</td><td>mm</td><td>1.05</td><td>0.171</td><td>0.171</td><td>1.62</td><td>1.62</td></tr><tr><td>Thickness</td><td>μm</td><td>185</td><td>85</td><td>15</td><td>145</td><td>20</td></tr></tbody></table>

<table><caption>**Table 4** Comparison of natural frequencies for the prototype of architecture 2 described in Table 3, between the real system, model D0, and the two equivalent multibody models: model D1<sub>1</sub> and model D1<sub>2</sub>. The errors are with respect to model D0</caption><thead><tr><th rowspan="2">Mode</th><th colspan="2">Model D0</th><th colspan="2">Model D1<sub>1</sub></th><th colspan="2">Model D1<sub>2</sub></th></tr><tr><th>Frequency (Hz)</th><th>Frequency (Hz)</th><th>Error (%)</th><th>Frequency (Hz)</th><th>Error (%)</th><th>Error (%)</th></tr></thead><tbody><tr><td>1</td><td>38.38</td><td>38.75</td><td>0.96</td><td>38.63</td><td>0.65</td><td>0.65</td></tr><tr><td>2</td><td>82.84</td><td>86.7</td><td>4.66</td><td>82.58</td><td>0.31</td><td>0.31</td></tr><tr><td>3</td><td>139.87</td><td>146.86</td><td>4.95</td><td>135.09</td><td>3.42</td><td>3.42</td></tr></tbody></table>

<table><caption>**Table 5** Comparison of the mode shape amplitudes for the prototype of architecture 2 described in Table 3, between the real system, model D0, and the two equivalent multibody models: model D1<sub>1</sub> and model D1<sub>2</sub></caption><thead><tr><th rowspan="2">Mode</th><th colspan="2">Model D0</th><th colspan="2">Model D1<sub>1</sub></th><th colspan="2">Model D1<sub>2</sub></th></tr><tr><th>Ratio r<sub>1/3</sub></th><th>Ratio r<sub>1/3</sub></th><th>Error (%)</th><th>Ratio r<sub>1/3</sub></th><th>Error (%)</th><th>Error (%)</th></tr></thead><tbody><tr><td>1</td><td>0.313</td><td>0.291</td><td><strong>7.03</strong></td><td>0.306</td><td>2.24</td><td>2.24</td></tr><tr><td>2</td><td>11.78</td><td>16.69</td><td><strong>41.7</strong></td><td>12.12</td><td>2.89</td><td>2.89</td></tr><tr><td>3</td><td>-2.145</td><td>-2.043</td><td>4.75</td><td>-2.102</td><td>2.00</td><td>2.00</td></tr></tbody></table><table><thead><tr><th rowspan="2">Mode</th><th colspan="2">Ratio r<sub>2/3</sub></th><th colspan="2">Ratio r<sub>2/3</sub></th><th colspan="2">Ratio r<sub>2/3</sub></th></tr><tr><th>Ratio r<sub>2/3</sub></th><th>Error (%)</th><th>Ratio r<sub>2/3</sub></th><th>Error (%)</th><th>Ratio r<sub>2/3</sub></th><th>Error (%)</th></tr></thead><tbody><tr><td>1</td><td>0.264</td><td>0.547</td><td><strong>107.2</strong></td><td>0.268</td><td>1.52</td><td>1.52</td></tr><tr><td>2</td><td>-451.54</td><td>-285.48</td><td><strong>36.78</strong></td><td>-462.13</td><td>2.35</td><td>2.35</td></tr><tr><td>3</td><td>-1.529</td><td>-3.3163</td><td><strong>117.4</strong></td><td>-1.521</td><td>0.52</td><td>0.52</td></tr></tbody></table>

Errors larger than 5% are shown in bold

The errors are with respect to model D0

with $\lambda_3 = L_3$ and:

$$ \lambda_1 = \frac{2}{3} L_1, \quad \lambda_2 = \frac{2L_2L_3}{3L_3 - 2L_2 \sin(\beta)}, \qquad (35a) $$

$$ m_1 = \frac{3}{4} M_1, \quad m_2 = \frac{M_2 (2L_2 \sin(\beta) - 3L_3)^2}{12L_3^2}, \qquad (35b) $$

$$ m_3 = \frac{1}{3} M_3 + \frac{1}{4} (M_1 + M_2), $$

$$ k_1 = \frac{\Gamma_1}{\lambda_1^2}, \quad k_2 = \frac{\Gamma_2}{\lambda_2^2}, \quad k_3 = \frac{\Gamma_3}{\lambda_3^2}. \qquad (35c) $$

To verify this equivalence with a real, fully flexible, system, model D0 shown in Fig. 17 is considered with a FEA (Abaqus software with B32 quadratic Timoshenko beam elements). It is identical to model B0 of Fig. 13 with an additional rib for the leading edge. Table 4 presents the eigenfrequencies of architecture 2 obtained for both the real system (model D0) and the two equivalent multibody models D1<sub>1</sub> and D1<sub>2</sub>. Model D1<sub>1</sub> is shown in Fig. 17, whereas model D1<sub>2</sub> is similar, with the center of hinged slightly shifted from $L_{li}/2$ (i=1,2,3), like what was done for model A1<sub>2</sub> of Fig. 9. The comparison are performed using the dimensions of Table 3. In addition, Table 5 presents a comparison of the mode shapes calculated for the same models. The first part of the table presents the ratio between the amplitudes of the
leading edge and the central beam, $r_{1/3} = \theta/\gamma$, whereas the second part of the table shows the ratio between the amplitudes of the trailing edge and the central beam, $r_{2/3} = \alpha/\gamma$.

In analysing these numerical results, one observes an excellent agreement between the models for the natural frequencies, with an error of less than 5% regardless of the torsional spring positions (models D1₁ or D1₂). However, analysis of the mode shapes reveals that employing the D1₂ model is crucial for achieving accurate results (with less than 5% of error). To fully understand these results, one has to remember that the difference between models D1₁ and D1₂ lies in whether or not the length of the compliant joint is taken into account in the equivalent model (Fig. 9). Thus, depending on the D1 model, the λ distance of the point along the rib where the equivalence with the D0 model is made is not the same. This is why the choice was made to compare, in Table 5, the ratio of the amplitudes of the angles rather than the displacements (which was done in Fig. 15). To obtain the value of the angle formed by the ribs for the D1 model, knowing that the calculation performed with Eq. (19) gives the displacement, Eq. (20) is used. It can therefore be noted that the value of the angle is inversely proportional to λ. However, as shown by the expressions of λ of the D1₂ model in Eq. (77), they are highly dependent on the lengths of the compliant links, especially λ₂, which depends on both $L_{l2}$ and $L_{l3}$. It is this dependence that is the origin of the differences observed in Table 5. For the mode shape amplitudes of model D1₁, by neglecting the length of the

compliant links, the equivalent distance (λ), and therefore the angle, is underestimated. It can also be noticed that the less the joint deforms, the less negligible its length becomes, and the greater the errors between model D1₁ and model D0 can be. So, it is preferable to consider the springs in the middle of the compliant links, i.e. model D1₂, such that 0D models could be used to simulate and optimise NAV architecture 2.

## 4 Optimisation of Wing Kinematics from 0D Models

As explained in Sect. 2.1, we are interesting in optimising the structure of our NAVs to naturally obtain a combination of flapping and twisting motions of its wings using a resonant motion. This requires to bring the flapping and twisting modes resonances as close together as possible, while ensuring comparable amplitudes for each of them, as proposed in previous studies [13, 14, 28]. To achieve this, the optimisation is based on the 0D models introduced in Sect. 3, which have the advantages of their reduced number of degrees of freedom. The optimisation criteria are (1) obtaining a periodic motion of the wing, at a particular driving frequency, equivalent to a combination of flapping and twisting motion with a quadrature phase shift and (2) having as much amplitude as possible at this frequency. In our study, since the wing deformation is defined by the positions of its leading and trailing edges, criteria (1) is achieved by looking for a particular motion for which the two edges move periodically with a quadrature phase shift.

Only the dynamic properties of the structure are optimised, to target the particular combined flapping/twisting kinematics. Several points are thus not directly considered in our process. The first one is the produced aerodynamic lift force. It is not computed, and thus not used for the optimisation, since we assume that it is maximal if the particular kinematics with combined flapping and twisting motions is realised, as shown in [36]. The evaluation of the lift force produced by vibrating wings is not an easy task and would require complex and time consuming computational fluid dynamics (CFD) simulations, coupled to a structure motion [42], largely out of the scope of the present article. The second one is the effect of geometrical nonlinearities, linked to the large amplitude motion of the wings, which is targeted and obtained (see Sect. 2.2). The main effect of geometrical nonlinearities is to give a dependence of the resonance frequency upon the amplitude of the motion, a so-called hardening/softening behaviour, that depends on the resonance considered as well as on the boundary conditions of the structure [43, 44]. In the present case of the wings that are hinged at one end and free at the other one, with deformed

<figure><img src="image_29.png" /><figcaption>**Fig. 18** 0D basic model (named E2) of a resonant wing with 2 degrees of freedom with damping, masses ($m_1, m_2$) representing the leading and trailing edges, respectively</figcaption></figure>
shapes close to those of a cantilever beam, the hardening/softening effect can be assumed very weak. For instance, the change of frequency of the first nonlinear bending mode of a cantilever beam is less than 4%, for very large oscillations with a tip rotation of 80° [23, 29, 30, 45]. The geometrical nonlinearities are thus neglected in our design process. Finally, our 0D model also partially neglect the aerodynamic influence of the fluid on the dynamic properties of the structure. Those effect can be viewed as an addition of damping and mass. The damping effect is due to the aeroelastic drag force, mainly proportional to the square of the local wing velocity [36]. This effect is approximated in our model with an equivalent linear damping term (the constants $c_1$, $c_2$ of model E2 of Sect. 4.1), equivalent to the aeroelastic drag at a given motion amplitude [23]. Then, the added mass effect of the fluid mainly affects the value of the resonances frequency, by slightly decreasing their values. In this study, the added mass effect is not taken into account, assuming that it affects both the bending and twisting resonance in the same way, such that our main optimisation parameter, the frequency distance between the two modes, is well predicted by our model.

## 4.1 A Basic Two Degrees of Freedom 0D Model

Let us first consider the simple model E2 of Fig. 18 (and studied in [18]). This architecture corresponds to a reduction of the model C2 (architecture 1) of Fig. 16, where only the leading and trailing edges are considered, with the central beam assumed to be rigid and with its driving functionality replaced by an external force $F(t)$ applied to the leading edge. Let $(m_1, k_1, c_1)$, $(m_2, c_2, k_2)$ be the mass, stiffness damping of the leading and trailing edge, respectively, and $(u_1(t), u_2(t))$ their displacements. Notice that the $c_2$ damper, placed between the $(m_2)$ trailing edge and the ground, models the air damping.

The equations of motion of model E2 are written in a dimensionless form to obtain a reduced number of independent parameters. One defines, $i = 1, 2$:

$$ 
\begin{aligned}
\mu &= \frac{m_2}{m_1}, \quad \omega_i = \sqrt{\frac{k_i}{m_i}}, \quad \xi_i = \frac{c_i}{2m_i\omega_i}, \\
\bar{\omega} &= \frac{\omega_2}{\omega_1}, \quad \bar{t} = \omega_1 t, \quad \bar{\Omega} = \frac{\Omega}{\omega_1},
\end{aligned}
\quad (36) $$

where $\Omega$ is an arbitrary frequency. This leads to the following scaled equations of motion:

$$ 
\begin{aligned}
\ddot{u}_1 + 2\xi_1\dot{u}_1 + (1 + \mu\bar{\omega}^2)u_1 - \mu\bar{\omega}^2 u_2 &= F_0 \cos(\bar{\Omega}\bar{t}), \\
\ddot{u}_2 + 2\xi_2\bar{\omega}\dot{u}_2 + \bar{\omega}^2(u_2 - u_1) &= 0.
\end{aligned}
\quad (37) $$

where $F_0$ is the amplitude of the harmonic forcing $F(t)$. By rewriting Eq. (37) in the frequency domain, one can express the complex amplitude of $(u_1, u_2)$, of the leading and trailing edge motions as:

$$ U_1(\bar{\Omega}) = \frac{\bar{\omega}^2 - \bar{\Omega}^2 + 2j\xi_2\bar{\omega}\bar{\Omega}}{D(\bar{\Omega})} F_0, \quad U_2(\bar{\Omega}) = \frac{\bar{\omega}^2}{D(\bar{\Omega})} F_0, \quad (38) $$

with:

$$ 
\begin{aligned}
D(\bar{\Omega}) = \bar{\Omega}^4 & - 2j\bar{\Omega}^3 (\bar{\omega}\xi_2 + \xi_1) \\
& - \bar{\Omega}^2 (\bar{\omega}^2 (\mu + 1) + 4\bar{\omega}\xi_1\xi_2 + 1) \\
& + 2j\bar{\Omega}\bar{\omega} (\xi_2 (\bar{\omega}^2\mu + 1) + \bar{\omega}\xi_1) + \bar{\omega}^2.
\end{aligned}
\quad (39) $$

We also define the phases of the two motions and their difference:

$$ 
\begin{aligned}
\phi_1 &= \arg U_1, \quad \phi_2 = \arg U_2, \\
\phi &= \phi_2 - \phi_1 = \arg(\bar{\omega}^2) - \arg(\bar{\omega}^2 - \bar{\Omega}^2 + 2j\xi_2\bar{\omega}\bar{\Omega}) \\
&= -\operatorname{atan}\frac{2\xi_2\bar{\omega}\bar{\Omega}}{\bar{\omega}^2 - \bar{\Omega}^2}.
\end{aligned}
\quad (40) $$

From Eqs. (38) and (40), the Frequency Response Function (FRF) for the leading edge $U_1(\bar{\Omega})/F_0$ and trailing edge $U_2(\bar{\Omega})/F_0$ are displayed in Fig. 19 by the solid orange and purple lines, respectively, for representative values of the parameters $\bar{\omega}$, $\mu$, $\xi_1$, $\xi_2$. It can be first observed that the FRF displays two resonances corresponding to a low frequency mode with the two motions in phase (at $\bar{\Omega} = \bar{\omega}_- \simeq 0.97$ with $\phi = 0$) and a high frequency mode with the two motions in phase opposition (at $\bar{\Omega} = \bar{\omega}_+ \simeq 1.08$ with $\phi \simeq -\pi$), identified as a flapping and a twisting mode, respectively. In between, at a particular frequency $\bar{\omega}_q \simeq 1.05$, the two motions $u_1(t)$ and $u_2(t)$ are in quadrature of phase ($\phi = -\pi/2$). Following the closed-form expression of $\phi$ in Eq. (40), the denominator must be zero to obtain $\phi = -\pi/2$, such that this phase quadrature frequency is found exactly as $\bar{\omega}_q = \bar{\omega} \Leftrightarrow \omega_q = \omega_2$. Choosing the driving frequency equal to this quadrature frequency, $\bar{\Omega} = \bar{\omega}_q$, leads to the optimal motion to create a lift force since the leading and the trailing edge are in phase quadrature (see Sect. 2.1). However, observing Fig. 19 shows that this working point appears to be close to an anti-resonance of the leading edge motion, which is an issue in the present case since a large amplitude motion of the wing is required to create a sufficient lift. This issue is addressed in the following, by identifying the optimal value of the independent parameters of Eq. (37) ($\mu$, $\bar{\omega}$, $\xi_1$ and $\xi_2$) to maximise the amplitude $|U_1(\bar{\omega}_q)|$ of the leading edge at the quadrature frequency. One may however notice that this antiresonance property is
<figure><img src="image_30.png" /><figcaption>**Fig. 19** FRF of the leading edge ($U_1(\bar{\Omega})/F_0$, orange solid line) and trailing edge ($U_2(\bar{\Omega})/F_0$, purple solid line), for model E2 in Fig. 18, with $\bar{\omega} = 1.05$, $\mu = 0.01$ and $\xi_1 = \xi_2 = 0.005$. (top) modulus in [dB] and (bottom) phase of ($U_1(\bar{\Omega}), U_2(\bar{\Omega}))/F_0$. The solid black line represents the phase difference $\phi$ between the leading and trailing edges, and the horizontal dotted black line indicates phase quadrature. The black vertical dashed lines indicate flapping and twisting modes natural frequencies ($\omega_-$, $\omega_+$) and the quadrature frequency $\bar{\omega}_2$.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>**Fig. 20** Modulus of the FRF of the leading edge $|U_1(\bar{\Omega})|$ of model E2 (Eq. 38) for several values of the parameters ($\mu, \xi_1, \xi_2, \bar{\omega}$), as specified in the legend. The phase quadrature frequencies are shown with a bullet point.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>Fig. 21 Evolution of the natural frequency distance $\bar{\omega}_+ - \bar{\omega}_-$, for model E2 in Fig. 18, as a function of $\mu$ and $\bar{\omega}$, from Eq. (41)</figcaption></figure>

sometimes an advantage, used for vibration reduction applications in certain tuned vibration dampers [46, 47].

Figure 20 shows the frequency response $|U_1(\bar{\Omega})|$ of the leading edge for several selected values of the system's parameters ($\mu, \xi_1, \xi_2, \bar{\omega}$). One can observe that the anti-resonance point is close to the phase quadrature point, as already noticed, but that the amplitude $|U_1(\bar{\omega}_q)|$ at this point can be increased for a particular choice of the parameters. Indeed, $\bar{\omega}$ impacts the symmetry of the frequency response, that has to be preserved to increase $|U_1(\bar{\omega}_q)|$. Consequently, $\bar{\omega} = 1$ seems a good choice. Moreover, decreasing $\mu$ leads to put closer the two resonances, increasing $|U_1(\bar{\omega}_q)|$ in the meantime since $\bar{\omega}_- < \bar{\omega}_q < \bar{\omega}_+$. Finally, increasing the damping ratio ($\xi_1, \xi_2$) also appears beneficial since it seems to increase $|U_1(\bar{\omega}_q)|$, under a certain limit above which $|U_1(\bar{\omega}_q)|$ starts to decrease, suggesting that there is probably an optimal value of the damping ratios. More quantitative results are proposed in the following.

<figure><img src="image_33.png" /><figcaption>Fig. 22 a Leading edge amplitude at the phase quadrature frequency $|U_1(\bar{\omega}_q)|$, for $\bar{\omega} = 1$ and $\mu = 0.01$, as a function of the damping ratios $\xi_1$ and $\xi_2$ for model E2 in Fig. 18; the black solid line shows the case with $\xi_1 = \xi_2$. b (blue dashed curves, left axis) leading edge amplitude at the phase quadrature frequency $|U_1(\bar{\omega}_q)|$, for $\bar{\omega} = 1$, as a function</figcaption></figure>

of the natural frequencies ($\bar{\omega}_-, \bar{\omega}_+$) of the two modes can be obtained as the zeros of $D(\bar{\Omega})$ in Eq. (38), by neglecting the damping ($\xi_1 = \xi_2 = 0$):

$$ \bar{\omega}_{\pm} = \frac{\sqrt{1 + (1 + \mu)\bar{\omega}^2 \pm \sqrt{(1 + (1 + \mu)\bar{\omega}^2)^2 - 4\bar{\omega}^2}}}{\sqrt{2}} \quad (41) $$

Figure 21 shows the difference $\omega_+ - \omega_-$ as a function of the two parameters ($\bar{\omega}, \mu$).

A minimum is indeed obtained for $\bar{\omega} = 1$, confirming this choice as an optimal one. Then, this plot also confirms that $\omega_+ - \omega_- \text{ decreases with } \mu$. For $\bar{\omega} = 1$ and small values of $\mu$, one can show that $\omega_+ - \omega_- \simeq \sqrt{\mu}$, based on Eq. (41). In conclusion, it is necessary to have $\bar{\omega} = 1$ ($\omega_2 = \omega_1$) and $\mu \ll 1$ ($m_2 \ll m_1$), in order to bring closer the flapping and twisting modes' resonances.

A more quantitative analysis is illustrated in Fig. 22a, that shows the leading edge amplitude at the phase quadrature frequency $|U_1(\bar{\omega}_q)|$, in the case $\bar{\omega} = 1$ and $\mu = 0.01$, as a function of the damping ratios $\xi_1$ and $\xi_2$, computed with Eq. (38). This plot shows that the influences of the two damping ratios on the leading edge amplitude are antagonist: $\xi_2$ must be increased whereas $\xi_1$ must be decreased, with a lesser influence of the latter. However, in practice, it is difficult to control the values of $\xi_1$ and $\xi_2$ separately, since they both mainly stem from the aerodynamic coupling. In addition, they depend nonlinearly on the amplitude of the motion [23, 30]. Consequently, it seems reasonable to consider that $\xi_1 \simeq \xi_2$, illustrated by the black solid line in Fig. 22a, a case considered in details in the following.

Figure 22b shows in dashed blue lines the leading edge amplitude at the phase quadrature frequency $|U_1(\bar{\omega}_q)|$ as a function of equal damping ratios $\xi_1 = \xi_2$, for several values of $\mu$. It is observed that all the curves present a maxima,

of $\xi_1 = \xi_2$, for several values of $\mu$, as indicated; (black solid curve, left axis) locus of the maxima of the dashed curves, i.e. maximum leading edge amplitude as a function of $\xi_1 = \xi_2$; (orange solid curve, right axis) value of $\mu_{opt}$ of $\mu$ at the maxima of the dashed curves, as a function of $\xi_1 = \xi_2$
identified as an optimal point. As a consequence, one can plot the curve that crosses all the maxima when $\mu$ is varied, which can be viewed as a three-dimensional (3D) curve in the space ($\xi_1 = \xi_2, \mu, |U_1(\bar{\omega}_q)|$). One obtains the solid black and orange curves of Fig. 22b, which are indeed the projections of this 3D curve in the planes ($\xi_1 = \xi_2, |U_1(\bar{\omega}_q)|$) and ($\xi_1 = \xi_2, \mu$). These curves show that for a given damping ratio $\xi_1 = \xi_2$, there exist an optimal value $\mu_{opt}$ (defined by the orange curve) of $\mu$ for which the leading edge amplitude is maximal. Recalling that $|U_1(\bar{\omega}_q)|$ is dimensionless (see Eq. 38) and that it gives the amplitude at the phase quadrature point relative to its static value (the value it has for $\Omega = 0$), we are interested in values of $|U_1(\bar{\omega}_q)| > 1$ (larger than one) to obtain an amplification effect. The black curve of Fig. 22b shows that provided the right $\mu_{opt}$ is guaranteed by the design, large amplification factors can be theoretically obtained in practice ($|U_1(\bar{\omega}_q)| = 9$ for $\xi_1 = \xi_2 = 0.025$ and $\mu = 0.03$ or $|U_1(\bar{\omega}_q)| = 5$ for $\xi_1 = \xi_2 = 0.05$ and $\mu = 0.01$).

## 4.2 0D Model of NAV Architecture 1

In this section, the aim is to optimise the model C2 of Fig. 16. It differs from the system in the previous section Sect. 4.1 by the addition of a third degree of freedom (mass 3) between the ground and the two—leading and trailing—edges (masses 1 and 2), that takes into account the central beam of our NAV architecture 2, that constitutes the motor.

In the same way as in the previous section Sect. 4.1, we introduce the following dimensionless parameters, to simplify the analysis of this system (i.e. reduce the number of parameters):

$$ \mu_2 = \frac{m_2}{m_1}, \quad \bar{\omega}_2 = \frac{\omega_2}{\omega_1}, \quad \mu_3 = \frac{m_3}{m_1}, \quad \bar{\omega}_3 = \frac{\omega_3}{\omega_1}. \qquad (42) $$

Notice that parameters ($\mu_2, \bar{\omega}_2$) are equivalent to parameters ($\mu, \bar{\omega}$) of Sect. 4.1.

Although it was shown in Sect. 3.3 that the equivalence between the 0D model (model C2) and the multibody model (model C1) is not perfect for this NAV architecture, we will here study the actual matrix of model C2, corresponding to $M = 0$ in Eq. (29). The accuracy of this approach is limited, but it has the benefits of simplifying the analysis and allows to understand the impact of the independent parameters of Eq. (42) on this system. Moreover, following the results of Sect. 4.1, a key point is to design the structure such that its flapping and twisting resonances are as close as possible. We therefore focus on a simple modal analysis, disregarding the effect of damping, that can be considered in a second step to increase the value of the motion amplitude at the phase quadrature frequency, as done in Sect. 4.1. The

dimensionless matrices of inertia and stiffness are obtained using Eq. (42) in Eq. (29), as:

$$ M = \begin{bmatrix} 1 + \mu_2 + \mu_3 & 1 + \mu_2 & \mu_2 \\ 1 + \mu_2 & 1 + \mu_2 & \mu_2 \\ \mu_2 & \mu_2 & \mu_2 \end{bmatrix}, \quad K = \begin{bmatrix} \mu_3 \bar{\omega}_3^2 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & \mu_2 \bar{\omega}_2^2 \end{bmatrix}. \qquad (43) $$

The three natural modes of this system are obtained by solving the eigenvalue problem (Eq. 25). Based on the conclusions of Sect. 4.1 about model E2, we first set $\bar{\omega}_2 = \bar{\omega} = 1$ and $\mu_2 = \mu = 0.01$. Typical mode shapes are represented in Fig. 23a, b for $\mu_3 = 0.01$ and $\bar{\omega}_3$ equal to 3.5 and 0.5 respectively. The triangle represents the masses that have been moved and the direction of the motion, the rectangle the initial positions and the lines the springs' motion. The colour code is the same as previously, the blue mass represents $m_3$, the orange $m_1$ and the purple $m_2$. Based on the relative motions of the three degrees of freedom, the three modes will be called the flapping mode, the twisting mode and the central beam mode. The flapping mode corresponds to the case where the trailing and leading edges are in phase, whereas they are in phase opposition for the twisting mode. The third mode is called the central beam mode since it involves a large motion of $m_3$, associated to the central beam, and also because it depends mainly on its parameters.

In Fig. 23c, the three natural frequencies are plotted as a function of $\bar{\omega}_3$ for several values of $\mu_3$. The zoom allows to illustrate the mode veering, i.e. the zone where the natural frequency of two modes veer apart instead of crossing when varying a parameter, here $\bar{\omega}_3$. This veering also changes the mode shape, as shown by the position of the central beam mode in Fig. 23a, b. Moreover, the variation of $\mu_3$ is represented using the variation of the line style: the solid line is for $\mu_3 = 0.01$, the dashed line for $\mu_3 = 0.1$ and the dotted line for $\mu_3 = 1$.

The effect of the central beams is clear in Fig. 23c, as it modifies the natural frequencies values around the (black thin) line $\bar{\omega}_3 = \bar{\Omega}$, giving the skeleton of the plot. In order to recover the behaviour of model E2, with close flapping and twisting resonances, the central beam mode natural frequency must be rejected in the high frequencies by choosing $\bar{\omega}_3 > 1$. Moreover, contrary to $\mu_2, \mu_3$ must be chosen large ($\mu_3 \approx 1$) to bring closer the twisting and flapping modes resonances. This means that the central beam must be far more rigid than the leading and trailing edges.

## 4.3 0D Model of NAV Architecture 2

Let us now study the NAV architecture 2 (Fig. 7b), where the trailing edge is attached to the central beam. In this case, there is a complete equivalence between the multibody model (model D1) and the 0D model (model D2). Applying the same methodology as Sect. 4.2, by introducing the
<figure><img src="image_34.png" /><figcaption>Fig. 23 Natural modes of model C2 for $\bar{\omega}_2 = 1$ and $\mu_2 = 0.01$, a, b View of the three mode shapes for $a\bar{\omega}_3 = 3.5$, $\mu_3 = 0.01$ and $b\bar{\omega}_3 = 0.5$ and $\mu_3 = 0.01$, where the square/triangular markers denote the reference/deformed position of the masses (the direction of the triangular marker also shows the direction of deformation) and the orange/pur-</figcaption></figure>

ple/blue colors denote masses $m_1, m_2, m_3$, respectively; c evolution of the three natural frequencies (x-axis), shown in blue, orange and yellow, as a function of $\bar{\omega}_3$ (y-axis), for several values of $\mu_3$ (see the legend)

dimensionless parameters of Eq. (42) into Eq. (34), we obtain the following inertia and stiffness matrices for this geometry,

$$M = \begin{bmatrix} 1 + \mu_2 + \mu_3 & 1 & \mu_2 \\ 1 & 1 & 0 \\ \mu_2 & 0 & \mu_2 \end{bmatrix}, \quad K = \begin{bmatrix} \mu_3 \bar{\omega}_3^2 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & \mu_2 \bar{\omega}_2^2 \end{bmatrix}. \quad (44)$$

The three mode shapes of this system are obtained by solving the eigenvalue problem (Eq. 25) and are illustrated in Fig. 24 for $\bar{\omega}_3 = 3$ and $\mu_3 = 1.3$ and several values of $(\bar{\omega}_2, \mu_2)$. One can first observe that the case for which $\bar{\omega}_2 = 1$ and $\mu_2 = 1$ (Fig. 24a), representative of a symmetric

system, must be avoided because the twisting mode is symmetric, making it unable to be excited by the central beam (blue mass).

Based on the conclusions of Sect. 4.2, it is necessary to have $\bar{\omega}_3 > 1$ to reject the central beam mode in the high frequencies and $\mu_3 \approx 1$ to bring the twisting and flapping modes closer together. Therefore, we set $\bar{\omega}_3 = 3$ and $\mu_3 = 1.3$, and determine the optimum value of $\bar{\omega}_2$ to have the maximum possible amplitude at the phase quadrature frequency. Because in this case the trailing edge is fixed on the central beam and no longer on the leading edge, the $\bar{\omega}_2$ optimisation of Sect. 4.1 is no more applicable.
<figure><img src="image_35.png" /><figcaption>**Fig. 24** Natural modes of model D2 for $\bar{\omega}_3 = 3$ and $\mu_3 = 1.3$. **a–d** view of the three mode shapes for **a** $\bar{\omega}_2 = 1$, $\mu_2 = 1$, **b** $\bar{\omega}_2 = 0.96$, $\mu_2 = 0.01$, **c** $\bar{\omega}_2 = 3.5$, $\mu_2 = 0.01$, **d** $\bar{\omega}_2 = 0.5$, $\mu_2 = 0.01$, where the square/triangular markers denote the reference/deformed position of the masses (the direction of the triangular marker also shows the direction of deformation) and the orange/purple/blue colors denote masses $m_1, m_2, m_3$, respectively; **e** evolution of the three natural frequencies (x-axis), shown in blue, orange and yellow, as a function of $\bar{\omega}_2$ (y-axis), for several values of $\mu_2$ (see the legend)</figcaption></figure>

Figure 24e illustrates that, as in the case of Sect. 4.1, $\mu_2 \ll 1$ must be verified to bring the flapping and twisting modes resonances close. Moreover, one can note that the skeleton of the plot is given by three lines: (i) the oblique line $\bar{\omega}_2 = 1$, corresponding to the natural frequency of oscillator ($m_2, k_2$) grounded (with $u_3 = 0$), shown with a thin black line, and (ii) two vertical lines (shown dash-dotted) corresponding to the two natural frequencies $\bar{\omega}_{1,3}^0$ of the system ($m_1, k_1, m_3, k_3$) uncoupled from ($m_2, k_2$) (i.e. with ($k_2, m_2$) = (0, 0) and with the last rows and column of ($M, K$) canceled in Eq. (44)), that take the values:

Results are illustrated in Fig. 24, which is built on the same principle as Fig. 23 but in (e) the modes frequencies are a function of $\mu_2$ and $\bar{\omega}_2$. In Fig. 24b–d, the three mode shapes are observed for $\bar{\omega}_2 \in \{0.96, 3.5, 0.5\}$. Because it is $\bar{\omega}_2$ that varies in this case, one can notice that the order of appearance in the spectrum of the modes (flapping, twisting and central beam, as defined earlier) does not change after the mode veering. Instead, it is the amplitude of the trailing edge ($m_2$, purple mass) that changes, as illustrated in Fig. 24b–d, where $m_2$ has the greatest amplitude for the central beam mode when $\bar{\omega}_2 = 3.5$ and the flapping mode when $\bar{\omega}_2 = 0.5$.
<figure><img src="image_36.png" /><figcaption>Fig. 25 a Optimal value of $\bar{\omega}_2 = \bar{\omega}_1^0$, as defined by Eq. (45), as a function of $(\bar{\omega}_3, \mu_3)$; b difference $\Delta = \bar{\omega}_+ - \bar{\omega}_-$ between the flapping and twisting natural frequencies for model D2, in the optimal case $\bar{\omega}_2 = \bar{\omega}_1^0$, as a function of $(\bar{\omega}_3, \mu_3)$. Those two surfaces are plotted for $\mu_2 = 0.03$</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>Fig. 26 Frequency difference between the twisting and flapping modes for an optimal value $\bar{\omega}_2$ with $\bar{\omega}_3 = 3$ and $\mu_3 = 1.3$ and for the two architectures (models C2 and D2), as specified in the legend</figcaption></figure>

$$ \bar{\omega}_{1,3}^0 = \sqrt{\frac{\bar{\omega}_3^2 + 1}{2} + \frac{1}{2\mu_3} \pm \frac{\sqrt{\mu_3^2 (\bar{\omega}_3^2 - 1)^2 + 2\mu_3 (\bar{\omega}_3^2 + 1) + 1}}{2\mu_3}} \quad (45) $$

One can notice that there is an optimal value of $\bar{\omega}_2$, detailed in the inset of Fig. 24e, that brings the flapping and twisting natural frequencies the closest to each other, which corresponds to the intersection of the oblique line and the first vertical one, such that $\bar{\omega}_2 = \bar{\omega}_1^0$, defined in Eq. (45). This optimal value of $\bar{\omega}_2$ is thus a function of $(\bar{\omega}_3, \mu_3)$ only, shown as a surface in Fig. 25a.

The difference between the natural frequency of the flapping and twisting mode, $\Delta = \bar{\omega}_+ - \bar{\omega}_-$, is shown in Fig. 25b as a function of $(\bar{\omega}_3, \mu_3)$ in the optimal case $\bar{\omega}_2 = \bar{\omega}_1^0$ and for $\mu_2 = 0.03$. It shows that the interesting regions to obtain the two natural frequencies close to each other, i.e. the low values of $\Delta$, are the yellow areas in the

left and upper right parts of the $(\bar{\omega}_3, \mu_3)$ space. Those two areas correspond to values of $\bar{\omega}_2$, readable in Fig. 25a, close to $\bar{\omega}_2 \simeq 0$ and $\bar{\omega}_2 \simeq 1$, respectively. For other values of $\mu_2$, the shape of the surfaces is the same. Between those two areas of low $\Delta$, the left one, obtained for small values of $\bar{\omega}_3$, is not acceptable in practice since it brings the central beam mode natural frequency below the flapping and twisting ones, as illustrated in Fig. 23b, c (these plots are in fact for model C2, but the effect is the same for model D2). Consequently, the optimal parameter area is the yellow part on the right of Fig. 25b, thus for high values of $\mu_3$ and $\bar{\omega}_3$. To illustrate this point, the level set for $\bar{\omega}_2 = 0.96$ is shown as a black curve on the surface of Fig. 25a, with the blue curve showing the corresponding line in the $(\bar{\omega}_3, \mu_3)$ plane. The corresponding curve on the surface of Fig. 25b, obtained for the same values of $(\bar{\omega}_3, \mu_3)$, is also shown. It is associated to very low values of $\Delta \simeq 0.002$.

In conclusion, to bring the flapping and twisting resonances as close together as possible and maximise the vibration amplitude at the phase quadrature frequency, one must have $\mu_2 \ll 1$. Then, since the optimal $\bar{\omega}_2 = \bar{\omega}_1^0$ defined by Eq. (45) depends on the two parameters $(\bar{\omega}_3, \mu_3)$, one needs to select a point lying in the optimal yellow area on the right part of the surface of Fig. 25b, associated with high values of $(\bar{\omega}_3, \mu_3)$. In practice, $\bar{\omega}_3$ must be chosen not too high to guarantee a not too high stiffness $k_3$ of the central beam, such that it remains flexible enough to effectively excite the NAV as an actuator, as explained in paragraph Sect. 2.1.

## 4.4 Comparison Between the Two Architectures 1 and 2

As a reminder, obtaining an optimal configuration means bringing the flapping and twisting modes close enough together. To compare the two architectures (see Fig. 7, models
C2 and D2), we select $\bar{\omega}_3 = 3$, $\mu_3 = 1.3$ in the optimal area of Fig. 25b, and take the optimal value of $\bar{\omega}_2$. For architecture 1, the optimal value is $\bar{\omega}_2 = 1$, as illustrated in Fig. 21, and for architecture 2, it is defined by Eq. (45) ($\bar{\omega}_2 = 0.96$ in this case). Figure 26 shows the natural frequency difference between the twisting and flapping modes ($\Delta = \bar{\omega}_+ - \bar{\omega}_-$) as a function of the mass ratio $\mu_2 = m_2/m_1$. One can notice that architecture 2, where both edges are attached to the central beam (model D2), is more efficient than architecture 1 (with the masses in series, model C2), because for a given value of $\mu_2$, the frequency difference $\Delta$ is much smaller.

Another remark is that a numerical analysis on the mode shapes shows that for $\bar{\omega}_2$ optimal and $\mu_2$ small, the amplitude of the leading edge is always very small as compared to the one of the trailing edge. This is illustrated for architecture 1 in Fig. 19 (notice that the $|U_1|$ curve (leading edge) is about 20 dB smaller, at the two resonances, than the $|U_2|$ curve (trailing edge)) and in Fig. 23b, and for architecture 2 in Fig. 24b (the orange arrows are much smaller than the purple arrows). This is a direct consequence of the small value of $\mu_2$. However, this effect is lessened in the case of a real NAV, thanks to the non-linear aeroelastic damping that increases with the motion amplitude [23].

## 5 Dimensioned Geometry

In the previous Sect. 4.3, a detailed parameter analysis was proposed to address the conditions to bring the flapping and twisting modes resonances and maximising at the same time the motion amplitude at the phase quadrature frequency for architecture 2, which was found superior to architecture 1. Since all the study was performed on a reduced set of independent dimensionless parameters, based on those guidelines, it is now necessary to define the real geometry of an optimal wing, thus giving values to all dimensional material and geometrical parameters. The case of a two wings NAV has also to be discussed as well. This section thus addresses the complete sizing of a NAV, followed by a comparison of FRFs obtained by model D2 and by a finite element model of a NAV with a real geometry. A detail is that the central

beam will be considered here as fully compliant, as it is the practical manufacturing solution chosen for our NAVs. It also reduces the number of parameters to be set for the sizing of this central beam, since only a single set of width/ thickness/length parameters is required since no compliant link is used.

### 5.1 A Single Compliant Central Beam

This section addresses the specific case of the central beam, which is fully compliant in the real architectures of Fig. 7, in order to model its vibrations by an equivalent mass/spring system, as shown in Fig. 27. Considering it as a cantilever beam of length $L$ and rectangular cross section, with (h, b) its thickness and width, its motion is governed by Eq. (14). Its vibrations onto the first eigenmode can be written:

$$ 
\begin{aligned}
w(x,t) &= \Phi(x)q(t), \quad \ddot{q} + \omega_1^2 q = Q(t), \\
Q(t) &= \frac{1}{\rho S} \int_0^L p(x,t)\Phi(x) \,dx / \int_0^L \Phi^2(x) \,dx,
\end{aligned} 
\quad (46) 
$$

with $\Phi(x)$ its first mode shape, defined by Eq. (5a). The above equations are obtained using the orthogonality property of $\Phi(x)$ [40]. With the normalisation, one has $\Phi(x = L) = 1$ such that the modal coordinate is directly the tip displacement of the cantilever, $q(t) = w_L(t)$. Moreover, applying a point force $F_L(t)$ at the tip of the cantilever leads to $p(x,t) = F_L(t)\delta(x - L)$, with $\delta$ the Dirac delta function, such that $Q(t) = 4F_L(t)/(\rho SL)$ since $\int_0^L \Phi^2(x) \,dx = L/4$. Since the first natural frequency of the cantilever is $\omega_1 = \beta_1^2/L^2\sqrt{EI/(\rho S)}$, the dynamics of this system reads:

$$ M\ddot{w}_L + Kw_L = F_L, \quad (47) $$

with

$$ M = \frac{\rho SL}{4}, \quad K = \frac{\rho SL\omega_1^2}{4} = \frac{\beta_1^4 E bh^3}{48L^3}, \quad (48) $$

<figure><img src="image_38.png" /><figcaption>Fig. 27 Equivalence between the central beam of the NAV and a mass/spring system</figcaption></figure>
and $\beta_1 = 1.8751$ as in Sect. 2.3. Equation (47) proves the equivalence with the motion of a mass/spring system illustrated in Fig. 27.

## 5.2 Dimension Analysis

The aim of this section is to resize an optimised system of architecture 2, Fig. 17, in order to validate the optimisation of Sect. 4.3 while also achieving dimensions of a NAV that can be manufactured. For this purpose, a NAV with 2 wings is studied. If we consider the NAV prototype to be perfectly symmetrical, studying 2 wings is equivalent to analysing model D2 with 4 mass/spring systems attached in parallel to the blue mass $m_3$ in Fig. 17, instead of 2. Two of these systems are identical and form the leading edges, while the other two, also identical, form the trailing edges. If the symmetry is perfect, the full system keeps the same flapping an twisting modes than model D2, with identical displacements on the two wings, but with doubled values for the leading and trailing edge stiffnesses ($k_1, k_2$) and masses ($m_1, m_2$). Considering two identical wings also adds in the spectrum purely asymmetric modes, for which the displacements on the two wings are exactly opposite, with a zero displacement of $m_3$. Because of this, those asymmetric modes are in practice not excited (or only slightly) by the motion of the central beam $m_3$ and are neglected in the present study. Consequently, according to Eq. (34), the following matrices are considered:

$$ M = \begin{bmatrix} m_3 + 2m_1 + 2m_2 & 2m_1 & 2m_2 \\ 2m_1 & 2m_1 & 0 \\ 2m_2 & 0 & 2m_2 \end{bmatrix}, \quad K = \begin{bmatrix} k_3 & 0 & 0 \\ 0 & 2k_1 & 0 \\ 0 & 0 & 2k_2 \end{bmatrix}, \quad (49) $$

with $m_1, k_1, m_2, k_2$ and $k_3$ the equivalent masses and stiffnesses of model D2 (Fig. 17b) defined by Eq. (35). Only the value of $m_3$ will change, as it depends on the sum of the masses of the system (M[1, 1] in Eq. (49) and Eq. (34)) and using the same approach than in Sect. 3.2 with the mass matrices of Eqs. (33) and (49), $m_3 = 1/3M_3 + 1/2(M_1 + M_2)$.

The sizing consists in choosing all the geometrical and material parameters defining the NAV. From Sect. 2.1, all the structural parts of the NAV are built in the same SU-8 material. The geometry is the one of model D0 of Fig. 17, defined by the angle $\beta$, 6 geometry parameters for the leading and the trailing edges and only 3 for the central beam

since it is made fully compliant (see Sect. 5.1), giving 16 geometrical parameters to be chosen. The optimisation study of Sect. 4.3 leads to the choice of the 4 dimensionless parameters ($\mu_2, \bar{\omega}_2, \mu_3, \bar{\omega}_3$), such that the problem is under constrained.

Thus, a number of parameters must be chosen, as explained hereafter and illustrated on an example.

We first choose $\beta = 35^\circ$ such the NAV wings share similarities to real insects wings. The SU-8 properties are $E = 4.88\,GPa$ and $\rho = 1290\,kg.m^{-3}$. Then, equal width for the compliant links and the ribs are chosen, to avoid the appearance of fillets at the extremity of the compliant links during manufacturing, leading to:

$$ b_{l1} = b_1, \quad b_{l2} = b_2. \quad (50) $$

To simplify the process by minimising the number of manufacturing layers, the thicknesses of the two compliant links $\Gamma_1$ and $\Gamma_2$ are made equal and the thicknesses of the leading edge and that of the central beam are also matched:

$$ h_{l2} = h_{l1}, \quad h_3 = h_1. \quad (51) $$

Finally, the two (flapping and twisting) modes natural frequencies must be around a chosen characteristic frequency $f_{char} = 100\,Hz$. Based on Eq. (13) with $\omega_1 = 2\pi f_{char}$, the length $L_{l1}$ can be determined as follows:

$$ L_{l1} = \frac{Eb_{l1}h_{l1}^3}{16L_1^3 f_{char}^2 \pi^2 \rho b_1 h_1}. \quad (52) $$

These five constraints, in addition to the choice of the four dimensionless optimised parameters ($\mu_2, \bar{\omega}_2, \mu_3, \bar{\omega}_3$), reduce the number of free parameter to $16 - 5 - 4 = 7$. The seven underlined parameters in Table 6 ($L_1, b_1, h_1, h_{l1}, L_2, h_2, L_3$) are chosen as leading ones, which condition the values of the other parameters ($b_{l1}, b_{l2}, h_{l2}, h_3, L_{l1}$) (see above equations) and ($b_2, L_{l2}, b_3$) thanks to Eqs. (79, 80, 81) in Appendix B.4. The process is illustrated as follows.

Based on the optimisation discussed in Sect. 4.3, $\mu_2$ is chosen small to ensure the flapping and twisting resonances closer together. A value of $\mu_2 = 0.03$ is selected. For the same reason, the optimal area of Fig. 25b is considered to choose ($\bar{\omega}_3, \mu_3$), which imposes the value of $\bar{\omega}_2$ with Eq. (45).

<table><caption>**Table 6** Dimensions of an optimised prototype of architecture 2 (model D0 of Fig. 17) with two wings, for $\beta = 35^\circ$ and having as dimensionless parameters $\bar{\omega}_3 = 3$, $\mu_3 = 1.3$, $\mu_2 = 0.03$ and $\bar{\omega}_2 = 0.96$ and the material parameters of SU-8 $E = 4.88\,GPa$ and $\rho = 1290\,kg.m^{-3}$</caption><thead><tr><th></th><th></th><th>Leading edge</th><th>Compliant link (Γ<sub>1</sub>)</th><th>Trailing edge</th><th>Compliant link (Γ<sub>2</sub>)</th><th>Central beam</th></tr></thead><tbody><tr><td>Length</td><td>mm</td><td><u>L<sub>1</sub></u></td><td><u>L<sub>l1</sub></u> Eq. (52)</td><td><u>L<sub>2</sub></u></td><td><u>L<sub>l2</sub></u> Eq. (80)</td><td><u>L<sub>3</sub></u></td></tr><tr><td>Width</td><td>mm</td><td><u>b<sub>1</sub></u></td><td><u>b<sub>l1</sub></u> = b<sub>1</sub></td><td><u>b<sub>2</sub></u>(Eq.79)</td><td><u>b<sub>l2</sub></u> = b<sub>2</sub></td><td><u>b<sub>3</sub></u>(Eq.81)</td></tr><tr><td>Thickness</td><td>µm</td><td><u>h<sub>1</sub></u></td><td><u>h<sub>l1</sub></u></td><td><u>h<sub>2</sub></u></td><td><u>h<sub>l2</sub></u> = h<sub>l1</sub></td><td><u>h<sub>3</sub></u> = h<sub>1</sub></td></tr></tbody></table>

The underlined parameters are the primary chosen parameters
The two values $\bar{\omega}_3 = 3$, $\mu_3 = 1.3$ are selected, leading to $\bar{\omega}_2 = 0.96$ from Eq. (45) and $\Delta = (\omega_+ - \omega_-)/\omega_1 = 0.013$. The central beam must remains flexible enough to effectively excite the NAV as an actuator, leading to a not too large $\bar{\omega}_3 = \omega_3/\omega_1$.

Then, the dimensions of the leading edge and its compliant link $\Gamma_1$ are first chosen, as the system has been made dimensionless relative to their properties (Eq. 42). $L_1 = 10$ mm is selected as a characteristic length of our NAV, $h_1$ is chosen as an integer multiple of 10 µm because of manufacturing constraints and $b_1$ (as well as $h_1$) are chosen to have a $b_2$ not too small ($> 0.1$ mm). $h_{l1}$ is chosen also as an integer multiple of 10 µm to have a not two small $L_{l1}$ ($> 0.1$ mm). Then, $L_2$ is chosen close to $L_1$, $h_2$ an integer multiple of 10 µm, both chosen to control the value of $h_{l2} \in [0.1, 0.5]$ mm. Finally, $L_3$ is selected so as to obtain $b_3 > 1.5$ mm as this represents the diameter of the magnet to be attached to the central beam. $L_3$ is also manually adjusted to obtain a $h_3$ value that verifies both $h_3 = h_1$ (Eq. (51), coming from manufacturing constraints) and Eq. (82), the latter being mandatory to have the right values of the dimensionless optimisation parameters. One has also to verify that the compliant links play their role with the ribs remaining rigid under deformation, by ensuring $h_i > 3h_{li}$ and $L_i > 15L_{li}$, $i = 1, 2$, as specified at the end of Sect. 3.1.

One has to note that the above sizing is based on Eq. (35), obtained from the simplified model $D1_1$, that does not include the correction from the length of the compliant link (compare models $A1_1$ and $A1_2$ in Fig. 9). To include it

and obtain a finely adjusted prototype, we use Eq. (77) (in Appendix B.3) of model $D1_2$ to compute $\mu_3$ and $\bar{\omega}_3$ to obtain $\bar{\omega}_2$ with Eq. (45) (in this case $\bar{\omega}_2 = 0.92$), which is different from the targeted optimal value of model $D1_1$ ($\bar{\omega}_2 = 0.96$). This is corrected in a second step by adjusting manually the value of $L_{l2}$ (which conditions the stiffness $\Gamma_2$).

Table 6 gathers all the numerical values of an optimised prototype of architecture 2, model D0 of Fig. 17. In practice, all the constraints and equations are entered in a computer spreadsheet and the optimisation is finely adjusted by trial and error.

Figure 28 shows the frequency response functions (FRFs) computed using the 0D equivalence model (model D2 of Fig. 17, solid line) and a FEA for the optimised prototype of architecture 2 (model D0, dashed line) with the parameters of Table 6. Precisely, the response of model D2 is the solution in the frequency domain of $M\ddot{q} + Kq = F$ with $\dot{q} = [u_3 \ u_1 \ u_2]^T$ and $F = [\cos(\Omega t) \ 0 \ 0]^T$ with $M$ and $K$ of Eq. (49). For the FEA simulation, $(u_1, u_2, u_3)$ are recovered at the point in the geometry defined by $(\lambda_1, \lambda_2, \lambda_3)$ of Eq. (35). A structural damping model is used in both simulations, by using a complex stiffness matrix in the frequency domain (i.e. replacing $K$ by $K(1 + j\eta)$ with $\eta = 0.01$). An excellent match is obtained on the resonance frequencies, thus validating our optimisation process. Some slight discrepancies are noticed on the amplitudes, that may be attributed to the differences between the two models. In practice, slight changes on the value of $L_{l2}$ enable to finely adjust the geometry of the prototype, directly checked on the FEA.

<figure><img src="image_39.png" /><figcaption>**Fig. 28** a FRF of architecture 2 with the geometry parameters of Table 6; solid line correspond to model D2 with 2 wings and dashed line to the FEA model. $(u_1, u_2, u_3)$ (leading edge, trailing edge, cen-</figcaption></figure>

tral beam) are shown in orange, purple and blue colors, respectively. (b,c) FEA model deformed shapes at the flapping and twisting resonance frequencies
<figure><img src="image_40.png" /><figcaption>Fig. 29 Compliant link and rib</figcaption></figure>

## 6 Conclusion

The goal of this work is to present a new design and micro-fabrication of a bioinspired NAV. As compared to previous fully flexible designs, we localise the compliance in small elastic regions connected to rigid bodies, which enables modeling of the complete structure as a three degree of freedom (0D) mass/spring system, provided that geometric criteria for the compliant links are respected. This simplifies the optimisation process, allowing us to optimise a NAV with only three parameters fixed by fabrication criteria and the actuator. This new design also allows for very high vibration amplitude without damaging the structure. The resulting prototypes, manufactured using microfabrication technologies, are capable of producing large-amplitude wing motions up to a 37° angle (Fig.5), while also capable of being driven in their flapping and twisting modes independently or in phase quadrature. Finally, for a NAV with two wings, an excellent correlation with FEA is achieved, with less than 5% error on resonance frequencies. These results validate the new design strategy for a NAV outlined in this work.

The methodology presented in this article is the basis of the design of a new generation of NAVs with two independently actuated wings, with the first results recently presented in [48, 49].

## Appendix A. Simple Model of a Compliant Link

This section aims at showing that the dynamics of a long rib of length L connected to a small compliant link of length Ll (such that $Ll \ll L$) can be approximated with a compliant link subjected to a bending moment only, Fig. 29.

The deformations of the compliant link are considered in statics. If ($F_L, M_L$) denote the shear force and the bending moment applied by the rib on the compliant link at its tip B, simple strength of material calculus leads to the following expressions for the transverse displacement $w_L$ of B and the rotation $\theta_L$:

$$ \theta_L = \frac{Ll}{EI} \left( M_L + \frac{Ll}{2} F_L \right), \quad w_L = \frac{Ll}{EI} \left( \frac{Ll}{2} M_L + \frac{Ll^2}{3} F_L \right). \quad (57) $$

Then, if the mass of the rib is m and we assume to be in the linear domain, its kinetic energy reads:

$$ \mathcal{T} = \frac{1}{2} [m\dot{w}_L^2 + mL\dot{w}_L\dot{\theta}_L + (I + mL^2/4)\dot{\theta}_L^2], \quad (58) $$

with $I = mL^2/12$ the moment of inertia of the rib around its center of inertia C. Using Lagrange's equations, one obtains the following equations of motion:

$$ m\ddot{w}_L + \frac{mL}{2}\ddot{\theta}_L = -F_L, \quad \frac{mL}{2}\ddot{w}_L + \frac{mL^2}{3}\ddot{\theta}_L = -M_L. \quad (59) $$

Inverting Eq. (57) and eliminating ($F_L, M_L$) with Eq. (59), the dynamics of the rib thus verifies:

$$ m \begin{bmatrix} 1 & L/2 \\ L/2 & L^2/3 \end{bmatrix} \begin{bmatrix} \ddot{w}_L \\ \dot{\theta}_L \end{bmatrix} + \frac{EI}{L_l^3} \begin{bmatrix} 12 & -6L_l \\ -6L_l & 4L_l^2 \end{bmatrix} \begin{bmatrix} w_L \\ \theta_L \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}. \quad (60) $$

These equations of motion are rewritten under a dimensionless form by defining:

$$ \bar{w}_L = \frac{w_L}{L_l}, \quad \bar{t} = \frac{t}{T_0}, \quad \bar{\omega} = T_0\omega, \quad T_0^2 = \frac{mL^2L_l}{EI}, \quad \bar{L} = \frac{L_l}{L}, \quad (61) $$

where $\bar{L}$ is a small parameter. One obtains:

$$ \begin{bmatrix} \bar{L}^2 & \bar{L}/2 \\ \bar{L}/2 & 1/3 \end{bmatrix} \begin{bmatrix} \ddot{w}_L \\ \dot{\theta}_L \end{bmatrix} + \begin{bmatrix} 12 & -6 \\ -6 & 4 \end{bmatrix} \begin{bmatrix} \bar{w}_L \\ \theta_L \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \quad (62) $$

where now $\dot{\phi} = d \circ d\bar{t}$.

This enables to compute the two dimensionless natural frequencies of this system as a function of $\bar{L}$. Then, by keeping the first order term of a Taylor expansion in $\bar{L}$ and using Eq. (61), one obtains the two natural frequencies (in rad/s):

$$ \omega_1^2 \approx \frac{3}{T_0^2} = \frac{3EI}{mL_lL^2}, \quad \omega_2^2 \approx \frac{48}{\bar{L}^2 T_0^2} = \frac{16}{\bar{L}}\omega_1^2. \quad (63) $$

One can observe that if the rib's length L is large with respect to the one of the compliant link $Ll$ ($\bar{L} \ll 1$), it leads to $\omega_2 \gg \omega_1$: the second natural frequency is much larger than the first one. Moreover, the first natural frequency coincides
with the one computed in Sect. 3.1. As a consequence, in the low frequency band under interest of our systems, if $\bar{L} \ll 1$, one can neglect the high natural frequency of the system.

This is equivalent to neglect all terms of the order $o(\bar{L})$ and above in Eq. (62), leading to neglect the translational inertia ($\ddot{w}_L \simeq 0$) of the rib. As a consequence, $F_L \simeq 0$ and the loading of the link is reduced to the bending moment $M_L$ only. The first Eq. (60) leads to $w_L = L_l \theta_L / 2$, and the second leads to $\ddot{\theta}_L + \omega_1^2 \theta_L = 0$ with $\omega_1 = \omega_{rb}$ of Eq. (9).

# Appendix B. Details of the Derivation of the Equations of Motion of Various Models

## Appendix B.1. Lagrange's Equation for Multibody and 0D Models of Architecture 1

For the model C2 in Fig. 16, the kinetic energies of the system of Fig. 16b:

$$ \mathcal{T}_{\Sigma / 0} = \frac{1}{2} (m_3 \dot{u}_3^2 + m_1 (\dot{u}_\theta + \dot{u}_f)^2 + m_2 (\dot{u}_\alpha + \dot{u}_f + \dot{u}_3)^2), \quad (64) $$

and the potential energy is,

$$ \mathcal{U}_{\Sigma / 0} = \frac{1}{2} (k_3 u_3^2 + k_1 u_\theta^2 + k_2 u_\alpha^2). \quad (65) $$

Based on Eqs. (64) and (65), the Lagrange equation [50] in a matrix form leads to the equations of motion of the form (19) with:

$$ M = \begin{bmatrix} m_3 + m_1 + m_2 & m_2 + m_1 & m_2 \\ m_2 + m_1 & m_1 + m_2 & m_2 \\ m_2 & m_2 & m_2 \end{bmatrix}, \quad K = \begin{bmatrix} k_3 & 0 & 0 \\ 0 & k_1 & 0 \\ 0 & 0 & k_2 \end{bmatrix}, \quad (66) $$

where $q = [u_3, u_\theta, u_\alpha]^T$ represent the three degrees of freedom of the studied mass/spring model.

For the multibody model, denoted model C1 in Fig. 16, knowing that the kinetic energy is written in the form:

$$ \mathcal{T}_{i/0} = \frac{1}{2} M_i V_{Gi/0}^2 + \frac{1}{2} (\mathcal{I}_{Gi} \omega_{i/0}) \omega_{i/0} \quad (67) $$

for a beam $i$, where $I_{Gi}$ is moment of inertia with $I_i = M_i L_i^2 / 12$ the diagonal term in the beam reference frame, $V_{Gi/0}^2$ is the velocity with respect to the origin at the point $G_i$ defined in Fig. 16 and $\omega_{i/0}$ is the angular velocity with respect to the base. By linearising kinetic energies around the equilibrium position, ($\gamma = \theta = \alpha = 0$), for the three beams represented in Fig. 16a this leads to:

$$ \begin{align} \mathcal{T}_{1/0} &= \frac{L_1^2 \dot{\theta}^2}{6} + \frac{L_1 L_3 \dot{\gamma} \dot{\theta}}{2} + \frac{L_3^2 \gamma^2}{2} \nonumber \\ \mathcal{T}_{2/0} &= \frac{L_2^2 \dot{\alpha}^2 M_2}{6} + \frac{L_2^2 \dot{\alpha} \dot{\theta} M_2 \cos(\beta)}{3} + \frac{L_2^2 \dot{\theta}^2 M_2 \cos^2(\beta)}{6} \nonumber \\ &\quad + \frac{L_2 \dot{\alpha} \dot{\gamma} M_2 (-2L_2 \sin(\beta) + 3L_3)}{6} + \frac{L_2 \dot{\gamma} \dot{\theta} M_2 (-L_2 \sin(2\beta) + 3L_3) \cos(\beta)}{6} \nonumber \\ &\quad + \frac{\dot{\gamma}^2 M_2 (L_2^2 \sin^2(\beta) - 3L_2 L_3 \sin(\beta) + 3L_3^2)}{6} \nonumber \\ \mathcal{T}_{3/0} &= \frac{L_3^2 \dot{\gamma}^2 M_3}{6}, \tag{68} \end{align} $$

with $\mathcal{T}_{1/0}$, $\mathcal{T}_{2/0}$, and $\mathcal{T}_{3/0}$, representing the kinetic energies of the leading edge, trailing edge, and central beam, respectively. The potential energies expression is:

$$ \mathcal{U}_{\Sigma / 0} = \frac{1}{2} (\Gamma_1 \gamma^2 + \Gamma_1 \theta^2 + \Gamma_2 \alpha^2). \quad (69) $$

Based on Eqs. (68) and (69), the Lagrange equation of this system leads to Eq. (33) of Sect. 3.3.

## Appendix B.2. Lagrange's equations for Multibody and 0D Model of Architecture 2

The kinematic energy of the mass/spring, model D2, represented in Fig. 17a, is as follows:

$$ \mathcal{T}_{\Sigma / 0} = \frac{1}{2} (m_3 \dot{u}_3^2 + m_1 (\dot{u}_\theta + \dot{u}_3)^2 + m_2 (\dot{u}_\alpha + \dot{u}_3)^2), \quad (70) $$

while the potential energy as the same form as in Eq. (65). By applying the Lagrange equations, this leads to the equations of motion of the form (19) with:

$$ M = \begin{bmatrix} m_3 + m_1 + m_2 & m_1 & m_2 \\ m_1 & m_1 & 0 \\ m_2 & 0 & m_2 \end{bmatrix}, \quad K = \begin{bmatrix} k_3 & 0 & 0 \\ 0 & k_1 & 0 \\ 0 & 0 & k_2 \end{bmatrix}, \quad (71) $$

where $q = [u_3, u_\theta, u_\alpha]^T$. It can be observed that the matrices of Eq. (71) correspond to a combination of two models: one with the central beam and the leading edge (Eq. 72):

$$ M = \begin{bmatrix} m_3 + m_1 & m_1 & 0 \\ m_1 & m_1 & 0 \\ 0 & 0 & 0 \end{bmatrix}, \quad K = \begin{bmatrix} k_3 & 0 & 0 \\ 0 & k_1 & 0 \\ 0 & 0 & 0 \end{bmatrix}, \quad (72) $$

the second being one studied on the Sect. 3.2, (Eq. 73) and correspond to the central beam and the trailing edge:

$$ M = \begin{bmatrix} m_3 + m_2 & 0 & m_2 \\ 0 & 0 & 0 \\ m_2 & 0 & m_2 \end{bmatrix}, \quad K = \begin{bmatrix} k_3 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & k_2 \end{bmatrix}. \quad (73) $$

For the multibody model of Fig. 13, model B1, by following the same approach as in Appendix B.1, linearising the kinematic energy around the equilibrium position ($\gamma = \theta = \alpha = 0$), one can obtain:
$$
\begin{align}
\mathcal{T}_{3/0} &= \frac{1}{2} M_3 \frac{L_3^2}{3} \dot{\gamma}^2 \nonumber \\
\mathcal{T}_{2/0} &= \frac{1}{2} M_2 \left( \frac{L_2^2 \dot{\alpha}^2}{3} + 2 \dot{\alpha} \dot{\gamma} \left( -\frac{L_2^2 \sin(\beta)}{3} + \frac{L_2 L_3}{2} \right) \right) \tag{74} \\
&\quad + \dot{\gamma}^2 \left( \frac{L_2^2 \sin^2(\beta)}{3} - L_2 L_3 \sin(\beta) + L_3^2 \right) \nonumber
\end{align}
$$

$$
\mathcal{T}_{1/0} = \frac{1}{2} M_1 \left( \frac{L_1^2 \dot{\theta}^2}{3} + L_1 L_3 \dot{\gamma} \dot{\theta} + L_3^2 \dot{\gamma}^2 \right) \quad (75)
$$

Notices that the potential energies have the same form then
in Eq. (69).

Using the expression of Eqs. (74) and (75), the Lagrange
equation of this system leads to Eq. (33) of Sect. 3.4.

For the case of the leading edge, we fixed $\beta = 0$, so it leads
us to :

for $\lambda_3 = L_3 + \frac{L_{l3}}{2}$. And for the case where $\lambda_2 = L_2 + \frac{L_{l2}}{2}$:

$$
\begin{align}
\lambda_3 &= \frac{(2L_2 + L_{l2}) (3L_3 (L_2 + L_{l2}) - 2 (L_2^2 + 3L_2L_{l2} + 3L_{l2}^2) \sin(\beta_2))}{4 (L_2^2 + 3L_2L_{l2} + 3L_{l2}^2)}, \nonumber \\
m_3 &= \frac{4L_3^2 (3M_2 (L_2 + 3L_{l2})^2 + 4M_3 (L_2^2 + 3L_2L_{l2} + 3L_{l2}^2)) (L_2^2 + 3L_2L_{l2} + 3L_{l2}^2)}{3 (2L_2 + L_{l2})^2 (3L_2L_3 + 3L_3L_{l2} - 2 (L_2^2 + 3L_2L_{l2} + 3L_{l2}^2) \sin(\beta_2))^2}, \nonumber \\
m_2 &= \frac{4M_2 (L_2^2 + 3L_2L_{l2} + 3L_{l2}^2)}{3 (2L_2 + L_{l2})^2}. \tag{78}
\end{align}
$$

**Appendix B.3. Lagrange’s Equations for the**
**Model with two Edges Fixed on the Central Beam**
**(Architecture 2) with Torsional Spring at L₁/₂**

Here, the same configuration as in Appendix B.2 is stud-
ied, but with the torsional spring placed in the middle of the

compliant link. In this case, the second moment of area is
written as $I_i = (L_i + 3L_{li})^2/12$ according to the Huygens-
Steiner theorem. Therefore, the kinetic energies of the three
beams are as follows:

$$
\begin{align}
\mathcal{T}_{1/0} &= \frac{1}{2} M_1 \left( \frac{\dot{\gamma}^2 (2L_3 + L_{l3})^2}{4} + \frac{\dot{\gamma} \dot{\theta} (L_1 + L_{l1}) (2L_3 + L_{l3})}{2} + \frac{\dot{\theta}^2 (L_1^2 + 3L_1 L_{l1} + 3L_{l1}^2)}{3} \right), \nonumber \\
\mathcal{T}_{2/0} &= \frac{1}{2} M_2 \left( \frac{\dot{\alpha}^2 (L_2^2 + 3L_2 L_{l2} + 3L_{l2}^2)}{3} + \dot{\alpha} \dot{\gamma} \left( \frac{(L_2 + L_{l2}) (2L_3 + L_{l3})}{2} - \frac{2 (L_2^2 + 3L_2 L_{l2} + 3L_{l2}^2) \sin(\beta)}{3} \right) + \right. \nonumber \\
&\qquad \left. \dot{\gamma}^2 \left( -\frac{(L_2 + L_{l2}) (2L_3 + L_{l3}) \sin(\beta)}{2} + \frac{(2L_3 + L_{l3})^2}{4} + \frac{(L_2^2 + 3L_2 L_{l2} + 3L_{l2}^2) \sin^2(\beta)}{3} \right) \right), \tag{76} \\
\mathcal{T}_{3/0} &= \frac{1}{2} M_3 \left( \dot{\gamma}^2 \left( \frac{L_3}{2} + \frac{L_{l3}}{2} \right)^2 + \frac{\dot{\gamma}^2 (L_3 + 3L_{l3})^2}{12} \right). \nonumber
\end{align}
$$

The potential energy is the same tan before, Eq. (69).
Following the same method as in Sects. 3.3 and 3.4, the

equivalence between this model and the 0D model is given
by,

$$
\begin{align}
\lambda_1 &= \frac{2 \left( \frac{L_1^2}{3} + L_1 L_{l1} + L_{l1}^2 \right)}{L_1 + L_{l1}}, & m_1 &= \frac{M_1 (L_1 + L_{l1})^2}{4 \left( \frac{L_1^2}{3} + L_1 L_{l1} + L_{l1}^2 \right)}, \nonumber \\
\lambda_2 &= \frac{6 (2L_3 + L_{l3}) \left( \frac{L_2^2}{3} + L_2 L_{l2} + L_{l2}^2 \right)}{3 (L_2 + L_{l2}) (2L_3 + L_{l3}) - 4 (L_2^2 + 3L_2 L_{l2} + 3L_{l2}^2) \sin(\beta)}, \nonumber \\
m_2 &= \frac{M_2 (3 (L_2 + L_{l2}) (2L_3 + L_{l3}) - 4 (L_2^2 + 3L_2 L_{l2} + 3L_{l2}^2) \sin(\beta))^2}{36 (2L_3 + L_{l3})^2 \left( \frac{L_2^2}{3} + L_2 L_{l2} + L_{l2}^2 \right)}, \tag{77} \\
m_3 &= \frac{M_1 (L_1 + 3L_{l1})^2}{4 (L_1^2 + 3L_1 L_{l1} + 3L_{l1}^2)} + \frac{M_2 (L_2 + 3L_{l2})^2}{4 (L_2^2 + 3L_2 L_{l2} + 3L_{l2}^2)} + \frac{4M_3 (L_3^2 + 3L_3 L_{l3} + 3L_{l3}^2)}{3 (2L_3 + L_{l3})^2}. \nonumber
\end{align}
$$
## Resizing Equation for a NAV of Architecture 2 with 2 Wings

To obtain the expressions for $b_2$ and $L_{l2}$ as functions of $\bar{\omega}_2$ and $\mu_2$, it is first necessary to replace the stiffness of torsional springs $\Gamma_1$ and $\Gamma_3$ in Eq. (35) with their detailed expressions from Eq. (12), as well as the beam masses $M_1 = b_1 L_1 h_1 \rho$ and $M_2 = b_2 L_2 h_2 \rho$, with $\rho$ the SU-8 density. Then, by introducing the detailed expressions from Eq. (34) into Eq. (42), one can obtain:

$$b_2 = \frac{9L_3^2 L_1 \mu_2 b_1 h_1}{L_2 h_2 (3L_3 - 2L_2 \sin(\beta))^2} \quad (79)$$

$$h_3 = \sqrt{\frac{9L_3^4 \mu_3 \bar{\omega}_3^2 b_1 h_1^3 (3L_3 - 2L_2 \sin(\beta))^2}{L_1^3 L_{l1} \beta_1^4 b_1 h_1 (18L_3^2 \mu_2 - 3\mu_3(3L_3 - 2L_2 \sin(\beta))^2 + 2(3L_3 - 2L_2 \sin(\beta))^2)}} \quad (82)$$

**Acknowledgements** This work was supported by ANR-ASTRID NANOFLY (ANR-19-ASTR-0023) and the French AID (Defense Innovation Agency).

**Author Contributions** All authors contributed to the study conception and design. Conceptualisation, methodology, formal analysis and investigation were performed by Marguerite de La Bigne, Mathieu Colin, Marie Zwingelstein, Sébastien Grondel and Olivier Thomas. The first draft of the manuscript was written by Marguerite de La Bigne, Mathieu Colin, Sébastien Grondel, Olivier Thomas and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

**Funding** Open access funding provided by Arts et Metiers Institute of Technology. Not applicable.

**Data Availability** Not applicable.

## Declarations

**Conflict of interest** On behalf of all authors, the corresponding author states that there is no conflict of interest.

**Ethics approval and consent to participate** Not applicable.

**Consent for publication** Not applicable.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright

$$L_{l2} = \frac{L_1^2 L_{l1} b_2 h_2^3 (3L_3 - 2L_2 \sin(\beta))^2}{9L_3^2 L_2^2 \mu_2 \bar{\omega}_2^2 b_1 h_1^3} \quad (80)$$

Similarly, for the central beam, by introducing the expression of $M_1$, $\Gamma_1$, from Eq. (12), that of $M_3$ and $K_3$ from Eq. (48) into Eq. (35) then injecting this result into Eq. (42), one can obtain:

$$b_3 = \frac{L_1(L_3)b_1h_1(-18L_3^2\mu_2 + (3L_3 - 2L_2\sin(\beta))^2(3\mu_3 - 2))}{L_3h_3(3L_3 - 2L_2\sin(\beta))^2}, \quad (81)$$

holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

1. Srygley, R., & Thomas, A. (2003). Unconventional lift-generating mechanisms in free-flying butterflies. *Nature*, **420**, 660–4. http://doi.org/10.1038/nature01223

2. Wootton, R. J. (1991). The functional morphology of the wings of odonata. *Advances in odonatology*, **5**(1), 153–169.

3. Dickinson, M., & Muijres, F. (2016). The aerodynamics and control of free flight manoeuvres in drosophila. *Philosophical Transactions of the Royal Society of London Series B, Biological Sciences*. https://doi.org/10.1098/rstb.2015.0388

4. Ansari, S., & Knowles, K. (2008). Insectlike flapping wings in the hover part i: Effect of wing kinematics. *Journal of Aircraft*, **45**, 1945–1954. https://doi.org/10.2514/1.35311

5. Ellington, C. (1999). The novel aerodynamics of insect flight: Applications to micro-air vehicles. *Journal of Experimental Biology*, **202**(23), 3439–3448. https://doi.org/10.1242/jeb.202.23.3439

6. Shengjie, X., Hu, K., Huang, B., Deng, H., & Ding, X. (2021). A review of research on the mechanical design of hoverable flapping wing micro-air vehicles. *Journal of Bionic Engineering*, **18**, 1235–1254. https://doi.org/10.1007/s42235-021-00118-4

7. Petricca, L., Ohlckers, P., & Grinde, C. (2011). Micro- and nano-air vehicles: State of the art. *International Journal of Aerospace Engineering*, **20**(1), 1–17. https://doi.org/10.1155/2011/214549

8. Hu, K., Deng, H., Shengjie, X., Yang, G., & Sun, Y. (2024). Design of a bio-inspired, two-winged, flapping-wing micro air vehicle with high-lift performance. *Journal of Bionic Engineering*, **21**, 1191–1207. https://doi.org/10.1007/s42235-024-00486-7

9. Zou, Y., Zhang, W., Zhou, S., Ke, X., Cui, F., & Liu, W. (2018). Monolithic fabrication of an insect-scale self-lifting flapping-wing robot. *Micro & Nano Letters*, **13**(2), 267–269. https://doi.org/10.1049/mnl.2017.0730

10. Wood, R. (2007). Liftoff of a 60 mg flapping-wing mav. In *International conference on intelligent robots and systems* (pp. 1889–1894). https://doi.org/10.1109/IROS.2007.4399502.
11. Meng, K., Zhang, W., Chen, W., Li, H., Chi, P., Zou, C., Wu, X., Cui, F., Liu, W., & Chen, J. (2012). The design and micromachining of an electromagnetic mems flapping-wing micro air vehicle. *Microsystem Technologies, 18*, 127–136. https://doi.org/10.1007/s00542-011-1388-6

12. Rashmikant, K., & Ishihara, D. (2023). Iterative design window search for polymer micromachined flapping-wing nano air vehicles using nonlinear dynamic analysis. *International Journal of Mechanics and Materials in Design, 19*, 407–429. https://doi.org/10.1007/s10999-022-09635-4

13. Faux, D., Thomas, O., Cattan, É., & Grondel, S. (2018). Two modes resonant combined motion for insect wings kinematics reproduction and lift generation. *Europhysics Letters, 121*(6), Article 114883. https://doi.org/10.1209/0295-5075/121/66001

14. Faux, D., Thomas, O., Grondel, S., & Cattan, É. (2019). Dynamic simulation and optimization of artificial insect-sized flapping wings for a bioinspired kinematics using a two resonant vibration modes combination. *Journal of Sound and Vibration, 460*, Article 114883. https://doi.org/10.1016/j.jsv.2019.114883

15. Shang, J., Combes, S., Finio, B., & Wood, R. (2009). Artificial insect wings of diverse morphology for flapping-wing micro air vehicles. *Bioinspiration & Biomimetics, 4*, Article 036002. https://doi.org/10.1088/1748-3182/4/3/036002

16. Bao, X. Q., Bontemps, A., Grondel, S., & Cattan, É. (2011). Design and fabrication of insect-inspired composite wings for may application using mems technology. *Journal of Micromechanics and Microengineering, 21*(12), Article 125020. https://doi.org/10.1088/0960-1317/21/12/125020

17. Dargent, T., Bao, X., Grondel, S., Brun, G., Paquet, J., Soyer, C., & Cattan, É. (2009). Micromachining of an su-8 flapping-wing flying micro-electro-mechanical system. *Journal of Micromechanics and Microengineering, 19*, 085028-1-10. https://doi.org/10.1088/0960-1317/19/8/085028

18. Colin, M. (2021). Conception d'un nano-drone à ailes vibrantes. PhD thesis, École Nationale Supérieure d'Arts et Métiers.

19. Wu, K. (2023). Design, modeling, optimization and control of compliant mechanisms. PhD thesis, Centrale Lille, France.

20. Howell, L.L., Magleby, S.P., & Olsen, B.M. (eds.) (2013). *Handbook of compliant mechanisms*. Wiley.

21. Samadikhoshkhaa, Z., Saive, E., & Lipsett, M. G. (2025). A review of compliant mechanisms for contact robotics applications. *Robotics and Autonomous Systems, 186*, Article 104902. https://doi.org/10.1016/j.robot.2024.104902

22. Platl, V., & Zentner, L. (2024). Optimizing natural frequencies in compliant mechanisms through geometric scaling. *Mechanism and Machine Theory, 203*, Article 105822. https://doi.org/10.1016/j.mechmachtheory.2024.105822

23. Colin, M., Thomas, O., Grondel, S., & Cattan, É. (2020). Very large amplitude vibrations of flexible structures: Experimental identification and validation of a quadratic drag damping model. *Journal of Fluids and Structures, 97*, Article 103056. https://doi.org/10.1016/j.jfluidstructs.2020.103056

24. Fry, S., Sayaman, R., & Dickinson, M. (2003). The aerodynamics of free-flight maneuvers in drosophila. *Science (New York, NY)*, 300, 495–498. https://doi.org/10.1126/science.1081944

25. Chin, D. D., & Lentink, D. (2016). Flapping wing aerodynamics: from insects to vertebrates. *Journal of Experimental Biology, 219*(7), 920–932. https://doi.org/10.1242/jeb.042317

26. Liu, H., & Aono, H. (2009). Size effects on insect hovering aerodynamics: An integrated computational study. *Bioinspiration & Biomimetics, 4*, Article 015002. https://doi.org/10.1088/1748-3182/4/1/015002

27. Wootton, R. (1999). How flies fly. *Nature, 400*(6740), 112–113. https://doi.org/10.1038/21998

28. Faux, D., Cattan, É., Grondel, S., & Thomas, O. (2017). Coupling of two resonant modes for insect wing mimicking in a

flexible-wing nav and generate lift. In *SMASIS 2017 (ASME 2017) Conference on smart materials adaptive structures and intelligent systems*, Snowbird, USA, paper 3770.

29. Debeurre, M., Grolet, A., Cochelin, B., & Thomas, O. (2023). Finite element computation of nonlinear modes and frequency response of geometrically exact beam structures. *Journal of Sound and Vibration, 548*, Article 117534. https://doi.org/10.1016/j.jsv.2022.117534

30. Debeurre, M., Benacchio, S., Grolet, A., Grenat, C., Giraud-Audine, C., & Thomas, O. (2024). Phase resonance testing of highly flexible structures: Measurement of conservative nonlinear modes and nonlinear damping identification. *Mechanical Systems and Signal Processing, 215*, Article 111423. https://doi.org/10.1016/j.ymssp.2024.111423

31. Zheng, L., Hedrick, T., & Mittal, R. (2013). Time-varying wing-twist improves aerodynamic efficiency of forward flight in buterflies. *PloS one, 8*, Article e53060. https://doi.org/10.1371/journal.pone.0053060

32. Wehmann, H. N., Heepe, L., Gorb, S. N., Engels, T., & Lehmann, F. O. (2019). Local deformation and stiffness distribution in fly wings. *Biology Open*. https://doi.org/10.1242/bio.038299

33. Truong, H., Engels, T., Wehmann, H., Kolomenskiy, D., Lehmann, F. O., & Schneider, K. (2022). An experimental data-driven mass-spring model of flexible calliphora wings. *Bioinspiration & Biomimetics, 17*, Article 026003. https://doi.org/10.1088/1748-3190/ac2f56

34. Cote, B., Weston, S., & Jankauski, M. (2022). Modeling and analysis of a simple flexible wing-thorax system in flapping-wing insects. *Biomimetics*. https://doi.org/10.3390/biomimetics7040207

35. Wang, L., Zhang, H., Zhang, L., Song, B., Sun, Z., & Zhang, W. M. (2024). Analysis and actuation design of a novel at-scale 3-dof biomimetic flapping-wing mechanism inspired by flying insects. *Bioinspiration & Biomimetics*. https://doi.org/10.1088/1748-3190/ad94c2

36. Colin, M., Fernandez, P., Godoy-Diana, R., Thiria, B., Grondel, S., Cattan, É., & Thomas, O. (2024). A simple aerodynamic model of the lift generation of hovering apping wings. *Physical Review Fluids* (submitted).

37. Saya, D., Dezest, D., Welsh, A. J., Mathieu, F., Thomas, O., Leïchlé, T., Trolier-McKinstry, S., & Nicu, L. (2020). Piezoelectric nanoelectromechanical systems integrating microcontact printed lead zirconate titanate films. *Journal of Micromechanics and Microengineering, 30*, Article 035004. https://doi.org/10.1088/1361-6439/ab60bf

38. Zhang, C., & Rossi, C. (2017). A review of compliant transmission mechanisms for bio-inspired flapping-wing micro air vehicles. *Bioinspiration & Biomimetics, 12*, Article 025005. https://doi.org/10.1088/1748-3190/aa58d3

39. Howell, L. L. (2001). *Compliant Mechanisms*. John Wiley and Sons.

40. Meirovitch, L. (1967). *Analytical methods in vibrations*. MacMillan.

41. Guillon, S., Saya, D., Mazenq, L., Perisanu, S., Vincent, P., Lazarus, A., Thomas, O., & Nicu, L. (2011). Effect of non-ideal clamping shape on the resonance frequencies of silicon nanopacntilevers. *Nanotechnology, 22*(24), Article 245501. https://doi.org/10.1088/0957-4484/22/24/245501

42. Engels, T., Wehmann, H. N., & Lehmann, F. O. (2020). Three-dimensional wing structure attenuates aerodynamic efficiency in flapping fly wings. *Journal of The Royal Society Interface, 17*, 20190804. https://doi.org/10.1098/rsif.2019.0804

43. Touzé, C., Vizzaccaro, A., & Thomas, O. (2021). Model order reduction methods for geometrically nonlinear structures: a review of nonlinear techniques. *Nonlinear Dynamics, 105*, 1141–1190. https://doi.org/10.1007/s11071-021-06693-9
44. Thomas, O. (2025). Understanding, computing and identifying the nonlinear dynamics of elastic and piezoelectric structures thanks to nonlinear modes. In C. Touzé & A. Frangi (Eds.), *Model order reduction for design*. Springer series CISM courses and lectures: Analysis and control of nonlinear vibratory systems.

45. Debeurre, M., Grolet, A., & Thomas, O. (2023). Extreme nonlinear dynamics of cantilever beams: Effect of gravity and slenderness on the nonlinear modes. *Nonlinear Dynamics*, **111**, 12787–12815. https://doi.org/10.1007/s11071-023-08637-x

46. Auleley, M., Thomas, O., Giraud-Audine, C., & Mahé, H. (2020). Enhancement of a dynamic vibration absorber by means of an electromagnetic shunt. *Journal of Intelligent Material Systems and Structures*, **32**(3), 331–354. https://doi.org/10.1177/1045389X20957097

47. Auleley, M., Giraud-Audine, C., Mahé, H., & Thomas, O. (2021). Tunable electromagnetic resonant shunt using pulse-width modulation. *Journal of Sound and Vibration*, **500**, Article 116018. http://doi.org/10.1016/j.jsv.2021.116018

48. de La Bigne, M., Cattan, É., Itawi, A., Ghenna, S., Grondel, S., & Thomas, O. (2024). Optimization of the design and microfabrication of a biologically inspired nano-aerial flapping wing vehicle.

In *Proceedings of IEEE 19th international conference on nano/micro engineered and molecular systems (NEMS)*, Kyoto, pp. 1–6.

49. de La Bigne, M., Cattan, É., Ghenna, S., Zwingelstein, M., Grondel, S., & Thomas, O. (2024). A new design and optimization method of a nano-aerial vehicle with resonant vibrating wings. In *Proceedings of ISMA 2024, International conference on noise and vibration engineering*, Leuven.

50. Inman, D. (2001). *Engineering vibration*. Prentice Hall.