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 St is a dimensionless parameter defined as,
where f is the flapping frequency, h is the peak to peak amplitude of the flapping stroke, and U is the forward velocity. This definition describes a ratio between the oscillating speed (fh) and the forward speed. Another dimensionless parameter that characterizes the aero/hydrodynamic performance and wake signature of flying and swimming animals is the reduced frequency k, which is a measure of the residence time of a vortex (or a particle) convecting over the wing/fin chord compared to the period of motion and is defined as,
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 ≤ St ≤ 0.5, and for two different reduced frequency values, one corresponding to high frequency and the other one to low frequency, this was done to study leading edge vortex shedding dependency.
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 (u, v, w, p), the governing equations of the initialboundaryvalue problem (IBVP) for the laminar incompressible NavierStokes equations can be written as
with the following boundary conditions and initial conditions,
In this IBVP, the vector x = (x, y, z) contains the Cartesian coordinates in physical space
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, c is the wing’s chord, and h is the heaving amplitude of the heaving wing. The background grid (BG) extends 4.0 × c away from the wing’s leading edge (LE), 10.0 × c away from the wing’s trailing edge (TE), 2.0 × c away from the left and right wing’s tips (LHWT and RHWT, respectively), and 4.0 × c + h away from the point of maximum thickness of the upper and lower surfaces. This overlapping grid system layout corresponds to the instant, when the wing is in the mid position of the heaving cycle (as illustrated in Figure 1 ). In this manuscript, all the base units are expressed in the international system.
In all the cases studied, a rectangular wing with an aspect ratio AR equal to 2 was used. The crosssection of the wing is an ellipse, with a corresponding major axis a = 0.25 and a minor axis b = 0.025. Therefore, the wing’s chord c is equal to 2 × a = 0.5.
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 (u = (1.0, 0.0, 0.0),
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 y axis in Figure 1 ) and according to the following function,
where y(t) is the heaving motion (and is defined positive upwards), h is the heaving amplitude, f is the heaving oscillating frequency, φ is the phase angle of the heaving motion (0 in this case), and t is the physical time.
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 St = 0.3 and at a reduced frequency k = 1.570795, and unsteadiness was observed to disappear typically after 4 to 5 cycles of wing heaving motion, and further calculations show negligible nonperiodicity. Each simulation was checked for acceptable iterative convergence.
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 r is equal to 2. In
Figure 2
, we illustrate a typical overlapping grid system
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 
BG stands for background grid, WGCS stands for wing grid centersection, WGTS stands for wing grid tipsection (only one tipsection), GSR stands for grid refinement spacing ratio, and 1NW stands for position of the first node normal to the wall.
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 T instead of drag force D. The thrust force is equal in magnitude but opposite in direction to the drag force (T = − D). To quantify the unsteady aerodynamics performance, we computed the lift coefficient c_{l} , the drag coefficient c_{d} , and thrust coefficient c_{t} as follows,
In Eq. (16), the lift force L and the thrust force T (where T = − D) are computed by integrating the viscous and pressure forces over the wing surface.
In
Table 2
, we present the average drag coefficient
Based on the outcome of the GCI study (refer to
Table 3
),
Outcome 



pobs  2.061650  1.947838 
χ _{ spacing = 0}  0.057708  0.908738 
GCI_{12} (%)  0.397203  0.402503 
GCI_{23} (%)  1.641617  1.567203 
(GCI_{23}/GCI_{12}) × (1/r_{pobs} )  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
In Table 4 , we can read that for values of St < 0.30, the heaving wing produces drag; for values of 0.30 < St < 0.35, the heaving wing produces little or no drag (or thrust), whereas for values of St > 0.35, the heaving wing produces thrust. These results suggest that there is a range of St values, where the heaving wing generates thrust. The results also point at the presence of a combination of St and k values, where propulsive efficiency peaks. In general, the results are inline with the hypothesis that flying and swimming animals cruise at Strouhal numbers corresponding to a regime of vortex growth and shedding in which propulsion efficiency is high.
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 St and k values studied and the simplified wing geometry and heaving kinematics covered in this study, all the qualitative and quantitative results presented are in close agreement with the experimental observations of Rohr and Fish [12], Triantafyllou et al. [13], Nudds et al. [14], Taylor et al. [15], and Parker et al. [32]; this supports the hypothesis that “flying and swimming animals cruise at a Strouhal number tuned for high power efficiency” [15].
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.