Effects of mass and chordwise flexibility on 2D self-propelled flapping wings

Journal of Fluids and Structures 64 (2016) 46–66

Contents lists available at ScienceDirect

Journal of Fluids and Structures journal homepage: www.elsevier.com/locate/jfs

Effects of mass and chordwise flexibility on 2D self-propelled flapping wings Mathieu Olivier n, Guy Dumas Laboratoire de Mécanique des Fluides Numérique, Department of Mechanical Engineering, Laval University, Quebec City, Quebec, Canada G1V 0A6

a r t i c l e in f o

abstract

Article history: Received 5 January 2016 Received in revised form 1 April 2016 Accepted 13 April 2016

A self-propelled flexible flapping wing 2D numerical model undergoing a combined pitching and heaving motion is presented. Since such freely moving foil experiences zero net thrust, a definition of efficiency for this kind of problem is proposed and discussed against other formulations found in the literature. It is also shown that the deviation motion of wings such as that found in natural flyers is likely a consequence of the fluid– structure dynamics of the wings. The passive deviation motion observed in numerical simulations is either a consequence of a feathering mechanism referred to as rigid feathering or of the inertial displacement caused by the wing deformation. The effects of flexibility on the performance of the wing are also presented. It is found that flexibility may significantly enhance the efficiency in pressure-driven deformation cases. The rigid feathering mechanism is found to have an effect similar to that of the feathering caused by wing flexibility on the performances of pressure-driven deformation cases. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Flexible flapping wings Fluid–structure interaction Self-propulsion Deviation motion

1. Introduction While most studies of oscillating wing are performed in an imposed upstream flow, the motion of flying animals or manmade flying devices is not constrained; it is a consequence of the propulsive force. In the case of flapping wings, the resulting velocity is most likely to fluctuate since the propulsive mechanism is unsteady. This remains true unless the mass of the body on which the wings are attached, is much larger than the wing themselves. Recently, the study of self-propelled flapping wings seems to have gained prevalence over “clamped” flapping wings (Zhang et al., 2009, 2010; Spagnolie et al., 2010; Thiria and Godoy-Diana, 2010; Ramananarivo et al., 2011; Alben et al., 2012; Hua et al., 2013; Lee and Lee, 2013; Yeh and Alexeev, 2014; Zhu et al., 2014a, 2014b). Letting the wing free to move in the direction transverse to the strokes inherently introduces some passive mechanisms on the wing. Yet, passive mechanisms are known to be used by flying and swimming animals to minimize their energy consumption. For example, some insects, such as dragonflies, benefit from fluid torque to achieve effortless wing rotation during the pronation (see Wang, 2005). As another example, fishes that swim upstream behind an obstacle not only do benefit from a low velocity wake region that reduces the effort needed to stay put, but actually experience thrust. Indeed, Beal et al. (2006) demonstrated that even a dead fish can experience thrust when its body resonates with the vortices in the wake of a bluff body. Further experiments showed that it is also possible to apply this principle to a high-aspect-ratio foil n

Corresponding author. E-mail addresses: [email protected] (M. Olivier), [email protected] (G. Dumas).

http://dx.doi.org/10.1016/j.jfluidstructs.2016.04.002 0889-9746/& 2016 Elsevier Ltd. All rights reserved.

M. Olivier, G. Dumas / Journal of Fluids and Structures 64 (2016) 46–66

Nomenclature

c c CP Cp CT (d x , d y ) E e e˜ F′ f f⁎ h h⁎ h0 I′ M′

N′ ^ n P p

relative convective velocity chord length power coefficient pressure coefficient thrust coefficient displacement vector field Young modulus wing thickness Green–Lagrange strain force vector per unit span frequency reduced frequency heaving displacement normalized heaving amplitude heaving amplitude beam cross-section area moment of inertia per unit span internal bending moment in a beam per unit span, moment per unit span internal normal force in a beam per unit span unit surface normal vector power pressure field

Re S St T t Un V v ^ v

α

αEff χ˜ δP⁎

η μ ν ω

ϕ˜

ρ Σ τ θ θ0

θEff

47

Reynolds number surface Strouhal number thrust force time normalized velocity volume velocity vector field control surface velocity angle of attack effective angle of attack curvature of the deformed beam normalized flexibility efficiency dynamic viscosity kinematic viscosity vorticity field (z-component) phase shift density pressure-to-inertia ratio shear load pitching angle pitching amplitude effective pitching angle

behind a long D-cylinder. A similar study was also performed by Eldredge and Pisani (2008) using a viscous vortex particles method to model three rigid bodies connected with hinges. It was found, surprisingly, that the mechanism works equally well on multi-elements rigid bodies whether the hinges are linked with torsional springs or locked. Using the same numerical method, Wilson and Eldredge (2011) studied the physics of bodies. These bodies were either actively or passively controlled. In the case of passively controlled bodies, the motion was specified on certain hinges while the other hinges were linked with torsional springs. On the other hand, actively controlled bodies used prescribed kinematics on all hinges. It was found that some configurations with passively controlled bodies provide optimal swimming speed and efficiency. On a more fundamental perspective, von Ellenrieder et al. (2008) suggested that the choice of wing frequency in swimming and flying animals may be the result of a limit cycle process. That is, the flapping wing dynamics of such species would rely significantly on passive mechanisms. In that context, a device that propels itself using a specific flapping wing configuration (geometry and kinematics) without active control should be designed so that it naturally reaches a limit cycle that corresponds to an optimal regime. In this paper, mass and chordwise flexibility effects on self-propelled 2D flapping wings performances are investigated with respect to a formal dimensionless parametric space previously introduced (Olivier and Dumas, 2016). The imposed kinematics consists of a combined pitching and heaving motion while the wing is let free to move in the horizontal direction. The imposed pitching motion forces the wing to travel in a specific direction (a zero pitching amplitude can produce a hysteretic regime involving either forward or backward motion, see Zhang et al. (2009, 2010) and Spagnolie et al. (2010)). In cases of strong fluid–structure interactions (i.e., light wings with respect to fluid), it will be shown that two feathering mechanisms act simultaneously and strongly affect the dynamics of the wing and of the flow. On the other hand, when weak fluid–structure interactions occur (i.e., heavy wings with respect to fluid), the rigid wing acts mostly as if it was moving at constant velocity in the x-direction while the flexible ones exhibit a small deviation motion. Moreover, while self-propelled wings reach a terminal average velocity from which a reduced frequency can be computed, the circumstances under which the optimal reduced frequency of the corresponding constrained scenario is reached will be discussed. Section 2 introduces the proposed flexible flapping wing problem definition and the corresponding performance metrics. Section 3 presents the mathematical model that describes the incompressible flow and the elastic wing section as well as the numerical methods used to solve the fluid–structure coupled problem. Section 4 presents numerical verifications and discusses the validity of the numerical method. Finally, Section 5 reports the results and discussions of 2D self-propelled rigid and flexible flapping wings.

48

M. Olivier, G. Dumas / Journal of Fluids and Structures 64 (2016) 46–66

Fig. 1. Geometry and motion description of the flat plate wing considered here. The effective pitching angle θ Eff is also shown.

