A modified ALE method for fluid flows around bodies moving in close proximity

Computers and Structures 145 (2014) 1–11

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

A modified ALE method for fluid flows around bodies moving in close proximity T. Gillebaart a,⇑, W.B. Tay a,b, A.H. van Zuijlen a, H. Bijl a a b

Faculty of Aerospace Engineering, Delft University of Technology, P.O. Box 5058, 2600GB Delft, The Netherlands Temasek Laboratories, National University of Singapore, 5A Engineering Drive 1, 117411, Singapore

a r t i c l e

i n f o

Article history: Received 17 April 2013 Accepted 23 July 2014

Keywords: Clap-and-fling Flapping wings Immersed boundary method Mesh topology changing method ALE Moving meshes

a b s t r a c t A modified ALE method combining a body conforming mesh and a topology changing algorithm is developed to simulate unsteady symmetric fluid problems with bodies in close proximity. By using this method, the advantages of the original ALE method are used, while the disadvantages are reduced or eliminated. A mesh topology changing algorithm is used to create a computational domain boundary, while still having a smoothly deforming mesh in the remaining part of the mesh. It is shown that at higher Reynolds numbers the newly developed method performs better in terms of spurious oscillations than an immersed boundary method (IBM). Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Bodies moving in close proximity of each other is a challenge in many engineering applications. For example cars including ground effects, vortex induced vibrations on piping installations and more recently flapping Micro Aerial Vehicles (MAVs). At Delft University of Technology a MAV using flapping wings called Delfly has been developed [1,2] and both numerical and experimental studies have been performed on flapping wings [3,4]. Its design is inspired by a dragonfly. The wing of the Delfly performs a clap-and-fling motion at both sides of the body. The aerodynamic mechanism present during such a motion is referred to as the Weis-Fogh mechanism (see [5]). In Fig. 1 the clap-and-fling motion is illustrated. During this motion the wings of the Delfly deform significantly, as illustrated in Fig. 12. In hovering conditions the typical Reynolds number is in the order of 104 , which is relatively high for the clap-and-fling motion. Many flapping wing studies have been performed in the past, but only few focus on the clap and fling mechanism as shown in [6], which reviews the recent progress in flapping wing aerodynamics. Experimental studies did reveal the complex aerodynamics around the wings during clap and fling, which is highly dependent on kinematics, Reynolds number and wing design [7–10]. Only a few computational studies have been performed. ⇑ Corresponding author. Tel.: +31 15 27 87527. E-mail address: [email protected] (T. Gillebaart). http://dx.doi.org/10.1016/j.compstruc.2014.07.016 0045-7949/Ó 2014 Elsevier Ltd. All rights reserved.

A two-dimensional clap-and-fling motion at low Reynolds number (Re < 128) with artificial flexibility has been studied [11,12], which showed the potential benefits on the force production of the clapand-fling kinematics. The relation between the vortex shedding and force production for clap-and-fling kinematics at different Reynolds numbers, ranging from 80 to 1400, has been investigated for both two- and three-dimensional computations [13,14]. High Reynolds number flows with bodies in close proximity is a challenge in a wide range of engineering applications. With a special interest in physics behind the Delfly kinematics there are two open issues: clap-and-fling at high Reynolds number (104 ) and influence of wing deformation for the Delfly wings. To numerically study these phenomena three methods can be used: Arbitrary Lagrangian Eulerian (ALE) method, immersed boundary method (IBM) and overset-grids. However, all of them have significant disadvantages when used for simulating flexible clap-and-fling at relatively high Reynolds numbers. In the ALE method the motion of the moving body is interpolated to the mesh points to deform the internal mesh accordingly [15]. An important advantage is a body conforming mesh for each computational time step. Due to this body conforming mesh the boundary layer can be solved accurately and efficiently for high Reynolds numbers, as cells with high aspect ratio aligned with the body can be used. A second important advantage of ALE over the other two methods is that no interpolation of flow variables between meshes or bodies is needed, while both IBM and the overset-grid method interpolate flow variables. However, when

2

T. Gillebaart et al. / Computers and Structures 145 (2014) 1–11

Fig. 1. Illustration of the two-dimensional clap-and-fling motion and the fundamental flow phenomena present. The illustration is obtained from [29].

multiple moving bodies are present, or a body is moving close to a static object (e.g. domain boundary), the ALE method cannot be applied. The mesh quality will decrease rapidly due to extreme deformation of cells between the bodies, causing a larger error, unstable calculations or even an invalid mesh. The extreme case of clap-and-fling motion is impossible to simulate with the standard ALE method, since the bodies are (almost) touching each other. In this cases highly stretched/compressed cells will arise together with negative volumes. Continuous remeshing would be needed to obtain high quality meshes at all times, which comes at the expense of additional computational costs and interpolation between meshes. At low Reynolds numbers IBM has been used for clap-and-fling simulations [12,14]. The most important disadvantage of the IBM is the need for interpolation of the body boundary conditions to the flow. Correctly interpolating the body boundary conditions becomes difficult at higher Reynolds numbers [16]. Numerical instabilities arise, causing oscillations in the pressure distribution on the body. In addition, a higher mesh resolution is needed to capture the boundary layer appropriately compared to a body conforming mesh. Due to the interpolation spurious oscillations occur, especially at higher Reynolds numbers. Two main sources of these spurious oscillations are identified [17]: pressure discontinuity across the immersed boundary and temporal velocity discontinuity near the immersed boundary. Recently, a Cartesian cut-cell method coupled with a virtual merging technique has been introduced [18], which reduces the spurious oscillations significantly at lower Reynolds numbers. However, enforcing the boundary conditions in the immersed boundary method for moving objects remains a cumbersome task, especially at higher Reynolds numbers. Taking into account the future application in incompressible fluid–structure-interaction (FSI), reduction of spurious oscillations is critical for performance, since FSI is very sensitive to (spurious) oscillations in the pressure. Among others, an important advantage of the immersed boundary method is the use of a Cartesian mesh [16]. However, higher Reynolds number flows are less suitable for the immersed boundary method, due to increasing sensitivity of the boundary interpolation at higher Reynolds numbers [16]. ALE should be preferred over IBM, according to [19], because of the body conforming mesh, robustness and accuracy as long as remeshing is not required.

