Grid dimensions used for the grid dependence study.
Simulating the three-dimensional flow features generated by heaving wings constitutes a great challenge due to the computational effort required to compute the complex three-dimensional flow produced as a function of the kinematics parameters, wing geometry, and Reynolds number. Hereafter, we study the wake topology generated by oscillating rigid wings and the validity of the Strouhal number as the fundamental parameters used to assess the aerodynamic performance of heaving wings. The unsteady laminar incompressible Navier-Stokes equations are solved on moving overlapping structured grids using a second-order accurate in space and time finite-difference numerical method. The numerical simulations are performed at a Reynolds number of Re = 250 and at different values of Strouhal number and heaving frequency.
- heaving wings
- aerodynamic performance
- wake topology
- Strouhal number
- overlapping grids
Flapping wings for flying and oscillating fins for swimming stand out as one of the most complex yet widespread propulsion methods found in nature. Natural flyers and swimmers (which have evolved over millions of years) represent illuminating examples of biokinetics, unsteady aerodynamics, high maneuverability, endurance, and large aero/hydrodynamics efficiency.
In the field of flapping flight, biologist, zoologist, and engineers are sharing findings and conducting research together. From the point of view of a biologist or zoologist, studying flapping flight in nature is of great importance for understanding the biology, allometry, flight patterns, flight skills, and the migratory habits of avian life. From an engineering point of view, the main reason for studying flapping flight is the use of animal locomotion as inspiration for improving existing applications or developing new technologies by just mimicking nature evolutionary-optimization process (biomimetics). Such applications may include drag reduction, noise reduction, and flow control by using feather-like structures  and flippers tubercles ; the development of new propulsion/lift generation systems for micro-air-vehicles (MAVs), nano-air-vehicles (NAVs), and autonomous-underwater-vehicles (AUVs) is inspired by flapping wings or oscillating fins [3, 4, 5, 6, 7, 8], energy harvesting applications , and even robotic extraterrestrial exploring missions .
An important aspect of flying using flapping wings or swimming by using oscillating fins or fanning is the ability to generate thrust with relatively high propulsive efficiency. Early attempts at building fish-inspired mechanisms achieved disappointingly low propulsive efficiencies . It was only through a deeper understanding of the vorticity and wake produced by swimming animals, significant progress was achieved .
Many researchers [13, 14, 15] have found that flying and swimming animals cruise in a narrow range of Strouhal numbers (between 0.2 and 0.4), corresponding to a regime of vortex growth and shedding in which the propulsion efficiency peaks. The Strouhal number
Hence, it becomes evident that gaining a better understanding of the wing/fin motion parameters driving forces generation, vortices generation and shedding, the manner in which the vortices interact with the moving surfaces and themselves, and how they contribute to lift and propulsion would aid in better understanding the propulsion mechanism of birds, insects, and fishes, independently of their possible practical applications.
In the current numerical study, we aim at performing a comprehensive analysis of the wake signature and aerodynamic performance of finite-span rigid wings undergoing pure heaving motion. The laminar incompressible Navier-Stokes equations are numerically approximated, and all unsteady, viscous, and three-dimensional effects are solved. The simulations are conducted for Strouhal numbers values between 0.15 ≤
The remainder of this paper is organized as follows. In Section 2, we give a brief description of the numerical method and gridding methodology. In Section 3, we present a description of the computational domain, case setup, and heaving kinematics. In Section 4, we present a short discussion of the quantitative and qualitative results obtained from a grid dependence study. In Section 5, we present a detailed discussion of the results. Finally, in Section 6, we present the conclusions and future developments.
2. Numerical method
Hereafter, we summarize the numerical method used to solve the governing equations on structured overlapping grids. For a complete description of the numerical method and gridding methodology, the interested reader should refer to the papers by Henshaw , Henshaw and Petersson , and Chesshire and Henshaw .
In primitive variables (
with the following boundary conditions and initial conditions,
In this IBVP, the vector
An alternative formulation of the system of Eqs. (3)–(6), called the velocity-pressure formulation, can be written as follows,
with the following boundary and initial conditions,
The system of Eqs. (7)–(11) is equivalent to the original formulation (Eqs. (3)–(6) [16, 19, 20]) and is the form of the equations that will be discretized. Hence, we look for an approximate numerical solution of Eqs. (7) and (8), in a given domain , with prescribed boundary conditions and given initial conditions (Eqs. (9)–(11)). Eqs. (7)–(11) are solved in logically rectangular grids in the transformed computational space (the interested reader should refer to the following papers for a detailed derivation [16, 18, 21, 22]), using second-order centered finite-difference approximations on structured overlapping grids.
The structured overlapping grids method consists in generating a set of conforming structured component grids that completely cover the domain that is being modeled in physical space and overlap where they meet. In this newly generated overlapping grid system , domain connectivity is obtained through proper interpolation in the overlapping areas, and in the case of moving bodies, grid connectivity information is recomputed at each time step.
As the problem of heaving wings is implicitly a problem with moving bodies, we need to solve the governing equations in a frame that moves with the component grids. For moving overlapping grids, Eqs. (7) and (8) can be expressed in a reference frame moving with the component grids as follows,
where represents the velocity of the components grids. It is worth noting that the new governing equations (Eqs. (12) and (13)) must be accompanied by proper boundary conditions. For a moving body with a corresponding moving no-slip wall (e.g., a heaving wing), we should impose the velocity at the wall as follows,
After spatial discretization of the governing equations, we can proceed with the temporal discretization. By proceeding in this way, we are using the method of lines (MOL) [23, 24]. The main advantage of the MOL method is that it allows to select numerical approximations of different accuracy for the spatial and temporal terms. Each term can be treated differently to yield different accuracies.
At this point, the discretized equations can be integrated in time by using any time discretization method. In this work, we used a semi-implicit multistep method, where the viscous terms are approximated using the Crank-Nicolson scheme, and the convective terms are approximated using the Adams-Bashforth/Adams-Moulton predictor-corrector approach (where the velocity is advanced in time by using a second-order Adams-Bahforth predictor step, followed by a second-order Adams-Moulton corrector step). We chose to implicitly treat the viscous terms because if they were treated explicitly, we could have a severe time step restriction proportional to the spatial discretization squared. The solution method previously described yields to a second-order accurate in space and time numerical scheme on structured overlapping grids.
To assemble the overlapping grid system and solve the laminar incompressible Navier-Stokes equations in their velocity-pressure formulation, the overture framework was used (http://www.overtureframework.org/). The large sparse system of linear algebraic equations arising from the discretization of the laminar incompressible Navier-Stokes equations is solved using the PETSc library (http://www.mcs.anl.gov/petsc/), which was interfaced with overture. The system of equations is then solved using a Newton-Krylov iterative method, in combination with a suitable preconditioner and in parallel computational architectures.
3. Geometry, boundary conditions, initial conditions, and wing kinematics
, we present an illustration of the overlapping grid system layout used to conduct this parametric study. In the figure,
In all the cases studied, a rectangular wing with an aspect ratio
The initial conditions used for all the heaving wings simulations are those of a fully converged solution of the corresponding fixed wing case. In
(left), the left boundary of the BG corresponds to an inflow boundary condition (
In all the simulations conducted, we assumed that the wing is undergoing pure heaving motion, wherein the wing cross-section heaves in the vertical direction (or
4. Grid dependence study and vortical structures visualization
When conducting numerical simulations, it is well known that the grid resolution can affect the results from a quantitative and qualitative point of view. In this section, we present the outcome of the grid dependence study used to determine the best overlapping grid system in terms of computing time, stability, and accuracy of the solution. To conduct this study, we used the grid convergence index method or GCI, as described by Roache in [25, 26].
During this study, fixed and moving wings were considered, but for simplicity, we will only present the results related to pure heaving motion (where the uncertainties are higher). Several simulations were run at a Strouhal number
The different grid sizes, layouts, and mesh stretching ratios considered during this study are shown in
, the grid spacing ratio (GSR) is the refinement ratio from the finer grid to the coarser grid and is expressed in reference to the position of the first node normal to the wing surface (1NW). Therefore, the grid spacing refinement ratio
|161 × 121 × 101||221 × 121 × 41||81 × 61 × 41||1||0:001 × 2c|
|161 × 121 × 101||201 × 101 × 31||61 × 51 × 31||2||0:002 × 2c|
|161 × 121 × 101||161 × 81 × 31||51 × 41 × 31||4||0:004 × 2c|
Since in the study of heaving wings propulsion, the main task is thrust production, and it is more convenient to think in terms of thrust force
In Eq. (16), the lift force
, we present the average drag coefficient
and the root mean square of the lift coefficient
for the grid dependence study. These values were used to compute the observed order of convergence
Based on the outcome of the GCI study (refer to
at zero grid spacing is estimated to be 0.057708, with an error band of 0.397203% for grid
and an error band of 1.641617% for grid
at zero grid spacing is estimated to be 0.908738 with an error band of 0.402503% for grid
and an error band of 1.567203% for grid
. The constancy of
|(GCI23/GCI12) × (1/
To supplement the quantitative GCI study, we conducted an additional grid dependence study but from the qualitative point of view (vortical structures resolution on the wake of the wing). The qualitative results for grids
) are illustrated in
, we use the iso-surfaces of vorticity magnitude (∣
Summarizing the quantitative and qualitative results previously presented, we can conclude that the solutions obtained by using the overlapping grid systems
are grid independent. Taking into account the computational resources available, CPU time restrictions, and solution accuracy,
with a value of 1 NW equal to 0.001 × 2
5. Simulation results
Hereafter, we carry out a comprehensive parametric study to assess the wake signature and aerodynamic performance of heaving rigid wings. In
, we present the kinematics parameters governing the heaving motion (described by Eq. (15)). In
, we also present the quantitative results obtained, where
is the average thrust coefficient and
is the maximum lift coefficient (which was measured during the downstroke). By inspecting these results, we observe that as we increase
, we can read that for values of
To understand how the heaving wing generates thrust, we need to take a closer look at the vortex shedding process. In Figures 5 and 6 , we illustrate the vortex topology for a heaving case in the thrust producing regime (case 3DH-11 in Table 4 ). These figures show that the downstream wake (or far field wake) of this case consists of two sets of doughnut-shaped vortex rings (VR) which convect at oblique angles about the centerline of the heaving motion. Thus, the flow induced by each vortex ring along its axis is expected to have a net streamwise component linked to thrust production. This net streamwise momentum excess in the wake is connected with the thrust production of heaving wings. The process by which the vortex rings are formed can be explained by examining the vortex formation and shedding close to the wing’s surface. During the heaving motion and close to the wing, four vortices are formed; namely, one leading edge vortex (LEV), one trailing edge vortex (TEV), and two wing-tip vortices or WTV (one WTV on the left wing-tip or LH-WTV and another WTV on the right wing-tip or RH-WTV). These four vortices are all connected and form a closed vortex loop (VL). During the heaving motion and as this vortex loop is convected downstream, it disconnects from the wing, creating in this way the doughnut-shaped vortex rings. It is also of interest in the fact that each vortex loop has two sets of thin contrails (TC) that connect the VL to the VR generated in the previous stroke; these structures are segments of the wing-tip vortices, and, as the vortex rings are convected downstream, they become weaker and ultimately disappear (as in vortex ring VR2 in Figure 5 ). During a complete heaving cycle, two VRs are formed, one at the end of the upstroke and the other one at the end of the downstroke.
In Figure 7 , we present the wake topology for a drag-producing case (which corresponds to case 3DH-3 in Table 4 ). It is clear from this figure that the wake topology is very different from the one of the thrust producing case. In this case, as the vortex loops are convected downstream, they do not morph into vortex rings. Instead, they keep their original shape, and as they are convected, they diffuse. We can also observe that the wake height is very compact, in comparison to that of the thrust production case (as depicted in Figure 8 ). Finally, notice how the flow induced by each vortex loop is inclined in the same direction of the wing’s travel direction, resulting this in a momentum surfeit linked to drag production. The momentum deficit and momentum excess scenarios can be better appreciated in Figure 8 , where A corresponds to drag production (momentum deficit) and B corresponds to trust production (momentum excess).
Animations of several simulation cases are available at the author’s web site: http://www.dicat.unige.it/guerrero/flapsim1.html (last accessed: Sep 2017).
6. Conclusions and perspectives
In this manuscript, we studied the unsteady aerodynamics of heaving rigid wings. The laminar incompressible Navier-Stokes equations were solved in their velocity-pressure formulation using a second-order accurate in space and time finite-difference numerical method, and to efficiently deal with moving bodies, we used overlapping structured grids. To study the dependence of the aerodynamic forces and wake topology on the wing kinematics, many simulations were conducted at different values of Strouhal number and at two reduced frequency values (low and high oscillating frequency).
The simulations show that the wake of thrust producing, rigid heaving wings is formed by two sets of interconnected vortex loops that slowly convert into vortex rings as they are convected downstream. It is also observed that the vortex rings are inclined with respect to the free-stream flow, whereas for thrust producing configurations, the angle of inclination of the vortex rings is in the same direction of their travel, and for drag-producing configurations, the angle of inclination of the vortex rings is opposite to the direction of their travel. The presence of thin contrails that link the vortex loops is of interest; these structures are segments of the wing-tip vortices, and as the vortex loops are convected downstream, they become weaker and ultimately disappear. In general, the observed structures are qualitatively similar to those observed in the experiments by Parker et al.  and Von Ellenrieder et al.  and the numerical simulations of Dong et al.  and Blondeaux et al. .
From the force measurement study, two different behaviors were observed for the average thrust coefficient
and maximum lift coefficient
. It seems that the reduced frequency
Finally, for the limited range of
The results presented in this manuscript are limited to laminar flow; nevertheless, they provide an excellent insight into the wake signature of the unsteady aerodynamics of heaving wings. We envisage to extend the current study to higher Reynolds numbers and turbulent cases and use more realistic wing geometries and kinematics. Finally, in this manuscript, we did not cover propulsive efficiency and optimal frequency selection, but we hope to address these issues in future studies.
The use of the computing resources at CINECA high performance computing center was possible thanks to the ISCRA grant 2016, project IsC45_DO4EnD. We also thank the DLTM (http://www.dltm.it) for letting us use their HPC facilities.
Favier J, Dauptain A, Basso D, Bottaro A. Passive separation control using a self-adaptive hairy coating. Journal of Fluid Mechanics. 2009; 627:451-483
Favier J, Pinelli A, Piomelli U. Control of the separated flow around an airfoil using a wavy leading edge inspired by humpback whale flippers. Comptes Rendus Mecanique. 2012; 340:107-114
Jones KD, Platzer MF. Experimental investigation of the aerodynamic characteristics of flapping-wing micro air vehicles. In: AIAA Paper 2003-0418-CP. 2003. pp. 1-11
Singh SN, Simha A, Mittal T. AUV maneuvering by pectoral fins: Inverse control design based on CFD parametrization. IEEE Journal of Oceanic Engineering. 2004; 29:777-785
McMichael JM, Francis MS. Micro Air Vehicles toward a New Dimension in Flight. Technical report, Defense Advanced Research Projects Agency (DARPA). 1997
Guerrero J. Wake signature and aerodynamic performance of finite-span root flapping rigid wings. Journal of Bionic Engineering. 2010; 7:S109-S122
Guerrero J, Pacioselli C, Pralits JO, Negrello F, Silvestri P, Lucifredi A, Bottaro A. Preliminary design of a small-sized flapping UAV I. Aerodynamic performance and static longitudinal stability. Meccanica. 2016; 51:1343-1367
Negrello F, Silvestri P, Lucifredi A, Guerrero J, Bottaro A. Preliminary design of a small-sized flapping UAV II. Kinematic and structural aspects. Meccanica. 2016; 51:1369-1385
Orchini A, Mazzino A, Guerrero J, Festa R, Boragno C. Flapping states of an elastically anchored plate in a uniform flow with applications to energy harvesting by fluid/structure interaction. Physics of Fluids. 2013; 25:097105–1/097105–17
Michelson R, Naqvi M. Extraterrestrail flight (entompter-based mars surveyor). Von Karman Institute for Fluid Dynamics. Lecture Series. November 2003; 2003:24-28
Triantafyllou GS, Triantafyllou MS. An efficient swimming machine. Scientific American. 1995; 272:40-48
Rohr J, Fish F. Strouhal number and optimization of swimming by odontocete cetaceans. The Journal of Experimental Biology. 2004; 207:1633-1642
Triantafyllou MS, Triantafyllou GS, Gopalkrishnan R. Wake mechanics for thrust generation in oscillating foils. Physics of Fluids. 1991; 3:2835-2837
Nudds RL, Taylor GK, Thomas AR. Tuning of strouhal number for high propulsive efficiency accurately predicts how wingbeat frequency and stroke amplitude relate and scale with size and flight speed in birds. Proceedings of the Biological Sciences. 2004; 7:2071-2076
Taylor GK, Nudds RL, Thomas AR. Flying and swimming animals cruise at a strouhal number tuned for high power efficiency. Letters to Nature. 2003; 425:707-711
Henshaw WD. A fourth-order accurate method for the incompressible Navier-Stokes equations on overlapping grids. Journal of Computational Physics. 1994; 113:13-25
Henshaw WD, Petersson NA. A split-step scheme for the incompressible Navier-Stokes equations. Numerical Solutions of Incompressible Flows. Preprint UCRL-JC-144. 2003. pp. 1-14
Chesshire G, Henshaw W. Composite overlapping meshes for the solution of partial differential equations. Journal of Computational Physics. 1990; 90:1-64
Petersson N. Stability of pressure boundary conditions for stokes and navier-stokes equations. Journal of Computational Physics. 2001; 172:40-70
Sani RL, Shen J, Pironneau O, Gresho PM. Pressure boundary condition for the time-dependent incompressible navier-stokes equations. International Journal for Numerical Methods in Fluids. 2006; 50:673-682
Guerrero J. Numerical Simulation of the Unsteady Aerodynamics of Flapping Flight PhD thesis. Italy: University of Genoa, Deparment of Civil, Enviromental and Architectural Engineering; 2009. pp. 1-220
Vinokur M. Conservation equations of gas-dynamics in curvilinear coordinate systems. Journal of Computational Physics. 1974; 14:105-125
Schiesser WE. The Numerical Method of Lines: Integration of Partial Differential Equations. San Diego, California: Academic Press; 1991. pp. 1-326
Wesseling P. Principles of Computational Fluid Dynamics. Berlin, Heidelberg: Springer; 2001. pp. 1-644
Roache P. Verification and Validation in Computational Science and Engineering. Albuquerque, New Mexico: Hermosa Publishers; 1998. pp. 1-464
Roache P. Quantification of uncertainty in computational fluid dynamics. Annual Review of Fluid Mechanics. 1997; 29:123-160
Jeong J, Hussain F. On the identification of a vortex. Journal of Fluid Mechanics. 1995; 285:69-94
Haimes R, Kenwright D. On the velocity gradient tensor and fluid feature extraction. In: AIAA Paper 1999–3288-CP. 1999. pp. 1-10
Wang ZJ. Vortex shedding and frequency selection in flapping flight. Journal of Fluid Mechanics. 2000; 410:323-341
Young J, Lai J. Oscillation frequency and amplitude effects on the wake of a plunging airfoil. AIAA Journal. 2004; 42:2042-2052
Guerrero J. Aerodynamic performance of cambered heaving airfoils. AIAA Journal. 2010; 48:2694-2698
Parker K, Soria J, von Ellenrieder K. Thrust measurements from a finite-span flapping wing. AIAA Journal. 2007; 45:58-70
Parker K, von K, Soria J. Flow structures behind a heaving and pitching finite-span wing. Journal of Fluid Mechanics. 2003; 490:129-138
Dong H, Mittal R, Najjar FM. Wake topology and hydrodynamic performance of low-aspect-ratio flapping foils. Journal of Fluid Mechanics. 2006; 566:309-343
Blondeaux P, Fornarelli F, Guglielmini L, Triantafyllou MS, Verzicco R. Numerical experiments on flapping foils mimicking fish-like locomotion. Physics of Fluids. 2005; 17:113601-113612
Young J, Lai J. Mechanisms influencing the efficiency of oscillating airfoil propulsion. AIAA Journal. 2007; 45:1695-1702