2. Numerical method

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), nanoair-vehicles (NAVs), and autonomous-underwater-vehicles (AUVs) is inspired by flapping wings or oscillating fins [3–8], energy harvesting applications [9], and even robotic extraterres-

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

Many researchers [13–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 dimen-

St <sup>¼</sup> f h

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 (f h) 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

<sup>k</sup> <sup>¼</sup> <sup>π</sup>fc

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

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,

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

over the wing/fin chord compared to the period of motion and is defined as,

fishes, independently of their possible practical applications.

this was done to study leading edge vortex shedding dependency.

<sup>U</sup> (1)

<sup>U</sup> (2)

by swimming animals, significant progress was achieved [12].

trial exploring missions [10].

120 Flight Physics - Models, Techniques and Technologies

sionless parameter defined as,

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 \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} = \frac{-\nabla p}{\rho} + \nu \nabla^2 \mathbf{u} \qquad \text{for} \quad \mathbf{x} \in \mathcal{D}, \qquad t \ge 0,\tag{3}$$

$$
\nabla \cdot \mathbf{u} = 0 \qquad \qquad \text{for} \quad \mathbf{x} \in \mathcal{D}, \qquad t \ge 0,\tag{4}
$$

with the following boundary conditions and initial conditions,

$$B(\mathbf{u}, p) = \mathbf{g} \qquad \text{for} \quad \mathbf{x} \in \partial \mathcal{D}, \qquad t \ge 0,\tag{5}$$

$$
\stackrel{\circ}{\mathbf{Q}}(\mathbf{x},0) = \mathbf{q}\_0(\mathbf{x}) \qquad \text{for} \quad \mathbf{x} \in \mathcal{D}, \qquad t = 0. \tag{6}
$$

In this IBVP, the vector x = (x, y, z) contains the Cartesian coordinates in physical space <sup>P</sup> <sup>¼</sup> <sup>P</sup>(x, <sup>y</sup>, <sup>z</sup>, <sup>t</sup>), <sup>D</sup> is a bounded domain in <sup>P</sup> <sup>∈</sup> <sup>R</sup><sup>N</sup> (where <sup>N</sup> is the number of space dimensions), ∂D are the boundaries of the bounded domain 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), Q ∘ is the initial condition, and q<sup>0</sup> 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 \nabla \mathbf{u} = \frac{-\nabla p}{\rho} + \nu \nabla^2 \mathbf{u} \quad \text{for} \quad \mathbf{x} \in \mathcal{D}\_{\prime} \quad t \ge 0,\tag{7}$$

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

with the following boundary and initial conditions,

$$B(\mathbf{u}, p) = \mathbf{g} \tag{9}$$

$$f \mathbf{e} \quad \qquad \qquad \qquad \qquad \qquad \qquad \quad f \mathbf{e} \in \partial \mathcal{D}, \qquad t \ge 0,\tag{9}$$

$$\nabla \cdot \mathbf{u} = 0 \qquad\qquad\qquad\text{for}\quad \mathbf{x} \in \partial \mathcal{D}, \qquad t \ge 0,\tag{10}$$

$$\stackrel{\circ}{\mathbf{Q}}(\mathbf{x},0) = \mathfrak{q}\_0(\mathbf{x}) \qquad\qquad\qquad\text{for}\quad \mathbf{x} \in \mathcal{D}, \qquad t = 0. \tag{11}$$

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 D, with prescribed boundary conditions ∂D and given initial conditions 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 G<sup>g</sup> that completely cover the domain D that is being modeled in physical space P and overlap where they meet. In this newly generated overlapping grid system 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} - \dot{\mathbf{G}} \right) \cdot \nabla \right] \mathbf{u} = \frac{-\nabla p}{\rho} + \nu \nabla^2 \mathbf{u} \quad \text{for} \quad \mathbf{x} \in \mathcal{D}, \quad t \ge 0,\tag{12}
$$

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

where 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}(\mathbf{x}\_{\text{wall}},t) = \dot{\mathbf{G}}\left(\mathbf{x}\_{\text{wall}},t\right), \quad \text{where} \quad \mathbf{x}\_{\text{wall}} \in \partial \mathcal{D}\_{\text{wall}}(t). \tag{14}$$

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.

Bð Þ¼ u; p g for x ∈∂D, t ≥ 0, (9)

∇ � u ¼ 0 for x ∈∂D, t ≥ 0, (10)

ð Þ¼ x; 0 q0ð Þx for x∈ D, t ¼ 0: (11)

(Eqs. (9)–(11)). Eqs. (7)–(11) are solved in logically

u for x ∈ D, t ≥ 0, (12)

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 D, with prescribed boundary

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

The structured overlapping grids method consists in generating a set of conforming structured component grids G<sup>g</sup> that completely cover the domain D that is being modeled in physical space P and overlap where they meet. In this newly generated overlapping grid system 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

<sup>þ</sup> <sup>ν</sup>∇<sup>2</sup>

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

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

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

þ ∇u � u<sup>x</sup> þ ∇v � u<sup>y</sup> þ ∇w � u<sup>z</sup> ¼ 0 for x ∈ D, t ≥ 0, (13)

xwall ð Þ ; t , where xwall ∈ ∂Dwallð Þt : (14)

∘

Q ∘

122 Flight Physics - Models, Techniques and Technologies

conditions ∂D and given initial conditions Q

component grids as follows,

yield different accuracies.

where G �

∂u

∇<sup>2</sup>p ρ

<sup>∂</sup><sup>t</sup> <sup>þ</sup> <sup>u</sup> � <sup>G</sup> � � �

should impose the velocity at the wall as follows,

u xwall ð Þ¼ ; t G

�

� ∇

<sup>u</sup> <sup>¼</sup> �∇<sup>p</sup> ρ

h i

finite-difference approximations 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.