Finally overset-grids have been used to simulate moving bodies. Overset-grids have been used to simulate a clap-and-fling motion [20]. As in the ALE method a big advantage is the body conforming mesh. However, interpolation is still needed at the interface between the overset-grids and the main grid. A comparison is made between the ALE method and the overset-grid method [21]. It is shown that the interpolation between the body mesh and the background mesh is important for accurate results. During clap-and-fling this interpolation is performed close to the wings, since the wings are in close proximity of each other. Although overset-grids uses a body conforming mesh, cross-mesh interpolation of flow variables near the body is undesired. All three methods have disadvantages and advantages when a high Reynolds number clap-and-fling motion needs to be simulated. However, none of the three methods provide the robustness and accuracy needed to simulate high Reynolds number clap-andfling kinematics with flexible wings. The ALE method provides important advantages: body conforming meshes (good for higher Reynolds numbers), low mesh requirements and no interpolation of flow variables required. However, the ALE method does not allow for creation or destruction of cells. Therefore, the goal of this study is to develop a modified ALE method combining a body conforming mesh with a local mesh topology changing algorithm near the symmetry plane to simulate clap-and-fling kinematics with deforming bodies at relatively high Reynolds numbers O(104). This modified ALE method still has the two most important advantages of the ALE method: (1) capturing the boundary layer on the body and (2) direct enforcement of all boundary conditions. Additionally, this modified ALE method uses a topology changing algorithm to directly impose the symmetry boundaries conditions via an actual domain boundary using projection of the cells near the symmetry plane and remapping of data from new to old meshes. As a results bodies can move in close proximity of each other while still using the basis of the ALE method. Three different wing kinematics are used to determine the validity and quality of the new method. The wing kinematics considered in this paper are (1) the low Reynolds number (8–128) clap-and-fling motion using an immersed boundary method [11], and (2) the plunging airfoils in biplane configuration [22] and the Delfly wing motion obtained experimentally in [4]. The low

T. Gillebaart et al. / Computers and Structures 145 (2014) 1–11

Reynolds number clap-and-fling kinematics are used to show that the new method provides reliable results. To show the performance at higher Reynolds numbers two plunging airfoils are simulated at a Reynolds number of 104 and compared to [22], which used overset-grids. Finally, the Delfly wing motion is simulated, which combines the high Reynolds numbers with the small distance between the wings during a clap-and-fling motion. All three kinematics are also simulated using an immersed boundary method for comparison. First the governing equations are stated. Secondly the working principle of the modified ALE method will be explained. Subsequently the basis of the boundary method is presented. After this the kinematics used in this paper are shown. The results for the three kinematics from [11,22,4] are discussed next. Finally the conclusions are drawn. 2. Governing equations For all simulations the unsteady incompressible Navier–Stokes (NS) equations are used and laminar flow is assumed. The equations are written in ALE form and discretized using the finite volume method. The equations are solved using the PISO algorithm [23]. The present method is implemented into OpenFOAM [24], which is used to calculate the unsteady incompressible laminar NS in finite volume formulation. For details on the discretization and implementation of the PISO algorithm see [24]. The unsteady incompressible Navier–Stokes momentum equation in integral ALE formulation is

@ @t

Z

udV þ

I

V C ðtÞ

¼

Z

SC ðtÞ

$p V C ðtÞ

n  ðu  um ÞudS 

q

Z

r  ðm$uÞdV

V C ðtÞ

dV:

ð1Þ

Here V C is the control volume, SC the surface of the control volume, u is the velocity vector, um the mesh velocity, m the kinematic viscosity, p the pressure and q the density. The discretized momentum equation using the finite volume method, a first order time discretization and Gaussian theorem, is

0 1 nþ1 nþ1 n n X X u V u V B C þ nf  ðmruÞf Sf @ /  /m Auf  |{z} |{z} Dt f f X p ¼ nf Sf : f

q

nf uf Sf

In case of incompressible flow the following boundary conditions hold for the symmetry boundary:

ns  u ¼ 0   ~ u  nT ¼ 0 ns  r s

ð4Þ

~ p ¼ 0; ns  r

ð5Þ

ð2Þ

f

Here / is the fluid flux at the face, /m is the mesh flux at the face from the ALE formulation, superscript n þ 1 indicates the current time, superscript n indicates the old time, subscript f indicates the face value, nf is the face normal and Sf is the face area. Two terms are influenced by the mesh motion. First, the temporal term includes the change in volume of the cell over time. Secondly, the convective term has a mesh flux term for each face. The mesh flux is calculated using the swept volume of each face. This is done to satisfy the discrete geometric conservation law (DGCL). Both large translations and rotations are present in the kinematics considered in this study. Two mesh movement strategies are available: using grid-connectivity, such as Laplacian smoothing successfully applied in [25]; or point-by-point deformation, such as Radial Basis Function (RBF) interpolation. It is shown that the RBF interpolation provides a robust and efficient method of deforming the mesh for ALE computations in case of large displacements and rotations of the body [26]. Therefore RBF interpolation with the Thin Plate Spline (TPS) is used for the basic mesh deformation.

ð3Þ

where ns is the normal vector of the symmetry plane. 3. Immersed boundary method All kinematics are also simulated using a typical immersed boundary method to provide a more detailed comparison between a typical immersed boundary method and the modified ALE method. The method by Ravoux et al. is chosen for the comparison, since it is a well known method. It combines ideas from the IBM and the volume-of-fluid (VOF) method [27] to simulate bodies on a Cartesian mesh. An external body force density (f c ), which signifies the presence of the solid body, is introduced into the flow governing equations:

@u 1 ¼ u  ru þ r2 u  rp þ f c : @t Re

ð6Þ

This method, which uses a similar form of fractional step method except for an additional step after solving the momentum equation, is outlined in three steps. First, the discretized momentum equation is solved without the forcing term. Hence, f c ¼ 0. The velocities obtained in this step will not be divergence free. The second step, which is an essential component to Ravoux’s IBM, is to consider the whole system as a ‘‘binary’’ fluid. Hence, one phase is the fluid, while the other is the body. As shown in Fig. 2, the volume fraction of the rigid body phase is unity whereas the ordinary fluid phase has a volume fraction (U) of zero. In between, the medium is partially made up of the fluid and rigid body phases and the volume fraction has a value between zero and unity. By imposing the forcing term (f c ) in cells that are partially or fully occupied by the body, the velocity field is modified to make it vanish in the body. This is achieved by setting

~ 00 ði; jÞ ¼ ½1  Uði; jÞu ~ 0 ði; jÞ; u ~0

nf ðum Þf Sf

3

ð7Þ

where u is the velocity vector obtained from solving the momen~ 00 is the velocity vector modified tum equation in the first step and u due to the forcing term f c , hence, in this scheme, f c is not explicitly defined. Rather, its effect is implicitly added into the overall scheme through the second step. The third step involves solving the Poisson equation to obtain the pressure and updating the velocities to divergence-free. The current solver is also slightly different from the implementation by Ravoux et al. [27] in that a fully implicit Crank Nicholson second

