EXPERIMENTAL OPTIMIZATION OF 2-D FLAPPING

PLATE KINEMATICS USING GLOBAL SURROGATE-

ASSISTED GENETIC ALGORITHM WITH KRIGING

THESIS

This paper is submitted

as partial fulfillment for Master degree
from Institut Teknologi Bandung

by
DEA DANIELLA
23613005
(Aeronautics and Astronautics Program)

INSTITUT TEKNOLOGI BANDUNG

2015
 ABSTRACT

EXPERIMENTAL OPTIMIZATION OF 2-D FLAPPING

PLATE KINEMATICS USING GLOBAL SURROGATE-
ASSISTED GENETIC ALGORITHM WITH KRIGING

by
Dea Daniella
23613005
(Aeronautics and Astronautics Program)

In this thesis, an optimization of 2-D flapping plate kinematics under hovering

condition was done to obtain optimum kinematics that maximizes lift and study
the effect of certain kinematics to lift generation. The optimization was done using
G-SAGA with Kriging, where Kriging method constructed surrogate model of
time-averaged lift coefficient as a function of flapping kinematics and Genetic
Algorithm optimized explicit function obtained from surrogate modelling. Fitness
evaluation was done using towing tank experiments.

Four kinematics parameters were split into two test cases of optimization. The
results show that lift generation and enhancement are heavily relied on rotational
timing and translational angle of attack. Force history and flow analysis from PIV
experiments showed that lift is enchanced if the plate rotates prior to the end of a
half-stroke. This mechanism, known as advanced rotation, further increases lift at
the beginning of the following half stroke as the plate hits shed vortices from
previous stroke. Furthermore, maximum lift of flapping plate is achieved in a
certain range of angle of attack between 35º ≤ α ≤ 55º due to delayed stall of LEV
that sustains lift. The effect of stroke amplitude and rotational duration, however,
are less apparent than the other two parameters.

Keywords: Genetic Algorithm, Kriging, surrogate model, experimental fluid

dynamics, hovering insect flight, flapping plate

i
 LEGALIZATION PAGE

EXPERIMENTAL OPTIMIZATION OF 2-D FLAPPING

PLATE KINEMATICS USING GLOBAL SURROGATE-
ASSISTED GENETIC ALGORITHM WITH KRIGING

by
Dea Daniella
23613005
(Aeronautics and Astronautics Program)

Institut Teknologi Bandung

THESIS APPROVAL
This thesis has been approved and legalized as partial fulfillment for Master
degree in Aeronautics and Astronautics Program.

Bandung, August 26, 2015

___________________________
(Dr. Lavi R. Zuhal)

ii
 THESIS USER GUIDE

This master degree thesis is registered and available in the library of Institut
Teknologi Bandung (ITB). This thesis is open for public with copyright on the
author and by following the rules of intellectual property rights in Institut
Teknologi Bandung (ITB). Literature references are allowed to be recorded but
citations are only allowed with permission from the author. Publishing or
reproducing a part or this entire thesis is prohibited and is only allowed with
consent from Director of Graduate Programme of Institut Teknologi Bandung
(ITB).

iii
 ACKNOWLEDGEMENT

I would like to thank my advisor, Dr. Lavi R. Zuhal for his invaluable guidance,
encouragement, and support during the process of research for this thesis. I would
also like to thank my thesis examiners: Dr. M. Agoes Moelyadi, Dr. Rianto Adhy
Sasongko, and Dr. Djoko Sarjadi, and also Dr. Leonardo Gunawan, who has
kindly helped me throughout seminar and thesis defense.

I am forever grateful for the support provided by my friends and relatives during
my study: my parents, my sister, friends from Department of Aeronautics and
Astronautics, and lab mates from Aerodynamics Laboratory. In particular, I wish
to thank my research partner, Yohanes Bimo Dwianto, for providing the
optimization algorithm, Stepen, whose image processing code was used for flow
analysis in this thesis, and Adityo Saputra, whose generous assistance helped me
understand experimental design, instrumentation, and electrical engineering aspect
of my thesis far better. I would also like to thank our lab technician: Pak Dali and
Mas Adi, who helped me solve several technical problems that arose while I was
working with this thesis.

This research was funded by Institute of Technology Bandung’s Riset KK 2015

program. My graduate student scholarship was provided by DIKTI (Directorate
General of Higher Education of Indonesia), to whom I am grateful.

Bandung, August 26, 2015

Dea Daniella

iv
 LIST OF CONTENTS

ABSTRACT ............................................................................................................. i
LEGALIZATION PAGE ........................................................................................ ii
THESIS USER GUIDE ......................................................................................... iii
ACKNOWLEDGEMENT ..................................................................................... iv
LIST OF CONTENTS ............................................................................................ v
LIST OF FIGURES .............................................................................................. vii
LIST OF TABLES ................................................................................................. ix
LIST OF SYMBOLS .............................................................................................. x
Chapter I Introduction ......................................................................................... 1
I.1 Historical Background .......................................................................... 1
I.2 Objectives .............................................................................................. 5
I.3 Scope of Work....................................................................................... 6
I.4 Research Methodology.......................................................................... 6
I.5 Thesis Outline ....................................................................................... 7
Chapter II Fundamentals of Flapping Flight .................................................... 9
II.1 Introduction ........................................................................................... 9
II.2 Terminology and Kinematics of Flapping Flight .................................. 9
II.3 Equation of Flow ................................................................................. 11
II.4 Scaling Parameters .............................................................................. 13
II.5 Attributes of Flapping Flight Aerodynamics ...................................... 14
Chapter III Fundamentals of Surrogate – Assisted Genetic Algorithm ........... 18
III.1 Introduction ......................................................................................... 18
III.2 Fundamentals of Natural Optimization ............................................... 18
III.3 Single-Objective Real Genetic Algorithm (SOGA) ............................ 19
III.4 Surrogate – Assisted Genetic Algorithm (SAGA) .............................. 21
III.5 Kriging ................................................................................................ 23
III.6 Global Surrogate – Assisted Genetic Algorithm with Kriging ........... 26
Chapter IV Experimental System Design and Setup ....................................... 27
IV.1 Introduction ......................................................................................... 27
IV.2 Towing Tank Experiment ................................................................... 27

v
 IV.3 Flapping plate model ........................................................................... 29
IV.4 Translational System ........................................................................... 29
IV.4.1 Two-Stage PID Controller ................................................................ 31
IV.4.2 Limit Switch and Logic Gate ............................................................ 33
IV.4.3 Ground Loop ..................................................................................... 34
IV.4.4 Encoder Pulse Divider....................................................................... 36
IV.5 Rotational System ............................................................................... 37
IV.6 Force Measurement System ................................................................ 37
IV.7 Particle Image Velocimetry................................................................. 39
Chapter V Experimental Optimization of Flapping Plate Kinematics using G-
SAGA with Kriging .............................................................................................. 43
V.1 Introduction ......................................................................................... 43
V.2 Preliminary Experiment and Validation.............................................. 43
V.3 Experimental Optimization of Flapping Plate Kinematics using G-
SAGA with Kriging.......................................................................................... 45
V.3.1 Case 1: The Effect of Plate Rotation ................................................. 47
V.3.2 Case 2: The Effect of Stroke Amplitude and Angle of Attack ......... 56
V.4 Result Analysis and Discussion .......................................................... 58
V.4.1 Case 1: The Effect of Plate Rotation ................................................. 65
V.4.2 Case 2: The Effect of Stroke Amplitude and Angle of Attack ......... 68
Chapter VI Conclusion and Future Works....................................................... 72
VI.1 Conclusion........................................................................................... 72
VI.2 Future Works ....................................................................................... 73
REFERENCES...................................................................................................... 74

vi
 LIST OF FIGURES

Figure I.1 (a) Ornithopter; (b) Lilienthal's glider, taken from Shyy (2013)............ 1
Figure I.2 Examples of biologically inspired MAV, taken from Shyy (2013) ....... 2
Figure II.1 The terminology in insect flight (Sane, 2003) ...................................... 9
Figure II.2 Kinematics of 3-D flapping flight (Shyy, et al., 2010) ....................... 10
Figure II.3 Kinematics of 2-D hovering flight (Sane, 2003)................................. 11
Figure II.4 Clap and fling mechanism (Shyy, et al., 2010) ................................... 15
Figure II.5 Illustration of wake capture (Trizila, Kang, Aono, & Shyy, 2011) .... 16
Figure II.6 Illustration of wake dynamics during one stroke ................................ 16
Figure III.1 Flowchart of (a) basic real, single objective Genetic Algorithm ....... 23
Figure IV.1 Towing tank in Laboratory of Aerogasdynamics .............................. 27
Figure IV.2 Schematic of towing tank experimental system ................................ 28
Figure IV.3 Two-stage PID system for AC motor control.................................... 31
Figure IV.4 (a) Logic gate for limit switches and relay (b) PCB containing ........ 34
Figure IV.5 Electrical schematic of voltage coupler ............................................. 35
Figure IV.6 Encoder waveform output (Autonics Corporation, 2015) ................. 36
Figure IV.7 Decomposition of normal force, measured by load cell .................... 38
Figure IV.8 PIV system for towing tank experiment ............................................ 40
Figure IV.9 Basic PIV scheme for flowfield measurement .................................. 41
Figure V.1 Schematic of the reproduction of Lisoski's normal flat plate ............. 44
Figure V.2 (a) Drag coefficient Cd comparison of present experiment to ............ 45
Figure V.3 Flowchat of flapping plate optimization algorithm using real ............ 46
Figure V.4 (a) Typical surface and (b) contour of time-averaged lift ................... 49
Figure V.5 Histogram of optimum set of rotational kinematics in Case 1 ........... 49
Figure V.6 Experimental result of first optimum solution of Case 1 .................... 50
Figure V.7 Experimental result of second optimum solution of Case 1 ............... 52
Figure V.8 Experimental result of suboptimum solution of Case 1 ...................... 54
Figure V.9 (a) Typical surface and (b) contour of time-averaged lift ................... 57
Figure V.10 Histogram of optimum set of flapping kinematics in Case 2 ........... 57
Figure V.11 Experimental result of first optimum solution of Case 2 .................. 59
Figure V.12 Experimental result of second optimum solution of Case 2 ............. 61

vii
Figure V.13 Experimental result of suboptimum solution of Case 2 .................... 63
Figure V.14 Illustration of lift generation due to plate pronation ......................... 65
Figure V.15 Velocity vector on the end of a half stroke in (a) Case 1.................. 67
Figure V.16 Result of parametric study of the effect of rotational kinematics ..... 67
Figure V.17 (a) Result of parametric study of the effect of stroke amplitude ...... 69

viii
 LIST OF TABLES

Table IV.1 Two-stage PID constants .................................................................... 33

Table IV.2 Karnaugh map for limit switch logic function .................................... 34
Table IV.3 Prosilica GX1050 setting for image acquisition ................................. 42
Table IV.4 Image processing parameter configuration in PIV_GUI.m ................ 42

ix
 LIST OF SYMBOLS

ABBREVIATION Name First appeared

in page

AC Alternating current 29
ADC Analog-to-digital converter 37
AR Aspect ratio 13
CCD Charge-coupled device 7
CDI Central Difference Interrogation 42
CFD Computational Fluid Dynamics 4
DAQ Data acquisition 37
DOF Degree of freedom 3
FFT Fast Fourier Transform 42
GA Genetic Algorithm 3
G-SAGA Global Surogate-Assisted Genetic 5
Algorithm
IA Interrogation area 42
IC Integrated circuit 29
IDE Integrated development 31
environment
ISR Interrupt Service Routine 36
LEV Leading edge vortex 17
MAV Micro Aerial Vehicle 2
MOSFET Metal-oxide-semiconductor field- 35
effect transistor
PCB Printed circuit board 29
PD Proportional-derivative 32
PID Proportional-integral-derivative 31
PIV Particle Image Velocimetry 7
P/R Pulse per revolution 32
PRS Polynomial Response Surface 4
PWM Pulse width modulation 30
RBNN Radial Basis Neural Network 4
RGA Real-coded Genetic Algorithm 6
RMS Root mean square 45
RPM Rotation per minute 30
SA Search area 42
SBX Simulated Binary Crossover 21
SOGA Single-Objective Genetic Algorithm 19
SQP Sequential Quadratic Programming 3
SVR Support Vector Regression 4
SWF Forward switch 33
SWR Reverse switch 33
TEV Trailing edge vortex 17
TiV Tip vortex 17

x
ROMAN SYMBOL Name First appeared
in page

b Wing span 29
C Aerodynamic force coefficient 38
c Chord length 13
D Drag 38
d Distance of stroke 13
e Error 31
F Force vector 12
f Frequency 13
K Control system constant 31
k Reduced frequency 14
L Likelihood function, lift 24, 38
N Number of stroke 47
n Number of poles 30
n Normal vector 12
P Pressure 11
p Number of pulse 32
p Exponent vector 24
r Distance 12
Re Reynolds number 11
S Wing area 29
St Strouhal number 13
T Stroke period 47
t Time 11
U Flow velocity 13
u Velocity vector 11
v Carriange’s translational velocity 30
v Induced velocity vector 12

GREEK SYMBOL

α Angle of attack 10
β Stroke plane 10
Γ Circulation 12
Δ Difference 32
θ Basis vector 24
θ Deviation angle 10
λ Regression constant 25
μ Mean 24
ν Kinematic viscosity 13
ρ Density 12
σ Square root of variance 24
φ Stroke angle 10
χ Body angle 10
Ψ Gram matrix 24
ψ Basis function 24

xi
GREEK SYMBOL Name First appeared
in page

ω Vorticity vector 12

ACCENT AND
SUBSCRIPT

~ Non-dimensional variable 11
^ Surrogate 22
¯ Time-averaged 44
∞ Free stream 13
a Amplitude 13
d Derivative, drag 31, 38
i Integral 31
l Lift 38
m Mean 13
n Normal 38
opt Optimum 26
out Output 31
p Proportional 31
pop Population 20
r Rotation 47
read Read or measured from sensor 31
ref Reference 13
subopt Suboptimum 56
var variable 20

