Nonlinear dynamics of a flapping rotary wing:Modeling and optimal wing kinematic analysis

Qiuqiu WEN,Shijun GUO,Ho LI,Wei DONG

aSchool of Aerospace,Beijing Institute of Technology,Beijing 100081,China

bCentre for Aeronautics,SATM,Cran field University,Cran field MK43 0AL,UK

1.Introduction

The Micro Air Vehicle(MAV)has become an active research area due to the potentiality for the civil and military application.1The typical characteristics of MAV are small dimension(wing spans within 15 cm),low weight(gross take-off weight ranging from 100 to 200 g)and low flight speed(between 10 and 15 m/s).In recent two decades,a variety of MAV layouts,which mainly include fixed wing,rotary wing,and flapping wing,had been put forward.However,due to the extremely small dimension and high lift and efficiency requirements at

Nomenclature

low Reynolds numbers Re,few practical MAVs with load carrying capabilities has been accomplished.Research efforts for new and practical designs of MAVs have never been stopped.

In 2004,Vandenberghe et al.2employed the experimental method and found that a pair of wing flapping up and down can freely rotate spontaneously around the horizontal shaft as a critical frequency was exceeded.Based on this discovery,Guo et al.3,4proposed the design of Flapping Rotary Wing(FRW) flight vehicle as a new con figuration of MAV.Similar concept was also proposed and applied in full-scaled helicopter rotor by Van Holten et al.5As shown in Fig.1,a pair of antisymmetrically mounted wings,which can flap along the vertical direction by a drive shaft,is fixed on the rotary rigid base.The thrust generated by the wings’vertically flapping motion drives them to rotate around the shaft,resulting in a flapping and simultaneously rotating kinematics.Combined with tuning the pitch angles of the wings asymmetrically in the upstroke and down-stroke,the high lift force is produced to make FRW take-off and hover.

Recently,experimental works6were used to measure the force produced and proved that the lift from flapping rotary wing was larger than that from conventional rotary wing in the range of Re from 2600 to 5000.Wu et al.7conducted a computational fluid dynamics method to research the unsteady aerodynamic behavior of FRW.It is observed that the leading edge vortex attached on the wing surface during the whole flapping period,which is the main reason for the high lift generation by FRW.Unlike the ordinary Flapping Wing(FW)flight vehicle which is driven by its own motor to rotate,the flapping rotary wing is driven by the aerodynamic force to rotate passively.Previous works on FRW have mostly assigned a constant rotation velocity by assuming an‘equilibrium’state.However,for a practical wing,the inertia forces associated with the complicated kinematics will essentially interact with the aerodynamic force production.The influence of the wing inertia and the dynamic process as the wing converges to the equilibrium status will necessarily have a nontrivial effect on the aerodynamic performance of FRW.The varying rotation velocity conversely affects the aerodynamics and flow structure of the flapping rotary wing,8resulting in a coupling between the passively rotary motion and aerodynamic force.Therefore,the nonlinear dynamic model,especially for FRW,is needed to analyze its aerodynamic performance.

To date,only a few studies have focused on the dynamics of FRW.However,many studies have been relevant to flapping wing flight vehicle.In these studies,dynamics of FW are generally investigated using standard aircraft equations with six degrees of freedom.9,10However,this approach neglects the inertial effects of the mass of the wings.Recently,some studies have investigated the dynamics from the aspect of multiple body nonlinear system,such as Gebert11and Sun12et al.Orlowski and Girard13modeled a flapping wing micro air vehicle as a system of three rigid bodies,a body and two wings,and studied the influence of the mass of the wings to the dynamics.Mahjoubi and Byl14developed the dynamic multi-body model using Lagrangian method and the proposed control approach to optimize the wings’mass and mechanical impedance properties of the joints.These studies have indicated that the multiple-body dynamic theory may be used to analyze the dynamics of FRW.

In this paper,a simplified FRW is modeled as three rigid bodies,one for the body of rotary base and others for each wing.The wing pitching motion is assumed to be actively driven through a control servo,as shown in Fig.1.Thereby,each flapping wing owns three degrees of freedom:the actively flapping,pitching and the passively rotating.Using the D’Alembert’s Principle given in Ref.15,a multi-body dynamic model is derived for FRW.In addition,a quasi-steady aerodynamic model is utilized for the calculation of the aerodynamic forces and moments.The motion process of wings is simulated in a selected typical parameter set to understand the coupling with the lift/thrust production.Finally,the effects of the kinematics of wings on the dynamic rotational equilibrium of FWR and the aerodynamic performances are presented.

