Home

Search

Collections

Journals

About

Contact us

My IOPscience

Numerical simulation of X-wing type biplane flapping wings in 3D using the immersed boundary method

This content has been downloaded from IOPscience. Please scroll down to see the full text. 2014 Bioinspir. Biomim. 9 036001 (http://iopscience.iop.org/1748-3190/9/3/036001) View the table of contents for this issue, or go to the journal homepage for more

Download details: IP Address: 130.207.50.37 This content was downloaded on 11/11/2014 at 23:26

Please note that terms and conditions apply.

Bioinspiration & Biomimetics Bioinspir. Biomim. 9 (2014) 036001 (21pp)

doi:10.1088/1748-3182/9/3/036001

Numerical simulation of X-wing type biplane flapping wings in 3D using the immersed boundary method W B Tay 1,2 , B W van Oudheusden 1 and H Bijl 1 1

Faculty of Aerospace Engineering, Delft University of Technology, Kluyverweg 1, 2629 HS Delft, The Netherlands 2 Temasek Laboratories, National University of Singapore, 5A Engineering Drive 1, 117411, Singapore E-mail: [email protected] Received 19 August 2013, revised 14 January 2014 Accepted for publication 31 January 2014 Published 3 March 2014 Abstract

The numerical simulation of an insect-sized ‘X-wing’ type biplane flapping wing configuration is performed in 3D using an immersed boundary method solver at Reynolds numbers equal to 1000 (1 k) and 5 k, based on the wing’s root chord length. This X-wing type flapping configuration draws its inspiration from Delfly, a bio-inspired ornithopter MAV which has two pairs of wings flapping in anti-phase in a biplane configuration. The objective of the present investigation is to assess the aerodynamic performance when the original Delfly flapping wing micro-aerial vehicle (FMAV) is reduced to the size of an insect. Results show that the X-wing configuration gives more than twice the average thrust compared with only flapping the upper pair of wings of the X-wing. However, the X-wing’s average thrust is only 40% that of the upper wing flapping at twice the stroke angle. Despite this, the increased stability which results from the smaller lift and moment variation of the X-wing configuration makes it more suited for sharp image capture and recognition. These advantages make the X-wing configuration an attractive alternative design for insect-sized FMAVS compared to the single wing configuration. In the Reynolds number comparison, the vorticity iso-surface plot at a Reynolds number of 5 k revealed smaller, finer vortical structures compared to the simulation at 1 k, due to vortices’ breakup. In comparison, the force output difference is much smaller between Re = 1 k and 5 k. Increasing the body inclination angle generates a uniform leading edge vortex instead of a conical one along the wingspan, giving higher lift. Understanding the force variation as the body inclination angle increases will allow FMAV designers to optimize the thrust and lift ratio for higher efficiency under different operational requirements. Lastly, increasing the spanwise flexibility of the wings increases the thrust slightly but decreases the efficiency. The thrust result is similar to one of the spanwise studies, but the efficiency result contradicts it, indicating that other flapping parameters are involved as well. Results from this study provide a deeper understanding of the underlying aerodynamics of the X-wing type, which will help to improve the performance of insect-sized FMAVs using this unique configuration. Keywords: biplane, flapping wings, immersed boundary method, numerical (Some figures may appear in colour only in the online journal)

maneuverability and speed. The popularity of flapping wing micro-aerial vehicles (FMAVs) has led to numerous studies [1, 2] focusing on different areas. We will briefly discuss some of these areas, which include unsteady aerodynamic effects,

1. Introduction Insects and birds have long been inspiring us to create flying vehicles which resemble them, due to their excellent 1748-3182/14/036001+21$33.00

1

© 2014 IOP Publishing Ltd

Printed in the UK

W B Tay et al

Bioinspir. Biomim. 9 (2014) 036001

(a)

(d)

(b)

(e)

(c)

(f )

the flapping motion of a foil with either chordwise or spanwise flexibility. The simulations were performed immersed in both low and high density fluids. In the first scenario, it was found that chordwise flexibility reduces both the thrust and the propulsion efficiency, while the spanwise flexibility increases the thrust without efficiency reduction, within a small range of structural parameters. In the second scenario, chordwise flexibility increases the efficiency, while spanwise flexibility reduces both the thrust and the efficiency. Another numerical investigation by Aono et al using a computational framework for FSI found that there are two key factors associated with spanwise wing deformation which affect thrust generation. They are the in-phase motion between the wing tip and root, and the increased effective angle of attack of the deformed wing. Under the conditions whereby the wing motions resulting from both prescribed motion and deformation are correlated, increasing the effective angle of attack at the tip enhances the aerodynamic performance. However, high wing flexibility decreases the aerodynamic performance. The results of Heathcote et al [8] and Zhu [7] clearly contradict one another in the high density (water) scenario. Hence, depending on different conditions, spanwise and chordwise flexibility may or may not be beneficial to the performance of the FMAV. Wing–wing interactions refer to the interactions when two wings or a wing and a tail are in close proximity of one another. The clap and fling motion [4] described earlier is the consequence of wing–wing interaction. In two-dimensional flow simulations, Tuncer and Kaya [10] used an overset grid method to simulate the unsteady viscous flow over flapping airfoils in a biplane configuration using a Reynolds-averaged Navier–Stokes solver. The Spalart–Allmaras turbulence model was used to simulate turbulent flow. The Reynolds number, Re, was 10 000 (10 k) and the airfoils were flapping at a heaving amplitude of 0.4 (non-dimensionalized with respect to the airfoil chord) and reduced frequency f varying from 0.08 to 0.48. The biplane configuration was found to produce up to 40% more thrust than a single (flapping) airfoil. Tay et al [11] used an immersed boundary method (IBM) [12] solver to simulate biplane flapping airfoils at Re of 1 k. It was found that the overall efficiency and average thrust per airfoil increase up to 17% and 126%, respectively, with respect to the single airfoil due to the formation of a strong momentum jet that results from the interaction of the two airfoils. With regards to the design of FMAV, novel materials and mechanisms have now been used in different aspects. Studies have been performed to investigate the use of smart materials, such as piezoelectric materials [13] or ionic polymer metal composite [14] for the FMAVs’ wing. Energy storage flapping mechanism, similar to that of an insect thorax, has been proposed for improved efficiency [15] as well. Besides FMAV with a single pair of wings [16, 17], there are other variants as well [18, 19]. One of these is the Delfly [20], a family of bio-inspired ornithopter MAV developed at the Delft University of Technology, which has two pairs of wings in a X-wing type flapping configuration (figure 2). Although it is unlike any particular species of flying insect or bird found in nature, simple experiments [20], using rubberband powered flying models to measure the flying distance,

Figure 1. Schematic illustration of the clap-and-fling/peel motion; streamlines are illustrated in black, light blue vectors are net forces and dark blue vectors represent the induced velocities. ‘Clap and fling’ adapted with permission from [6] © 2003 Journal of Experimental Biology. ‘Peel’ added in by the authors.

wing flexibility, wing–wing interaction, smart materials, novel mechanisms and flapping configurations. One of the unsteady aerodynamic effects is the leading edge vortex (LEV), which was first visualized in an experiment by Ellington [3], who used a 10 × scale mechanical flapping model. The LEV, created by dynamic stall, was used to explain the high lift generated by insects. It was also postulated that the LEV managed to stabilize throughout the up or downstroke due to the positive spanwise velocity removing vorticity from the LEV. Another unsteady aerodynamic effect is the clapand-fling motion [4], which comprises two phases, as shown in figure 1. In the ‘clap’ phase, the leading edges of the wings come together rotating about the leading edges, occurring until the gap between the wings disappears. As the gap between the wings closes progressively, both wings’ circulations cancel each other out, resulting in air being expelled out in the form of a momentum jet. In the ‘fling’ phase, the wings rotate about their trailing edges forming a gap in between, sucking air into the gap. This suction causes an increase in circulation and hence an increase in lift. A variation of the clapand-fling motion is the clap-and-peel motion (figure 1). The difference between the classic clap-and-fling and the clap-andpeel motion is that the ‘fling’ is replaced by the ‘peel’, whereby the wings peel apart due to the fluid–structure interaction (FSI) between the air and the flexible membrane wings, resulting in a more gradual build-up of the circulation, which can prevent an unstable LEV from shedding into the wake [5]. This motion can be found in the Lepidoptera and Drosophila. Insect wings are flexible and studies [7–9] have shown that spanwise and chordwise flexibility can be beneficial or detrimental to force production and efficiency. Heathcote et al performed a water tunnel study to investigate the effect of spanwise flexibility on the thrust, lift and propulsive efficiency of a rectangular wing oscillating in pure heave. A small increase in thrust coefficient and higher efficiency are achieved for a degree of spanwise flexibility when the Strouhal number is greater than 0.2. However, a further increase in the flexibility decreases the thrust and efficiency significantly. Zhu [7] performed a fully coupled FSI numerical study to investigate 2

W B Tay et al

Bioinspir. Biomim. 9 (2014) 036001

X-wing type flapping configuration in 3D in forward flight condition, which motivates this study to be performed under this condition. A reference simulation running at the original Re of 36 k is attempted too. However, as will be explained in the results section, the current laminar numerical solver is not suitable and the results obtained do not reflect the performance of the actual Delfly. Due to the complexity in the design of the Delfly X-wing type flapping configuration, it will be difficult to perform an exact simulation, which requires the capability of FSI and wing material modeling. Hence, the main objective of this study is to perform a simplified analysis of the X-wing type flapping configuration in 3D in forward flight condition which captures the essential flow features and allows us to estimate the aerodynamic performance. The simplifications used in this simulation study include using thin plates, increasing the distance between the upper and lower wings at the root to prevent wing intersection. An IBM [12] laminar Navier–Stokes solver 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 of one another, which is a situation for which the IBM approach is very suited, as will be discussed below. There are four objectives to be achieved in this study. They are as follows.

Figure 2. Different versions of the Delfly family: I (right), II (left), micro (middle).

flight time, and required rubber band winding, have shown that it has lower power requirement compared to FMAVs with a single pair of wings or two wings in tandem. Moreover, it gives minimal rocking amplitude (vertical oscillation, perpendicular to the FMAVs flight path), which is a beneficial property for FMAVs to be used for sharp image capture and recognition. Numerical and experimental studies have been performed on the Delfly to improve our understanding on its unique X-wing type flapping configuration. De Clercq et al [21] performed particle image velocimetry (PIV) measurements and simultaneous force measurements on the hovering Delfly II to investigate the flow-field behavior and its relation to the aerodynamic forces generated. Similarly, Groen et al [22] performed a parametric wing geometry study in combination with a detailed aerodynamic and aeroelastic investigation. In both Delfly experimental studies, the clap-and-peel motion was observed, due to the flexibility of the wings. In the numerical aspect, Nakata et al [23] performed simulations on a hummingbird-inspired, flexible wing FMAV (very similar to the Delfly micro, the smallest FMAV in the Delfly family) using an overset solver under hovering conditions. A strong negative pressure region accompanied by a LEV was found on the upper and lower wings during both of the half strokes. The excellent properties of the Delfly X-wing type flapping configuration motivate us to perform numerical simulation on an insect-sized FMAV using the Delfly wing configuration in the forward flight condition. As mentioned earlier, the advancement in smart materials, together with MEMS technology enables the fabrication of smaller and smaller FMAVs, such as the flapping wing microrobot [17] from the Harvard Microrobotics Laboratory. We would like to assess its aerodynamic performance when the original Delfly II FMAV is reduced to the size of an insect, which corresponds to the Reynolds number (Re) between 1 k and 5 k. The Delfly II (hereafter referred to simply as the Delfly) is specifically chosen because it is an established and welldocumented platform on which a number of experiments [21, 22, 24] have been performed to characterize it. The current smallest X-wing type FMAV is the Delfly micro, with a wingspan of 10 cm and a weight of 3 g, corresponding to the Re of around 5 k. Moreover, to the authors’ best knowledge, there has not been any numerical study performed on the

(1) To assess the possible aerodynamic advantage of an insect-sized X-wing type flapping configuration at Re = 1 k, as compared with the more common single wing flapping configurations. As mentioned in the previous sections, simulations of biplane configurations in 2D [10, 11] showed that there is thrust increase compared to the single airfoil case. Hence, we would like to simulate and prove that there is similar benefit in the 3D case as well at Re = 1 k. (2) To evaluate and compare the effects of Re at 1 k and 5 k. The Re of 1 k and 5 k correspond approximately to the flight envelope of honeybee and hawkmoth, respectively [2]. Although both Re still belong to the laminar regime, it is expected that finer vortical structures will appear for the higher Re case. Previous studies [25–27] have shown that output forces generally increase as Re increases. Therefore, it is expected that output forces at Re = 5 k will be higher than their counterparts at Re = 1 k. (3) During forward flight of FMAVs like the Delfly, the body inclination angle, which is the angle of attack of the body with respect to the incoming flow (figure 3), is always nonzero to vector the forces, such as to produce sufficient thrust and lift at the same time. However, most biplane simulations [10, 11] have been performed for a zero body angle. We would like to evaluate if this has a large influence on the aerodynamic forces and flow structures. (4) Lastly, the current numerical solver does not have FSI capability. The addition of prescribed spanwise wing deformation will allow us to study and compare the effect of flexibility with respect to rigid wings. Past studies [7–9] have found that depending on different conditions, spanwise flexibility may or may not be beneficial to the performance of the FMAV. Hence, this motivates 3

W B Tay et al

Bioinspir. Biomim. 9 (2014) 036001

validated to give good accuracy compared to grid conforming solvers and experiments in laminar flow from Re = 40 to 800 [28–30] and turbulent flow at Re = 4 k and 10 k. It has also been shown to give better accuracy at lower grid resolution compared to the continuous forcing approach [31]. A more detailed description will be given in the later section. The main drawback of the IBM is that the size of the Cartesian grid used has to increase much faster compared to structured grids when the Re increases, in order to resolve the flow gradients near the surface [12]. As a result, the computational resource requirement is much higher than using conforming grids. 2.1. Discrete forcing IBM

The current discrete forcing IBM is based on a combination of the methods developed by Yang and Balaras [28], Kim et al [29] and Liao et al [30]. A summary of the current method is given below. The reader is referred to the above-mentioned papers for more detailed descriptions. In the discrete forcing IBM, the momentum equation is modified with the addition of the forcing term (fc) to become ∂u 1 2 = −u · ∇u + ∇ u − ∇ p + f c, (1) ∂t Re where u is the velocity vector, t is the time, p is the pressure and Re is the Reynolds number. Equation (1) has been nondimensionalized using the far field incoming velocity (U∞) and root chord length (c) as the reference velocity and length. Accordingly, the calculation of Re is based on U∞ and c. In the paper by Liao et al [30], the forcing term is provisionally calculated explicitly using the second-order Adam–Bashforth (AB2) schemes for both the viscous and convective terms. However, in this study, the forcing term is provisionally calculated explicitly using the first-order forward Euler and second-order AB2 schemes for the viscous and convective terms, respectively. This is because the forcing term will be equivalent to the right-hand side of the linearized momentum equation, which is the negative of the sum of the convective, viscous and pressure terms in equation (1). This will reduce the computational cost and it is found that there is no observable difference in the results compared to the original scheme, when using the forward Euler/AB2 explicit combination scheme. It is given by   u f − un 3 n 1 n−1 f cn+1 = − Dn + ∇ pn , (2) + C − C t 2 2

Figure 3. Explanation of the body inclination angle.

us to investigate the effects of flexibility on the X-wing 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 vorticity contour plots. It is expected that the results obtained from this study will further enhance the understanding of the aerodynamic behavior of the insect-sized X-wing type FMAVs and as such be beneficial to their further improvement. 2. Numerical method The unsteady viscous flow generated by the flapping wings 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 small. Moreover, the moving wings require the re-calculation of the grid at every timestep, which will be expensive in a 3D simulation. 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 relatively straightforward. However, since the grid does not conform to the outlines of the body, the boundary conditions around the grid require special modifications. In fact, it is the different modifications used by the various research groups that distinguish them. They are generally classified into the continuous and discrete forcing approach [12]. The difference between the two is that the modification or forcing term is added into the governing equation before and after discretization in the continuous and discrete forcing approach, respectively. In this study, the discrete forcing approach [28–30] is used because it enables a sharp representation of the body contour, which is important when simulating higher Re cases of 5 k. The discrete forcing approach has been

where C and D are spatial operators containing the convective 1 ∇ 2 u terms, respectively. ∇ · uu and viscous Re Here, uf, the velocity of the interface point, is obtained through interpolation of the boundary point (u1) plus three surrounding points (u2–4), using the triangle linear interpolation in 3D. An example in 2D (without u4 in the third dimension) is shown in figure 4. The 3D interpolation scheme is given by u f = a1 u1 + a2 u2 + a3 u3 + a4 u4 ,

(3)

where a1–4 are the coefficients calculated using the boundary point plus three surrounding points. 4

W B Tay et al

Bioinspir. Biomim. 9 (2014) 036001

 Fi = − solid

f cn+1 dV + i

solid

∂ui u j ∂ui + ∂t ∂x j

 dV,

(5)

where V is the volume of the body. Equation (5) is based on the force exchange between control volumes of the solid and fluid. More details about its derivation can be found in [35]. This method is much simpler and faster than the first method because the forcing term has already been calculated and stored in the computer’s memory. Single and biplane simulations using the force calculation methods have been performed and there is less than 1% difference between them. Hence, in this study, all subsequent simulations and the results presented use the second method to calculate the forces. The power input for the flapping is calculated based on the angular velocity of the wing and the torque it produces. The torque is given by  Tr = (r × F ) dx, (6)

Figure 4. Triangle linear interpolation in 2D, u1 = boundary point, u2–3 are the surrounding points.

2.2. Fractional step method

A finite volume fractional step method, which is based on an improved projection method [32], is used to solve the modified non-dimensionalized incompressible Navier–Stokes equation (momentum equation (1) and mass continuity equation (4)): ∇ · u = 0.





where r and F are the distance of the point from the rotating axis and the output forces. Power input is then given by

(4)

Pn = Tr · ω,

The time integration scheme is a semi-implicit scheme, such that the convective and viscous terms use the second-order AB2 and Crank Nicolson (CN2) discretization, respectively. All spatial derivatives are discretized using central differencing and the QUICK [33] scheme on a staggered grid for Re = 1 k and 5 k, respectively. The general steps of the fractional step method are outlined briefly. The discretized momentum is first solved to obtain a velocity field which is non-divergence free. Using the non-divergence free velocity, the Poisson equation is then solved once to obtain the pressure field, which in turn updates the velocity to be divergence free. It is similar to the SIMPLE [34] algorithm, except that the SIMPLE algorithm solves the momentum and Poisson equations iteratively many times within a time step until convergence. The solver is laminar and no turbulence models have been added.

(7)

where ω is the angular velocity of the wing. Propulsive efficiency η is defined as a measure of the energy lost in the wake versus energy used in creating the necessary thrust, and is given by c¯t η= , (8) P¯in where ct is the output thrust. The overbar indicates averaged values. 2.4. Validation

The accuracy of the IBM solver has been validated by comparison with the experiment performed by Calderon et al [36] at the University of Bath. The experiment involves a plunging NACA0012 rectangular wing with a semi aspect ratio of 2, in a water tunnel. The wing, at a constant angle of attack of 20◦ , undergoes pure plunging motion at an amplitude of 0.15c while the reduced frequency f varies from 0.0 to 2.0 at a Re of 10 k. The outputs from the experiment include vorticity contour plots as well as forces measured using a twocomponent aluminum binocular strain gauge force balance. The forces included both time-dependent aerodynamic forces and inertia forces. By doing time-averaging, the inertia force during the plunging up cycle cancels with the inertia force during the plunging down cycle, leaving only the averaged aerodynamic forces. The vorticity contour plots have also been phase-averaged. For the validation, two frequencies of 0.45 and 0.9 are selected. The inflow velocity is in the z direction. A grid convergence study has been performed using three different grid sizes, as given in table 1. Only the ct at different resolutions is shown in figure 5 because the grid requirement for the lift and side forces is much lower than that for the thrust. Results show that the minimum grid spacings of 0.018, 0.006 and 0.006

2.3. Force and power calculation

Since the grid is not conformal to the body, the force coefficients3 on the body are calculated differently compared to that of a body conforming grid. When using a body conforming grid, the boundary grid coincides with the body’s surface. Hence, numerical integration can be applied directly to obtain the pressure and viscous forces. On the other hand, in the IBM, the forces are calculated in a different way. There are two methods to obtain the forces when using the IBM. The first method explicitly calculates the pressure and viscous stress using the coefficients arising from the tri-linear interpolation stencil, as given in (3). More details about this force calculation can be found in [28]. The second method uses the forcing term fcn+1 obtained earlier. In this case, the total force is calculated as 3

For simplicity, the force coefficients ct and cl are assumed to be equal to the forces Fi. This does not pose any problem because we only compare between results within this study. 5

W B Tay et al

Bioinspir. Biomim. 9 (2014) 036001

Figure 5. cd at f = 0.45 using different grid resolutions with zoom in.

Figure 6. PIV and IBM vorticity contour plots comparison at f = 0.45.

Figure 7. PIV and IBM vorticity contour plots comparison at f = 0.90. Table 1. Minimum grid lengths in x, y and z directions and total grid

Table 2. Force coefficient comparison between experimental and numerical results.

size.

Grid

Minimum Minimum Minimum grid grid grid lengths (x) lengths (y) lengths (z) Total grid size

Coarse 0.036 Medium 0.018 Fine 0.018

0.012 0.006 0.003

0.012 0.006 0.003

f

Average Average cd (Calderon) cd (IBM)

0.45 0.44 0.90 0.34

101 × 221 × 421 177 × 368 × 734 177 × 651 × 1338

Average Average cl (Calderon) cl (IBM)

0.39 (−11%) 1.36 0.24 (−29%) 1.60

1.25 (−7%) 1.31 (−18%)

and z directions. The simulation results agree reasonably with the experimental results, as shown in table 2, figures 6 and 7, although it under-predicts all the force coefficients. The discrepancy in values between the numerical and experimental results is much smaller at f = 0.45, compared to f = 0.9. It is most likely that the flow is transiting from laminar to

in the x, y and z directions, respectively, which corresponds to a grid resolution of 177 × 368 × 734, are sufficient for accurate results. 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 6

W B Tay et al

Bioinspir. Biomim. 9 (2014) 036001

(a)

(b)

(c)

Figure 8. Normal, close-up isometric and front view of the 3D Cartesian grid with the biplane wings in red.

turbulence, especially at f = 0.9, which appears in the form of fine vortical structures, as shown in figure 7. Nevertheless, the trends in forces from f = 0.45 to 0.9 are predicted correctly. In an earlier numerical/experimental comparison [31] at the same Re, it is found that higher frequencies ( f = 1.3 to 2.0) in the same experiment result in completely turbulent flow. The current simulations run at Re = 1 k and 5 k, which are still considered laminar flows.

Table 3. Convergence test using ct for one period.

dx

Root mean |rms| Order square (rms) = |rmsfine−rmscoarse|a = |rms|fine/|rms|coarse

0.018 0.012 0.009 0.006

0.566 0.624 0.597 0.613

0.0588 0.0277 0.0158

2.13 1.76

a

Fine and coarse refer to the two corresponding dx, for example, 0.012 and 0.018.

2.5. Simulation setup Table 4. Average time taken to run one simulation cycle in parallel at different dx.

This section gives the setup of the simulations for this study. Figure 8 shows the 3D Cartesian grid with the biplane wings. The orientation of the wings is such that the incoming flow is along the z axis direction. 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 use of a staggered grid and numerical truncation error may prevent the full wing from being captured. 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 figure 8(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 × 441 × 405.

dx

Processor type

Number of processors

Time taken (h)

0.018 0.012 0.009 0.006

AMD opteron 6136 AMD opteron 6136 AMD opteron 6234 AMD opteron 6234

32 64 96 144

2.5 3.25 30 107

Re = 5 k because convergence study at Re = 1 k already shows that dx = 0.018 is not sufficiently accurate. Only the comparison of ct at different resolutions is shown because the grid requirement for lift and side forces is much lower than that for the thrust. Figure 9 shows that a minimum grid length of 0.009 gives sufficient accuracy and is computationally cheaper than the minimum grid length of 0.006 at Re = 1 k and 5 k. Moreover, the convergence test using ct for one period is also given in table 3 and it shows second-order convergence. There are some small spikes in the figure, which can either be due to the IBM [35] or to the actual flow aerodynamics. Even at dx = 0.006, the spikes are present, albeit to a lesser extent. However, the number of spikes is very small and they do not affect the overall accuracy of the forces. Hence, dx = 0.009 will be used for all simulations in this study. The average time taken to run one simulation cycle in parallel at different dx is given in table 4.

2.6. Grid convergence study

A grid convergence study is performed using the basic configuration (Xw) given in the research methodology section at Re = 1 k and 5 k. Due to the non-conformal nature of the IBM solver, the grid requirement is higher than that of grid conforming solvers. The convergence study at Re = 1 k uses minimum grid lengths (dx) of 0.006, 0.009, 0.012 and 0.018 chord lengths, corresponding to grid sizes of 381 × 629 × 556, 266 × 441 × 405, 207 × 346 × 327 and 148 × 249 × 245, respectively. Subsequently, the convergence study at Re = 5 k uses minimum grid lengths of 0.006, 0.009, 0.012. The dx = 0.018 case is omitted at

3. Research methodology The main objective of this study is to gain a deeper understanding of the underlying aerodynamics of the X-wing type biplane flapping configuration. This section will first 7

W B Tay et al

Bioinspir. Biomim. 9 (2014) 036001

(a)

(b)

Figure 9. ct of the upper wing at minimum grid lengths of 0.006, 0.009, 0.012 and 0.018c at Re = (a) 1 k and (b) 5 k, with zoom in.

varying pitching of 10◦ is added to the wings during flapping. Another concern is that the design of the biplane wings of the Delfly 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, a small amount of clearance (equivalent to 0.3c) is added between the two wings to prevent the wings from intersecting one another. Despite these simplifications, it is expected that the simulation is capable of capturing the essential flow behavior and trends. Therefore, valuable insights can still be obtained from the simulation results.

discuss the simplifications used in this study, followed by the details on how the simulations are performed to achieve the objectives. 3.1. Simplifications

Due to the complexity in the design of Delfly and its wings, some simplifications have been used in the simulations. Some of these have been briefly mentioned earlier in the paper but they are repeated here for completeness, while the additional modifications discussed below are based on the experimental observations on the actual Delfly. The Delfly uses thin membrane Mylar sheets for its wings, which flex heavily depending on the FSI. In comparison, while insects’ wings are also flexible, they are much more rigid compared to the thin membrane Mylar wings of the Delfly. Hence, although the current solver does not have an FSI capability, by representing the wings by rigid thin plates (thickness = 0.06c) with or without prescribed deformation, it is not too far off from the insects’ wings approximation. Moreover, to better simulate the chordwise pitching of the wings due to FSI, a sinusoidally

3.2. Simulation parameters

In this section, the simulation parameters used in the simulations are discussed. Figure 10 shows the planform and stroke angle of the wing. It uses the same shape and stroke angle as the Delfly. As mentioned in the introduction, the Delfly is specifically chosen because it is an established and well-documented platform on which a number of experiments [21, 22, 24] have been performed to characterize it. For the 8

W B Tay et al

Bioinspir. Biomim. 9 (2014) 036001

(a)

(b)

(c)

Figure 10. (a) Planform, (b) stroke angle, layout of the X-wing (c) stroke angle, layout of the actual Delfly [37].

configurations, as given in table 6. For the comparison between the single and X-wing (Sw versus Xw) configurations, two types of single wing configurations are used. The first is the Sw-u configuration, which refers only to the upper wing of the Xw configuration. In other words, the Sw-u configuration is exactly the same as that of the Xw, minus the lower wing. The lower wing is not simulated because both the upper and lower wing configurations are very similar. Another single wing configuration is the Sw-full configuration. In this case, the single wing undergoes the full stroke angle of θ o = 44◦ , similar to that of the Xw configuration. In other words, the tip of the single wing at its top-most and bottom-most flapping positions coincides with the tip of the upper and lower wing of the Xw configuration at the end of its outstroke, respectively. Its flapping parameters are given in table 5. The Xw-5 k and Xw-36 k configurations have exactly the same parameters as the Xw configuration except that they are simulated at Re = 5 k and 36 k, respectively. Moreover, the Xw-36 k configuration has a reduced frequency4 of 0.126, which is calculated based on the actual flapping frequency of 11 Hz, forward flight velocity of 7 m s−1 and chord length of 0.08 m. For the body inclination comparison, the Xw configuration is rotated along the x axis to give body inclination angle α b of 15◦ or 30◦ (Xw-15◦ or 30◦ ). Lastly, the Xwdeform-0.15/0.30 configurations have flexible wings along the spanwise direction with the maximum deformed length (dlmax) at 0.15/0.30c. The next section will describe the spanwise flexibility of the wing in more detail.

Table 5. Wing angles of the upper, lower and full wings in degrees.

θ mid (degrees)

Wing angles Upper Lower Full

30 −8 14

θ o (degrees) 22 22 44

model wings, the flapping stroke and pitch angles are given respectively by θ = θ0 sin(2π f (t − t0 ) − φ) + θmid α = α0 sin(2π f (t − t0 )) + αb ,

(9) (10)

where θmid , θ0 , θ , αb , α0 , α, f and φ are the midpoint stroke angle, stroke amplitude, stroke angle, body inclination angle, pitch amplitude, pitch angle, reduced frequency and phase angle, respectively. Their values for the upper and lower wings of the X-wing are given in table 5, which are based on the Delfly II specification given in appendix A of [38]. The pitch amplitude and the center of rotation of the pitching are at 10◦ and 0.5c, respectively. The selection of the reduced frequency f is based on the criterion that average thrust > 0. Frequencies of 0.15, 0.3 and 0.6 are tested and it is found that simulations at f = 0.15 and 0.3 produce drag. At f = 0.6, thrust is generated and hence it is selected as the flapping frequency in this study. The phase angle φ is fixed at 90◦ . 3.3. Wing configurations

In a previous section, four objectives to be achieved in this study were mentioned. These include single/X-wing configuration comparison, Re number 1 k/5 k comparison, body inclination effect and spanwise flexibility effect. In order to achieve the objectives, the wings are simulated in different

3.4. Spanwise flexibility

A front view of the prescribed deformation of the wings (Xw-deform) at a particular instant in time is shown in 4

The reduced frequency f is defined as actual frequency × c/U∞ .

Table 6. Different wing simulation configurations.

Configuration type

Simulation description

Sw-u Sw-full Xw (-u/d) Xw-5 k Xw-36 k Xw-15◦ /30◦ Xw-deform-0.15/0.30

Single upper wing at Re = 1 k without wing deformation Single wing equivalent to full stroke angle θ o = 44◦ of Xw at Re = 1 k without wing deformation Basic X-wing without body inclination or wing deformation at Re = 1 k (upper/lower) Basic X-wing without body inclination or wing deformation at Re = 5 k Basic X-wing without body inclination or wing deformation at Re = 36 k, with f = 0.126 X-wing with body inclination angle α b = 15◦ or 30◦ at Re = 1 k without wing deformation X-wing with wing deformation dlmax = 0.15 or 0.30c at Re = 1 k 9

W B Tay et al

Bioinspir. Biomim. 9 (2014) 036001

Table 7. Total ct , Pin and η of the Sw-u, Sw-full and Xw

configurations. Configuration Total ct

Total cl Total Pin

0.58 (2 × 0.29)a −0.094 1.95 (+236%) 0.62 0.77 (+33%) 0.045

Sw-u Sw-full Xw

Total η

6.8 (2 × 3.4) 0.085 28.3 0.069 8.82 0.087

a The ct of Sw-u has been multiplied by a factor of two to ensure a fair comparison between Sw-u, Sw-full and Xw.

4. Results and discussions The simulations are run for five periods. Results show that the forces becomes 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 [41] vorticity plots. These plots enable one to identify and visualize the 3D vortex structures. The results and discussions will be presented in the following four sections according to the objectives of the study.

Figure 11. Front view of the right side of the X-wing configuration 5

at a particular instant in time .

figure 11. In this study, the upper and lower wings undergo sinusoidal spanwise deformation with a maximum deformed length (dlmax) of 0.15c or 0.30c. A maximum value of 0.3c is chosen because past studies on chordwise prescribed flexibility of airfoils [39, 40] show that dlmax = 0.3c gives the optimum configuration. However, the value of 0.3c is used only as a reference comparison since the deformations in this study are spanwise. Moreover, the objective of the deformation study is to compare the differences between rigid and flexible wings, instead of obtaining the optimum configuration. The spanwise flexibility equation is given by  2 dc dl = dlmax sin(2π f (t − t0 )), (11) dsemispan

4.1. Single and X-wing flapping configuration comparison

The vortical flow structures are visualized through the Q criterion [41] iso-surface. Figures 12, 13 and 14 show the evolution of the vortex structures of Sw-u, Sw-full and Xw at Q = 250 over one period (1 T), respectively. The upstroke and downstroke start from a to c and d to f , respectively. Figure 15 shows the force coefficients outputs ct and cl of the Sw-u, Swfull and Xw configurations over one period. Figure 16 shows the moment about the z axis (momz) and power input (Pin) of the Sw-u, Sw-full and Xw configurations over one period. The labeled vertical lines represent the same time instant as in figures 12, 13 and 14. Table 7 shows the total ct , Pin and η of the Sw-u, Sw-full and Xw configurations. In figure 12 of the Sw-u configuration, at the start of the upstroke at time t = 0.12 T, a spiral conical LEV, with increasing strength from root to tail, starts to grow at the lower surface of the wing. The reason for the varying vortex strength

where dl, dlmax, dc and dsemispan are the perpendicular deformed length magnitude, maximum deformed length, distance from a point on the wing to the root and semi wingspan of the wing respectively. dl always acts opposite to the direction of flapping. The f used is the same as the reduced flapping frequency of 0.6. 5

The front view of the wings appear much thicker in this case due to their pitching.

(a)

(b)

(c)

(e)

(d)

(a) (f)

Figure 12. Vortex structures of Sw-u at Q criterion (Q) = 250 over one period (1 T). The double black arrow represents the inflow direction. The green, purple and blue circles represent the LEV, TV and TEV respectively. 10

W B Tay et al

Bioinspir. Biomim. 9 (2014) 036001

(a)

(b)

(c)

(f)

(e)

(d)

Figure 13. Vortex structures of Sw-full at Q = 250 over one period.

is due to the increase in the flapping velocity of the wing along the span, which results in a spanwise variation of the effective angle of attack. This LEV is accompanied by the appearance of the tip vortex (TV) and partial trailing edge vortex6 (TEV) as well. The TEV is only present at the slanted edge of the wing, near the tip and it is much weaker in strength. The growth of these vortices increases until the mid-upstroke, after which they break down and shed, due to the decrease in the flapping velocity and the increase in the negative spanwise velocity near the tip. The shedding of the vortices causes a peak in the thrust, as seen around the b instant of figure 15(a). As mentioned in the introduction, Ellington [3] postulated that the conical LEV managed to stabilize throughout the up or downstroke due to the positive spanwise velocity removing vorticity from the conical LEV. However, in this case, the spanwise velocity becomes increasingly negative near the tip, as shown in figure 17. The vorticity is not removed from the conical LEV and the buildup of vorticity reaches a point whereby it becomes unstable, resulting in the conical LEV breakup. The same event happens at the downstroke (starting from figure 12(d)), whereby a conical LEV and TV start to grow in strength at the upper surface of the wing, accompanied by a weak TEV. In figure 12(e), the LEV is at its maximum strength, which corresponds to the instant of highest lift in e of figure 15(b). The vortices shed near the end of the downstroke.

The vortical iso-surface result of figure 12 is very similar to the simulation of a thin flapping wing in forward flight at Re = 10 k by Gopalakrishnan and Tafti [42]. In their study, the LEV also did not extend up to the base and breaks up before the end of the up or downstroke, due to the negative spanwise velocity. On the other hand, the vortex structure of the Swfull configuration in figure 13 shows a more chaotic vortex shedding, although the general features of the vortex structure such as the conical LEV are similar to that of the Sw-u configuration. Due to the higher stroke amplitude of 44◦ , the flapping velocity of the wings is much higher, resulting in the generation of stronger conical LEV, as shown in figure 13(e). Moving on to the X-wing (Xw) configuration, figure 14(a) denotes the start of the outstroke of the Xw configuration, with the conical LEV appearing at the lower surface of the upper wing and at the upper surface of the lower wing. At the same time, the TV and TEV are also growing in strength for both upper and lower wings. The conical LEV tube grows in magnitude as one moves from the root to the tip of the wing, similar to the Sw-u and Sw-full configurations. Moreover, ringlike LEV tubes can be seen in Xw, nearer the wing tip region of the wings. They are clearly missing in the Sw-u and Swfull configurations. The tip vortices in the Xw configuration are stronger compared to that of Sw-us, but weaker than that of Sw-fulls. The two differences can be seen more clearly in figure 18 at two particular snapshots, using a lower Q criterion of 50. The green and purple circles denote the differences in the LEV. The bottom view of the surface pressure of the

6 The TEV is only visible from another perspective angle and it is not shown here.

11

W B Tay et al

Bioinspir. Biomim. 9 (2014) 036001

(a)

(b)

(c)

(f)

(e)

(d)

Figure 14. Vortex structures of Xw at Q = 250 over one period. The yellow circles represent the ring-like LEV.

Figure 15. ct and cl of the Sw-u, Sw-full, Xw-u and Xw configuration over one period.

wing of the Sw-u, Sw-full and Xw-u (upper wing of Xw) configuration over one period is shown in figure 19. The top view is not shown because the difference is more pronounced in the bottom view. The lower pressure between the upper and lower wings of Xw around the instroke period provides higher suction and hence stronger LEV and TV compared to the Sw-u configuration. However, the Sw-full configuration shows the largest variation in the pressure due to having a flapping stroke angle (θ 0) twice (44◦ ) that of the Sw-u and Xw configurations (22◦ ). Figure 15 shows the ct and cl of the Sw-u, Sw-full and Xw configuration over one period. Between t = 0.12 T and 0.24 T, in figures 12(a) and (b), 13(a) and (b), and 14(a) and

(b), the strength of LEV at the (upper) wing’s lower surface is growing stronger, which is evident from the corresponding minimum cl and maximum ct peak in figure 15. After this instant (figure 15(b)), the vortical structures shed, resulting in a decrease in the ct. At the start of the instroke (after t = 0.5 T), the LEV and TV once again build up, leading to an increase in the cl and ct, followed by the shedding of the vortices after figure 15(d). At figure 15(d), the LEV, TV and TEV then shed, giving a second peak in the ct. This second peak of the Xw configuration is higher than the first, unlike in the Sw-u and Sw-full configuration, where both peaks are equal. This is the result of the upper and lower wings clapping together. As observed in figure 19(c), at the start of the instroke, the 12

W B Tay et al

Bioinspir. Biomim. 9 (2014) 036001

Figure 16. momz and Pin of the Sw-u, Sw-full, Xw-u and Xw configuration over one period.

similar pressure change throughout the flapping stroke, their changes in magnitudes are symmetrical in the instroke and outstroke phases. It is interesting to note that in the hovering simulation of the Delfly micro-variant by Nakata et al [23], there are also two unequal thrust peaks. However, in their case, the thrust at the outstroke is higher than that at the instroke. The two possible reasons are different flapping parameters and flight conditions. The Delfly micro-variant wings clap together both during the outstroke and the instroke, unlike that of the Delfly, which only clap together during the instroke. Moreover, the flight condition in this study is forward flight while in Nakata et al [23], it is hovering flight. Similarly, Lehmann et al [43, 44] performed clap and fling experiments using a dynamically scaled electromechanical flapping device in the quiescent condition. By varying the kinematic pattern, the higher lift peaks can either happen at the instroke or outstroke. Moreover, in the clap-and-fling (or clap-and-peel) motion, the center of rotation changes from the leading edge to the trailing edge while transiting from clap to fling, but in this study, the wings rotate about the mid chord position. The shedding of the vortices is shown in figure 20 at two radial slices of 0.3 and 1.35. The shedding of a pair of opposite vortices can be seen in both instants for the Xw configuration and it is not observed in the Sw-u and Sw-full configurations. A similar vortex pair is also observed in the 2D biplane airfoil simulations performed by Tay et al [11] at the same Re of 1 k. Hence, the biplane simulation results in 2D match the 3D ones very closely. This is also true with respect to the total average thrust (c¯t ) comparison. There is an increase in the total average thrust (c¯t ) from 0.58 (2 × 0.29) (Sw-u) to 0.77 (Xw) (+33%), similar to the thrust improvement in the 2D biplane simulations by Tuncer [10] and Tay et al [11]. The improvement in thrust is even greater for the Sw-full configuration, which gives ct = 1.95, 236% more than that of the Sw-u configuration. The strong conical LEV seen in figure 13 generates high thrust, as a result of flapping at twice the stroke angle (θ 0 = 44◦ ). This result tallies with that of the simple experiment using rubberband powered flying models, whereby the average flight speed of the single wing model (2.35 m s−1) is higher than that of the X-wing model (1.40 m s−1), due to the higher wing loading. The total cl of Sw-u, Sw-full and Xw is very small ( 0. This gives a Strouhal number (St) of 0.78, 1.44 and 0.72/0.788 for the Sw-u, Sw-full and Xw configurations. According to Taylor et al [45], the optimum St at cruising speed for birds and insects is between 0.2 and 0.4. The St of the current study is much higher than 0.4 and hence it is expected to be low for the η. However, given that the c¯t of the Xw configuration is higher than that of the Sw-u configuration, while the η remains almost constant for both configurations, this indicates that the Xw configuration performs better than the Sw-u configuration in this study. Nevertheless, it should be emphasized that more research using different flapping parameters must be performed to give a more conclusive statement since the study by Nakata et al [23] showed that either the single or X-wing configurations can give higher lift forces depending on the mylar membrane

(a)

(b)

(c) Figure 20. Vortices shedding at the radius r = 0.3 and r = 1.35 for

(a) Sw-u, (b) Sw-full and (c) Xw configurations at different instants in time. Green circles denote regions of differences.

8 Depending on how the St is calculated, it can have different values. The calculation of St can be found in the appendix.

14

W B Tay et al

Bioinspir. Biomim. 9 (2014) 036001

(a)

Figure 21. Drag output of the lower and upper wing of the Xw-36 k

configuration.

thickness and the body inclination angle of the FMAV. On the other hand, despite the much higher ct of the Sw-full configuration, its η is lower than that of the Sw-u and Xw configurations. Its momz and Pin are also higher than the other two configurations. This means that in an actual FMAV, to achieve the higher thrust, a more powerful motor must be used. Moreover, the overall structure of the FMAV must be able to withstand the higher vibration. The comparison between single wing and X-wing configurations shows that at Re = 1 k, there are also performance advantages in the X-wing configuration, hence making it an attractive alternative design for insect-sized FMAVs.

(b) Figure 23. X vorticity slice at x = 0.1 and 1.0 for Re of (a) 1 k and (b) 5 k. Green circles denote regions of differences.

actual performance of the Delfly. Hence, the result from the Xw-36 k configuration is not used in the comparison. Figure 22 shows the vortex structures of Xw at Re = 1 k and 5 k (Xw-5 k) over one period with Q criterion at 250. On comparison, the predominant features in Xw, such as the conical LEV, are also present in Xw-5 k. However, the vortex structures at Xw-5 k are smaller and finer. This is expected due to the lower viscosity at higher Re. The same variation of vortex structures across increasing Re has also been found in many other studies, such as Garmann et al [27], who performed simulation from Re = 200 to 60 k of a rotating plate. X vorticity slices at x = 0.1 and 1.0 for Re of 1 k and 5 k in figure 23 show more diffusive wake

4.2. Re = 1 k/5 k comparison

This section describes the differences between the simulation results at Re = 1 k and 5 k. A reference simulation running at the original Re of 36 k is attempted. However, the flow field has now transitioned to fully turbulent and the current laminar solver is no longer suitable. As shown in figure 21, there are many spikes present in the drag output. A small amount of drag instead of thrust is produced and does not reflect the

(a)

(b)

Figure 22. Vortex structures of Xw at (a) Re = 1 k (b) Re = 5 k (Xw-5 k) over one period at Q = 250. 15

W B Tay et al

Bioinspir. Biomim. 9 (2014) 036001

(a)

(b) Figure 24. ct and cl for the (a) upper and (b) lower wing at Re = 1 k and 5 k.

(a)

Table 8. Total ct , ct and η of the Xw and Xw-5 k configurations.

(b)

Configuration

Total ct

Total cl

Total η

Xw Xw-5 k

0.77 0.89 (+16%)

0.045 0.046

0.087 0.098 (+13%)

do not defer much in terms of force production and efficiency. In fact, several studies show that force output variation is very similar when using a laminar model compared to a turbulent model at Re < 60 k [46]. This means that it may be possible to estimate the actual force generated by FMAV at higher Re by doing simulation at a lower Re. The general important features found in the vorticity contour plots of the Xw-5 k configuration are also present in the Xw-1 k configuration, although it is smaller and finer. 4.3. Body inclination angle comparison Figure 25. Surface pressure of the wings at Re = 1 k and 5 k at

The body inclination angle α b of all FMAVs in forward flight is almost never zero. The slower the forward speed, the higher the body inclination angle. Previous 2D biplane simulations [10, 11] were performed with a zero angle of attack. Hence, in this section, we would like to evaluate if there is a large difference in the aerodynamic forces and flow structures when α b = 15◦ and 30◦ . Figure 26 shows the ct and cl of the upper wing of the Xw configuration (Xw-u) at α b = 0◦ , 15◦ and 30◦ over one period. The ct and cl of the lower wing of the Xw configuration are similar and hence they are not shown here. The change in the cl graphs is expected; the graphs shift up as α b increases, resulting in an increasingly positive cl . This is because part of the thrust component has been vectored to become lift. However, at first glance, the shape of ct graphs looks very different as α b increases. This is because it is not possible to match the shape of the graphs at α b = 15◦ and 30◦ to that of α b = 0◦ by a simple shift. However, by rotating the ct and cl of the α b = 0◦ case by 30◦ and resolving them into components along the vertical and horizontal axes (figure 27), one can see that the resultant graph’s shape becomes very similar to that of the α b = 30◦ case, as shown in figure 28. Despite the similarity now in the graphs’ shapes, there are still some differences in the peak magnitudes of ct and cl. The resulting ct of the Xw-u configuration at α b = 0◦ rotated 30◦ and α b = 30◦ are 0.32 and −0.16, respectively. On the other hand, the cl at α b =

time = (a) 0.2 T and (b) 0.8 T. Green circles denote differences.

at Re = 5 k. Moreover, at x = 1.0, the LEV at Re = 5 k is less well defined, due to mixing of positive and negative vorticity. Figure 24 shows the ct and cl for the (a) upper and (b) lower wing at Re = 1 k and 5 k. The shapes of the graphs are generally similar. The total ct at Re = 5 k is larger than that of Re = 1 k by 16% (0.89 compared to 0.77). Increasing ct as a result of increasing Re was also observed by Garmann et al [27] in the rotating plate simulation. However, in this case, the ct increase is limited to the region from the start of the outstroke to the mid-instroke. These differences are reflected more clearly in the surface pressure of the wing, as shown in figure 25. The distribution of the surface pressure is very similar for both Re, except for the higher suction (lower pressure) at the wing tip for the Xw-5 k configuration at t = 0.2 T. This may have explained the ct increase during this instant. At t = 0.8 T, the pressure distribution on the wing surface for both Re are almost the same, corresponding to similar ct for both Re configuration at that instant. There is also no difference in the cl at both Re. As a result of this and the higher ct for the Xw-5 k configuration, the η is also slightly higher at 0.098, compared to Xw’s 0.087. The total ct , ct and η of the Xw and Xw-5 k configurations are summarized in table 8. These results show that FMAVs in the Re range between 1 k and 5 k 16

W B Tay et al

Bioinspir. Biomim. 9 (2014) 036001

(a)

(b) Figure 26. (a) ct and (b) cl of the Xw-u configuration at α b = 0◦ , 15◦ and 30◦ . Table 9. Total ct , cl and force magnitude of the Xw, Xw-15◦ /30◦

and Xw rotated 15◦ /30◦ configurations. Configuration

Total ct

Xw 0.77 0.25 Xw-15◦ −0.72 Xw-30◦ Xw rotated 15◦ 0.73 Xw rotated 30◦ 0.66

Total cl

Resultant force magnitude

0.045 1.94 3.30 0.15 0.33

0.77 1.96 3.38 (with drag) 0.75 0.74

Figure 27. Rotation of Xw configuration by 30◦ and resolving the original ct and cl into new components.

instead of thrust, despite the high lift coefficient. Similar to the upper wing discussion earlier, the total ct and cl of the rotated configurations are much higher and lower compared to the Xw-15/30◦ configurations, respectively. The total force magnitude for the Xw-15◦ configurations is also higher than that of the Xw configuration. Figure 29 shows the vortex structures of Xw-30◦ at Q = 250 over one period. The labeled vertical lines in figure 28 represent the same time instants as in figure 29.

0◦ rotated 30◦ and α b = 30◦ are 0.18 and 1.32, respectively. This shows that increasing α b to 30◦ causes a reduction in the expected thrust, but an improvement over the expected lift. This result also holds true for the α b = 15◦ case, although it is not shown here. The total ct , cl and resultant force magnitude of the Xw, Xw-15/30◦ and Xw rotated 15/30◦ configurations are shown in table 9. The Xw-30◦ configuration produces drag

(a)

(b) Figure 28. (a) ct and (b) cl of the Xw-u configuration at α b = 0◦ , 30◦ and 0◦ rotated to 30◦ .

(a)

(b)

(c)

(d)

Figure 29. Vortex structures of Xw-30◦ at Q = 250 over one period. 17

W B Tay et al

Bioinspir. Biomim. 9 (2014) 036001

(a)

(b) Figure 30. (a) LEV and (b) TV differences between Xw and Xw-30◦ configurations.

(a)

Table 10. Total ct , cl , Pin and η of the Xw-u, Xw-deform-0.15 and

(b)

Xw-deform-0.30 configurations. Configuration

Total ct

Total Total cl Pin Total η

Xw 0.77 0.045 8.82 0.087 Xw- deform-0.15 0.78 (+1.3%) 0.067 9.92 0.078 (−10.3%) Xw- deform-0.30 0.80 (+3.9%) 0.015 11.34 0.071 (−18.4%)

to the Xw configuration. This additional knowledge regarding the aerodynamics performance of the Xw configuration at 15◦ and 30◦ will allow FMAV designers to better adjust the body inclination angle α b to optimize the thrust and lift ratio for higher efficiency under different operational requirements.

Figure 31. Vertical velocity v comparison at z = 0.15 c from the

wing’s leading edge at α b = (a) 0◦ and (b) 30◦9.

Maximum thrust and lift now occur at the start of the outstroke (figures 28(a) and 29(a) and the mid-instroke10 (c) respectively. In comparison with the vortex structures of Xw in figure 14, both configurations are similar in terms of overall vortex structures. However, the Xw-30◦ configuration shows more chaotic structures. The vortices are also more unstable and shed earlier, compared to the Xw configuration, as observed in the upper wing of figures 14(a) and 29(a). It is also observed in the lower wing of figures 14(d) and 29(c). This is expected since it is easier for flows to be separated more easily at higher angles of attack. Two other key differences lie in the LEV and TV, as shown in figure 30. In figure 30(a), the size of the LEV for the Xw-30◦ configuration is uniform across the entire wingspan, unlike the conical LEV in figure 14. As mentioned in the Swu/Xw configuration comparison, the increasing strength of the conical LEV along the wingspan is due to the increasing flapping velocity. In the current case, the wing is constantly at a high angle of attack (30◦ ) and hence flow will be separated even at the root, similar to flow past a stationary wing at a high angle of attack. However, unlike the stationary wing, the LEV tube is formed throughout the wingspan and it remains stable for approximately 0.2 T with the help of the spanwise velocity. Moreover, it can be observed that the LEV is always formed at the top surface of the upper and lower wings. This asymmetry therefore resulted in the higher than expected lift, as shown in table 9. The high angle of attack also increases the strength of the TV, as shown in figure 31. In the Xw-30◦ configuration, the v velocity around the tip is higher (green circle), compared

4.4. Spanwise flexibility comparison

Figure 32 shows the vortex structures of the Xw-deform-0.15 and Xw-deform-0.30 configuration at Q = 250 over one period. In comparison to figure 14 at the same time instants, the most obvious difference is that the formation of the conical LEV is delayed when the wing is spanwise flexible. Moreover, higher flexibility causes more delay. In figure 33, the ct and cl graphs of the original and deformed configurations for both the upper and lower wings show that the shapes of the different configuration are generally similar, except for a slight shift horizontally to the right. The ct graph of the upper wing for the Xw-deform0.30 configuration also has a slight dip at around t = 0.25 T. The reason for the shift is due to the delay in the conical LEV generation, which is in turn due to the reduction in the velocity of the wing perpendicular to the spanwise direction. This can be seen more clearly by differentiating the deformation equation (11), after which one will obtain the deformation velocity as  2 dc vdl = dlmax (2π f ) cos(2π f (t − t0 )). (12) dsemispan Equation (12) acts in the opposite direction to the wing flapping direction, hence reducing the final velocity in the deformed configuration. From table 10, one can observe that the ct increases by a small amount (up to 3.9%) as dlmax increases. The magnitude of the lift changes as dlmax increases but it is too small relative to its peak value ( 0.30 but 19

W B Tay et al

Bioinspir. Biomim. 9 (2014) 036001

efficiency remains constant as spanwise flexibility increases, contrary to the current results. The flapping parameter involved in the two studies are different, and hence, it is not possible to do a direct comparison between these different studies.

force output variation is very similar when using a laminar model compared to a turbulent model at Re < 60 k [46], indicating that it may be possible to estimate the actual force generated by FMAV at higher Re by doing simulation at a lower Re. Together with experiments and perhaps further numerical optimization study, the reliance on trial and error in the design of FMAV can be greatly reduced and improvement in the performance of insect-sized FMAVs can be achieved in a shorter period.

5. Conclusion Simulations have been performed using an IBM solver to investigate the underlying aerodynamics of an insect-sized Xwing type (Xw) flapping configuration in 3D at Re = 1 k and 5 k. The Q criterion iso-surface visualization of the Xw configuration reveals interesting features such as the ring-like LEV tubes not found in the single wing (Sw-u and Sw-full) configurations. The average thrust ct of the Xw configuration is more than twice that of the Sw-u configuration. Similarly, given that the Xw’s efficiency η is similar to that of Sw-u’s, this proves that the Xw configuration is superior in performance compared to the Sw-u configuration. It must however be noted that these conclusions are based on the particular set of flapping parameters used in this study. Hence, additional research is required to reach a more conclusive statement. On the other hand, although the Xw configuration provides lower thrust compared to the Sw-full configuration, it gives higher η and lower vibration and power input, making it more suited for sharp image capture and recognition. Hence, the performance advantages in the X-wing configuration make it an attractive alternative design for insect-sized FMAVS compared to the single wing configuration. Q criterion iso-surface comparison between the Xw configuration at Re of 1 k and 5 k shows that both configurations show similar predominant features, but the vortex structures at Xw-5 k are smaller and finer due to vortex breakup. ct and η of the Xw-5 k configuration are also slightly higher, due to the higher suction at its wing tip region. It is not possible to do a comparison with the actual Delfly configuration at Re = 36 k because the flow field is completely turbulent. When the body inclination angle α b increases from 0◦ to 30◦ , the shape of the thrust and lift curves at α b = 15/30◦ are very similar to that of the zero body inclination angle configuration rotated over 15/30◦ . However, the values of ct and cl are very different in these two cases. This is because at high α b, instead of conical LEV, a uniform LEV tube is formed on the wing. Moreover, the LEV is only formed at the top surface of the upper and lower wings, resulting in much higher lift. Understanding the force variation as the body inclination angle increases will allow FMAV designers to optimize the thrust and lift ratio for higher efficiency under different operational requirements. Lastly, comparison of rigid and spanwise flexible wings in this study shows that increasing the spanwise flexibility of the wings increases the thrust slightly but decreases the efficiency, similar partly to the results of Zhu’s spanwise study [7]. However, many parameters are at play and they are interlinked together, resulting in conflicting results from different studies. The design of an actual FMAV is a very complex process. In the past, most groups rely on trial and error in their FMAV design. The current study shows that we can achieve a deeper understanding in the underlying aerodynamics of the X-wing type configuration through simulation. Studies have shown that

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. Appendix. Strouhal number (St) calculations The St is calculated according to Taylor et al [45], whereby St is given by fA , (A.1) St = U∞ where f , A and U∞ are the flapping frequency, wing amplitude and far field incoming velocity. The computation of A, the wing amplitude of Xw is more complicated because of the unique X-wing configuration. In this case, A can be defined in two ways—the distance from the upper most position of the upper wing to the lowest most position of the lower wing, or simply the amplitude of the upper or lower wing. In the first method, A has to be divided by two since it involves two wings. In the first method, 1.75 × (sin 58◦ ) + 1.75 × (sin 30◦ ) = 1.2 A= 2 1.2 × 0.6 = 0.72. 1.0 In the second method, St =

A = 2 × 1.75 × sin 22◦ = 1.3 St =

1.3 × 0.6 = 0.78. 1.0

References [1] Ho S, Nassef H, Pornsinsirirak N, Tai Y and Ho C 2003 Unsteady aerodynamics and flow control for flapping wing flyers Prog. Aerosp. Sci. 39 635–81 [2] Shyy W, Aono H, Chimakurthi S K, Trizila P, Kang C-K, Cesnik C E S and Liu H 2010 Recent progress in flapping wing aerodynamics and aeroelasticity Prog. Aerosp. Sci. 46 284–327 [3] Ellington C P 1999 The novel aerodynamics of insect flight: applications to micro-air vehicles J. Exp. Biol. 202 3439–48 http://jeb.biologists.org/content/202/23/3439.abstract [4] Weis-fogh T 1973 Quick estimates of flight fitness in hovering animals, including novel mechanisms for lift production J. Exp. Biol. 59 169–230 http://jeb.biologists.org/content/59/1/169.abstract 20

W B Tay et al

Bioinspir. Biomim. 9 (2014) 036001

[27] Garmann D J, Visbal M R and Orkwis P D 2012 Three-dimensional flow structure and aerodynamic loading on a low aspect ratio, revolving wing 42nd AIAA Fluid Dynamics Conf. and Exhibit (New Orleans, Louisiana) pp 1–22 [28] Yang J and Balaras E 2006 An embedded-boundary formulation for large-eddy simulation of turbulent flows interacting with moving boundaries J. Comput. Phys. 215 12–40 [29] Kim J, Kim D and Choi H 2001 An immersed-boundary finite-volume method for simulations of flow in complex geometries J. Comput. Phys. 171 132–50 [30] Liao C-C, Chang Y-W, Lin C-A and McDonough J M 2010 Simulating flows with moving rigid boundary using immersed-boundary method Comput. Fluids 39 152–67 [31] Tay W B, van Oudheusden B W and Bijl H 2012 Validation of immersed boundary method for the numerical simulation of flapping wing flight Int. Micro Air Vehicle Conf. and Flight Competition (Braunschweig) pp 1–21 [32] Kim D and Choi H 2000 A second-order time-accurate finite volume method for unsteady incompressible flow on hybrid unstructured grids J. Comput. Phys. 162 411–28 [33] Leonard B P 1979 A stable and accurate convective modelling procedure based on quadratic upstream interpolation Comput. Methods Appl. Mech. Eng. 19 59–98 [34] Patankar S V and Spalding D B 1972 A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows Int. J. Heat Mass Transf. 15 1787–806 [35] Lee J, Kim J, Choi H and Yang K-S 2011 Sources of spurious force oscillations from an immersed boundary method for moving-body problems J. Comput. Phys. 230 2677–95 [36] Calderon D E, Wang Z and Gursul I 2010 Lift enhancement of a rectangular wing undergoing a small amplitude plunging motion 48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition (Orlando, FL) pp 1–18 [37] Bruggeman B 2010 Improving flight performance of DelFly II in hover by improving wing design and driving mechanism MSc Thesis Delft University of Technology [38] Groen M A 2010 PIV and force measurements on the flapping-wing MAV DelFly II MSc Thesis Delft University of Technology [39] Miao J-M and Ho M-H 2006 Effect of flexure on aerodynamic propulsive efficiency of flapping flexible airfoil J. Fluids Struct. 22 401–19 [40] Tay W B and Lim K B 2010 Numerical analysis of active chordwise flexibility on the performance of non-symmetrical flapping airfoils J. Fluids Struct. 26 74–91 [41] Hunt J C R, Wray A A and Moin P 1988 Eddies, streams, and convergence zones in turbulent flows Proc. Summer Program pp 193–208 http://adsabs.harvard.edu/abs/1988stun.proc..193H [42] Gopalakrishnan P and Tafti D K 2009 Effect of rotation kinematics and angle of attack on flapping flight AIAA J. 47 2505–19 [43] Lehmann F-O, Sane S P and Dickinson M 2005 The aerodynamic effects of wing-wing interaction in flapping insect wings J. Exp. Biol. 208 3075–92 [44] Lehmann F-O and Pick S 2007 The aerodynamic benefit of wing-wing interaction depends on stroke trajectory in flapping insect wings J. Exp. Biol. 210 1362–77 [45] Taylor G K, Nudds R L and Thomas A L R 2003 Flying and swimming animals cruise at a Strouhal number tuned for high power efficiency Nature 425 707–11 [46] Lian Y, Broering T, Hord K and Prater R 2013 The characterization of tandem and corrugated wings Prog. Aerosp. Sci. 65 41–69

[5] Ellington C P 1984 The aerodynamics of hovering insect flight: III. Kinematics Philos. Trans. R. Soc. B 305 41–78 [6] Sane S P 2003 The aerodynamics of insect flight J. Exp. Biol. 206 4191–208 [7] Zhu Q 2007 Numerical simulation of a flapping foil with chordwise or spanwise flexibility AIAA J. 45 2448–57 [8] Heathcote S, Wang Z and Gursul I 2008 Effect of spanwise flexibility on flapping wing propulsion J. Fluids Struct. 24 183–99 [9] Aono H, Chimakurthi S, Cesnik C E S, Liu H and Shyy W 2009 Computational modeling of spanwise flexibility effects on flapping wing aerodynamics 47th AIAA Aerospace Sciences Meeting Including The New Horizons Forum and Aerospace Exposition (Orlando, FL) pp 1–18 [10] Tuncer I H and Kaya M 2003 Thrust generation caused by flapping airfoils in a biplane configuration J. Aircr. 40 509–15 [11] Tay W B, Bijl H and van Oudheusden B W 2012 Analysis of biplane flapping flight with tail 42nd AIAA Fluid Dynamics Conf. and Exhibition (New Orleans, Louisiana) pp 1–17 [12] Mittal R and Iaccarino G 2005 Immersed boundary methods Annu. Rev. Fluid Mech. 37 239–61 [13] Mukherjee S and Ganguli R 2010 Non-linear dynamic analysis of a piezoelectrically actuated flapping wing J. Intell. Mater. Syst. Struct. 21 1157–67 [14] Mukherjee S and Ganguli R 2010 Ionic polymer metal composite flapping actuator mimicking dragonflies Comput. Mater. Contin. 19 105–34 [15] Madangopal R, Khan Z A and Agrawal S K 2005 Biologically inspired design of small flapping wing air vehicles using four-bar mechanisms and quasi-steady aerodynamics J. Mech. Des. 127 809–16 [16] Pornsin-sirirak T N, Tai Y C, Nassef H and Ho C M 2001 Titanium-alloy MEMS wing technology for a micro aerial vehicle application Sensors Actuators A 89 95–103 [17] P´erez-Arancibia N O, Ma K Y, Galloway K C, Greenberg J D and Wood R J 2011 First controlled vertical flight of a biologically inspired microrobot Bioinspir. Biomim. 6 036009 [18] Jones K D, Bradshaw C J, Papadopoulos J and Platzer M F 2003 Development and flight testing of flapping-wing propelled micro air vehicles 2nd AIAA ‘Unmanned Unlimited’ Systems, Technologies, and Operations (San Diego, CA) pp 1–9 [19] Festo 2013 BionicOpter–Inspired by dragonfly flight [20] De Croon G C H E, de Clercq K M E, Ruijsink R, Remes B and de Wagter C 2009 Design, aerodynamics, and visionbased control of the DelFly Int. J. Micro Air Veh. 1 71–97 [21] Clercq K M E, De Kat R, Remes B, van Oudheusden B W and Bijl H 2009 Aerodynamic experiments on DelFly II: unsteady lift enhancement Int. J. Micro Air Veh. 1 255–62 [22] Groen M, Bruggeman B, Remes B, Ruijsink R, van Oudheusden B and Bijl H 2010 Improving flight performance of the flapping wing MAV DelFly II Int. Micro Air Vehicle Conf. and Competitions (Braunschweig) pp 1–17 [23] Nakata T, Liu H, Tanaka Y, Nishihashi N, Wang X and Sato a 2011 Aerodynamics of a bio-inspired flexible flapping-wing micro air vehicle Bioinspir. Biomim. 6 045002 [24] Eisma J 2012 Flow visualization and force measurements on a flapping-wing MAV DelFly II in forward flight configuration MSc Thesis Delft University of Technology [25] Birch J M 2004 Force production and flow structure of the leading edge vortex on flapping wings at high and low Reynolds numbers J. Exp. Biol. 207 1063–72 [26] Lentink D and Dickinson M H 2009 Biofluid dynamic scaling of flapping, spinning and translating fins and wings J. Exp. Biol. 212 2691–704

21