AIAA 2015-3155

AIAA Aviation
22-26 June 2015, Dallas, TX
33rd AIAA Applied Aerodynamics Conference

Numerical Simulation of a Flapping Wing in Ansys-Fluent

Based on the Fruit Fly Wing Kinematics

Brayant Nieto*, **, Xiolesmy Alfonso† and Karen M. Basto*.

Universidad de San Buenaventura, Bogota D.C., Colombia, 110911

Jaime A. Escobar‡.
Universidad de San Buenaventura, Bogota D.C., Colombia, 110911

**
Author for correspondence (e-mail: bnieto@[Link])
Downloaded by CARLETON UNIVERSITY LIBRARY on July 17, 2015 | [Link] | DOI: 10.2514/6.2015-3155

This document describes the aerodynamic analysis of fruit fly ‘Drosophila Melanogaster’ flapping wings
for a hovering configuration and symmetric movement of trajectories including phase difference.
Computational Fluid Dynamics CFD is applied to understand the kinematics behavior and the flapping
phenomenon. The simulation was performed in 2D employing the ANSYS® FLUENT, Academic
Research software, replicating, in standard conditions, an elliptical geometry of 3% thickness in an O-type
computational grid; a laminar, incompressible and unstable fluid was used, with a Reynolds number of
151, evaluating cases with angles of attack of 25°, 35° and 45°. Once results were obtained, the physical
effects present in this type of flight such as the Kramer and Wagner effects, Karman and Leading Edge
vortices, and wake capture, were identified using detailed vortex contours throughout the flight cycles
studied, giving insight of lift and drag coefficients tendency as well as the increase in vortex generation as
the angle of attack is augmented, followed by an analysis and comparison with experimental and
simulation data obtained from previous researches in order to validate the methodology applied in this
study.

Nomenclature
𝐴𝑜𝐴 Angle of Attack
𝐿 Plate longitude
𝑅𝑒 Reynolds number
𝛼 Rotational movement
𝛼𝑎 Movement amplitude
𝛼0 Phase difference between velocity [x, y] and omega [z]
𝛽 Stroke plane angle
𝐶𝑙 Lift coefficient
𝐶𝑑 Drag coefficient
𝑓 Flapping frequency
ℎ Translational movement
ℎ𝑎 Movement amplitude in a half stroke
Θ Translational velocity
∅ Phase difference between velocity [x, y] and omega [z]
𝑥 Horizontal axis
𝑦 Vertical axis

*
Undergraduate Student, Department of Aeronautical Engineering, kbasto@[Link],
†
Undergraduate Student, Department of Aeronautical Engineering, malfonso@[Link], AIAA
Student Member.
‡
Associate Professor, Department of Aeronautical Engineering, jescobar@[Link], AIAA Member.

Copyright © 2015 by Brayant Nieto; Xiolesmy Alfonso; Karen Basto; Jaime Escobar. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.
 𝑦+ Dimensionless wall distance
Ω Rotational velocity
∆𝑡 Time increment
𝑡 Time
I. Introduction

T he complex flapping mechanism, surface roughness and design of the wings, apart from other factors, are the
main cause that insects have such an efficient, stable flight with great maneuverability compared with other flying
animals; reason why researchers around the world such as Michael Dickinson, J. Birch, J. Wang1 S. N. Fry, R.
Sayaman2, W. Shyy, Y. Lian, J. Tang, D. Viieru, and H. Liu3 have studied, in the last decades, insect´s flapping flight
aerodynamics. The technological and social contribution of this study in the development of a M.A.V. prototype
depends of a deep understanding of the kinematics and aerodynamic forces associated, in this case, with the fruit fly´s
‘Drosophila Melanogaster’ flapping flight.
With today’s tendency of Computational Fluid Dynamics (CFD) simulations, which allow us to recreate accurately
Downloaded by CARLETON UNIVERSITY LIBRARY on July 17, 2015 | [Link] | DOI: 10.2514/6.2015-3155

the physical conditions of fluids, the goal of this study is to simulate in ANSYS® FLUENT, Academic Research,
Release 15.0 the kinematics of a fruit fly´s wings, in order to predict its aerodynamic characteristics and provide
information adaptable to future M.A.V. prototype designs as well as an appropriate implementation of the
methodology in order to retrieve accurate and valid results when CFD is being applied. Throughout the development
of this research, various similar studies are referenced, taking into consideration fundamental parameters, to establish
the CFD simulation and the necessary conditions needed to obtain the results presented. We also explain in a simplified
manner the basic characteristics of the wing movement and the physical theories implied in it.

II. Fruit Fly Flapping Kinematics

The fruit fly depends on different trajectories of wing flapping in order to produce the aerodynamic forces necessary
to maintain its position in all of its flight configurations. Fig. 1 shows an example of trajectories for the hovering
flight, which were used for this parametric study.

