Wake Topology and Aerodynamic Performance of Heaving Wings

Abstract

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.

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 evolutionary-optimization process (biomimetics). Such applications may include drag reduction, noise reduction, and flow control by using feather-like structures [1] and flippers tubercles [2]; 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 [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 fish-inspired 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,

$$\mathit{St}=\frac{\mathit{fh}}{U}$$ E1

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,

$$k=\frac{\mathit{\pi fc}}{U}$$ E2

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 ≤ 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 initial-boundary-value problem (IBVP) for the laminar incompressible Navier-Stokes equations can be written as

$$\frac{\partial \mathrm{u}}{\partial t}+\mathbf{u}\cdot \mathbf{\nabla }\mathbf{u}=\frac{-\nabla p}{\rho }+\nu {\mathbf{\nabla }}^{\mathbf{2}}\mathbf{u}\mathit{\text{for}}\mathbf{x}\in \mathcal{D},t\ge 0,$$ E3

$$\mathbf{\nabla }\cdot \mathbf{u}=0\mathit{\text{for}}\mathbf{x}\in \mathcal{D},t\ge 0,$$ E4

with the following boundary conditions and initial conditions,

$$B\left(\mathbf{u},p\right)=\mathbf{g}\mathit{\text{for}}\mathbf{x}\in \partial \mathcal{D},t\ge 0,$$ E5

$$\stackrel{\circ }{\mathbf{Q}}\left(\mathit{x},0\right)={\mathbf{q}}_0\left(x\right)\mathit{\text{for}}\mathbf{x}\in \mathcal{D},t=0.$$ E6

In this IBVP, the vector x = (x, y, z) contains the Cartesian coordinates in physical space $\mathcal{P}=\mathcal{P}$ (x, y, z, t), $\mathcal{D}$ is a bounded domain in $\mathcal{P}\in {\mathfrak{R}}^{\mathbb{N}}$ (where $\mathbb{N}$ is the number of space dimensions), $\text{∂}\mathcal{D}$ are the boundaries of the bounded domain $\mathcal{D}$ , t is the physical time, u is a vector containing the velocity components (u, v, w), the scalar p is the pressure, and the constants ν and ρ are the kinematic viscosity and density. In Eq. (5), B is a boundary operator (that can leave to a Dirichlet or Neumann boundary condition), and g is the boundary condition input value. In Eq. (6), ${\displaystyle \stackrel{\circ }{\mathbf{Q}}}$ is the initial condition, and q ₀ is the input value of the initial conditions (which can be a uniform or nonuniform field).

An alternative formulation of the system of Eqs. (3)–(6), called the velocity-pressure formulation, can be written as follows,

$$\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \mathbf{\nabla }\mathbf{u}=\frac{-\nabla p}{\rho }+\nu {\mathbf{\nabla }}^{\mathbf{2}}\mathbf{u}\mathit{\text{for}}\mathbf{x}\in \mathcal{D},t\ge 0,$$ E7

$$\frac{{\mathbf{\nabla }}^{2}p}{\rho }+\nabla u\cdot {\mathbf{u}}_x+\nabla v\cdot {\mathbf{u}}_y+\nabla w\cdot {\mathbf{u}}_z=0\mathit{\text{for}}\mathbf{x}\in \mathcal{D},t\ge 0,$$ E8

with the following boundary and initial conditions,

$$B\left(\mathbf{u},p\right)=\mathit{g}\mathit{\text{for}}\mathbf{x}\in \partial \mathcal{D},t\ge 0,$$ E9

$$\nabla \cdot \mathbf{u}=0\mathit{\text{for}}\mathbf{x}\in \partial \mathcal{D},t\ge 0,$$ E10

$$\stackrel{\circ }{\mathbf{Q}}\left(\mathbf{x},0\right)={\mathit{q}}_0\left(\mathbf{x}\right)\mathit{\text{for}}\mathbf{x}\in \mathcal{D},t=0.$$ E11

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 $\mathcal{D}$ , with prescribed boundary conditions $\text{∂}\mathcal{D}$ and given initial conditions ${\displaystyle \stackrel{\circ }{\mathbf{Q}}}$ (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 ${\mathcal{G}}_g$ that completely cover the domain $\mathcal{D}$ that is being modeled in physical space $\mathcal{P}$ and overlap where they meet. In this newly generated overlapping grid system $\mathbb{G}$ , 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,

$$\frac{\partial \mathbf{u}}{\partial t}+\left[\left(\mathbf{u}-{\displaystyle \stackrel{\cdot }{\mathbf{G}}}\right)\cdot \nabla \right]\mathbf{u}=\frac{-\nabla p}{\rho }+\nu {\nabla }^2\mathbf{u}\mathit{\text{for}}\mathbf{x}\in \mathcal{D},t\ge 0,$$ E12

$$\frac{{\nabla }^{2}p}{\rho }+\nabla u\cdot {\mathbf{u}}_x+\nabla v\cdot {\mathbf{u}}_y+\nabla w\cdot {\mathbf{u}}_z=0\mathit{\text{for}}\mathbf{x}\in \mathcal{D},t\ge 0,$$ E13

where ${\displaystyle \stackrel{\cdot }{\mathbf{G}}}$ 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,

$$\mathbf{u}\left({\mathbf{x}}_{\mathit{\text{wall}}},t\right)=\stackrel{\cdot }{\mathbf{G}}\left({\mathbf{x}}_{\mathit{\text{wall}}},t\right),\mathit{\text{where}}{\mathbf{x}}_{\mathit{\text{wall}}}\in \partial {\mathcal{D}}_{\mathit{\text{wall}}}\left(t\right).$$ E14

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

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 (LH-WT and RH-WT, 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.

Figure 1.

Left: computational domain layout in the xy plane. Right: computational domain layout in the zy plane. The figure is not to scale.

In all the cases studied, a rectangular wing with an aspect ratio AR equal to 2 was used. The cross-section 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), ${\text{∂}}_{\widehat{\mathbf{n}}}p=0$ ), and the top, bottom, and right boundaries of the BG are outflow boundaries (velocity extrapolated from the interior points). In Figure 1 (right), all the boundaries of the BG correspond to outflow conditions. On the wing surface (which is a moving body), we impose a no-slip boundary condition for moving walls ( $\mathbf{u}={{\displaystyle \stackrel{\cdot }{\mathbf{G}}}}_x,{{\displaystyle \stackrel{\cdot }{\mathbf{G}}}}_y,{{\displaystyle \stackrel{\cdot }{\mathbf{G}}}}_z$ ). The rest of the boundaries is interpolation boundaries, where we used a nonconservative Lagrange interpolation scheme. The Reynolds number (defined as Re = U × c/ν) is equal to 250 for all the simulations.

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 y axis in Figure 1 ) and according to the following function,

