Grid dimensions used for the grid dependence study.
Abstract
Simulating the threedimensional flow features generated by heaving wings constitutes a great challenge due to the computational effort required to compute the complex threedimensional 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 NavierStokes equations are solved on moving overlapping structured grids using a secondorder accurate in space and time finitedifference numerical method. The numerical simulations are performed at a Reynolds number of Re = 250 and at different values of Strouhal number and heaving frequency.
Keywords
 heaving wings
 aerodynamic performance
 wake topology
 Strouhal number
 overlapping grids
1. Introduction
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 evolutionaryoptimization process (biomimetics). Such applications may include drag reduction, noise reduction, and flow control by using featherlike structures [1] and flippers tubercles [2]; the development of new propulsion/lift generation systems for microairvehicles (MAVs), nanoairvehicles (NAVs), and autonomousunderwatervehicles (AUVs) is inspired by flapping wings or oscillating fins [3, 4, 5, 6, 7, 8], energy harvesting applications [9], and even robotic extraterrestrial exploring missions [10].
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 fishinspired mechanisms achieved disappointingly low propulsive efficiencies [11]. It was only through a deeper understanding of the vorticity and wake produced by swimming animals, significant progress was achieved [12].
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
where
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 finitespan rigid wings undergoing pure heaving motion. The laminar incompressible NavierStokes equations are numerically approximated, and all unsteady, viscous, and threedimensional 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 [16], Henshaw and Petersson [17], and Chesshire and Henshaw [18].
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 velocitypressure 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
The structured overlapping grids method consists in generating a set of conforming structured component grids
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
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 semiimplicit multistep method, where the viscous terms are approximated using the CrankNicolson scheme, and the convective terms are approximated using the AdamsBashforth/AdamsMoulton predictorcorrector approach (where the velocity is advanced in time by using a secondorder AdamsBahforth predictor step, followed by a secondorder AdamsMoulton 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 secondorder accurate in space and time numerical scheme on structured overlapping grids.
To assemble the overlapping grid system and solve the laminar incompressible NavierStokes equations in their velocitypressure 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 NavierStokes 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 NewtonKrylov iterative method, in combination with a suitable preconditioner and in parallel computational architectures.
3. Geometry, boundary conditions, initial conditions, and wing kinematics
In
Figure 1
, 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
Figure 1
(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 crosssection heaves in the vertical direction (or
where
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
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
Table 1
. In
Table 1
, 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
Grid

BG  WGCS  WGTS  GSR  1NW 


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
In
Table 2
, we present the average drag coefficient
Grid

GSR 




1  0.057892  0.905822 

2  0.058476  0.897486 

4  0.060914  0.865326 
Based on the outcome of the GCI study (refer to
Table 3
),
Outcome 




2.061650  1.947838 

0.057708  0.908738 
GCI_{12} (%)  0.397203  0.402503 
GCI_{23} (%)  1.641617  1.567203 
(GCI_{23}/GCI_{12}) × (1/ 
0.990012  1.009249 
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
Summarizing the quantitative and qualitative results previously presented, we can conclude that the solutions obtained by using the overlapping grid systems
5. Simulation results
Hereafter, we carry out a comprehensive parametric study to assess the wake signature and aerodynamic performance of heaving rigid wings. In
Table 4
, we present the kinematics parameters governing the heaving motion (described by Eq. (15)). In
Table 4
, we also present the quantitative results obtained, where
Case number  St  k 


Regime 

3DH1  0.15  1.570795  −0.122675  0.560613  Drag production 
3DH2  0.15  0.785397  −0.111949  0.422249  Drag production 
3DH3  0.20  1.570795  −0.106101  0.784628  Drag production 
3DH4  0.20  0.785397  −0.087877  0.597156  Drag production 
3DH5  0.25  1.570795  −0.086136  1.040502  Drag production 
3DH6  0.25  0.785397  −0.068118  0.798169  Drag production 
3DH7  0.30  1.570795  −0.058294  1.238780  Drag production 
3DH8  0.30  0.785397  −0.052475  1.024220  Drag production 
3DH9  0.35  1.570795  0.022962  1.651670  Thrust production 
3DH10  0.35  0.785397  −0.036191  1.257870  Drag production 
3DH11  0.40  1.570795  0.061370  2.030980  Thrust production 
3DH12  0.40  0.785397  0.008296  1.539740  Thrust production 
3DH13  0.45  1.570795  0.143680  2.417032  Thrust production 
3DH14  0.45  0.785397  0.010296  1.891345  Thrust production 
3DH15  0.50  1.570795  0.199366  2.944020  Thrust production 
3DH16  0.50  0.785397  0.024617  2.150710  Thrust production 
In
Table 4
, 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 3DH11 in Table 4 ). These figures show that the downstream wake (or far field wake) of this case consists of two sets of doughnutshaped 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 wingtip vortices or WTV (one WTV on the left wingtip or LHWTV and another WTV on the right wingtip or RHWTV). 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 doughnutshaped 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 wingtip 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 dragproducing case (which corresponds to case 3DH3 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 NavierStokes equations were solved in their velocitypressure formulation using a secondorder accurate in space and time finitedifference 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 freestream flow, whereas for thrust producing configurations, the angle of inclination of the vortex rings is in the same direction of their travel, and for dragproducing 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 wingtip 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. [32] and Von Ellenrieder et al. [33] and the numerical simulations of Dong et al. [34] and Blondeaux et al. [35].
From the force measurement study, two different behaviors were observed for the average thrust coefficient
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.
Acknowledgments
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.
References
 1.
Favier J, Dauptain A, Basso D, Bottaro A. Passive separation control using a selfadaptive hairy coating. Journal of Fluid Mechanics. 2009; 627 :451483  2.
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 :107114  3.
Jones KD, Platzer MF. Experimental investigation of the aerodynamic characteristics of flappingwing micro air vehicles. In: AIAA Paper 20030418CP. 2003. pp. 111  4.
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 :777785  5.
McMichael JM, Francis MS. Micro Air Vehicles toward a New Dimension in Flight. Technical report, Defense Advanced Research Projects Agency (DARPA). 1997  6.
Guerrero J. Wake signature and aerodynamic performance of finitespan root flapping rigid wings. Journal of Bionic Engineering. 2010; 7 :S109S122  7.
Guerrero J, Pacioselli C, Pralits JO, Negrello F, Silvestri P, Lucifredi A, Bottaro A. Preliminary design of a smallsized flapping UAV I. Aerodynamic performance and static longitudinal stability. Meccanica. 2016; 51 :13431367  8.
Negrello F, Silvestri P, Lucifredi A, Guerrero J, Bottaro A. Preliminary design of a smallsized flapping UAV II. Kinematic and structural aspects. Meccanica. 2016; 51 :13691385  9.
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  10.
Michelson R, Naqvi M. Extraterrestrail flight (entompterbased mars surveyor). Von Karman Institute for Fluid Dynamics. Lecture Series. November 2003; 2003 :2428  11.
Triantafyllou GS, Triantafyllou MS. An efficient swimming machine. Scientific American. 1995; 272 :4048  12.
Rohr J, Fish F. Strouhal number and optimization of swimming by odontocete cetaceans. The Journal of Experimental Biology. 2004; 207 :16331642  13.
Triantafyllou MS, Triantafyllou GS, Gopalkrishnan R. Wake mechanics for thrust generation in oscillating foils. Physics of Fluids. 1991; 3 :28352837  14.
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 :20712076  15.
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 :707711  16.
Henshaw WD. A fourthorder accurate method for the incompressible NavierStokes equations on overlapping grids. Journal of Computational Physics. 1994; 113 :1325  17.
Henshaw WD, Petersson NA. A splitstep scheme for the incompressible NavierStokes equations. Numerical Solutions of Incompressible Flows. Preprint UCRLJC144. 2003. pp. 114  18.
Chesshire G, Henshaw W. Composite overlapping meshes for the solution of partial differential equations. Journal of Computational Physics. 1990; 90 :164  19.
Petersson N. Stability of pressure boundary conditions for stokes and navierstokes equations. Journal of Computational Physics. 2001; 172 :4070  20.
Sani RL, Shen J, Pironneau O, Gresho PM. Pressure boundary condition for the timedependent incompressible navierstokes equations. International Journal for Numerical Methods in Fluids. 2006; 50 :673682  21.
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. 1220  22.
Vinokur M. Conservation equations of gasdynamics in curvilinear coordinate systems. Journal of Computational Physics. 1974; 14 :105125  23.
Schiesser WE. The Numerical Method of Lines: Integration of Partial Differential Equations. San Diego, California: Academic Press; 1991. pp. 1326  24.
Wesseling P. Principles of Computational Fluid Dynamics. Berlin, Heidelberg: Springer; 2001. pp. 1644  25.
Roache P. Verification and Validation in Computational Science and Engineering. Albuquerque, New Mexico: Hermosa Publishers; 1998. pp. 1464  26.
Roache P. Quantification of uncertainty in computational fluid dynamics. Annual Review of Fluid Mechanics. 1997; 29 :123160  27.
Jeong J, Hussain F. On the identification of a vortex. Journal of Fluid Mechanics. 1995; 285 :6994  28.
Haimes R, Kenwright D. On the velocity gradient tensor and fluid feature extraction. In: AIAA Paper 1999–3288CP. 1999. pp. 110  29.
Wang ZJ. Vortex shedding and frequency selection in flapping flight. Journal of Fluid Mechanics. 2000; 410 :323341  30.
Young J, Lai J. Oscillation frequency and amplitude effects on the wake of a plunging airfoil. AIAA Journal. 2004; 42 :20422052  31.
Guerrero J. Aerodynamic performance of cambered heaving airfoils. AIAA Journal. 2010; 48 :26942698  32.
Parker K, Soria J, von Ellenrieder K. Thrust measurements from a finitespan flapping wing. AIAA Journal. 2007; 45 :5870  33.
Parker K, von K, Soria J. Flow structures behind a heaving and pitching finitespan wing. Journal of Fluid Mechanics. 2003; 490 :129138  34.
Dong H, Mittal R, Najjar FM. Wake topology and hydrodynamic performance of lowaspectratio flapping foils. Journal of Fluid Mechanics. 2006; 566 :309343  35.
Blondeaux P, Fornarelli F, Guglielmini L, Triantafyllou MS, Verzicco R. Numerical experiments on flapping foils mimicking fishlike locomotion. Physics of Fluids. 2005; 17 :113601113612  36.
Young J, Lai J. Mechanisms influencing the efficiency of oscillating airfoil propulsion. AIAA Journal. 2007; 45 :16951702