xii
 Chapter I Introduction

I.1 Historical Background

Since the beginning of century, humans have been fascinated by biological flyers.
They longed to fly like birds and tried to imitate them, but to no avail. Icarus, a
Greek mythological character, was depicted as a man with a pair bird-like wings
strapped on his back and his arms. A sketch by Leonardo da Vinci in 16th century
showed a human-powered “flying machine” or “ornithopter”, a bat-like wing
attached to one’s back with pedals to generate thrust. Three centuries later, a stork
inspired Otto Lilienthal’s infamous glider.

Figure I.1 (a) Ornithopter; (b) Lilienthal's glider, taken from Shyy (2013)

Biological flyers indeed possess desirable flight characteristics; the one that stands
out most is probably their high maneuverability. While a highly aerobatic aircraft
has a roll rate of approximately 720°/s, a swallow (Hirundo rustics) has a roll rate
in excess of 5000°/s. Some birds can withstand 10 to 14 G-forces, whereas the
maximum G-forces in some military aircrafts are limited to 8 or 10 G (Shyy, Lian,
Tang, Vieru, & Liu, 2008). To perform those sophisticated maneuvers, biological
flyers have complex wing mechanism that is difficult to be imitated by man-made
machines. They flex, flap, twist, sweep, and plunge their wings in order to adjust
their flying speed while sustaining their body weight, generating lift and
propulsive force altogether. These motions create unsteady and complicated fluid
flow around their wings, involving vortex formation and interaction, yet another
challenging topic to be studied in aerodynamics.

1
One of the pioneers of quantitative study of biological flight was Weis-Fogh
(1973), who indentified clap-and-fling mechanism in lift generation of wasps
(Encarsia formosa). Similar mechanism in insect lift generation was also
observed by Ellington (1984), who proposed “vortex models” to study the
aerodynamics of insect wing analytically. In this model, insect wing was
approximated with propeller blade in quasi-steady condition. Aside from clap-
and-fling mechanism, Ellington also observed a lift enhancement mechanism due
to leading edge vortex formation called “delayed stall”, which was later confirmed
by Dickinson (1999) from his experiment in fruit fly. Along with Sane (2001),
Dickinson identified another specific kinematics of flapping wing called rotational
forces, which increases lift when the wing rotates at the end of each half-stroke.
Dickinson also identified “wake capture”, which increases lift as the wing
undergoes stroke reversal and hits shed vortex from previous stroke.

Figure I.2 Examples of biologically inspired MAV, taken from Shyy (2013)

In engineering, study of biological flyers often goes hand in hand with its
application: MAVs (Micro Aerial Vehicles). MAV was originally defined as
light-weight aerial vehicle which has maximum dimension of 15 cm, flight speed
between 10-20 m/s, flies at low Reynolds number region (similar to biological
flyers), and is equipped with a camera and sensor to gather information (Shyy,

2
Lian, Tang, Vieru, & Liu, 2008). Today, MAVs are developed for surveillance,
reconnaissance, homeland security, and sensing at hazardous environments. Based
on its wing, MAV can be divided into 3 categories: fixed wing, rotary wing, and
flapping wing. Compared with the other two, MAV with flapping wing still leaves
many rooms to be explored, mainly because of its lack of effectiveness compared
to natural flyers. Although biological flyers have complex wing kinematics and
structure as aforementioned, it is quite possible for MAV to copy wing kinematics
and shape from them using advanced technology in MAV design at present.
However, MAV is still proved to be less effective than those created by Mother
Nature throughout years of evolution (Jones & Yamaleev, 2013). Therefore, in
order to obtain a highly maneuverable, effective, and agile MAV, optimization is
needed after preliminary conceptual design.

The study about optimization of flapping wing mostly concerns itself with
maximizing lift, thrust, or propulsive efficiency by varying design variables: wing
kinematics parameters (pitch angle, phase angle, flapping frequency), wing shape
parameters (wingspan, chord length, wing thickness, dihedral angle, twist angle),
and global flow parameters (Reynolds, Strouhal, and Mach number). Because of
the massive number of design variables, finding a function that relates
aerodynamic forces and/or propulsive force to those variables will be very
difficult. Therefore, conventional gradient-based optimization is unsuitable for
this case. The proposed optimization method for MAV design is Genetic
Algorithm (GA), a global non-gradient method based on the principle of genetics,
evolution, and natural selection.

The utilization of GA for flapping flight aerodynamic optimization has been done
in several previous researches. Ito (2002) optimized pitching and plunging
frequency, amplitude, and phase shift of a rectangular rigid flat plate with 2 DOF
to obtain maximum propulsive efficiency. In his work, Ito applied GA and
combined it with SQP (Sequential Quadratic Programming), a gradient-based
optimization method. Similar object was also studied by Milano and Gharib
(2005), albeit with different kind of motion and kinematic parameters. Milano

3
used GA to optimize translational velocity, travelling distance, angular velocity,
and angle of attack of plate which was translated and rotated back and forth,
mimicking flapping motion, in a towing tank experiment. Coupling between GA
and experimental method was also done by Salles (2004) to optimize kinematics
parameters of 3 DOF robotic model of hawkmoth (Manduca sexta) wing used for
planetary exploration. Both of these works optimized flapping kinematics in order
to obtain maximum lift.

Similar to Ito’s work, Choi (2011) applied GA to optimize flapping frequency,

angle of attack, and phase shift of a pitching and plunging 2D flat plate for
maximum thrust. In his work, Choi combined GA with a surrogate modelling
technique called Kriging method which served the purpose of reducing
computational cost. This, however, was not the first time surrogate models were
used in flapping wing research. Surrogate models were previously incorporated in
the researches of Trizila (2008), where several surrogate models, such as PRS,
RBNN, Kriging, and SVR, were used to map the effects of 2-D airfoil kinematics
parameters (angle of attack, stroke amplitude, and phase lag) on lift and power
generation. Same research was conducted and comparison between 2-D and 3-D
wing was drawn as a continuation of this research (Trizila P. , Kang, Visbal, &
Shyy, 2008). All of the works mentioned above optimized symmetric hovering
kinematics of insect flight and used CFD method to evaluate fitness function.

In Laboratory of Aerogasdynamics, Institute of Technology Bandung, GA has

been developed and studied in several cases of aerodynamic optimization. GA
was employed in the works of Adhynugraha (2009) and Dwianto (2014) to
optimize airfoil geometry for design purpose. Another case is flapping flight
optimization, which is based on previous researches in the same laboratory about
flow characteristics around a 2-D flapping plate by Nguyen (2006), Wuwung
(2008), and Nguyen (2009). Their works showed that flow characteristics (e.g.
velocity and vorticity field) around flapping plate are heavily influenced by
flapping plate kinematics. Furthermore, it was also demonstrated by Nguyen
(2009) that certain kinematics are able to enchance force generated by flapping

4
plate. Using this insight, Palar (2014) developed Micro-GA, which was coupled
with experimental method, to optimize flapping flight kinematics for maximum
thrust generation. Furthermore, Palar suggests that optimization alone is not
enough: parametric study is needed to study the effect of flapping kinematics on
force production adequately (Palar, 2012).

It was suggested by previous researches in insect flight that hovering insect has
peculiar, yet sophisticated kinematics (Dickinson, Lehmann, & Sane, 1999).
Furthermore, Ellington (1984) theorizes that lift generated solely by flapping
during hover must be the highest since there is no additional lift from insect
body’s forward velocity. Therefore, in this thesis, I wish to study how a set of
kinematics parameters affect lift generation of a flapping plate and investigate
optimum kinematics that generate maximum lift during hover. Kinematics
parameters that were studied in this thesis are rotational parameters, stroke
amplitude, and angle of attack. Kriging method was employed to construct
surrogate model of lift generated by flapping plate. Finally, Genetic Algorithm
(GA) optimization with the assistance of surrogate modelling (G-SAGA with
Kriging) was done to find a set of optimum kinematics that generates maximum
lift. It is hoped that the combination of GA and surrogate modeling will increase
searching efficiency of GA, so that optimum solution can be obtained with lower
number of fitness evaluation. Computational efficiency was increased further by
combining optimization algorithm with experimental method (instead of CFD) for
fitness evaluation. All in all, this research is intended to cast some light in
flapping wing design for MAV application, especially from aerodynamics’ point
of view and its correlation to kinematics of flapping wing, and also to improve the
efficiency of an engineering design tool using Genetic Algorithm coupled with
surrogate modelling.

I.2 Objectives
The objectives of this thesis are:
(1) To investigate the effect of rotational parameters, stroke amplitude, and
angle of attack of a 2-D flapping plate on lift generation.

5
 (2) To obtain a set of optimum kinematics parameters that maximizes lift of 2-
D flapping plate under hovering condition.
(3) To study the flow physics of 2-D flapping plate under optimum condition
and its correlation to lift generation and enhancement mechanism.

I.3 Scope of Work

To limit the scope of work, optimization of 2-D flapping plate kinematics using
G-SAGA with Kriging were done as follows:
(1) The optimization was done for 2-D rectangular flapping plate model under
symmetric hovering condition, where Re = 1200 and stroke plane was
horizontal, in two test cases, where each case consists of two kinematics
parameters.
(2) The 2-D flapping mechanism consists of translational (downstroke and
upstroke) and rotational (pronation and supination) mechanism.
(3) The optimization was done using single objective Global Surrogate-
Assisted Genetic Algorithm (G-SAGA) with Kriging as surrogate
modeling method and real GA (RGA) as optimization algorithm to
maximize time-averaged lift coefficient of 2-D flapping plate.
(4) Fitness evaluation during optimization was done in towing tank
experiments.
(5) Force measurement system used for fitness evaluation does not separate
inertial force from plate’s acceleration and plate’s aerodynamic force

I.4 Research Methodology

This research is divided into 3 major parts: experimental system development,
optimization, and PIV experiments. Experimental system consists of 3
subsystems: translational, rotational, and force measurement system. A
rectangular flat plate was translated and rotated using a 3-phase AC motor and a
servo, respectively, to mimic 2-D flapping motion in a water-filled towing tank. A
load cell was used to measure force experienced by plate. A code was developed
in MATLAB 2010a to control translational and rotational mechanism via Arduino
UNO board and to obtain force data from load cell via National Instruments data

6
acquisition device according to kinematics parameter input to the said code. This
code is able to work simultaneously whenever fitness evaluation is needed.

After the experimental system for fitness (force) calculation had been ready,
optimization was done using single objective G-SAGA with Kriging. The code
was first developed by Dwianto (2015) with several adjustments to connect the
algorithm with experimental system. The optimization was done on two test cases,
each producing a set of optimum kinematics solution. The first test case optimized
and examined the effect of rotational parameters (time when rotation starts and the
duration of rotation) on lift generation. The second test case optimized and
examined the effect of stroke amplitude and angle of attack. For comparison, a
suboptimum kinematics was picked from solution space of each test case.

Two-dimensional flow measurement and visualization was conducted using PIV

experiments for both optimum and suboptimum kinematics obtained from each
test case to draw connection between force history and the physics of flow. Image
acquisition was done using a CCD camera. A code was created in MATLAB
2010a to control image acquisition parameters simultaneously with experimental
system. Image processing was done using cross correlation algorithm enhanced
with advanced processing method made by Stepen (2011).

I.5 Thesis Outline

This thesis consists of six chapters. The following describes the general outline of
this thesis and a brief description about each chapter:

CHAPTER I: INTRODUCTION
This chapter contains historical background, objectives, scope of work, research
methodology, and the outline of the thesis.

CHAPTER II: FUNDAMENTALS OF FLAPPING FLIGHT

This chapter contains brief explanation about theoretical background of insect’s
flapping wing, including some useful terminology and scaling parameters

7
commonly used in flapping flight research, derivation of Navier-Stokes equation
to obtain force from unsteady flapping kinematics, and some unique attributes of
flapping aerodynamics.

CHAPTER III: FUNDAMENTALS OF SURROGATE-ASSISTED GENETIC

ALGORITHM
This chapter contains brief explanation about theoretical background of GA and
surrogate modelling technique, especially Kriging method. The explanation
continues with coupling between GA and a global surrogate modelling technique
(G-SAGA) that was used in this thesis.

CHAPTER IV: EXPERIMENTAL SYSTEM DESIGN AND SETUP

This chapter contains detailed explanation about the design process and setup of
experimental system, including: translational, rotational, force measurement, and
image acquisition system. Some design choices, design problems that occurred
during design process, and their solutions are discussed briefly in this chapter.

CHAPTER V: OPTIMIZATION OF 2-D FLAPPING PLATE USING G-SAGA

WITH KRIGING
This chapter begins with explanation about preliminary towing tank experiment in
order to validate force measurement system for fitness calculation during
optimization process. Coupling between G-SAGA and experimental method,
elaboration about each test case and its optimum solution, a closer look on flow
physics of optimum solutions obtained from PIV experiments, and lengthy
discussion about the results fill the latter part of this chapter.

CHAPTER VI: CONCLUSION AND FUTURE WORKS

This chapter contains conclusion of present research and some recommendation
regarding future research.

8
 Chapter II Fundamentals of Flapping Flight

II.1 Introduction
In this chapter, the fundamentals of flapping flight will be explained. Due to the
limitation of research scope in this thesis, this chapter will focus solely on 2-D
insect flight kinematics. This chapter consists of four subchapters. This chapter
will be started with the terminology and kinematics of real-life, 3-D insect flight
and its simplification to 2-D kinematics. Formulation of aerodynamic forces,
generated from flight kinematics and derived from Navier-Stokes equation, and
some scaling parameters commonly used in flapping wing research will follow
next. This chapter will be ended by several attributes of flapping flight that, as
shown by previous researches, contribute to lift enhancement.

II.2 Terminology and Kinematics of Flapping Flight