$$y\left(t\right)=h\times \mathit{sin}\left(2\times \pi \times f\times t+\phi \right)$$ E15

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 $\mathbb{G}$ 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 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 $\mathbb{G}$ used in this study. In Figure 2 , the background grid is shown in images a, b, and c (rectangular grid). The wing grid is made up of three component grids (image d), namely, center section, left wing grid tip-section, and right wing grid tip-section.

Grid ${\mathbb{G}}_g$	BG	WG-CS	WG-TS	GSR	1NW
${\mathbb{G}}_1$	161 × 121 × 101	221 × 121 × 41	81 × 61 × 41	1	0:001 × 2c
${\mathbb{G}}_2$	161 × 121 × 101	201 × 101 × 31	61 × 51 × 31	2	0:002 × 2c
${\mathbb{G}}_3$	161 × 121 × 101	161 × 81 × 31	51 × 41 × 31	4	0:004 × 2c

Table 1.

Grid dimensions used for the grid dependence study.

BG stands for background grid, WG-CS stands for wing grid center-section, WG-TS stands for wing grid tip-section (only one tip-section), GSR stands for grid refinement spacing ratio, and 1NW stands for position of the first node normal to the wall.

Figure 2.

Typical overlapping grid system $\mathbb{G}$ used in this study. (A) Side view. (B) Front view. (C) Bottom view. (D) Perspective view.

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,