2. Problem definition The problem definition used for the study of self-propelled flexible flapping wings is very similar to that of the constrained flapping wings except that the wing is let free to move in the x-direction. The wing geometry used for this study is illustrated in Fig. 1. The wing is actually a flexible flat plate with massless and rigid rounded edges. The chord length c is given by the distance between points LE and TE on the undeformed wing and the wing thickness e is set to 0.01 c throughout the study. The wing motion is imposed at point P, which is located right behind the front rounded edge, and is defined by a combined heave and pitch motion given respectively by:

h = h0 cos(2πft ),

(1)

θ = θ0 sin(2πft ),

(2)

where h0 and θ0 are respectively the heaving and pitching amplitudes and f is the motion frequency. The initial x-velocity of the wing is zero. It is the thrust generated by the combined pitching and heaving motion that accelerates the wing up to a constant cycle-averaged velocity. In the case of flexible wings, the deformation results in an effective pitching angle which is defined in Fig. 1. The corresponding angle of attack and effective angle of attack are respectively given by:

⎛ ḣ ⎞ α(t ) = arc tan⎜ − ⎟ − θ (t ), vx ⎠ ⎝

(3)

⎛ ḣ ⎞ αEff (t ) = arc tan⎜ − ⎟ − θ Eff (t ), vx ⎠ ⎝

(4)

where θEff is the effective pitching angle that accounts for the wing deformation (see Fig. 1) and vx is the x-velocity of the wing. Table 1 Dimensionless parameters for the freely moving flexible flapping wing problem. Name Reynolds number

Definition

Ref =

Normalized heaving amplitude

h⁎ =

Pitching amplitude Pressure-to-inertia ratio

θ0

Flexibility

h0fc ν

h0 c

Σ=

ρf h0 ρs e

δP⁎ =

ρf (h0f )2c 3 EI′

M. Olivier, G. Dumas / Journal of Fluids and Structures 64 (2016) 46–66

49

2.1. Parametric space The dimensionless parameters of this problem are chosen based on those of the constrained flapping wing presented by Olivier and Dumas (2016). However, since the free-stream velocity is now a resulting metric rather than an actual parameter, the reference velocity used in the dimensionless parameters is chosen to be h0f , a characteristic velocity that scales with the heaving velocity. The resulting set of dimensionless parameters is summarized in Table 1 in which ν is the kinematic viscosity of the fluid, ρf and ρs are the density of the fluid and the solid respectively, E is the Young modulus of the wing, and I′ = e 3/12. When comparing this set with the one presented by Olivier and Dumas (2016), one notes that the Reynolds number Re is replaced by Ref which is simply the Reynolds number based on the vertical velocity and the reduced frequency f ⁎ = fc /U∞ is no more a parameter. In fact, the latter is now a result as it is being related to the a priori unknown freestream velocity. The parameter δP⁎, which is called flexibility, is chosen to provide an estimation of the plate deflection with respect to the fluid pressure load caused by the heaving motion. This parameter could alternatively have been based on a pressure load estimate caused by the pitching motion. However, because both the pitching and the heaving amplitudes remain fixed in this study, this choice does not have any impact on the interpretation of results. The parameter Σ, which is called pressure-to-inertia ratio, provides an estimation of pressure forces over inertia forces by taking into account the geometry of the structure since pressure forces act on surfaces while inertia acts on volumes. Here again, the load estimations are chosen to be based on the heaving motion. These last two parameters provide an insight into the strength of the fluid–structure interaction. For the interaction to be strong, two conditions must be met:


δP⁎ must be high so that the structure deformation cannot be ignored when computing the flow fields.
Σ must be high so that pressure forces cannot be ignored when computing the structure deformation. 2.2. Performance metrics The performance metrics of the self-propelled flapping wing problem differ from the case where the x-velocity is imposed. Following the fact that the wing is not constrained in the x-direction, the instantaneous horizontal reaction force on the driving mechanism, hence the actual instantaneous thrust, is zero. In fact, when horizontal aerodynamic forces are not null, the wing accelerates instead of providing a driving force. Thus, if the wing reaches a constant cycle-averaged velocity, namely a periodic regime, the acceleration is periodic and the cycle-averaged aerodynamic thrust is zero (this assessment is actually supported by the numerical results presented in Section 5). Note that, in a real system, the driving mechanism would carry a payload that produces drag and the reaction force would compensate that drag. Although the mean thrust coefficient is systematically zero in the present case, two metrics are readily available: the mean velocity and the required power input to sustain a constant mean velocity. By normalizing these quantities, one obtains the normalized mean velocity and the power coefficient given respectively by:

U⁎ =

CP =

vx , fh0 F′P·y^ ḣ 1 ρ (h0f )3c 2 f

(5)

+

MP′θ ̇ , 1 ρf (h0f )3c 2

(6)

^ is where F′P and MP′ are respectively the reaction force and moment per unit span on the driving mechanism at point P and y a unit vector pointing in the y-direction. On the other hand, the definition of efficiency requires more discussion in the context of self-propelled flexible flapping foils since the thrust is zero. Indeed, the traditional definition of efficiency for propulsive flapping wing systems is typically expressed as the thrust power over the power input. For a flapping wing with a varying periodic x-velocity, this efficiency definition is thus given by:

T vx , P

(7)

where T and P are respectively the thrust acting on the driving mechanism and power input per unit span. Obviously, a zero thrust results systematically in a zero efficiency. While one could attempt to isolate thrust and drag from the total aerodynamic force acting along the flight direction, such a separation remains arbitrary in the context where there is no payload in the model. To circumvent this problem, alternate efficiency definitions are proposed in the literature. Among these definitions, one that seems to gain acceptance from the community is the ratio of kinetic energy of the foil over the input energy during one period of oscillation (Kern and Koumoutsakos, 2006; Zhang et al., 2010; Hua et al., 2013; Zhu et al., 2014a, 2014b). Using the parameters defined in Table 1, this definition would scale as η ∼ U⁎/( CPΣ ). While other definitions exist, Zhu et al. (2014a) discuss these other alternatives and highlight the fact that the latter is the only one that has an actual physical interpretation. Yet, this definition has two drawbacks. Firstly, the thrust is not directly responsible for the actual

50

M. Olivier, G. Dumas / Journal of Fluids and Structures 64 (2016) 46–66

kinetic energy of the foil, but rather for a kinetic energy variation in a context where there is no drag from an external source such as a body. Indeed, in most cases, the foil requires more or less than exactly one single cycle to gain its terminal average velocity. In that context, comparing the kinetic energy with the input energy of a single cycle is meaningless. Secondly, when the mass of the foil tends toward zero, the efficiency also tends toward zero (because Σ is then very large). While the concept of a massless foil is purely theoretical, it represents the limiting case of very strong fluid–solid interactions and it is, as such, of great interest. A similar argument can be made about foils with very large masses that would produce high efficiencies. Another alternative would be to simply use the ratio U⁎/ CP . As pointed out by Yeh and Alexeev (2014), this intuitive ratio can be obtained by using Eq. (7) and by assuming that the mean thrust is estimated by T = 21 ρf (h0f )2c . However, this thrust estimation does not take into account the effect of the specific shape of the thruster nor possible Reynolds number effects (that would change the thrust coefficient which is approximated to 1 in this definition). Therefore, while this definition has some physical background, the actual efficiency values remain somewhat arbitrary. The efficiency definition used in this study is inspired on the definition used by Thiria and Godoy-Diana (2010) and Ramananarivo et al. (2011). The idea is to assume that the thrust compensates the drag that would act on the flapping wing if it were moving at its terminal averaged speed without flapping. Therefore, the thrust is given by T = (1/2)CDρf vx 2c in which the drag coefficient can be obtained by a correlation based on independent numerical simulations of the wing in a steady configuration (i.e., a flat plate of finite thickness parallel to the flow in laminar conditions, see Appendix A for the actual correlation). Using Eq. (7), it follows that the efficiency is defined as: 3