Fig. 2. Definition of volume fraction of a cell used in the in-house IBM.

4

T. Gillebaart et al. / Computers and Structures 145 (2014) 1–11

order (CN2) scheme is used instead of the original fully explicit scheme to improve the stability of the solver. More details about the solver and its algorithm can be found in [27]. 4. Modified ALE method The modified ALE method combines the ALE approach with a mesh changing algorithm near the symmetry plane to overcome the problems of deteriorating mesh quality near the symmetry plane. This method is proposed to simulate a moving body that is in close proximity of a symmetry plane. Three important properties of the modified ALE method are direct enforcement of boundary conditions into the spatial discretization schemes, compliance with the DGCL, mass conservation, but no momentum conservation near the symmetry plane. The modified ALE method is described in 5 steps, which are explained in detail in Section 4.1. The first step is to define the mesh deformation due to the wing motion, as in a regular ALE method. By ignoring the symmetry plane, such that the mesh can move unconstrained through this symmetry plane, a smoothly deforming mesh is created with a high quality mesh near the body. Secondly, the symmetry plane is created by means of projection of points on the symmetry plane. As a third step the mesh topology at the previous time is redefined to ensure a similar mesh topology as the mesh at the current time. After this the velocity and pressure at the previous time are interpolated for the newly entered cells based on surrounding cells and the symmetry boundary conditions. Finally the regular ALE system of equations can be constructed and solved for. Each of these steps are explained detailed in Section 4.1, after which the limitations and performance of the modified ALE method are discussed in Sections 4.2 and 4.3. 4.1. Method The modified ALE method is described in five steps. In Fig. 4 the first four steps are shown schematically and these steps will be discussed in this section. 1. Creation of mesh at t nþ1 (a) RBF mesh interpolation and face identification: The mesh is deformed based on the motion of the wing from tn to t nþ1 while ignoring the symmetry plane (see Fig. 4(a) and (c)) to ensure a smoothly deformed mesh is created with a high quality mesh near the body:

dinternal ¼ Hinterpolation dcontrol

ð8Þ

Here dinternal are the internal point displacements, Hinterpolation is the interpolation matrix from control points to internal points and dcontrol are the control point displacements, which consists out of the farfield point displacements and the wing point displacements. At this stage the symmetry plane is still detached from the mesh and thus not connected to the mesh. Based on the new mesh at t nþ1 the cells belonging to the computational domain are identified. When the cell centre is located to the left of the symmetry plane (cðxcc Þ > 0) it is considered to be in the computational domain, while cells located to the right of the symmetry plane are ignored for this time step:





cðxcc Þ ¼ xcc  xsp  nsp

ð9Þ

where xcc are the cell centre coordinates, xsp is a point on the symmetry plane and nsp is the symmetry plane normal. Faces at the edge of the computational domain are identified and will be used to construct the actual symmetry plane. In Fig. 4(c) these are indicated by the blue dashed dotted line.

(b) Constructing the symmetry plane: The symmetry plane could be treated as a regular immersed boundary, but the close proximity between the wing and symmetry plane still causes the boundary condition interpolation to have a significant undesired influence on this pressure distribution. Therefore the symmetry plane is constructed every time, by projecting points towards the symmetry plane. All the points associated with faces identified as symmetry faces are projected orthogonally to the symmetry plane:

   xprojected ¼ xinternal  nisp xinternal  xisp  nisp ;

ð10Þ

where xprojected are the projected point coordinates associated with the symmetry plane faces and xinternal are the original point coordinates. In Fig. 4(d) the resulting mesh is shown (t nþ1 sym ). An additional algorithm is present during this step to ensure that no zero or negative volume cells are created. The mesh after this step (at tnþ1 sym displayed in Fig. 4(d)) is used in step 3 to solve the laminar incompressible Navier–Stokes equations with the symmetry boundary conditions imposed directly (without additional interpolation) in the numerical schemes.  2. Redefining mesh and flow parameters from tn to tn (a) Redefining the previous mesh: In order to satisfy the DGCL n at tnþ1 sym the mesh at t sym is redefined to obtain a similar mesh nþ1 topology as at tsym . Only the points associated to the symmenþ1 try plane at tnþ1 sym are set to the position at t sym , while all remaining points stay at their original location (t n ):

     xinternal t nsym ¼ xinternal t nsym      xprojected t nsym ¼ xprojected tnþ1 sym

ð11Þ n

In Fig. 4(e) this mesh is shown and is named tsym . This operation changes the mesh topology at the symmetry plane to ensure a similar mesh topology as at t nþ1 sym . (b) Interpolation of velocity and pressure: Cells entering the computational domain due to the mesh movement performed in step 1 (see illustrational case in Fig. 3), need a value for the pressure and velocity at the previous time (t n sym ). RBF interpolation based on the symmetry boundary conditions and the surrounding cell values is used to determine the values of velocity and pressure at tn sym for the newly entered cells. In Eq. (12) the interpolation functions is stated.

sðxÞ ¼

N¼N n þN b X

0

1

ai /i @kx  xi kA þ |fflfflfflfflffl{zfflfflfflfflffl}

i¼1

r

M¼M Xg j¼1

0

1

bj hj @kx  xj kA: |fflfflfflfflffl{zfflfflfflfflffl}

ð12Þ

r

Here N n is the number of neighboring cells used for interpolation, N b is the number of fixed value boundary conditions, Mg is the number of gradient boundary conditions, ai are the weighting coefficients for the basis functions, /i ðr Þ are the radial basis functions, bj are the weighting coefficients for the second set of basis functions and hj ðr Þ are the second set of basis functions. To solve for all weighting coefficients M þ N boundary conditions are needed. In Eq. (13) the fixed value boundary condition is stated:

sðxk Þ ¼ hk ¼

N¼N n þN b X

ai /i ðkxk  xi kÞ þ

i¼1

M¼M Xg

  bj hj kxk  xj k

j¼1

for k ¼ 1; . . . ; N;

ð13Þ

and in Eq. (14) the gradient boundary condition is stated: N¼N n þN b X  d   dr d  0 s xq ¼ hq ¼ ai /i kxq  xi k dg dg dr i¼1

þ

M¼M Xg j¼1

bj

 dr d  hj kxq  xj k for q ¼ 1;. .. ;M; ð14Þ dg dr

T. Gillebaart et al. / Computers and Structures 145 (2014) 1–11

5

(b) t n*sym

(a) t nsym