$$c_l=\frac{L}{\frac{1}{2}\rho U^{2}A},c_d=\frac{D}{\frac{1}{2}\rho U^{2}A},c_t=\frac{T}{\frac{1}{2}\rho U^{2}A}$$ E16

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 $\overline{c_d}$ and the root mean square of the lift coefficient $c_l^{\mathit{rms}}$ for the grid dependence study. These values were used to compute the observed order of convergence p_obs , the extrapolated value of the observed quantity at zero grid spacing χ _{spacing = 0}, the fine grid convergence index GCI₁₂ and GCI₂₃, and the constancy of GCI ₂₃ = r^p_obs  × GCI ₁₂.

Grid ${\mathbb{G}\mathrm{g}}_g$	GSR	$\overline{c_d}$	$c_l^{\mathit{rms}}$
${\mathbb{G}}_1$	1	0.057892	0.905822
${\mathbb{G}}_2$	2	0.058476	0.897486
${\mathbb{G}}_3$	4	0.060914	0.865326

Table 2.

Observed values of $\overline{c_d}$ and $c_l^{\mathit{rms}}$ for the grid dependence study.

Based on the outcome of the GCI study (refer to Table 3 ), $\overline{c_d}$ at zero grid spacing is estimated to be 0.057708, with an error band of 0.397203% for grid ${\mathbb{G}}_1$ and an error band of 1.641617% for grid ${\mathbb{G}}_2$ , whereas $c_l^{\mathit{rms}}$ at zero grid spacing is estimated to be 0.908738 with an error band of 0.402503% for grid ${\mathbb{G}}_1$ and an error band of 1.567203% for grid ${\mathbb{G}}_2$ . The constancy of GCI ₂₃ = r^p_obs  × GCI ₁₂ indicates that the results for grids ${\mathbb{G}}_1$ and ${\mathbb{G}}_2$ are within the asymptotic range of convergence. Finally, the observed order of convergence p_obs represents a direct indication of the accuracy of the numerical method. In this study, we obtained a value of p_obs close to 2 for both quantities of interest, which indicate that the method is second-order accurate in space and time.

Outcome	$\overline{{\mathbf{c}}_{\mathbf{d}}}$	$c_l^{\mathbf{rms}}$
pobs	2.061650	1.947838
χ spacing = 0	0.057708	0.908738
GCI12 (%)	0.397203	0.402503
GCI23 (%)	1.641617	1.567203
(GCI23/GCI12) × (1/rpobs )	0.990012	1.009249

Table 3.

GCI study results.

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 ${\mathbb{G}}_1$ and ${\mathbb{G}}_2$ (refer to Table 1 ) are illustrated in Figure 3 . In Figure 3 , we use the iso-surfaces of vorticity magnitude (∣ω∣-criterion), and the iso-surfaces of Q-criterion [27, 28], to capture the vortices and their corresponding cores. As depicted in this figure, there are no discernible differences between the solutions, indicating that both grids are adequate for vortical structures resolution. It can also be seen that the ∣ω∣-criterion, although capable of capturing the general vortical structures, has the disadvantage of also showing the shear layers near the wing surface and between the vortices. On the other hand, the Q-criterion shows the details of the vortical structures more clearly as it does not show the shear layers close to the wing (as illustrated in Figure 4 ). Apart from this, both methods provide nearly identical structures in the far wake. Based on these qualitative results, we chose the Q-criterion as the main criterion for wake topology characterization.

Figure 3.

Vortex topology at the beginning of the upstroke (t = 5.0). Iso-surfaces of ∣ω∣-criterion are shown in light gray, and iso-surfaces of Q-criterion are shown in dark gray. Simulation parameters: Re = 250, St = 0.3, and k = 1.570795. A corresponds to overlapping grid system ${\mathbb{G}}_1$ and B corresponds to overlapping grid system ${\mathbb{G}}_2$ .

Figure 4.

Shear layers close to the wing at the beginning of the upstroke (t = 5.0). Left column corresponds to iso-surfaces of Q-criterion, and the right column corresponds to iso-surfaces of ∣ω∣-criterion. Simulation parameters: Re = 250, St = 0.3, and k = 1.570795.

