domain: spi phys

text:
Anh Le

Damien Doan

Samuel Faux

Eric Dupont

Sébastien Cattan

Grondel

Cattan

Doan

Damien Faux

Samuel Dupont

Eric Cattan

Sebastien Grondel

Modeling and simulation of the vertical take off and energy consumption of a vibrating wing nano air vehicle

Keywords: MEMS, Vibrating-wing nano air vehicle, PID, Bond Graph. I

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.

INTRODUCTION

Inspire of the attractive aerial capacities such as hovering, backward fly and recovering after sock, very few recent prototypes have the same geometric characteristics and masses as an insect [START_REF] Wood | Microrobot Design Using Fiber Reinforced Composites[END_REF][START_REF] Dargent | Micromachining of an SU-8 flapping-wing flying microelectro-mechanical system[END_REF][START_REF] Groen | Improving flight performance of the flapping wing MAV DelFly II International Micro Air Vehicle conference and competitions (IMAV 2010)[END_REF][START_REF] Park | Development of bio-mimetic composite wing structures and experimental study on flapping characteristics Robot[END_REF][START_REF] Lee | MEMS wing technology for a battery-powered ornithopter Micro Electro[END_REF]. This could be due to undesirable surface effects such as adhesion, friction and wear [START_REF] Wood | Microrobot Design Using Fiber Reinforced Composites[END_REF] which makes the traditional mechanisms and macro-scale manufacturing methods become inefficient for an insect scale. To overcome the problem, a prototype with kinematics similar to those of insects using MEMS manufacturing processes has been constructed [START_REF] Dargent | Micromachining of an SU-8 flapping-wing flying microelectro-mechanical system[END_REF][START_REF] Bontemps | Design and performance of an insect-inspired nano air vehicle Smart[END_REF][START_REF] Tai | Titanium-alloy MEMS wing technology for a micro aerial vehicle application Sensors[END_REF]. Large wing displacement and high power efficiency are two main objectives. In a previous paper [START_REF] Bontemps | Modeling and evaluation of power transmission of flapping wing nano air vehicle[END_REF], the study was focused only on the tergum actuation. This study extends the investigations on a full wings prototype providing a view of energy transmission in company with the dynamic analysis of the vibrating wing. The aerodynamic forces calculated from the wing torsion and bending motions could be used in a closed-loop vertical fly. Since the studied system combines the mechanical, the electromagnetic and the aerodynamics fields, a unique model Bond Graph formalism has been set up. This kind of formalism is widely used to model multiphysics systems and energetics efficiency. The first part of the paper is an overview of the prototype including actuation, transmission and wings. The second part answers to the question of the wing model. The aerodynamic forces are then established from retrieved wing motion. Finally, the global model is constructed with a Proportional-Integral-Derivative (PID) controller and the efficiency of the nano air vehicle is determined and discussed.

II. PRINCIPLE

A. Prototype structure

The prototype studied here is mainly composed of a flexible beam (tergum) with a wing at each extremity. An electromagnetic actuator glued to the thorax, positioned at the center of tergum plays a role of generating and maintaining the vibration. The compliant links transmit the vibration to the wings and affect the system efficiency. The prototype is realized in SU-8 EPOXY, which is commonly used material in MEMS technologies. Figure 1 presents an overall schematic view of the system.

B. The Word Bond Graph of the prototype

Since the prototype is considered as an incorporation of components from different engineering disciplines, Bond Graph formalism has the crucial advantage to represent it as a whole, figure 2. The model is similar to the well-known block diagram, with the major difference that the "bonds" which link the elements together represent bi-directional exchange of physical energy. Each bond depicts instantaneous flow of energy or power denoted by a pair of power variables called flow and effort. For example, the bond next to the battery block would present the flow of electrical energy and the power variable would be the current (i a ) and the voltage (U a ), whose product is power (P a ). In like manner, the magnet displacement (z ) and the electromagnetic force (F e ) are flow and effort of the corresponding power (P mechanic ). The model of first two blocks could be found in the previous work [START_REF] Bontemps | Modeling and evaluation of power transmission of flapping wing nano air vehicle[END_REF]. By understanding the dynamic phenomenon of the existing artificial wing, a modal model based on mechanical beam has been developed in the third block. The skeleton of wing is presented by an Euler Bernoulli beam as depicted in figure 3. The compliant links are as a simple supports in flexion and a springs in torsion. The coupling between flexion and torsion is simulated via the continuous conditions at the junction of the trailing edge and the leading edge. Bending and passive torsion motions extracted from the third block comes into the last one as signals to create the aerodynamic forces. Assumption is taken here that there is a low coupling between aerodynamic forces and structural forces, which means that the kinematics is essentially due to the structural part whereas the membrane induces the aerodynamic forces. This assumption has already been validated in a previous paper [START_REF] Dargent | Micromachining of an SU-8 flapping-wing flying microelectro-mechanical system[END_REF]. In the next part, we will take a closer look at the wing model.