η=

CD U⁎ . CP

(8)

In the previous equation, the drag coefficient is defined with the characteristic velocity vx such that the drag is given by D = (1/2)CDρf vx 2c (as opposed to CP, which is based on the characteristic velocity h0f ).

3. Governing equations and numerical methods 3.1. Governing equations Low Reynolds number Newtonian incompressible flows are considered in this study. Because the wing considered moves and deforms, the governing equations are expressed in the integral form over an arbitrarily moving volume, sometimes referred to as the arbitrary Lagrangian Eulerian approach (see Donea et al., 2004). The equations for the fluid flow are:

∂ ∂t

∫V (t) ρf dV + ∫S(t) ρf n^ ·c dS = 0,

∂ ∂t

∫V (t) ρf v dV + ∫S(t) ρf v(n^ ·c) dS = ∫V (t) f dV − ∫S(t) p n^ dS + ∫S(t) n^ ·μ∇v dS,

(9)

(10)

^ is the control surface normal, p is the pressure field, v is the velocity field, v ^ is the where μ = ρf ν is the fluid viscosity, n

Fig. 2. Representation of the beam model. The dashed line represents the beam in the undeformed configuration. The variables dx, dy, p, and ⎛ ∂dy ∂d ⎞ represented at specific locations, but exist everywhere in the beam. The angle ψ can be obtained locally with − ⎜1 + x ⎟ tan(ψ ) = 0 . ∂x ∂x ⎠ ⎝

τ

are

M. Olivier, G. Dumas / Journal of Fluids and Structures 64 (2016) 46–66

51

^ is the relative convective velocity, and f is a generic body force term used to velocity of the control surface, c = v − v introduce a source term that allows computations to be made in an accelerating reference frame. The specific reference frame will be described in more detail later. In this study, the flexible foils are assimilated to a geometrically nonlinear beam model that represents a thin plate section. The beam model is based on the Euler–Bernoulli theory and supports large displacements. As opposed to the fluid equations, the solid equations are based on a Lagrangian formulation reported in the initial configuration. As such, the introduction of a moving reference frame does not need to be taken into account in the solid equations which are (see Epstein and Murray, 1976):

ρs e

ρs e

∂ 2dx ∂t 2 ∂ 2dy ∂t

2

=

=

2 ⎤ ⎡ ⎛ ∂dy ∂dx ⎞⎤ ∂dx ⎞ ∂ ⎡ ⎛ ∂ ⎢ ∂ dy ⎥ ∂ 2 ⎡ ∂dy ⎤ ⎥−p + τ⎜ 1 + ⎢ N′⎜ 1 + M′ 2 + 2 ⎢ M′ ⎟, ⎟⎥ + ⎝ ∂x ⎣ ⎝ ∂x ⎠⎦ ∂x ⎢⎣ ∂x ⎦ ∂x ∂x ⎠ ∂x ⎥⎦ ∂x ⎣

(11)

⎛ ∂dy ∂dx ⎞⎤ ∂dx ⎞ ∂ ⎡ ∂dy ⎤ ∂2 ⎡ ⎛ ∂ ⎡ ∂ 2dx ⎤ ⎢ M′ 2 ⎥ − p⎜ 1 + ⎢ N′ ⎥ − 2 ⎢ M′⎜ 1 + , ⎟⎥ − ⎟−τ ⎝ ∂x ⎣ ∂x ⎦ ∂x ⎣ ⎝ ∂x ⎠⎦ ∂x ⎣ ∂x ⎠ ∂x ∂x ⎦

(12)

where dx and dy are respectively the x- and y-components of the displacement field, p is the pressure (taken from the flow) on the deformed configuration, and τ is the shear (also taken from the flow) on the deformed configuration (see Fig. 2). The internal normal force N′ as well as the internal bending moment M′ is given by the Euler–Bernoulli beam theory:

⎡ ⎛ 2 ⎛ ∂dy ⎞2⎞⎤ ∂d 1 ⎛ ∂d ⎞ ⎟ ⎟⎥, N′ = Eee˜ = Ee⎢⎢ x + ⎜ ⎜ x ⎟ + ⎜ ∂x 2 ⎜⎝ ⎝ ∂x ⎠ ∂x ⎠ ⎟⎠⎥ ⎝ ⎦ ⎣

(13)

⎤ ⎡ ∂ 2d ⎛ ∂dx ⎞ ∂dy ∂ 2dx ⎥ y M′ = EI′χ˜ = EI′⎢ 2 ⎜ 1 + , ⎟− 2 ∂x ⎠ ∂x ∂x ⎥⎦ ⎢⎣ ∂x ⎝

(14)

where e˜ is the Green–Lagrange strain and χ˜ is the curvature of the deformed beam reported in the undeformed configuration. While Eqs. (11)–(14) allow the beam to deform axially, axial deformations are not significant in this study since thin beams (e/c = 0.01 here) are much more sensitive to transverse loads. At the fluid–solid interface, momentum and mass conservation yields the following conditions:

vf = vs,

(15)

^ ·σ = n ^ ·σ , n f s

(16)

where σ is the Cauchy stress field and the subscripts f and s represent respectively the fluid and the solid.

3.2. Numerical methods The numerical code used to solve the flow equations is based on the OpenFOAM library (OpenFOAM Foundation, 2013) and, as such, it uses a second order finite-volume discretization. The temporal discretization is handled by a second-order backward scheme and the pressure–velocity coupling is achieved by an iterative PISO algorithm (also called PIMPLE in the OpenFOAM documentation). The discretization takes into account the geometric conservation law to ensure proper conservation of mass (Demirdžić and Perić, 1988; Donea et al., 2004). The mesh motion is handled by an inverse distance weighting interpolation technique based on Witteveen (2010). The interpolation is made by using only the nodes on the wing to compute the mesh motion and radial damping is used to damp the motion in the far field. Regarding the solid, the discretization is done with the Galerkin finite-element method along with the Hermite polynomials shape functions. The Newton–Raphson method is used to handle nonlinear terms. The fluid–solid coupling is made by using radial basis function interpolation to transfer the loads and the displacements at the fluid–solid interface (Smith et al., 1995; Beckert and Wendland, 2001; de Boer et al., 2007). An artificial compressibility term is introduced in the pressure equation in order to stabilize the partitioned coupling scheme (Järvinen et al. (2008); Degroote et al. (2010)). More details on the numerical methods are provided in Olivier (2014) as well as in Olivier and Dumas (2016). The problem of the self-propelled wing can be set in the numerical code by applying the following boundary conditions on the leading edge in the beam solver:

52

M. Olivier, G. Dumas / Journal of Fluids and Structures 64 (2016) 46–66

⎛ ∂ 2dy ∂dx ⎞ ∂ ⎛ ∂dy ⎞ ⎜ M′ ⎟ = Fx′,LE, N′⎜ 1 + ⎟ + M′ 2 + ⎝ ∂x ⎠ ∂x ⎝ ∂x ⎠ ∂x d y = h, ⎛ ∂dx ⎞ − ⎜1 + ⎟ tan(θ ) = 0, ⎝ ∂x ∂x ⎠