2.Reference definition

To describe the motion of rotary base in the FWR body frame,and the motion of wings with respect to rotary base,four reference frames are used.The body frame Obxbybzbis attached to the center of the body of FWR.As shown in Fig.2,the positive xbaxis is along the longitudinal axis of the central body.The ybaxis locus in the vertical symmetry plane of body and is perpendicular to the xbaxis with a positive upward.The zbaxis is perpendicular to the xOy plane.The unit vectors of the body frame are presented by ex，b，ey，band ez，b.

After rotating an angle ψrabout the ybaxis of the body frame for the rotary base,it becomes the rotary plane frame Orxryrzr(shown as the subscript ‘r”).The rotate plane frame de fines the rotary motion of two connected wings.

The wing- fixed frames Owxwywzware two fixed frames attached to the wings.The initial orientation of the wing fixed frames is parallel to the rotate plane frame with an origin coincident with the rotation of the wings joint.The orientation of the wings with respect to the rotary plane is determined by the pitch angle ϑwand flap angle γwof the wings.Here we use the subscripts L,R to represent the left and right wings,respectively.

The rotation matrix from the body frame to rotation plane frame is

As shown in Fig.2,the wings successively rotate about the xrand zwaxis with the angles of γwand ϑwto reach the ultimate position.The rotation matrices for the right wing are

The rotation matrices for the left wing,with respect to the rotary base,are combined in the same manner as in Eq.(2).The only difference that the signs of γwand ϑware interchanged for the right and left wings.

Combining Eqs.(1)and(2),the rotation matrices from the body frame to right and left wing- fixed frames are

3.Dynamic model of FRW

3.1.Method and assumption

The flapping wing vehicle is modeled as a system of three rigid bodies:a central body of rotary base with two rigid wings attached at ideal hinges.The method chosen to derive the equations of motion is D’Alembert’s Principle Extended to Multiple Rigid Bodies.15The functions of generalized inertia force are described as

where i presents the number of rigid bodies and j denotes the number of generalized coordinates.γijrepresents velocity coefficients.βijare angular velocity coefficients.The inertia forceand momentof the ith rigid body are given as

where mi，vi，ρi，ωiand Iidenote the mass,velocities,reference vectors,angular velocities and the resulting mass moments of inertia matrices of the ith rigid body.

In this study,the two wings are assumed to be attached to the rotary base body by joints that allow two degrees of freedom respectively with a common rotation degree of freedom.To simplify the derivation, firstly,the inertia tensors for the individual bodies are calculated with respect to the reference point and they do not need to be calculated at the timevarying center of mass of the system.Then,the body of FRW is assumed to be always fixed on the ground,thereby its motions relative to inertial space are neglected and the body frame is equal to inertial frame.As a result,a dynamic system of three rigid bodies:one for the rotary base,the other two for each wing,are considered.The five degrees of freedom are selected to be described by the generalized coordinates xi,listed together as

The related quasi-velocities of coordinates,expressed in inertia frame,are

The variables ωx，ωyand ωzdescribe the angular velocity of the each selected center rigid bodies in the body frame.Especially, ωy，rdenotes the rotation angular velocity of the rotary base with two wings,and the flapping and pitching angular velocities of each wing are expressed as[ωz，wRωx，wR]and[ ωz，wLωx，wL],respectively.

3.2.Velocities and reference vectors

The angular velocity vector of the wing related to rotary base and expressed in rotation plane frame can be obtained by the time derivative of the two Euler angles˙γwand˙ϑw,which are assumed to be known as the command input.For right wing,the equation is defined as

Related to body mass center,the joint point owns an angular velocity

With the combination of Eqs.(8)and(9),the angular velocity of the right and left wings with respect to the body frame,and expressed in the body frame,are

The reference vectors denote the position of the center of mass of the ith body with respect to the reference point.For the rotary base,the reference point is chosen to be its respective center of mass,thereby the reference vector ρ1equals zero.In each wing- fixed frame,the positon of mass center owns two components along xwaxis and zwaxis directions:

The related reference vectors are transformed from the wing- fixed frames according to