III. MODAL MODEL OF THE WINGS

As mentioned above, the skeleton of wings is modeled by a set of Euler-Bernoulli beams numbered from one to five as in figure 4. The vibration in bending and torsion in each beam could be expressed by the two following equations:

+ =0 - =0 (1)
where w is the transverse displacement due to flexion, is the angle of torsion, EI and GI 0 are the bending and torsion stiffnesses, is the mass per unit length of the beam and I 0 is the moment inertia of rotation. The lengths of leading edge, trailing edge, and the distance between two supports are L, H, and l respectively. The method of separation of variables in time and space can be used to solve equation ( 1) by separating the spatial and temporal function as:

( , ) = ( ) ( ) ( , ) = ( ) ( ) (2)
which can be substituted into the flexion part of (1) to give:

( ) ( ) =- ( ) ( ) with = (3)
The left hand side of (3) depends on x alone while the right hand side depends on t alone. Since x and t are independent variables, both sides of (3) must be equal to the same positive constant :

( ) ( ) =- ( ) ( ) = (4)
yielding:

( ) - ( ) =0 with = (5) ( ) + ( ) =0 (6)
the solution forms of ( 5) and ( 6) are then:

( ) = + + ℎ + ℎ (7) ( ) = + (8)
where A, B, C, D, a and b are unknown constants determined by the boundary and initial conditions respectively.

Inserting (2) into the torsion part of (1) we get:

( ) ( ) = ( ) ( ) =- (9)
yielding:

( ) + ( ) =0 with = (10)
the solution form of (10) is:

( ) = ( ) + ( ) (11)
in which E and F are unknown constants determined by the boundary conditions. The following analysis applies the above solutions for each beam of the wing of type "L".

Suppose that all the segments have the same section and mechanical properties. Since the torsion has very little effect on the trailing edge (beam number 1 and 5), there are five spatial equations for flexion but only three for torsion:

′′′′( ) - ( ) =0 ′′′′( ) - ( ) =0 ′′′′( ) - ( ) =0 ′′′′( ) - ( ) =0 ′′′′( ) - ( ) =0 (12) ′′( ) + ( ) =0 ′′( ) + ( ) =0 ′′( ) + ( ) =0 (13)
It is a general mathematical principle that the number of boundary conditions necessary to determine a solution to a differential equation matches the order of the differential equation; therefore twenty-six limiting conditions are needed to solve (12) and (13). As can be seen that this set of equations has the form of the homogeneous system equations. Solve them for the value of β and the non-trivial solutions (A ,B ,C ,D ,E ,F ), we now obtain the spatial equations of bending and torsion:

( ) = + + ℎ + ℎ ( ) = ( ) + ( ) (14)
for 1≤i≤5 and 2≤j≤4. Figure 5 shows the plot of ( 14) for the first mode of beam model. The membrane of the wing has been added but has no influence on the kinematics as mention previously The experiments in [START_REF] Bontemps | [END_REF] demonstrate the rightness of this model. The bending and torsion angles and also the resonance frequency derived from the simulation are close to those of the experiment. However, this system is not totally compatible with Bond Graph formalism and need a little tweaking to make itself fit. In other word, the L-shaped wing is now identified as a spring-mass-damper system. The spring stiffness (1/k eq ) and the equivalent mass (m eq ) are characterized in bond graph formalism by a C-element (potential energy) and a I-element (kinetic energy) respectively. The value of a R-element (dissipated energy) presents the system damping (R eq ) could be found in [START_REF] Bontemps | [END_REF]. The spting stiffness k eq and the equivalent mass m eq are determined using the modal superposition [START_REF] Balachandran | Vibrations[END_REF] as the following integrals taken along all segments of the wing:

= ( ) (15) = ( ) ( ) (16)
The results are summarized in the table 1 for clarity sake. In the mechanical domain, the 1-junction represents force balance and is a generalization of the Newton's second law. All flows of the connecting bonds at a 1-junction have the same value and the efforts sum is equal to zero:

( ) + ( ) + ( ) = ( 2 ) (17)
Note that ( /2) is the electromagnetic force acting at the middle of two supports computed directly using the modal superposition principle and q(t) is the function of time as mentioned in [START_REF] Dargent | Micromachining of an SU-8 flapping-wing flying microelectro-mechanical system[END_REF]. Finally, a full prototype in Bond Graph language is shown in the figure 6. Recall that the electromagnetic force formula has the form:

= ( ).2 . . . = ( ). ( 18
)
where k em (z) depicts the output signal of a block diagram MultilyDivide element which is the multiplication of the mean magnetic flux density B m (z) with the coil radius r the coil number of turn n and the 2 value. A modulated gyrator (MGY) with its gyration ratio (k em (z)) symbolizes the transformation of electrical energy into mechanical ones. Battery is presented via a modulated source of flux i a (MSf) and a R element stands for by coil resistor (R coil ). To connect the electromagnetic part with the L-shaped beam part we also need another type of transformer which is a TF element with transform ratio ( ( /2)). For more detail of the electromagnetic actuator refers to [START_REF] Bontemps | Modeling and evaluation of power transmission of flapping wing nano air vehicle[END_REF].

IV. ENERGY EVALUATION

When the current source is connected to the prototype, energy is provided to compensate the energy dissipation and maintain the wing vibration. In fact, for an electrical input power of 1mW, corresponding to a current (rms) of 80 mA and a voltage (rms) of 12.5mA, the supplying energy in one cycle is:

= ×( ) (19)
The power transformation in the electromagnetic actuator costs a major part of the input energy because of the Joules effect on the coil resistor (E Rcoil ). The remaining part comes into the wings and will dissipate (E d ) or store (E s ) under different ways resulting in a mechanical energy (E mechanic ). In physic and engineering, a measure of quality of an oscillatory system is its Q, defined as:

=2 (20)
The Q factor is a dimensionless parameter that indicates the energy losses within a resonant element of a mechanical structure or an electronic circuit such as a spring or an inductor. The higher the value of Q the lower the rate of energy loss and hence the better the system is. In order to evaluate the prototype efficiency, we need to identify the two sources of energies mentioned in (20). The energy stored (E s ) corresponds to the one stored under kinematics (E I ) and potential energy (E C ). The energy dissipated (E d ) is the loss due to the resistor of coil in the electromagnetic actuator (R coil ) and the damping value (R eq ). The same analysis as in [START_REF] Bontemps | [END_REF] revealed that E Rcoil is responsible for 98% of energetics losses and therefore this work is mainly focus on the efficiency of new wing model based on remaining energy.

It is important to realize that the system is always excited to its first resonant frequency by an alternating current (i a =I 0 sin

). In this mode, the stored energy does not change and simply moves between the stored C and I Bond graph elements. When the flow is maximum across the I bond graph element (

), all the energy is in fact stored in the I element (kinetic energy):

= 1 2 (21)
Now let's introduce the dissipated energy in the system. At each cycle, the resistor R eq will dissipate energy:

= 1 2 2 (22)
Substitute the results of ( 22) and ( 21) into (20), we get an encouraging value 11.5 of Q. In other words, mechanical damping (R eq ) dissipates only 35% of the energy rest. It also means that almost two-third of the mechanical energy is useful for generating the wings motion. The energy analysis will be the key answer for the future question of battery's capacity and mass. The compromise between them both ensures a stable and long duration of fly. Figure 7 brings to light the energy distribution.

V. AERODYNAMIC AND VERTICAL FLY

In nature, the flying insects achieve flight through the use of flapping wings. Robot mimicking them flap and rotate their wing periodically with sufficient large angular amplitude to perform aerial maneuvers. Flight forces with great precision are required to remain hover. The wings model developed in the section III is adapted to this requirement. In hover, the ambient air velocity is negligible, the aerodynamic forces developing only through the wing relative motion. The bending and torsion angle velocities are determined as:

( , ) = ( ) * ( )

( , ) = ( ) * ( ) (23)
The aerodynamic forces is commonly composed of lift and drag components, F L and F D . Both have the same expression, the only difference being the coefficient, C L and C D .

= 1 2 ( ) ( ) = 1 2 ( ) ( ) (24)
where = -, r is radial distance from the wing root, and c is the wing chord length. Since we suppose that there is a low coupling between the aerodynamic and mechanical part, the aerodynamic block used here is instead a block diagram with the equations of aerodynamic forces inside and not a Bond Graph block. With the available wing's dimensions and the value of C L and C D found in [START_REF] Samuel | Bond graph model of a flapping wing micro-air vehicle[END_REF], the output lift and drag forces are calculated by the input q (t) as in equations ( 23) and (24). Integrating this block to the previous Bond Graph ones; we obtain a mean lift force of 650μN which is enough to elevate the prototype (26 mg). One actuator prototype shows the advantage of large and symmetrical wings motion, however, that means the lift force between each wing cannot be varied for a high maneuverability of a fly. In this case, a Bond Graph model of a stable vertical takeoff flight will be demonstrated as in figure 8. In order for a fly vehicle to rise into air, the lift force (F Lift ) created must be at least greater than or equal to the force of gravity (F Gravity ) as shown in equation ( 25):

+ = - ( 25
)
where z is the prototype's altitude and m prototype is the mass of prototype. Suppose that the prototype moves at relatively low speeds through the air, the viscous resistance f viscous is approximately proportional to its velocity:

=- ( 26
)
where b is the viscous coefficient. A 1-junction is employed to present the equation (25) as in the Vertical takeoff block. A control system is modeled for this prototype using PIDcontroller, where the input parameter is the unstable altitude signal of vertical takeoff block (z). The altitude control result and input current are shown in figure 9. The coefficients for the proportional, integral, and derivative term (K P , K I and K D ) are 9, 8 and 3 respectively. For a simpler model, the mean magnetic flux density B m (z) is approximated by linear signal and the lift forces is calculated with mean value of torsional angle. Remember that the prototype always runs in its resonant mode, controlling the amplitude of input current one can have stabilized system as well as a desired altitude. The prototype reaches its desired position of 50cm after 2.5 seconds.

VI. CONCLUSION AND PERSPECTIVES

This paper proposes a Bond Graph model of a VWNAV wings which has been experimentally validated and used to evaluate the whole system energy consumption. The low percentage of dissipated energy (E d ) due to mechanical damping (R eq ) comparing to one of the stored energy (E S ) shows a reasonable system efficiency. The vertical fly simulation reveals the prototype ability of flying up and being stable at a desired position in the vertical direction. In the next step, the strategy will be developed to reach 6 Degree Of Freedom and therefore allowing the control of VWNAV in the real space.

Fig. 1 .

 1 Fig. 1. Presentation of the prototype architecture.

Fig. 2 .

 2 Fig. 2. Word Bond Graph of the prototype.

Fig. 3 .

 3 Fig. 3. Diagram of the wing skeleton.

Fig. 4 .

 4 Fig. 4. Diagram of a wing of type "L"

Fig. 5 .

 5 Fig. 5. First mode of beam model.

Fig. 6 .

 6 Fig. 6. Bond Graph formalism of L-shaped wing.

Fig. 7 .

 7 Fig. 7. Energy analysis.

Fig. 8 .

 8 Fig. 8. Bond Graph model of closed loop fly.

Fig. 9 .

 9 Fig. 9. Altitude control and input current.

TABLE I

 I

	.	MECHANICAL PARAMETERS
		Used values
	k eq	4e-7 N/m
	m eq	1.4e-11 kg
	R eq	2e-10 kg/s
	Resonance frequency	27Hz

ACKNOWLEDGMENT

This work was realized within the ANR CLEAR-Flight project 13-0012-001 (Controlled lift for efficient flight of an Artificial Insect) which is supported by the French program ANR ASTRID.


Samuel Dupont
email: samuel.dupont@univ-valenciennes.fr

Sébastien Grondel

Alexandre Bontemps

Eric Cattan