Fig. 3. Schematic representation of a single cell entering the computational domain using the following symbols: dark grey area are the cells not considered for the current time step, white area is computational domain, dashed purple line is (immersed) symmetry plane, the red square is the newly entered cell with unknown flow values at tn sym , purple triangles are face centers associated to symmetry faces, green triangles are face centers associated to symmetry faces that are used for interpolation, green crosses are cell centers used for the interpolation and black dots are other cell centers. In the green cell centers both velocity and pressure are known and they are used as fixed value condition. For the green triangles on the symmetry plane both the flow values as well as gradient conditions are known. In addition to the fixed value condition also the gradient boundary conditions stated in Eqs. (4) and (5) are used for the green triangles. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Snap

(a) tn

Step 2.(b): Unew & pnew interpolation

(b) tnsym

Step 2.(a): Internal mesh points Step 1.(a): RBF deformation

tn

Step 2.(a): Symmetry points

(e) tn*sym

Step 1.(b): Snap

t n+1

(c) tn+1

(d) tn+11sym

Fig. 4. Schematic representation of the five steps used to implement the modified ALE method using the following symbols: dark grey area are the cells not considered for the current time step, white area is computational domain, dashed purple line is immersed symmetry plane, blue dashed-dotted line represents the symmetry faces, purple triangles are points associated to symmetry faces, red squares are cell centers of the first row of cells outside the computational domain, green crosses are cell centers of the first row of cells in the computational domain and black dots are other cell centers. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

where xk are the positions of the known fixed values, hk are the known fixed values, xq are the positions of the known 0 gradient values, hq are the known gradient values and with

g the direction in which the gradient is taken. The positions of the gradient conditions coincides with some of the fixed value condition positions. Due to this, the basis functions

6

T. Gillebaart et al. / Computers and Structures 145 (2014) 1–11

used for the gradient conditions (hi ðrÞ) must be different from the basis functions used for the fixed value conditions (/i ðr Þ). This is to prevent a singular system of equations. For /i ðr Þ the thin-plate-spline (TPS) is used, while for hi ðrÞ the Gaussian basis function is used (see [26]). For all the newly entered cells this interpolation is performed for the pressure and the 2 (in 2-D) velocity directions. 3. Solving system of equations Finally the system of discrete equations (Eq. (16)) based on the n mesh at tnþ1 sym and t sym can be solved in the computational domain for the velocity and pressure:

X nf  uf Sf ¼ 0

ð15Þ

f

u

nþ1

¼

V

nþ1

n

u V Dt

n

0

X p nf Sf : f

q

4.3. Performance of modified method

1

XB X C þ nf  ðmruÞf Sf @ /  /m Auf  |{z} |{z} f f nf uf Sf

n f ðu m Þf S f

ð16Þ

f

These equations are solved using the PISO algorithm. The symmetry boundary condition is enforced implicitly in the finite volume framework without interpolation because of the previous steps. Thus, the boundary is treated as a regular symmetry boundary and equations are solved using the standard linear solvers available in OpenFOAM as if the mesh only made a regular mesh nþ1 deformation step from t n sym to t sym . 4.2. Limitations of the method Even though the modified ALE method combines the good properties of the ALE approach with a mesh topology changing method, there are three limitations to the current implementation. First, the redefinition of the mesh causes an error in the momentum conservation. Secondly, the relatively large displacement of the cell centers near the symmetry plane during the redefinition step changes the pressure gradient field near the symmetry plane. Finally, the current implementation is limited to two-dimensional and symmetric cases. Redefining the mesh at the previous time step causes a local error in the momentum conservation. During this step the shape, volume and location of specific cells at the symmetry boundary change, while their values for pressure and velocity remain unchanged. Additionally, the flow variable values for the newly entered cells are interpolated. Considering the general formulation of the discrete conservation law, as stated in Eq. (17), it can be concluded that if the volume changes, while the physical quantity (/) remains unchanged, conservation for the momentum may not be preserved: n

n

Nt Nt X X n n /i V ti – /j V tj : i

at a single previous time step is redefined, only time discretization schemes based on a single previous time are available. Additionally, the method is currently only applicable for two-dimensional (2-D) symmetric calculations. However, by applying the method on both sides of the immersed plane and introducing a sliding interface also non-symmetric flows and motions can be simulated. Extension to three-dimensions is possible, but it is a cumbersome task because of the increased complexity of the projection algorithm. Finally, a combination of quadrilateral and pentagon cells are used to create meshes during this study. Some additional work might be needed to ensure a correct projection of the faces to the symmetry plane for meshes with more complex mesh topologies.

ð17Þ

j

However mass is conserved, since the density is constant and the total volume of the computational domain remains unchanged. Due to the redefinition of the mesh at the previous time, cell positions (and thus cell centre locations) change, which alters the pressure gradient field used in the momentum predictor step of the PISO loop. The error in momentum conservation and pressure gradient influences the pressure field locally, creating small jumps in the pressure. Due to the use of the incompressible N–S these small jumps affect the pressure field instantly, which result in small jumps in the force values. Lastly, there are limitations to the applicability of the method and the choice of time discretization schemes. Since only the mesh

In order to determine the performance of the modified ALE method two indicators are used: order behavior and error. By introducing components of the modified ALE method step by step into the solver the effect of the different algorithms on the order behavior and the error is determined. Five intermediate cases are used: non moving wing, moving wing (ALE), moving wing with projection of faces (no symmetry conditions), moving wing with projection and redefinition of previous mesh and finally the complete modified ALE method (moving wing with projection, redefinition and symmetry conditions). These five intermediate cases are named steady, moving, projection, redefinition and modALE, respectively. 4.3.1. Test case The performance is analyzed using the clap-and-fling motion described by the equations from [11]. The following parameters are used: Reynolds number of 10, chord length of 0.1 m, frequency of 2.5 Hz, maximum angle deflection of 30 degrees, translational amplitude of 0.5 times the chord, rotation time is 0.25 of period (summing up to 1 with four rotation phases in a period), both acceleration and deceleration time of 0.25 of the period (summing up to 1 with two acceleration and two deceleration phases), symmetry plane location at 12.5% from initial wing position and 0.4 s simulation time (one period of motion). In Fig. 5 a schematic representation of the domain and the movement is given. For steady the domain is kept the same, but the wing remains static. For moving, projection and redefinition the wing movement is imposed but without enforcing the symmetry conditions. Finally, modALE uses the same movement but with symmetry boundary conditions enforced. 4.3.2. Temporal order Each of the five solvers use five different normalized time steps (normalized by time step of 103 s): 1, 0.5, 0.25, 0.125 and 0.0125. The results obtained with the final time step are used as a reference solution. The unsteady vertical velocity in the region around the symmetry plane location is used to calculate the error relative to the reference solution. The RMS is taken of these errors on which the order is determined. A first order backward-Euler scheme is used in all simulations. The resulting temporal orders are shown in Table 1. Based on this test case the correct time order is present in all parts of the modified ALE method. The order for the full method (last column in Table 1) is most important and is approximately equal to 1 as it should be for the first order backward-Euler time integration scheme. 4.3.3. Spatial order Determining the spatial order in this method is not straightforward, since the modified ALE method can only be applied to unsteady simulations. However the influence of the mesh quality on the spatial order can be assessed. Therefore steady state

