Optimal pitching axis location of flapping wings for efficient hovering flight
Wang, Qi; Goosen, Hans; van Keulen, Fred DOI
10.1088/1748-3190/aa7795 Publication date
2017
Document Version
Accepted author manuscript Published in
Bioinspiration & Biomimetics: learning from nature
Citation (APA)
Wang, Q., Goosen, H., & van Keulen, F. (2017). Optimal pitching axis location of flapping wings for efficient hovering flight. Bioinspiration & Biomimetics: learning from nature, 12(5), [056001].
https://doi.org/10.1088/1748-3190/aa7795 Important note
To cite this publication, please use the final published version (if applicable). Please check the document version above.
Copyright
Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy
Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.
This work is downloaded from Delft University of Technology.
efficient hovering flight
Q. Wang, J.F.L. Goosen & F. van Keulen
Department of Precision and Microsystems Engineering, Structural Optimization and Mechanics Section, Delft University of Technology, Mekelweg 2, 2628CD, Delft, The Netherlands
E-mail: q.wang-3@tudelft.nl
Abstract. Flapping wings can pitch passively about their pitching axes due to their flexibility, inertia, and aerodynamic loads. A shift in the pitching axis location can dynamically alter the aerodynamic loads, which in turn changes the passive pitching motion and the flight efficiency. Therefore, it is of great interest to investigate the optimal pitching axis for flapping wings to maximize the power efficiency during hovering flight. In this study, flapping wings are modeled as rigid plates with non-uniform mass distribution. The wing flexibility is represented by a linearly torsional spring at the wing root. A predictive quasi-steady aerodynamic model is used to evaluate the lift generated by such wings. Two extreme power consumption scenarios are modeled for hovering flight, i.e., the power consumed by a drive system with and without the capacity of kinetic energy recovery. For wings with different shapes, the optimal pitching axis location is found such that the cycle-averaged power consumption during hovering flight is minimized. Optimization results show that the optimal pitching axis is located between the leading edge and the mid-chord line, which shows close resemblance to insect wings. An optimal pitching axis can save up to 33% of power during hovering flight when compared to traditional wings used by most of flapping wing micro air vehicles (FWMAVs). Traditional wings typically use the straight leading edge as the pitching axis. With the optimized pitching axis, flapping wings show higher pitching amplitudes and start the pitching reversals in advance of the sweeping reversals. These phenomena lead to higher lift-to-drag ratios and, thus, explain the lower power consumption. In addition, the optimized pitching axis provides the drive system higher potential to recycle energy during the deceleration phases as compared to their counterparts. This observation underlines the particular importance of the wing pitching axis location for energy-efficient FWMAVs when using kinetic energy recovery drive systems.
Keywords: flapping wings, hovering flight, passive pitching, pitching axis location, efficiency
Abbreviations
AOA angle of attack
BPDF Beta probability density function
CFD computational fluid dynamics
CWAD chordwise area distribution
CWMD chordwise mass distribution
FWMAV flapping wing micro air vehicle
HM-wing hawkmoth wing
KERS kinetic energy recovery drive system
LE leading edge
PA pitching axis
QE-wing quarter-ellipsoidal wing
Rect-wing rectangular wing
SWAD spanwise area distribution
SWMD spanwise mass distribution
TE trailing edge
Tri-wing triangular wing
1. Introduction
There exists a great diversity of wing morphology among insects, especially from the wing shape perspective. Compared to the wing outlines of various winged-insects, the wings of most flapping wing micro air vehicles (FWMAVs) (de Croon et al., 2009; Bolsman et al., 2009; Keennon et al., 2012; Ma et al., 2013; Nguyen et al., 2015) have much simpler shapes (e.g., quarter-ellipse, polygon). Although this simplification helps to decrease the design complexity and fabrication difficulty, the effect of wing shape on the flapping flight performance, particularly on the energy-efficiency, can not be completely ignored.
Wing shape determines the area distribution which can be further decomposed into distributions in the span and chord directions. Aerodynamic loads on flapping wings roughly increase quadratically with the velocity and, thus, are sensitive to the spanwise area distribution (SWAD). As a result, its influence on flapping wing performance has been extensively studied in the past. For instance, Ansari et al. (2008) investigated the influence of the spanwise wing area distribution and aspect ratio on the lift generation and power efficiency using an unsteady aerodynamic model (Ansari, 2004). Recently, Shahzad et al. (2016) conducted similar investigations using a computational fluid dynamics (CFD) method. They both showed that wings with more outboard area and higher aspect ratio generate higher lift but also consume more power. This is because the lift and power scale with the square and cube of flapping velocity, respectively. Optimal wing shapes have also been studied in literature using quasi-steady aerodynamic models. Typically, the total lift force is constrained in these studies. By keeping the wing aspect ratio constant, Wang et al. (2013) showed that flapping wings with more outboard area
can have better performance, i.e., higher lift generation or less power consumption. Stewart et al. (2013) found that the optimal wing shape which provides maximal thrust for forward flight tends to push much of the planform away from the wing root. Wang et al. (2014) investigated the optimal wing shape for energy-efficient hovering flight in passive pitching motion while keeping the total wing area constant. A wing planform with a low aspect ratio, similar to butterfly wings, was obtained.
In contrast, the influence of chordwise area distribution (CWAD) on flapping flight performance has been less studied. CWAD determines the inertia and aerodynamic loads a wing needs to overcome when pitching about its pitching axis. Particularly when the wing pitches passively due to its flexibility and the inertial and aerodynamic loads CWAD directly determines the passive pitching behavior, including the amplitude, the phase lag from sweeping motion, and its deviation from a harmonic motion. Passive pitching can be found for both insects (Ennos, 1989; Combes and Daniel, 2003; Bergou et al., 2007; Ishihara et al., 2009) and FWMAVs (Wood, 2008; de Croon et al., 2009; Keennon et al., 2012; Ma et al., 2013). Particularly for FWMAVs, the passive pitching motion is preferred in order to simplify the drive mechanism. Consequently, the location of the pitching axis becomes an important design aspect for artificial wings in order to achieve energy-efficient flight for FWMAVs.
This work aims to find the optimal pitching axis location for flapping wings with passive pitching motion. Flapping wings are modeled as rigid plates with non-uniform mass distribution. The wing flexibility is represented by a linearly torsional spring at the wing root, which is commonly used for the wing design (Whitney and Wood, 2010; Wang et al., 2014). A predictive quasi-steady aerodynamic model (Wang et al., 2016) is used to evaluate the lift generated by such wings. A short description of this model is provided in Appendix A. Two extreme power consumption scenarios are modeled for hovering flight, i.e., the power consumed by a drive system with or without the capacity of kinetic energy recovery. In this study, hovering flight is considered not only because it is generally more energy-consuming compared to forward flight (Dudley, 2002), but also it is a required and an important capability for FWMAVs.
The following sections start with the flapping wing modeling in Sec. 2 which consists of the wing area and mass distribution models as well as the kinematics model. In Sec. 3, the aerodynamic model and the power consumption model are introduced. Based on the optimization model described in Sec. 4, we investigate the optimal location of pitching axis for different wings and analyse its influences on the power efficiency during hovering flight in Sec. 5. Conclusions are presented in Sec. 6.
2. Flapping wing modeling 2.1. Area distribution
Wing area distribution is decomposed into spanwise area distribution (SWAD) and chordwise area distribution (CWAD).
Ellington (1984) proposed that the SWAD of insect wings can be approximately described by Beta probability density function (BPDF). The mean and standard variance of the BPDF are represented by the dimensionless radii of the first and second spanwise moment of area, which are denoted as ˆrs1 and ˆrs2, respectively. They are
defined as ˆ rs1 = 1 SR Z R 0 crdr and rˆs2 = s 1 SR2 Z R 0 cr2dr, (1)
where R, c and S are the wing span, chord length and the total area, respectively. Furthermore, ˆrs1 was correlated with ˆrs2 by the relation ˆrs2 = 0.929ˆr
0.732
s2 by Ellington
(1984). In contrast, for wings of FWMAVs, the SWAD is often simplified to reduce the modelling and fabrication complexity. For instance, flapping wings with rectangular (Seshadri et al., 2013), quarter-ellipsoidal (Keennon et al., 2012), triangular (Chaudhuri et al., 2013) and polygonal (de Croon et al., 2009) shapes have been adopted.
Compared to the SWAD, the CWAD was less studied for two reasons. First,
successful take-off and stable flight are still the research objectives of most FWMAVs. Thus, people are more interested in the lift and thrust generation which are more sensitive to the wing SWAD. Second, the influence of the CWAD on flight performance is more pronounced for wings that pitch passively when compared to wings with fully prescribed kinematics, and the latter is more frequently used. In order to model the
CWAD, we introduce a dimensionless parameter ˆd which is defined as the
local-chord-normalized distance from the leading edge (LE) to the pitching axis (PA). The parameter ˆ
d is formulated as a linear function of the spanwise radius r, i.e., ˆ
d(r) = r
R( ˆdt− ˆdr) + ˆdr, 0 6 r 6 R, (2)
where ˆdrand ˆdtrepresent the values of ˆd at the wing root and the tip, respectively. When ˆ
dr = ˆdt = 0, the wing takes its straight LE as the pitching axis, which has been widely adopted for most FMWAVs (de Croon et al., 2009; Bolsman et al., 2009; Keennon et al., 2012; Ma et al., 2013; Nguyen et al., 2015). When ˆdr = ˆdt = 0.5, the wing planform is symmetric about its PA, which has been used as a simplified wing model for studying the optimal kinematics by Berman and Wang (2007). In contrast, the PA of insect wings is generally located between the LE and the mid-chord line. The difference of the PA location between artificial and insect wings inspired us to investigate the influence of the PA on the flight efficiency and find the optimal PA location for artificial wings.
In this work, four commonly used wing shapes with straight leading edges (i.e., ˆ
dr = ˆdt = 0) but different spanwise area distributions (SWADs) are studied. As shown in Figs. 1 (a) - (d), the SWADs of these wings are identical to the rectangular wing (Rect-wing), the hawkmoth wing (HM-(Rect-wing), the quarter-ellipsoidal wing (QE-wing) and the triangular wing (Tri-wing), respectively. However, the shapes of these wings can vary significantly if the CWAD is changed. The resulting effect on the flight performance is still unclear and will be studied in this work. To facilitate the comparison, all wings have the same wing span and surface area S.
root LE tip TE
(a) Rect-wing (b) HM-wing (c) QE-wing (d) Tri-wing
Figure 1: Four traditional wings with the same straight leading edge (LE) but different trailing edges (TEs). The shape of the TE is determined by the spanwise area distributions (SWADs). The SWADs of these wings correspond to the rectangular wing (Rect-wing), the hawkmoth wing (HM-wing), the quarter-ellipsoidal wing (QE-wing) and the triangular wing (Tri-wing). The values of ˆrs1 are 0.5, 0.45, 0.42
and 0.33 for wings from left to right.
LE tm TE LE tm TE
(a) kite profile (b) uniform profile
Figure 2: Schematic diagram of two cross-sectional profiles that are used to describe the chordwise mass distribution. The maximal thickness of the cross-section is denoted as tm. LE and TE represent the leading edge and trailing edge, respectively.
2.2. Mass distribution
Similar to the wing area distribution, the wing mass distribution is also characterized by the spanwise mass distribution (SWMD) and the chordwise mass distribution (CWMD). Span wise, the dimensionless radii of the first and second moment of inertia, which are denoted as ˆrm1 and ˆrm2, respectively, are generally used to reflect the mass
distribution. According to studies on different species of insects (Ellington, 1984; Betts and Wootton, 1988; Dudley and Ellington, 1990; Willmott and Ellington, 1997b), ˆrm1
typically ranges from 0.2 to 0.5, and ˆrm2 from 0.3 to 0.6. Combes and Daniel (2003)
measured the spanwise bending stiffness variation from the root to the tip of a hawkmoth wing (Manduca sexta) and a dragonfly wing (Aeshna), and found that the variation can be approximated by an exponential decline. The spanwise stiffness is related to the cross-sectional profile in chordwise direction and proportional to the term c(r)t3
m(r),
where c and tm are the chord length and the maximal thickness of the cross-section.
Therefore, if the variation of the chord length is ignored, the change of tm along the span can be approximately described by an exponential function
tm(r) = λ1eλ2r, 0 6 r 6 R, (3)
where λ1 and λ2 are two parameters to be determined. In this work, we assume that
the wing is made of the same material with a density ρw, which implies that the mass distribution is equivalent to the wing thickness variation. Therefore, the function of tm can reflect the SWMD if the wing has the same cross-sectional profile along the span.
Table 1: Comparison of the calculated and measured values of ˆrm1 for different insect wings.
The value of ˆrm1 of the hawkmoth wing was measured by Willmott and Ellington
(1997b). The values of ˆrm1 of the other three wings were measured by Ellington
(1984). The measurement errors are not taken into account. The model error is defined as rˆcalm 1− ˆr meas m1 /ˆrmeasm 1 .
insect wing mass measured ˆrm2 measured ˆrm1 calculated ˆrm1 model error
species [mg] [-] [-] [-] (%)
hawkmoth 44.79 0.38 0.29 0.31 6.90
hoverfly 0.21 0.44 0.36 0.37 2.78
dronefly 0.53 0.40 0.32 0.33 3.13
bumblebee 0.54 0.44 0.36 0.37 2.78
The CWMD is also important since it determines the pitching inertia and, thus, the passive pitching motion. Due to the lack of knowledge on the CWMD of insect wings, we assume two different cross-sectional profiles, as shown in Fig. 2. The first is a kite profile, i.e., the thickness (or mass) of an arbitrary cross-section increases linearly from zero at both the leading and trailing edges to the maximal value tm at the pitching axis location. The second is a uniform thickness along the chord.
Based on the proposed SWMD and CWMD models, the mass distribution of an
arbitrary wing can be quantified. When the wing mass mw and the moment of inertia
are known, the unknown parameters λ1 and λ2 in Eq. 3 can be determined by solving
following system of equations
Z R 0 ρwγλ1eλ2rcdr− mw = 0, (4a) Z R 0 ρwγλ1eλ2rcr2dr− mwˆrm22 = 0, (4b)
where the value of γ depends on the CWMD (γ = 1/2 for the kite profile and γ = 1 for the uniform profile). It should be noted that the spanwise radii of the first and second moment of inertia (i.e., ˆrm1 and ˆrm2) do not change with the value of γ.
To validate the proposed mass distribution model, we use it to describe the wing mass distribution of four insects, including a hawkmoth, a hoverfly, a bumblebee and a dronefly. The wing mass and the dimensionless spanwise radii of the first and second moment of inertia (i.e., ˆrm1 and ˆrm2) of these insects have been measured in literature
(Ellington, 1984; Willmott and Ellington, 1997b). The values of λ1 and λ2 for each
insect wing are determined using Eqs. 4. Thus, the mass distribution over the entire wing is known. Then, the values of ˆrm1 are calculated and compared with the measured
values, as shown in Table 1. The comparison shows that the calculated ˆrm1 deviates
from the measured value by less than 6.9% for all four wings, which implies that the proposed model provides a good representation of the SWMD of insects wings.
2.3. Kinematics
Two frames are introduced to describe the flapping kinematics of a rigid wing, i.e., the inertial frame xiyizi and the co-rotating frame xcyczc, as shown in Fig. 3. The inertial frame is fixed at the wing joint. The xi axis coincides with the projection of the leading edge on the horizontal stroke plane while the wing is at rest at the middle stroke. The zi axis is perpendicular to the stroke plane, and the yi axis is defined by the
right-hand-rule. The co-rotating frame rotates with the wing, whose xc axis coincides with the
PA, and the yc and zc axes are perpendicular and parallel to the wing platform. The
sweeping angle φ and pitching angle η can be described by these two frames, as shown in Fig. 3.
A harmonic function
φ(t) = φmsin(2πf t) (5)
is used to prescribe the sweeping motion with an amplitude of φm and a drive frequency f . The passive pitching motion is determined by the wing inertia, the rotational stiffness of the linear torsional spring and the aerodynamic torque about the PA (i.e., the xcaxis in Fig. 3). Based on Euler’s second law of motion for a rigid body, the equation of motion for the pitching motion can be derived from
τxappliedc + τxinerc = 0, (6)
where τapplied
xc includes the elastic torque τ
elas
xc due to the resistance from deformed spring
and the aerodynamic torque τxaeroc that will be discussed in next section. The inertial torque τiner can be given in the co-rotating frame by
τiner=−Iαc− ωc× (Iωc) , (7)
where I is the inertia matrix, and ωc and αc are the wing angular velocity and
acceleration in the co-rotating frame, respectively. ωc and αc can be obtained by
ωc= [ ˙η, ˙φ sin η, ˙φ cos η]T, (8)
αc= [¨η, ˙φ ˙η cos η + ¨φ sin η, ¨φ cos η− ˙φ ˙η sin η]T. (9) The inertia matrix I can be calculated by
I = Ixcxc Ixcyc Ixczc Iycxc Iycyc Iyczc Izcxc Izcyc Izczc = R Mz 2 cdm 0 − R M xczcdm 0 RM(x2c+ zc2)dm 0 −R Mxczcdm 0 R Mx 2 cdm , (10)
where the cross term Ixczc results from the asymmetry of the wing mass distribution
about the PA. The inertia matrix expressed in the co-rotating frame does not change with the flapping motion, and this facilitates the easy derivation of the equation of motion. The inertial terms related to the wing thickness are ignored since the thickness is negligible when compared to the span and average chord length.
xc yc zc xi yi zi LE TE ϕ η
Figure 3: Two frames used for describing the motion of flapping wings. The inertial frame xiyizi is fixed to the wing joint, while the co-rotating frame xcyczcco-rotates with
the wing. φ and η represent the sweeping and pitching angles, respectively.
Substituting the xc component of τiner and the formula of the elastic torque (i.e., τelas
xc =−kηη) into Eq. 6, the equation of motion can be rewritten as
Ixcxcη + k¨ ηη = τ aero xc + 1 2Ixcxcφ˙ 2sin(2η) − Ixczcφ cos η,¨ (11)
where the last two inertial terms represent the centrifugal force and Euler force induced
by the sweeping motion, respectively. It can be seen that both the aerodynamic
and inertial terms introduce nonlinearity to the system. The influence of the inertial nonlinearity to the pitching motion depends on Ixcxc, Ixczc and φ(t). The values of Ixcxc
and Ixczc depend on the PA location, while φ(t) depends on the flapping frequency.
In this study, the influence of the PA location and flapping frequency on the pitching motion will be studied, and the optimal wing design will be investigated to maximize the energy-efficiency during hovering flight.
3. Aerodynamic and power consumption modelling 3.1. A quasi-steady aerodynamic model
The aerodynamic term τaero
xc in the equation of motion (Eq. 11) is calculated by a
predictive quasi-steady aerodynamic model as proposed by Wang et al. (2016) (see also Appendix A). The quasi-steady assumption means that the transient aerodynamic loads on a flapping wing are equivalent to that on the wing undergoing a steady motion at the same transient velocity and angle of attack. The viscous drag at the thin boundary layer and the suction load at the leading edge are ignored since they are negligible compared to the load resulting from the pressure difference between two sides of flapping wings (Sane, 2003). Consequently, the resultant aerodynamic force can be assumed to be
perpendicular to the wing platform, i.e., aligned to the yc axis in the co-rotating frame, over the entire stroke. Therefore, it is convenient to formulate the aerodynamic loads in the co-rotating frame. The resultant aerodynamic force and torque are calculated by Fyaeroc = Fytransc + Fyrotc + Fycouplc + Fyamc , (12a) τxaeroc = τxtransc + τxrotc + τxcouplc + τxamc , (12b) τzaeroc = τztransc + τzrotc + τzcouplc + τzamc , (12c)
where the four loading terms at the right-hand side of each equation originate from the wing translation, the rotation, the couping between them and the added mass effect. The translation-induced loading terms (i.e., Ftrans
yc , τ
trans
xc and τ
trans
zc ) represent the loads
when a wing undergoes pure translation at a certain angle of attack. The rotation-induced loading terms (i.e., Frot
yc , τ
rot xc and τ
rot
zc ) reflect the aerodynamic damping loads
on the wing that purely rotates about its PA. The coupling terms (i.e., Fycouplc , τxcouplc and τzcouplc ) represent the loads due to the coupling effect between the wing translation and rotation. The terms Fyamc , τxamc and τzamc reflect the loads due the added-mass effect as a result of the accelerated or decelerated fluid surrounding the wing. More details on the formulations of each term can be found in Appendix A. The lift L, which is in the direction of the zi axis, can be calculated by transforming Fyaeroc from the co-rotating
frame to the inertial frame.
The model is used in this paper for its advantages over other quasi-steady models in terms of two aspects. First, in addition to the aerodynamic forces, the model provides a good prediction for the aerodynamic torque about the pitching axis, which is essential for calculating the passive pitching motion of flapping wings. Second, the model does not rely on any empirical parameters, which provides more freedom to set the dimension of flapping wings.
3.2. Power consumption
Due to the limited on-board power supply, the flight endurance of existing FWMAVs is generally less than 20 minutes (Floreano and Wood, 2015). This issue motivates us to improve the power efficiency of flapping wing designs.
For flapping wings with passive pitching, the total power consumption consists
of three components: (1) aerodynamic power Paero used to overcome the aerodynamic
drag, (2) inertial power Pinerto accelerate the wing and surrounding fluid, and (3) elastic power Pelas due to the resistance from the elastic spring. These can be formulated as
Paero(t) =−τxaero c ωxc − τ aero zc ωzc, (13a) Piner(t) =−τxinerc ωxc − τ iner yc ωyc− τ iner zc ωzc, (13b) Pelas(t) = kηηωxc, (13c)
respectively. Due to aerodynamic drag, Paero > 0 during most parts of a flapping
the sweeping motion due to the wing rotational effect. For such cases, Paero < 0. The kinetic energy of the wing increases during the accelerating phase of each half-stroke and, thus, Piner > 0. During deceleration, the wing starts to lose its kinetic energy, which implies Piner < 0. The lost kinetic energy can be dissipated, used to compensate the energy consumed by the drag, or stored in the kinetic energy recovery drive system (KERS). The KERS can be in different forms, e.g., insect thoraxes (Dudley, 2002) and elastic structures of FWMAVs (Bolsman et al., 2009; Lau et al., 2014).
The uncertainty of the kinetic energy transformation complicates the modelling of the exact power consumed by flapping wings. Instead, two extreme power consumption scenarios are modeled. First, a KERS is used, thus, the cycle-averaged elastic and inertia power will be zero. Therefore, the total average power consumption ¯PKERStotal is equal to the average aerodynamic power, i.e.,
¯ PKERStotal = 1 T Z T Paerodt, (14)
where T is the period of a flapping cycle. Second, a drive system which can not recover kinetic energy is used. In this case, the kinetic energy and the elastic energy will be first used to compensate the energy consumed by drag. The extra energy will be dissipated. For convenience, this type of drive system will be referred to as non-KERS in this work. The total average power consumption ¯Pnon-KERStotal for this extreme case can be calculated by ¯ Pnon-KERStotal = 1 T Z T
Ξ Paero+ Piner+ Pelas dt, (15)
where Ξ(•) is an operator to set negative values (excess stored energy) to zero and keep positive values the same.
The power consumption by flapping wings should lie between ¯Ptotal
KERSand ¯Pnon-KERStotal . To facilitate the comparison of energy-efficiency between different wing designs, the average power is normalized by the mass that can be lifted by the corresponding wing design and denoted by adding a cap (e.g., ˆP¯total
KERS and ˆP¯non-KERStotal ).
4. Optimization model
To investigate the influence of the pitching axis (PA) location on the power efficiency, ˆ
dr and ˆdt are optimized. Meanwhile, the kinematics is also optimized by setting
the sweeping frequency f and the rotational stiffness kη as design variables. The
minimization of the cycle-averaged power consumption ˆP¯total
KERS and ˆP¯non-KERStotal are set as the optimization objective for KERS and non-KERS, respectively. The lift generated by optimized wings needs to be equal to the required lift (e.g., half of the body weight for two-winged FWMAVs).
It is interesting to compare the energy efficiency between wings with an optimal PA location and traditional wing designs which have a straight LE as the PA (see Fig. 1). For a fair comparison, the kinematics, i.e. the values of f and kη, are also optimized for traditional wings subjected to the same objective and constraint.
For a fair comparison of the power consumption, all the wings studied in this work
share the same span R, aspect ratio A, wing mass mw, spanwise radius of the second
moment of inertia ˆrm2, sweeping amplitude φm and the lift generation requirement ¯Lreq.
These parameters are prescribed by referring to an adult female hawkmoth (HMF2) studied by Willmott and Ellington (1997b), as listed in Table 2.
Table 2: Prescribed parameters of flapping wings. A is the aspect ratio of a single wing and defined as R2/S. ¯L
req is the lift generation required for a single wing.
parameters R A mw rˆm2 φm L¯req
units [×10−3m] [-] [×10−6Kg] [-] [degrees] [×10−3N]
HMF2 51.00 2.77 46.01 0.38 61.64 9.21
model wings 50.00 2.50 50.00 0.38 60.00 9.80
Due to the nonlinearity existing in the equation of motion for the passive pitching motion, both the lift and power are non-convex functions in the design domain expanded by four design variables (i.e., ˆdr, ˆdt, f and kη). Therefore, a stochastic global optimization method (Li and Au, 2010) is used to obtain the rough optima. Then, the rough solution is taken as the initial value for a gradient-based optimization
method (fmincon from MATLAB R) which uses the sequential quadratic programming
algorithm.
In this work, we study the optimal PA for eight wings resulting from the combination of four spanwise area distributions (SWADs) (see Fig. 1) and two chordwise
mass distributions (CWMDs) (see Fig. 2). For each wing, both the two extreme
power consumption cases are considered as the optimization objective. Therefore, 16 optimization cases are studied in total.
5. Results and analysis
5.1. Optimal pitching axis location
For different wings, the shapes with the optimal PA location are compared to traditional wings which use a straight LE as the PA in Fig. 4. To facilitate the comparison, all the optimal pitching axes coincide with the straight LE of corresponding traditional wings. It can be seen that the pitching axes of all the optimized wings are located behind the LE. These optimized wings are quite different from traditional wings but can be very close to insect wings, as shown by the comparison between optimized HM-wings with the real hawkmoth wing in Fig. 4 (b).
Optimized wings with different CWMDs and subjected to different power
consumption scenarios show close PA locations. However, there are recognizable
differences between optimized wing shapes. First, when the same CWMD is used, the PA of optimized wings for KERS is further behind the LE as compared to the optimal PA for non-KERS. The difference is more prominent at the wing root than at the wing
0 0.02 0.04 −0.04 −0.02 0 0.02 r [m] zc [m] (a) Rect-wings 0 0.02 0.04 hawkmoth wing r [m] (b) HM-wings 0 0.02 0.04 r [m] (c) QE-wings 0 0.02 0.04 r [m] (d) Tri-wings KERS, CWMD: kite KERS, CWMD: uniform non-KERS, CWMD: kite non-KERS, CWMD: uniform
Figure 4: Optimal pitching axis locations for different types of wings. Gray areas represent traditional wings with a straight leading edge (LE) as the PA. To facilitate the comparison, all the optimal PAs coincide with the straight LE of corresponding traditional wings.
tip. Second, the optimal location of PA is more sensitive to the objective function than CWMD. It should be noted that aforementioned observations hold for all the four types of wings.
In order to interpret the optimal PA location, we plot the corresponding power consumption, pitching amplitudes, flapping frequency and the rotational stiffness for all the cases. Here, the pitching amplitude is defined as the maximum pitching angle during the stroke. Inter-comparisons are used to help the analysis.
5.2. Analysis
Figure 5 compares the power consumed by different optimized wings. It is clearly seen that the wings with optimal PA location show a dramatical reduction of the power consumption as compared to traditional wings. For KERS, up to 33% of power can be saved by the wing with the optimized PA when compared to the traditional HM-wing. The minimum amount of power that can be saved by the wing with optimized PA location is still more than 21%, as shown by the QE-wings for non-KERS.
It can be also observed that, for traditional wings, the use of KERS does not reduce the power consumption significantly as compared to non-KERS. In contrast, when the PA is at the optimal location, the power consumption shows a considerable drop for KERS. For instance, more than 13% of energy has been saved via recovering kinetic energy by the HM-wing with the optimal PA. The influence of PA location on the usefulness of kinetic energy recovery capacity can be better understood from the power plot in Fig. 6. The study of power history of the HM-wings with a kite profile CWMD is taken as an example. It can be seen that the aerodynamic power before pitching reversals is negligible for the HM-wings with optimal PA. Thus, the negative power
due to the wing deceleration can not be completely offset by the aerodynamic power. Such a significant amount of unnecessary power consumption can be avoided by using KERS. In contrast, the traditional HM-wing still experiences a significant amount of aerodynamic power consumption before reversals which is almost equal to the negative inertial power. Thus, the remainder of the kinetic energy that can be stored is very limited.
From Fig. 5, we can also see that the Tri-wings show the lowest power efficiency. This is due to the fact that the Tri-wings have a smaller ˆrs1 as compared to other wings
(see Fig. 1), which implies that the area of the Tri-wings is closer to the root on average. As a consequence, the aerodynamic efficiency of Tri-wings is lower due to the lower average flapping velocity.
The pitching amplitudes ηm corresponding to different optimized wings are
compared in Figs. 7 (a) and (b). The plots show that the pitching amplitudes of wings with a straight LE are lower than wings with an optimal PA, and the pitching amplitudes of Tri-wings are lower than other wings. By comparing these results with the power plots in Fig. 5, it can be concluded that a higher pitching amplitude leads to lower
power consumption. The conclusion can be explained by the relation between the
aerodynamic power and the pitching amplitude. This is because, for KERS and non-KERS, the power consumption is equal to or dominated by the aerodynamic power, respectively. We know that the aerodynamic power is proportional to the average drag, while the average drag is inversely proportional to the average lift-to-drag ratio since all the wings need to generate the same average lift. Furthermore, the value of the
lift-Rect-wings HM-wings QE-wings Tri-wings 30 40 50 60 p ow er [W/Kg] (a) CWMD: kite
Rect-wings HM-wings QE-wings Tri-wings 30 40 50 60 p ow er [W/Kg] (b) CWMD: uniform optimized wing, KERS traditional wing, KERS optimized wing, non-KERS traditional wing, non-KERS
Figure 5: Comparison of power consumption of different wings for both KERS and non-KERS. The horizontal axis ticks Rect-wing, HM-wing, Rect-wing and Tri-wing represent the wings with the chordwise area distribution of the rectangular wing, the hawkmoth wing, the quarter-ellipsoidal wing and the triangular wing as shown in Fig.1, respectively.
0 0.25 0.5 0.75 1 −50 0 50 100 0.25 0.75 0
upstroke downstroke upstroke
t/T [-]
p
ow
er
[W/Kg]
(a) HM-wing with optimal PA
0 0.25 0.5 0.75 1 −50 0 50 100 0.25 0.75 0
upstroke downstroke upstroke
t/T [-]
p
ow
er
[W/Kg]
(b) HM-wing with straight LE ˆ
Paero Pˆiner Pˆelas Pˆtotal
KERS Pˆnon-KERStotal
Figure 6: Comparison of the power histories for the HM-wings with a kite profile CWMD. (a) The wing has a optimal PA location. (b) The straight LE is taken as the PA. Filled areas between the curves for ˆPtotal
KERS and ˆPnon-KERStotal indicate the kinetic
energy that can be potentially recovered by drive system.
to-drag ratio is inversely proportional the angle of attack based on the condition that the resultant aerodynamic force is perpendicular to the chord. We also know that the average angle of attack is inversely proportional to the pitching amplitude. Therefore, the aerodynamic power is inversely proportional to the pitching amplitude.
The lift generation is related to both the angle of attack and the flapping velocity. According the aerodynamic model proposed by Wang et al. (2016), the lift coefficient of a translational wing reaches its maximal value at the angle of attack (AOA) of 45◦. In practice, due to the wing rotational and added-mass effects, a flapping wing
reaches the highest lift coefficient at an AOA that slightly deviates from 45◦. The
corresponding AOA will be referred to as lift-maximizing AOA for convenience. The maximal cycle-averaged lift coefficient is attainable when the cycle-averaged absolute deviation of the AOA from the lift-maximizing AOA is minimized. The deviation can be partly reflected by the pitching amplitude. For example, if the pitching motion is
assumed to be harmonic and the lift-maximizing AOA is 45◦, the pitching motion has
the minimal deviation when the pitching amplitude is 55◦. Therefore, for the wing
with a pitching amplitude much higher than the lift-maximizing AOA, the average lift coefficient is inversely proportional to the pitching amplitude. We can see from Fig. 7 (a) to Fig. 7 (d) that, for the wing with a higher pitching amplitude (i.e., a smaller average lift coefficient), a higher flapping frequency is required to increase the flapping velocity such that the lift constraint can be satisfied. The frequencies with respect to the optimized HM-wings with a optimal PA location (see Fig. 7 (c)) are 27.04 Hz and 26.46 Hz for KERS and non-KERS, respectively. In contrast, the average flapping frequency of the reference hawkmoth (Willmott and Ellington, 1997a) is 25.40 Hz. The
Rect-wings HM-wings QE-wings Tri-wings 40 50 60 70 80 90 ηm [degrees] (a) CWMD: kite
Rect-wings HM-wings QE-wings Tri-wings 40 50 60 70 80 90 ηm [degrees] (b) CWMD: uniform optimized wing, KERS traditional wing, KERS optimized wing, non-KERS traditional wing, non-KERS
Rect-wings HM-wings QE-wings Tri-wings 24 26 28 30 optimal f [Hz] (c) CWMD: kite
Rect-wings HM-wings QE-wings Tri-wings 24 26 28 30 optimal f [Hz] (d) CWMD: uniform
Rect-wings HM-wings QE-wings Tri-wings 10−5 10−4 10−3 wing types optimal kη [Nm/Rad] (d) CWMD: kite
Rect-wings HM-wings QE-wings Tri-wings 10−5 10−4 10−3 wing types optimal kη [Nm/Rad] (f) CWMD: uniform
Figure 7: Comparison of optimal results for four types of wings with two different chordwise mass distribution (CWMD) for both KERS and non-KERS.
slightly overestimated frequency can be attributed to the smaller span and the higher required lift production of the model wing as compared the hawkmoth wing (see Table 2). When the PA moves from the LE to the optimal location, the wing pitching inertia
decreases. Meanwhile, the PA gets closer to the center of pressure (CP) which is
normally located between the 1/4 chord to the mid-chord (Zhao et al., 2010; Wang et al., 2016). Both the change of the pitching inertia and the distance between the CP and PA tends to reduce the external torque applied on the PA. Therefore, in order to maintain the expected pitching amplitude, the wing rotational stiffness kη needs to decrease simultaneously. This is why the optimal stiffness used by wings with a straight LE is much larger than that for the wings with a optimal PA location, as shown in Figs. 7 (e) and (f). This also provides an explanation to the wing shape and flexibility of insect wings. With a properly located PA, an insect wing can generate sufficient lift in the most efficient way. Meanwhile, the wing can be very flexible which implies that
the structure can be light-weight. The light-weight wing design is also preferred for flapping wings of FWMAVs.
−60 −30 0 30 60 −90 −60 −30 0 30 60 90 ϕ [degrees] η [degrees] (CWMD: kite) (a) Rect-wings −60 −30 0 30 60 ϕ [degrees] (b) HM-wings −60 −30 0 30 60 ϕ [degrees] (c) QE-wings −60 −30 0 30 60 ϕ [degrees] (d) Tri-wings optimized wing, KERS optimized wing, non-KERS traditional wing, KERS traditional wing, non-KERS −60 −30 0 30 60 −90 −60 −30 0 30 60 90 ϕ [degrees] η [degrees] (CWMD: uniform) (e) Rect-wings −60 −30 0 30 60 ϕ [degrees] (f) HM-wings −60 −30 0 30 60 ϕ [degrees] (g) QE-wings −60 −30 0 30 60 ϕ [degrees] (h) Tri-wings
Figure 8: Phase portrait of optimal kinematics. The cross-line with a positive/negative slope means that the corresponding pitching motion is ahead/behind of the sweeping motion.
To gain more insight into the pitching motion of wings with the optimal PA, Figure. 8 shows the phase portraits of their kinematics and the comparison of the optimal kinematics between optimized and traditional wings. It can be seen that the reversals of pitching motion with respect to the wing with a optimal PA always take place in advance to the sweeping reversals. In contrast, the pitching reversals of traditional wings are generally delayed compared to sweeping reversals. The advanced pitching reversal can enhance the lift generation during the reversal phases (Dickinson et al., 1999; Nakata and Liu, 2012), and it has also been observed on the hovering hawkmoth (Nakata and Liu, 2012). The PA location might be also one of the explanations to the advanced pitching reversals observed for insect flight.
In Figs 5, 7 and 8, we showed the optimal power consumption, design variables and kinematics for wings with both the kite and uniform profile CWMDs. It can be
observed that the performance difference between the wings with these two CWMDs is marginal. This implies that the optimal PA location is not sensitive to the CWMD. For convenience, the analysis presented in this section was primarily based on the wings with the kite profile CWMD.
5.3. Influence of lift constraints
0 0.01 0.02 0.03 0.04 0.05 −0.02 −0.01 0 0.01 pitching axis r [m] zc [m] (a) KERS 0 0.01 0.02 0.03 0.04 0.05 pitching axis r [m] (b) non-KERS ¯ Lreq=0.5 g ¯ Lreq=0.75 g ¯ Lreq=1.25 g ¯ Lreq=1.5 g
Figure 9: Comparison of wing shapes with optimal location of pitching axes for the HM-wing with the kite profile CWMD when subjected to different lift constraints ( ¯Lreq). The
gray areas correspond to the lift constraint of 1 g.
In order to investigate the dependence of the optimal PA location on the lift constraint, we further optimize the PA location for different lift productions, including 50%, 75%, 125% and 150% of the original value (1 g). Taking the HM-wing with a kite profile CWMD as an example, the optimal shapes for different lift constraints are quite close to the original optimal wing shape, as shown in Fig. 9.
Figure. 10 plots the optimal frequency and torsional stiffness as well as the corresponding power consumption and pitching amplitudes with respect to different lift constraints. The power plots in Fig. 10 (a) confirm that the wings with the optimal PA location are more energy-efficient than traditional wings regardless of the lift constraint. The optimal pitching amplitudes as shown in Fig. 10 (b) are approximately equal to the amplitude of the original optimal wing. In order to satisfy different lift constraints, the optimal frequency increases with the required lift production, as shown in Fig. 10 (c). Figure. 10 (d) shows that the optimal rotational stiffness for both the wings with the optimal PA location and traditional wings increases with the required lift production
approximately at the same rate. However, the absolute change for wings with the
optimal PA location is one order smaller than traditional wings.
Combining the results in Figs. 9 and 10, we can see that the wing with a properly located PA can generate different lift while maintaining a high power efficiency. This can be realized by changing the flapping frequency and a fine tuning of the wing stiffness. For insects, there are several ways to achieve minor stiffness changes, including the vein blood circulation (Hou et al., 2014), wing warping (Ristroph and Childress, 2014). For FWMAVs, there also exist many approaches to tune the wing stiffness,
0.5 0.75 1 1.25 1.5 20 30 40 50 60 lift production [×10−3 Kg] p ow er [W/Kg] (a) 0.5 0.75 1 1.25 1.5 50 60 70 80 90 lift production [×10−3 Kg] optimal ηm [degrees] (b) optimized wing, KERS traditional wing, KERS optimized wing, non-KERS traditional wing, non-KERS 0.5 0.75 1 1.25 1.5 15 20 25 30 35 lift production [×10−3 Kg] optimal f [Hz] (c) 0.5 0.75 1 1.25 1.5 10−5 10−4 10−3 lift production [×10−3 Kg] optimal kη [Nm/Rad] (d)
Figure 10: Comparison of optimal designs of the HM-wings with the kite profile CWMD for different lift constraints. The HM-wings are optimized for both KERS and non-KERS.
e.g., piezoelectric polymers, electrorheological fluids, and electrostatic softening (Peters et al., 2015; Peters, 2016).
6. Conclusions
This paper studied the optimal pitching axis (PA) location for different flapping wings to maximize the energy efficiency during hovering flight. We found that the optimal PA is located between the leading edge (LE) and the mid-chord line, which can result in shapes very similar to insect wings. In contrast, traditional wings used by most flapping wing micro air vehicles (FWMAVs) simply pitch about their straight LEs. The comparison of power consumption shows that wings with optimal PA can save up to 33% of power as compared to traditional wings with optimized kinematics. The PA location also influences the usefulness of the kinetic energy recovery capacity of a drive system. More than 13% of the power consumption can be saved via recovering kinetic energy for wings with an optimal PA while the saved power is negligible for traditional wings. Furthermore, the wing with a properly located PA can generate different lift levels while maintaining a high energy efficiency. This can be realized by simultaneously changing
the flapping frequency and fine tuning of the wing stiffness.
It is nontrivial to give a general optimal PA location for flapping wings considering the diversity of flight conditions and the nonlinear relation between the PA location and the flight performance. However, for flapping wings pitching passively, the PA location should be carefully designed.
Acknowledgments
This work is financially supported by China Scholarship Council (CSC NO. 201206290060).
Appendix A. Quasi-steady aerodynamic model
As a consequence of the wing flapping motion, the translational velocity varies from the wing root to tip and from the LE to TE. Hence, the blade element method (Osborne, 1951) is applied for discretizing the wing planform in both chordwise and spanwise directions. The resultant aerodynamic load can be formulated as a superposition of four loading terms: translation-induced load, rotation-induced load, coupling load between wing translation and rotation as well as added-mass load. More specifically, the resultant translation-induced force in ycdirection and corresponding torques about the xc and zc axes can be formulated as
Fytransc =− sgn(ωzc) 1 2ρf(ω 2 yc+ ω 2 zc)C trans Fyc Z R 0 x2ccdxc, (A.1) τxtransc = −sgn (ωzc) 1 2ρf ω 2 yc+ ω 2 zc C trans Fyc ˆ zcptrans− ˆd RR 0 x 2 cc2dxc, ωyc ≤ 0 −sgn (ωzc) 1 2ρf ω 2 yc+ ω 2 zc C trans Fyc 1− ˆztrans cp − ˆd RR 0 x 2 cc2dxc, ωyc > 0 (A.2) and τztransc =−sgn (ωzc) 1 2ρf ω 2 yc + ω 2 zc C trans Fyc Z R 0 x3ccdxc, (A.3)
where sgn(·) is the signum function, ρf is the fluid density. CFtransyc is the translational force
coefficient which is formulated as a function of the angle of attack ( ˜α) and wing aspect ratio (A) (Taha et al., 2014), i.e., 2πA sin (˜α) 2 +√A2+ 4−1
. ˆztrans
cp is the center
of pressure due to the wing translation. It is calculated with the formula (2π)−1|˜α|, which is an analytical approximation of data from experimental measurements (Wang et al., 2016). The negative and positive values of ωyc mean that the translational velocity
component vzc(= −xcωyc) points at the LE and TE, respectively. Generally, ωyc changes
signs during wing reversal phases. For instance, ωyc becomes positive before the sweeping
reversal if the pitching reversal happens earlier.
as Fyrotc =1 2ρfωxc|ωxc|C rot D Z R 0 Z dcˆ ˆ dc−c zc|zc| dzcdxc, (A.4) τxrotc =−1 2ρfωxc|ωxc|C rot D Z R 0 Z dcˆ ˆ dc−c |zc|3dzcdxc, (A.5) and τzrot c = 1 2ρfωxc|ωxc|C rot D Z R 0 Z dcˆ ˆ dc−c zc|zc| xcdzcdxc, (A.6) where Crot
D is the rotational damping coefficient. The value of CDtrans when the AOA is equal to 90◦ is adopted as the rotational damping coefficient.
The resultant coupling force between the wing translation and rotation and corresponding torques are determined by
Fycoup c = πρfωxcωyc h RR 0 3 4 − ˆd c2x cdxc+ RR 0 1 4c 2x cdxc i , ωyc ≤ 0 πρfωxcωyc h RR 0 ˆd− 1 4 c2x cdxc+ RR 0 1 4c 2x cdxc i , ωyc > 0. (A.7) τxcoupc = πρfωxcωyc h RR 0 3 4 − ˆd 1 4 − ˆd c3xcdxc+ RR 0 1 4 3 4 − ˆd c3xcdxc i , ωyc ≤ 0 πρfωxcωyc h RR 0 3 4 − ˆd ˆd− 1 4 c3x cdxc+ RR 0 1 4 ˆd− 1 4 c3x cdxc i , ωyc > 0 (A.8) and τzcoupc = πρfωxcωyc h RR 0 3 4 − ˆd c2x2 cdxc+ RR 0 1 4c 2x2 cdxc i , ωyc ≤ 0 πρfωxcωyc h RR 0 ˆd− 1 4 c2x2 cdxc+ RR 0 1 4c 2x2 cdxc i , ωyc > 0. (A.9)
Finally, the resultant force and torques due to the added-mass effect are formulated as Fyamc =− Z R 0 m22aycdxc− αxc Z R 0 m24dxc (A.10) τxamc =− Z R 0 m42axcdxc− αxc Z R 0 m44dxc, (A.11) and τzam c =− Z R 0 m22aycxcdxc− αxc Z R 0 m24xcdxc. (A.12)
Here we denote the directions of translational motions along axes ycand zcof a wing strip as the “2” and “3” directions and the rotation about the xc axis as the “4” direction. The parameter mij is used to represent the load induced by the added-mass effect in the i direction due to a unit acceleration in the j direction. axc and ayc are the translational
With the resultant aerodynamic force in the co-rotating frame, the lift L and drag D, which are defined as the resultant force components with respect to zi and yi axes in the inertial frame, can be calculated by
L = (Fytransc + Fyrotc + Fycouplc + Fyamc ) sin η (A.13)
and
D = (Fytransc + Fyrotc + Fycouplc + Fyamc ) cos η, (A.14) respectively. For the harmonic sweeping motion as defined in Eq. 5, the average lift force ¯L is perpendicular to the stroke plan (i.e., xiyi plane), and the average drag force
¯
References
S. A. Ansari. A nonlinear, unsteady, aerodynamic model for insect-like flapping wings in the hover with micro air vehicle applications. Phd thesis, Cranfield University, 2004. S. A. Ansari, K. Knowles, and R. Zbikowski. Insectlike Flapping Wings in the Hover
Part II: Effect of Wing Geometry. Journal of Aircraft, 45(6):1976–1990, nov 2008. A. J. Bergou, S. Xu, and Z. J. Wang. Passive wing pitch reversal in insect flight. Journal
of Fluid Mechanics, 591:321–337, oct 2007.
G. J. Berman and Z. J. Wang. Energy-minimizing kinematics in hovering insect flight. Journal of Fluid Mechanics, 582:153–168, jun 2007.
C. R. Betts and R. J. Wootton. Wing shape and flight behaviour in butterflies
(Lepidoptera: Papilionoidea and Hesperioidea): a preliminary analysis. Journal of Experimental Biology, 288(1):271–288, 1988.
C. T. Bolsman, J. F. L. Goosen, and F. van Keulen. Design overview of a resonant wing actuation mechanism for application in flapping wing MAVs. International Journal of Micro Air Vehicles, 1(4):263–272, dec 2009.
A. Chaudhuri, R. Haftka, P. G. Ifju, D. Villanueva, K. Chang, and J. Rue. Experimental Optimization and Uncertainty Quantification of Flapping Wing of a Micro Air Vehicle. In 10th World Congress on Structural and Multidisciplinary Optimization, pages 1–10, 2013.
S. A. Combes and T. L. Daniel. Flexural stiffness in insect wings. II. Spatial distribution and dynamic wing bending. The Journal of experimental biology, 206:2989–2997, sep 2003.
G. C. H. E. de Croon, K. M. E. de Clercq, R. Ruijsink, B. Remes, and C. de Wagter. Design, aerodynamics, and vision-based control of the DelFly. International Journal of Micro Air Vehicles, 1(2):71–98, 2009.
M. H. Dickinson, F. O. Lehmann, and S. P. Sane. Wing rotation and the aerodynamic basis of insect flight. Science, 284(5422):1954–1960, jun 1999.
R. Dudley. The biomechanics of insect flight: form, function, evolution. Princeton University Press, 2002.
R. Dudley and C. P. Ellington. Mechanics of Forward Flight in Bumblebees: I.
Kinematics and Morphology. J. Exp. Biol., 148(1):19–52, jan 1990.
C. P. Ellington. The Aerodynamics of Hovering Insect Flight. II. Morphological
Parameters. Philosophical Transactions of the Royal Society B: Biological Sciences, 305(1122):17–40, feb 1984.
A. R. Ennos. Inertial and aerodynamic torques on the wings of diptera in flight. J. Exp. Biol., 142(1):87–95, mar 1989.
D. Floreano and R. J. Wood. Science, technology and the future of small autonomous drones. Nature, 521(7553):460–466, 2015.
D. Hou, Y. Yin, H. Zhao, and Z. Zhong. Effects of blood in veins of dragonfly wing on the vibration characteristics. Computers in biology and medicine, 58C:14–19, dec 2014.
D. Ishihara, T. Horie, and M. Denda. A two-dimensional computational study on
the fluidstructure interaction cause of wing pitch changes in dipteran flapping flight. Journal of Experimental Biology, 212:1–10, 2009.
M. Keennon, K. Klingebiel, H. Won, and A. Andriukov. Development of the nano hummingbird: A tailless flapping wing micro air vehicle. AIAA, (January):1–24, 2012.
G. K. Lau, Y. W. Chin, J. T. W. Goh, and R. J. Wood. Dipteran-Insect-Inspired Thoracic Mechanism With Nonlinear Stiffness to Save Inertial Power of Flapping-Wing Flight. IEEE TRANSACTIONS ON ROBOTICS, pages 1–11, 2014.
H.-S. Li and S.-K. Au. Design optimization using Subset Simulation algorithm.
Structural Safety, 32(6):384–392, nov 2010.
K. Y. Ma, P. Chirarattananon, S. B. Fuller, and R. J. Wood. Controlled flight of a biologically inspired, insect-scale robot. Science, 340(6132):603–607, may 2013. T. Nakata and H. Liu. Aerodynamic performance of a hovering hawkmoth with flexible
wings: a computational approach. Proceedings. Biological sciences / The Royal
Society, 279(1729):722–31, mar 2012.
Q.-V. Nguyen, W. L. Chan, and M. Debiasi. An insect-inspired flapping wing
micro air vehicle with double wing clap-fling effects and capability of sustained hovering. In A. Lakhtakia, M. Knez, and R. J. Mart´ın-Palma, editors, Bioinspiration, Biomimetics, and Bioreplication 2015, pages 1–11, mar 2015.
M. F. M. Osborne. Aerodynamics of flapping flight with application to insects. Journal of Experimental Biology, 28:221–245, 1951.
H. J. Peters. A Controllability Approach for Resonant Compliant Systems: Applied to a Flapping Wing Micro Air Vehicle. Doctoral thesis, Delft University of Technology, mar 2016.
H. J. Peters, Q. Wang, J. F. L. Goosen, and F. Van Keulen. Active control of the passive pitching of a flapping wing with electrostatic clamping. In 7th ECCOMAS Thematic Conference on Smart Structures and Materials (SMART 2015), pages 1–19, Azores, Portugal, 2015.
L. Ristroph and S. Childress. Stable hovering of a jellyfish-like flying machine. Journal of the Royal Society, Interface / the Royal Society, 11(20130992), 2014.
S. P. Sane. The aerodynamics of insect flight. Journal of Experimental Biology, 206 (23):4191–4208, dec 2003.
P. Seshadri, M. Benedict, and I. Chopra. Understanding Micro Air Vehicle Flapping-Wing Aerodynamics Using Force and Flowfield Measurements. Journal of Aircraft, 50(4):1070–1087, jul 2013.
A. Shahzad, F.-B. Tian, J. Young, and J. C. S. Lai. Effects of wing shape, aspect ratio and deviation angle on aerodynamic performance of flapping wings in hover. Physics of Fluids, 28(11):111901, nov 2016.
E. Stewart, M. Patil, R. Canfield, and R. Snyder. Shape Optimization of a Flapping Wing. In 54th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, pages 1–11, Boston, Massachusetts, 2013.
H. E. Taha, M. R. Hajj, and P. S. Beran. State-space representation of the unsteady aerodynamics of flapping flight. Aerospace Science and Technology, 34(1):1–11, apr 2014.
Q. Wang, J. F. L. Goosen, and F. van Keulen. Optimal hovering kinematics with respect to various flapping-wing shapes. In International Micro Air Vehicle Conference and Flight Competition (IMAV2013), Toulouse, France, 2013.
Q. Wang, J. F. L. Goosen, and F. van Keulen. Study of design parameters of
flapping-wings. In International Micro Air Vehicle Conference and Flight Competition (IMAV2014), Delft, The Netherlands, 2014.
Q. Wang, J. F. L. Goosen, and F. van Keulen. A predictive quasi-steady model of aerodynamic loads on flapping wings. Journal of Fluid Mechanics, 800:688–719, aug 2016.
J. P. Whitney and R. J. Wood. Aeromechanics of passive rotation in flapping flight. Journal of Fluid Mechanics, 660:197–220, jul 2010.
A. P. Willmott and C. P. Ellington. The mechanics of flight in the hawkmoth Manduca sexta. I. Kinematics of hovering and forward flight. The Journal of experimental biology, 200(21):2705–2722, nov 1997a.
A. P. Willmott and C. P. Ellington. The mechanics of flight in the hawkmoth Manduca sexta. II. Aerodynamic consequences of kinematic and morphological variation. Journal of Experimental Biology, 200(21):2723–2745, nov 1997b.
R. Wood. The First Takeoff of a Biologically Inspired At-Scale Robotic Insect. IEEE Transactions on Robotics, 24(2):341–347, apr 2008.
L. Zhao, Q. Huang, X. Deng, and S. P. Sane. Aerodynamic effects of flexibility in flapping wings. Journal of the Royal Society, Interface / the Royal Society, 7(44): 485–97, mar 2010.