Since the translational velocity of the rotary base,the reference velocity v1equals zero.And for each of the wings,the reference point of translational velocity is its joint point.The vectors from the center of rotary base to wing joint points,are expressed as rwRand rwL.As shown in Fig.1,the vectors defined in body frame own two components along ybaxis and zbaxis directions:

The reference velocity,for each of the wings,is the velocity of the respective wing joint in the inertia frame.The velocities of the wings are

The acceleration can be derived by differentiating the above equation

3.3.coefficients

The angular velocity coefficients βijare necessary for the derivation of dynamic model,which arise from the calculation of virtual work performed by moments.Each coefficient is vector and is determined for each rigid body and velocity combination.The angular velocity coefficients are defined as

The coefficients of center body of FWR are

The angular coefficients of right and left wings are

3.4.Mass moments of inertia

For the rotary base,the mass symmetry for xOy and xOz planes is assumed.No planes of mass symmetry are assumed for either wing during the model development.As a result,the resulting mass moments of inertia matrices for each rigid body are

3.5.Motion equations of rotary base

The derived equations of passive rotation motion,with all of the individual pieces together,are presented in vector notation.

The rotations of the right wing and the left wings are described by Eqs.(18)and(19),respectively.

Here,Q1is the rotational moment acted on rotary base.The generalized forces Q2，Q3are the control moment for the right wing,and Q4，Q5are the control moments for the left wing.

The rotation moment,expressed in body frame,can be dived as aerodynamic moments Maeroproduced by flapping wings and gravity moments Mmassdue to the mass of wings.In this study,the quasi-steady theory is used to calculate the aerodynamic moments Mwproduced by flapping wings.The calculation model is given in the following chapter.For each wing,the Mmassis calculated according to

where RIbdenotes the transfer matrix from inertial frame to body frame.As the assumption of this study,we have RIb=I.As a result,the Mmassalong the ybaxis equals zero.That means the gravity of wings will not produce the rotation moment,if the body of FWR does not have angular motion in inertial frame.Then,the rotation moment has an expression as

4.Aerodynamic model

In the numerical study of Wu et al.7on FRW,a strong spanwise flow on the wing was observed,and the LEV on the FRW wing merged with the tip vortex and the Trailing Edge Vortex(TEV),forming a vortex ring structure that stayed attached on the wing throughout the flapping cycle.These findings suggest that the quasi-steady model used in this study is applicable for modeling the aerodynamic forces of FRW.As a result,in this study,we extended the quasi-steady aerodynamic model to the application of the flapping and simultaneously rotating wing kinematics of FRW.

Firstly,a geometric model of the FWR wing is chosen and the detailed shape and definition of geometric parameters of the wing are given in Appendix A.For blade element analysis,it is convenient to write down the velocity and acceleration of a 2D wing chord due to the gyration of the wing at span-wise location r.The resultant velocity and acceleration vector when expressed in the wing- fixed frame are planar vectors with only two nontrivial in dices,i.e.the xwand ywcomponents:

and

where ex，w，ey，wand ez，wdenote the unit vectors of right wing frame;ωwis the angular velocity vector of the right wing related to inertia frame and expressed in wing frame.For two wings,the related ωwRand ωwLare obtained as

Since the velocity and acceleration of wing are expressed in wing frame,the effective Angle of Attack(AOA)of the wing αecan be easily found by inverse trigonometric function of the velocity components ratio of the wing:

In quasi-steady theory of flapping wing,this relationship holds,except that the force vector acts perpendicular to the wing chord.16The quasi-steady forces are divided by translational forces,rotational forces and virtual mass forces,and the corresponding coefficients are experimentally measured.17–19Here,we will use this definition.The corresponding equations for translational force Ftand rotational force Frfor two wings are expressed as

And

where ρ is the density of the surrounding air.Cris the rotational force coefficient due to wing pitching,and the value of this coefficient is chosen as Cr=1.6 in our calculation.The vector Ctis the translational force coefficient and can be treated as a unit force vector acting on the wing.In the velocity direction,Ctis expressed as life coefficient CLand drag coefficient CD,and can be approximated by the following equations20:

where the constant coefficients CLmax，CDmaxand CD0at the specific Reynolds number(Re≈4000)are valued from 3D CFD case calculations result.The values are given as:CLmax=0.18,CDmax=3.4,CD0=0.05.

Since we calculate the force and moment in the wing frame,the translational force coefficient Ctcan be obtained as