T. Gillebaart et al. / Computers and Structures 145 (2014) 1–11 u = 0, v = 0, ∇p = 0

1.5 chord

u=1 v=0 ∇p = 0

2.5 chord

1.0 chord

∇u = 0 ∇v = 0 p=0

1.5 chord

u = 0, v = 0, ∇p = 0

Fig. 5. Schematic representation of domain and wing movement used to analyze the performance of the modified ALE method. Red plate is the wing, green dashed line is the symmetry plane and grey area is the area not considered for computations if the symmetry conditions are applied. Blue arrows indicate movement of the wing and black arrows the size of the computational domain. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

calculations are performed on 10 deformed meshes obtained using the meshes from the modified ALE method results. The result from the finest mesh is used as reference solution. Cell values are interpolated, using a second order method, to a line close to the symmetry plane and a line in the wake. Based on these interpolated values the RMS of the error is calculated. The order is determined based on these errors. Steady, moving, projection and modALE are used to obtain the influence of the different algorithms on the spatial order. Redefinition is not used, since for steady state calculations there is no difference compared to projection. In Table 2 the results are stated. All calculations are performed with second order spatial discretization schemes. Although the order differs slightly for the different cases, there is no clear trend. It can be concluded that changes in the mesh and the mesh quality near the symmetry plane caused by the different algorithms have little influence on the spatial order. 4.3.4. Force error The error in the force is the final indicator of the performance. For moving, projection and redefinition the resulting forces are ideally the same. However, due to the projection of faces and the redefinition of the previous mesh, differences are present between moving and the other two. In Table 3 the normalized RMS of the error with respect to the moving results in the drag is stated for different meshes and time steps for both projection and redefinition. From these results three important conclusions can be derived. First, the error increases significantly when redefining the previous mesh. The larger error in redefinition compared to projection is expected to be caused by the error in the momentum conservation. The error present in the projection is caused by the discontinuous mesh movement due to the projection of faces. This error is mainly caused by the error in the pressure gradient. From this it can be

Table 1 Influence of the modified ALE method on temporal order. The time step sizes are normalized by a time step of 103 s.

D~t ¼ 1 : 0:5 D~t ¼ 0:5 : 0:25 D~t ¼ 0:25 : 0:125

Steady

Moving

Projection

Redefinition

modALE

1.021 1.039 1.079

1.013 1.035 1.076

1.012 1.034 1.071

1.011 1.031 1.069

1.015 1.033 1.075

7

concluded that the most significant error in the drag is caused by the error in the momentum conservation due to the redefinition. Secondly, the error decreases when a smaller time step size is used for both cases. Even though the overall error decreases, it must be stated that the amplitude of the oscillations present in the results increases when the time step is decreased. This is caused by the projection of faces leading to a discontinuous mesh movement. Since the amplitude of the discontinuous mesh movement is independent of the time step size, a smaller time step size will cause a larger mesh velocity and acceleration. This results in a larger amplitude of the oscillations. Finally, the error reduces when a higher spatial resolution is used as well. Since the amplitude of the discontinuous mesh movement scales with the spatial resolution, a higher resolution causes a smaller amplitude and thus a lower error. The final two conclusions (larger amplitude of numerical oscillations when reducing time step and smaller error for increasing mesh size) are comparable with the behavior in immersed boundary methods [17]. Since this method uses a similar approach as the cut-cell finite volume immersed boundary method (see [16]), this behavior is expected. 5. Results Three different wing kinematics are considered to validate the modified ALE method and compare with other methods: low Reynolds number clap-and-fling, high Reynolds number plunging airfoil and high Reynolds number clap-and-fling. The first two cases [11,22] are chosen to make a comparison with literature. However, [11] has a low Reynolds number (128), while in [22] wings are not in close proximity (minimum gap of 0.6 times the chord). Therefore the kinematics from [4], which describe clap-and-fling at high Reynolds number, are used to determine the performance of the modified ALE method compared to the earlier described IBM. In addition, oscillations in the forces for different methods are discussed in detail since they are an important factor for future fluid–structure-interaction simulations. 5.1. Clap and fling at low Reynolds number In order to validate the modified ALE method the low Reynolds number (128) clap-and-fling wing kinematics described in [11] are used and compared to their results, which are obtained using an IBM. The minimum gap between the wings is 1=6 times the chord. In Fig. 6 the normalized translational and rotational velocities are shown. The vertical force coefficient (C Y ) and horizontal force coefficients (C X ) from [11], the IBM discussed previously and the modified ALE method are presented in Figs. 7 and 8. The overall correspondence between the three different methods is good, however, there are differences in the peak amplitudes, phase and numerical oscillations. First, oscillations are found in C X (see Fig. 7 at 55% of the period), which are caused by cells entering the domain. However, this is expected due to the local error in the momentum conservation, change in pressure gradient field, interpolation of flow variables and simplicity of the projection algorithm near the symmetry plane. At the same time this is the most critical part of the period considering the numerical method, because many enter the computational domain and the wing is in close proximity to the symmetry plane. In these new cells interpolation errors are made causing small pressure jumps, which consequently affect the forces on the wing due to its close proximity. These oscillations reduce when the new cell interpolation algorithm is enhanced or a finer mesh is used, as is demonstrated for a different test case in Figs. 9 and 10. The interpolation technique can be enhanced by

8

T. Gillebaart et al. / Computers and Structures 145 (2014) 1–11

Table 2 Influence of the modified ALE method on spatial order. First value is the order based on interpolated values close to symmetry plane location. The value between brackets is the order based on the interpolated values in the wake. N is the total number of cells in the mesh. Spatial order is determined by calculating the error between the results of two consecutive meshes.

N = 2576:10304 N = 10304:41126

Steady

Moving

Projection

modALE

1.402 (1.784) 1.633 (1.907)

1.220 (1.281) 1.565 (1.713)

1.120 (1.247) 1.681 (1.665)

1.169 (1.664) 1.656 (1.740)