Before getting into the formulation of flapping flight kinematics, it is necessary to
define the terminology of flapping flight first. The majority of flapping flight
terminology is conveniently borrowed from fixed wing aerodynamics. As shown
in Figure II.1, “Wing span” refers to the length of straight line from tip to tip of
the wing. “Wing length” refers to the base-to-tip length of the wing. Usually, wing
span is twice as long as wing length, therefore ignoring the width of thorax.
“Wing chord” refers to the length of straight line that connects leading edge and
trailing edge of the wing, where the length might vary along wing span (Sane,
2003).

Figure II.1 The terminology in insect flight (Sane, 2003)

9
While in reality the kinematics of insect flight is complex and might be different
for every species, the flapping motion of insect can be simplified into rotational
motion in three axes. Consider a 3-D rigid body consists of wing and body of an
insect in Figure II.2. Body angle χ refers to the inclination of insect body relative
to the ground. Stroke plane angle β refers to the inclination of stroke plane relative
to the ground. In this plane, the three rotational motions are defined as follows:
stroke angle φ, which refers to the angle of rotation along x-axis; deviation angle
θ, which refers to the angle of rotation along z-axis; and angle of attack α, which
refers to the angle of rotation along y-axis (Shyy, et al., 2010). In most studies
regarding flapping flight, the term angle of attack itself normally refers to
“geometric angle of attack”, or the angle between wing chord and free-stream
velocity vector, instead of “aerodynamic angle of attack”, where the influence of
downwash velocity is included in the measurement (Sane, 2003).

Figure II.2 Kinematics of 3-D flapping flight (Shyy, et al., 2010)

The kinematics of 3-D insect flight can be simplified further into 2-D kinematics,
shown in Figure II.3. Consider the kinematics of insect flight where the wing
chord rotates along z-axis and wing chord axis. The rotation along z-axis can be
seen as a back-and-forth translational motion in the direction of x-axis and the
rotation about wing chord can be seen as the change of angle of attack along wing
trajectory. This kinematics consists of two motions per stroke: downstroke and

10
upstroke. “Downstroke” refers to dorsal-to-vental wing motion, whereas
“upstroke” refers to ventral-to-dorsal wing motion.

Figure II.3 Kinematics of 2-D hovering flight (Sane, 2003)

Between transitions of each half-stroke, the wing rotates along the axis of the
chord, therefore changes the angle of attack. The term “supination” refers to a
rapid wing rotation during downstroke-to-upstroke transition, whereas the term
“pronation” refers to a rapid wing rotation during upstroke-to-downstroke
transition.

In this thesis, 2-D insect kinematics is applied to “infinite wing” in order to

simulate chord-wise flow in a 2-D wing section, instead of span-wise flow in 3-D
wing. This condition can be achieved experimentally through “grazing condition”,
where the tip of finite wing is placed as close as possible to the wall. This setup
limits span-wise flow, creating a “nominally” 2-D flow about the wing chord
(Lisoski, 1993). It is also worth noting that while infinite wing cannot actually
perform flapping motion, the study of animal flight has benefitted from 2-D flow
assumption, particularly in the case where wing aspect ratio is high (Sane, 2003).

II.3 Equation of Flow

Consider a 2-D body submerged in fluid. Assuming incompressible 2-D flow,
non-dimensional Navier-Stokes equation gives
~
u
t~  u 
~~
~   ~ ~ 1 ~ 2~
u  P 
Re
u

(II.1)
where the right hand side terms represent pressure and viscous term, respectively.
Reynolds number Re is defined as the ratio of fluid inertia to viscous dissipation.

11
Using this non-dimensional formulation, it is easier to see that as the Re,
depending on flight condition, becomes smaller, viscous term becomes more
important than pressure term, and vice versa.

While Navier-Stokes equation is the basis of flow equation, it is not easy, from
experimental fluid mechanics point of view, to measure pressure distribution
around an object. Therefore, it is more convenient to rewrite Navier-Stokes
equation in terms of vorticity
~ ~
ω
   u ~  1 
~ ω ~ 2~
ω
~
t Re
(II.2)
where the vorticity ω is defined as
~ ~
~ 
ω u
(II.3)

To formulate force equation in insect wing, thin airfoil theory will be used.
Integrating small vorticity elements over the surface area S of an airfoil, bounded
by closed contour Σ
~ ~
 u~  dl   ω~  n~dS
 S

(II.4)
Term Γ, known as “circulation”, is obtained. Using Kutta-Jokowski theorem with
the addition of circulation term, force can be defined as

d d
F   
dt R
r  ωdR    vdA
dt S

(II.5)
where the first and second right hand term represent the force risen from vorticity
created by movement of airfoil and inertial force of fluid displaced by airfoil
body, respectively. The term v is defined as induced velocity, which is caused by
the presence of non-zero vorticity ω of a small fluid element dV at a distance r

12
from the region of interest. According to Biot-Savart law, this velocity can be
defined as
1 r
v 
4 V
ω  3 dV
r
(II.6)
Detailed derivation of aerodynamic force in flapping wing is elaborated in Sane’s
work (2003).

II.4 Scaling Parameters

The study of physical phenomena in nature is often conducted with scaling
parameters. Scaling is necessary to reduce the number and the complexity of
parameters involved by combining those parameters into a set of dimensionless
parameters (Shyy, et al., 2010). In the study of flapping wing aerodynamics, there
are 3 commonly used scaling parameters, such as:
(1) Reynolds number Re, which is defined as the ratio of fluid inertia to
viscous dissipation,
U ref Lref U cm
Re  
 
(II.7)
where free stream velocity U∞ and mean chord cm are commonly used as
reference velocity and length, respectively, and kinematic viscosity of
fluid flown is denoted by ν. For physically scaled objects with identical
Re, the value of non-dimensional fluid properties and force is the same
(Sane, 2003). Reynolds number of insect flight typically varies from
O(101) to O(104) (Shyy, Lian, Tang, Vieru, & Liu, 2008).
(2) Strouhal number St, which is defined as the ratio of forward flight speed
U∞ to flapping speed,
fLref fARc m d a
St  
U ref 2U 

(II.8)
where f, AR, and da represent flapping frequency, wing aspect ratio, and
full stroke amplitude (in radians). Strouhal number characterizes vortex

13
 formation and shedding during forward flight. Typical Strouhal number of
insect flight ranges from 0.2-0.4 (Shyy, et al., 2010).

(3) Reduced frequency k, which is a flapping flight parameter similar to

Strouhal number with better representation of flow unsteadiness. Reduced
frequency is defined as
fcm fcm
k 
U ref U
(II.9)

II.5 Attributes of Flapping Flight Aerodynamics

As mentioned in the previous chapter, force generation in flapping wing is very
much related to complex aerodynamics, i.e. vortex formation and interaction.
Therefore, in this part, we will discuss several physical phenomena in flapping
flight caused by vortex dynamics, such as:
(1) Clap and fling
Clap and fling motion, illustrated in in Figure II.4, was found by Weis-
Fogh (1973) as lift enhancement mechanism in wasp (Encarsia formosa).
During “clap” phase, a pair of wing rotates about leading edge, closing the
gap between wings. Consequently, during “fling phase”, a pair of wing
rotates about trailing edge, forming a gap between them (Percin, Hu, van
Oudheusden, Remes, & Scarano, 2011). Lift enhancement is generated
from flow circulation between wings. This mechanism is commonly found
in wasp, fruitfly, butterfly, and hawkmoth, to name a few (Shyy, et al.,
2010).

(2) Rapid pitch rotation

Rapid pitch rotation is another lift enhancement mechanism in hovering
flight. This theory suggests that wing’s own rotation will act as a source of
flow circulation around the wing, thus generates additional lift. At the end
of each half stroke, wing rotates rapidly, causing a phenomenon similar to

14
 Figure II.4 Clap and fling mechanism (Shyy, et al., 2010)

Magnus effect, a curved trajectory on spinning baseball. Baseball rotation

forms circulation around the ball, increasing fluid velocity at one surface
and simultaneously decreasing fluid velocity at the other surface.
“Backspin” is the term used if the velocity is higher at the top surface, thus
the lower pressure will generate upward force component or the lift.
“Topspin” is the term used if the velocity is higher at the lower surface,
thus higher pressure on top surface will generate downward force
component.

The same analogy can be applied in the case of insect wing rotation.
Dickinson (1999) observed that rotational pattern (direction and timing)
heavily influences force generation. Dickinson’s research shows that the
timing of stroke reversal contributes greatly to the aerodynamic force peak
at each stroke. “Advanced” rotation is the term used if the wing rotates
shortly before the end of half-stroke, where the leading edge rotates
backward relative to stroke direction (“backspin”), thus increasing
resultant force in the direction of translation and creating force peak at the
end of each stroke. “Delayed” rotation is the term used if the wing rotates
shortly after half-stroke starts, which, as opposed to advanced rotation,
reduces lift because the wing rotates forward relative to the stroke
direction (“topspin”).

15
 (3) Wake capture
Shortly after the wing rotates and half-stroke ends, another force peak
occurs at the beginning of new stroke (see Figure II.5 and Figure II.6).
This phenomenon is characterized as “wake capture” by Dickinson (1999),
or “wing-wake interactions” (Sane, 2003). During stroke reversal, wing
hits shed vortex from previous stroke, recovering energy from the fluid
that was lost during previous stroke and therefore generating higher force.

Figure II.5 Illustration of wake capture (Trizila, Kang, Aono, & Shyy, 2011)

Figure II.6 Illustration of wake dynamics during one stroke (Birch &
Dickinson, 2003)

(4) Delayed stall of LEV (Leading Edge Vortex)

The presence of LEV is a common trait of flapping flight with Re ≤ O(104)
(Liu & Aono, 2009). During translation of a high angle of attack flapping
wing, flow separation occurs, forming LEV and TEV (Trailing Edge
Vortex) at both edges of wing, and also TiV (Tip Vortex) at wingtip.

16
However, in 2-D wing translation, only LEV and TEV are formed.
Assuming a 2-D insect wing can be modeled with a thin, sharp-edged
airfoil, vortex generated at the leading edge has lower pressure on its core,
resulting in an additional suction force perpendicular to the wing.
Ellington (1996) suggested that this mechanism, due to the presence of
LEV, increases lift at high angle of attack, which is the characteristic of
insect flight. Further research (Liu, Ellington, Kawachi, Van den Berg, &
Willmott, 1998) provides better understanding towards the effect of LEV
in lift enhancement and the term “delayed stall”. Delayed stall refers to the
condition where LEV is still attached to the wing (stall is delayed) and lift
increases. As the wing translates, the LEV is shed to wake and stall occurs,
resulting in a drop in lift (Sane, 2003).

17
Chapter III Fundamentals of Surrogate – Assisted Genetic
Algorithm

III.1 Introduction
This chapter contains fundamentals of optimization algorithm used in flapping
plate kinematics optimization problem: Genetic Algorithm (GA) coupled with
surrogate modelling. This chapter starts with fundamentals of natural
optimization, with focus on real, single objective GA. Fundamentals of surrogate
modelling and the construction of surrogate model with Kriging method follow in
next subchapters. Finally, explanation about coupling between GA and surrogate
modeling closes this chapter.

III.2 Fundamentals of Natural Optimization

Optimization is defined as the process of adjusting a set of input variables to a
mathematical process or function in order to find maximum or minimum output or
result (Haupt & Haupt, 2004). The input of optimization algorithm is a set of
variables, the process is called “fitness function”, and the output is called
“fitness”.

Among many methods to do optimization, analytical optimization is probably the

most popular one. Based on calculus, the maximum or minimum of a function can
be obtained by equating its first derivative (the gradient) with zero. Unfortunately,
many real-life problems, even more so in engineering, are complex and not
straight forward. Sometimes the fitness itself is a function of many variables;
while in other cases, the fitness function is not easily obtained or explicit or
continuous. Another problem arises as the algorithm of analytical method tends to
move towards local optimum solution instead of global one. An alternative
approach called natural optimization method has been developed to overcome this
problem.

Natural optimization method differs from analytical optimization method in the

sense that, instead of finding the first derivative of a fitness function and heading

18
towards possibly local optimum solution, natural optimization constructs a large
finite search space and statistically moves points in the space towards more
optimal places by employing intelligent, nature-based search algorithm. Using this
approach, finding the first derivative, or even the fitness function itself, is no
longer necessary and the existence of discrete fitness function is no longer a
problem. One of the examples of natural optimization, Genetic Algorithm (GA),
was used as optimization algorithm in this thesis. A brief explanation about GA
will follow in the next subchapter.

III.3 Single-Objective Real Genetic Algorithm (SOGA)

Genetic Algorithm (GA) is an optimization and search technique based on the
principle of genetics and natural selection. In GA, a set of variables is seen as an
individual composing a population, which is allowed to evolve, according to the
principles of genetic recombination and natural selection, under specified
selection rules to obtain maximum or minimum “fitness”. GA was first developed
by John Holland and his student, David Goldberg, in 1970s to solve gas-pipeline
transmission control problem (Haupt & Haupt, 2004).

Based on the objectives, GA is divided into two types: single-objective and multi-
objectives. In single-objective GA (SOGA), the objective of GA is to find the
maximum or minimum fitness of a single fitness function; while in multi-
objectives GA, the objective is to find the maximum or minimum fitnesses of
several fitness functions from similar input variables. Since the scope of this
thesis is limited to single-objective optimization, the process of SOGA will be
explained shortly below.

Based on how the input variables are formulated, GA is divided into two types:
binary and real GA. Binary representation, used in the earlier version of GA,
represents each individual as string of binary numbers that needs to be decoded
before GA process. This method is known for its simplicity and similarity to the
workings of chromosomes in the nature. However, for continous problems,
problems with many variables, or problems that require high precision, it is more

19
convenient to represent the each individual with floating-point numbers (Haupt &
Haupt, 2004). This method is called real GA (RGA).