where CHis the translational force coefficient along xwaxis,and CVis the translational force coefficient along ywaxis.

For the calculation of the aerodynamic torque,the location of the Centre of Pressure(CP)at a chord-wise location r is defined as rCP=xCPex，w+rez，w.As a result,the aerodynamic torque of the above two forces can be decided by the following equation:

The virtual mass force and moment are calculated using Sedov’s formula,18which is suitable for our coordinate definition:

where vyand ωzcan be obtained from the vectors ωw，vwof the each wing. λaand λaωare the added mass force coefficients,which are obtained as

After integrating Eqs.(22),(28),(31)and(32)along the wing span orientation,as a result,the total aerodynamic forces and moments,expressed in body frame,are obtained as

The forces and moments for each wings are calculated based on its velocity and angular velocity respectively.Then,the necessary aerodynamic moments Maero Rand MaeroLin Eq.(22)are obtained.

5.Simulation conditions

5.1.Kinematic functions of wings

In this study,we use simple harmonic functions to describe the flapping and pitching motion of the wing,as previous studies for insects flight.21,22The kinematic functions of the wing is specified by giving the variation functions:

where fFis the flapping frequency,Δγwis the flapping amplitude angle,and Δα is the pitching amplitude;for modeling the asymmetric pitching,the angle α0is introduced.By this definition,the calculation functions between Δα，α0with the geometric AOA of the wing at mid up-stroke angle αUand mid down-stroke angle αDare given as

The time history of flapping motion and pitching motion is plotted as an example in Fig.3 to illustrate the relationship between flapping motion γwand pitching motion

5.2.Nondimensional coefficients

We used the mean chord length of the wingand the mean flapping velocity at the wingtip vt=2ΔγwfFR as the reference length and reference velocity.The aerodynamic lift and rotation moment coefficients are thus defined as

Here,Fyand Mxaeromean aerodynamic lift and rotational moment,¯CLandare the averaging aerodynamic lift and rotational moment coefficients during one flapping period,which are obtained by

where t0is the initial time at the beginning of one flapping period and TFmeans the flapping period.If the effect of geometric shape on aerodynamics is not considered,then the flapping period TFonly depends on frequency fF.The energetic cost of the FRW’s wings can be calculated by the time averaged power efficiency coefficient over a flapping period TF.For a practical MAV design,elastic storage is desirable for energy efficiency,of which the order is decided by the design property of the mechanical system.

In the current study,we consider that the mechanical system of the FRW can fully store the input power.The instantaneous aerodynamic power efficiency coefficient for hovering flight due to gyration equals directly minus dot product of the angular velocity vector ωwwith the aerodynamic torque

then,the average power efficiency coefficient is given by

Anondimensional variable¯t=t/TFis defined to describe clearly the time courses of wing motion during a flapping period.As shown in Fig.3,¯t∈[0，0.5]indicates that the wing’s motion is in the phase of up-stroke,whereas ¯t∈ [0.5，1.0]indicates that the wing’s motion is in the phase of down-stroke.

5.3.Validation

To validate the compatibility of the aerodynamic model,we compare the calculation results of our model with 3D CFD results presented by Wu et al.,7who studied the aerodynamic characteristic of FRW at a low Reynolds Number.In Wu’s work,7the CFD mode employed a boundary fitted dynamic grid to orientate the wing boundary at a different time with the prescribed kinematics.An OH type mesh was used for flow simulation.Grid 1 has dimensions of 31×33×37(in normal,chord wise,and span wise directions,respectively);and grids 2 and 3 have dimensions of 51×57×61 and 81×81×91,respectively.The outer boundary for these grids is located 30c away from the wing surface and 15c away from the wingtip.The first grid spacings from the wing surface of the three grids are 0.002,0.001,and 0.0005.

In the benchmark case,the flapping and pitching motions may be defined by the previous descriptions.In the example in CFD results of Wu et al.,7the dynamic of rotation motion is ignored and the rotation speed is assumed to be constant,and we use the same model to describe the rotation motion in this case.

where ψj0is the rotating speed.

Since the kinematics of FWR is combined by steady rotation and reciprocal flapping motion,a nondimensional rotation speed kRis defined to measure the deflection of the effective AOA:

Based on Eq.(43),ψj0can be obtained.