Table 3 Errors in drag for different solvers. Three normalized time steps and three mesh sizes are used. Time steps are normalized by a time step of 103 s. Errors are the 2 norms of the transient forces normalized by number of data points and the band of the drag coefficient (maximum amplitude – minimum amplitude). N = 10,304

D~t ¼ 0:1

D~t ¼ 1

D~t ¼ 0:5

D~t ¼ 0:1

D~t ¼ 0:1

0.030E3 0.152E3

0.014E3 0.133E3

0.010E3 0.097E3

0.007E3 0.052E3

0.002E3 0.022E3

Dimensionless velocities

Translation velocity

Miller−Peskin 2005

5 4 3

Angular velocity

2 1

1.5

0

1

−1

0

2

4

6

8

10

τ [−]

0.5

Fig. 8. Vertical force coefficient obtained using the mesh topology changing algorithm for the clap and fling motion described in [11]. ‘mALE’ indicates the result obtained with the modified ALE algorithm, ‘IBM’ is the immersed boundary method and ‘Miller–Peskin’ is the result presented in [11].

0 −0.5 −1 −1.5

0

2

4

6

8

10

τ [−] Fig. 6. Dimensionless translational and angular velocity for clap-and-fling motion described in [11].

mALE

IBM

Miller−Peskin 2005

12 10 8

CX [−]

IBM

N = 41,126

CY [−]

Projection Redefinition

N = 2576

mALE

6 4 2 0 0

2

4

6

8

10

τ [−] Fig. 7. Horizontal force coefficient obtained using the mesh topology changing algorithm for the clap and fling motion described in [11]. ‘mALE’ indicates the result obtained with the modified ALE algorithm, ‘IBM’ is the immersed boundary method and ‘Miller–Peskin’ is the result presented in [11].

considering mass and momentum conservation, while incorporating the discontinuous change in cell centre position due to the projection. The results from the immersed boundary method used in [11] do not show oscillations. However, the result presented in [11] use the same immersed boundary method as presented in [28]. In [28] all results were filtered to suppress the ‘noise’ in the numerical data and it is assumed that in [11] also the same filtering is used. Due to this filtering the frequency and amplitude of the oscillations cannot be compared accurately. Secondly, the absolute value of the largest peak in C X (Fig. 7) differs significantly (10%) and a small shift in time is present (1% of the period) in the results. A similar shift in time is found in the C Y (Fig. 8). These differences can be expected due to the difference in method and mesh. In Figs. 9 and 10 the results for the two-wing fling motion is shown for a coarse and fine mesh compared to mesh study results from [11]. From these results it can be seen that the mesh used to obtain the results presented in [11] is a relative coarse mesh (1230  1230 cells), which shows a maximum peak difference of approximately 8% with the finer mesh (2460  2460 cells). In the present study a similar mesh study, as done in [11], is performed for the fling motion. The resulting C X and C Y are presented in Figs. 9 and 10, respectively. A good correspondence with the fine mesh results from [11] is found for both the coarse and fine mesh used during this study. Differences between the results obtained in this study and [11] can thus partly be explained by differences in mesh topology and mesh resolution. A clear explanation for the shift is not found. However, the results of C X , obtained with both the modified ALE method and the IBM, show clear changes in the force data at the exact moment the kinematics (and thus the acceleration) change, while in [11] there is a small delay before the effect of the change

9

T. Gillebaart et al. / Computers and Structures 145 (2014) 1–11

Coarse

Fine

MP Coarse

MP Fine

6 5

CX [−]

4 3 2 1 0

0

1

2

3

4

5

τ [−] Fig. 9. Horizontal force coefficient for the fling motion described in [11] using a coarse (16,144 cells) and fine mesh (61,506 cells). The unsteady force coefficient presented in [11] for both their coarse (1230  1230) and fine mesh (2460  2460) is included for comparison. These results are indicated by ‘MP’.

Coarse

Fine

MP Coarse

MP Fine

3 2.5

CY [−]

2

5.2. Plunging airfoils at high Reynolds number Secondly, a higher Reynolds number flow is used for comparison to determine the performance of the modified ALE method at higher Reynolds numbers. In [22] plunging airfoils in biplane configuration are simulated with the unsteady laminar incompressible Navier–Stokes equations at a Reynolds number of 104 by using overset-grids. The plunging motion is described by a sinusoidal function with an amplitude of 0.4 times the chord, period of p seconds and a minimum gap of 0.6 times the chord. Both the IBM and modified ALE method are used to simulate this motion. As mentioned in [16] the mesh size needed for immersed boundary methods is considerably larger when compared to the ALE methods. For the IBM method a mesh of 1972  1824 (3,596,928) cells is used, while for the modified ALE method a mesh of 71,068 cells is used. The mesh used in [22] has a size of 44,112 cells. The lift and drag coefficients are shown in Fig. 11 including the results from [22]. At this higher Reynolds number (104 ) the modified ALE method performs similar as the overset-grid method from [22] in terms of absence of numerical oscillations. Compared to the results obtained with the immersed boundary method, the modified ALE method outperforms the current implementation of this method at this relatively high Reynolds number. After the start-up effects have reduced (after 3 periods) an almost periodic solution is found for both the results from [22] and the modified ALE method. The amplitudes of the peaks do correspond well for these two results, which indicates that the method performs well at higher Reynolds numbers. In the IBM result numerical oscillations are present. Even though the peak values of the immersed boundary method do correspond reasonably well with the results of [22], the significant numerical oscillations in the IBM results make the method less suitable for future fluid–structure-interaction computations and high Reynolds number calculations.

1.5

5.3. DelFly II wings in hovering

1 0.5 0

0

1

2

3

4

5

τ [−] Fig. 10. Vertical force coefficient for the fling motion described in [11] using a coarse (16,144 cells) and fine mesh (61,506 cells). The unsteady force coefficient presented in [11] for both their coarse (1230  1230) and fine mesh (2460  2460) is included for comparison. These results are indicated by ‘MP’.