The process of RGA can be divided into five steps, as shown by flowchart in
Figure III.1(a). Similar to the process of natural selection, a loop in GA to find the
maximum or minimum of fitness function f(x), where input variable x = {x1, x2,
x3, …, xN}T, consists of:
(1) Initial population generation
GA begins with generating a random Npop x Nvar matrix, where Nvar is the
number of variables and Npop is the number of individuals inside the
population on each generation. For the sake of convenience, all variables
are normalized into real numbers xnorm between [0, 1]. The unnormalized
variable x is obtained by the following equation:
x  xl  xh  xl xnorm
(III.1)
where xl and xh are the lowest and highest number in the variable range, or
the lower and upper bound, respectively.
(2) Fitness function evaluation
After initial population is generated, all individuals inside the population
matrix are evaluated using known fitness function f(x) to obtain the
“fitness”. If the fitness function is not known, the evaluation can be done
by creating a solver from numerical simulations or experiments. For the
sake of convenience, the objective of SOGA in this thesis is to find the
minimum fitness. However, when maximum fitness is the objective, the
algorithm can be adjusted by simply multiplying the fitness with -1.
(3) Selection
Selection is the process where individuals from a population are selected,
in a similar way to Darwinian natural selection, to form a mating pool.
There are several known methods to perform selection, such as: top-to-
bottom pairing, random pairing, weighted pairing (also known as roulette
wheel method), cost pairing, and tournament selection (Haupt & Haupt,
2004). In this thesis, tournament selection is chosen to be selection

20
 method. According to this approach, two individuals are randomly
selected from population and then undergo “mating competition”. The
individual with lower fitness will enter the mating pool. Selection process
will produce an Npop x Nvar mating pool matrix for crossover.
(4) Crossover
Similar to reproduction, crossover is the process in which two individuals
from mating pool (the parents) swap or blend their “genes” to create two
new individuals (the offsprings). As with selection, there are many
methods to perform crossover. In this thesis, the method for crossover is
SBX (Simulated Binary Crossover). SBX is an adaptive search method
where the crossover probability of distribution depends on the how close
the offsprings are to their parents (Deb & Kumar, 1995). As the solution
converges, the offsprings become more and more similar to the parents.
(5) Mutation
As found in the nature, mutation in GA is a genetic alteration introduced
on the offspings. However, unlike the nature, mutation in GA is needed to
prevent it from converging too quickly to local optimum. In GA, random
mutation on the value of offsprings occurs constantly according to a small
mutation probability.

A single loop in GA ends after mutation, in which an Npop x Nvar offspring matrix
is generated. This matrix will serve as the new population for next generation, and
so on. The iteration will end when the optimum solution converges or the number
of maximum generation is reached.

III.4 Surrogate – Assisted Genetic Algorithm (SAGA)

During the process of engineering design, decision making with regard to design
outcome must be done through analysis. The most common method for analysis in
engineering design is either experiment or numerical simulation. However, the
number of data available for analysis is very much limited to the computing power
or experimental capability. Surrogate or meta model acts as curve fit to available
data so that the result or outcome of an expensive simulation code may be

21
predicted. In addition to that, surrogate model can also provide calibration for
codes with limited accuracy, deal with noisy or missing data, and perform data
mining to obtain functional relationships between available design variables and
the result of interest (Forrester, Sobester, & Keane, 2008).

The core principal of surrogate construction is to construct y = f(x), where f is the

function that maps input vector x to scalar output y. As f is unknown, the solution
is to find a number of output values Y = {y(1), y(2), …, y(n)}T that correlate to input
vector X = {x(1), x(2), …, x(n)}T and construct the surrogate model f̂ x  for f based
on this observation (Forrester, Sobester, & Keane, 2008).

Due to its advantages, surrogate model can be included in GA codes. The purpose
of this is to increase the searching efficiency of GA, especially for complex
optimization problem (Le, Ong, Menzel, Jin, & Sendhoff, 2013). Coupling
between GA and surrogate model is known as Surrogate – Assisted Genetic
Algorithm (SAGA). Based on how the surrogate model is applied in GA, there are
two approaches of SAGA: local and global. In global SAGA (G-SAGA), the
surrogate model f̂ x  replaces the “true” fitness function f(x) for optimization;
while in local SAGA (L-SAGA), the surrogate model is constructed for every
individual in the population. It is found in the work of Dwianto (2015) that while
G-SAGA can reduce the number of fitness function evaluation significantly
compared to L-SAGA, this method becomes inaccurate for fitness function with
complex surface (e.g. surface with many local minima or maxima) and problem
with high dimensional optimization parameters.

From previous researches about flapping kinematics effect on force generation, it

is found that flapping plate kinematics optimization problem similar to this thesis
has a simple surface. Based on this finding and the scope of research, in which the
dimension of optimization parameters is two, the SAGA method chosen for
flapping plate kinematics optimization was G-SAGA with Kriging method for
surrogate model construction. Details about Kriging method and its coupling with
GA are explained in the next subchapters.

22
Figure III.1 Flowchart of (a) basic real, single objective Genetic Algorithm
(GA) (b) Global Surrogate – Assisted GA (G-SAGA) with Kriging
method

III.5 Kriging
Consider a set of n noise-free sample data, X = {x(1), x(2), …, x(n)}T, which
correlates to responses Y = {y(1), y(2), …, y(n)}T. The relation between sample set
and the responses can be denoted using a set of random vectors
 
 Y x 1

Y  
  
 Y x n  

(III.2)
Random vectors Y are connected with each other by a basis function

23
  k pj 
    
cor Y x i  , Y x l   exp    θ j x ji   x jl  
 j 1 
(III.3)
The right hand side term is a particular basis function for a surrogate model
construction method called Kriging. In Kriging method, the basis function ψ is in
the form of
 k 
 t   exp    θ j x jt   x j
pj


 j 1 
(III.4)
As shown by Equation (III.4), basis function in Kriging method has a varying

basis vector θ = {θ1, θ2, …, θk}T and exponent pj = {p1, p2, …, pk}T, pj ∈ [1, 2], for

each x. The objective of Kriging method is to construct surrogate model by

solving the equations for basis vector θ and exponent vector pj.

The construction of surrogate model by Kriging method begins with maximization

of likelihood function L, where L is expressed in natural logarithm form

n n
 
1
ln L    ln 2   ln   ln Ψ 
2 y  1  Ψ 1 y  1 
T

2 2 2 2
(III.5)
where σ2 is the variance and 1μ is the mean of random vectors Y (1 is a 1 x n
column vector consists of ones). Ψ is the so-called Gram matrix, which is defined
as
    
 cor Y x 1 , Y x 1

    
 cor Y x 1 , Y x n 
Ψ    
     
 cor Y x n  , Y x 1     
 cor Y x , Y x n 
n  

(III.6)
By finding the first derivative of Equation (III.5) and equating it with zero, the
MLEs (maximum likelihood estimates) for variable σ2 and is μ obtained,

24
 1T Ψ 1 y
ˆ 
1T Ψ 1 1
ˆ 2 
y  1 T Ψ 1 y  1 
n
(III.7)
Substituting the MLEs to Equation (III.5) and removing constant terms from the
same equation, a simplified function, called the “concentrated ln-likelihood
function”, is obtained
n
 
1
ln L    ln  2  ln Ψ
2 2
(III.8)
Through a lengthy derivation of likelihood function L, as elaborated by Forrester,
et. al. (2008), the surrogate model of f(x) is
fˆ x  ˆ  ψ T Ψ 1 y  1ˆ 
(III.9)
where ψ is the basis function vector, ψ = {ψ (1), ψ (2), …, ψ (n)}T. Equation (III.9)
shows that surrogate model f̂ x  is actually a function of basis vector θ and
exponent vector pj inside ψ and Ψ. This method is called Krigging interpolation.

Kriging interpolation, as mentioned previously, is designed for noise-free sample

data. In this thesis, however, random noises in GA fitness existed from force
measurement during towing tank experiment. To filter noisy sample data, a
regression constant λ is added to surrogate model in Equation (III.9) so that it
becomes
fˆ x  ˆ r  ψ T Ψ  I  y  1ˆ r 
1

(III.10)
where
1T Ψ  I 1 y
ˆ r 
1T Ψ  I 11
(III.11)

25
In this thesis, the exponent vector was set to a constant value, pj = 2, while the

basis vector θ was set in the range of θ ∈ [1 x 10-3, 1 x 102]. The value of

regression constant was set in the range of λ ∈ [1 x 10-6, 1].

III.6 Global Surrogate – Assisted Genetic Algorithm with Kriging

The coupling between GA and Kriging begins with population generation. The
population consists of Npop samples generated by sampling plan algorithm called
Halton sequence (Dwianto, 2015). Each individual inside population is then
evaluated using the “true” fitness function f(x) to obtain responses Y = {y(1), y(2),
…, y(Npop)}T that correlates to sample set X = {x(1), x(2), …, x(Npop)}T. Applying
Kriging method to the result of fitness evaluation, the surrogate model for true
fitness function f̂ x  is obtained. Surrogate model f̂ x  is optimized using RGA,
generating a set of optimum solution xopt. The optimum solution xopt is finally
evaluated using the “true” fitness function, generating the true, best fitness value
yopt = f(xopt). Flowchart of G-SAGA algorithm is shown in Figure III.1(b). This
thesis used G-SAGA code created by Dwianto (2015) for single objective case
with few adjustments to connect the algorithm with towing tank experimental
system.

26
 Chapter IV Experimental System Design and Setup

IV.1 Introduction
This chapter contains explanation about towing tank experimental system design
and setup. The first two subchapters explain towing tank experiment in general
and plate model that was used in the experiment. Since the plate must mimic 2-D
hovering insect kinematics, the next subchapters explain about 2-DOF robotic
system, which consists of translational (downstroke and upstroke) and rotational
(pronation and supination) mechanism. Some electrical and instrumentation
problems encountered during design process were handled and explained in these
subchapters. To measure force for optimization algorithm, force measurement
system was designed and explained in this chapter. Finally, explanation about PIV
and image acquisition system design for flow visualization around plate model
closes this chapter.

IV.2 Towing Tank Experiment

Figure IV.1 Towing tank in Laboratory of Aerogasdynamics, Department of

Aeronautics and Astronautics, Institute of Technology Bandung

Flapping plate optimization experiment was conducted in 100 cm x 60 cm x 60

cm glass-walled towing tank (Figure IV.1) in Laboratory of Aerogasdynamics,
Department of Aeronautics and Astronautics, Institute of Technology Bandung.
The tank was filled up with water (density ρ = 1000 kg/m3 and kinematic viscosity
ν = 1 x 10-6 m2/s). A metal carriage was attached to a threaded shaft above towing
tank to carry plate model, a servo for rotational plate motion, and a load cell for

27
Figure IV.2 Schematic of towing tank experimental system

28
aerodynamic force measurement. The threaded shaft was connected to a 3-phase
industrial AC motor for translational plate motion (see Figure IV.2).

As shown in Figure IV.2, the system for optimization can be divided into 3 parts:
translational motion, rotational motion, and force measurement system. For flow
visualization using PIV, the system can also be divided into 3 parts: translational
motion, rotational motion, and image acquisition system. All of the systems
mentioned above will be explained in details in the following subchapters.

IV.3 Flapping plate model

Flapping plate model was created from a 50 cm x 4 cm rectangular acrylic plate
with 4 mm thickness. The plate chord is equal to the width of the plate (4 cm) and
the span is equal to the typical water height inside the tank (43 cm). The aspect
ratio AR of the plate can be calculated using the following equation,
b2 b
AR  
S c
(IV.1)
Hence, the typical plate aspect ratio for this experiment is 10.75.

To ensure the nature of two-dimensional flow in this experiment, the plate was
placed within a very small gap from the bottom of the tank (see Figure IV.2). This
condition is called “grazing” (Lisoski, 1993).

IV.4 Translational System

As shown by Figure IV.2, translational system in this towing tank experiment
consists of a 3-phase industrial AC motor, an inverter, an optical encoder for
system feedback, and a PCB containing relay, logic gate for limit switch, and an
optocoupler IC for ground loop prevention. The system was designed to move the
plate back and forth within a straight trajectory, mimicking downstroke and
upstroke translation of an insect wing. To accomplish that motion, the plate
carriage was attached to a threaded shaft with 2 mm pitch. This shaft was
connected to 3-phase industrial AC motor. The correlation between rotational

29
speed (in rpm) of an ideal, no-load AC motor with its frequency f (in Hz) and
number of pole n is
120
rpm  f
n
(IV.2)
The 3-phase industrial AC motor used in this experiment has 4 poles with the
maximum rotational speed of 1450/1720 rpm for 50/60 Hz motor frequency.
Hence, the translational velocity of the carriage v (in cm/s) can be defined as
rpm
v d
60
(IV.3)
where d is the pitch of the threaded shaft, d = 2 mm (0.2 cm).

To control motor frequency and the direction of translation (or actually the
direction of shaft rotation), a Toshiba VF S-11 industrial inverter was connected
to the AC motor. Using this inverter, a direct, linear correlation between input
voltage and motor frequency can be obtained (60 Hz from 10 Vdc voltage input)
and the direction of shaft rotation can be set from F and R pin of the inverter,
where F and R denote the forward direction (shaft rotating clockwise) and reverse
direction (shaft rotating counter-clockwise), respectively (see Figure IV.2).

One of the objectives of experimental system design was to create an automatic

system. This means that, all parts of the system must be connected to a computer
and controlled using a software; in this case the software is MATLAB R2010a. In
translational system, an Arduino Uno R3 board was used to connect AC motor
and inverter with the computer. Arduino is essentially a microcontroller board
based on ATmega328, an AVR 8-bit processor. Arduino Uno has 14 digital I/O
pins, of which 6 pins can be used as PWM outputs, and 6 analog input pins. The
frequency of AC motor is controlled from VIA pin on inverter panel. Arduino
converts an 8-bit integer PWM output to input voltage for the inverter. Since the
operating voltage of Arduino Uno is limited to 5V, the maximum motor frequency
is ~30 Hz, which translates into ~3 cm/s maximum translational velocity.