∂dy

(17)

′ is the x-component of the aerodynamic force on the leading edge (the rounded part) which is applied as a nodal where Fx,LE force while h and θ are respectively given by Eqs. (1) and (2). These conditions correspond respectively to a natural condition acting mainly in the x-direction and to two essential conditions (on dy and ∂dy /∂x ) acting mainly in the y-direction. However, because of the geometrically nonlinear nature of the beam equations used here, the first and third conditions in Eq. (17) actually act on both displacement components. On the trailing edge, the aerodynamic force is applied as a nodal force such that the following natural conditions apply:

⎛ ∂ 2dy ∂dx ⎞ ∂ ⎛ ∂dy ⎞ ⎜ M′ ⎟ = Fx′,TE, N′⎜ 1 + ⎟ + M′ 2 + ⎝ ∂x ⎠ ∂x ⎝ ∂x ⎠ ∂x N′

∂dy ∂x

− M′

∂ 2dx ∂x

2

−

∂dx ⎞⎞ ∂ ⎛ ⎛ ⎜ M′⎜ 1 + ⎟⎟ = F′y,TE, ∂x ⎝ ⎝ ∂x ⎠⎠

′ , M′(2e˜ + 1) = MTE

(18)

′ is the ′ and F′y,TE are respectively the x- and y-component of the aerodynamic force on the trailing edge, and MTE where Fx,TE aerodynamic moment on the trailing edge. These boundary conditions on the beam model are similar to those used by Zhu et al. (2014b) except that here, the pitching amplitude is not zero and the fluid load on the leading and trailing edges is taken into account since the thickness of the plate is finite in the fluid domain. At the inlet of the fluid domain, a Dirichlet condition specifying an arbitrary constant velocity is applied along with a zero normal gradient for the pressure. At the outlet, a zero normal gradient is imposed on the velocity field and the pressure is set to 0. Symmetry conditions are applied at the top and bottom boundaries. Using this numerical setup, the flapping wing effectively moves freely. However, after a few oscillation cycles, the wing gets closer to the domain boundary resulting in the following consequences: the inlet uniform boundary condition is no longer physical as the wing is too close and the mesh deteriorates. Setting the inlet velocity equal to the wing average x-velocity would solves the problem of reaching the mesh boundary, but the latter is not known prior to run the simulation. Alternatively, the numerical code can be adapted to handle the problem in a moving frame of reference. The moving frame is chosen to move only in the x-direction, following the displacement of the wing leading edge position at the previous time-step. The reason why the frame position is chosen to lag one time-step is to avoid convergence deterioration. Indeed, the position of the wing is a part of the solution and it is not available at the beginning of a time-step. Using the value of the previous time-step thus avoids the need to use iteratively changing source terms and boundary conditions. It follows that the x-displacement of the wing leading edge with respect to its initial position in the reference frame is only the difference in displacement from two successive time-steps, which is always very small. The velocity of the moving frame is used as an inlet boundary condition and its acceleration is added as a source term in the momentum equation of the fluid domain. Thus, at time index n + 1, the specific boundary condition is vinlet = − vxn,P and the source term is set as f = ( − axn,P , 0, 0) where vxn,P and axn,P are respectively the velocity and the acceleration of point P on the wing taken at the previous time-step. The beam solver does not need to be modified as the computations are performed in the initial (undeformed) configuration of the beam. However, when it comes to computing fluid–solid interface conditions, structural grid points need to be temporarily moved to match the deformed wing position in the fluid domain reference frame in order to perform the interpolation correctly. The wing displacement field interpolated and applied on the mesh motion algorithm is the actual wing displacement minus the displacement of point P taken at the previous time-step, dxn,P . The boundary condition on the wing surface does not need to be modified since it is automatically computed using the fluid mesh node displacements. The fluid load on the structure is interpolated in the usual way. For the sake of simplifying the imposition of the outlet condition, the mesh external boundary is a square. That way, a constant uniform pressure can still be imposed on the outlet along with a zero normal gradient for the velocity. Otherwise, the hydrostatic pressure gradient produced by the source term would have made the pressure non-uniform on a curved outlet. Symmetry planes are used on the top and bottom boundaries and the moving frame velocity is imposed on the left side of the domain along with a zero normal gradient for the pressure.

4. Verification and validation The numerical code has been thoroughly validated and some validation tests involving a similar numerical setup as well as comparisons with experiments by Heathcote and Gursul (2005, 2007) have been presented in a previous paper (Olivier and Dumas, 2016). Nevertheless, the numerical resolution used for this specific problem and the strategy of employing a

M. Olivier, G. Dumas / Journal of Fluids and Structures 64 (2016) 46–66

53

Fig. 3. Mesh used in the fluid domain along with the far-field boundary conditions. The mesh contains 62,801 cells (medium resolution).

source term and time-varying boundary conditions based on the wing frame of reference needs to be verified. For that purpose, a specific flexible flapping wing scenario is used on different grid resolutions and time-step sizes, and the corresponding case is also run in the inertial reference frame for a short period of time (up to the point where the wing gets close to the inlet boundary). The parameters used for this assessment are:

Ref = 200,

h⁎ = 1,

θ0 = 30°,

δP⁎ = 0.324,

Σ = 5.

(19)

Numerical results are obtained using 3 different meshes and time-step sizes. The coarse mesh contains 15,376 cells with f Δt = 5 × 10−4 , the medium mesh contains 62,801 cells with f Δt = 2.5 × 10−4 , and the fine mesh contains 235,998 cells with f Δt = 1.25 × 10−4 . All three meshed domains have the same geometry which consists of a large square of length 80 c centered around the flat plate (see Fig. 3). The solid domain contains respectively 25, 50, and 100 beam elements uniformly distributed. The effects of the grid resolution and time-step size on the velocity of the wing and on the instantaneous vorticity field are illustrated respectively in Figs. 4 and 5. The results show that using the medium mesh is sufficient for the sake of the proposed study. Moreover, the solution obtained in the inertial reference frame confirms that the moving reference frame implementation is correct.

Fig. 4. Comparison of the instantaneous velocity of the wing obtained with different mesh resolutions and time-steps. The result obtained with the inertial reference frame is also shown.

54

M. Olivier, G. Dumas / Journal of Fluids and Structures 64 (2016) 46–66

Fig. 5. Vorticity field at ft ¼1.2 with the (a) coarse, (b) medium, and (c) fine resolution (mesh and time-step). The scenario run in the inertial reference frame (medium resolution) is shown in (d).

5. Results and discussion Numerical simulations have been carried out with the following dimensionless parameters:

Ref = 200,

h⁎ = 1,

θ0 = 30°.

(20)

The medium mesh and time-step have been used for these simulations which are run up to a physical time where the cycleaveraged metrics are periodic (constant from one cycle to another). 5.1. Deviation motion Numerical results show that, in all cases presented, the wing reaches a periodic motion in which the time-averaged xvelocity stabilizes to a constant value. Nevertheless, the wing still oscillates in the x-direction. In fact, this motion in the xdirection distinguishes the self-propelled wing from oscillating wings with an imposed constant x-velocity. This oscillation transverse to the stroke plane is called deviation. It is worth pointing out that, in the context of the self-propelled wing model used in this study where the wing is not attached to a body, the deviation motion is inherently passive. To illustrate the deviation motion, the harmonic motion of Eq. (1) is used along with a harmonic deviation given by:

dx,P =

h0 cos(4πft + ϕ), 4

(21)

where ϕ is an arbitrary phase shift. Note that the amplitude h0 /4 has been chosen arbitrarily for illustration purposes, but the fact that the deviation oscillates at twice the frequency of the heaving motion is inspired from actual insect wing kinematics (Fry et al., 2005). The resulting motion shown in Fig. 6a is similar to the harmonic model proposed by Bos et al.

Fig. 6. Deviation motion of rigid wings. (a) Harmonic heaving motion (Eq. (1)) combined with an imposed harmonic deviation (Eq. 21)). (b) Actual motion of rigid wings for different values of Σ obtained from numerical simulations.

M. Olivier, G. Dumas / Journal of Fluids and Structures 64 (2016) 46–66

55

(2008). It is shown that for ϕ = 0, the upstroke and the downstroke follow the same path and that for nonzero phase shifts, a figure-of-eight appears. The openness of the figure-of-eight reaches a maximum at ϕ = ± π /2. Positive phase angles (not shown) produce the same motion as their corresponding negative phase angles, but the motion direction is reversed. While Bos et al. (2008) propose that the deviation motion is used by insects to level the forces throughout an oscillation cycle, another interpretation is presented here. Indeed, numerical results suggest that the specific kinematics related to the deviation motion may be a consequence of passive mechanisms caused by the fluid–structure dynamics rather than a deliberate action. Three fluid–structure passive mechanisms can be identified from the present results: two aerodynamic feathering mechanisms identified as rigid feathering and flexible feathering as well as a structural mechanism called the inertial displacement which is caused by the wing deformation. 5.1.1. Rigid feathering Rigid feathering is an aerodynamic feathering mechanism that acts on self-propelled flapping wings. This feathering effect is caused by the wing horizontal acceleration resulting from aerodynamic forces which causes the wing to move back and forth. As results show, this effect is dominant on rigid wings, hence the name rigid feathering. Rigid feathering occurs when the wing reaches a configuration (i.e., velocity and angle of attack) that would produce a thrust increase. Since the wing is let free to move horizontally, it accelerates which has the immediate effect of reducing the effective angle of attack (see Eq. (3)) and, subsequently, the thrust itself. The high-thrust features such as moderate leading edge vortices are also much less intense when substantial rigid feathering occurs. Since the acceleration is inversely proportional to the mass for a given force, the rigid feathering effect is more important for light wings (high Σ). Trajectories obtained from numerical simulations of rigid wings are presented in Fig. 6b. As expected, the heavier wing (Σ = 0.1) has a nearly vertical trajectory with a very slight figure-of-eight openness. This means that the deviation is almost in phase with the heaving motion and that the magnitude of the deviation is small. Therefore, the kinematics of freely moving heavier wings is effectively very similar to that of fully constrained wings with an imposed x-velocity. This statement can be assessed further by observing the strong vortex shedding in Fig. 7 and the corresponding pressure fields in Fig. 8 which are similar to what is observed when the wing x-velocity is imposed for the same Σ value (see Olivier and Dumas, 2016). When Σ is high (pressure-driven deformation cases), the wing tends to follow a path that minimizes the aerodynamic force in the x-direction. In that sense, this scenario is analogous to a free-falling piece of paper that actually follows a path where the aerodynamic force tends to compensate the gravity (see Wang, 2005). Fig. 9 presents the instantaneous aerodynamic thrust coefficient for cases involving rigid wings and different Σ. It is found that the rigid feathering effect actually reduces the magnitude of the thrust. Moreover, the rigid feathering effect also has the consequence of leveling the thrust force. This is in good agreement with an analogous finding related to the deviation motion reported by Bos et al. (2008). 5.1.2. Flexible feathering The impact of wing flexibility on the figure-of-eight pattern of high Σ cases is shown in Fig. 10a. It appears that increasing wing flexibility has two consequences: it reduces the magnitude of the deviation motion and it modifies the openness of the figure-of-eight and eventually reverses it. These consequences can be explained by another feathering mechanism that affects flexible wings and which will therefore be referred to as flexible feathering. Flexible feathering occurs when a flexible wing deforms in a way that realigns its profile with the local flow, which reduces the effective angle of attack (see Zhu, 2007). As it was discussed earlier, the deviation motion of light wings is caused by rigid feathering which results from the fact that the wing accelerates along the x-axis in a way that reduces its angle of attack (see Eq. (3)). However, when Σ is high, flexible wings also undergo flexible feathering which has the same consequence on the effective angle of attack. The reduction of the effective angle of attack is then caused by an increased effective pitching angle θEff (see Eq. (4)). Since both feathering effects cannot alter the angle of attack beyond the feathering limit, if flexible feathering occurs, rigid feathering, and thus the magnitude of the deviation motion, will unavoidably be less important. The second consequence on the deviation motion, that is the reversal of the figure-of-eight, can be explained by the effect of flexible feathering on the evolution of the thrust coefficient which is reported in Fig. 11. Effectively, the flexible feathering effect has the global effect of reducing the aerodynamic thrust coefficient. Moreover, it appears that rigid cases experience a drag peak near ft¼0 and ft¼ 0.5. This drag peak results in a significant backward acceleration and, subsequently, a significant backward motion when the wing is light (Σ = 50 for instance). The flexible feathering effect that occurs when the wing flexibility is increased causes this drag peak to almost vanish, which has the effect of reversing the figure-ofeight path. 5.1.3. Inertial displacement Decreasing Σ (e.g., by increasing the wing mass), has the effect of increasing the relative inertia of the wing. This has the consequence that aerodynamic forces have much less immediate impact on the kinematics of the wing. That is, it takes longer for the wing to reach its terminal average velocity and the velocity fluctuations in the x-direction caused by the inherent unsteadiness of the flow around the flapping wing are much less important. This means that the deviation motion is reduced as depicted in Fig. 10b. Yet, increasing the flexibility in low Σ cases produces an increase in the deviation amplitude. This is explained by the geometrically nonlinear deformation of the wing. Indeed, while the inertia of the wing causes the deformation by pulling the wing upward (or downward) at the extremities of its course, the fact that the

56

M. Olivier, G. Dumas / Journal of Fluids and Structures 64 (2016) 46–66

Fig. 7. Vorticity field for various flexibilities involving weak interaction (Σ = 0.1).

deformation is large causes the trailing edge of the wing to go slightly forward. This has the consequence of pulling the leading edge backward so that the center of mass position remains mostly unchanged (except for the slight acceleration caused by unsteady aerodynamic forces).

M. Olivier, G. Dumas / Journal of Fluids and Structures 64 (2016) 46–66

Fig. 8. Pressure coefficient field for various flexibilities involving weak interaction (Σ = 0.1) .

57

58

M. Olivier, G. Dumas / Journal of Fluids and Structures 64 (2016) 46–66

Fig. 9. Aerodynamic thrust coefficient of rigid wings for various

Σ.

Fig. 10. Wing displacement relative to a point that moves at the wing average velocity for various flexibilities δP⁎ with (a) Σ = 50 and (b) Σ = 0.1.

