Journal of Fluids and Structures: C.M. Hoke, J. Young, J.C.S. Lai
Journal of Fluids and Structures: C.M. Hoke, J. Young, J.C.S. Lai
Journal of Fluids and Structures: C.M. Hoke, J. Young, J.C.S. Lai
a r t i c l e in f o abstract
Article history: Research into the effectiveness of flapping foil propulsion and power extraction systems
Received 25 June 2014 has typically focused on the kinematic parameters defining the pitching and plunging
Accepted 3 May 2015 motion, and used simple basic rigid airfoil shapes. In this paper, the effects of a time-
Available online 15 June 2015
varying deformable foil shape on the power extraction and propulsive efficiency of
Keywords: flapping foil systems are studied. Two-dimensional Navier–Stokes solutions using the
Flapping propulsion commercial flow solver Fluent are performed with a deformable mesh to alter the camber
Foil shape of the foil during the flapping cycle. All simulations are assumed to be fully laminar, with a
Flapping power extraction free stream Reynolds number of Re ¼ 1100 for the power extraction cases and Re ¼ 20 000
Low Reynolds number
for the propulsion cases, chosen to match representative cases in the existing literature.
Camber deformation
The shape of the foil is constructed by deforming the camber line via a circular arc
Navier–Stokes simulation
centered at the mid-span, with the magnitude of the circular arc deformation varying
sinusoidally in time. The phase angle of this sinusoidal camber variation is the primary
independent variable, resulting in a broad range of foil interactions with the shed vortex
flow fields as the shape varies over the flapping cycle. This camber deformation is
superimposed on kinematically constrained sinusoidal motions for both the power
extraction and propulsion regimes, as well as semi-passive motions for power extraction.
The results show that the efficiency of these systems can be increased by judiciously
deforming the foil shape to interact with the resulting leading and trailing edge vortex
structures. For all cases, the power required to deform the foil surface is included in the
overall efficiency calculation, and for several cases is a significant factor in the final
efficiency. The performance of the power extraction cases is primarily determined by the
interactions between the shed LEV and the foil horizontal surface during the plunging
stroke and the interaction of the trailing edge as it passes through the vortex during the
pitch reversal portion of the cycle. The semi-passive cases for power extraction are also
strongly affected by the vortex interactions, and the resulting forces and moments can
significantly alter the flapping frequency, further increasing their impact. The propulsion
cases show a high dependence on the interaction of the leading edge vortex with the foil
curvature to provide high leading suction and greater propulsive forces. For all three
regimes, regions of phase angle where the efficiency is increased by the deformation of
the foil are shown.
& 2015 Elsevier Ltd. All rights reserved.
n
Corresponding author. Tel.: þ61 26268 8209.
E-mail address: [email protected] (C.M. Hoke).
http://dx.doi.org/10.1016/j.jfluidstructs.2015.05.001
0889-9746/& 2015 Elsevier Ltd. All rights reserved.
C.M. Hoke et al. / Journal of Fluids and Structures 56 (2015) 152–176 153
1. Introduction
Flapping wings, such as those used in nature by birds, insects, and aquatic animals, show great promise in both
propulsive applications and as novel energy extraction devices. However, due to the complexity of the flow physics and
fluid–structure interactions involved in such flapping flight, there is still room to explore how to make such systems more
efficient. Flapping wing systems have the potential to create much higher lift and thrust forces than more typical fixed-wing
systems for the appropriate flow regimes (Shyy and Liu, 2007; Bandyopadhyay, 2002). Devices that take advantage of these
forces show great promise to revolutionize both the micro-UAV and energy extraction fields.
These flapping airfoil systems are difficult to study experimentally, and therefore much work has been done to simulate the
physics of flapping wing flight with Computational Fluid Dynamics (CFD). CFD is useful due to the ease of computing a variety of
different flow regimes and Reynolds numbers, altering the airfoil shape, and implementing a variety of airfoil motions that would
be difficult to recreate in an experimental setting. However, limitations exist to the computational approach, particularly relating
to accurate computation of transitional and turbulent flows at moderate Reynolds numbers (Peng and Zhu, 2009) and the
extreme computational cost of 3-dimensional solutions (Visbal, 2009; Ashraf et al., 2012; Kinsey and Dumas, 2012). Due to these
limitations, much of the computational work has been performed at low Reynolds numbers, thus precluding the possibility of
transitional flow and reducing the effect of that 3D flow structures would have in a higher Re 3D simulation.
154 C.M. Hoke et al. / Journal of Fluids and Structures 56 (2015) 152–176
The power extraction problem presents a unique challenge because biological examples of flapping wings are exclusively
of the propulsive variety. Therefore, a great deal of work has been done to determine the optimum sinusoidal kinematic
parameters for power extraction, such as flapping frequency, the phase shift between pitch and plunge, the amplitude of
plunging motion, and the maximum pitch angle of the foil during the flapping cycle (McKinney and DeLaurier, 1981; Young
and Lai, 2004; Tuncer and Kaya, 2005; Kinsey and Dumas, 2005, 2008; Ashraf et al., 2012; Young et al., 2014; Xiao and Zhu,
2014). In addition to these sinusoidal cases, there has been work into non-sinusoidal flapping motions, including pitch-angle
control and angle-of-attack control schemes(Hover et al., 2004; Ashraf et al., 2011a; Platzer et al., 2009; Xiao et al., 2012).
Many of these different schemes show promise to produce high efficiencies, and alternative parameterizations warrant
more study as more efficient motions are found. Semi-passive and fully kinematically passive cases where either the plunge
or pitch motions (or both) are allowed to develop naturally in response to the oncoming flow have also shown the potential
for high efficiency, as they do not require any power input into the system to maintain a prescribed motion profile (Young
et al., 2010; Zhu, 2012; Young et al., 2013; Kinsey and Dumas, 2014).
The flapping foil propulsion problem is a relatively mature field due to the copious examples in nature and the
applicability to micro-UAVs. Several papers cover the history of flapping wing propulsion quite well (Lai and Platzer, 1999;
Shyy et al., 1999; Rozhdestvensky and Ryzhov, 2003; Platzer et al., 2008). There has been extensive computational work, and
indeed many of these have also focused on the sinusoidal kinematic parameters of flapping foil propulsion (Streitlien and
Triantafyllou, 1998; Anderson et al., 1998; Young and Lai, 2004, 2007a; Tuncer and Kaya, 2005; Shyy and Liu, 2007; Bos et al.,
2008; Kinsey and Dumas, 2014), as well as non-sinusoidal motions (Hover et al., 2004; Zhu and Peng, 2009; Xiao and Liao,
2010; Ashraf et al., 2012; Lu et al., 2013).
There are very few studies on shape effects in the power extraction problem, with the vast majority of cases using
symmetric foils. There has been considerable study into the effects of thickness on both symmetric NACA airfoils (Kinsey and
Dumas, 2008) and flat plates (Usoh et al., 2012) or even scallop shell shapes (Le et al., 2013), lending weight to the argument
that the ideal shape for power extraction is yet to be determined. All of these studies focus on static shapes, with very little
work on shapes that actively change shape during the flapping cycle. A major exception to this is the work by Liu et al.
(2013), which used prescribed deformations of both the leading and trailing edges of the foil. However, deflections of the
leading and trailing edges produce artificially altered pitching angles, making it difficult to isolate the effects of the shape
alone. In addition, the power required to deform the foil is not calculated.
For the propulsion problem, shape effects have been investigated more thoroughly. Recently, several studies using CFD
have looked at the effect of foil thickness (Wang, 2000; Rozhdestvensky and Ryzhov, 2003; Cebeci et al., 2005; Ashraf et al.,
2012) and a few have focused on the effect on static camber deformations (Lentink and Gerritsma, 2003; Ashraf et al.,
2011b). There are a few studies on the effect of an actively varying camber line on a purely plunging foil (Miao and Ho, 2006;
Tay and Lim, 2010). However, the camber deformation in both cases was achieved by deflecting the trailing edge, which
artificially pitches the foil, so comparing the results to other purely plunging cases, or even cases with the same nominal
pitch angle, is problematic because the effect of the shape is difficult to separate from the effect of the resultant pitch angle.
Passive deformation, or flexibility, has been more extensively studied, especially for the propulsion regime. Flapping
wings in nature, especially insect wings, are often passively flexible, and indeed this flexibility has been found to be quite
important in the efficiency of these systems (Miao and Ho, 2006; Heathcote and Gursul, 2007; Hoogedoorn et al., 2010;
Unger et al., 2012; Tian et al., 2013). All of these studies determined that passive flexibility was very important for the
generation of thrust, and indeed was necessary for the kind of performance characteristics seen in biological examples.
Similar studies for flapping power extraction, however, are sparse, as most work has focused on small-scale fluttering
piezioelectric devices such as the “energy harvesting eel” (Allen and Smits, 2001; Taylor et al., 2001) and vortex induced
vibration of cylinders or other shapes (Bernitsas et al., 2008; Oh et al., 2010).
In this study, a deforming foil undergoes a prescribed deformation during the flapping cycle for three representative
motions, as shown in Table 1, and the resulting efficiencies are compared with published rigid foil results. For all of the flow
regimes considered, the Reynolds number at which the simulations are performed makes a large difference in the resulting
flow physics. In this paper, the power extraction results are chosen to run at Re ¼ 1000, while the propulsion cases are run at
Re ¼ 20 000. Both of these values are chosen primarily to match with previous studies as denoted in the table to give a
readily available basis for comparison. These previous studies all make a fully laminar boundary layer assumption that is
justifiable for these low Reynolds numbers, albeit marginally for the propulsion cases. However, as many researchers have
noted (Ashraf et al., 2011b; Kinsey and Dumas, 2014; Young and Lai, 2007b; Young et al., 2013), these Reynolds numbers do
not accurately match any macro scale devices that could be built easily. This is a known limitation of these results.
Table 1
Representative rigid foil motion cases used for baseline comparisons.
1. Kinematically constrained sinusoidal pitch and plunge Power extraction Kinsey and Dumas (2008)
2. Semi-passive flow driven α control Power extraction Young et al. (2013)
3. Kinematically constrained sinusoidal pure plunging Propulsion Miao and Ho (2006)
C.M. Hoke et al. / Journal of Fluids and Structures 56 (2015) 152–176 155
In all cases, the camber deformation is superimposed on a best-case condition, which is a combination of flapping
parameters identified in the corresponding reference that produced the highest efficiency for the flow regime of interest. All
three of the reference cases (Kinsey and Dumas, 2008; Young et al., 2013; Miao and Ho, 2006) performed parameter searches
to identify these best case conditions, so any increase in efficiency can be attributed to the deformation of the foil alone.
2. Methods
2.1. Kinematics
Fig. 1. Pitch/plunge relationship at ϕ ¼ 90, reproduced from Young et al. (2014) with permission. (a) Propulsion mode and (b) power Extraction mode.
Fig. 2. Schematic of the semi-passive flapping wing turbine, reproduced from Young et al. (2013) with permission.
156 C.M. Hoke et al. / Journal of Fluids and Structures 56 (2015) 152–176
top and the bottom of the pitching sinusoid and has a force vector that acts in the direction of travel during the stroke, as shown in
Fig. 1. These cases are all calculated at a Reynolds number of 1100. The sinusoidal pitching and plunging motions are defined by
V y ¼ H 0 θ_ cos ðθt
_ þ ϕÞ; ð1Þ
_
θ ¼ θ0 sin ðθtÞ: ð2Þ
Table 2
Motion breakpoints for semi-passive angle of attack controlled motion (0 r β r 2πÞ.
Table 3
Equations for angle of attack linkage function.
β gðβÞ
β o β1 αmax
β β1
β o β2 αmax þ αmax ð1 cos ðπzÞÞ, z ¼
β2 β1
β o β3 αmax
β β3
β o β4 αmax αmax ð1 cos ðπzÞÞ, z ¼
β4 β3
β 4β4 αmax
C.M. Hoke et al. / Journal of Fluids and Structures 56 (2015) 152–176 157
Fig. 3. Angle of attack linkage function f ðβÞ and resulting pitch angle θ, Δβ^ R ¼ 0:2.
In the general case, the equation of motion for this system is modeled as a flywheel with a viscous damper attached, as
well as linear and rotational spring–damper systems attached to the foil itself. Ignoring the mass of the moving linkage
elements, the equation of motion is determined via conservation of energy as
mfoil y_ y€ þ I foil θ_ θ€ þ I fly β_ β€ ¼ Ly_ þM θ_ cplunge y_ 2 cfoil θ_ cfly β_ kplunge yy_ kfoil θθ:
_
2 2
ð5Þ
Consistent with the work by Young et al. (2013) in which this semi-passive configuration is introduced, the system is
simplified by assuming that the flywheel mass is an order of magnitude greater than that of the airfoil, and that the only
non-zero spring or damper element is the power-extracting flywheel itself. An actual turbine system, of course, will have
non-zero values for these elements. However, as we are comparing these results with other idealized 2D cases that use
proscribed motions that do not specify the mechanism by which the foil is moved, these idealizations are necessary to make
meaningful comparisons of their efficiencies. The generalized equation as well as the simplified equation used in this work
is given in the following equations:
1 h 2i
β€ ¼ 2
Lf β þMg β cfly β_ mfoil f β f ββ þ I foil g β g ββ β_ ; ð6Þ
mfoil f β þ I foil g β þI fly
2
1h i
β€ ¼ Lf β þMg β cfly β_ : ð7Þ
I fly
The other kinematic parameters of the system used for all of the semi-passive runs here are a non-dimensional flywheel
inertia I 0 ¼ 5:67 (Eq. (8)), a chord length and plunging amplitude of 0:15c, an assumed span length of 1c, a pivot location at
1=2c, and a nondimensional flywheel damper strength c0fly ¼ 3:0 (Eq. (9)). These parameters are all reproduced from the
optimum flow driven reference case. These non-dimensionalizations are derived in Young et al. (2013), which is derived by
analogy with the approach used by Zhu and Peng (2009) and Zhu et al. (2009):
I fly
I 0fly ¼ ; ð8Þ
1 2
πρh H o c2 s
4
cfly
cfly 0 ¼ : ð9Þ
1
πρc3 U 1 s
16
All simulations are done using the commercial CFD code Fluent (version 14.5) with an unsteady Navier–Stokes flow
solver and a fully laminar boundary layer assumption. The convective terms are discretized with a second-order upwind
scheme and the diffusion terms are discretized with a second-order central differencing scheme. Time discretization is done
via a second-order implicit scheme. The plunging motion is modeled with a time-dependent momentum source term
introduced into the governing equations applied over the entire computational domain, which rotates the freestream flow
158 C.M. Hoke et al. / Journal of Fluids and Structures 56 (2015) 152–176
Fig. 4. Foil structured mesh details and sliding mesh motion setup.
Table 4
Details of grid and time-step refinement studies for sinusoidal motion cases.
direction to the proper angle corresponding to the instantaneous vertical plunge velocity, hence the simulation is performed
in the heaving reference frame.
The boundary conditions utilize a velocity inlet condition on the left, top, and bottom farfield boundaries that specify a
time-varying vertical velocity to be consistent with the momentum source term. A two-zone mesh with a circular structured
inner portion and a rectangular unstructured triangular outer region is used to accommodate the pitching motion of the foil,
with a sliding mesh interface between the two sections (Fig. 4). The circular mesh zone has a radius of four chord lengths,
which puts it far enough away from the airfoil surface to minimize interactions with the flow structures but keeps the
relative grid motions between the two zones low. The circular pitching motions and the deformable mesh used to modify
the foil shape are implemented by iteratively specifying the grid point locations of the circular moving mesh zone explicitly
for every time step, described in detail in Section 2.3. The mesh for the flow driven motion uses a slightly different grid that
is further refined at the circular interface zone to better resolve the grid motion for increased rotational velocities, but is
otherwise similar.
Table 5
Details of grid and time-step refinement studies for semi-passive cases.
Fig. 7. Comparison of semi-passive grid and time step studies to Kinsey and Dumas baseline lift and power coefficients.
The time-varying camber line variation is implemented with a deforming mesh and a simple circular arc based camber
line definition about the mid-span of the foil. This approach maintains the location of the leading and trailing edges of the
foil and therefore maintains a constant chord length and chord line so that the pitch angle and angle of attack are not
functions of the deformation. This ensures that the effects of the deformation can be isolated from the kinematic motion of
the foil. The maximum chord line deflection is H c ¼ 0:1c and the camberline is recomputed at every time step via Eq. (10) for
the sinusoidal motions, or Eq. (11) for the semi-passive cases. The grid deformation is implemented by storing the original
grid point locations of the circular rotating mesh zone at case initialization, and explicitly recalculating and updating their
160 C.M. Hoke et al. / Journal of Fluids and Structures 56 (2015) 152–176
Fig. 8. Definition of circular arc camber line and resulting deformed grid.
Fig. 9. Foil deformation and motion example, sinusoidal motion (ϕc ¼ 01).
coordinates at each time step. The camberline deformation and an example of the deformed grid are shown in Fig. 8:
hcam; sinusoidal
_ þ ϕc Þ
¼ H c sin ðθt ð10Þ
The mesh surrounding the foil is updated by rigidly displacing all grid points within 0:1c in the y direction to maintain
the proper boundary layer spacing and normal direction growth rates next to the foil surface. At a distance from 0:1c to 1:1c,
the grid deformation is a simple linear function of distance from the airfoil. For clarity the deformation superimposed on the
sinusoidal motion case is displayed in Fig. 9.
The phase angle of the camber deformation, ϕc, is the primary independent variable, and both the instantaneous
deformation at a given point in the flapping cycle as well as the time history of the deformation leading up to that point has
a strong influence on the resulting flow field and power extraction efficiency. Table 6 details the relationship of the camber
phase angle to the motion of the foil.
C.M. Hoke et al. / Journal of Fluids and Structures 56 (2015) 152–176 161
Table 6
Relationship of ϕc with the foil pitch/plunge motion ( þ θ¼ pitch down).
In the case of a proscribed deformation, additional power is required to deform the foil. This quantity may be positive or
negative depending on the pressure distribution around the foil, as the localized pressure gradient may assist or resist the
deformation at a given point in the cycle. The power due to deformation at each time step is calculated by integrating the
dot product of the local face velocity vector with the pressure and shear forces at each face of the foil, Eq. (14). The face
velocity vector is calculated in the foil reference frame, and this quantity is added to the pitching and plunging power:
I I
! !
P def ¼ ^ face P face dAface þ
V face n V face !
τ face Þ dAface : ð14Þ
The lift force, L(t), and moment about the pivot point, M(t) both contribute to the generation of power as the foil moves,
and the total aerodynamic power extracted is the time average of these two quantities integrated over the flapping cycle,
plus the deformation power, shown in non-dimensional form in Eq. (15). The power extraction efficiency, η, is defined as the
ratio of the power extracted over one flapping cycle to the power available in the flow as defined by the Betz criterion,
Eqs. (16) and (17):
Z 1 _
V y ðtÞ θðtÞc t
C Ptotal ¼ C Paero þC Pdef ¼ C P y þC P θ þC P def ¼ C y ðt Þ þC M ðt Þ þ C Pdef d ; ð15Þ
0 U1 U1 T
P P y þP θ c
ηaero ¼ ¼ ¼ C Paero ; ð16Þ
P a 1ρU 31 d d
2
P P y þ P θ þ P def c
ηtotal ¼ ¼ ¼ C Ptotal : ð17Þ
Pa 1
2ρU 31 d d
where d is the total swept distance of the foil including the trailing edge:
For the semi-passive cases, the power extracted is related to the time-average rate of energy dissipation in the flywheel
damper. This provides a physical mechanism for the power extraction, similar to a rotary system. The power required to
deform the foil is presumed to be extracted or added from the flywheel power, depending on whether it is a positive or
negative contribution, and is included in the efficiency calculation:
Z 1 t
cfly β_ þP def d
2
P¼ : ð19Þ
0 T
The propulsive efficiency is defined as the ratio of the propulsive power to the power consumption:
CT
ηpropulsion;total ¼ : ð21Þ
C Ptotal
3. Results
Foil camber phase angles from 180 r ϕc r 1801 in 22.51 intervals are simulated for 12 flapping cycles each to ensure
that the flow is fully developed and all startup transients have dissipated. The effect of camber on the flapping foil efficiency
is surprisingly strong, resulting in lower aerodynamic efficiencies over a large range of camber phase angles, but a significant
region of increased aerodynamic efficiency with a maximum of 37.9% at ϕc ¼ 1351 as shown in Fig. 10. The deformation
power was also a strong contributor to the overall efficiency, as shown by the CPtotal line on the phase efficiency plot. For
these kinematically constrained cases, the deformation power contribution is negative for phase angles less than zero
degrees, as in general the grid motion is opposed by the pressure gradient in those regions. This shifted the overall optimum
case from ϕc ¼ 1351 to ϕc ¼ 1351. Clearly, in any scheme that includes an active deformation of the airfoil this power
contribution can drastically change the ideal operating parameters. However, for the flow visualizations in the remainder of
this paper we will focus on the best aerodynamic efficiency to illustrate the interaction of the foil shape with the flow fields.
Clearly, the shape of the foil can have a strong effect on the overall aerodynamic efficiency of the flapping cycle. There are
several phenomena that determine the overall cycle efficiency, including most importantly the interaction between the foil
shape and the leading and trailing edge vortex structures, the shape of the airfoil during the plunging stroke, and
the effective angle of attack during the stroke. The power extraction history over a flapping cycle reveals the areas of
favorable interactions, shown in Fig. 11, by areas where the instantaneous power coefficient is higher than the non-
deforming baseline case. Interestingly, the best case does not produce the most aerodynamic power at any given portion of
the cycle, but the time integration of power production over the entire cycle is maximum. This demonstrates the difficulty in
optimizing the foil shape – maximizing the power output over any given portion of the flapping cycle will not necessarily
Fig. 11. Instantaneous coefficient of power for a cycle for all kinematically constrained cases.
C.M. Hoke et al. / Journal of Fluids and Structures 56 (2015) 152–176 163
Fig. 12. Power extraction parameters for four representative sinusoidal cases.
Fig. 13. Comparison of best and worst cases for sinusoidal power extraction, t=T ¼ 0:25.
Fig. 14. Comparison of best and worst cases for sinusoidal power extraction, t=T ¼ 0:45.
Fig. 15. Comparison of best and worst cases for sinusoidal power extraction, t=T ¼ 0:65.
of the trailing edge with this vortex has a large effect on the pitching motion power extraction, and a favorable interaction is
what determines the overall contribution of the C Pθ term to the overall efficiency. If the foil shape interacts with the flow to
form the LEV at the right time in the cycle so that it has translated to the trailing edge at the time of the pitch reversal
maneuver, the low pressure helps us to “flip” the foil, resulting in a high efficiency as the pressure field resultant force
vectors act in the direction of the prescribed motion.
During the power stroke at t=T ¼ 0:25, the best case at ϕc ¼ 1351 exhibits the weakest LEV. All of the other cases show a
more coherent LEV vortex at this instant, but the low pressure associated with the vortex has minimal horizontal surface to
act on, which results in low lift and power extraction. The best case, conversely, is curved in such a way that there is more
horizontal surface to act on, and indeed this case has the highest value of CPy and CP at this instant. This is a surprising result,
C.M. Hoke et al. / Journal of Fluids and Structures 56 (2015) 152–176 165
as the large coherent vortex typically associated with flapping foil flow dynamics is delayed in this case. This portion of the
flapping cycle is very important to the overall efficiency, as the sinusoidal cases get the majority of their power during the
plunging stroke. This is clearly exhibited in the vorticity magnitude and pressure coefficient contours shown in Fig. 13,
which compares the best case to the worst case at the midpoint of the plunging cycle.
At t=T ¼ 0:45, just prior to the lower apex of the flapping cycle (Fig. 14), the LEV has translated to the trailing edge, and
the foil is rotating through it. This vortex, in the best case, has translated aft the least of all the cases, due primarily to the
delay in formation discussed at t=T ¼ 0:25. Therefore, later in the stroke at t=T ¼ 0:50, the vortex is still strong and coherent
when the foil is rotating through it. By comparison, in the low efficiency case, the vortex has moved off the trailing edge and
is not inducing a favorable pressure gradient to assist with the rotational motion. The ϕc ¼ 01 case has the best C Pθ at this
instant, as the trailing edge is actually in the process of moving through the vortex core and beginning to break it up.
Even though the best case exhibits the highest overall efficiency, there are still regions where the pressure fields are
unfavorable. Just after t=T ¼ 0:50, there is a negative spike in C Pθ for the best case as the break up of the vortex when the TE
passes through it temporarily forms a secondary vortex on the opposite side of the foil that resists the pitching maneuver. Of
the four representative cases, the best case actually has the lowest C Pθ at this instant. The increased efficiency during the
plunging stroke makes up for this brief period of low power generation, however.
Finally, at t=T ¼ 0:65 (Fig. 15), the beginning of the upward plunging stroke, it is again observed that the baseline and best
cases have the most horizontal surface for the lifting forces to act upon, producing the highest CPy and delaying the LEV
formation.
This work elucidates several important flow features that are important to the generation of power. For a sinusoidally flapping
foil where the pitch rates are relatively low, the plunging power generation term, CPy, is the primary contributor to the overall
power extraction. Therefore, deforming the foil to provide the most lift at the beginning and middle of the power cycle is
important, but the pitch angles must still be sufficient to generate a strong leading edge vortex so that the trailing edge can
interact with it during the pitch reversal portion of the cycle. It is possible, as exhibited by the low efficiency cases, to generate the
strong LEV too soon in the cycle, disrupting the complex flow interaction needed for high efficiency at pitch reversal.
Clearly, the vortex interactions can be influenced by the shape of the foil to increase the efficiency of the system, but
determining the optimum shape histories is a difficult prospect. Using this simple sinusoidal variation of the camber line,
the aerodynamic efficiency in the best case has been increased from η ¼ 0:329 to η ¼ 0:381, a relative 15.8% increase.
However, this is only a refinement of a single case, and the optimum shape history will change drastically with any change
in Reynolds number, foil shape, or motion kinematics.
It should be emphasized again that the inclusion of the deformation power has a large effect on the overall efficiency of
the system. The preceding paragraphs focus on the aerodynamic efficiency in isolation to investigate the interaction of the
foil surface with the flow structures, as that was the original focus of the research. However, for any real system that uses a
proscribed motion, the contribution of the power required to deform the foil will clearly determine the overall ideal
operating conditions.
All the semi-passive cases are initialized with a non-deformable startup case designed to get the flywheel spinning and
eliminate initial transients. Once this startup case has reached a steady state oscillation (approximately 12 flapping cycles),
the deformation is ramped in and the oscillation allowed to stabilize. The number of cycles needed to reach a steady
oscillation varies with camber phase angle, but to ensure that any restart transients have dissipated all cases are run for at
least 12 more flapping cycles. The phase angle is varied from 1801 r ϕc r 1801, at intervals of 22.51.
Similar to the kinematically constrained sinusoidal cases, the results show a region where the efficiency is increased, and
in fact the variations in power extraction with camber deformation are far more pronounced than for the sinusoidal cases.
This is primarily a result of the flow-driven nature of the motion, because any unfavorable interaction with the pressure field
not only has an instantaneous reduction in the aerodynamic power output, it also results in a reduction in the rate of
rotation of the flywheel. In fact, the highest efficiency cases correspond exactly to the cases that produce the shortest
Fig. 16. Efficiency and flapping cycle duration dependence on camber phase angle – semi-passive motion.
166 C.M. Hoke et al. / Journal of Fluids and Structures 56 (2015) 152–176
Fig. 17. Instantaneous power coefficients over a cycle for all semi-passive cases.
Fig. 18. Instantaneous power coefficients over a cycle for four semi-passive cases.
C.M. Hoke et al. / Journal of Fluids and Structures 56 (2015) 152–176 167
average cycle duration, which is a clear consequence of the semi-passive power equation, Eq. (19). The overall efficiency and
flapping cycle durations are shown in Fig. 16.
The deformation power again was a large contributor to the overall efficiency, although in these semi-passive cases the contribution
does not change the optimal phase angle. Note that the power required to deform the airfoil did not affect the flapping frequency or
flywheel rate in this simulation, but is considered to be a net gain or loss in power that is supplied after the flywheel power is considered.
This would require some sort of driving mechanism that would deform the foil to be powered by the flywheel output power, but the
details of the machinery are left abstracted here as the primary focus is on the foil shape interactions with the flow fields.
The range of favorable interactions centers around ϕc ¼ 157:51 for these cases, with a maximum relative efficiency
increase of 18.7%. Fig. 17 shows the power extraction coefficient of all the semi-passive cases over one flapping cycle. As for
the kinematically constrained cases, the best case does not have the highest power coefficient for most of the flapping cycle,
but the time integration over the entire cycle yields the highest overall efficiency.
In Fig. 18, four representative cases are again compared with the power contributions broken up into plunging, pitching,
aerodynamic, deformation, and total contributions. Four representative cases again considered, consisting of the baseline case,
the best case at ϕc ¼ 157:51, the worst case at ϕc ¼ 901, and a medium efficiency case at ϕc ¼ 451. Again, the terms “best,
medium, and worst” refer specifically to the aerodynamic efficiencies. The main feature of these plots when compared to the
kinematically constrained motions of Fig. 12 is the dominance of the C Pθ term. This arises from the much higher pitch rate
during the stroke reversal portions of the flapping cycle that are a consequence of the angle of attack controlled flapping
motion. This in combination with the fact that the highest efficiency case has the highest flapping frequency means that the
most important factor in the semi-passive case is the interaction of the foil with the shed vortex structure at the time of pitch
reversal. Several of the non-optimum cases actually exhibit stronger lift characteristics during the power stroke, denoted by
higher CPy values, and indeed the best case, ϕc ¼ 157:51, is among the lowest CPy at certain instants in the cycle. However, the
large increase in efficiency during the stroke reversal, clearly seen in the plot of C Pθ , results in a larger overall efficiency for a
complete flapping cycle. The details of the flow structure and the foil shape interactions for these four cases are visualized in
the Appendix, Figs. A3 and A4. Figs. 19–21 are extracted from the appendix data here to elucidate some key points.
During the middle of the plunging stroke at β^ ¼ 0:50 (Fig. 19), both the baseline case and the best case show a LEV that
has just shed and remains relatively close to the foil surface. Both foils are nearly vertically orientated as the pitch increases
Fig. 19. Comparison of best and worst cases for semi-passive power extraction, β^ ¼ 0:50.
Fig. 20. Comparison of best and worst cases for semi-passive power extraction, β^ ¼ 0:70.
168 C.M. Hoke et al. / Journal of Fluids and Structures 56 (2015) 152–176
to maintain αmax ¼ 401, which accounts for their low CPy during the stroke. In contrast, the worst case experiences its
maximum positive camber at this point in the cycle, and it is clear that the shape of the foil preceding this point in the cycle
causes the LEV to form much earlier, and it has shed well off the surface of the foil by this time. This is primarily because
more time has elapsed during the stroke due to the longer flapping cycle of the low-efficiency case. While this has a small
effect on the CPy at this point, the position of the vortex is of paramount importance as the flapping cycle continues.
The pitching maneuver begins at β^ ¼ 0:65, and in the flow visualization immediately following its initiation at β^ ¼ 0:70 (Fig. 20)
again some similarity between the baseline and best cases is observed, as the LEV has propagated downstream and is positioned to
interact with the TE. The low efficiency case, however, has taken much longer to reach this point in the cycle and the initial strong LEV
has induced a new strong secondary vortex underneath the foil and a new leading edge vortex has appeared. These two vortices are
rotating in opposite directions and interact during the pitching maneuver. The medium efficiency case, ϕc ¼ 451, exhibits a mix of the
high and low efficiency behaviors, with the initial vortex closer to the foil and the secondary vortex structures not as fully developed.
At the apex of the pitching motion, β^ ¼ 0:75 (Fig. 21), the angular rates are the highest. It is evident here that the trailing edge
of the highest efficiency case has rotated into the shed LEV, which greatly assists the pitching maneuver and is the reason for the
quite high C Pθ and resulting high efficiencies. The baseline case trailing edge also approaches the vortex, although has not
entered the vortex core as quickly as the best case. However, in the low efficiency case, the secondary vortex and the new leading
edge vortex have interacted and largely broken each other up, resulting in a much weaker pressure gradient in the direction of
the pitching maneuver. The medium case, again, exhibits a mix of these two flow patterns, and gets some assistance from the
initial LEV but still has a complex secondary vortex interaction that does not significantly impact the overall efficiency.
Finally, during the beginning of the plunging stroke, shown at β^ ¼ 0:90 in Fig. 22, the baseline case and best cases have
not yet formed a strong leading edge vortex, and in fact the best case at ϕc ¼ 157:51 has a positive camber that, like the
kinematically constrained cases, is curved in the proper direction to provide more horizontal surface area to provide lift.
However, the angle of attack control has resulted in a higher overall pitch angle, so the effect on CPy is less pronounced. The
low efficiency case, however, has already shed a strong LEV at this point in the cycle. As seen earlier, this LEV propagates too
far from the foil to be of much use to generate power.
Fig. 21. Comparison of best and worst cases for flow-driven power extraction, β^ ¼ 0:75.
Fig. 22. Comparison of best and worst cases for semi-passive power extraction, β^ ¼ 0:90.
C.M. Hoke et al. / Journal of Fluids and Structures 56 (2015) 152–176 169
The pure plunging propulsion cases also exhibited an increase in efficiency for a narrow range of camber deformation
phase angles, with a maximum propulsive efficiency of 22.2% as compared to the baseline case of 16.2%, a relative increase of
37.3%. The effects of the deformation power are quite small in these propulsion cases. The increases in efficiency center
around a camber phase angle of ϕc ¼ 01, which corresponds to maximum positive camber at the middle of the downward
stroke when the foil is plunging downward the fastest, maximum negative camber during the upward stroke, and a
symmetric profile at the upper and lower points of the cycle. The worst case is at ϕc ¼ 1351, which reverses the camber
deflection, i.e. has maximum positive camber during the upstroke. Fig. 23 shows the dependence of the propulsive efficiency
on the phase angle of deformation.
Similar to the power extraction cases, the best case solution does not produce the maximum propulsive force over the
entire flapping cycle, as clearly displayed in Fig. 24. The effect of the camber is quite pronounced, with the high efficiency
cases producing over twice the propulsive force for certain portions of the cycle.
Comparing the baseline, the best case at ϕc ¼ 01, the worst case at ϕc ¼ 1351, and a medium case at ϕc ¼ 901 reveal the
trends in the thrust and power coefficient variations over the flapping cycle. The best case solution produces high thrusts
during large portions of the flapping cycle, particularly during the second half of each plunging stroke, from t=T ¼ 0:25 to
t=T ¼ 0:50 or so. The lowest efficiency case, at ϕc ¼ 1351 by contrast, is producing net drag over most of the flapping cycle.
The medium case is interesting, however, as it starts to produce high thrust coefficients earlier in the plunging stroke than
the best case. These four cases are compared in Fig. 25.
To fully understand the physical flow structures that are causing these behaviors, the contours of vorticity magnitude,
instantaneous streamlines, and pressure coefficients for these representative cases have been included in fully in the
appendix (Figs. A5 and A6) with Figs. 26–28 extracted from the appendix data to elucidate some key points. The best case
produces its high thrust during the plunging cycle after the LEV vortex has formed, and the deformation of the camber line
has curved the foil in such a way as to provide a large area for the low pressure in the vortex to act upon, inducing a large
amount of leading edge suction in the horizontal thrusting direction. This effect is clearly seen in Fig. 26, which occurs just
after the midpoint of the plunging stroke.
Interestingly, it is clear that the leading edge vortex is much more developed at this point in the cycle for the baseline
case, and the pressure coefficients are lower, but the interaction of the vortex with the foil is not favorable. Similar to the
power extraction cases, it appears that there is benefit to delaying the formation of the vortex if the shape interactions
170 C.M. Hoke et al. / Journal of Fluids and Structures 56 (2015) 152–176
Fig. 25. Instantaneous coefficient of power over a cycle for four representative propulsion cases.
Fig. 26. Comparison of baseline and best cases for propulsion during the plunging stroke, t=T ¼ 0:30.
warrant it, even if the overall vortex strength is diminished. However, for the “medium” efficiency case at ϕc ¼ 901, there is
high efficiency earlier in the plunging stroke than there is in the best case, shown in Fig. 27.
This suggests that there may be even more favorable shape interactions that could be conceived through more complex
methods of foil deformation that produces an early, coherent vortex and maintains a shape that continually provides a leading
edge suction. The worst case, predictably, alters the shape in the least favorable direction and in some cases results in increased
drag. However, the worst case also exhibits the strongest vortex structures, which is a surprising result, shown in Fig. 28. It is quite
clear that the vortex is much more strongly developed for the worst case due to the high relative curvature of the leading edge to
the oncoming free stream flow, but it is also clear that the shape of the foil does not take advantage of this vortex structure.
C.M. Hoke et al. / Journal of Fluids and Structures 56 (2015) 152–176 171
Fig. 27. Comparison of best and medium cases for propulsion early in the plunging stroke, t=T ¼ 0:20.
Fig. 28. Comparison of best and worst cases for propulsion, showing disparity in vortex strength, t=T ¼ 0:20.
The shape of the foil used in a flapping foil power generation system has a large effect on the overall efficiency of the
system, and careful manipulation of the foil shape during the flapping cycle may yield significant benefits. We found that for
well-studied kinematically constrained sinusoidal flapping motions, judicious variation of the circular camber line can
increase the relative aerodynamic efficiency of the system by 15.8%. The timing of the formation of the leading edge vortex
during the high pitch angle plunging strokes is a key factor in the efficiency of these cases, both due to the interaction of this
vortex during the plunging stroke to create lift and the position of the vortex during the pitching maneuver to induce
favorable pressure gradients in the direction of rotation during the pitch reversals. The power required to deform the foil in
these cases was a large contributor to the overall efficiency, where it changed the phase angle for maximum efficiency
considerably.
Semi-passive cases that connect a foil to a flywheel exhibit a similar benefit over the non-flexible optimum, increasing
relative efficiency by 18.7%. The complex flow physics and interactions between the shed vortex structures with the foil
surface are the primary contributors to the final efficiencies. The overall dependence on the camber is more pronounced due
to the effect of the pressure fields on the overall flapping cycle frequency. Additionally, the higher pitch angles and faster
pitch rates in the semi-passive cases mean that the rotational power term is more important, and therefore the vortex
interaction with the trailing edge is more important than they are for the kinematically constrained cases.
The propulsion cases showed a similar region of increased efficiency, with the best case increasing the relative efficiency
by 37.3%. The primary shape interaction in these cases is the interaction of the shed vortices with the leading edge curvature.
In cases where the airfoil is curved such that the leading edge has a lot of vertical surface to act on as the shed LEV passes
over it, the propulsive efficiency is quite high. Similar to the power extraction cases, the highest efficiency cases often delay
the formation of the LEV and reduce its strength, but the positive interactions of the shape with the vortex that does form
are enough to offset the weaker vortex.
These insights warrant much further investigation, as it is quite possible that more complex functions for both the pitch
and plunge motions as well as camber deformation could increase the advantageous flow interactions further. Flat plates
and scalloped foil configurations have shown promise in other works, and deformation of these basic structures could also
172 C.M. Hoke et al. / Journal of Fluids and Structures 56 (2015) 152–176
Fig. A1. Instantaneous streamlines and vorticity magnitude during the flapping cycle for four kinematically constrained sinusoidal motion cases.
Fig. A2. Contours of pressure coefficient during the flapping cycle for four kinematically constrained sinusoidal motion cases.
C.M. Hoke et al. / Journal of Fluids and Structures 56 (2015) 152–176 173
Fig. A3. Instantaneous streamlines and vorticity magnitude during the flapping cycle for four semi-passive motion cases.
Fig. A4. Contours of pressure coefficient during the flapping cycle for four semi-passive motion cases.
174 C.M. Hoke et al. / Journal of Fluids and Structures 56 (2015) 152–176
Fig. A5. Instantaneous streamlines and vorticity magnitude during the flapping cycle for four plunging propulsion cases.
Fig. A6. Contours of pressure coefficient during the flapping cycle for four plunging propulsion cases.
C.M. Hoke et al. / Journal of Fluids and Structures 56 (2015) 152–176 175
produce high efficiency interactions. Of course, fully coupled fluid structure interaction solvers with shapes that deform in
response to the pressure field are a logical next step, as they more closely mirror insect wings and avoid any sort of
mechanical structure needed to enforce a proscribed deformation.
References
Allen, J.J., Smits, A.J., 2001. Energy harvesting eel. Journal of Fluids and Structures 15 (3-4), 629–640.
Anderson, J.M., Streitlien, K., Barrett, D.S., Triantafyllou, M.S., 1998. Oscillating foils of high propulsive efficiency. Journal of Fluid Mechanics 360 (1), 41–72.
Ashraf, M.A., Young, J., Lai, J.C.S., Platzer, M.F., 2011a. Numerical analysis of an oscillating-wing wind and hydropower generator. AIAA Journal 49 (7),
1374–1386.
Ashraf, M.A., Young, J., Lai, J.C.S., 2011b. Reynolds number, thickness and camber effects on flapping airfoil propulsion. Journal of Fluids and Structures 27
(2), 145–160.
Ashraf, M.A., Young, J., Lai, J.C.S., 2012. Oscillation frequency and amplitude effects on plunging airfoil propulsion and flow periodicity. AIAA Journal 50 (11),
2308–2324.
Bandyopadhyay, P.R., 2002. Maneuvering hydrodynamics of fish and small underwater vehicles. Integrative and Comparative Biology 42 (1), 102–117.
Bernitsas, M.M., Raghavan, K., Ben-Simon, Y., Garcia, E.M., 2008. Vivace (vortex induced vibration aquatic clean energy): a new concept in generation of
clean and renewable energy from fluid flow. Journal of Offshore Mechanics and Arctic Engineering 130 (4), 041101.
Bos, F.M., Lentink, D., Van Oudheusden, B.W., Bijl, H., 2008. Influence of wing kinematics on aerodynamic performance in hovering insect flight. Journal of
Fluid Mechanics 594 (1), 341–368.
Cebeci, T., Platzer, M., Chen, H., Chang, K.C., Shao, J.P., 2005. Analysis of Low-Speed Unsteady Airfoil Flows. Horizons Publishing Inc, Springer, Long Beach,
California, USA.
Heathcote, S., Gursul, I., 2007. Flexible flapping airfoil propulsion at low reynolds numbers. AIAA Journal 45 (5), 1066–1079.
Hoogedoorn, E., Jacobs, G.B., Beyene, A., 2010. Aero-elastic behavior of a flexible blade for wind turbine application: a 2D computational study. Energy 35
(2), 778–785.
Hover, F.S., Haugsdal, O., Triantafyllou, M.S., 2004. Effect of angle of attack profiles in flapping foil propulsion. Journal of Fluids and Structures 19 (1), 37–47.
Kinsey, T., Dumas, G., 2005. Aerodynamics of oscillating wings and performance as wind turbines. In: The 23rd AIAA Applied Aerodynamics Conference,
Toronto, Ontario, Canada.
Kinsey, T., Dumas, G., 2008. Parametric study of an oscillating airfoil in a power-extraction regime. AIAA Journal 46 (6), 1318–1330.
Kinsey, T., Dumas, G., 2012. Computational fluid dynamics analysis of a hydrokinetic turbine based on oscillating hydrofoils. Journal of Fluids Engineering
134 (2). 021104 (16 pp).
Kinsey, T., Dumas, G., 2014. Optimal operating parameters for an oscillating foil turbine at reynolds number 500,000. AIAA Journal 52 (9), 1885–1895,
http://dx.doi.org/10.2514/1.J052700.
Lai, J.C.S., Platzer, M.F., 1999. Jet characteristics of a plunging airfoil. AIAA Journal 37 (12), 1529–1537.
Le, T., Ko, J., Byun, D., 2013. Morphological effect of a scallop shell on a flapping-type tidal stream generator. Bioinspiration and Biomimetics 8 (3), 11.
Lentink, D., Gerritsma, M., 2003. Influence of airfoil shape on performance in insect flight. In: The 33rd AIAA Fluid Dynamics Conference and Exhibit,
Orlando, Florida.
Liu, W., Xiao, Q., Cheng, F., 2013. A bio-inspired study on tidal energy extraction with flexible flapping wings. Bioinspiration and Biomimetics 8 (3), 036011.
Lu, K., Xie, Y.H., Zhang, D., 2013. Numerical study of large amplitude, nonsinusoidal motion and camber effects on pitching airfoil propulsion. Journal of
Fluids and Structures 36 (1), 184–194.
McKinney, W., DeLaurier, J., 1981. Wingmill: an oscillating-wing windmill. Journal of Energy 5 (2), 109–115.
Miao, J.M., Ho, M.H., 2006. Effect of flexure on aerodynamic propulsive efficiency of flapping flexible airfoil. Journal of Fluids and Structures 22 (3), 401–419.
Oh, S.J., Han, H.J., Han, S.B., Lee, J.Y., Chun, W.G., 2010. Development of a tree-shaped wind power system using piezoelectric materials. International Journal
of Energy Research 34 (5), 431–437.
Peng, Z., Zhu, Q., 2009. Energy harvesting through flow-induced oscillations of a foil. Physics of Fluids 21 (12), 1–9.
Platzer, M.F., Jones, K.D., Young, J., Lai, J.C.S., 2008. Flapping wing aerodynamics: progress and challenges. AIAA Journal 46 (9), 2136–2149.
Platzer, M., Young, J., Lai, J.C.S., Ashraf, M., 2009. Development of a new oscillating-wing wind and hydropower generator. In: The 47th AIAA Aerospace
Sciences Meeting, Orlando, FL.
Rozhdestvensky, K.V., Ryzhov, V.A., 2003. Aerohydrodynamics of flapping-wing propulsors. Progress in Aerospace Sciences 39 (8), 585–633.
Shyy, W., Liu, H., 2007. Flapping wings and aerodynamic lift: the role of leading-edge vortices. AIAA Journal 45 (12), 2817–2819.
Shyy, W., Berg, M., Ljungqvist, D., 1999. Flapping and flexible wings for biological and micro air vehicles. Progress in Aerospace Sciences 35 (5), 455–505.
Streitlien, K., Triantafyllou, G.S., 1998. On thrust estimates for flapping foils. Journal of Fluids and Structures 12 (1), 47–55.
Tay, W.B., Lim, K.B., 2010. Numerical analysis of active chordwise flexibility on the performance of non-symmetrical flapping airfoils. Journal of Fluids and
Structures 26 (1), 74–91.
Taylor, G.W., Burns, J.R., Kammann, S.A., Powers, W.B., Welsh, T.R., 2001. The energy harvesting eel: a small subsurface ocean/river power generator. IEEE
Journal of Oceanic Engineering 26 (4), 539–547.
Tian, F.B., Luo, H., Song, J., Lu, X., 2013. Force production and asymmetric deformation of a flexible flapping wing in forward flight. Journal of Fluids and
Structures 36, 149–161.
Tuncer, I.H., Kaya, M., 2005. Optimization of flapping airfoils for maximum thrust and propulsive efficiency. AIAA Journal 43 (11), 2329–2336.
Unger, R., Haupt, M.C., Horst, P., Radespiel, R., 2012. Fluid–structure analysis of a flexible flapping airfoil at low reynolds number flow. Journal of Fluids and
Structures 28, 72–88.
Usoh, C.O., Young, J., Lai, J.C.S., Ashraf, M.A., 2012. Numerical analysis of a non-profiled plate for flapping wing turbines. In: The 18th Australasian Fluid
Mechanics Conference, Launceston, Australia.
Visbal, M.R., 2009. High-fidelity simulation of transitional flows past a plunging airfoil. AIAA Journal 47 (11), 2685–2697.
Wang, Z.J., 2000. Vortex shedding and frequency selection in flapping flight. Journal of Fluid Mechanics 410 (1), 323–341.
Xiao, Q., Liao, W., 2010. Numerical investigation of angle of attack profile on propulsion performance of an oscillating foil. Computers and Fluids 39 (8),
1366–1380.
Xiao, Q., Zhu, Q., 2014. A review on flow energy harvesters based on flapping foils. Journal of Fluids and Structures 46, 174–191.
Xiao, Q., Liao, W., Yang, S., Peng, Y., 2012. How motion trajectory affects energy extraction performance of a biomimic energy generator with an oscillating
foil? Renewable Energy 37 (1), 61–75.
Young, J., Lai, J.C.S., 2004. Oscillation frequency and amplitude effects on the wake of a plunging airfoil. AIAA Journal 42 (10), 2042–2052.
Young, J., Lai, J.C.S., 2007a. Vortex lock-in phenomenon in the wake of a plunging airfoil. AIAA Journal 45 (2), 485–490.
Young, J., Lai, J.C.S., 2007b. Mechanisms influencing the efficiency of oscillating airfoil propulsion. AIAA Journal 45 (7), 1695–1702.
176 C.M. Hoke et al. / Journal of Fluids and Structures 56 (2015) 152–176
Young, J., Ashraf, M.A., Lai, J.C.S., Platzer, M.F., 2010. Numerical simulation of flow-driven flapping-wing turbines for wind and water power generation. In:
The 17th Australian Fluid Mechanics Conference, Auckland, New Zealand.
Young, J., Ashraf, M.A., Lai, J.C.S., Platzer, M.F., 2013. Numerical simulation of fully passive flapping foil power generation. AIAA Journal 51 (11), 2727–2739.
Young, J., Lai, J.C.S., Platzer, M.F., 2014. A review of progress and challenges in flapping foil power generation. Progress in Aerospace Sciences 67, 2–28.
Zhu, Q., 2012. Energy harvesting by a purely passive flapping foil from shear flows. Journal of Fluids and Structures 34 (1), 157–169.
Zhu, Q., Peng, Z., 2009. Mode coupling and flow energy harvesting by a flapping foil. Physics of Fluids 21 (3), 10.
Zhu, Q., Haase, M., Wu, C.H., 2009. Modeling the capacity of a novel flow-energy harvester. Applied Mathematical Modelling 33 (5), 2207–2217
doi: http://dx.doi.org/10.1016/j.apm.2008.05.027.