30
To control the direction of translation, a relay was used. The relay works in such
way that the carriage moves in forward direction if F pin is shorted with CC pin,
reverse direction if R pin is shorted with CC pin on inverter panel. Arduino Uno
controls the relay using digital output 0 or 1 for F or R direction, respectively.

IV.4.1 Two-Stage PID Controller

To ensure that the translational velocity output from AC motor is equal to the
input or reference velocity, a feedback system is required. PID controller, the most
common control system with feedback for industrial application (Astrom &
Murray, 2008), was designed and coded in MATLAB R2010a and Arduino IDE.
In this thesis, PID controller for AC motor is split into 2 stages: during transient
and steady state response, where PID constants are automatically changed.

Figure IV.3 Two-stage PID system for AC motor control

As indicated by its name, PID controller consists of 3 parts: proportional (P),

integral (I), and derivative (D). In this experimental system, as shown by a
simplified block diagram in Figure IV.3, a reference velocity vref acts as the input
of the system (in experiments, vref is equal to translational plate velocity U). PID
system acts as compensator for AC motor response. This output is measured as
vread and fed back to the system by an optical encoder attached to the end of the
shaft (see Figure IV.3). Error e between reference and output velocity serves as
the input for PID transfer function, where the output velocity from PID controller
to the AC motor is
t
de
vout t   K p et   K i  et dt  K d
0
dt

(IV.4)

31
where Kp, Ki, and Kd are proportional, integral, and derivative constant,
respectively.

The towing tank experiment used Autonics E50S8 rotary encoder with 1024 P/R
resolution. This is an incremental rotary encoder that converts angular position or
motion of the shaft to number of pulses p. According to its specifications, the
encoder sends 1024 pulse per shaft revolution to Arduino. Hence, AC motor
response can be converted to measured translational velocity vread using equation
d  pt  pt 1 
vread  
t 1024
(IV.5)

All three constants of PID in equation (4.4) will affect system response.
Proportional constant Kp acts as response amplifier with adjustable gain (Ogata,
1997). As integral constant increases, steady-state error will decrease, but the
response has higher tendency towards oscillation. The derivative constant acts as
damper to oscillating response up to a certain value, then it becomes ineffective
(Astrom & Murray, 2008). Also, the transient response will improve as the Kd
increases.

There are several ways to determine the value of PID constants for a system, for
example: Ziegler-Nichols method, determining critical gain experimentally, et
cetera (Ogata, 1997). In towing tank experiments, however, the response of AC
motor itself was irregular and unpredictable at most times, therefore could not be
modeled. Another approach chosen to determine the value of PID constants was
fine tuning. This means that the value of Kp, Ki, and Kd are determined by trial-
and-error. It was found during fine tuning process that increasing Ki would create
oscillation in transient response, which translated into rapid acceleration-
deceleration in motor speed, alternatingly. This result is undesireable since the
translational system is required to reach reference velocity as rapid as possible and
then moves at constant velocity. Therefore, as the result of fine tuning, Ki was set
to zero (the system became PD) and the value of Kp and Kd are shown in Table
IV.1.

32
 Table IV.1 Two-stage PID constants
PID constants First stage Second stage
Kp 0.20 0.040
Kd 0.03 0.003
Ki 0 0

Another control engineering problem that arose during system design was
irregularities in steady-state response. Due to the compensating nature of PID
controller, sometimes obstacles along carriage path caused motor speed
deceleration, which would be followed by sudden motor acceleration from PID
response, creating “spikes” along steady-state response of the motor. This
response is, of course, undesireable. Therefore, a two-stage PID controller was
designed. During transient response, the controller works with prescribed PID
constants. When velocity error is less than 10% of its reference, Kp is dropped to
20% while both Ki and Kd are dropped to 10% of their prescribed values (Table
IV.1). With this design, obstacles or disturbance along carriage path will not cause
sudden acceleration or “spikes” in PID response.

IV.4.2 Limit Switch and Logic Gate

Another issue that also arose while conducting this research was error or
inaccuracy in control system. In several occasions, this error or inaccuracy could
cause the motor to rotate until the carriage hit the edge of tank. Since this
condition is undesireable, two limit switches (denoted by SWF and SWR) were
mounted at both edges of tank (see Figure IV.2). These switches work in such
way that motor’s PWM signal from Arduino will be cut off when the carriage hit
either of them and the carriage will move in reverse direction of the present relay
condition whenever the relay gives signal to do so.

To accomplish those tasks, a logic gate was designed using Karnaugh map, or K-
map for short. The map has three inputs: relay, forward switch (SWF), and reverse
switch (SWR) state. The state of relay is either F (0) or R (1). Since both switches
were normally closed type, the state is either closed (switch is not pressed, 1) or
open (switch is pressed, 0). The statement is true for the following conditions:
(1) Both switches are closed and the carriage moves in either F or R direction;

33
 (2) SWF is open and the carriage moves in F direction;
(3) SWR is open and the carriage moves in R direction.

Table IV.2 Karnaugh map for limit switch logic function

SWR, SWF
00 01 11 10
Relay 0 (F) 0 1 1 0
1 (R) 0 0 1 1

The K-map (Karnaugh map) for logic function above is shown in Table IV.2.
Using this table, a logic function for towing tank limit switch case is also
obtained,
out  SWR  R  SWF  F
(IV.6)

Figure IV.4 (a) Logic gate for limit switches and relay (b) PCB containing
relay, limit switches logic gate, and voltage coupler

The function can be translated into logic gate in Figure IV.4(a). The gate consists
of AND, OR, and NOT gates. A PCB was designed for limit switch logic gate
with 4049 hex inverting buffer IC for NOT gate, 4071 OR gate IC, and 4073
AND gate IC (see Figure IV.4(b)).

IV.4.3 Ground Loop

Ground loop was another instrumentation problem that occured during towing
tank experiments. Ground loop occurs when devices that are supposed to be in the
same potential, usually common ground, are actually at different potentials. In

34
towing tank experiments, ground loop occurred whenever inverter was plugged
into the system, causing noises in force measurement.

Ground loop problem in this case can be overcome by simply separating inverter
and Kyowa CDV-700A signal conditioner’s power supply, which are connected
to the load cell. However, noises from ground loop were once again apparent
since inverter was connected to computer via Arduino.

Figure IV.5 Electrical schematic of voltage coupler

To overcome ground loop problem, a PCB, shown in Figure IV.5, was designed to
isolate inverter from the rest of the system, including the ground. An optocoupler
IC, 4N35, acts as voltage coupler between Arduino board and inverter. A 4N35
optocoupler IC consists of two essential parts: an infrared LED and NPN
phototransistor (see Figure IV.5). LED’s brightness changes according to PWM
signal’s voltage from Arduino. NPN phototransistor acts as voltage switch, with
operating voltage between ground and PP10 pin (10 Vdc) on inverter panel.
However, NPN transistor switches PWM input from HIGH to LOW and vice
versa; hence 4N35 IC must be connected to MOSFET to flip the logic. Finally, the
output voltage ranges from 0 to ~5 Vdc through VIA pin on inverter panel. This
way, the inverter is isolated from the rest of the system; hence the noise during
force measurement due to ground loop (~60 mV) is no longer apparent.

35
IV.4.4 Encoder Pulse Divider
As explained previously in subchapter IV.4.1, an incremental roraty encoder acts
as translational system feedback by converting angular position of shaft to number
of pulses. For this experiment, encoder must be able to recognize whether the
carriage moves in forward (shaft rotating clockwise, pulse increasing), or reverse
direction (shaft rotating counter-clockwise, pulse decreasing).

An incremental rotary encoder used in this experiment is designed with six output
phases A, A, B, B, and Z, , also known as index line. For this experiment, only
outputs from line A and B were used. Waveform output from line A and B has
similar amplitude and period, but different phase (see Figure IV.6). The direction
of rotation can be sensed by comparing the state of waveform output between
these two lines using ISR (Interrupt Service Routine) in Arduino.

Figure IV.6 Encoder waveform output (Autonics Corporation, 2015)

For encoder application, ISR is programmed to interrupt Arduino serial

communication whenever it senses voltage changes (rise, drop, or both) in line A.
Once the interrupt is initiated, state between line A and B is compared to
determine the direction of rotation, whether the number of pulses is increasing or
decreasing.

36
For this experiment, the system is designed to connect input from MATLAB to
Arduino via serial communication, also an ISR. Unfortunately, ATmega328 is
unable to process serial communication and encoder ISR simultaneously. As an
illustration, a 1024 P/R rotary encoder sends roughly 15 pulse/ms if the shaft is
rotated at 900 rpm, which translates into maximum translational velocity 3 cm/s.
Servo, on the other hand, requires 1-2 ms to set angle between 0-180º and some
more time to receive command from serial communication. From this illustration,
we can conclude that while setting servo to a particular angle, encoder ISR will be
skipped, along with at least 15-30 pulses from encoder, causing inaccuracy in
translational motion.

To overcome this problem, a pulse divider was designed to divide encoder

resolution by 8, hence dropping the resolution from 1024 P/R to 128 P/R, before
the encoder sends output to Arduino. Pulse divider was designed using another
Arduino Uno R3 board (see Figure IV.2). The divider works in such way that it
receives pulse from pin A and pin B, which are attached to interrupt pin 2 and 3 in
Arduino, and repeatedly sends pulse to main Arduino only after the counter hits 8.
To shorten computing cycle, pulse divider board was programmed in Arduino
IDE using DPM (direct pin manipulation) and bitwise operation. It was also found
to be a good practice to attach pull-down resistors to encoder output cables.

IV.5 Rotational System

To mimic pronation and supination, a rotational sytem was designed to rotate the
plate and load cell. This motion is accomplished using Futaba S3003 servo,
powered by 5V pin on Arduino (Figure IV.2). Servo sets angle from 0-180º by
responding to pulse length modulation with short duration between 1-2 ms. In this
experiment, servo is controlled using <Servo.h> library in Arduino.

IV.6 Force Measurement System

To measure force generated by plate during flapping motion, force measurement
system was designed. The system consists of Scaime BE 1 kg single point load
cell, Kyowa CDV-700A signal conditioner as signal amplifier and low pass filter,

37
and National Instruments DAQ (data acquisition) USB-6211 as ADC (analog to
digital converter), which transmits analog signal from force measurement to
computer.

Figure IV.7 Decomposition of normal force, measured by load cell

The objectives of the system is to measure aerodynamic forces (lift and drag)
generated by plate in flapping motion. Plate model is attached perpendicularly to a
binocular, single point load cell to measure normal force experienced by plate,
denoted by Fn (see Figure IV.7). This normal force consists of inertial force from
plate’s acceleration and aerodynamic force from fluid interaction around the plate.
Normal force is decomposed into lift L and drag D according to convention in
Figure IV.7 using equation
L  Fn cos 
D  Fn sin 
(IV.7)
Lift and drag coefficient are defined as
2L
Cl 
U 2 S
2D
Cd 
U 2 S
(IV.8)

Scaime BE 1 kg is a binocular, single point load cell that contains a strain gauge
as force sensor. The gauge is arranged to form a quarter-bridge strain gauge
circuit, where it acts as the active arm, depicted as a variable resistor in Figure

38
IV.2. Tension and compression experienced by the gauge changes gauge’s
resistance and changes Wheatstone bridge’s output voltage, simultaneously.
Before conducting the experiment, load cell was calibrated using a set of
calibration weights to obtain linear relationship between output voltage and
weight (gravitational force).

As with any other measurement system, noise also occurs during flapping force
measurement process. In the towing tank experiment, this is a problem since the
force measured from flapping plate is very small, in the order of O(102) mN,
compared to load cell’s capacity. Noise problem can be overcome by setting the
SNR (signal-to-noise ratio) to the highest possible value and/or filtering the data
obtained from measurement. This can be achieved by increasing bridge excitation
voltage, increasing gain during signal processing, and creating lowpass filter.
During the experiment, Kyowa CDV-700A signal conditioner was operated in 4V
bridge excitation voltage, highest possible gain (or smallest RANGE in the value
of 100), and 10 Hz cutoff frequency for analog low pass filter. Using this setting,
bridge output voltage was typically measured in the order of O(101) -O(102) mV.
Another noise source around 60 mV from ground loop was diminished by
isolating the inverter, as explained in subchapter IV.4.3.

IV.7 Particle Image Velocimetry

Flow measurement and visualization was the next step of this research after
optimization and towing tank experiments produced optimized plate kinematics.
The purpose of flow measurement and visualization is to obtain quantitative
values of flow properties and a clear view of flow physics around a flapping plate
during hovering flight, to be analyzed and studied further.

Due to the complexity of fluid flow around a flapping plate, it is difficult to

capture the flow and measure properties of fluid using conventional experiment
method, e.g. wind tunnel test. A proposed experimental method to answer this
problem is PIV (Particle Image Velocimetry). PIV is an experimental method for
fluid mechanics research using optical principles to obtain fluid properties and

39
Figure IV.8 PIV system for towing tank experiment

40
flow visualization. Compared to well-known conventional approaches, PIV has
several advantages. First, it is non-intrusive. This means that PIV does not require
a sensor or other equipments that may disrupt the flow. A wind tunnel test usually
uses pitot tube to measure fluid pressure at one point only and needs other method
(e.g. dye flow visualization) to visualize fluid flow. PIV, on the other hand, is able
to measure fluid properties globally and capture flow visualization at the same
time.

Figure IV.9 Basic PIV scheme for flowfield measurement (Raffel, Willert,
Wereley, & Kompenhans, 2007)

The principals of PIV will be explained shortly as follows. Hollow spherical

particles are immersed in fluid that will be the object of the research (see Figure
IV.9). A laser is shot through concave and cylindrical lens to form a surface which
illuminates fluid particles. A CCD camera is used to capture two images of the
flow within a short time step. By measuring the displacement of each particle on
the field from two consecutive images using cross correlation algorithm, other
derivatives such as velocity vector and vorticity can be acquired. The final result
of PIV is a velocity field surrounding an object, observed at a certain time. Taking
several images within a time span, the evolution of velocity field (as a function of
time) can be observed.