a b
)

c d
Figure 1. Fruit fly flapping trajectories. a) U-shaped trajectory of Drosophila in free flight, b) 8-shaped trajectory
without correction at the end of each half stroke, c) Linear trajectory showing the amplitude of the stroke without
phase difference, d) Elliptical trajectory of the wing 15.

A. Magnus Effect
The Magnus effect produces a lift force, perpendicular to that of the fluid current, in a body of smoothened surface
that rotates and moves at the same time. This is due to the induced velocity and pressure differences between both the
upper and lower surfaces of an airfoil, generated by the vortices that flow with and against the stream flow trajectory
evidenced by Jost Seifert4. In insect flapping, the wing´s direction of rotation influences the total force orientation,
using the Magnus effect, along with the wake capture, to produce a maximum lift peak during the rapid pitch-up
maneuver at the rotations of the wing as described by W. Shyy et al.3.

B. Leading Edge Vortex – Delayed Stall

Sanjay Sane5 mentions a lift enhancement mechanism known as the Leading Edge Vortex (LEV), which produces
a delayed stall. This is due to the rapid acceleration that the airfoil experiences at the beginning of each half stroke.
The LEV remains attached to the wing during most of the half stroke and is separated at the end when the rotational
phase begins, it is produced and fed by the detachment of the fluid from the leading edge and it moves throughout the

2
 wing span starting at the root of the wing, working its way to the tip where it meets the wing tip vortex TiV, becoming
a shedding wake at the end, as has been previously suggested by M. H. Dickinson et al. 6.

C. Kramer Effect
In order to analyze insect flying maneuvers, at the end
of each upstroke and down stroke, the wings perform
pronation and supination respectively, taking its span as
the axis of rotation, which allows the wing to maintain a
positive angle of attack and generate lift as it moves as it
is shown in Fig. 27.
According to Dickinson et al.6, it is due to the
rotational circulation that the total forces in a wing act
normal to its chord, not perpendicular to the direction of
Downloaded by CARLETON UNIVERSITY LIBRARY on July 17, 2015 | [Link] | DOI: 10.2514/6.2015-3155

movement; likewise, the viscous forces of the fluid make

it flow smoothly along the trailing edge. This effect is
known as the Kutta condition, and fixes a stagnation point
for the fluid in the trailing edge and as consequence, the
circulation quantity and the force produced depends
mainly on the position of the rotation axis. S. Sane and M. Figure 2. Kramer Effect. Schematic representation of
Dickinson8 detect that during the rotation of the wing, delayed stall (2-1) and rotational lift (2-2). Busra Akay7.
additional flow circulation is created, which produces
rotational forces that increase or decrease the total
translation force, effect known as the Kramer effect.

D. Wagner Effect
Insect flight depends on repetitive cycles
made of translation and rotational phases,
simplified as start, end and reverse which are
equal to one complete flapping cycle. Ansari,
Bikowski and Knowles 9 consider each rotation
at the start and/or end of a half stroke is
composed by three stages: deceleration, reverse
and acceleration as shown by parts a, b and c of
Fig. 3. This behavior ‘start-end-reverse’ is a key
feature for the aerodynamics enabling insect
flight; at the start of the movement, the wing
c sheds a starting vortex that flows against the
direction of the fluid present around the airfoil,
augmenting the lift produced as it is shown in
figure 3, which is known as the Wagner effect.
Figure 3. Wagner effect and Wake Capture. Vortex in grey When the wing decelerates at the end of each
rotate clockwise. (a) As the wing decreases its speed at the end of half stroke, a stopping vortex is shed, which
the half cycle, a stopping vortex is shed. (b) During the rotation, a doesn`t allow the lift produced to be reduced.
vortex of opposite circulation is shed, result of the pitch up When the wing rotates and starts the following
performed. (c) At the start of the second half stroke, a starting half stroke at high AoA, a new starting vortex is
vortex is shed rotating in the same direction as the previous produced above the leading edge that interact
stopping vortex generated. Ansari, Bikowski, and Knowles9. with the stopping vortex shed during the
preceding half cycle, thus increasing the lift and
the Wagner effect, as was documented by W.
Shyy, Y. Lian, J. Tang, D. Viieru, and H. Liu 3.

E. Wake Capture

Based on experimentation, F. O. Lehman10 saw that the energy, in form of momentum and heat, released from the
flapping motion is regained and reused due to the interaction between the wing and the wakes created; this process is

3
 known as wake capture, which is applied to
increase the aerodynamic performance. Sane5
emphasizes that this effect is observed mainly
during hovering, making it necessary to perform
various similar flapping cycles in order to achieve
a periodic state.
This phenomenon is produced immediately
with the wing`s change of direction during the
pronation and supination at the start and end of
each half cycle9. Fig. 4 shows a deceleration of the
wing (A, B) in order to change direction and
perform the rotation, process in which the LEV is
shed and a Trailing Edge Vortex TEV is created.
These vortices of opposite rotation create a flow
Downloaded by CARLETON UNIVERSITY LIBRARY on July 17, 2015 | [Link] | DOI: 10.2514/6.2015-3155