5.1.4. Deviation motion of insect wings Using insect wing characteristics provided by Combes and Daniel (2003a) and references herein (wing density of 1200 kg/m3, wing thickness of about 45 μm ), it is possible to estimate the value of Σ for insect wings which is about 0.2. This estimate is obtained by using a heaving amplitude of 0.01 m which corresponds to the average chord length of the species reported by Combes and Daniel (2003a). This estimate confirms the trend reported in the literature that insect wings deformation is mainly governed by inertia forces (Daniel and Combes, 2002; Combes and Daniel, 2003b; Thiria and GodoyDiana, 2010). However, considering the uncertainties and the fact that the actual heaving amplitude varies along the wing (because the wing undergoes a rotational stroke motion as it is held at its root), it is not excluded that the actual value of Σ of insect wings becomes sufficiently high to increase the effects of pressure forces over inertia forces and eventually to produce some feathering effects. This is in good agreement with the findings of Ramananarivo et al. (2011), who point out that insect wing deformation is dominated by inertia effects, but nonlinear damping produces significant effects on the phase lag of the deformation. Improvements in performance are obtained for large phase lags. As such, even if the wing deformation is mainly governed by its inertia, the fluid damping can have an important impact on the performances.

M. Olivier, G. Dumas / Journal of Fluids and Structures 64 (2016) 46–66

59

Fig. 11. Aerodynamic thrust coefficient of flexible wings for Σ = 50 .

The kinematics of the fruit fly wings is reported by Fry et al. (2003) and the corresponding 2D models are proposed by Bos et al. (2008). Moreover, Liu and Sun (2008) measured the wing kinematics of drone flies. In both species, the wing trajectory can be approximated fairly well by harmonic functions so that the figure-of-eight pattern presented in Fig. 6a is similar to the wing trajectory of an actual insect when ϕ = − π /20. Moreover, while the motion proposed here has a vertical stroke plane as opposed to a horizontal stroke planes for hovering insects, it is worth pointing out that some insects actually control their velocity by inclining the stroke plane at various angles (Willmott and Ellington, 1997; Azuma, 2006). As such, the wing motion used in this study is representative of that of insect wings. By comparing this proposed harmonic trajectory using ϕ = − π /20 with Fig. 6b, it seems that the wings of fruit flies and drone flies follow a path that is similar to wings undergoing a passive rigid feathering mechanism with low-to-moderate Σ values. Indeed, if the deviation was caused by inertial displacement, the wing would follow its path in the reversed direction. Therefore, it seems likely that some insects actually rely on passive feathering to level the forces during wing strokes. 5.2. Performances of freely moving wings Before discussing performances of self-propelled flexible flapping wings, it is worthwhile to look at the results of rigid wings. Indeed, even if parameters of Eq. (20) are held fixed, the flow field depends on Σ even for rigid wings (δP⁎ = 0). This is not the case for rigid wings when the x-velocity is imposed. Fig. 12 reports the effect of varying Σ on the performance metrics. It is found that both the optimal mean normalized velocity and the optimal efficiency are obtained when Σ is about 10. The power coefficient also reaches a local maximum at that point and a local minimum is observed near Σ = 0.6. Interestingly, the latter point is very close to the estimation of Σ of insect wings discussed earlier. The local power coefficient maximum can be explained by the larger amplitude of the deviation motion when Σ is high. Indeed, these cases exhibit a larger x-velocity at the ends of the vertical course. This results in a larger effort required to pitch the wing. Increasing Σ beyond 10 has a similar effect, but the wing reacts more rapidly to the forces coming from the flow so that the power peak required to rotate the wing is still significant but has a shorter duration. Also, it is worth mentioning that the fact that some curves cross each other in Figs. 13a and 14 shows that optimal performances are not necessarily always obtained when

Fig. 12. Propulsion performance metrics for rigid freely moving wings with respect to (c) efficiency.

Σ:

(a) mean normalized velocity; (b) mean power coefficient;

60

M. Olivier, G. Dumas / Journal of Fluids and Structures 64 (2016) 46–66

Fig. 13. Efficiency of the freely moving flexible flapping wing as a function of (a) the flexibility and (b) the frequency ratio. Various pressure-to-inertia ratios are shown.

Σ = 10 if the flexibility is increased. This suggests that the optimal value of Σ found here does not necessarily represent a general result since it is likely to change with the problem configuration (i.e., parameters given by Eq. (20)). This remains to be investigated. The mean normalized velocity can be related to the reduced frequency by f ⁎ = 1/( U⁎h⁎) while the Strouhal number can be approximated by St ≃ 2/ U⁎ if one assumes the wake width is given by twice the heaving amplitude h0. By considering results of Fig. 12a, it appears that the reduced frequencies of the rigid cases studied span the range 0.092–0.12, which corresponds to a Strouhal number in the range 0.184–0.24. This is in good agreement with the Strouhal number range reported in flying and swimming animals which is between 0.2 and 0.4 (see Triantafyllou et al., 2000). Moreover, rigid wing results from Olivier and Dumas (2016) report that the optimal efficiency for cases with an imposed freestream velocity is obtained when the reduced frequency is near 0.15 for θ0 = 30°. Before discussing this matter further, it is worth recalling that numerical simulations with an imposed x-velocity (or experiments in wind or water tunnels with a fully constrained flapping wing) actually represent a scenario of a self-propelled flapping wing carrying a payload that would have the following characteristics:


The interaction of the payload with the flow surrounding the flapping wing is minimal.
The mass of the payload is significantly higher than the mass of the wing. This would have the effect of canceling the back and forth motion (deviation) of the wing.


The drag of the payload would correspond to the thrust of the constrained flapping wing at the given flapping configuration. While it is not strictly possible to mimic a fully constrained configuration by increasing Σ (because the resulting thrust of the self-propelled wing will remain zero), the trend of Fig. 12a suggests that reducing Σ further may increase the reduced frequency further to a value that may be close to 0.15. Therefore, the present results support the proposition of von Ellenrieder et al. (2008) that optimal oscillating wings in nature reach an optimal operation by means of a locking phenomenon. Performance metrics of self-propelled flexible wings are reported in Figs. 13–15. The trends found in these figures are analogous to the ones reported by Olivier and Dumas (2016) when the x-velocity is imposed. Indeed, increasing the wing flexibility reduces the efficiency for inertia-driven deformation cases (low Σ) while a moderate amount of flexibility increases the efficiency of pressure-driven deformation cases (high Σ). However, the increase in efficiency in this case is slightly less significant than it was with cases where the x-motion of the wing was imposed. This can be readily explained by the additional rigid feathering experienced by the freely moving wing. Indeed, Fig. 16 shows that the additional aerodynamic feathering effect suppresses the leading edge vortex for all flexibilities. However, high-efficiency cases reported by Olivier and Dumas (2016) always benefit from moderate leading edge vortices, and thus delayed stall, that usually increase aerodynamic forces without penalizing the power input. Increasing the wing flexibility also reduces the mean normalized velocity of all cases because of the so-called flexible feathering effect. This can be directly related to the diminution of the pressure field variation around the wing depicted in Fig. 17. This behavior is analogous to the effect of flexibility on the thrust coefficient in Olivier and Dumas (2016). This trend

M. Olivier, G. Dumas / Journal of Fluids and Structures 64 (2016) 46–66

61

Fig. 14. Mean normalized velocity of the freely moving flexible flapping wing as a function of the flexibility. Various pressure-to-inertia ratios are shown.

Fig. 15. Mean power coefficient of the freely moving flexible flapping wing as a function of the flexibility. Various pressure-to-inertia ratios are shown.