In this research, an AVT Prosilica GX1050 mono CCD camera, along with Cosina
lens (AF 28-80 mm, F 3.5-5.6) and Nikon F-C lens adapter, was mounted under

41
the tank (see Figure IV.8) to capture 2-D chordwise flow around flapping plate. A
Viasho VA-II greenlight laser was shot through a cylindrical planar concave lens
(focal length -20 mm) to form an illuminated surface on plate’s midspan inside the
tank where Dantec Dynamics’ hollow glass sphere particle (10 μm diameter) was
immersed. Image acquisition was conducted using MATLAB 2010a, where the
camera was triggered and the images were acquired via Intel PCI-E ethernet card.
Frame rate and shutter speed were set in such way that particle elongation, a
condition where the particle appears in an image as white streaks instead of dots,
was prevented. These settings, along with other camera settings during image
acquisition are shown in Table IV.3.

Table IV.3 Prosilica GX1050 setting for image acquisition

Parameter Unit Value
Image type - Mono, 8-bit
Image resolution pixels 1024 x 1024
Frame rate fps 90
Shutter speed μs 11111
Gain - 2

Table IV.4 Image processing parameter configuration in PIV_GUI.m

Parameter Configuration
Grid size 32 pixels
IA size 32 pixels
SA size 32 pixels
Time step 0.0111 s
Cross correlation method FFT
Subpixel interpolation method Compensated Gaussian
Outlier test method Local median
Outlier removal 3x
Window shifting technique CDI
Window deformation technique Bilinear Interpolation

To obtain velocity vector and vorticity magnitude around plate model from
images acquired during image acquisition, image processing was conducted using
cross correlation algorithm in PIV_GUI.m created by Stepen (2011). This
algorithm is enhanced with advanced processing technique called window shifting
and window deformation to prevent in-plane loss of pair. The parameter
configuration of image processing in this thesis is shown in Table IV.4.

42
 Chapter V Experimental Optimization of Flapping Plate
Kinematics using G-SAGA with Kriging

V.1 Introduction
In this chapter lays the core of this thesis. It will be begun by a preliminary towing
tank experiment in order to validate force measurement system for fitness
calculation during optimization process. The optimization of flapping plate was
conducted within 2 test cases to maximize time-averaged lift coefficient in
hovering condition using G-SAGA with Kriging for surrogate modeling method.
In each test case, 2 kinematics parameters were optimized and their effect on lift
generation was also plotted and studied further. First test case was intended to
optimize plate rotational kinematics and study their effect on lift generation,
whereas the objective of the second test case was to optimize and study the effect
of stroke amplitude and angle of attack during hovering condition. Furthermore, it
is also desired to study the physics of flow under optimum condition. The flow
physics of optimum solutions and their comparison with suboptimum solutions
was measured and visualized using PIV (see subchapter IV.7). Finally, analysis
and discussion about the results of flapping plate optimization in this thesis will
close this chapter.

V.2 Preliminary Experiment and Validation

Before conducting optimization and experiments, it is necessary to validate force
measurement system, which serves as the fitness in optimization algorithm. The
validation was done by recreating an existing experiment and comparing the drag
coefficient obtained from those two similar experiments. The experiment chosen
to be the reference is an experiment to study nominally 2-D flow around a normal
flat plate, which was done by Lisoski (1993).

In his experiment, Lisoski accelerated a rectangular flat plate, angled 90º with
respect to direction of translation, from rest until it reached maximum velocity U
within the distance of 1.31 times the chord length. The flat plate then moved in
constant velocity U until it traveled 64 times the chord length from start. Lisoski

43
found that the drag peaked to the value of Cd = +3 during acceleration. To ensure
the nature of 2-D flow, the plate was placed in grazing condition from the bottom
of the tank (Lisoski, 1993). The experiment was conducted in several Re: 1000,
2500, and 5000.

Figure V.1 Schematic of the reproduction of Lisoski's normal flat plate

experiment (1993)

Similar experiment was conducted in towing tank in Laboratory of

Aerogasdynamics, Institute of Technology Bandung. An acrylic flat plate with 4
cm chord length was used for this experiment. Similar to Lisoski’s, the plate was
angled 90º with respect to direction of translation and accelerated from rest to
maximum velocity U = 2.5 cm/s (Re = 1000) within the distance of 1.31 times the
chord length (see Figure V.1). At the end of experiment, the plate was decelerated
to rest within the distance of 1.31 times the chord length. Due to the limitation in
towing tank’s length, the plate only traveled 14 times chord length from start until
it stopped, instead of 64. Like Lisoski’s, the experiment was repeated several
times (in this case, 24 times) and the drag coefficient result was ensembled and
averaged.

Comparison between present towing tank experiment and Lisoski’s normal flat
plate experiment is shown in Figure V.2(a). As expected, Cd peaked at the
beginning of the experiment, during acceleration phase. However, while the peak
in Lisoski’s experiment reaches 3.16, the peak value in present experiment is
about 2.9. This is the result of typical run, ensembled from 24 runs to ensure the
repeatability of the experiment and shown in Figure V.2(b). Acceleration phase,
characterized by Cd peak, is immediately followed by constant velocity motion,
where initial symmetric recirculating bubble occurs, corresponding to Cd drop to
approximately 1.3, followed by organized vortex shedding, raising the d to ~1.6.

44
Figure V.2 (a) Drag coefficient Cd comparison of present experiment to
Lisoski's experiment, (b) Ensemble average and individual run to
ensure repeatability of the experiment

The RMS error of time-averaged drag coefficient d, compared to the reference, is

10.77%. The error is caused, particularly, by load cell’s limited precision: the
capacity of load cell is 1 kg, but the measurement has the order of O(10 1) to
O(102) gr. However, the trend of Cd measured from towing tank experiment in
Figure V.2(a) is considered close enough to the reference. Therefore, the present
force measurement system is sufficient to measure aerodynamic force in towing
tank experiment for fitness calculation.

V.3 Experimental Optimization of Flapping Plate Kinematics using G-

SAGA with Kriging
As stated previously, the goals of this research are to find optimized 2-D
kinematics of a flapping plate model and to study the effect of various kinematics
to flow field, and consequently, aerodynamic force generation of flapping plate
during hover. Hovering is a particular condition where the forward velocity of
insect body is zero (reduced frequency k goes to infinity); hence the insect body is
supended in the air (Ellington, 1984), hence lift generated solely by flapping
during hover must be the highest since there is no additional lift from insect
body’s forward velocity. In this thesis, two test cases, each with different
kinematics parameters, were optimized for maximum time-averaged lift

45
coefficient. Brief explanation about the flowchart of flapping plate optimization
algorithm will be given below.

Figure V.3 Flowchat of flapping plate optimization algorithm using real,

single objective G-SAGA with Kriging, coupled with towing tank
experiment

Optimization was conducted using real, single objective G-SAGA with Kriging as
surrogate construction method. Optimization using G-SAGA began by generating
Npop number of population which consisted of a set of kinematics parameter x.
The algorithm was coupled with towing tank experiment for fitness evaluation
(see flowchart in Figure V.3). For flapping case, the objective was to find
maximum fitness, where fitness was the measured time-averaged aerodynamic lift
coefficient l(x) of plate model, whose motion was dictated by kinematics
parameter input x. The force was measured by load cell during towing tank
experiment. Kriging method was employed to construct surrogate model of fitness
function, which would be optimized using RGA. Since RGA was programmed to

46
find minimum fitness instead of maximum, it was necessary to multiply all fitness
values by -1. The optimum set of flapping kinematics xopt obtained from RGA was
studied further using PIV for flow visualization.

In this thesis, optimization was done for a 2-D flapping plate in hovering
condition. To model hovering condition, plate model is programmed using
Arduino Uno R3 to conduct 2-DOF motion: downstroke and upstroke translation,
pronation and supination (see Figure II.3). The stroke plane was horizontal (β =
0), as observed in the kinematics of flies, bees (Bombus terrestris), and wasps
(Encarsia formosa) (Wang, 2004). This simplification of flapping kinematics is
deemed sufficient since most insects flap their wings in a horizontal or inclined
stroke plane (Ellington, 1984). Furthermore, the 2-D flapping plate model is rigid;
therefore the passive pitching/ rotation phenomenon due to wing flexibility, which
is actually closer in mimicking real insects, is outside the scope of this thesis and
can be neglected.

Optimization was conducted in constant Reynolds number Re = 1200 (constant

translational velocity U = 3 cm/s), which was within the range of Reynolds
number of honeybees (Apis mellifera), hoverflies (Eristalis tenax), and hawker
dragonflies (Aeshna juncea) (Ellington, 1984). The fitness was calculated from
the first four strokes of each evaluation (N = 4). To prevent flow disturbances
inside the towing tank, a 30-seconds pause was given between each fitness
function evaluation. The number of population Npop was 100 for each G-SAGA
run; hence the total number of fitness evaluation is 100 plus 1 final evaluation on
optimum result.

V.3.1 Case 1: The Effect of Plate Rotation

It was suggested by Dickinson (1999) that at the end of each half-stroke, insects
rotate their wings rapidly to create flow circulation that generates lift.
Furthermore, Dickinson observed that lift enhancement during stroke reversal is
heavily affected by rotational timing (see II.5). Through experiments, it is found
that the wing must rotate prior to stroke reversal, a kinematics known as

47
“advanced rotation”, in order to enhance lift at the end of half-stroke. Hence,
previous researches suggested that rotation must be rapid and advanced in order to
enhance lift. Now the questions are how these two rotational parameters
simultaneously affect lift generation and which set of rotational parameters
generates maximum lift.

To answer the question, an optimization using G-SAGA with Kriging was done
for two kinematics parameters x = {tr, Δtr}T, where tr is the time when rotation
starts and Δtr is the duration of rotation, both are relative to stroke period T. The
objective of optimization was to maximize time-averaged lift coefficient. The path
of optimization is as follows:
tr , t r   0.5 0.5
Find
t r , t r  0 0.5
4T
1
to maximize Cl 
4T  c t dt
0
l

The minus sign in rotation time indicates that the rotation starts before downstroke
ends (signifying advanced rotation), whereas plus sign indicates that the rotation
starts after downstroke ends (signifying delayed rotation). When tr = 0, rotation
starts exactly at the end of downstroke phase. The optimization was conducted in
Re = 1200, stroke amplitude da = 4 times the chord length, and under symmetric
hovering condition: the downstroke and upstroke angle of attack of the plate
model were equal, α = 45º.

Figure V.4 shows typical flapping plate’s time-averaged lift coefficient surface
and contour plot of Case 1 as a function of rotation time tr and duration Δtr
constructed by Kriging, and also maximum time-averaged lift coefficient from
optimization using RGA. As shown by Figure V.4, the surface has two valleys:
one of global minimum at tr > 0 (delayed rotation) and one of global maximum at
tr < 0 (advanced rotation).

It is important to keep in mind that the result pictured in Figure V.4 is a typical
result. Due to the presence of random error in experimental system, different G-

48
Figure V.4 (a) Typical surface and (b) contour of time-averaged lift coefficient
of flapping plate of Case 1. The surface and contour shows the
effect of rotational parameters. Blue dot indicates maximum
fitness.

Figure V.5 Histogram of optimum set of rotational kinematics in Case 1,

ensembled from 25 G-SAGA run

SAGA optimization run may result in different optimum kinematics and different
maximum fitness time-averaged lift coefficient. Figure V.5 shows the histogram
of optimum set of kinematics parameter xopt = {tr, Δtr}T ensembled from 25 G-
SAGA run. Two sets of optimum kinematics parameter that have the highest bar,

x opt,1   0.25 0.3

T

x opt, 2   0.15 0.2

T

were chosen as the optimum solutions. Illustration of flapping kinematics at one

full stroke in shown in Figure V.6(a) and Figure V.7(a). Figure V.6(b) and Figure
V.7(b) show the force history of flapping plate during 4 strokes. To study the
physics of flow around plate model and how it affects lift generation, an image
sequence of flow field evolution, in terms of velocity vector and vorticity, during

49
Figure V.6 Experimental result of first optimum solution of Case 1: (a) an
illustration of flapping motion of plate model during one full
stroke, black dots indicate leading edge; (b) Force and angle of
attack history of 4 flapping strokes; (c) the evolution of fluid flow
during the first 2 strokes (black arrow indicates velocity vector and
colorbar depicts magnitude of vorticity ω). Due to perspective
error, flow properties at grey area near plate are unavailable,
especially when α ≠ 90 º

50
Figure V.6 (c) the evolution of fluid flow during the first 2 strokes (black
arrow indicates velocity vector and colorbar depicts magnitude of
vorticity ω). Due to perspective error, flow properties at grey area
near plate are unavailable, especially when α ≠ 90 º

51
Figure V.7 Experimental result of second optimum solution of Case 1: (a) an
illustration of flapping motion of plate model during one full
stroke, black dots indicate leading edge; (b) Force and angle of
attack history of 4 flapping strokes; (c) the evolution of fluid flow
during the first 2 strokes (black arrow indicates velocity vector and
colorbar depicts magnitude of vorticity ω). Due to perspective
error, flow properties at grey area near plate are unavailable,
especially when α ≠ 90 º

52
Figure V.7 (c) the evolution of fluid flow during the first 2 strokes (black
arrow indicates velocity vector and colorbar depicts magnitude of
vorticity ω). Due to perspective error, flow properties at grey area
near plate are unavailable, especially when α ≠ 90 º

53
Figure V.8 Experimental result of suboptimum solution of Case 1: (a) an
illustration of flapping motion of plate model during one full
stroke, black dots indicate leading edge; (b) Force and angle of
attack history of 4 flapping strokes; (c) the evolution of fluid flow
during the first 2 strokes (black arrow indicates velocity vector and
colorbar depicts magnitude of vorticity ω). Due to perspective
error, flow properties at grey area near plate are unavailable,
especially when α ≠ 90 º