current between them (C, D) which increase the

relative velocity going into the wing`s airfoil (E,
F), generating some lift and drag maximum peaks Figure 4. Wing-Wake Interactions. Parts A-F show a wing
as Ansari, Bikowski, Knowles 9, Sanjay Sane5, D. section during a flapping cycle5.
Lentink and M. H. Dickinson11 noted.

III. Methodology and Processing

Taking into consideration the study performed by Dilek F. Kurtulus 12, who indicates that the best results are of the
1% thickness elliptical airfoil, when compared to 9% and 12 % thickness using a Digital Particle Image Velocimetry
(DPIV), the present study considers an elliptical airfoil geometry of 3% thickness and a chord equal to 1 mm, based
also on the analysis performed by Sebastian Leguizamon 13, noting that a 1% thickness airfoil would be too thin for
practical applications and would affect the LEV behavior along with the flow structure.
Angles of attack of 25°, 35°, and 45° were simulated in order to coincide with the regions in which two-dimensional
and three-dimensional time averaged aerodynamic forces results are comparable according to W. Shyy et al.14. The
stroke amplitude chosen was 5 times the chord due to the affirmations from Dickinson et al.1 that results obtained by
changing the translation of the airfoil from 2 to 5 times the chord do not vary significantly. In order to identify the
correct time step size, preliminary simulations were carried out using sizes of 1e-05 s and 1e-06 s, showing similar
results, thus deciding on a 1e-05 s time step analyzed by 20 iterations, to perform the simulation process, reaching a
total of 5004 time steps per case to replicate 12 flapping cycles with accurate results in the flight configuration
discussed. A convergence criteria was applied evaluating the tendency of the aerodynamic coefficients all through the
12 cycles, taking critical points in each of the cycles, valuing the 30%, 60% and 80% of the stroke in the lift coefficient
representations, and the 5%, 35% and 90% of the stroke in the drag coefficient curves, which generated a trend
condition to a unique value allowing the condition of convergence.
Comparable results are to be obtained by varying parameters such as trajectory and rotational timing, in order to
examine comparable aerodynamic forces and the structures of vortices created.

A. Flapping Wing Kinematics

For this study, the flapping cycle is divided into a
translation and rotational phase. The translation consists
of two half strokes: down-up strokes. At the end of each
of these phases starts the rotation of the wing, which is
divided into supination and pronation respectively.
In order to define the wing kinematics adequately, a
UDF was used, programmed in C language including
the sinusoidal equations of velocity with the same
translation and rotation frequency, described in all the
different trajectories, establishing the center of rotation
at 50% of the chord and having a phase difference of
rotation timing depending of the value given to 𝛼0 and
∅.
Figure 5. Normal hovering representation in
symmetric flapping flight 3.

4
 Positioning Equations

ℎ [𝑥, 𝑦] = (ℎ𝑎 )𝑠𝑒𝑛(2𝜋𝑓𝑡 + 𝛼0 ) (1)

𝛼[𝑧] = 𝛼𝑎 𝑠𝑒𝑛(2𝜋𝑓𝑡 + ∅) (2)

Velocity Equations

Θ [𝑥, 𝑦] = (ℎ𝑎 )(2𝜋𝑓)𝑐𝑜𝑠(2𝜋𝑓𝑡 + 𝛼0 ) (3)

Ω[𝑧] = 𝛼𝑎 (2𝜋𝑓)𝑐𝑜𝑠(2𝜋𝑓𝑡 + ∅) (4)

where ℎ𝑎 is given a value of 2.5 mm, setting the flapping frequency to 240 Hz and the stroke amplitude to 5 times the
chord.
Downloaded by CARLETON UNIVERSITY LIBRARY on July 17, 2015 | [Link] | DOI: 10.2514/6.2015-3155

B. Meshing Theory and Computational Configuration.

A dynamic, hybrid, O-type mesh was designed using ANSYS® ICEM, Academic Research, Release 15.0,
composed by an internal section of quadrilaterals with a radius of 4 times the chord. Surrounding this internal section,
a tetrahedral mesh was made, extending up to a boundary limit or radius equivalent to 20 times the chord. The
computational domain is re-meshed and smoothened as the airfoil recreates the cycles throughout the flow of time,
making use of a finer mesh close to its surface with a ‘y’ of 0.001 m based on a Y+ of 0.01 in order to obtain a better
definition of the vortices shed. A concentric mesh was the result of this design process, starting with a down stroke
and having a total of 72430 cells.

Figure 6. Mesh created. a) Hybrid mesh showing connection between inner and outer mesh. b) Close-up of airfoil
surrounded by quadrilateral cell mesh.
The software ANSYS® FLUENT, Academic Research, Release 15.0 was used to perform the simulations,
applying the Laminar model, coupled pressure-based solver, and pressure-velocity SIMPLE C method, ad transient
flow. In order to perform the re-meshing, the dynamic mesh was deformed with a skewness factor of 0.8.

IV. Results
When analyzing the two-dimensional air flow around the wing section of the “Drosophila Melanogaster” and
comparing the results with experimental data, graphs were created in order to observe the behavior of the lift and drag
coefficients according to the three angles of attack previously established, and instantaneous snapshots of certain time
sequences during the 12th flapping cycle were gathered in order to perform a qualitative study of vortex generation
and interaction.

A. Analysis of Aerodynamic Forces

Three different cases were studied varying the angle of attack during the translational phase. The aerodynamic
force coefficients results for cases 1, 2 and 3 that were analyzed with AoA of 25°, 35°, and 45° respectively are
presented in Fig. 7, which yielded aerodynamic coefficients of different magnitude, with case 3 having a wider range
and higher set of values. Point 1 in the graphs of this figure represents the moment in which the Wagner effect and a

5
 pitch-up rotation performed at the end of each pronation and supination, increase the value of lift to a first peak and
the drag to its maximum value during the upstroke. Followed by an automatic decrease in the coefficient values shown
by point 2, produced by the acceleration of the fluid towards the trailing edge, decreasing simultaneously the high-
pressure stagnation area on the lower surface of the airfoil, as observed by Jian Tang, D Viieru, and Wei Shyy15.
The three cases present a discrepancy between them regarding the time in which coefficient values increase and
decrease, as seen especially in point three of Fig. 7, in which there is a drag increment and the lift reaches its maximum
value due to wake capture of vortices shed during the previous half stroke as part of the translational phase. After that,
due to the speed reduction and increase of AoA during the beginning of the rotational phase, there is an abrupt fall of
forces down to point 4, where the minimum values are reached.
With the analysis of Fig. 7, it can be observed that the highest and lowest values of lift coefficients are present
during the upstroke, whereas the highest and lowest values of drag coefficients occur in the downstroke. However, in
case 1 the values are the same throughout both of the half cycles, indicating that the aerodynamic forces and their
behavior all through the stroke depend considerately on the angle of attack.
Downloaded by CARLETON UNIVERSITY LIBRARY on July 17, 2015 | [Link] | DOI: 10.2514/6.2015-3155

Figure 7. Lift and Drag Coefficients comparison for different angles of attack. Coefficient values during 12th
flapping cycle.

B. Analysis of flow field and vorticity contours.

The 12th flapping cycle, starting from rest, is analyzed (Fig. 8) in a similar manner in which Jian Tang, D Viieru,
and Wei Shyy15 evaluate their results, observing, first of all, an airfoil acceleration and reduction of the AoA during
the culmination of the pronation at the beginning of the downstroke t/T=4.0; this rotation accelerates the flow
surrounding the leading edge and trailing edge, due to vortices A and B present in t/T=4.1 for 𝛼=45°, creating a
suction zone (zone 1 at t/T=4.1 for 𝛼=35°) on the upper surface of the wing section. The high-pressure stagnation area
in the lower surface of the airfoil (zone 2 at t/T=4.1 for 𝛼=35°) is augmented due to the fluid drawn in by the vortices
formed previously (t/T=4.1, vortices C and D for 𝛼=45°). This combination of high and low pressure areas produces
a higher lift prior to starting the translational phase of the downstroke as shown by point 1 in figure 7. As the airfoil
rotates decreasing its AoA, the flow speed is raised near to the trailing edge, diminishing the high pressure stagnation
area at the lower surface of the airfoil at t/T= 4.2, reducing thus the lift reaching a minimum value pictured by point 2
in figure 7.
About halfway of both down and upstroke, the airfoil travels at a constant AoA and velocity creating, on the upper
surface of the airfoil (zone 3 at t/T= 4.3 for 𝛼=45°), a vortex that is dragged along with the wakes captured along the
translation in order to raise the lift and drag to its highest value as it is shown by point 3 of Fig. 7. After reaching the
AoA adequate for a constant translational velocity, the wing section reduces its speed and increases the AoA due to a
pitch up motion, resulting in flow separation from the upper airfoil surface (t/T 4.4 and 4.5 for 𝛼=45°, zone 4). This
shedding of vortices, combined with the fast translational deceleration, reduce the vortex flow and the drag coefficient
reaches its lowest peak at the start and end of each half stroke, as seen in point 4 of figure 9. This analysis agrees with
the annotations made by W. Shyy, Y. Lian, J. Tang, D. Viieru, and H. Liu3 who describe with precision the phenomena
present throughout the flapping cycle.
The description of the vorticity contours is performed only from t/T 4.0 to 4.5 corresponding to a downstroke
section in which the moments prior to it (t/T 4.6 to 4.9) present the same aerodynamic effects with opposite vortex
flow circulation, giving a mirror illusion when comparing A-C versus D-F in Fig. 8.

6
 25 DEGREES 35 DEGREES 45 DEGREES 25 DEGREES 35 DEGREES 45 DEGREES
Downloaded by CARLETON UNIVERSITY LIBRARY on July 17, 2015 | [Link] | DOI: 10.2514/6.2015-3155

Figure 8. Vorticity Contours for angles of attack of 25°, 35° and 45°. Red represents the counterclockwise rotating
vortices. Blue represents the clockwise rotating contours. A-C presents the downstroke and D-F presents the upstroke.

V. Comparison with Experimental and Computational Data

The comparison of instantaneous vortex circulation contours was performed by sectioning the 12th flapping cycle
every 4.2 e-04 s. When comparing the results obtained from the present study with the 2D simulation results gathered
by Z. J. Wang et al.1 shown in the first column, the experiment performed by Dickinson et al. 6 analyzing a wing section
65% of the span with DPIV monitors in the second column, and the 2D numerical simulation by Busra Akay 7 in the
third column, a great similitude is evident in the magnitude and circulation of the vortex structures.

7
 COMPUTATIONAL EXPERIMENTAL COMPUTATIONAL PRESENT COMPUTATIONAL EXPERIMENTAL COMPUTATIONAL PRESENT
WANG DICKINSON BUSRA AKAY STUDY WANG DICKINSON BUSRA AKAY STUDY
Downloaded by CARLETON UNIVERSITY LIBRARY on July 17, 2015 | [Link] | DOI: 10.2514/6.2015-3155

Figure 9. Comparison of computational and experimental vortex circulation at 45° AoA during the downstroke
and upstroke with previous research data.
The color scales of the vortices do not correspond to the exact values in the contours scale for all the columns,
resulting in different magnitudes of wakes around the airfoil. Based on this fact, the comparison was executed in a
qualitative manner more than a quantitative kind of way, as suggested by Busra Akay 7. When performing the mirror
illusion analysis, the first two columns show some discrepancies between the vortices created during the down stroke
and upstroke; whereas the computational results of the third column along with our results present the same effects of
opposite circulation (Fig. 9).

Figure 10. Lift and Drag Coefficients results comparison with Wang et al. 1 and Dickinson et al. 6.

The aerodynamic coefficients were compared and analyzed in Fig. 10, noticing a high resemblance in the curve
tendency and range of values obtained, especially with the preceding numerical simulation 6. Fluctuations are also
present when strictly comparing experimental and computational data, due to the transition from three dimensional to
two dimensional results, the environmental parameters of different flow properties used, and the sensor precision used
to monitor the moment and force during the experiments. As the Wang et al. 1 computational simulations, it can be
observed that the results of the present study capture and recreate the principal characteristics of the vortex dynamics
throughout the flapping cycle.

8
 VI. Effects of Rotation timing
The relationship between the translational and rotational phases establishes the wing motions into advanced,
symmetric, and delayed rotation patterns (Fig. 11). As W. Shyy et al.3 explains it, the advanced rotation is when the
wing rotates before it reverses its direction at the end of each half stroke, rotating the leading edge backwards relative
to the translation, producing an upward component of lift. In the symmetrical rotation, the wing experiences an angle
of attack of 90º at the end of each half stroke, while in the delayed rotation, the wing rotates after it reverses its
direction at the end of each half stroke, making the leading edge rotate forward relative to translation, producing a
downward lift 6.
Downloaded by CARLETON UNIVERSITY LIBRARY on July 17, 2015 | [Link] | DOI: 10.2514/6.2015-3155

Figure 11. Force vectors based on rotation timing. The red vectors describe the direction of the lift generated by
the airfoil through the entire stroke for the advanced, symmetrical and delayed rotation 6.

An advanced rotation results in a positive lift peak at the end of the half stroke, while a delayed rotation results in
a negative lift at the beginning of each half stroke. A symmetrical rotation produces a positive and a negative peak
before and after the stroke reversal respectively, according to Dickinson et al.6, as the angle of attack exceeds 90º, the
force vector dips below the stroke plane and the sign of lift changes from positive to negative. In Fig. 12, both the
advanced and delayed rotation show an offset in the graph, making the beginning point of the wing stroke different
Figure 11.
for each rotation pattern; the advanced rotation begins the stroke at the 0.7 part of the symmetric cycle, while the
delayed rotation starts the stroke at the 0.3 part of that same symmetrical rotation cycle.
Dickinson et al.6 showed that advanced rotation produced large lift peaks at the beginning and near end of each
stroke. In symmetrical rotation, the lift peak at the beginning became smaller and was followed by a dip, and the lift
peak at the end of each stroke also became smaller, as for delayed rotation, no lift peak appeared and a large dip
occurred at the beginning of the stroke 16.

Figure 12. Lift and Drag Coefficients results comparison between symmetric, advanced and delayed rotation
in the straight line trajectory.
9
 The peaks observed in the curves of the types of rotation patterns mentioned, are created due to the interaction of
the wing with the vortices shed throughout the previous wing strokes, affected as well by the angle of attack of the
airfoil regarding the translational movement.
Fig. 13 compares the three rotation timing schemes of the wing during the last cycle, analyzed in the linear
trajectory, beginning with a downstroke at t/T=4.0; it can be observed that during the advanced and delayed rotation
movements, the interaction of the wing with the vortices is accelerated, affecting the loss and gain of energy
respectively.
During the delayed rotation mode, the interaction between vortices increases as a result of the delay in the wing´s
position at the beginning of the stroke, cancelling the suction zone created by vortices A-B in Fig. 8, because the
position of the airfoil and its interaction with the vortices don´t give enough time for such use of energy.
Meanwhile, over the course of the advanced rotation stroke, the LEV is not shed from the leading edge after the
wing has reversed its direction during the translational phase, benefiting from the kinetic energy in the wakes, shown
by the lift peaks in Fig. 12 during the beginning of each half stroke.
Downloaded by CARLETON UNIVERSITY LIBRARY on July 17, 2015 | [Link] | DOI: 10.2514/6.2015-3155

Figure 13. Vorticity Contours for a straight line flapping trajectory. Red represents the counterclockwise rotating
vortices. Blue represents the clockwise rotating contours.

10
 VII. Effects of Flapping Trajectories
Based on the flapping trajectories described in Fig. 1, the wing´s translational pathways were varied in order to
compare and analyze the effects that each trajectory has on the force generated by the 2D wing model.
The four trajectories mentioned were simulated maintaining an angle of attack of 45° and changing another
parameter such as the rotation timing, yielding force coefficient results for delayed, symmetric and advanced rotation.

A. Analysis of Aerodynamic Forces

In order to analyze the behavior of the curves obtained in Figs. 14, 15 and 16, the results of the linear trajectory
are taken as reference. In Fig. 14, the four trajectories have a similar tendency throughout the entire flapping cycle,
varying only the aerodynamic coefficient values, points 1, 2 and 3 in Fig. 7. During the down stroke, the ‘8’ and ‘oval’
or ‘elliptical’ trajectories do not exhibit the first lift peak shown as point 1 in Fig. 7, as they do during the upstroke,
thus reducing the mean aerodynamic coefficients. The other trajectories such as straight line and ‘U’ shaped present
Downloaded by CARLETON UNIVERSITY LIBRARY on July 17, 2015 | [Link] | DOI: 10.2514/6.2015-3155

these peaks throughout the entire flapping cycle, based on the fact that the results and the behavior of the trajectories
are identical during both the downstroke and upstroke. The highest coefficient value is present in the ‘U’ trajectory,
along the entire flapping cycle, mainly because of the new angle formed between the inclination of the stroke plane
and the wing at the end of the half stroke and during the translational phase.

Figure 14. Lift and Drag Coefficient results for linear, oval/ elliptical, U shaped and 8 shaped trajectories with
symmetrical rotation, maintaining 45° angle of attack during the translation phase.

When comparing the curves in Fig. 15, the same tendency as in Fig. 14 can be seen; however, when analyzing the
results between the symmetric and advanced rotation, an increment is evident in the magnitude of the lift and drag
coefficients, most specifically in the ‘8-shaped’ trajectory, consequence of the complete shedding of the LEV halfway
the translational phase, and when the airfoil returns in the upstroke, it recovers the previously shed vortex, increasing
the energy usage and thus increasing the lift. The positive increment of the coefficients can be explained as an effect
of the angle made with regard to the stroke plane due to the advanced rotation, generating the lowest peak of the drag
coefficient. Advanced rotation in the ‘U’ trajectory shows two positive lift peaks during each half stroke, maintaining
their range of magnitude during the upstroke, and acting in form of mirror when compared with the downstroke.

Figure 15. Lift and Drag Coefficient results for linear, oval, U shaped and 8 shaped trajectories with advanced
rotation, maintaining 45° angle of attack during the translation phase.

11
 Regardless of the trajectory performed, at the beginning of each half stroke, the lift will tend to zero, while the
drag will show a negative value as consequence of the interaction between the vortices shed in the wake trail of the
previous strokes and the energy surrounding the airfoil at all times. Generally, how fast the increase in lift is generated
depends of the trajectory and the immediate angle of attack. As a clear example, the ‘8-shaped’ trajectory with
advanced rotation takes the longest for the wing to generate a positive value of the coefficient during the upstroke.

The four trajectories analyzed present similar trends and tendencies shown in Fig. 16, with the exception of the
elliptical trajectory for both lift and drag coefficients due to the surrounding environment of the wing during the
translational phase, which can be considered close to a clean environment where the previous vortices shed are not
used nor taken advantage to increase the aerodynamic forces. Generally, the behavior of the curves is similar to the
ones presented in advanced and symmetric rotation, only differing in the slight decrease of lift and increase of drag,
showing an inefficiency in the translational phase when a delayed rotation is performed.
Downloaded by CARLETON UNIVERSITY LIBRARY on July 17, 2015 | [Link] | DOI: 10.2514/6.2015-3155

Figure 16. Lift and Drag Coefficient results for linear, oval, U shaped and 8 shaped trajectories with delayed
rotation, maintaining 45° angle of attack during the translation phase.

B. Analysis of flow field and vorticity contours

For the vortex analysis between trajectories, as for the one performed in section 5, images of detailed vortex
contours for the 12th flapping cycle were gathered, taking into consideration a consistent 45° of angle of attack for
the translation phase, maintaining a symmetrical rotation. Taking as reference the behavior of the lineal trajectory, the
‘U’ movement is the one with the most similar conduct as a result of having a least significant deviation in the y axis.
The major difference is shown in the oval trajectory, due to the nature of the movement, in which the energy of the
vortices is not harnessed as all other trajectories do.
As seen in Fig. 17, during the first half of the downstroke, the wing doesn´t interact with the wakes left by the
previous flapping cycle in the elliptical and ‘8-shaped’ trajectories, because it starts its movement upwards, leaving
the vortices behind not encountering them again in future cycle repetitions, making this a big difference compared
with the linear and ‘U-shaped’ trajectories, which start their movement in a horizontal manner and downwards
respectively, where the wing interacts more frequently with the vortex wake than in any other trajectory analyzed;
however, the ‘U’ shaped trajectory is the only one that doesn´t keep the recirculation bubble on the upper surface of
the airfoil, described in section B of chapter 4 (Fig. 8, zone 4), during the second half of the downstroke, thus
decreasing the aerodynamic coefficients.

12
Downloaded by CARLETON UNIVERSITY LIBRARY on July 17, 2015 | [Link] | DOI: 10.2514/6.2015-3155

Figure 17. Vorticity Contours comparison for a straight line, oval, U shaped and 8 shaped flapping trajectories.

During the upstroke, it can be seen how the wing-wake interaction is different in the oval trajectory, mainly because
the change in direction is more abrupt than in other trajectories, seeing the wing move upwards, resulting in a small
increase of the lift coefficient, reaching higher values than the ones reached during the linear trajectory, but not as
high as the peaks reached in the drag coefficient. The linear, ‘8’ and ‘U’ shaped trajectories interact with the vortex
wakes, achieving its greatest wing-wake interaction during the downstroke of the ‘U-shaped’ trajectory. For the ‘8-
shaped’ trajectory, the aerodynamic coefficients are lower, consequence of not interacting with the vortices, which are
detached quickly during the second half of the upstroke.
Finally, it can be considered that apart from the linear movement, the ‘U-shaped’ trajectory is the only one that
shows the mirror effect when both upstroke and downstroke are compared, similar to that described in the linear
trajectory, given the dynamics of the translational movement. From these results, the ‘U-shaped’ flapping trajectory
could present the optimal conditions and characteristics in the fruit fly´s translational movement, given the force
patterns and aerodynamic effects.

C. Conclusion
The mathematical model that describes the movement of the fruit fly wing was defined based on the previous
research leading to the numerical simulation. Even though the equations implemented represent the kinematics of the
four trajectories and do not describe 100% accurately the wing movement of the Drosophila, taking into consideration
parameters of frequency, Reynolds number and environment conditions, the flow structure results obtained are valid
and comparable with previous computational and experimental studies.
In the course of the analysis, the main aerodynamic effects were visible throughout the flapping cycle, recognizing
the Wagner effect, which produces the LEV remaining attached to the surface of the airfoil for an estimated time
interacting with the wakes left by the previous cyclic movements. This condition is most notorious in case 3 of a 45°
AoA, due to the fact that with an increase of the angle of attack, the magnitude and vortices interaction is increased.
Similarly, we can observe a mirror effect in the behavior of the fluid during the upstroke and downstroke with the

13
 instantaneous vortex contours analyzed. 2D results simplifying the fruit fly wing kinematics are completely valid and
indeed comparable with experimental and computational data for a hovering configuration.
With the variation of the rotation timing, a control of the aerodynamic effects was possible, due mainly to the
change of the angle of attack of the airfoil, with respect to the translational trajectory, at the beginning and end of each
half stroke. An increase of lift coefficient was observed in the advanced rotations, compared with the symmetrical
rotation, while the delayed rotation experienced a decrease of lift, and a rise of drag before or after the pronation and
supination of the wing.
By adjusting this rotation timing, insects can generate lift making use of the vorticity field from previous strokes.
Within the nature of the fruit fly´s flapping motion, a combination of rotation timings can be used to perform pitch
and yaw turns, adjusting a delayed rotation in one wing and an advanced rotation in the other, increasing the
maneuverability in flight.
Along with the variation of the rotation timing, the change in flapping trajectory of the wing allows the insect to
increase the aerodynamic forces produced, taking advantage of the energy dissipated after every stroke still
Downloaded by CARLETON UNIVERSITY LIBRARY on July 17, 2015 | [Link] | DOI: 10.2514/6.2015-3155

surrounding the wing during hovering flight. An appropriate combination of advanced rotation and U-shaped
trajectory can improve the flight´s efficiency during hovering and steering maneuvers.
Future work implies performing computational simulations in 3D with all the degrees of freedom necessary to
analyze all the parameters in the fruit fly wing kinematics and evaluate all the aerodynamic effects present in this type
of flight.

Acknowledgments
The authors would like to thank Jaime A. Escobar of the Universidad de San Buenaventura for his support and
orientation to fulfill this study. We are also thankful to the Universidad de San Buenaventura for sponsoring this
research.

References
1
M. H. Dickinson, J. Birch and J. Wang, "Unsteady forces and flows in low Reynolds number hovering flight: two-
dimensional computations vs robotic wing experiments," The Journal of Experimental Biology 207, pp. 449-460,
2004.
2
S. N. Fry., R. Sayaman and M. Dickinson, "The aerodynamics of hovering flight in Drosophila," The Journal of
Experimental Biology 208, pp. 203-2318, 2005.
3
W. Shyy, Y. Lian, J. Tang, D. Viieru and H. Liu, Aerodynamics of Low Reynolds Number Flyers, Cambridge
University Press, 2008.
4
J. Seifert, A review of a Magnus effect in aeronautics, S. S. Direct, Ed., Journal Progress in Aerospace Science,
2012.
5
S. P. Sane, "The aerodynamics of insect flight," The Journal of Experimental Biology 206, pp. 4191-4208, 2003.
6
M. Dickinson., F. O. L. and S. P. Sane., "Wing rotation and the Aerodynamic basis of insect flight," Science 284,
p. 1954, 1999.
7
Busra Akay., "Unsteady aerodynamic Calculations of flapping wing motion," Middle East Technical University,
Ankara, 2007.
8
S. P. Sane and M. Dickinson, "The aerodynamic effects of wing rotation and a revised quai-steady model of
flapping," The Journal of Experimental Biology 205, pp. 1087-1096, 2002.
9
S. Ansari., R. Bikowski and K. Knowles., "Aerodynamic modelling of insect - like flapping flight for micro air
vehicles.," Sci Verse Science Direct, 2006.

14
 10
F. O. Lehmann., "When wings touch wakes: Understanding locomotor force control by wake-wing interference in
insect wings," The journal of Experimental Biology 211, pp. 224-233, 2008.
11
D. Lentink and M. H. Dickinson, "Biofluid dynamic sacling of flapping spinning and translating fins and wings,"
The Journal of Experimental Biology 212, pp. 2691-2704, 2009.
12
D. F. Kurtulus, "Numerical and experimental analysis of flapping motion in hover applications to micro air
vehicle." Middle East Technical University, Ankara, 2005.
13
S. Leguizamón, "Parametric simulation of the non - stationary Aerodynamics of flapping wing at low Reynolds
numbers," 2012.
14
W. Shyy et al, “Recent progress in flapping wing aerodynamics and aeroelasticity,” Progress in Aerospace
Sciences, Science Direct, 2010.
Downloaded by CARLETON UNIVERSITY LIBRARY on July 17, 2015 | [Link] | DOI: 10.2514/6.2015-3155

15
J. T. D. Viieru and W. Shyy, "Effects of Reynolds Number and Flapping Kinematics on Hovering Aerodynamics,"
45th AIAA Aerospace Sciences Meeting and Exhibit, 2007.
16
S. P. Sane and M. Dickinson., "The control of Flight Force by a Flapping Wing: Lift and Drag Production" The
Journal of Experimental Biology 204, pp. 2607-2626, 2001.

15