of the terminal velocity may seem in contradiction with what is often reported in the literature, i.e., that moderate flexibility usually implies a significant increase of the terminal averaged velocity or thrust in the case of fully constrained flapping wings (Heathcote and Gursul, 2007; Hua et al., 2013; Yeh and Alexeev, 2014; Zhu et al., 2014b). However, in these studies involving flexible flapping plates, only a heaving motion is imposed. In this context, flexibility helps in re-orienting the thrust force in the x-direction, which is highly beneficial, as shown by Olivier and Dumas (2016). In this study, this effect of force re-orientation is mitigated by the imposed pitching motion, which explains why flexibility does not increase the terminal velocity. Lastly, the trend of the power coefficient is also very similar to what was found by Olivier and Dumas (2016) for fully constrained wings. Indeed, wing flexibility reduces the power consumption because aerodynamic forces that act against the imposed motion are reduced by flexible feathering when Σ is high. On the contrary, when Σ is low, the power consumption increases with flexibility because the wing deformation causes strong boundary layer separations (see Fig. 7) that result in higher aerodynamic forces that are in opposition with the prescribed motion.

62

M. Olivier, G. Dumas / Journal of Fluids and Structures 64 (2016) 46–66

Fig. 16. Vorticity field for various flexibilities involving strong interaction (Σ = 50) .

M. Olivier, G. Dumas / Journal of Fluids and Structures 64 (2016) 46–66

Fig. 17. Pressure coefficient field for various flexibilities involving strong interaction (Σ = 50) .

63

64

M. Olivier, G. Dumas / Journal of Fluids and Structures 64 (2016) 46–66

0.35

Simulations Polynomial fit Blasius solution

0.3 0.25

CD

0.2 0.15 0.1 0.05 0

0

1000

2000

Re

3000

4000

5000

Fig. 18. Drag coefficient of the non-flapping wing against the Reynolds number obtained from numerical simulations. The polynomial fit and the Blasius solution are also shown.

In order to assess the effects of resonance, Fig. 13b introduces the frequency ratio which is the frequency of the imposed motion over the natural frequency of the first vibration mode of a beam (see Young and Felgar, 1949). The frequency ratio is given by:

⎛ δ⁎ ⎞ f ≃ 1.787 ⎜ P ⁎ ⎟ ⎝ Σh ⎠ f1

(22)

It is shown that most cases have a frequency ratio below 1, which means that structural resonance does not occur, at least when Σ is low. Indeed, when Σ is high, the fluid added-mass effect lowers the natural frequency of the coupled fluid–solid system which could enable resonance. However, the actual wing deformation observed in most cases does not exhibit higher modes unless the wing is very flexible, which corresponds to inefficient configurations. A similar investigation on fully constrained flapping wings showed, however, that for similar configuration at high Σ values, the second and third deformation modes have local effects on the efficiency (Olivier and Dumas, 2016). Specifically, the appearance of the second and third modes produces respectively a slight efficiency drop and increase. Moreover, it appears that the impact of resonance is more obvious when θ0 = 0. Note that for similar parameters, the second mode occurs at f /f1 ≃ 0.9 for Σ = 10 and at f /f1 ≃ 0.5 for Σ = 50. In the present study, it appears that letting the wing free to move in the x-direction, along with the fact that θ0 = 30°, inhibits this resonance effect. This can be assessed in Fig. 13 where no obvious local variations are observed near these frequency ratios.

6. Conclusion The self-propelled 2D flexible flapping wing problem is presented with an unconstrained motion in the x-direction while a combined pitching and heaving motion is imposed. That is, the wing accelerates itself until it reaches a constant cycleaveraged velocity. The simulation code performs the computations in a reference frame that moves with the wing in the xdirection. This variant of the problem has the same dimensionless parameters as the problem of the driven flapping wing except that the Reynolds number is now defined with the heaving velocity h0f and that the reduced frequency is no more a parameter, but an actual result. Regarding the performance metrics, the situation is different from the case of the fully constrained flexible flapping wing. Indeed, for this specific problem, the mean thrust coefficient is zero when the periodic regime is reached. Thus, the metrics used to characterize the performance of the system are: the cycle-averaged normalized velocity, the mean power coefficient, and an efficiency definition similar to that of Thiria and Godoy-Diana (2010) and Ramananarivo et al. (2011). Simulation results of self-propelled flapping wings exhibit periodic deviation motion. Comparing the overall wing motion with that of insects reported in the literature, it can be inferred that the deviation motion of insect wings benefits from a passive reaction that results from unsteady forces. In fact, the deviation motion appears to be the result of what was defined as a rigid feathering effect. This also implies that aerodynamic forces cannot completely be ignored when evaluating the deformation of insect wings during flight. Moreover, the effect of wing inertia and flexibility on the deviation motion is presented. Increasing wing inertia results in a reduction in the amplitude of the deviation motion. Increasing the flexibility in the case of pressure-driven deformation configurations also results in a reduction of the amplitude of the deviation and eventually has the effect of inverting the figure-of-eight pattern because of additional aerodynamic feathering caused by flexibility. In the case of high-inertia configurations, increasing the flexibility increases the amplitude of the deviation

M. Olivier, G. Dumas / Journal of Fluids and Structures 64 (2016) 46–66

65

because the location of the center of mass changes with respect to the leading edge position. This effect is referred to as inertial displacement. Regarding the influence of wing inertia and flexibility on performance metrics, it is found that the main trends are similar to the case of the fully constrained flexible flapping wing that undergoes the same pitching and heaving motion. That is, there is an optimal flexibility that maximizes the efficiency when Σ is high. On the other hand, when Σ is low, flexibility deteriorates the efficiency systematically. Increasing the flexibility also decreases the cycle-averaged velocity for any Σ value. While the physical mechanisms at stakes are similar to that of the constrained wing, an additional feathering effect, namely the rigid feathering, occurs when Σ is high. Acknowledgments The authors would like to acknowledge financial support from the Natural Sciences and Engineering Research Council (NSERC) of Canada and from the Fonds de Recherche du Québec – Nature et Technologies (FRQNT) as well as the access to high-performance computer facilities of Compute Canada.

Appendix A. Thrust of the self-propelled flapping wing The thrust of the self-propelled flapping wing is determined by its drag when the latter is not flapping. The drag is obtained by running numerical simulations of the steady wing at different constant terminal velocities vx (see Fig. 18). A correlation providing the drag coefficient CD = 2D/ρf vx 2c against the Reynolds number Re = vxc/ν is obtained with a polynomial fit such that:

CD =

102.57 1.0130 + + 0.013861. Re Re1/2

(A.1) 2

the determination coefficient of this regression is r = 0.99998. it is worth noting that, for thinner plates at a slightly higher Reynolds number, the Blasius solution cd = 2 × 1.328 Re−1/2 could be used. indeed, at Re¼4800, the viscous contribution of the drag coefficient (i.e., by removing the pressure contribution on the leading and trailing edges) is 0.040, which is in very good agreement with the Blasius solution.