Finally, a high Reynolds number flow in combination with a small minimum gap is used to determine the performance of the modified ALE method in an extreme case. The two-dimensional chord-wise deforming wing shapes from the Delfly performing clap-and-fling have been obtained using PIV experiments [4]. During the experiments the following parameters were used: flapping frequency of 11 Hz, chord length of 7.4 cm at 71% of the span and a Reynolds number (based on the frequency and the stroke path mALE

IBM

Tuncer 2003

0.2 0 −0.2 −0.4

CD [−]

in the kinematics is seen. Therefore it is expected that a small shift in time is present in the results of [11], which might be caused by their filtering technique. Considering all differences between the studies (mesh, mesh convergence, method) the results do correspond well and the modified ALE method does perform well during clap-and-fling at this Reynolds number. Additionally the IBM is compared with the modified ALE method for the force coefficients in Figs. 7 and 8. At this low Reynolds number (128) the oscillations in the IBM are expected to be small. In both force coefficients some small oscillations are present during the first and last phase of the motion, as is shown in Figs. 7 and 8. Nevertheless, the amplitudes of the oscillations are small for both the IBM and modified ALE method. Hence the results for both methods are acceptable when comparing the overall correspondence with the results from [11]. For the IBM a mesh of 2681  1456 is used. It is clear that for similar results the mesh size needed for each of the two immersed boundary methods is significantly larger compared to the modified ALE method.

−0.6 −0.8 −1 −1.2 −1.4 0

5

10

15

Time [s] Fig. 11. Drag coefficient for plunging airfoil motion described in [22]. Also included are the results from [22] using the an overset-grid method and results obtained using the IBM.

T. Gillebaart et al. / Computers and Structures 145 (2014) 1–11

length) in the order of 104 . This resulted in 50 wing shapes at 71% of the span for a single period. A continuous function of the wing shape in time is created by means of Fourier series interpolation. In Fig. 12 the interpolated wing shapes are shown for several instances during the clap-and-peel cycle. The mesh sizes used in the present method were 246,024 and 59,906 cells, while for the IBM method meshes of 2,182,304 cells and 545,576 cells were used. Simulations are performed on a fixed Courant number of 0.4. In Figs. 13 and 14 the force coefficients during the first period, respectively C Y and C X , are presented. A grid independent solution is found for the modified ALE method for the first part of the motion. However, at the end of the motion the chaotic effects due to the circular pattern of the flow become apparent and the results deviate. For the IBM the initial response is not shown, because large oscillations obscure the result. However, after the start up phase the two results are similar in global behavior except for the final part where again the chaotic effects dominate the flow. The details of this chaotic behavior will be treated in future studies focussing on the physics of the clap-and-peel motion. Despite this chaotic behavior and resulting differences, the performance of the two methods in high Reynolds number cases are illustrated well within these results.

Fig. 12. Interpolated chordwise clap (right) and peel (left) motion of the DelFly II wing at 71% of the span. Only the motion and deformation of the left wing is recorded and shown in this figure. This varying wing shape is mirrored to obtain the motion of the right wing. Due to this symmetric assumption the modified ALE method can be used. The red dashed line represents the virtual symmetry plane.

mALE

mALE− coarse

IBM

IBM − coarse

mALE− coarse

IBM

IBM − coarse

30 20 10 0 −10 −20 −30

0

0.2

0.4

0.6

0.8

1

Time [−] Fig. 14. Horizontal force coefficient for the chordwise deforming wing profile of the DelFly II making a clap-and-fling motion. Results are shown for both the modified ALE method and the immersed boundary method. The two meshes used for the modified ALE method have 246,024 cells (fine) and 61,506 (coarse), while the two meshes used for IBM have 2,182,304 cells (fine) and 545,567 cells (coarse).

Due to the small gap between the wings and a high Reynolds number differences between the two method become apparent. The modified ALE method shows a clear advantage over the IBM: Even though small force jumps are observed in the results obtained with the modified ALE method, the results obtained using the IBM show larger and more oscillations at all phases of the motion. Especially when fluid–structure-interaction would be applied, these oscillations would have a devastating effect on the stability, since the thin foil wings of the Delfly will respond directly to the pressure fluctuations. Besides the stability problems also the accuracy of the resulting foil is affected by these oscillations. Considering this future application of these methods with flexible wings, the oscillations should be reduced as much as possible. For the modified ALE method oscillations are present at the most critical part of the motion: at the start of out-stroke. During this phase many cells enter the computational domain and cause small pressure fluctuations directly causing small force oscillations. However, on a global scale these oscillations are small and the results from the modified ALE method show a smoother result when compared to the IBM results.

6. Conclusions

20

15

CY [−]

mALE

CX [−]

10

10

5

0 0

0.2

0.4

0.6

0.8

1

Time [−] Fig. 13. Vertical force coefficient for the chordwise deforming wing profile of the DelFly II making a clap-and-fling motion. Results are shown for both the modified ALE method and the immersed boundary method. The two meshes used for the modified ALE method have 246,024 cells (fine) and 61,506 (coarse), while the two meshes used for IBM have 2,182,304 cells (fine) and 545,567 cells (coarse).

A modified ALE method is introduced based on the ALE method and a mesh topology changing algorithm near a symmetry plane for two-dimensional flow cases with moving objects in close proximity of a symmetry plane. The numerical performance of the modified ALE method is analyzed by determining the spatial order, temporal order and the error in the drag by introducing the method step-by-step, creating five intermediate cases. A small reduction is found for the spatial order when compared to a single moving wing (ALE). First order temporal accuracy is maintained in all intermediate cases. The observed error in the drag in the projection part is caused by the discontinuous mesh movement resulting from the projection of faces to the symmetry plane. Projecting faces towards the symmetry plane has a limited influence in the final modified ALE method. The largest contribution to the error in the drag is due to the redefinition of the mesh, velocity and pressure at the previous time step. The error made in the momentum conservation is most probably the cause of this. Resolving this issue is part of ongoing studies. A higher mesh resolution or a

T. Gillebaart et al. / Computers and Structures 145 (2014) 1–11

smaller time step reduces the error. However, it must be stated that the amplitude of the numerical oscillations in the forces increases when smaller time steps are used, while the amplitudes decrease for higher mesh resolutions. Since the amplitude of the discontinuous mesh movement scales inversely with mesh resolution, but does not scale with time step size, a smaller time step causes the induced mesh velocities (amplitude over time) to increase and thus the amplitude of the oscillations to increase. First, the modified ALE method is validated for low Reynolds number (128) clap-and-peel kinematics. Secondly, a symmetric plunging airfoil motion at Reynolds number 104 is used to validate the method at higher Reynolds numbers. Based on the results of both these motions, the modified ALE method performs well in terms of accuracy and validity. Differences between the results are explained by differences in mesh, mesh convergence and method. In the full modified ALE method large reductions in oscillations are found, especially compared to IBM. At lower Reynolds number (128) there is no significant difference between the methods. However, the modified ALE method is able to handle deforming object in close proximity of a symmetry plane at moderate high Reynolds numbers (104 ), demonstrated by the Delfly wing motion, while the IBM shows significantly larger oscillations in the forces. For the purpose of simulating the chord-wise deforming crosssection of the Delfly II in hover the performance is good and no extension is needed since the problem is symmetric and twodimensional. Thus this method will be used to perform further studies into the physics of the deforming flapping wings performing clap-and-peel as seen in the Delfly II during hovering. Nevertheless, the present applicability is limited to 2D symmetric problems and further development and improvement is possible. Extension to three-dimensional simulations is possible even though this increases the complexity in the projecting algorithm. The conditions at the detached plane can be changed to function as a sliding interface with two computational domains at each side projecting faces on the this plane. Cases with objects moving nonsymmetric in close proximity of each other in a non-symmetric flow could then also be simulated.