The necessary parameters in Eqs.(35),(37)and(43)for validation case are Δγw=30°,fF=22 Hz,αU=-30°,αD=0°,kR=0.25.The Reynolds number for flapping flight may be defined by Re=vt¯c/v,where v is the kinematic viscosity of the air.In this case,we have Re=4058.As shown in Fig.4,a reasonable agreement is achieved between the lift and rotation moment trends obtained from the present model and CFD results.The averaging coefficients of the present model areandwhileand0.168 in CFD results.Thereby,the case study shows that the model employed in the present study is credible.

6.Results and analysis

6.1.Analysis of a typical case

A typical case(Δγw=30°,fF=22 Hz,αU=-30°,αD=20°,Re=4058)is selected and discussed in this section to investigate the rotation performance of the FRW with the structural parameters listed in Table 1.

The result of the rotation angular velocity varying with time is presented in Fig.5,and the rotational moment coefficient CRand its averaging value¯CRin each flapping periods are given in Fig.6.It can be clearly seen that the aerodynamic rotational moment in the up-stroke decreases to negative value with the increase of the rotational velocity,and the aerodynamic average rotational moment finally decreases close to zero after fifteen flapping periods,and then maintains at a small value continually.It differs from other existing kinematics with prescribed wing motion,such as rotary wing and insects flapping wing.

However,note that¯CRof each period cannot converge to zero strictly and ωyroscillates continuously with an amplitude of±0.5 r/s,even after a lot of flapping periods.It is caused by inertia coupling phenomenon of two wings,which are not symmetrical about the axis of rotation.According to Table 1,if the inertial products Ixy，w，Ixz，wand Iyz，ware small and ignored,according to Eqs.(17)and(22),the inertia coupling rotation moment Mcouplcan be expressed as

Table 1 Basic parameters of FRW.

Similar to coefficients CRand,the coefficients of Mcouplare defined as

Fig.7 gives the result of rotation coefficients CRcouplandin each flapping period.After the comparison with Fig.7,it can be seen that the inertia coupling rotation moment always exists even because of the high-speed flapping and pitching motions of two wings.As a result,the rotational moment is balanced with not only the resistance drag of the fluid,but also inertia coupling moment.However,the inertia coupling moment is small compared with the aerodynamic moment,and the total moment is still zero.Consequently,an FWR would reach and stay in an equilibrium rotation speed.In this study,we define the motion status ofnear zero as the Equilibrium Rotational Status(ERS)for FRW.

Fig.8 gives the results of lift coefficients CLandin each flapping period.As a comparison,constant rotational velocity model with kR=0.47 is used to calculate the same case.It can be seen from Fig.8(a)that the lift coefficient vastly increases with the increasing rotational velocity,especially in the upstroke,where the large negative lift becomes positive at the rotational equilibrium state.As a result,as shown in Fig.8(b),the variance ofis increasing gradually and reaches the equilibrium value finally.Only in this time,the lift force of FRW keeps stable with certain kinematic parameters input.Notably,the lift generation of FRW presents first-order inertia system characteristics for a new kinematics of wings input because of the inertia damping in passive rotation motion.That is different from the motions of rotary wing and insect flapping wing,the life force of whom changes and reaches stable immediately with the variance of wings kinematics.4,13Thereby,the dynamic inertia time-lag phenomena of lift generation due to passively induced rotational velocity is a unique feature of the FRW con figuration.

In order to assess the dynamic system for lift generation,a time-lag constant parameter τaerois defined to present the time cost of arriving ERS.In this case,the τaerois obtained as 0.68 s,which equivalent to 15 flapping period.

The following figures present a comparison between the simulation results ofandwith the full,2 times,4 times and 1/2 times rotational moment of inertia of wings Iy，w.Fig.9 shows that as the rotational moment of inertia of the wings relative to the central body increase,the inertia time-lag phenomena of lift generation becomes serious.τaerochanges from 0.41 s,corresponding to 1/2 times Iy，w,to 2.72 s,corresponding to 4 times Iy，w.However,once FWR reaches and stays in an equilibrium rotational status,the lift generation is stable and not affected by the variance of Iy，w.Thus,a small rotational moment of inertia of wings is useful to decrease time-lag constant.

6.2.Effects of the kinematics of wings

As given in Eq.(36),the pitching kinematics of wings is described as sinusoidal wave functions.For a certain flapping frequency fF,the pitch angle of the wing at any instantaneous time is decided by AOAs parameters:mid up-stroke αUand mid down-stroke αD.The two parameters,donated pitching kinematics,may be selected to analyze the in fluence of wings kinematics to aerodynamic performance.