References Alben, S., Witt, C., Baker, T.V., Anderson, E., Lauder, G.V., 2012. Dynamics of freely swimming flexible foils. Phys. Fluids 24 (5). Azuma, A., 2006. The Biokinetics of Flying and Swimming, 2nd edition. AIAA, Reston, Virginia. Beal, D., Hover, F., Triantafyllou, M., Liao, J., Lauder, G., 2006. Passive propulsion in vortex wakes. J. Fluid Mech. 549, 385–402. Beckert, A., Wendland, H., 2001. Multivariate interpolation for fluid–structure interaction problems using radial basis functions. Aerosp. Sci. Technol. 5 (2), 125–134. Bos, F., Lentink, D., van Oudheusden, B., Bijl, H., 2008. Influence of wing kinematic on aerodynamic performance in hovering insect flight. J. Fluid Mech. 594, 341–368. Combes, S., Daniel, T., 2003a. Flexural stiffness in insect wings. I. Scaling and the influence of wing venation. J. Exp. Biol. 206 (17), 2979–2987. Combes, S., Daniel, T., 2003b. Into thin air: contributions of aerodynamic and inertial-elastic forces to wing bending in the hawkmoth Manduca Sexta. J. Exp. Biol. 206, 2999–3006. Daniel, T., Combes, S., 2002. Flexible wings and fins: bending by inertial or fluid-dynamic forces?. Integr. Comparat. Biol. 42, 1044–1049. de Boer, A., Van Zuijlen, A., Bijl, H., 2007. Review of coupling methods for non-matching meshes. Comput. Methods Appl. Mech. Eng. 196, 1515–1525. Degroote, J., Swillens, A., Bruggement, P., Haelterman, R., Segers, P., Vierendeels, J., 2010. Simulation of fluid–structure interaction with the interface artificial compressibility method. Int. J. Numer. Methods Biomed. Eng. 26, 276–289. Demirdžić, I., Perić, M., 1988. Space conservation law in finite volume calculations of fluid flow. Int. J. Numer. Methods Fluids 8, 1037–1050. Donea, J., Huerta, A., Ponthot, J.-P., Rodríguez-Ferran, A., 2004. Arbitrary Lagrangian–Eulerian Methods. Encyclopedia Comp. Mec. 1, 14. Eldredge, J., Pisani, D., 2008. Passive locomotion of a simple articulated fish-like system in the wake of an obstacle. J. Fluid Mech. 607, 279–288. Epstein, M., Murray, D., 1976. Large deformation in-plane analysis of elastic beams. Comput. Struct. 6, 1–9. Fry, S.N., Sayaman, R., Dickinson, M.H., 2003. The aerodynamics of free-flight maneuvers in drosophila. Science 300 (5618), 495–498. Fry, S.N., Sayaman, R., Dickinson, M.H., 2005. The aerodynamics of hovering flight in drosophila. J. Exp. Biol. 208 (12), 2303–2318. Heathcote, S., Gursul, I., 2005. Flexible flapping airfoil propulsion at low Reynolds numbers. In: 43th AIAA Aerospace Sciences Conference and Exhibit, paper AIAA-2005-1405, Reno, Nevada (January). Heathcote, S., Gursul, I., 2007. Flexible flapping airfoil propulsion at low Reynolds numbers. AIAA J. 45 (5), 1066–1079. Hua, R.-N., Zhu, L., Lu, X.-Y., 2013. Locomotion of a flapping flexible plate. Phys. Fluids 25 (12). Järvinen, E., Råback, P., Lyly, M., Salenius, J., 2008. A method for partitioned fluid–structure interaction computation of flow in arteries. Med. Eng. Phys. 30, 917–923. Kern, S., Koumoutsakos, P., 2006. Simulations of optimized anguilliform swimming. J. Exp. Biol. 209 (24), 4841–4857. Lee, J., Lee, S., 2013. Fluid–structure interaction for the propulsive velocity of a flapping flexible plate at low Reynolds number. Comput. Fluids 71 (0), 348–374. Liu, Y., Sun, M., 2008. Wing kinematics measurement and aerodynamics of hovering drone flies. J. Exp. Biol. 211 (13), 2014–2025. Olivier, M., 2014. A Fluid–Structure Interaction Partitioned Algorithm Applied to Flexible Flapping Wing Propulsion (Ph.D. thesis). Université Laval. URL 〈http://www.theses.ulaval.ca/2014/30609/30609.pdf〉. Olivier, M., Dumas, G., 2016. A parametric investigation of the propulsion of 2D chordwise-flexible flapping wings at low Reynolds number using numerical simulations. J. Fluids Struct. 63, 210–237, http://dx.doi.org/10.1016/j.jfluidstructs.2016.03.010. OpenFOAM Foundation, 2013. OpenFOAM – The Open Source CFD Toolbox – User Guide. URL 〈http://www.openfoam.org/〉. Ramananarivo, S., Godoy-Diana, R., Thiria, B., 2011. Rather than resonance, flapping wing flyers may play on aerodynamics to improve performance. Proc. Natl. Acad. Sci. 108 (15), 5964–5969.

66

M. Olivier, G. Dumas / Journal of Fluids and Structures 64 (2016) 46–66

Smith, J., Hodges, D., Cesnik, C., 1995. An evaluation of computational algorithms to interface between CFD and CSD methodologies. Technical Report WLTR-96-3055, Wright Laboratory (November). Spagnolie, S.E., Moret, L., Shelley, M.J., Zhang, J., 2010. Surprising behaviors in flapping locomotion with passive pitching. Phys. Fluids 22 (4). Thiria, B., Godoy-Diana, R., 2010. How wing compliance drives the efficiency of self-propelled flapping flyers. Phys. Rev. E 82 (July), 015303. Triantafyllou, M.S., Triantafyllou, G.S., Yue, D.K.P., 2000. Hydrodynamics of fishlike swimming. Annu. Rev. Fluid Mech. 32 (1), 33–53. von Ellenrieder, K., Parker, K., Soria, J., 2008. Fluid mechanics of flapping wings. Exp. Therm. Fluid Sci. 32 (8), 1578–1589. Wang, Z., 2005. Dissecting insect flight. Annu. Rev. Fluid Mech. 37, 183–210. Willmott, A.P., Ellington, C.P., 1997. The mechanics of flight in the hawkmoth Manduca sexta. I. Kinematics of hovering and forward flight. J. Exp. Biol. 200 (21), 2705–2722. Wilson, M., Eldredge, J., 2011. Performance improvement through passive mechanics in jellyfish-inspired swimming. Int. J. Non-Linear Mech. 46 (4), 557–567. Witteveen, J., 2010. Explicit and robust inverse distance weighting mesh deformation for CFD. In: 48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, Paper 2010-165, Orlando, Florida (January). Yeh, P.D., Alexeev, A., 2014. Free swimming of an elastic plate plunging at low Reynolds number. Phys. Fluids 26 (5). Young, D., Felgar, R., 1949. Tables of characteristic functions representing normal modes of vibration beam, vol. 44. The University of Texas Publication, p. 4913. Zhang, J., Liu, N.-S., Lu, X.-Y., 2010. Locomotion of a passively flapping flat plate. J. Fluid Mech. 659 (September), 43–68. Zhang, X., Ni, S., Wang, S., He, G., 2009. Effects of geometric shape on the hydrodynamics of a self-propelled flapping foil. Phys. Fluids 21 (10). Zhu, Q., 2007. Numerical simulation of a flapping foil with chordwise or spanwise flexibility. AIAA J. 45 (10), 2448–2457. Zhu, X., He, G., Zhang, X., 7 2014a. How flexibility affects the wake symmetry properties of a self-propelled plunging foil. J. Fluid Mech. 751, 164–183. Zhu, X., He, G., Zhang, X., 2014b. Numerical study on hydrodynamic effect of flexibility in a self-propelled plunging foil. Comput. Fluids 97, 1–20.