References [1] de Croon G, de Clercq K, Ruijsink R, Remes B, de Wagter C, et al. Design, aerodynamics, and vision-based control of the DelFly. Int J Micro Air Veh 2009;1(2):71–97. [2] De Clercq K, de Kat R, Remes B, van Oudheusden B, Bijl H. Flow visualization and force measurements on a hovering flapping-wing MAV DelFly II. In: 39th AIAA fluid dynamics conference, no. June, 2009. p. 1–10. [3] Bos F, Van Oudheusden B, Bijl H. Wing performance and 3D vortical structure formation in flapping flight. J Fluid Struct 2013;42:130–51. [4] Groen M, Bruggeman B, Remes B, Ruijsink R, van Oudheusden B, Bijl H. Improving flight performance of the flapping wing MAV DelFly II. In: International micro air vehicle conference and flight competition. Braunschweig, Germany; 2010. p. 1–17. [5] Weis-Fogh T. Quick estimates of flight fitness in hovering animals, including novel mechanisms for lift production. J Exp Biol 1973;59:169–230.

11

[6] Shyy W, Aono H, Chimakurthi S, Trizila P, Kang C-K, Cesnik C, et al. Recent progress in flapping wing aerodynamics and aeroelasticity. Prog Aerosp Sci 2010;46(7):284–327. http://dx.doi.org/10.1016/j.paerosci.2010.01.001. [7] Lehmann F-O, Sane S, Dickinson M. The aerodynamic effects of wing–wing interaction in flapping insect wings. J Exp Biol 2005;208(Pt 16):3075–92. http://dx.doi.org/10.1242/jeb.01744. [8] Lehmann F-O, Pick S. The aerodynamic benefit of wing–wing interaction depends on stroke trajectory in flapping insect wings. J Exp Biol 2007;210(Pt 8):1362–77. http://dx.doi.org/10.1242/jeb.02746. [9] Lehmann F-O. When wings touch wakes: understanding locomotor force control by wake wing interference in insect wings. J Exp Biol 2008;211(Pt 2):224–33. http://dx.doi.org/10.1242/jeb.007575. [10] Sohn M, Chang J. Flow visualization and aerodynamic load calculation of three types of clap-fling motions in a Weis-Fogh mechanism. Aerosp Sci Technol 2007;11(2–3):119–29. http://dx.doi.org/10.1016/j.ast.2006.10.003. [11] Miller L, Peskin C. A computational fluid dynamics of ‘clap and fling’ in the smallest insects. J Exp Biol 2005;208(Pt 2):195–212. http://dx.doi.org/ 10.1242/jeb.01376. [12] Miller L, Peskin C. Flexible clap and fling in tiny insect flight. J Exp Biol 2009;212(19):3076–90. http://dx.doi.org/10.1242/jeb.028662. [13] Kolomenskiy D, Moffatt H, Farge M, Schneider K. Vorticity generation during the clap-fling-sweep of some hovering insects. Theoret Comput Fluid Dyn 2009;24(1–4):209–15. http://dx.doi.org/10.1007/s00162-009-0137-2. [14] Kolomenskiy D, Moffatt H, Farge M, Schneider K. Two- and three-dimensional numerical simulations of the clap-fling-sweep of hovering insects. J Fluids Struct 2011;27(5-6):784–91. http://dx.doi.org/10.1016/j.jfluidstructs.2011. 05.002. [15] Donea J, Giuliani S, Halleux J. An arbitrary Lagrangian–Eulerian finite element method for transient dynamic fluid-structure interactions. Comput Methods Appl Mech Eng 1982;33:689–723. [16] Mittal R, Iaccarino G. Immersed boundary methods. Annu Rev Fluid Mech 2005;37(1):239–61. http://dx.doi.org/10.1146/annurev.fluid.37.061903. 175743. [17] Lee J, Kim J, Choi H, Yang K. Sources of spurious force oscillations from an immersed boundary method for moving body problems. J Comput Phys 2011;230:2677–95. [18] Seo J, Mittal R. A sharp-interface immersed boundary method with improved mass conservation and reduced spurious oscillations. J Comput Phys 2011; 230:7347–63. [19] van Loon R, Anderson P, Vandevosse F, Sherwin S. Comparison of various fluidstructure interaction methods for deformable bodies. Comput Struct 2007;85(11–14):833–43. http://dx.doi.org/10.1016/j.compstruc.2007.01.010. [20] Mao S, Xin Y. Flows around two airfoils performing fling and subsequent translation and translation and subsequent clap. Acta Mech Sin 2003;19(2):103–17. http://dx.doi.org/10.1007/BF02487671. [21] Sengupta T, Suman V, Singh N. Solving Navier–Stokes equation for flow past cylinders using single-block structured and overset grids. J Comput Phys 2010;229:178–99. [22] Tuncer I, Kaya M. Thrust generation caused by flapping airfoils in a biplane configuration. J. Aircraft 2003;40:509–15. [23] Issa R. Solution of the implicitly discretised fluid flow equations by operator splitting. J Comput Phys 1985;62:40–65. [24] Jasak H. Error analysis and estimation for the finite volume method with applications to fluid flows, Ph.D. thesis. Imperial College of Science, Technology and Medicine; 1996. [25] Bathe KJ, Zhang H, Ji S. Finite element analysis of fluid flows fully coupled with structural interactions. Comput Struct 1999;72(1–3):1–16. [26] de Boer A, van der Schoot M, Bijl H. Mesh deformation based on radial basis function interpolation. Comput Struct 2007;85:784–95. http://dx.doi.org/ 10.1016/j.compstruc.2007.01.013. [27] Ravoux J, Nadim A, Haj-Hariri H. An embedding method for bluff body flows: interactions of two side-by-side cylinder wakes. Theoret Comput Fluid Dyn 2003;16:433–66. [28] Miller L, Peskin C. When vortices stick: an aerodynamic transition in tiny insect flight. J Exp Biol 2004;207(Pt 17):3073–88. http://dx.doi.org/10.1242/ jeb.01138. [29] Sane S. The aerodynamics of insect flight. J Exp Biol 2003;206(23):4191–208. http://dx.doi.org/10.1242/jeb.00663.