In order to present the aerodynamic performance of FRW in the equilibrium rotational status,the period average lift coefficientpower efficiency coefficientand nondimensional rotational velocityare defined as

where istabmeans the flapping period while FRW has been in the ERS,j is the total number of flapping period used to calculate coefficients,andP¯fstabis defined as positive always.In this case,we let j=10 to acquire accurate description of aerodynamic performance for FRW.

In a typical case of Δγw=30°,fF=22 Hz,Re=4058,Fig.10 presents the results of coefficientsand τaero,while the mid-up stroke αUvaries from 0°to-90°with increments of-2.5°,and the mid-down stroke αDis fixed to αD=10°,20°,30°.

As shown in Fig.10(a),notice that the lift of FRW is near zero when the αUequals-αD,e.g.the lift produced in up-stroke and down-stroke may cancel each other.In order to acquire a high positive lift generation of FRW,we need to increase αUand decrease αD.Since in the up-stroke,the high value of αUmakes the wing surface close to the air flow direction decided by rotation motion and flapping motion together.As a result,the negative effective AOA of wing αedecreases,or even changes to be positive.The negative lift generated in upstroke phase decreases consequently.On the contrary,the small value of αDincreases the αeat down-stroke phase,which produces more positive lift force.The power efficiency does not increase with the lift performance improving,while compared with Fig.10(a)and(b).Since the high drag force follows as high life generation,the optimal lift generation does not correspond to the optimal power efficiency.In addition,a high equilibrium rotation velocity occurs at small AOA parameters while both αUand αDare within 20°to-20°.However,as shown in Fig.10(c)and(d),highmeans more time is required to arrive ERS,thereby the τaerobecomes large for a small αUand αD.

Table 2 Maximum nondimensional parameters and the corresponding geometric AOAs.

In order to present the effect of pitching kinematic on FRW,three zones are defined,which include high lift force zonehigh power efficiency zoneand low time-lag zone(τaero≤ 0.6 s).The maximum nondimensional parameters and the corresponding geometric AOAs are listed in Table 2.Fig.11 shows three zones distribution as αUvaries from-80°to 0°and αDvaries from 0°to50°.The maximum lift coefficientresponses the point at αU=-20.0°, αD=6°,and the maximum power efficiency coefficient happens at αU=-53°, αD=34°.If taking the high lift,high power efficiency and low time-lag into consideration,as given in Fig.11,an optimal AOAs parameters occur at the point where αDis relative small(within 7–15°),and αUis in the range-35°to-40°.

7.Conclusions

The passive rotation motion and aerodynamic performance of a rotary base with two flapping wings as a simplified model of FRW flight vehicle are studied.The nonlinear,multiple body equations of motion for an FRW is derived using D’Alembert’s Principle for Multiple Rigid Bodies.In addition,a quasi-steady aerodynamic model is utilized for the calculation of the aerodynamic forces and moments at a low Reynolds number(Re≈4000).

The simulation of typical case shows that the passive rotation motion of FRW is a continuous dynamic process of convergence into rotary velocity equilibrium status due to an interaction between aerodynamic thrust and rotation velocity.That causes the unique time-lag of stable lift generation.Even in the rotation equilibrium status,the lift force is still oscillating with small amplitude.

The pitching kinematics of wings greatly affects the equilibrium rotational characteristics,thus the aerodynamic performance of FRW.The lift force generation,power efficiency,equilibrium rotational velocity and dynamic time-lag are studies for various AOA parameters of wings.The result shows that in order to acquire a high positive lift generation with high power efficiency and small dynamic time-lag,a relative high mid up-stroke αUand low mid down-stroke αDare necessary.In the zone where αDis within 7–15°and αUis within-35°to-40°,the performance of FRW is optimal.

Appendix A.

The geometry of the FRW wing is modeled by keeping the morphological parameters of quasi-static analysis similar with available insects’data.In generation,we assume that the thickness of the wing is small enough and has little effect on the wing’s aerodynamic.Thereby,a 2D wing with a span length R is shown in Fig.A1.Along the span direction,the wing is dived by in finite small strip with width dr.The chord length c(r)is shown at a chord-wise location r;h refers to the coordinate of the major axis of the ellipse in xwaxis.The Center of Gravity(CG)of wing has two components xCG，zCGin the wing frame.The AR is within 3–5;the shape parameters,including wing aspect ratio AR,the first,second and third radii of nondimensional moment of wing areaandare:AR=3.6,0.63.