Summarizing the quantitative and qualitative results previously presented, we can conclude that the solutions obtained by using the overlapping grid systems ${\mathbb{G}}_1$ and ${\mathbb{G}}_2$ are grid independent. Taking into account the computational resources available, CPU time restrictions, and solution accuracy, ${\mathbb{G}}_2$ with a value of 1 NW equal to 0.001 × 2c is used as the base grid to perform all further computations. In the case of a smaller or bigger computational domain, the grid dimensions are scaled in order to keep the same grid spacing as for this domain.

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 $\overline{c_t}$ is the average thrust coefficient and ${\widehat{c}}_l$ is the maximum lift coefficient (which was measured during the downstroke). By inspecting these results, we observe that as we increase St and k, the values of $\overline{c_t}$ and ${\widehat{c}}_l$ also increase. This result is expected because as we increase St and k (therefore, the oscillating frequency and heaving amplitude), the vertical velocity of the wing is higher. Hence, the forces excerpted on the wing’s surface are larger. We also observe two different behaviors of the aerodynamic forces for high and low reduced frequencies k values. Hence, it seems that for heaving wings, the oscillating frequency plays an important role in the vortex generation and shedding and, henceforth, on the aerodynamic forces. These frequency dependence observations are similar to those of Wang [29], Young and Lai [30], and Guerrero [31], but here, we extend them to three-dimensional cases.

Case number	St	k	$\overline{c_t}$	${\widehat{c}}_l$	Regime
3DH-1	0.15	1.570795	−0.122675	0.560613	Drag production
3DH-2	0.15	0.785397	−0.111949	0.422249	Drag production
3DH-3	0.20	1.570795	−0.106101	0.784628	Drag production
3DH-4	0.20	0.785397	−0.087877	0.597156	Drag production
3DH-5	0.25	1.570795	−0.086136	1.040502	Drag production
3DH-6	0.25	0.785397	−0.068118	0.798169	Drag production
3DH-7	0.30	1.570795	−0.058294	1.238780	Drag production
3DH-8	0.30	0.785397	−0.052475	1.024220	Drag production
3DH-9	0.35	1.570795	0.022962	1.651670	Thrust production
3DH-10	0.35	0.785397	−0.036191	1.257870	Drag production
3DH-11	0.40	1.570795	0.061370	2.030980	Thrust production
3DH-12	0.40	0.785397	0.008296	1.539740	Thrust production
3DH-13	0.45	1.570795	0.143680	2.417032	Thrust production
3DH-14	0.45	0.785397	0.010296	1.891345	Thrust production
3DH-15	0.50	1.570795	0.199366	2.944020	Thrust production
3DH-16	0.50	0.785397	0.024617	2.150710	Thrust production

Table 4.

Simulation results for the pure heaving parametric study (positive $\overline{c_t}$ indicates thrust production whereas negative $\overline{c_t}$ indicates drag production).

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 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.

Figure 5.

Vortex wake topology at the beginning of the upstroke for case 3DH-11 (t = 5.0). Heaving parameters: Re = 250, St = 0.4, and k = 1.570795 (thrust producing wake).

Figure 6.

Vortex wake topology at the beginning of the upstroke for case 3DH-11 (t = 5.0). Heaving parameters: Re = 250, St = 0.4, and k = 1.570795 (thrust producing wake).

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).

Figure 7.

Vortex wake topology at the beginning of the upstroke for case 3DH-3 (t = 5.0). Heaving parameters: Re = 250, St = 0.2, and k = 1.570795 (drag-producing wake).

Figure 8.

Wake comparison of a drag-producing case (A) and thrust producing case (B). The wake height was measured approximately at 3 × c behind the trailing edge.

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. [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 $\overline{c_t}$ and maximum lift coefficient ${\widehat{c}}_l$ . It seems that the reduced frequency k (and hence the oscillating frequency) plays an important role in the vortex generation and shedding frequency, therefore, on the aerodynamic forces. Thus, LEV convection and separation introduce a frequency dependence into the results. This provides a mechanism of optimal selection of heaving frequency (in the sense of propulsive efficiency) as discussed by Wang [29], Guerrero [31], and Young and Lai [36]. It is worth mentioning that the results presented in the previous references were obtained for two-dimensional airfoils. The results presented in this manuscript extend these observations to three-dimensional wings.

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.

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.