54
Figure V.8 (c) the evolution of fluid flow during the first 2 strokes (black
arrow indicates velocity vector and colorbar depicts magnitude of
vorticity ω). Due to perspective error, flow properties at grey area
near plate are unavailable, especially when α ≠ 90 º

55
the first 2 strokes is shown in Figure V.6(c) and Figure V.7(c). Detailed analysis
and further discussion on these results will be presented in subchapter V.4. Note
that the result from PIV experiments contains error due to perspective error and
mechanical vibration in the plate’s and camera’s carriage. Wake cannot be shown
clearly because the chord is too long for camera’s ROI. This is the consequence of
maintaining Re to the value of 1200.

To obtain clear comparison between optimum solution and the rest of the
solutions in the search space, one suboptimum solution was picked from the
search space by simply delaying the rotation to the very end of downstroke phase
(or at the beginning of upstroke that follows it),
x subopt  0 0.3
T

Illustrations of flapping kinematics at one full stroke, force history of flapping

plate during 4 strokes, and the evolution of flow field during flapping are shown
in Figure V.8. Comparison between suboptimum and optimum result will be
discussed in detail in subchapter V.4.

V.3.2 Case 2: The Effect of Stroke Amplitude and Angle of Attack

Previous researches by Sane (2001), Trizila (2008), and Trizila (2008) regarding
the effect of various flapping kinematics on aerodynamic force generation show
that angle of attack affects the magnitude of flapping force greatly. Surprisingly,
all researches mentioned above concluded that the effect of stroke or plunging
amplitude is insignificant. Similar to their researches, the effect of angle of attack
and stroke amplitude is also studied in this thesis with the addition of G-SAGA
optimization to discover optimum flapping plate kinematics.

G-SAGA with Kriging was done for two kinematics parameters x = {da, α}T,
where da is the stroke amplitude, or the distance traveled by plate model during
downstroke or upstroke, relative to chord length c and was varied from 3 to 6
times the chord length. The optimization was done under symmetric hovering
condition, where angle of attack α, which ranges from 0-90 º, during downstroke

56
and upstroke, is equal. The objective of optimization was to maximize time-
averaged lift coefficient. The path of optimization is as follows:
da , d a  3 6
Find
 ,  0 90
4T
1
to maximize Cl 
4T  c t dt
0
l

Other parameters were kept constant during optimization: Reynolds number of the
plate Re = 1200 and the rotation was symmetrical (tr = -0.1, Δtr = 0.2).

Figure V.9 (a) Typical surface and (b) contour of time-averaged lift coefficient
of flapping plate of Case 2. The surface and contour shows the
effect of stroke amplitude and angle of attack. Blue dot indicates
maximum fitness.

Figure V.10 Histogram of optimum set of flapping kinematics in Case 2,

ensembled from 25 G-SAGA run

Figure V.9 shows typical flapping plate’s time-averaged lift coefficient surface
and contour plot of Case 2 as a function of stroke amplitude da and angle of attack

57
α constructed by Kriging, and also maximum time-averaged lift coefficient from
optimization using RGA. As shown by Figure V.9, the surface has one global
optimum at high angle of attack.

As with Case 1, this optimization problem also has random error. Figure V.10
shows the histogram of optimum set of kinematics parameter xopt = {da, α}T
ensembled from 25 G-SAGA run. Two sets of optimum kinematics parameter that
have the highest bar,

x opt,1  3 40
T

x opt, 2  3 46
T

were chosen as the optimum solutions. It is interesting to note that despite the
presence of random error, all G-SAGA run produced da = 3 as optimum solution
because of reason that will discussed later in subchapter V.4. Similar to Case 1,
illustration of flapping kinematics at one full stroke, force history of flapping plate
during 4 strokes, and image sequence of flow field evolution, in terms of velocity
vector and vorticity, during the first 2 strokes are presented in Figure V.11and
Figure V.12. One suboptimum solution was picked from the search space by
increasing the angle of attack,
x subopt  3 80
T

Illustrations of flapping kinematics at one full stroke, force history of flapping

plate during 4 strokes, and the evolution of flow field during flapping are shown
in Figure V.13. Detailed analysis and further discussion on both optimum and
suboptimum result will be presented in subchapter V.4.

V.4 Result Analysis and Discussion

As shown in Figure V.4 and Figure V.9, the G-SAGA with Kriging code has
successfully predicted lift of 2-D flapping plate as a function of flapping
kinematics and optimized the kinematics to obtain maximum lift. Whether the
results are acceptable or not will be discussed in the following subchapters.

58
Figure V.11 Experimental result of first optimum solution of Case 2: (a) an
illustration of flapping motion of plate model during one full
stroke, black dots indicate leading edge; (b) Force and angle of
attack history of 4 flapping strokes; (c) the evolution of fluid flow
during the first 2 strokes (black arrow indicates velocity vector and
colorbar depicts magnitude of vorticity ω). Due to perspective
error, flow properties at grey area near plate are unavailable,
especially when α ≠ 90 º

59
Figure V.11 (c) the evolution of fluid flow during the first 2 strokes (black
arrow indicates velocity vector and colorbar depicts magnitude of
vorticity ω). Due to perspective error, flow properties at grey area
near plate are unavailable, especially when α ≠ 90 º

60
Figure V.12 Experimental result of second optimum solution of Case 2: (a) an
illustration of flapping motion of plate model during one full
stroke, black dots indicate leading edge; (b) Force and angle of
attack history of 4 flapping strokes; (c) the evolution of fluid flow
during the first 2 strokes (black arrow indicates velocity vector and
colorbar depicts magnitude of vorticity ω). Due to perspective
error, flow properties at grey area near plate are unavailable,
especially when α ≠ 90 º

61
Figure V.12 (c) the evolution of fluid flow during the first 2 strokes (black
arrow indicates velocity vector and colorbar depicts magnitude of
vorticity ω). Due to perspective error, flow properties at grey area
near plate are unavailable, especially when α ≠ 90 º

62
Figure V.13 Experimental result of suboptimum solution of Case 2: (a) an
illustration of flapping motion of plate model during one full
stroke, black dots indicate leading edge; (b) Force and angle of
attack history of 4 flapping strokes; (c) the evolution of fluid flow
during the first 2 strokes (black arrow indicates velocity vector and
colorbar depicts magnitude of vorticity ω). Due to perspective
error, flow properties at grey area near plate are unavailable,
especially when α ≠ 90 º

63
Figure V.13 (c) the evolution of fluid flow during the first 2 strokes (black
arrow indicates velocity vector and colorbar depicts magnitude of
vorticity ω). Due to perspective error, flow properties at grey area
near plate are unavailable, especially when α ≠ 90 º

64
V.4.1 Case 1: The Effect of Plate Rotation
Lift approximation from surrogate modelling (Figure V.4) shows the tendency of
time-averaged lift coefficient to be higher if the rotation is advanced and lower
during delayed rotation, denoted by global maximum valley at tr < 0 and
minimum valley at tr > 0, respectively. Thus, this result confirms Dickinson’s
conclusion about the importance of rotation timing: that wing rotation must
supinate or pronate before a half stroke ends (Dickinson, Lehmann, & Sane,
1999). While the timing of rotation is critical in lift generation, the effect of
rotational duration is less pronounced than its timing, as shown in Figure V.5.

Figure V.14 Illustration of lift generation due to plate pronation at the first
downstroke, followed by upstroke in Case 1

Force history on optimum kinematics exhibits 2 peaks during a half-stroke, one at

the beginning and one at the end of each half-stroke (Figure V.6(b) and Figure
V.7(b)). The peak at the end suggested a lift enhancement kinematics known as
advanced rotation, where the plate supinates or pronates rapidly before a half
stroke ends. Rapid plate rotation in Case 1 is comparable to another phenomenon
called Magnus effect (see II.5) in spinning baseball trajectory. Both optimum
solutions in Figure V.14 exhibit “backspin”, where the plate pronates backward
relative to the direction of translation before a downstroke ends. Higher velocity at
the leading edge generates upward force component (positive lift, L > 0). On the
contrary, suboptimum solution in Figure V.14 exhibits “topspin” (the plate
pronates forward relative to the direction of translation), and thus generates

65
downward force (negative lift, L < 0). As the result of this kinematics, negative
lift at the end of every half stroke is apparent in force history of suboptimum
solution (see Figure V.8(b)).

A small force peak at the end of each half stroke due to rapid plate rotation is
immediately followed by negative lift, as shown in Figure V.6(b) and Figure
V.7(b). According to researches by Dickinson (1999) and Sane (2001), the force
acting on insect wing is roughly normal to the wing surface, even during rotation.
Hence, as the angle of attack exceeds 90º during rotation, the force vector changes
sign from positive to negative, resulting in lift drop at the very end of half stroke.

Wake capture and plate’s inertial force due to sudden acceleration are the causes
of the peak at the beginning of a half-stroke. Wake capture is a mechanism that
enhances lift solely from translational kinematics; hence it does not depend on
rotational timing. As the plate changes from downstroke to upstroke translation
and vice versa, the plate hits shed vortices from previous stroke, as shown by
Figure V.6(c) and Figure V.7(c) at t/T = 0.55, 1.05, and 1.55. Essentially, wake is
energy loss, but when the plate changes direction of translation, energy loss can be
recovered as the plate hits shed vortices from previous stroke. Dickinson (1999)
called this “a remarkable feature” of insect flight that insect is able to extract
energy from its own wake. Force peak due to wake capture, however, is not
apparent at the beginning of very first downstroke. The reason is quite obvious: at
the beginning of first downstroke, there is no shed wake from previous stroke.
Thus, small peak at that time is caused by plate’s inertial force during impulsive
start, similar to Lisoski’s case (see V.2).

Further observation on force history in both optimum and suboptimum solution

shows that the peak from wake capture on suboptimum solutions is not as high as
the peak on optimum solutions. Dickinson (1999) suggested that although wake
capture mechanism is generated by translational kinematics, its magnitude and
direction depend on interaction between translational-rotational kinematics. If the
rotation is advanced, there exists a large concentration of velocity magnitude at

66
the back of the plate, indicating large vorticity magnitude from shed vortices on
previous stroke (Figure V.15). As the plate travels in reverse direction, it captures
the wake at at a certain angle of attack that results in highly positive lift (Figure
V.6(b) and Figure V.7(b)). This large concentration of velocity magnitude,

Figure V.15 Velocity vector on the end of a half stroke in (a) Case 1 and (b)
Dickinson’s research (1999); colorbar depicts magnitude of
velocity. Both results are generated from PIV experiments in
towing tank.

Figure V.16 Result of parametric study of the effect of rotational kinematics on

lift generation by Sane (2001)

however, is not apparent in suboptimum solution, where the rotation is delayed

(Figure V.15). As the plate travels in reverse direction, the plate hits the flow at an
angle that produces negative lift. Hence, this mechanism explains why the lift is
highly negative between strokes in the case of delayed rotation (see Figure
V.8(b)).

Similar conclusion about wing rotation timing is also drawn by Sane (2001) in a
parametric study about the effect of flapping kinematics on force production (see
Figure V.16). In this research, Sane mapped the lift, drag, and L/D production of a

67
3-D insect wing model with 25-cm wingspan and Re = O(102), as a function of
time when rotation starts and the duration of rotation, using towing tank
experiments. Sane used 99 points across the map, where the starting time of the
rotation is limited to -0.5 ≤ tr ≤ 0 or half the search space investigated in this
thesis. The result is roughly similar to the advanced rotation part of the contour
(see Figure V.4(b)), where red color dominates that particular area, suggesting
highly positive lift, and blue area at the tip of the search space, suggesting
negative lift if the rotation is very rapid and starts at the very beginning of each
half stroke. In Sane’s work, it is noticeable that lift generally increases when the
rotation is faster and starts closer to the end of half stroke. That is not the case
with the result obtained from G-SAGA optimization (see Figure V.4(b)). Random
error is apparent in the histogram in Figure V.5, where 25 G-SAGA run result in
different optimum kinematics with large range of optimum solution, especially on
the duration of rotation. Random error can be further detected by comparison
between 2 optimum solutions chosen from histogram, where second optimum
solution, which has lower bar in the histogram, generates higher time-averaged lift
than the first solution. The presence of random error, plus the small range of lift in
advanced rotation area as shown by contour plot in Figure V.4(b), contributes to
the error in optimization and surrogate modelling result.

V.4.2 Case 2: The Effect of Stroke Amplitude and Angle of Attack

Lift approximation from surrogate modelling (Figure V.9) shows the peak of
time-averaged lift coefficient ranges between 20º ≤ α ≤ 65º and narrows down to
35º ≤ α ≤ 60º as the stroke amplitude increases. It can also be observed from the
contour plot that the gradient of time-averaged lift coefficient in the direction of
stroke amplitude is much smaller than angle of attack. This means that stroke
amplitude has little effect on lift generation and thus confirms the result of
previous researches by Sane (2001), Trizila (2008), and Trizila (2008).

As shown in Figure V.17(a), the lift tends to increase as the stroke amplitude
increases. This result is contradictory to result obtained from G-SAGA, where the
algorithm is always converged to smallest stroke amplitude (3 times the chord

68
length) as the optimum solution (see Figure V.10). There are two main reasons
that cause this result. First, the stroke period decreases with stroke amplitude.
Therefore, same non-dimensional rotational kinematics will result in faster
rotation as the stroke amplitude decreases, resulting in higher force peak. After
peak from wake capture, lift decreases into a plateau that ends during advanced
rotation before a half-stroke ends. When stroke period is lower, the length of the
plateau is shorter. Thus, when time-averaged lift coefficient is taken, kinematics
with lower stroke amplitude will generate higher lift. Second, although surrogate
modelling result by Trizila (2008) shows that stroke amplitude has insignificant
effect to 2-D airfoil lift, CFD analysis on several stroke amplitudes shows that
wake capture peak becomes smaller as stroke amplitude increases. If stroke
amplitude increases, the intensity of shed vortices decays before they interacted
with wing during the start of the following stroke (Trizila P. , Kang, Aono, &
Shyy, 2011).