The chord-wise location of Centre of Pressure(CP)is in the xwaxis.Referring to Dickinson et al.,17the xCPwith respect to location r varies linearly with the change of the effective AOA αe(given in Eq.(26)and the linear relation can be expressed as

References

1.Pines DJ,Bohorquez F.Challenges facing future micro-air-vehicle development.J Aircr 2006;43(2):290–305.

2.Vandenberghe N,Zhang J,Childress S.Symmetry breaking leads to forward flapping flight.J Fluid Mech 2004;506:147–55.

3.Guo S,Yang D,Huang Z.A smart material aeroelastic flapping wing micro rotorcraft.Int Forum Aeroelasticity Struct Dynam;2009.

4.Guo S,Li D,Wu J.Theoretical and experimental study of a piezoelectric flapping wing rotor for micro aerial vehicle.Aerosp Sci Technol 2012;23(1):429–38.

5.Van Holten T,Heiligers M,Kuiper R.Forced flapping mechanism designs for the Ornicopter:a single rotor helicopter without reaction torque.In:European rotorcraft forum;2004.

6.Zhou C,Wu J,Guo S.Experimental study on the lift generated by a flapping rotary wing applied in a micro air vehicle.Proc Inst Mech Eng Part G – J Aerospace Eng 2014;228(11):2083–93.

7.Wu J,Wang D,Zhang Y.Aerodynamic analysis of a flapping rotary wing at a low Reynolds number.AIAA J 2015;53(10):1–16.

8.Li H.Aerodynamic analysis and experiment of a micro flapping wing rotor[dissertation].Bedford:Cran field University;2015.

9.Etkin B,Reid L.Dynamics of flight:stability and control.New York:John Wiley&Sons Inc;1996.

10.Khan Z,Agrawal S.Modeling and simulation of flapping wing micro air vehicles.In:Proceedings of IDETC/CIE’2005 2005 ASME international design engineering technical conferences;2005 September 24–28;Long Beach,CA,USA.New York:ASME;2005.

11.Gebert G,Gallmeier P,Evers J.Equations of motion for flapping flight.In:Proceedings of the AIAA atmospheric flight mechanics conference.Reston:AIAA;2002.

12.Sun M,Wang J,Xiong Y.Dynamic flight stability of hovering insects.Acta Mech Sin 2007;23(3):231–46[Chinese].

13.Orlowski CT,Girard AR.Modeling and simulation of nonlinear dynamics of flapping wing micro air vehicles.AIAA J 2011;49(5):969–81.

14.Mahjoubi H,Byl K.Dynamics of insect-inspired flapping-wing MAVs:multibody modeling and flight control simulations.In:Proceedings of 2014 American control conference(ACC);2014 June 4–6;Portland,Oregon,USA.2014.

15.Greenwood DT.Advanced dynamics.New York:Cambridge University Press;2006.

16.Usherwood JR,Ellington CP.The aerodynamics of revolving wingsI.Modelhawkmoth wings.J Exp Biol2002;205(11):1547–64.

17.Dickinson MH,Lehmann FO,Sane SP.Wing rotation and the aerodynamic basis of insect flight. Science 1999;284(5422):1954–60.

18.Wang ZJ,Birch JM,Dickinson MH.Unsteady forces and flows in low Reynolds number hovering flight:two-dimensional computations vs robotic wing experiments.J Exp Biol 2004;207(3):449–60.

19.Sane SP,Dickinson MH.The control of flight force by a flapping wing:lift and drag production.J Exp Biol 2001;204(15):2607–26.

20.Whitney JP,Wood RJ.Aeromechanics of passive rotation in flapping flight.J Fluid Mech 2010;660(1):197–220.

21.Wu JH,Sun M.Unsteady aerodynamic forces of a flapping wing.J Exp Biol 2004;207(7):1137–50.

22.Bin Abas MF,Bin Mohd Ra fie AS,Bin Yusoff H,Bin Ahmad KA.Flapping wing micro-aerial-vehicle:kinematics,membranes,and flapping mechanisms of ornithopter and insect flight.Chin J Aeronaut 2016;29(5):1159–77.