Numerical simulation of a flapping four-wing micro-aerial vehicle

Journal of Fluids and Structures 55 (2015) 237–261 Contents lists available at ScienceDirect Journal of Fluids and Structures journal homepage: www.elsevier.com/locate/jfs Numerical simulation of a flapping four-wing micro-aerial vehicle W.B. Tay a,b,n, B.W. van Oudheusden a, H. Bijl a a b Faculty of Aerospace Engineering, Delft University of Technology, Kluyverweg 1, 2629 HS Delft, The Netherlands Temasek Laboratories, National University of Singapore, 5A Engineering Drive 1, 117411, Republic of Singapore a r t i c l e i n f o abstract Article history: Received 4 June 2014 Accepted 3 March 2015 Available online 21 April 2015 A three-dimensional numerical simulation of a four-wing (two wings on each side, one on top of another) flapping micro-aerial vehicle (FMAV), known as the Delfly micro, is performed using an immersed boundary method Navier–Stokes finite volume solver at Reynolds numbers of 5500 (forward flight condition). The objective of the present investigation is to gain an insight to the aerodynamics of flapping wing biplane configuration, by making an analysis on a geometry that is simplified, yet captures the major aspects of the wing behavior. The fractional step method is used to solve the Navier–Stokes equations. Results show that in comparison to the Delfly II flapping kinematics (a similar FMAV configuration but smaller flapping stroke angles), the Delfly-Micro flapping kinematics provides more thrust while maintaining the same efficiency. The Delfly-Micro biplane configuration generates more lift than expected when the inclination angle increases, due to the formation of a uniform leading edge vortex. Estimates of the lift produced in the forward flight conditions confirm that in the current design, the MAV is able to sustain forward flight. The potential effect of wing flexibility on the aerodynamic performance in the biplane configuration context is investigated through prescribed flexibility in the simulations. Increasing the wing' spanwise flexibility increases thrust but increasing chordwise flexibility causes thrust to first increase and then decrease. Moreover, combining both spanwise and chordwise flexibility outperforms cases with only either spanwise or chordwise flexibility. & 2015 Elsevier Ltd. All rights reserved. Keywords: Flapping wings MAV Immersed boundary methods Delfly Biplane 1. Introduction The Delfly-Micro (DFM) is a flapping wing micro-aerial vehicle (FMAV) developed at the Delft University of Technology (De Croon et al., 2009). It is the latest addition to the Delfly family of ornithopters, which so far consists of the Delfly I and II, with a wingspan of 50 and 28 cm respectively, as shown in Fig. 1.1 The DFM, the smallest of the three, weighs only 3 g and has a wingspan of 10 cm, making it the smallest flying ornithopter carrying a camera. The unique feature of the Delfly MAVs, which differentiate them from other FMAVs (see e.g. Keennon et al., 2012; Pornsin-sirirak et al., 2000) is that they have two pairs of wings, one on top of the other, instead of one (biplane flapping configuration). It therefore generates more thrust compared to a single pair of wings, and gives minimal rocking amplitude (vertical oscillation, perpendicular to the FMAVs flight path), which is a beneficial property for FMAVs to be used as a camera n Corresponding author. Tel.: þ65 6516 7330. E-mail address: tsltaywb@nus.edu.sg (W.B. Tay). 1 More information can be obtained from http://www.delfly.nl. http://dx.doi.org/10.1016/j.jfluidstructs.2015.03.003 0889-9746/& 2015 Elsevier Ltd. All rights reserved. 238 W.B. Tay et al. / Journal of Fluids and Structures 55 (2015) 237–261 Fig. 1. Different versions of Delfly: I (right), II (left), Micro (middle). Table 1 Parameters of the DF2 and DFM under typical forward flight conditions. DF2 DFM Wingspan (cm) Wing root chord (cm) Forward speed (m/s) Re Stroke angle (deg) Flapping frequency (Hz) 28.0 10.0 8.1 2.89 7.0 3.0 36,000 5500 22.0 25.5 11.0 37.5 platform. Due to the close distance between the upper and lower wings during the instroke, an important lift enhancing mechanism known as the clap-and-peel motion (Ellington, 1984) is active. This motion, which is actually a variant of the clap-and-fling motion (Weis-Fogh, 1973) known from insect flight, further helps to improve the thrust generation. Of the three Delfly models, the Delfly II (DF2) is the most established and well-documented platform as a substantial number of experiments and simulations has been performed upon this configuration. Tay et al. (2014) used the immersed boundary method (IBM) (Mittal and Iaccarino, 2005) to perform numerical simulations on a simplified DF2 model in forward-flight configuration with either rigid wings or with a prescribed spanwise deformation at Reynolds number (Re) of 1000 and 5000. It was shown that the biplane wing configuration produces more than twice the average thrust of only the equivalent upper wing of DF2, confirming the thrust enhancing effect of the wing interaction. However, DF2's average thrust is only 40% of that of the upper wing when it is made to flap at twice the stroke angle amplitude (θ0), which indicates the potential for thrust enhancement by increasing the stroke angle. Although the latter result indicates that larger thrust can be generated by using, instead of the two pair of wings, only a single pair of wings flapping over the full stroke angle, the increased stability due to the DF2's smaller lift and moment variation makes it more suited as a camera platform. For that reason, the DFM retains the biplane configuration, but with an increased stroke angle for increasing the lift. The study further revealed that increasing the body inclination angle generates a spanwise uniform LEV instead of a conical one along the wingspan, which is accompanied by higher lift. Lastly, increasing the spanwise flexibility of the wings increases the thrust slightly but decreases the efficiency. In comparison to the extensive analysis (numerical and experimental) of the DF2 configuration, to the best of the authors' knowledge, there are only two studies specific to a flapping wing MAV of a configuration comparable to the DFM. A hummingbird-inspired, flexible wing FMAV developed by Nakata et al. (2011), which is very similar to the DFM, formed the subject of a combined numerical and experimental study. In the numerical flow simulation, an overset solver was applied for simulating the FMAV under hovering conditions. A strong negative pressure region associated to a leading edge vortex (LEV) was found on the upper and lower wings during both of the half strokes. In the experimental investigation, Nakata et al. performed wind tunnel experiments on their FMAV under forward flight conditions. It was found that the biplane wing configuration may or may not be better compared to the single wing pair configuration. The former generated twice as much lift compared to the latter configuration at body inclination angles larger than 401. But if softer Mylar wings were used in the single wing pair configuration, it generated more lift at body inclination angles larger than 301 compared to the biplane wing configuration. Deng et al. (2013) performed experimental investigation on the DFM by measuring the force generated and studying the effect of wing flexibility. Results show that while a relatively rigid wing can produce more force, but at a lower efficiency. Sustained flight at 2–3 m/s is possible at a flapping frequency 28–34 Hz and angle of attack between 251 and 361. As mentioned earlier, despite similarities in terms of appearance between DF2 and DFM, there are several differences between the two FMAVs, as documented in Table 1. Being much smaller, DFM's Re is only 5500 (based on forward speed and chord length), indicating that its flow field is mostly laminar, whereas for DF2, Re is about 36 000 and the flow is likely to be transitional. Due to the larger flapping stroke amplitude, the wings of the DFM touch one another not only during the W.B. Tay et al. / Journal of Fluids and Structures 55 (2015) 237–261 239 instroke of the flapping cycle at the sides (as in the DF2), but also during the outstroke at the top. Hence, the clap-and-peel mechanism occurs at the sides during the instroke and on top during the outstroke. These differences show that their aerodynamics behaviors are different and this motivates to perform numerical investigation specifically on the DFM. The exact numerical simulation of the DFM would be extremely challenging, as it requires the capability of FSI and membrane modeling. An alternative analysis approach is to capture the actual geometric deformation and kinematics of the DFM's wings using multiple cameras and introduce this information as boundary conditions in the numerical simulation, as done by Nakata et al. (2011). However, the analysis of the simulation results will become even more complex and it will be difficult to assess the separate influence of individual variables (such as wing kinematics, spanwise/chordwise deformation) on the efficiency, lift and thrust of the DFM. In view of these considerations, we propose to analyze the effects of different variables on the aerodynamics of the DFM in a more “controlled” manner, and from a more fundamental viewpoint. Moreover, numerical simulations on the effect of spanwise/chordwise flexibility have been mostly limited to the standard FMAV with a single pair wings. In this study, we will investigate the effect of kinematics and flexibility in the context of the biplane DFM configuration. Hence, the main objective of this study is to obtain an estimate of the aerodynamic performance of the DFM from 3D flow simulations on a geometrically simplified model, which still allows us to capture the essential flow features. The simplifications used in this simulation study include using thin plates instead of flexible membrane wings. Also, the distance between the upper and lower wings at the root is slightly increased, in order to prevent wing intersection in the simulations. By comparing the lift generated by the simplified model with the actual DFM, we will show that our model can be considered aerodynamically similar to the actual DFM. However, for easy reference, we will still refer to the simplified model as the DFM. An IBM Navier–Stokes solver (Tay et al., 2012) is used for the simulations. The motivation for choosing IBM is that in the simulations, besides having two wings moving independently, they also come into close proximity to each other, which is a situation for which the IBM approach is very well suited, as will be discussed below. In our study, we intend to characterize not only the current MAV geometry, but also to explore variations that may assist in further improvement of the design. Investigating each variation (such as spanwise flexibility) individually also helps to single out the exact cause and effect. In this context, the following objectives are pursued: 1. To assess the possible aerodynamic benefit of the larger flapping stroke angle of DFM, as compared to that of the DF2. One expects that the thrust will improve with the larger flapping stoke angle, but how will it affect the efficiency and other factors? 2. To evaluate and compare the effects of the body inclination angle, since during forward flight of FMAVs like the DFM, the body inclination angle is always non-zero to vector the forces, such as to produce sufficient thrust and lift at the same time. 3. The addition of prescribed spanwise and chordwise wing deformation will allow us to study and compare the effect of each type of flexibility with respect to rigid wings without the complexity of the fluid structure interaction (FSI). Past studies (Aono et al., 2009; Heathcote et al., 2008; Tay and Lim, 2010; Yang et al., 2012; Zhu, 2007) have found that depending on different conditions, spanwise or chordwise flexibility may or may not be beneficial to the performance of the FMAV. Hence, this motivates to investigate the effects of flexibility in the context of the biplane configuration. The result will be analyzed in terms of the thrust (ct), lift (cl) coefficient, propulsive efficiency (η) and power input, while the flow itself will be visualized by means of pressure and Q criterion iso-surface contour plots. It is expected that the results obtained from this study will further enhance the understanding of the aerodynamic behavior of DFM and as such be beneficial to their further improvement. 2. Numerical method The unsteady viscous flow generated by the flapping wings at relatively low Reynolds number is computed by solving the non-dimensional Navier–Stokes equations using the fractional step method together with the IBM approach. The IBM approach is chosen for a number of reasons. First of all, the biplane wings move in close proximity with respect to one another. Hence, it is problematic to use conforming grids with an Arbitrary Lagrangian–Eulerian (ALE) formulation. Maintaining grid quality with the ALE method is challenging, especially when the distance between the wings becomes progressively smaller. Moreover, the moving wings require the re-calculation of the grid at every timestep, which is expensive in 3D simulations. In contrast, the IBM approach employs a stationary Cartesian grid, eliminating the need to move or deform the grid. The body outlines simply cut through the grid, therefore simulating a complicated or moving body is not a problem. A summary of the numerical method is given in the next section. The reader may refer to the papers by Yang and Balaras (2006), Kim et al. (2001) and Liao et al. (2010) for a more detailed description of the IBM algorithm. 2.1. Discrete forcing IBM and fractional step method The discrete forcing (Mittal and Iaccarino, 2005) IBM is based on a combination of the methods developed by Yang and Balaras (2006), Kim et al. (2001) and Liao et al. (2010). In this method, the momentum equation is modified with the 240 W.B. Tay et al. / Journal of Fluids and Structures 55 (2015) 237–261 addition of the forcing term (fc) to become: ∂u ¼ ∂t u U ∇u þ 1 2 ∇ u Re ∇p þ f c; ð1Þ where u is the velocity vector, t is the time, p is the pressure and Re is the Reynolds number. Eq. (1) has been nondimensionalized using the far field incoming velocity (in forward flight, Uref) and root chord length (c) as the reference velocity and length respectively. Accordingly, the calculation of Re is based on Uref and c. The forcing term is added to simulate the presence of a flow boundary inside the computational domain (Mohd-yusof, 1997). Rewriting Eq. (1), we can obtain fc as fc ¼ ∂u þ u U ∇u ∂t 1 2 ∇ u þ∇p: Re ð2Þ One way to evaluate fc is to take a preliminary explicit “predictor” step by using values from the previous time step (Kim et al., 2001). In the current study, we discretize the right hand side of Eq. (1) explicitly using the forward Euler and 2nd order Adam–Bashforth (AB2) schemes for the viscous and convective terms respectively. A finite volume fractional step method, which is based on an improved projection method (Kim and Choi, 2000), is then used to solve the modified nondimensionalized incompressible Navier–Stokes equations (1) and (3): ∇ U u ¼ 0: ð3Þ The time integration scheme uses a semi-implicit scheme, such that the convective and viscous terms use the second order AB2 and Crank Nicolson (CN2) discretization respectively. The convective and viscous spatial derivatives are discretized using QUICK (Leonard, 1979) and second order central differencing respectively on a staggered grid. Hence, for the calculation of fc, an explicit scheme is used. On the other hand, for the time integration, a semi-implicit scheme is used. The equations to be solved based on the fractional step method are given in discrete form from Eqs. (4)–(7) for completeness. No turbulence model has been added: δu^ i 1X v F ðδu^ i Þ ¼ Vol þ 2 Δt uni u^i ¼ uni þ 1 pn ni A Δt X n p ni A; Vol X ∂pn þ 1 ∂xi X F v ðuni Þ 3X c n n 1X c n F ðui ; U Þ þ F ðui 2 2 1 ; Un 1 Þ; ð4Þ ð5Þ 1 X n ui ni A; Δt ð6Þ Δt X n þ 1 p ni A; Vol ð7Þ ni A ¼ uni ¼ X where δu^ i ¼ u^ i uni , u^ i and uni are the intermediate velocities, uni , Δt, A and Vol are velocities at the n time step, computational time step, cell face area and cell volume respectively. Fv and Fc are the viscous and convective flux respectively and they are defined as F v ðuÞ ¼  1 ∇u U n^ Af Re ð8Þ and F c ðu; U f Þ ¼ U f Af u; ð9Þ ^ Uf and Af are the unit normal, outward face velocity perpendicular to the cell face and cell face area respectively. where n, The discretized momentum equation is first solved to obtain a velocity field which is non-divergence free. Using the nondivergence free velocity, we solve the Poisson equation to obtain the pressure field, which in turn updates the velocity to be divergence free. The momentum and Poisson equations are solved using the open source PETSc (Balay et al., 1997) and HYPRE (Falgout et al., 2006) respectively. 2.2. Force and power calculation Since the grid is not conformal to the body, the force coefficients on the body are calculated differently compared to that of a body conforming grid, by using the forcing term fc obtained earlier. In this case, the total force, in non-dimensional terms, based on Kim et al. (2001), is calculated as   Z Z ∂ui ∂ui uj þ Fg ¼ f cg dVol þ dVol; ð10Þ ∂t ∂xj solid solid where Vol is the volume of the wing. All the grid cells' forces within the volume of the wing are summed up to give the total force (Fg). W.B. Tay et al. / Journal of Fluids and Structures 55 (2015) 237–261 241 Table 2 Force coefficient comparison between experimental and numerical results. f Average cd (Calderon) Average cd (IBM) Average cl (Calderon) Average cl (IBM) 0.45 0.44 0.39 ( 1.36 1.25 ( 11%) 7%) Fig. 2. PIV and IBM phase-averaged vorticity contour plots comparison at f=0.45. PIV results extracted from Calderon et al. (2010). The force coefficients ct and cl are defined as cg ¼ Fg 1=2ρv2 S ð11Þ where g, ρ, ν and S represents lift (l) or thrust (t), density of the fluid, reference velocity and the wing planform area respectively For simplicity, the force coefficients ct and cl are assumed to be equal to the forces Fg. In other words, we let 1=2ρv2 S ¼ 1. This does not pose any problem because we only compare between results within this study. The (nondimensional) power input for the flapping is calculated based on the angular velocity of the wing and the torque it produces. Torque (Tr) is given by  Z ð12Þ Tr ¼ ðr  F ÞdVol;  ~ ~ where r and F are the distance of the grid cell from the rotating axis and the grid cells' forces at that cell respectively. All the ~ ~ grid cells' torque within the volume of the wing are summed up to give the total torque (Tr).  Power input is then given by P in ¼ Tr Uω;  ~ ð13Þ where ω is the angular velocity of the wing. ~ Propulsive efficiency η, defined as a measure of the energy lost in the wake versus energy used in creating the necessary thrust, is given by η¼ ct ; P in ð14Þ where ct is the output thrust coefficient. The overbar indicates averaged values. 2.3. Solver validation The IBM solver has been validated against a 3D PIV plunging wing experiment at a Re of 10 000 (Calderon et al., 2010), as described in the DF2 simulation paper by Tay et al. (2014). One of the cases (f ¼0.45) is shown here again for completeness. The wing, with a semi-aspect ratio of 2, was placed in a water tunnel and plunged with an amplitude of 0.15c. It was held at a constant angle of attack of 201 with a non-dimensional frequency of 0.45. The grid used is 177  368  734, with minimum grid spacings of 0.018, 0.006 and 0.006 in the x, y and z directions respectively. The minimum grid spacing used in the x direction is lower in resolution because the velocity along the span in the x direction is much lower compared to the y and z directions. The simulation results (force coefficients, Table 2 and vorticity contour plots, Fig. 2) agree reasonably with the experimental results. A further validation presented here is the comparison against the experiment by Lua et al. (2014) involving the simultaneous sweeping and pitching motion of a hawkmoth-like wing in a water tunnel at a Re of 10 000, based on the 242 W.B. Tay et al. / Journal of Fluids and Structures 55 (2015) 237–261 Fig. 3. Force coefficient comparison between experimental and numerical results. (The paper by Lua et al. (2014) only contains the lift coefficient result. The drag coefficient result was obtained from Lua through personal communication.) Fig. 4. (a) Normal, (b) close-up isometric and (c) front view of the 3D Cartesian grid with the biplane wings in red. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.) average wing velocity at the location of the second moment of wing area. The wing motion is simple harmonic under quiescent flow. An in-house-developed waterproof force transducer measured the forces acting on the wing. The experiment was simulated using the IBM solver at a minimum grid spacing of 0.018c, which corresponds to a grid resolution of 321  298  574. Due to the large aspect ratio of the hawkmoth wing, the computational requirement is high and hence higher resolutions were not attempted. Nevertheless, the results between the experiment and simulation at the current resolution (Fig. 3) agree reasonably, with the peak lift coefficient of the simulation slightly lower than that of the experiment by 10%. The validations show that the IBM solver is able to accurately simulate moving wings up to Re¼10 000. 2.4. Simulation setup This section describes the setup of the simulations in the present study. Fig. 4 shows the 3D Cartesian grid with the biplane wings. The orientation of the wings is such that the incoming flow is along the positive z axis direction in forward flight condition. The wings are mirrored along the yz plane such that a symmetry boundary condition at x¼ 0 is used to reduce computational cost. A small portion of the wing (0.1c) is shifted inwards along the x direction to enable the IBM to capture the shape of the wing fully. If the root of the wing is exactly at the yz plane, the full geometry of the wing will not be captured. Instead, some parts of the wing near the root will be missing. Hence, the length of the wing in the simulation has a span of 1.85c. The extra 0.1c portion of the wing is purely fictitious and does not participate in the simulation, as shown in Fig. 4(c). The total size of the computational domain in chord lengths is 8  16  20 in the x, y and z-directions, respectively, with a minimum grid length of 0.009c; refinement is used in the region near the wings. The resultant total number of cells for the domain is 266  530  405, based on the grid convergence study performed by Tay et al. (2014). The average time taken for one simulation cycle using 96 AMD Opteron 6234 processors running in parallel (MPI) is 30 h. 2.5. Grid convergence study Grid convergence study was performed in Tay et al. (2014) based on the DF2 configuration. The convergence study is repeated for the DFM configuration using similar minimum grid lengths (dx) of 0.018, 0.009 and 0.006. Of the three force coefficients (thrust, lift, side), only the thrust coefficient is compared because the grid requirement for the lift and side forces is much lower than that of the thrust. Besides comparing the thrust coefficient, we also compare the pressure contour plots at different resolutions to ensure convergence. Figs. 5 and 6 shows that a minimum grid length of 0.009 gives sufficient W.B. Tay et al. / Journal of Fluids and Structures 55 (2015) 237–261 243 Fig. 5. ct of the lower wing of the DFM at dx¼0.006, 0.009 and 0.018. Fig. 6. Pressure contour plot at time ¼ 0.76 T for dx: (a) 0.018, (b) 0.009 and (c) 0.006. Green circles denote differences. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.). accuracy and is computationally cheaper than the minimum grid length of 0.006. Hence, dx¼0.009 will be used for all subsequent simulations in this paper. 3. Research methodology The main objective of this study is to gain a deeper understanding of the aerodynamics of the DFM biplane flapping configuration, also in comparison to the DF2 configuration, from which it distinguishes itself mainly by the larger stroke angle and associated additional clap-and-fling moment in the stroke cycle. The particular focal points of the study have been stated in the introduction. This section will first describe the geometrical simplifications that were introduced in the simulations, followed by the details on how we perform the simulations to achieve the objectives. 3.1. Model simplifications Due to the complexity in the design of the DFM and its wings, we have introduced some simplifications in the simulations. The actual DFM wings consist of a stiff carbon rod as leading edge to which a thin membrane Mylar sheet is attached supported with thin carbon rod stiffeners. This flexible wing construction deforms during the flapping cycle under the influence of the fluid and structure interaction (FSI). Since the current solver does not have FSI capability, the wings are represented by thin plates (thickness ¼0.06c) that are either treated as rigid surfaces, or that are given a prescribed deformation, the latter as an approach to study the effects of wing flexibility. To simulate the chordwise pitching of the wings due to FSI, a sinusoidally varying pitching of α0 ¼ 101 is added to the wings during flapping. Moreover, the design of the biplane wings of the DFM is such that the roots of the upper and lower wings are very close to one another. However, due to the chordwise pitching and thickness of the wing, we add a small amount of clearance (equivalent to 0.3c) between the two wings to prevent the wings from intersecting. An indication that the simplified model is a reasonable approximation of the actual DFM is to compare the lift generated by both of them. The DFM typically flies with a non-zero body inclination angle αb (approximately from 101 to 301, based on captured videos2), as shown in Fig. 7, to vector the force generation by the wings into lift and thrust. As will be explained later, this correspond to rotating the body reference frame XbYbZb along the x axis by an angle of αb. The average total cl from 2 The captured videos can be found on http://www.youtube.com by using the keyword “Delfly micro”. 244 W.B. Tay et al. / Journal of Fluids and Structures 55 (2015) 237–261 αb ¼01 to 201 is obtained from the body inclination angle comparison results in Table 5. We choose αb ¼151 for the lift approximation because it can give the highest lift with a very small amount of drag (0.02). The total cl generated by the right hand side of DFM-151's wings is 1.63. Hence the total cl produced by both sides of the wings is 3.26. The actual lift force is given by Lif t ¼ cl ðρU 2ref c2 Þ ¼ ð3:26Þð1:18Þð3:0Þ2 ð0:0289Þ2 ¼ 0:029N ¼ 3:0g; ð15Þ where ρ, Uref and c are the density of air (1.18 kg/m3), reference forward velocity (3.0 m/s) and chord length (0.0289 m) respectively. The current weight of the DFM is approximately 3 g, equivalent to 0.029 N. Hence, under the current simulation, the lift generated according to the numerical simulations is equivalent to the force required for sustained flight. In actual flight conditions, due to parasitic drag, the total drag experienced by the DFM will be larger, indicating that αb must be reduced to increase its thrust to overcome the drag. In correspondence, the lift will be lower than 0.029 N. However, the current configuration uses rigid thin plate as the wing instead of flexible membrane, as in the actual DFM. We will show in following sections that the addition of flexibility can improve the performance of the DFM and provide higher lift. Nevertheless, this current calculation serves as a useful gauge regarding the generation of lift of DFM in the forward flight condition. Based on this validation, it is expected that the simulation is capable of capturing the essential flow behavior and that therefore valuable insights can still be obtained from the simulation results. 3.2. Simulation parameters This section discusses the parameters used in the simulations. Fig. 8 shows the planform, stroke angle and pitch angle of the DFM wing. It has the same shape as that of DF2 but its stroke angle is larger than that of the DF2, as given in Table 1. For the model wings, the flapping stroke (θ) and pitch (α) angles are given respectively by θ ¼ θ0 sin ð2πf t ϕÞ þθmid ; α ¼ α0 sin ð2πf tÞ þ αb ; ð16Þ ð17Þ where θmid ; θ0 ; θ; αb ; α0 ; α; f ; ϕ are the midpoint stroke angle, stroke amplitude, stroke angle, body inclination angle, pitch amplitude, pitch angle, frequency and phase angle respectively. The pitch amplitude and the center of rotation of pitching are at 101 and 0.5c respectively. The phase angle ϕ is fixed at 901 in all simulations. With reference to Fig. 9, XYZ, XbYbZb and xwywzw refers to the global, body and wing fixed axes respectively. The flapping stroke (θ) angle is attached to the wing xwywzw axis and it makes an angle θ with respect to the z axis of the XbYbZb coordinate system. The wing pitches with an angle α with respect to the x axis of the wing xwywzw axis. Hence, the mean value of the pitch angle is equal to the body inclination angle αb described earlier. The Reynolds number is determined using: Re ¼ U ref c ; ν ð18Þ where Uref, c, ν are the reference velocity, chord length and kinematic viscosity. The forward flight velocity or the mean wing tip velocity is used for Uref in the forward flight conditions, respectively. With a forward flight speed of 3 m/s and chord length of 0.0289 m, the Re for DFM is approximately 5500 in forward flight condition. The non-dimensional frequency f is given as f¼ f actual c ; U ref ð19Þ where factual is the actual flapping frequency. With the actual flapping frequency at 37.5 Hz, the non-dimensional frequency f is calculated to be 0.36. In the first objective, a comparison is made between the performance of the DFM and DF2 flapping configurations. To isolate pure size effects between the two MAVs, all parameters in the simulations are taken the same for the two configurations except for the stroke and midpoint stroke angles, which are given in Table 3. This will allow a fair comparison to determine if the stroke and dihedral angles of the DFM are indeed more advantageous under otherwise equal conditions. In the second objective involving body inclination comparison, the DFM is rotated along the x axis in the global reference frame to give body inclination angle αb from 101 to 201. 3.3. Prescribed spanwise and chordwise flexibility As mentioned earlier in the introduction, prescribing spanwise and chordwise wing deformation will allow us to study and compare the effect of each flexibility with respect to the rigid wings. The standard DFM configuration is used. Hence, αb is 01 and the wing angles are as given in Fig. 8. The deformation of the wings is modeled based on slow motion video capture of the actual DFM wings, as shown in Fig. 10(a). Looking at the wing deformation, there is almost no spanwise deformation W.B. Tay et al. / Journal of Fluids and Structures 55 (2015) 237–261 245 Fig. 7. Explanation of the body inclination angle. Fig. 8. (a) Planform. (Only the actual planform of the wing is shown. The additional 0.1c portion of the wing explained in the simulation setup is not shown.) (b) Stroke angle of the DFM wing (front view). (c) Pitch angle of the DFM wing (side view). The thick solid lines represent the possible positions of the wings. present. Nevertheless, we impose a quadratic variation along the span to investigate the effect of spanwise deformation on the performance of the DFM. So in this study, the upper and lower wings undergo a sinusoidally varying spanwise deformation with a maximum spanwise deformed length (dlsw-max) from 0.15c to 0.45c. A maximum value of 0.45c is chosen because the wings will intersect one another when dlsw-max 40.45c. The spanwise flexibility distribution, with reference to Fig. 11(a), is given by dlsw ¼ dlsw max dcsw  dsemispan 2 sin ð2πf tÞ ð20Þ where dlsw, dlsw-max, dcsw, dsemispan are the perpendicular spanwise deformed length magnitude, maximum spanwise deformed length, distance from point on wing to the root along spanwise direction and semi wingspan of the wing respectively. In Eq. (20), dlsw increases from the root to the tip of the wing, as shown in Fig. 11(a). The deformation dlsw always acts opposite to the direction of flapping. The f used is the same as the non-dimensional flapping frequency. On the other hand, in Fig. 10, one can see the chord deformation is increasing along the span from the root to approximately its mid-section (red dotting line), after which the deformation becomes almost uniform. For simplicity, we take the distance from root to wing's mid-section droot_to_mid to be 0.5c. Based on this information, we modeled the chord deformation equation, with referenced to Fig. 11(b), as For dccw odistance from root to wing's mid-section droot_to_mid: dlcw ¼ dlcw max  1 4:0dccw c 0:5c 2  dcsw droot_to_mid 2 sin ð2πf tÞ ð21Þ 246 W.B. Tay et al. / Journal of Fluids and Structures 55 (2015) 237–261 Fig. 9. Isometric view of one of the wings with the enclosed yellow region as the body. Second wing remove for clarity. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.) Table 3 Wing angles of the upper and lower wings of the DFM and DF2 in degrees. Wing angles Upper Lower Upper Lower (DFM) (DFM) (DF2) (DF2) θmid (deg) 53.5 3.5 30.0 8.0 θo (deg) 25.5 25.5 22.0 22.0 For dccw Zdistance from root to wing's mid-section:  2 4:0dccw 0:5c dlcw ¼ dlcw max 1 sin ð2πf tÞ c ð22Þ where dlcw, dlcw-max, dccw, c are the perpendicular chordwise deformed length magnitude, maximum chordwise deformed length, distance from point on wing to the leading edge along chordwise direction and chord length of the wing respectively. Eqs. (21) and (22) have been formulated such that dlcw has its maximum length at the wing's mid chord location, as shown in Fig. 11(b). This is based on the slow motion DFM wing capture and close visual observation of the DFM's wing.3 Similarly, dlcw always acts opposite to the direction of flapping and dlcw-max varies from 0.15c to a maximum of 0.35c, to prevent wing intersection. There are still some minor discrepancies between the actual and modeled wings due to model simplifications mentioned in Section 3.1 and the limitation of the above deformation formulations. However, it should be emphasized that the purpose of this deformation study is not to simulate an almost exact DFM spanwise/chordwise wing deformation using empirical formulation. Rather, we sought to obtain a modeled wing which gives a reasonable deformation approximation, and from there study the deformation effects. From the spanwise and chordwise results, an additional simulation involving both spanwise and chordwise flexibility (superposition of both) will be included to assess the effect of combining both forms of flexibility. 3.4. Different wing simulation configurations In order to achieve the objectives mentioned in the introduction, we simulate the wings in different configurations, as given in Table 4. We will first start with the simulation of DFM and DF2 without body inclination or wing deformation in forward flight. It allows us to compare the performance between the two configurations. Next, we move on to the body 3 The deformed shape was given by S. Deng (who performed close visual observation) through personal communication. W.B. Tay et al. / Journal of Fluids and Structures 55 (2015) 237–261 247 Fig. 10. Front view of the actual video (left) and modeled (right) DFM right wings at time: (a) 0.28 T and (b) 0.84 T. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.) Fig. 11. (a) Front view of wing undergoing spanwise deformation and (b) isometric view of wing undergoing chordwise deformation. inclination test by simulating DFM at varying αb of 101, 151 and 201. This is followed by the DFM spanwise and chordwise deformation investigation, whereby dlsw-max and dlcw-max varies from 0.15c to 0.45c and 0.15c to 0.35c respectively. Lastly, a case which combines spanwise dlsw-max ¼0.25c and chordwise deformation dlcw-max ¼0.25c is run. 4. Results and discussions The simulations were run for five flapping periods. Results show that the forces become almost periodical after the first cycle. The aerodynamic performance is analyzed in terms of efficiency, thrust, and lift coefficients, while the flow behavior is visualized with the help of Q criterion (Q) (Hunt et al., 1988) iso-surface contour plots. The results and discussions will be presented in the following sections, addressing the four focal objectives of the study. 248 W.B. Tay et al. / Journal of Fluids and Structures 55 (2015) 237–261 Table 4 Different wing simulation configurations. Configuration type Simulation description DFM DF2 DFM-10/15/201 DFM-s_df-0.15/0.30/0.45 DFM-c_df-0.15/0.25/0.35 DFM-s/c_df-0.25 Delfly-Micro without body inclination or wing deformation in forward flight Delfly II, same as DFM except for the stroke angles DFM with body inclination angle αb ¼101/151/201 DFM with spanwise deformation dlsw-max ¼ 0.15c/0.30c/0.45c DFM with chordwise deformation dlcw-max ¼0.15c/0.25c/0.35c DFM with spanwise dlsw-max ¼0.25c and chordwise deformation dlcw-max ¼0.25c Fig. 12. ct and cl of the DFM and DF2 upper and lower wing configurations over one period (1 T). The total cl of the DFM and DF2 are also shown. Fig. 13. Mz and Pin of the DFM and DF2 upper and lower wing configurations over one period. 4.1. DFM and DF2 comparison Figs. 12 and 13 show the variation in the ct/cl and Mz (moment about the z axis)/Pin (power input) over one period for the DFM and DF2 upper and lower wing configurations respectively. There are some spikes in these figures, which happens at times when IBM is applied to moving body problems like in the current scenarios (Lee et al., 2011). The two reasons for the spikes are pressure spatial discontinuity when a grid point in the solid body changes to a fluid, and velocity temporal discontinuity when a grid point in the fluid changes to that of a solid body. As shown in both grid convergence study (current paper and the paper by Tay et al. (2014)), increasing the grid resolution reduces the number of spikes. The current results only have a small number of spikes, and more importantly, these spikes do not influence the assessment of the results. Hence, to reduce computational cost, the minimum grid length of 0.009 is used. Fig. 14 shows a visualization of the vortex structures of DFM and DF2 by means of Q-criterion iso-surfaces with Q¼250 at three instants in time.4 The labeled vertical lines in Figs. 12 and 13 represent the same time instants as in Fig. 14. Due to the higher stroke angle for DFM, the ct for both its upper and lower wings are higher than that of DF2. Another difference lies in the twin peaks of each configuration. In Fig. 12, the ct of DF2 for both upper and lower wings at position a (corresponding to the mid-outstroke) is 4 Note that the lateral viewing direction for the last instant is opposite to that of the previous two instants, in order to properly visualize the vortex formation during the instroke. W.B. Tay et al. / Journal of Fluids and Structures 55 (2015) 237–261 249 Fig. 14. Vortex structures of DFM and DF2 at Q criterion (Q) ¼250 at three instants in time. The double black arrow represents in the inflow direction. (a) and (b) represent the front isometric view while (c) represents the back isometric view. Green circles denote differences. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.) lower than position c (mid-instroke, where the upper and lower wings clap together). This difference in peak height is also observed in the DF2 simulation at f¼ 0.6 in Tay et al. (2014) at Re¼1000 and 5000. The reason for this is that when transiting from instroke to outstroke, the upper and lower wings are in close proximity to one another, resulting in the clap and fling mechanism (Weis-Fogh, 1973) as described in the introduction. This produces higher ct, also evident from the high pressure in the region between the upper and lower wings during the mid-instroke, as shown in Fig. 15(a) at t¼0.32 T. However, when transiting from outstroke to instroke, there is no clap and fling mechanism involved in the DF2 simulation, and hence ct and the pressure in the region between the right and left (mirrored) upper wings are lower. On the other hand, in the DFM case, the higher stroke angle causes the upper and lower wings to be in close proximity to one another during the instroke and outstroke. Hence, the clap and fling mechanism is present and this produces almost equal ct peaks during both instants. As shown in Fig. 15(b) at t ¼0.32 T and 0.80 T, the pressure in the region between the upper and lower wings and in the region between the left and right upper wings are similar in magnitude. On the other hand, Nakata et al. (2011), who investigated a comparable MAV configuration with approximately the same stroke angles, obtained uneven peaks for both the thrust and lift output. In their simulation, the outstroke generates higher ct than the instroke for both the upper and lower wings. The reasons for this difference could be due to: 1. The motion and flexibility of the wings used in the simulation of Nakata et al. (2011) had been reconstructed based on images captured using high-speed cameras. The current study uses rigid wings. 2. The Re and f of the current study's and Nakata et al.'s simulation is 5500/0.36 and 3500/0.62 respectively. 250 W.B. Tay et al. / Journal of Fluids and Structures 55 (2015) 237–261 Fig. 15. Surface pressure of wings and pressure at a distance z¼ 0.25 slice from the wings' leading edge for the (a) DF2 and (b) DFM. Green circles denote differences. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.) Fig. 16. Total ct and cl of the DFM at αb ¼ 01, 101, 151, 201 over one period. In the lift comparison, the main difference between DFM and DF2 lies in their peak magnitude of their respective upper and lower wings. The peak magnitude of DFM's lower wing is higher than its upper wing while in DF2, the peak magnitudes of the upper and lower wings are approximately the same value. To investigate this difference, we simulate only the DFM's individual upper or lower wing (not shown) separately. We find that the peak magnitude of the lower wing is also higher than its upper wing, albeit with a smaller magnitude due to the absence of the clap and fling motion. Hence, the reason for the unequal lift peak in the DFM's case is due to the different θmid of the upper and lower wings, which is 53.51 and 3.51 respectively. One disadvantage of the unequal peak is the increased rocking amplitude during flight, as shown in the larger total lift peak of DFM, as compared to DF2's in Fig. 12. The rocking may affect the mounted camera's image clarity, although the amplitude of the total cl is still approximately three times smaller than that of a single DFM upper wing. If necessary, one can make use of suitable image stabilization algorithm to improve the image quality. Because of the higher stroke angles of the DFM, its peak magnitudes of the Mz and Pin (2.9 ( þ26%), 2.9 (þ 38%)) are higher than that of the DF2 (2.3, 2.1). Despite the increase, the total efficiency η, obtained using Eq. (14), gives approximately the same value of 0.11 for both cases. Considering the larger ct produced by the DFM (0.30 ( þ58%)), the DFM flapping configuration is in general more superior to the DF2's (0.19) under the current Re and flapping frequency. However, one practical disadvantage is the requirement of a more powerful motor for the actual DFM FMAV due to the higher power input. Fig. 14 visualizes the vortex structures for both DFM and DF2. The relatively high Re of 5500 results in the breakup of the vortex patterns into small and fine vortex structures, which complicates the analysis of the vortical structure evolution during flapping. The spiral conical leading edge vortex (LEV) can be observed in both cases along the wingspan, similar to the DF2 simulation by Tay et al. (2014) at Re¼1000 and 5000. The higher stroke angles of the DFM impart more energy into the flow compared to DF2 and hence the size and density of the vortical structures are also larger. The observation of the trailing edge vortex (TEV) in Fig. 14c for the DFM case (which is not observed in DF2's) is most likely due to the same reason. The TEV in DF2 can be observed at low Q criterion values (Qo50, not shown), indicated that its TEV is much weaker compared to that of DFM. W.B. Tay et al. / Journal of Fluids and Structures 55 (2015) 237–261 251 Table 5 Average total ct, cl and η of DFM/-10/-15/-20o. Configuration type Average total ct Average total cl Total η DFM DFM-101 DFM-151 DFM-201 0.30 0.14 0.02 0.31 0.03 1.11 1.63 2.05 0.11 0.05 Not applicable Not applicable Fig. 17. (a) Rotation of the DFM configuration by 201, resolving the original ct/cl into new components and (b) the ct of the lower wing of DFM and DFM-201. Fig. 18. Change from conical to uniform LEV for DFM at αb: (a) 01, (b) 101 and (c) 201. Green circles denote differences. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.) 4.2. Body inclination angle comparison The body inclination angle αb of the DFM in forward flight has been varied from 01 to 201. Besides comparing the vortical structure, we also calculate the lift force in the earlier section to prove that the current simplified DFM model is a good approximation of the actual DFM. Fig. 16 shows the total ct and cl of the DFM at αb ¼01, 101, 151, 201 over one period while Table 5 shows the average total ct, cl and η. Since the current focus is in the change in total force of the upper and lower wings, we do not show the individual forces generated by individual wings. As αb increases from 01 to 201, ct changes its shape from equal to unequal twin peaks, and ct from 0.30 to 0.31. In other words, at αb ¼201, the DFM is incapable of providing sufficient thrust to sustain forward flight. In the current simulation, the maximum αb for forward flight should lie around 151, where the ct is close to zero. Note that the thrust as defined here is the net forward aerodynamic force opposite to the flight direction and, hence, takes the drag force into account. The change in the ct 's shape is due to part of the thrust being vectored to become lift. By rotating the ct and cl of the αb ¼01 case by the corresponding αb and resolving them into components along the vertical and horizontal axes one can see that the resultant graph's shape becomes very similar to that of the αb case, as shown in Fig. 17. The graph of cl shift up as αb increases, producing as much as cl ¼ 2:05 at αb ¼201, albeit 252 W.B. Tay et al. / Journal of Fluids and Structures 55 (2015) 237–261 Fig. 19. Total ct and cl of the DFM from dlsw-max ¼ 0 to 0.45c over one period. Fig. 20. Vortex structures of DFM and DFM-s_df-0.45 at Q ¼250 at three instants in time. Fig. 14(a) of DFM is shown here again for ease of comparison. Green and blue circles denote differences. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) drag is also produced. Besides the vectoring of forces, another reason for the increase in lift is due to the change from the conical to the stronger uniform LEV, as shown in Fig. 18. Unlike the conical LEV, whereby its cross-sectional area decreases as it gets closer to the wing root, the uniform LEV has the same cross sectional area of the vortex along its span. The larger size of the uniform LEV indicates stronger vortex strength. Besides LEV size, the strength of the vortex also depends on the vorticity in its core. However, in this case, this should be similar between DFM at different body inclination angle because the flapping stroke, pitch angles and frequency are the same. The key difference is the body inclination angle, which cause the LEV to remain attach to the wing, thereby increasing the lift. The uniform LEV had also been observed in the simulation of the DF2 in forward flight at αb ¼301 (Tay et al., 2014). The Re used in the simulation was 1000 and the uniform LEV was in W.B. Tay et al. / Journal of Fluids and Structures 55 (2015) 237–261 253 Table 6 Average total ct, cl, Pin and η of DFM,DFM-s_df-0.15/0.30/0.45. Configuration type Average total ct DFM DFM-s_df-0.15 DFM-s_df-0.30 DFM-s_df-0.45 0.30 0.36 0.39 0.45 Average total cl 0.03 0.00 0.00 0.07 Average total Pin Total η 2.68 3.10 3.40 3.83 0.11 0.12 0.11 0.12 Fig. 21. Pressure iso-surface of (a) DFM and (b) DFM-s_df-0.45 at radius¼ 1.0c at different time instants. The surface pressure of the wings is also shown. Green circles denote differences. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.) comparison larger than the current simulation, which was due to the higher αb and lower Re. The appearance of the uniform LEV in these two simulations despite having different flapping parameters show that it may not be a unique feature which happens only during very specific flapping motion. It will also be interesting to perform simulations in single wing configuration to check if uniform LEV is present at non-zero inclination angle. 4.3. Spanwise flexibility comparison In this section, the flexibility along the span of the wing varies from zero to a maximum wing tip deflection of 0.45c. Figs. 19 and 20 show the total ct and cl of the DFM from dlsw-max ¼0 to 0.45c over one period and the vortex structures of DFM and DFM-s_df-0.45 at three instants in time respectively. We only show the outstroke duration because at this time interval, the difference in the LEV is most obvious under the current isometric view. From the force plots in Fig. 19, there is an obvious phase shift in both the ct and cl plots as dlsw-max increases. The shifting or phase delay is due to the spanwise wing deformation, which causes the same point on the wing to reach its new position at a later time, compared to a 254 W.B. Tay et al. / Journal of Fluids and Structures 55 (2015) 237–261 Fig. 22. Pressure contours of (a) DFM and (b) DFM-s_df-0.45, obtained from plane projection of previous figure. Green circles denote differences. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.) non-deformed wing. The difference in position between the deformed and non-deformed case is greatest during the midstroke. The resultant velocity of the point is smaller because the deformation acts in the opposite direction to the wing's sinusoidal motion. Similar phase shifts have been observed in the spanwise deformation studies by Tay et al. (2014) and Yang et al. (2012). Fig. 20 shows that the conical LEVs (circled green) on the DFM-s_df-0.45's wings stay much strongly attached compared to that of the DFM, which contributes to the larger thrust, as shown in Table 6. In other words, the LEV of DFM-s_df-0.45 stays attached to the leading edge wing surface for a much longer time along the span. In comparison, the LEV of DFM detaches and dissipates away, especially in the region near the wing root. The reason for the higher thrust can be understood by observing the pressure iso-surface at a radius of 1.0c from the location mid-way between the axis of rotation of the upper and lower wings in Fig. 21. For ease of analysis, the pressure iso-surface is “un-wrapped” and projected onto a plane in Fig. 22. The time instants shown correspond to the instants of maximum thrust. The pressure on the pressure and suction side of the wings of the DFM-s_df-0.45 are higher and lower than that of the DFM respectively, therefore giving more thrust. There is a time difference of approximately 0.04 T between the DFM and DFM-s_df-0.45 due to the phase shift explained earlier. This strong LEV attachment mentioned earlier was not observed in the DF2 spanwise deformation study (Tay et al., 2014), although the DF2 and the current DFM have the biplane wing configuration and both undergo similar spanwise deformation. Their key differences lie in the flapping kinematics and flapping frequency. It was originally postulated by Ellington et al. (1996) that the reason for the attachment of the LEV is due to the positive spanwise velocity removing vorticity from the conical LEV. This allows the LEV circulation to remain sufficiently small for attachment. Another difference W.B. Tay et al. / Journal of Fluids and Structures 55 (2015) 237–261 Fig. 23. Total ct and cl of the DFM from dlcw-max ¼ 0 to 0.35c over one period. Fig. 24. Vortex structures at Q ¼250 for time ¼0.2 T, 0.4 T and 0.76 T of (a) DFM, (b) DFM-c_df ¼0.25 and (c) DFM-c_df¼ 0.35. 255 256 W.B. Tay et al. / Journal of Fluids and Structures 55 (2015) 237–261 Fig. 25. z slice at distance 0.1c from the wing's leading edge, showing strength of axial velocity along spanwise direction for (a) DFM and (b) DFM-c_df¼ 0.25. between the DFM and DFM-s_df-0.45 is the stronger tip vortices (circled blue) for the latter configuration. This observation coincides with the results by Birch and Dickinson (2001), who used a dynamically scaled robotic model of a Drosophila wing flapping in mineral oil to show that it is the downward flow induced by tip vortices which limits the growth of the LEV and improves its stability. Similarly, Beem et al. (2012) performed PIV experiments to investigate the effect of spanwise flow on LEV stabilization. It was found that no LEV stabilization is present for varying degrees of spanwise flow tested. The author mentions that it is possible that LEV stabilization is due to a combination of spanwise flow, three-dimensional kinematics, rotation, and/or flexibility. Hence, in this case, the unique combination of spanwise flexibility and the DFM stroke angles show improvement in the LEV stabilization. The spanwise deformation introduces additional constraints which also increases the input power. However, due to the corresponding increase in average thrust, the η remains almost the same at increasing flexibility. 4.4. Chordwise flexibility comparison In this section, the wings of the DFM undergo chordwise deformation. As mentioned earlier in the research methodology section, the wings deform along the wing chord. At the same time, there are variations in the chordwise deformation along the wingspan. This is unlike the spanwise deformation case earlier, whereby there is uniform spanwise deformation along the wing chord. Fig. 23 shows the variation of the total ct and cl of the DFM from dlcw-max ¼ 0 to 0.35c over one period while Fig. 24 shows the vortex structures at Q¼250 for time ¼0.2 T, 0.4 T and 0.76 T of the DFM from dlcw-max ¼0 to 0.35c. The DFM-c_df ¼0.15 case is not shown in Fig. 24 because it is similar to that of DFM-c_df¼0.25. The phase shift of the ct and cl plots also happens, similar to that of the spanwise deformation. There are three differences between the deformed and non-deformed cases in Fig. 24. Firstly, the conical LEV is now much weaker (appearing for a short duration and thereafter detaches from the leading edge) and strong tip vortex (TV) emerges in the deformed case. At time ¼0.76 T, the TEV is also stronger in the deformed case. This is the result of the difference in the wing twist angle at the tip, due to the varying degrees of chordwise deformation along the span for the deformed case. Next, there is more shed vortices as dlcw-max increases, resulting in a more “messy” wake region. Hence, unlike the spanwise deformation case earlier, the strong TV does not accompany the strong LEV. It has been argued that the stability of the LEV is due to strong axial velocity along the spanwise direction (Aono et al., 2008; Ellington, 1999). Here, we look at the z slice at 0.1c from the wing's leading edge in Fig. 25 for DFM and DFM-c_df¼0.25. With reference to Fig. 24 at time¼0.2 T, the strong LEV of DFM is accompanied by strong axial velocity. On the other hand, the axial flow speed is much lower for the DFM-c_df¼0.25. This shows that in this case, the stability of the LEV directly relates to the axial velocity. W.B. Tay et al. / Journal of Fluids and Structures 55 (2015) 237–261 257 Table 7 Average total ct, cl and η of DFM,DFM-c_df-0.15–0.35. Configuration type Average total ct Peak |cl| Average total Pin Total η DFM DFM-c_df-0.15 DFM-c_df-0.25 DFM-c_df-0.35 0.31 0.67 0.74 0.56 2.85 3.04 3.54 3.87 2.75 3.32 3.64 4.03 0.11 0.20 0.20 0.14 From the average total ct, cl and η of DFM,DFM-c_df-0.15/0.25/0.35 in Table 7, the current preferred configuration type is DFMc_df-0.25, which give both maximum ct and η. Increasing dlcw-max beyond 0.25c to 0.35c instead reduces the ct and η, although it is still higher than the non-deformed case. S. Deng et al. (2013) also show in their DFM experiments that high flexibility reduces the thrust produced by the wings. To explain the variation in ct , we look at the pressure contours of the chordwise deformed cases in Fig. 26, obtained from the contours of pressure on a cylindrical surface of radius¼1.5c. We choose a larger radius of 1.5c (instead of 1.0c as in Fig. 22) because the difference in chordwise flexibility increases as the spanwise distance increases. At time¼0.24 T and 0.80 T, we can observe the large pressure difference between the pressure and suction side of the wings. As dlcw-max increases, the pitching angle of each section of the wing along the spanwise direction increases due to the chordwise deformation. These two reasons together with the vectoring of the forces therefore increases the thrust. However, at time ¼0.16 T, we can observe the TEV appearing for both the dlcw-max ¼0.25c and 0.35c. This creates low pressure around the rear part of the wing, reducing the thrust. In fact, if we refer to the force graphs in Fig. 23, time¼0.16 T (dotted line “a”) corresponds to drag production for dlcw-max ¼0.35c. The same event (appearance of TEV at the rear) happens at around time ¼0.62 T as well, although it is not shown here. The TEV is much stronger in the dlcw-max ¼0.35c case and hence the resultant ct at dlcw-max ¼0.35c is smaller than that of dlcw-max ¼0.25c. On the other hand, η increases as dlcw-max increases, until a maximum of 0.21, which represents an increase of 90% over the original non-deformed case. At DFM-c_df-0.35, η drops to 0.14. η depends on both thrust and power input, as given in Eq. (14). ct increases and then decreases as dlcw-max increases. However, P in increases as dlcw-max increases. As a result, η drops. One can also look from another perspective through the vortical shedding in Fig. 24, whereby the DFM-c_df-0.25 case has cleaner, more compact and coherent vortex shedding as compared to that of DFM-c_df-0.35. Less energy is wasted in creating the vortices and hence η increases (Pedro et al., 2003). 4.5. Simultaneous spanwise and chordwise flexibility Results from the previous two sections show that increasing the spanwise flexibility increases the thrust. On the other hand, thrust and efficiency first increases and then decreases as the chordwise flexibility increases. Based on these results, the preferred dlsw-max and dlcw-max are 0.45 and 0.25 respectively. However, combining both dlsw-max ¼0.45 and dlcw-max ¼0.25 causes the wings to intersect each other during flapping. Hence, a simulation based on both dlsw-max and dlcw-max ¼0.25 is attempted. Fig. 27 shows the total ct and cl of the DFM with dlsw-max ¼0.30c, dlcw-max ¼0.25c and dls/cw-max ¼ 0.25c over one period. The labeled vertical dotted lines a and b represent the time instants when the instantaneous thrust are at their maximum and minimum respectively. The contours of pressure on a cylindrical surface of radius¼ 1.5c for DFM-s_df¼0.30, DFM-c_df-0.25 and DFM-s/c_df-0.25 at these two time instants are given in Fig. 28. The Q criterion iso-surface contour plots in Fig. 29 shows that the DFM-s/c_df-0.25 has both relatively strong LEV and TV, essentially a combination of DFM-s_df¼0.30 and DFM-c_df-0.25. For DFM-s/c_df-0.25, the appearance of the LEV at time¼0.36 T (line “a”) further increases the pressure difference between the pressure and suction side of the wing. This, together with the high pitching angle due to the chordwise deformation, generates the highest thrust among the four cases compared, as shown in Table 8. The proportional increase of P in for the DFM-s/c_df-0.25 helps to maintain the η at 0.20, similar to that of the DFM-c_df-0.25. The combined spanwise/chordwise flexibility indeed gives better performance compared to either only spanwise or chordwise flexibility of the same degree. However, similar to the chordwise deformation study, appearance of TEV reduces the thrust. In another numerical study on tethered desert locust (Schistocerca gregariaby) by Young et al. (2009) at Re¼5000, it was found that that time-varying wing twist (due to combined spanwise and chordwise flexibility) and camber were important for the maintenance of attached flow. Hence, in their study, wing twisting also improve the LEV attachment. In Section 3.3, we mention that the slow motion video capture of the actual DFM wings (Fig. 10) showed almost no spanwise deformation. From the simultaneous spanwise and chordwise deformation results, we see that adding spanwise deformation further improves the wing's performance. Hence, in the design of the actual DFM wings, slightly thinner (or hollow) carbon rods can be chosen for the leading edge, to increase the spanwise deformation. 5. Conclusion Simulations have been performed using an IBM solver to investigate the underlying aerodynamics of the flapping-wing MAV with a Delfly-Micro (DFM) biplane configuration in 3D under different conditions. These include: 1. Forward flight condition 2. Varying inclination angles 3. Prescribed span/chordwise flexibility 258 W.B. Tay et al. / Journal of Fluids and Structures 55 (2015) 237–261 Fig. 26. Pressure contours of (a) DFM, (b) DFM-c_df-0.25 and (c) DFM-c_df-0.35, obtained from the contours of pressure on a cylindrical surface of radius¼1.5c. Black directional arrows indicate in flow velocity. W.B. Tay et al. / Journal of Fluids and Structures 55 (2015) 237–261 259 Fig. 27. Total ct and cl of the DFM with dlcw-max ¼0.25c, dlsw-max ¼ 0.30c and dls/cw-max ¼0.25c over one period. Fig. 28. Pressure contours of (a) DFM-s_df¼ 0.30, (b) DFM-c_df-0.25 and (c) DFM-s/c_df-0.25, obtained from the contours of pressure on a cylindrical surface of radius¼ 1.5c at time ¼0.36 T and 0.68 T. Green circles denote differences. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.) 260 W.B. Tay et al. / Journal of Fluids and Structures 55 (2015) 237–261 Fig. 29. Vortex structures at Q¼ 250 for time ¼ 0.36 T of (a) DFM-s_cf¼ 0.30, (b) DFM-c_df ¼ 0.25 and (c) DFM-s/c_df ¼ 0.25. Table 8 Average total ct, Pin and η of DFM, DFM-s_df ¼0.30, DFM-c_df-0.25 and DFM-s/c_df-0.25. Configuration type Average total ct Average total Pin Total η DFM DFM-s_df-0.30 DFM-c_df-0.25 DFM-s/c_df-0.25 0.31 0.39 0.74 0.87 2.75 3.40 3.64 4.24 0.11 0.11 0.20 0.20 Results show that aerodynamic performance of the DFM flapping configuration is in general superior to that of DF2 at a Re of 5000. The former provides more thrust while maintaining the same efficiency. However, to achieve this performance, the DFM requires a more powerful motor to flap its wings due to the higher power input, compared to a similar sized FMAV using the smaller flapping stroke angles of the DF2. The higher rocking amplitude of the DFM (compared to DF2) during flight may affect the images captured, although this problem can be relieved by using a suitable image stabilization algorithm. The DFM generates more lift than expected when the inclination angle increases, due to the generation of a spanwise uniform leading edge vortex. At a body inclination angle of 151, the DFM generates lift equivalent to its weight, indicating that sustained forward flight is indeed possible for the DFM. Increasing the spanwise flexibility of the wings causes the LEV to attach more strongly to the wings, which in turn increases the thrust. This improvement in the leading edge vortex stabilization is due to a unique combination of spanwise flow, three-dimensional kinematics, rotation, and/or flexibility. On the other hand, increasing the chordwise flexibility of the wings from dlcw-max ¼0 to 0.35c increases the efficiency and thrust, with maximum thrust occurring at dlcw-max ¼0.25c. This is due to two opposing factors – the appearance of trailing edge vortex which decreases the thrust, and the increased wing pitching angle due to chordwise flexibility which increases the thrust. However, the leading edge vortex has almost disappeared. An additional simulation with both spanwise and chordwise flexibility shows that this combination outperforms cases undergoing either only spanwise or chordwise flexibility. As mentioned earlier, reducing the rigidity of the leading edge carbon rod may further improve the wing's performance, as indicated in the chordwise deformation results. Currently, this study investigates the kinematics and wing flexibility from a more fundamental aspect. Future research will shift towards simulating the actual DFM, which include the modeling of the wing as thin membranes instead of plates and inclusion of fluid structure interaction to replicate the actual wing deformation during flapping. All these results provide a further understanding in the underlying aerodynamics of the DFM, which can help to improve the performance of FMAVs that are employing this unique propulsion configuration. Acknowledgments The work is carried out under support of the Dutch Technology Foundation (STW). The authors wish to thank STW for the financial support, Grant number 11023. The authors also wish to thank Asst. Prof. K.B. Lua and S. Deng for providing the drag coefficient data for the solver validation and the slow motion videos of the DFM's wing respectively. Part of this work is supported by the AnSTAR Computational Resource Centre through the use of its high performance computing facilities. W.B. Tay et al. / Journal of Fluids and Structures 55 (2015) 237–261 261 References Aono, H., Liang, F., Liu, H., 2008. Near- and far-field aerodynamics in insect hovering flight: an integrated computational study. Journal of Experimental Biology 211, 239–257, http://dx.doi.org/10.1242/jeb.008649. Aono, H., Chimakurthi, S., Cesnik, C.E.S., Liu, H., Shyy, W., 2009. Computational modeling of spanwise flexibility effects on flapping wing aerodynamics. In: Proceedings of the 47th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition. Orlando, Florida, pp. 1–18. Balay, S., Gropp, W.D., McInnes, L.C., Smith, B.F., 1997. Efficient management of parallelism in object oriented numerical software libraries. In: Arge, E., Bruaset, A.M., Langtangen, H.P. (Eds.), Modern Software Tools in Scientific Computing. pp. 163–202. Beem, H.R., Rival, D.E., Triantafyllou, M.S., 2012. On the stabilization of leading-edge vortices with spanwise flow. Experiments in Fluids 52, 511–517, http: //dx.doi.org/10.1007/s00348-011-1241-9. Birch, J.M., Dickinson, M.H., 2001. Spanwise flow and the attachment of the leading-edge vortex on insect wings. Nature 412, 729–733, http://dx.doi.org/ 10.1038/35089071. Calderon, D.E., Wang, Z., Gursul, I., 2010. Lift enhancement of a rectangular wing undergoing a small amplitude plunging motion. In: Proceedings of the 48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition. Orlando, Florida, pp. 1–18. De Croon, G.C.H.E., de Clercq, K.M.E., Ruijsink, R., Remes, B., de Wagter, C., 2009. Design, aerodynamics, and vision-based control of the DelFly. International Journal of Micro Air Vehicles 1, 71–97. Deng, S., van Oudheusden, B.W., Remes, B.D.W., Percin, M., Bijl, H., Ruijsink, H.M., 2013. Experimental investigation of the flapping performance on “ Delfly Micro,” In: Proceedings of the International Micro Air Vehicle Conference and Flight Competition. Toulouse, pp. 17–20. Ellington, C.P., 1984. The aerodynamics of hovering insect flight. III: kinematics. Philosophical Transactions of the Royal Society B 305, 41–78. Ellington, C.P., 1999. The novel aerodynamics of insect flight: applications to micro-air vehicles. Journal of Experimental Biology 202, 3439–3448. Ellington, C.P., van den Berg, C., Willmott, A.P., Thomas, A.L.R., 1996. Leading-edge vortices in insect flight. Nature 384, 626–630. Falgout, R.D., Jones, J.E., Yang, U.M., 2006. The design and implementation of hypre, a library of parallel high performance preconditioners. In: Bruaset, A.M., Tveito, A. (Eds.), Numerical Solution of Partial Differential Equations on Parallel Computers, Springer-Verlag, pp. 267–294. Heathcote, S., Wang, Z., Gursul, I., 2008. Effect of spanwise flexibility on flapping wing propulsion. Journal of Fluids and Structures 24, 183–199, http://dx. doi.org/10.1016/j.jfluidstructs.2007.08.003. Hunt, J.C.R., Wray, A.A., P. Moin, 1988. Eddies, streams, and convergence zones in turbulent flows. In: Proceedings of the Summer Program. pp. 193–208. Keennon, M., Klingebiel, K., Won, H., Andriukov, A., 2012. Tailless flapping wing propulsion and control development for the nano hummingbird micro air vehicle. In: Proceedings of the American Helicopter Society Future Vertical Lift Aircraft Design Conference. San Francisco, pp. 1–24. Kim, D., Choi, H., 2000. A second-order time-accurate finite volume method for unsteady incompressible flow on hybrid unstructured grids. Journal of Computational Physics 162, 411–428, http://dx.doi.org/10.1006/jcph.2000.6546. Kim, J., Kim, D., Choi, H., 2001. An immersed-boundary finite-volume method for simulations of flow in complex geometries. Journal of Computational Physics 171, 132–150, http://dx.doi.org/10.1006/jcph.2001.6778. Lee, J., Kim, J., Choi, H., Yang, K.-S., 2011. Sources of spurious force oscillations from an immersed boundary method for moving-body problems. Journal of Computational Physics 230, 2677–2695, http://dx.doi.org/10.1016/j.jcp.2011.01.004. Leonard, B.P., 1979. A stable and accurate convective modelling procedure based on quadratic upstream interpolation. Computer Methods in Applied Mechanics and Engineering 19, 59–98. Liao, C.-C., Chang, Y.-W., Lin, C.-A., McDonough, J.M., 2010. Simulating flows with moving rigid boundary using immersed-boundary method. Computers & Fluids 39, 152–167, http://dx.doi.org/10.1016/j.compfluid.2009.07.011. Lua, K.B., Lim, T.T., Yeo, K.S., 2014. Scaling of aerodynamic forces of three-dimensional flapping wings. AIAA Journal 52, 1095–1101, http://dx.doi.org/ 10.2514/1.J052730. Mittal, R., Iaccarino, G., 2005. Immersed boundary methods. Annual Review of Fluid Mechanics 37, 239–261, http://dx.doi.org/10.1146/annurev. fluid.37.061903.175743. Mohd-yusof, J., 1997. Combined immersed-boundary/B-spline methods for simulations of flow in complex geometries. CTR Annual Research Briefs, Center for Turbulence Research, NASA Ames/Stanford Univ., 317–327. Nakata, T., Liu, H., Tanaka, Y., Nishihashi, N., Wang, X., Sato, a, 2011. Aerodynamics of a bio-inspired flexible flapping-wing micro air vehicle. Bioinspiration and Biomimetics 6, 045002, http://dx.doi.org/10.1088/1748-3182/6/4/045002. Pedro, G., Suleman, A., Djilali, N., 2003. A numerical study of the propulsive efficiency of a flapping hydrofoil. International Journal for Numerical Methods in Fluids 42, 493–526. Pornsin-sirirak, T.N., Lee, S.W., Nassef, H., Grasmeyer, J., Tai, Y.C., Ho, C.M., Keennon, M., 2000. Mems wing technology for a battery-powered ornithopter. In: Proceedings of the 13th IEEE Annual International Conference on MEMS. Miyazaki, pp. 709–804. Tay, W.B., Lim, K.B., 2010. Numerical analysis of active chordwise flexibility on the performance of non-symmetrical flapping airfoils. Journal of Fluids and Structures 26, 74–91, http://dx.doi.org/10.1016/j.jfluidstructs.2009.10.005. Tay, W.B., van Oudheusden, B.W., Bijl, H., 2012. Validation of immersed boundary method for the numerical simulation of flapping wing flight. In: Proceedings of the International Micro Air Vehicle Conference and Flight Competition. Braunschweig, pp. 1–21. Tay, W.B., van Oudheusden, B.W., Bijl, H., 2014. Numerical simulation of X-wing type biplane flapping wings in 3D using the immersed boundary method. Bioinspiration and Biomimetics 9, 036001, http://dx.doi.org/10.1088/1748-3182/9/3/036001. Weis-Fogh, T., 1973. Quick estimates of flight fitness in hovering animals, including novel mechanisms for lift production. Journal of Experimental Biology 59, 169–230. Yang, J., Balaras, E., 2006. An embedded-boundary formulation for large-eddy simulation of turbulent flows interacting with moving boundaries. Journal of Computational Physics 215, 12–40, http://dx.doi.org/10.1016/j.jcp.2005.10.035. Yang, W., Song, B., Song, W., Wang, L., 2012. The effects of span-wise and chord-wise flexibility on the aerodynamic performance of micro flapping-wing. Chinese Science Bulletin 57, 2887–2897, http://dx.doi.org/10.1007/s11434-012-5249-1. Young, J., Walker, S.M., Bomphrey, R.J., Taylor, G.K., Thomas, A.L.R., 2009. Details of insect wing design and deformation enhance aerodynamic function and flight efficiency. Science 325, 1549–1552, http://dx.doi.org/10.1126/science.1175928. Zhu, Q., 2007. Numerical simulation of a flapping foil with chordwise or spanwise flexibility. AIAA Journal 45, 2448–2457, http://dx.doi.org/10.2514/ 1.28565.