Figure V.17 (a) Result of parametric study of the effect of stroke amplitude and
angle of attack on lift generation by Sane (2001); (b) result of
surrogate modelling of the effect of stroke amplitude and angular
amplitude (inversely proportional to angle of attack) on lift
generation of 2-D airfoil and 3-D wing by Trizila (2008)

In contrast to stroke amplitude, the effect of angle of attack on lift generation is

quite obvious. We have covered explanation about wake capture and rapid pitch
rotation which generate force peak at the beginning and the end of a half stroke,
respectively. What happens in a plateau between those two peaks will be
explained shortly. At 30º ≤ α ≤ 50º, plate generates force from the delayed stall

69
mechanism of the LEV (Sane & Dickinson, 2001). It was suggested several
previous researches that the presence of LEV is especially important to sustain lift
(and thus delay stall) in insect flight during translation, where the angle of attack
is high (see II.5). However, LEV in 2-D insect wing is only effective to delay stall
when it is attached to the leading edge. When LEV is detached from the wing and
shed, stall ocurrs, indicated by drop in the magnitude of lift (Sane, 2003).

The histogram in Figure V.10(b) shows that the optimum solution ranges between
35º ≤ α ≤ 55º. The dynamics of LEV will be examined through two sets of
optimum solution taken from Figure V.10(b). As shown by Figure V.11(c) and
Figure V.12(c), small LEV and TEV appear at the start of the first downstroke (t/T
= 0.05) and they continue to grow in size until t/T = 0.25 (the plate travels mid-
stroke, da = 1.5), where LEV is strong and about to detach and TEV is shed to the
wake. Following wake capture peak in the beginning of the first upstroke (t/T =
0.5), the LEV grows stronger until mid-stroke (t/T = 0.75), where the LEV is shed
to the wake (not pictured) and TEV grows stronger. As LEV is detached from the
plate, stall occurs, and thus explains the significant lift drop at mid-stroke (see
Figure V.11(c) and Figure V.12(c)). At the beginning of second downstroke, LEV
and TEV grows until TEV is shed at t/T = 1.1. At mid-stroke (t/T = 1.25), LEV is
shed to the wake, resulting in lift drop, and small TEV begins to form. Just before
rotation at t/T = 1.4, TEV is shed and a small LEV forms at the leading edge of
the plate. The process of alternating LEV-TEV formation and shedding generates
a set of counter-rotating vortices to the wake, a phenomenon known as “von
Karman vortex street” (Sane, 2003). Unfortunately, PIV experiments done for this
thesis are unable to capture the presence of von Karman vortex street due to the
limitation of camera’s ROI.

While the optimum solution exhibits delayed stall mechanism on a certain range
of high angle of attack, it is interesting to study what happens at low and very
high angle of attack (near 90º), where, as shown by G-SAGA result and Sane
(2001), the time-averaged lift coefficient is low. At low angle of attack, Sane
suggested that LEV and TEV generated during translation is small, thus the

70
contribution of lift is mostly given by rotation, since the plate must travel through
longer arc if the angle of attack is low.

At near-90º angle of attack, the resultant of force acting on plate is nearly parallel
to the direction of translation, thus the drag is very high and the lift is nearly zero,
as exhibited by force history of optimum solution (α = 80º) in Figure V.13(b).
Very high wake capture peak can be observed on drag history at every beginning
of half stroke, along with small peak at the end of half stroke due to advanced
plate rotation. This high wake capture peak is the result of plate hitting particulary
strong vortices at the end of half stroke (see Figure V.13(c) at t/T = 0.5). This
result comes as no surprise since the strength of vortices grow with increasing
angle of attack, as also found by Sane (2001).

71
 Chapter VI Conclusion and Future Works

VI.1 Conclusion
From the research done to accomplish the objectives in I.2, it can be concluded
that G-SAGA with Kriging developed for surface prediction and optimization has
successfully predicted lift of 2-D flapping plate as a function of flapping
kinematics and optimized the kinematics to obtain maximum lift. Coupling
between Kriging for surrogate modelling and GA for surrogate optimization
proves to be beneficial because it reduces number of fitness evaluation to a
hundred instead of thousands, as in the case of basic GA (Milano & Gharib,
2005). Furthermore, time required to evaluate fitness is reduced by evaluating
time-averaged lift coefficient in towing tank experiments instead of using
numerical simulations.

The result of first test case, which examines the effect of rotation on time-
averaged lift, shows that advanced rotation at the end of each half-stroke will
result in higher lift, thus moving towards optimum solution. Flow measurement
and visualization using PIV show that advanced rotation results in upward force
component from higher flow velocity at leading edge caused by backward plate
rotation, relative to the direction of translation. The peak is further increased at the
following stroke as the plate hits shed vortices from previous stroke, a mechanism
known as wake capture. While the timing of rotation is critical in lift generation,
the effect of rotational duration is less pronounced than its timing.

The result of second test case, which examines the effect of stroke amplitude and
translational angle of attack, shows that stroke amplitude has little effect on time-
averaged lift generation. However, maximum lift is achieved in a certain range of
angle of attack between 35º ≤ α ≤ 55º due to delayed stall. While LEV is still
attached to the plate at high angle of attack, lift is sustained. This mechanism is
not apparent if the angle of attack is low or nearly perpendicular to the direction of
translation.

72
VI.2 Future Works
From the research done to accomplish the objectives in I.2, it is hoped that in the
future:
(1) The research can be improved futher to optimize and study the effect of
several kinematics parameters at once on several aerodynamic design
objectives, such as lift, drag, L/D, or propulsive efficiency.
(2) The research can be improved improved futher to study the effect of
kinematics parameters on force generation of a 3-D insect wing model
with more complex kinematics that resembles real insect.
(3) The force measurement system can be further developed to separate
aerodynamic force from inertial force, so that the force generation from
fluid dynamics can be investigated thoroughly.

73
 REFERENCES

Adhynugraha, M. I. (2009). Optimization of Airfoil of Wing in Ground Effect

Using Genetic Algorithm, Master Thesis. Institut Teknologi Bandung.

Astrom, K. J., & Murray, R. M. (2008). Feedback Systems: An Introduction for

Scientists and Engineers. Princeton: Princeton University Press.

Autonics Corporation. (2015, April 28). Rotary Encoder (Incremental Type)

E50S/E50SP Series Manual. Retrieved from Autonics Corporation
Website:
http://download.autonics.com/upload/data/1420685595_E50S_E50SP_EN
_EP-KE-09-0140A_20141111.pdf

Birch, J. M., & Dickinson, M. H. (2003). The Influence of Wing-Wake

Interactions on The Production of Aerodynamics Forces in Flapping
Flight. Journal of Experimental Biology, 206, 2257-2272.

Choi, J. S., Kim, J. W., Lee, D. H., & Park, G. J. (2011). Optimization of the
Flapping Motion for the Hovering Flight. IEEE International Conference
on Robotics and Biomimetics, (pp. 1297-1301). Phuket.

Deb, K., & Kumar, A. (1995). Real-Coded Genetic Algorithms with Simulated
Binary Crossover: Studies on Multimodal and Multiobjective Problems.
Complaex Systems, 9, 431-454.

Dickinson, M. H., Lehmann, F., & Sane, S. P. (1999). Wing Rotation and the
Aerodynamic Basis of Insect Flight. Science, 284, pp. 1954-1960.

Dwianto, Y. B. (2015). Surrogate Assisted Genetic Algorithm for Expensive Fluid

Flow Optimization, Master Thesis. Institute of Technology Bandung.

Dwianto, Y. B., Palar, P. S., & Zuhal, L. R. (2014). Evolutionary Algorithm

Based Multi-Objective Aerodynamic Optimization Method for Low
Reynolds Number Airfoil. Applied Mechanics and Material, 660, 487-491.

Ellington, C. P. (1984). The Aerodynamics of Hovering Insect Flight I. The

Quasi-Steady Analysis. Philosophical Transactions of the Royal Society of
London. Series B, Biological Sciences, 305, 1-15.

Ellington, C. P., Van den Berg, C., Willmot, A. P., & Thomas, A. L. (1996).
Leading-edge Vortices in Insect Flight. Nature, 384, 626-630.

Forrester, A. I., Sobester, A., & Keane, A. J. (2008). Engineering Design via
Surrogate Modelling: A Practical Guide. Wiltshire: John Wiley & Sons
Ltd.

Haupt, R. L., & Haupt, S. E. (2004). Practical Genetic Algorithms. Hoboken,

New Jersey: John Wiley and Sons.

74
Ito, K. (2002). Optimization of Flapping Wing Motion. ICAS Congress, (pp.
814.1-814.7).

Jones, M., & Yamaleev, N. K. (2013). Adjoint-Based Shape and Kinematics

Optimization of Flapping Wing Propulsive Efficiency. 43rd AIAA Fluid
Dynamics Conference. San Diego: American Institute of Aeronautics and
Astronautics.

Le, M. N., Ong, Y. S., Menzel, S., Jin, Y., & Sendhoff, B. (2013). Evolution by
Adapting Surrogates. Evolutionary Computation, 21, 313-340.

Lisoski, D. L. (1993). Nominally 2-Dimensional Flow About a Normal Flat Plate,

Doctorate Dissertation. California Institute of Technology.

Liu, H., & Aono, H. (2009). Size Effects on Insect Hovering Aerodynamics: An
Integrated Computational Study. Bioinspiration and Biomimetics, 4, 1-13.

Liu, H., Ellington, C. P., Kawachi, K., Van den Berg, C., & Willmott, A. P.
(1998). A Computational Fluid Dynamic Study of Hawkmoth Hovering.
Journal of Experimental Biology, 201, 461-477.

Milano, M., & Gharib, M. (2005). Uncovering The Physics of Flapping Flat Plates
with Artificial Evolution. Journal of Fluid Mechanics, 534, 402-409.

Nguyen, H. L. (2009). Vortex Interaction and Force Generation in Flapping

Flight, Master Thesis. Institut Teknologi Bandung.

Nguyen, Q. N. (2006). Drag Due to Vortex Formation on Normal Flat Plate,

Master Thesis. Institut Teknologi Bandung.

Ogata, K. (1997). Modern Control Engineering (3rd ed.). Upper Saddle River,
New Jersey: Prentice-Hall Inc.

Palar, P. S. (2012). Experimental Optimization of Asymmetric Hovering Flapping

Flat Plate using Hybrid Artificial Evolution, Master Thesis. Bandung:
Institute of Technology Bandung.

Palar, P. S., & Zuhal, L. R. (2014). Flow Field Around Asymmetric Flapping Flat
Plate Optimized Using Micro Genetic Algorithm. 52nd AIAA Aerospace
Sciences Meeting - AIAA Science and Technology Forum and Exposition
(SciTech). National Harbor, MD: American Institute of Aeronautics and
Astronautics, Inc.

Percin, M., Hu, Y., van Oudheusden, B. W., Remes, B., & Scarano, F. (2011).
Wing Flexibility Effects in Clap-and-Fling. International Micro Air
Vehicles Conference.

Raffel, M., Willert, C., Wereley, S., & Kompenhans, J. (2007). Particle Image
Velociemtry: A Practical Guide (2nd ed.). Berlin: Springer.

75
Salles, R., & Schiele, A. (2004). Analysis of Flapping Wing Robots for Planetary
Exploration: An Evolutionary Approach. The 8th ESA Workshop on
Advanced Space Technologies for Robotics and Automation. Noordwijk.

Sane, S. P. (2003). The Aerodynamics of Flapping Flight. Journal of

Experimental Biology, 206, 4191-4208.

Sane, S. P., & Dickinson, M. H. (2001). The Control of Flight Force by a

Flapping Wing: Lift and Drag Production. Journal of Experimental
Biology, 204, 2607-2626.

Shyy, W., Aono, H., Chimakurthi, S. K., Trizila, P., Kang, C. K., Cesnik, C. E., &
Liu, H. (2010). Recent Progress in Flapping Wing Aerodynamics and
Aeroelasticity. Progress in Aerospace Sciences, 46, 284-327.

Shyy, W., Aono, H., Kang, C. K., & Liu, H. (2013). An Introduction to Flapping
Wing Aerodynamics. New York: Cambridge University Press.

Shyy, W., Lian, Y., Tang, J., Vieru, D., & Liu, H. (2008). Aerodynamics of Low
Reynolds Number Flyers. Cambridge: Cambridge University Press.

Stepen. (2011). Development, Evaluation, and Application of Image Processing

Code for Digital Particle Image Velocimetry (DPIV), Master Thesis.
Insitut Teknologi Bandung.

Trizila, P., Kang, C. K., Visbal, M. R., & Shyy, W. (2008). Unsteady Fluid
Physics and Surrogate Modeling of Low Reynolds Number, Flapping
Airfoils. 3th Fluid Dynamics Conference and Exhibit. Seattle,
Washington: AIAA 2008-3821.

Trizila, P., Kang, C. K., Visbal, M., & Shyy, W. (2008). A Surrogate Model
Approach in 2D versus 3D Flapping Wing Aerodynamic Analysis. 12th
AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference.
Victoria, British Columbia: AIAA-2008-5914.

Trizila, P., Kang, C.-K., Aono, H., & Shyy, W. (2011). Low Reynolds Number
Aerodynamics of a Flapping Rigid Flat Plate. AIAA Journal, 49, 806-823.

Wang, Z. J. (2004). The Role of Drag in Insect Hovering. Journal of Experimental

Biology, 207, 4147-4155.

Weis-Fogh, T. (1973). Quick Estimates of Flight Fitness in Hovering Animals,

Including Novel Mechanisms for Lift Production. Journal of Experimental
Biology , 59, 169-230.

Wuwung, V. S. (2008). Two Dimensional Particle Image Velocimetry

Measurement on The Downstroke Flapping Wing, Master Thesis. Institut
Teknologi Bandung.

76