Daniel Coutellier

Dupont Samuel

Grondel Sébastien

Bontemps Alexandre

Cattan Eric

Coutellier Daniel

Bond graph model of a flapping wing micro-air vehicle

Keywords: Bond Graph, MAV, dynamic model, flapping wings I

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Unmanned Aerial Vehicles (UAV's) enjoy the latest technology and seem ready to conquer the civilian and defence markets. In particular, a great interest is focused on micro-air vehicles (MAV's) and nano-air vehicles (NAV's) which expand the applications of more traditional UAV's due to their high manoeuvrability, portability and aeromechanical propulsion efficiency. The development of miniaturized structures, propulsion and navigation systems, aerodynamics, flight controllers, actuators and sensors are just a few of the current research topics in the field of MAV and NAV as shown in the last edition of the International Competition of Micro Aerial Vehicle and Flight competition (IMAV 2013) [1]. Among those topics, flapping MAVs are subject to a growing interest. They are inspired by insects and birds and fly by generating lift through oscillation or flapping of wings [START_REF] Tanaka | Flight dynamics of a butterfly-type ornithopter[END_REF]. The first active flapping MAV reviewed was the Terrestrial Entomopter: it was constructed of the Reciprocating Chemical Muscle (RCM) actuator and stereolithographic wings. Flight control was achieved through the actuation of the x-wing design of the system. The University of Delft [START_REF] Bolsman | Design Overview of a Resonant Wing Actuation Mechanism for Application in Flapping Wing MAVs[END_REF] proposed a MAV called the DelFly micro which was able to fly in outdoor and indoor for several minutes. This MAV only weighed 3 grams and had a size of 10cm from wing tip to wing tip. This makes it the smallest flying ornithopter carrying a camera in the world. From its side, Defence Advanced Research Projects Agency (DARPA) contributed to create a prototype "hummingbird-like" aircraft for a "NAV" program [START_REF] Keennon | Developpment of the nano hummingbird : a tailless flapping wing micro air vehicle[END_REF]. The result was called the Nano Hummingbird which can fly at 18km/h and move in three axes of motion. The artificial hummingbird manoeuvres using its flapping wings for propulsion and attitude control. It includes the systems required for flight: batteries, motors, and communications systems as well as the video camera payload. By scaling down prototypes around a wingspan of 3-5 cm, NAV's are scarce [START_REF]Pines and DARPA/DSO, Nano Air Vehicle Program[END_REF] and only the tiny fly by Wood et al. [START_REF] Wood | Microrobot Design Using Fiber Reinforced Composites[END_REF] is currently capable of flying under tethered power supply. Novel fabrication methods, design strategies as well as actuation concepts are here mandatory and micro-electromechanical systems (MEMS) technologies offer a very attractive solution by providing reliable, accurate and efficient microstructures such as required by a flapping-wing NAVs. To our knowledge, IEMN is the only ones to propose an entire MEMS structure called OVMI [START_REF] Bao | Design and fabrication of insect-inspired composite wings for MAV application using MEMS technology[END_REF][START_REF] Bontemps | Design and performance of an insect-inspired nano air vehicle[END_REF] able to produce large displacement and lift through the use of an electromagnetic actuator. The modelling and control of insect-like flapping wing MAV and NAV is widely studied and continuing to expand. Many of the research present in the literature focus on the standard aircraft six degrees of freedom equations found in many flight dynamics textbooks, e.g., the treatment by Etkin and Reid [START_REF] Etkin | Dynamics of Flight[END_REF]. For example, Doman et al. [START_REF] Doman | Dynamics and control of a minimally actuated biomimetic vehicle: Part i -aerodynamic model[END_REF][START_REF] Oppenheimer | Dynamics and control of a minimally actuated biomimetic vehicle: Part ii -control[END_REF] present modelling and control of a flapping wing MAV based on the RoboFly presented in [START_REF] Wood | The first takeoff_of a biologically inspired at-scale robotic insect[END_REF]. The aerodynamic model used in the simulations is developed in [START_REF] Hedrick | Flight Control in the Hawkmoth Manduca Sexta: The Inverse Problem of Hovering[END_REF] and based on the work of Sane and Dickinson in [START_REF] Dickinson | Wing rotation and the aerodynamic basis of insect flight[END_REF]. The standard aircraft equations of motion are justified here by assuming that either acceleration or the gyroscopic effects of the wings are small, compared to the body effects, or that the effects averaged over one flapping cycle are identically zero. An alternative approach is to use a combination of analysis and experimentation. For example, Perez et al [START_REF] Perez-Arancibia | First controlled vertical flight of a biologically inspired microrobot[END_REF] propose a two-stage methodology for tackling the problem: in the first stage, substantial a priori information about the system is gathered through system identification; in the second stage, the information collected through the static flapping experiments is employed to design the controller used in vertical flight. The design of the proposed controller relies on the idea of treating an exciting signal as a subsystem of the NAV. Since the controlled flight of our NAV [START_REF] Bontemps | Design and performance of an insect-inspired nano air vehicle[END_REF] requires optimizing the actuation in terms of power consumption, this paper introduces Bond Graph modelling to the field of flapping UAV's. Indeed, we chose to use language software adapted to mechatronic systems and energetic study. Our long-term goal is to build a full system description model, including the power actuation-transmission system [START_REF] Bontemps | Modeling and evaluation of power transmission of flapping wing nano air vehicle[END_REF] and the dynamic behaviour of the OVMI. As a first step, we focus in this paper on the Bond graph dynamics modelling of a MAV similar to the RoboFly in order to validate our modelling approach. To derive the equation system of the flapping wing MAV, we assume that the wings are rigid and that the MAV structure is symmetric. Figure 1 shows the coordinate system for the two models. The earth frame is denoted: E (E; x E , y E , z E ) and the body frame is denoted B (B; x, y, z). The right (left) wing frame is denoted W R (W R ; x R , y R , z R ) (W L (W L ; x L , y L , z L )), W R and W L being the articulation points of the right and left wings respectively. Another important frame is the aerodynamic frame A L,i ( A R,i ), which corresponds to the right (left) aerodynamic frame of the i th wing element with origin the aerodynamic centre of pressure. Assuming that the mass and inertial of the wings as well as the aerodynamic forces on the body are negligible, the equations of motion of the flapping wing MAV developed in the Newton-Euler formalism can be expressed as:

II. Flapping Wing MAV

A. Global dynamics equations E

B EB B ext V Ω V F m m ⊗ + = • (1) EB EB EB ext Ω Ω Ω M b b I I ⊗ + = • (2)
with m, I b being the body mass and inertia tensor about its respective centre of mass, whereas V B and Ω EB are the translational velocity with respect to the earth frame and the angular velocity of the body with respect to the earth frame expressed in the body frame respectively. Moreover, F ext and M ext stand for the sum of external forces and moments applied at the centre of mass, namely: gravitational and aerodynamic efforts.

g F F aero ext m + = ( 3
)
aero ext M M = (4)
To express the dynamic resultant F aero and moment M aero of the aerodynamic forces, the local aerodynamic forces and moments must be converted from the aerodynamic frame to the earth frame. [START_REF]Pines and DARPA/DSO, Nano Air Vehicle Program[END_REF] where F aeroR and F aeroL denote the aerodynamic forces within the earth frame, whereas D R,i , D L,i and L R,i , L L,i represent the drag and lift components on the i th chord-length element of the right and left wings respectively. Here R EB (resp. R BW ) is the rotation matrix between the earth frame and the body frame (resp. between the body frame and the wing frame) whereas R WAi corresponds to the rotation matrix between the wing frame and the aerodynamic frame of the i th wing element.

∑ ∑ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ + ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ = + = i i L i L i L A L W L BW EB i i R i R i R A R W R BW EB L D L D , , , , , , 0 0 R R R R R R F F F aeroL aeroR aero
⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ ⊗ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ + ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ + ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ ⊗ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ + ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ = ⊗ + ⊗ = + = ∑ ∑ → → i L i L i L A L W i i AL i AL WL WL WL B L W L BW i R i R i R A R W i i AR i AR WR WR WR B R W R BW L R L D y x z y x L D y x z y x BA BA , , , , , , , , , , 0 0 0
0 R R R R R R F F M M M
aeroL aeroR aeroL aeroR aero [START_REF] Wood | Microrobot Design Using Fiber Reinforced Composites[END_REF] Here R WB is the rotation matrix between the wing frame and the body frame whereas (x W , y W , z W ) and (x A,i , y A,i , 0 ) correspond to the coordinates of the articulation wing point and the aerodynamic pressure center of the i th wing element. Note that the indices R and L for lift and right wings respectively are deliberately removed in the following to simplify the expressions.

B. Aerodynamic forces

As shown on figure 1, the model of aerodynamic forces is based upon a local bi-dimensional approach: each wing is decomposed into a given number of parallel, chord-length elements, and the aerodynamic characteristics such as velocity V and angle of attack α are locally calculated for each element using aerodynamic 2D models. For a given i th wing element, the expression of the local aerodynamic velocity V i is given by [START_REF] Rakotomamonjy | Evaluation of flexoble flapping concept[END_REF]:

⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ ⊗ - ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ + ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ ⊗ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ - ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ - = ⊗ - ⊗ - ⊗ - - = → → → 0 0 , , , , i A i A i A i A W W W WB WB WE i i y x y x z y x r q p w v u WA WA BW BW BW EB EB B i Ω R R R Ω Ω Ω V V (7)
Where A i is the aerodynamic centre of the i th element and Ω BW the rotation vector between the wing and the body defined with the derivatives of the kinematic angles ξ, λ, µ, ν. More precisely, ξ is the angle of the stroke plane, λ the flapping angle, µ the deviation angle from the stroke plane and ν the pitch angle of the wing around its longitudinal axis. They correspond to the four possible rotations of the wing with respect to the body. Then the local aerodynamic angle of attack α i can be expressed as a function of the components of V i within (

)

) 1 j ( j 2 , , - = - - ∠ = z i x i i V V α ( 8
)
Both components D i , L i of the aerodynamic forces can be represented within the aerodynamic frame of the i th element as:

() () () [ ] () () () [ ] ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ + + = + + = + + = + + = i mas x i i rot x i i sta x i i mass i rot i st i i i mas z i i rot z i i sta z i i mass i rot i sta i i C C C S D D D D C C C S L L L L α α α ρ α α α ρ , , , 2 , , , 2 2 1 2 1
i i V V (9)
with the following aerodynamic coefficient functions respectively for the stationary, rotational, mass lift and drag:

( ) ( ) () ( ) () ( ) ( ) ( ) ( ) ( ) ⎪ ⎩ ⎪ ⎨ ⎧ = = ⎩ ⎨ ⎧ + = + = ⎪ ⎩ ⎪ ⎨ ⎧ = - = sin 2 cos 1 1 0 i i i i V V V V λ α π λ α π ν α π ν α π α α c y C c y C c c l C c c l C k C k k C Ai i mas x z Ai i mas x i x i rot z i x i rot x i i Cz sta z i i Cx Cx sta
x i [START_REF] Doman | Dynamics and control of a minimally actuated biomimetic vehicle: Part i -aerodynamic model[END_REF] Where S i , c, and l x represent the surface, the chord and distance between the leading edge and the axis of rotation. k Cxo , k Cx1 , k Cz1 , are constants based on measurements. The elementary aerodynamic forces are then summed up and integrated within the earth frame, in order to obtain the translational and rotational velocities of the global MAV.

C. Bond Graph modeling

Figure 2 displays the bond graph schematic for the flapping wing MAV. In the actuation block, a mathematical description of the rotations angles of the wings are programmed (ξ, λ, µ, ν). It feeds the kinematics blocks for the right and left wings respectively. This block allows computing the rotation matrix, angular and translational velocities. For example, the rotation matrix R WB between the wing frame and the body frame, as well as the rotation matrix R EW between the earth frame and the wing frame which are inputs of the aerodynamic velocity block. Then the velocity block solves equation 7 whereas the aerodynamic forces block computes equation 5, 6 and 8-10. Note that only the left side of the model is depicted in figure 2 as it is a symmetric model. Aerodynamic forces and velocities are also side dependant, whereas the dynamic block where the Newton-Euler equations are solved is common to both wings as all moments and forces (right and left) are taken into account. Then the Bond graph representation of the Newton Euler equations (1)-( 2) is shown in Figure 3 where the upper triangle relates the energy flow in translational dynamics and the lower triangle defines rotational dynamics. The three 1-junctions both in the upper and lower triangles provide the nodes for external forces and moments respectively.

III. Results and Discussion

To validate the bond graph model, open loop simulations have been performed. It is worth mentioning that the open loop system is unstable, as a consequence, in the model development, we have introduced the possibility to lock each of the 6 degrees of freedom of the system (3 in translation, 3 in rotation). The first step is to lock all of them as in this static configuration; a comparison with experimental results on flapping wing robots is possible [START_REF] Dickinson | Wing rotation and the aerodynamic basis of insect flight[END_REF]. Then, the study of motion along one axe is possible, such as elevation / descent when the z E translational degree of freedom is unlocked.

A. Validation of the aerodynamic model

Aerodynamic velocities and aerodynamic efforts have been calculated and compared to measurements available in the literature. A mathematical description of flapping wing insects was programmed to mimic wing kinematics, namely: flapping and pitch angle. In the case of wings insects µ = 0, moreover, the case considered here for validation, a flapping robot, requires to take ξ = 0. A typical evolution of the derivatives: dλ/dt and dν/dt is displayed Figure 4. The periodic laws define the input signals for our simulations; the pitch angle is a square like function as it has been shown that incidence is quite constant during stroke. They have been parameterised in order to be able to compare our simulations with the Robotfly experimental results [START_REF] Dickinson | Wing rotation and the aerodynamic basis of insect flight[END_REF]. Figures 5 and6 represent the evolution of the local aerodynamic horizontal (V xi ) and vertical (V zi ) velocities as a function of time, each curve being the result obtained for one of the 10 slices the wing has been devised into (i=1 represents the slice of the wing the closer to the wing hinge point). It is observed that aerodynamic speeds are increasing as a function of distance along the wing as expected. This results in the evolution of the aerodynamic velocity norm, figure 7. The evolution of lift forces are depicted by figure 8, it represents the steady lift force and the rotational lift forces. Both forces depend on the local aerodynamic norm. Moreover, it can be observed that the rotational lift force evolution depends also on the pitch angle derivative. In figure 9 evolution of the total lift force is plotted. Both qualitative evolution and quantitative value of experimental results are well reproduced [START_REF] Dickinson | Wing rotation and the aerodynamic basis of insect flight[END_REF].

B. Dynamic model

One degree of freedom has then been released: translation along z E -axis. This allows us to study the displacement, position, velocity, along this axis. In particular we want to verify the influence of the forces produced by the flapping of the wing on the MAV motion. Figure 10 displays the evolution of the velocity w along z Eaxis, it corresponds to a free fall of the MAV: only the weight is taken into account. Figure 11 illustrates that, when the flapping frequency is no more null, the total lift force produces a force that is sufficient to slow down the falling of the MAV.

IV. Conclusion

In this paper we have introduced bond graph modelling to the field of flapping wing MAV aerodynamics. This model has been validated thanks to experimental published results. Bond graph modelling enables an efficient approach of complex and coupled dynamics modelling, as illustrated by this paper. Flapping wing MAV is a mechatronic system, as such, the concern about different energies domains is crucial: in our approach, next step is to merge this model with the model of power transmission of a flapping wing nano air vehicle [START_REF] Bontemps | Modeling and evaluation of power transmission of flapping wing nano air vehicle[END_REF].

It is well know that the bond graph approach also facilitates extension of the model to include control concerns, future extensions of the model will be connected with regulation laws. Moreover, thanks to the discrete nature of the wing model, this model can be improved by taking into account in a future development the flexibility of a real wing.

Fig. 1 .

 1 Fig. 1. Flapping wing MAV schematic with reference frames.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Bond graph schematic of the flapping wing MAV (left side).

Fig. 4 .

 4 Fig. 4. Flapping and pitch angular velocities.

Fig. 5 .

 5 Fig. 5. Local aerodynamic horizontal velocity (Vxi)

Fig. 6 .Fig. 7 .Fig. 8 .Fig. 9 .

 6789 Fig. 6. Local aerodynamical horizontal velocity (Vzi)

Fig. 10 .

 10 Fig. 10. Evolution of the ascention velocity, open Loop Free Fall Simulation for zero thrust.

Fig. 11 .

 11 Fig. 11. Evolution of the ascention velocity, open Loop Free Fall Simulation with sinusoïdal flapping.