SUMMARY
A threedimensional numerical avian model was developed to investigate the unsteady and turbulent aerodynamic performance of flapping wings for varying wingbeat frequencies and flow velocities (up to 12 Hz and 9 m s^{–1}), corresponding to a reduced frequency range of k=0.22 to k=1.0 and a Reynolds number range of Re=16,000 to Re=50,000. The wings of the birdinspired model consist of an elastic membrane. Simplifying the complicated locomotion kinematics to a sinusoidal wing rotation about two axes, the main features of dynamic avian flight were approximated. Numerical simulation techniques of fluid–structure interaction (FSI) providing a fully resolved flow field were applied to calculate the aerodynamic performance of the flapping elastic wings with the Reynolds averaged Navier–Stokes (RANS) approach. The results were used to characterize and describe the macroscopic flow configurations in terms of starting, stopping, trailing and bound vortices. For high reduced frequencies up to k=0.67 it was shown that the wake does not consist of individual vortex rings known as the discrete vortex ring gait. Rather, the wake is dominated by a chain of elliptical vortex rings on each wing. The structures are interlocked at the starting and stopping vortices, which are shed in pairs at the reversal points of the wingbeat cycle. For decreasing reduced frequency, the results indicate a transition to a continuous vortex gait. The upstroke becomes more aerodynamically active, leading to a consistent circulation of the bound vortex on the wing and a continuous spanwise shedding of small scale vortices. The formation of the vortices shed spanwise in pairs at the reversal points is reduced and the wake is dominated by the tip and root vortices, which form long drawnout vortex structures.
INTRODUCTION
The biomechanics of flying birds and insects plays a crucial role in the domain of aeronautics. In view of the complexity of wingbeat kinematics, numerical studies of bird aerodynamics are limited. This is also the case for flow and solid simulation of simplified flapping wing motion, calculating aerodynamics and corresponding wing deformations at high Reynolds numbers^{a} (Re) and flapping frequencies, which has not been performed before. In the last few decades, there have been growing efforts and successes in describing the wake configuration of flying birds using theoretical models and experimental visualization techniques. Two idealized wake configurations for different flight conditions have been established: the continuous vortex gait (Rayner, 1986; Rayner and Gordon, 1998; Spedding, 1987) and the discrete vortex ring gait (Rayner, 1979a; Rayner, 1979b; Spedding, 1986). The latter considers the wake as a series of separate vortex rings as a consequence of an inactive upstroke and an active downstroke. In this case, the starting and stopping vortices shed at the beginning and end of each downstroke are connected via tip vortices yielding a separate elliptical or circular vortex ring. This flow configuration represents the wake for slowflying birds and was the basis for a wellknown theoretical model to determine the power output of flying birds (Rayner, 1979a; Rayner, 1979b). It was experimentally validated by visualizing fragments of ring vortex structures in the wake of slowflying passerines (Kokshaysky, 1979) and pigeons and jackdaws (Spedding et al., 1984; Spedding, 1986).
For higher flow velocities (equal flight speeds), a second type of wake configuration has been observed for bats (Rayner, 1986) and kestrels (Spedding, 1987). Instead of discrete vortex rings, the wake is dominated by tip vortices of nearly constant circulation following the path of the wing tips. Both the downstroke and upstroke are aerodynamically active and, thus, continuous spanwise shedding of small scale vortex structures occurs throughout the wingbeat cycle. By means of wind tunnel tests using particle image velocimetry (PIV), it was shown that the shifting between the two idealized wake configurations is not an abrupt process, but rather a continuous transition including intermittent vortex structures and a change of corresponding aerodynamic quantities (Spedding et al., 2003). Instead of a discrete vortex gait, doubleringed wakes were found for small flow velocities. The same visualization technique was used to investigate the wake configurations of birds and bats for a wide range of flow conditions (Spedding et al., 2003; Hedenström et al., 2006; Hedenström et al., 2007; Hennigson et al., 2008; Rosén et al., 2004). Twodimensional velocity distributions were used to reconstruct threedimensional wake models which differ enormously from the idealized wakes. For birds with an inflexible wing motion and during acceleration at high and medium flight speeds, a ladder wake is assumed (Pennycuick, 1988; Hedrick et al., 2002). In this case, the wake is dominated by the tip vortices, which are connected by spanwise vortices shed at the reversal points. The circulation of the shed spanwise vortices alternates during the upstroke and downstroke motions.
Instead of living animals, models of insects (Van Den Berg and Ellingtion, 1997a; Van Den Berg and Ellingtion, 1997b) and birds (Hubel and Tropea, 2005), and model wings (Videler et al., 2004; Dickinsion and Götz, 1993; Poelma et al., 2006) have been used to describe aerodynamic performance for low Reynolds number ranges or quasi stationary flow conditions. The present paper reports on the macroscopic flow configuration of an avian flight model for a wide range of flow velocities and wingbeat frequencies representing high Reynolds numbers and a high degree of unsteadiness. The results were obtained from interacting fluid–structure simulations and should help us to understand the complex flow mechanism behind flapping wings.
MATERIALS AND METHODS
Avian flight model
A threedimensional model was developed to simulate the flapping flight of birds in a simplified manner (Fig. 1). The model consists of two elastic wings, a rigid body and a mounting suspension. The modelling was influenced by specifications of FSI simulations (e.g. topologically identical model surfaces, convergence) and corresponding grid generation aspects (e.g. avoiding sharp edges, large cell skewness). To facilitate wind tunnel tests with an identical mechanical model, restrictions in manufacturing and mechanical realization were taken into account. The body and wings were inspired by birds, but do not represent real geometrical entities. The wingspan is b=0.543 m and the corresponding wing area A=313.8 cm^{2}, leading to an aspect ratio of AR=6.57. To avoid geometrical and numerical interference, the wings were placed sideways, creating a small gap between the wing and the body. The body was extended by semicircular protrusions to prevent flow through this gap. The airfoil of the wing is characterized by a round leading edge which fades into an asymmetric and elongated chamber. Its length starts to decline linearly at the halfwing distance and it has a mean chord of c=0.08 m. The wing consists of a stiff adapter and an elastic membrane. The adapter forms the front of the wing and covers about 20% of the projected wing area. It has a round leading edge and a rectangular support for the elastic membrane on the opposite side. Furthermore, a stiff frame was added sideways to the adapter at the wing base to stabilize the trailing edge, as required for the simulation. The adapter is the active element of the wing, transferring a predefined wingbeat motion to the front side of the elastic membrane, which was mounted between the adapter and the frame. Except in these wing areas, the membrane acts passively and its movement depends on the initial and aerodynamic forces occurring. Geometrical data of the model are set out in Table 1.
Downstroke kinematics of the wing are independent of the reduced frequency^{b} k, whereas upstroke kinematics differ enormously with k (Tobalske, 2000; Tobalske, 2007; Park et al., 2001; Oertel, 2008). For simplification, the wing motion can be reduced to a rotation about a single hinge (Azuma, 2006). Neglecting the forward velocity of the wing relative to the body, two degrees of freedom are used to realize a simplified wingbeat motion (Hubel and Tropea, 2005). In the real wingbeat kinematics of birds and bats (Rosén et al., 2004; Tian et al., 2006), the flapping motion can be approximated as a periodic function. Therefore, a sinusoidal rotation about an offcentre axis parallel to the body axis provides the flapping motion of the model. The angle of attack was controlled by phasedelayed pitching about the moving lateral wing axis. The flapping angle was γ=sin(ωt) and the angle of attack was θ=cos(ωt) with the angular frequency ω=2πf and the wingbeat and pitch frequency f. For the preformed simulations, the maximum and minimum angles were γ_{min/max}=+/–30 deg and θ_{min/max}=+/–20 deg.
Numerical models
The technique applied for simulation of the FSI requires the generation of two numerical models with topologically identical surfaces for the coupling approach: a flow model and a structural model.
The structural model includes the geometrical entities of the wing and its corresponding mechanical properties. Using tetra cells the grid generation was performed automatically resulting in a good quality mesh (cell average shape factor ASF=0.7). The components were modelled as homogeneous materials with isotropic material behaviour, constant densities ρ_{m} and linear Young's modulus E. For all simulations, the mechanical properties of the adapter, frame and membrane are consistent. The various mechanical properties are given in Table 2. To implement the wingbeat motion in the structural model, the predefined periodic functions were applied, transferring the wing kinematics to the stiff adapter.
The flow model represents the fluid domain and contains the wings and the body. An unstructured mesh of tetra cells for the first time step of the simulation was built by applying an automatic grid generation technique. Because of the moving and deforming wing surface, the grid needed to be remeshed close to the model at each calculation time step. The remeshing and smoothing algorithms for the dynamic mesh generation were provided by the FLUENT solver (ANSYS Germany GmbH, Darmstadt, Germany). In the wake the grid was not affected by the moving geometry and consequently no rebuilding algorithm was applied. Grid sensitivity studies were performed for an identical geometry setup with cell numbers ranging from 240,000 to 2,480,000. Fig. 2 shows polar curves for eight different meshes at reduced frequencies of k=0.67 and k=0.22. In the range 1,020,000 to 2,480,000 cells, a nearly unchanged force distribution and consistent vortex shedding occur. The mean force coefficient m_{max} for the maximum standard deviation SD_{max}, the maximum standard error of the mean SEM_{max} and the timeaveraged standard error of the mean SEM_{ta} for one wingbeat cycle are given in Table A1 in Appendix 2.
For grids with cell numbers above 310,000, the large scale vortices produced by the flapping motion dominate the wake topology and do not differ strongly. However, cell numbers below 720,000 indicated a low resolution of the near wake and led to an unsatisfactory analysis of the vortices structures. Fig. 3 shows the velocity magnitude for increasing cell numbers at two points in the wake. Point p_{1} was placed 0.3 m downstream from the trailing edge of the wing and point p_{2} 0.6 m downstream at the midwing position. As can be seen in Fig. 3, a converging solution was obtained for cell numbers above 720,000. Statistics for the velocity magnitude v and pressure difference Δp=p–p_{∞} (where p is pressure and p_{∞} is ambient pressure) at these two points are given in Table A2 in Appendix 2. Therefore, a grid with 1,780,000 cells was used to simulate the aerodynamics and to provide a fully resolved flow field for the analysis.
ALE formulation and coupling strategy
For numerical FSI simulations containing moving grids, either a Lagrangian or a Eulerian configuration may be applied to the computational domain. For structural problems the conservation equations for solid mechanics are primarily solved in the Lagrangian configuration with respect to a moving reference frame. Determination of the corresponding fluid mechanical conservation equations requires consideration of a hybrid approach taking the grid motion into account and, consequently, a blending of the Lagrangian and Eulerian configurations. Therefore the equations for mass and momentum are given in the arbitrary Lagrangian–Eulerian formulation (ALE). Derived from the conservation equations of continuum mechanics, the ALE formulation for mass and momentum of an incompressible medium are as follows: (1) (2) where ρ is the density, is the velocity vector in the fixed reference frame, is the local mesh velocity, σ is the Cauchy stress tensor and is the volume force.
The ALE formulation enables coupling of the fluid and solid mechanics via the Lagrangian configuration on the moving boundary surface S. Using partitioned code coupling, the numerical calculations of fluid and solid mechanics were considered separately, by solving successively the conservation equations of both systems and realizing the interaction by a surface coupling approach providing data exchange via the topologically identical surfaces S of the two computational domains. The required kinematic and dynamic coupling conditions for the data exchange can be formulated as follows.
Kinematic coupling condition – the velocity and position of the coupling boundary surface S must be the same for the fluid and solid mechanics: (3)
Dynamic coupling condition – the product of stress vector acting on the coupling boundary surface S and the corresponding local normalized vector must be the same for the fluid and solid mechanics: (4)
For the structural stiffness and the ratio of fluid density to structural density, ‘addedmass’ effects causing instabilities (Oertel, 2008) can be neglected for aeroelastic problems. Therefore, an explicit serial coupling scheme was applied.
Software and numerical methods
The numerical simulations of the threedimensional flow were performed by the finite volume method in the ALE formulation as implemented in the FLUENT solver. Using the Reynoldsaveraged Navier–Stokes (RANS) approach (Oertel et al., 2008), turbulence viscocity is modelled by the shear–stress transport (SST) kω model (Menter, 1993; Menter et al., 2003). The model combines the robust formulation of the kωmodel of Wilcox (Wilcox, 1988) at regions near the wall with a kϵmodel (Launder and Spalding, 1972) in the outer flow. A basic description of the model can be found in Appendix 1. For the spatial and time discretization, standardized approximation schemes were applied (2nd order upwind and euler implicit schemes, respectively) and the SIMPLE algorithm was used to solve the resulting algebraic equation system. For the structural calculations the finite element solver ABAQUS (SIMULIA Inc., Providence, RI, USA) was used. Nonlinear dynamic equilibrium equations were solved at each time increment by applying implicit integration schemes. The coupling procedure for interaction of the fluid and structure was realized by MpCCI (mesh based parallel code coupling interface, Fraunhofer Institute, Sankt Augustin, Germany) with an explicit serial coupling scheme, in which node positions and relative wall forces are exchanged at the boundary surface S for each coupling time step.
The numerical flow simulations were validated by means of wind tunnel tests (Ruck and Tischmacher, 2010). In this case, the flow behind the mechanical representation of the avian flight model was measured using a TRPIV system. The spatial and temporal development of the macroscopic vortex structures induced by the wing motion was in a good agreement with the results of the FSI simulations. Furthermore it was shown that the numerical simulations underestimated the maximum and minimum velocity values a little causing a reduction of the calculated aerodynamic coefficients (Oertel, 2008). Small scale flow structures which were additionally resolved by the experiments showed marginal influence on macroscopic vortices. Note, small scale flow structures were filtered by the mesh size; a much finer grid connected to anisotropic turbulence models would be needed to resolve a wider flow and turbulence spectrum.
RESULTS
The simulations were performed for three different flow velocities, u=3, 6 and 9 m s^{–1} and three wingbeat frequencies, f=8, 10 and 12 Hz. The flow conditions generated correlate with those determined for birds (Spedding et al., 2003; Rosén et al., 2004; Tobalske and Dial, 1996; Tobalske et al., 2003) and bats (Hedenström et al., 2007; Tian et al., 2006). The corresponding reduced frequencies range from k_{min}=0.22 to k_{max}=1.00 and the Strouhal numbers^{c} based on the wing amplitude range from St_{min}=0.12 to St_{max}=0.54. The aerodynamic quantities of the entire fluid domain and the corresponding mechanical quantities of the wings were recorded at a write frequency of f_{w}=200 Hz (for f=8 and 10 Hz) and f_{w}=330 Hz (for f=12 Hz). Furthermore, lift and drag forces were recorded at a frequency of f_{w}=1 kHz. For postprocessing one steadystate wingbeat cycle was considered providing data of a fully developed flow field. Hence, 21 (10 Hz), 26 (8 Hz) and 29 (12 Hz) solved data sets were obtained.
Wake analysis and topology
Reconstruction of the vortex configuration was performed for the model at three reduced frequencies. In this case, the flow was visualized on the isosurface of λ_{2} (Jeong and Hussain, 1993) to highlight vortex structures by displaying regions of pressure minima. The vortices were superimposed with the velocity magnitude. Figs 4 and 5 show the results for three conditions at the midupstroke and middownstroke wing positions, respectively. In addition, a schematic depiction was derived by analyzing the vorticity distribution, vortex cores and streamlines plotted on the vortex surfaces. The spatial resolution of the analysis was restricted by the cell size of the grid.
Considering the wake pattern over a wingbeat cycle at high reduced frequencies (k=0.67, f=8 Hz and u=3 m s^{–1}), the wake is dominated by a chain of separated vortices behind each wing. The time and spatial development of the wake is strongly influenced by the shedding of starting and stopping vortices which occurs at the reversal points. The vortices interlock with each other forming a continuous wake as shown in Fig. 4A and Fig. 5A. At the beginning of the downstroke, a clockwise rotating bound vortex is produced around each wing. Because of the wellknown vortexforming mechanism of classical airfoil theory (Oertel, 2009), corresponding anticlockwise rotating vortices (starting vortex) are generated and shed to the wake. Instead of uniform spanwise shedding of a compact starting vortex, separation begins by the formation a vortex sheet at the upper reversal point of the wingbeat cycle. The sheet extends widely over the wing and consists of the starting vortex structure and clockwise and anticlockwise rotating small scale vortices which are shed continuously into the wake. Separation of the vortex sheet starts at the wing tip and shifts to the wing base while moving downwards. The generated vortex structure is stretched streamwise, generating an elliptical vortex tube. The bound and starting vortex are connected via vortices at the tip (tip vortex) and base (root vortex) of the wing. After the midposition of the downstroke, the tip vortex is compact and well defined, whereas the root vortex is still overlaid by small scale structures. At the end of the downstroke, during the maximum wing pitching velocity, the bound vortex (stopping vortex) rolls to the wake, forming a closed vortex ring (DVS – downstroke vortex structure) behind each wing (see Fig. 5A). Simultaneously, the upstroke begins and a clockwise rotating starting vortex is generated behind each wing. This is interlaced with the synchronously rotating stopping vortex of the downstroke and is shed to the wake during the upstroke. As with the vortexforming mechanism of the downstroke, separation of the starting vortex is accompanied by the shedding of a vortex sheet. The sheet declines and starts to separate at the wing tip and is fully disconnected at the end of the upstroke. The vortices at the end and base of the wing change their rotational direction after passing the lower reversal point and generate a connection between the starting vortex and bound vortex of the upstroke. The bound vortex is not fully developed at the wing base and its rotational direction is indifferent. Therefore, the root vortex and the bound vortex close to the body are weak and have already begun to diffuse while the wing is moving upwards. At the end of the upstroke the wing pitches anticlockwise and the bound vortex (stopping vortex) moves to the wake and detaches completely, forming an elliptical semicircular vortex tube (UVS – upstroke vortex structure) behind each wing. The structure of the stopping vortex vanishes quickly as it moves further downstream. The starting vortices of the upstroke and the stopping vortices of the downstroke connect spanwise with each other and form heartshaped vortex loops while the connection to the root vortices vanishes. The wellpronounced tip vortex and the residual structures of the stopping vortex of the UVS are interlocked with the simultaneously shedding and synchronously rotating starting vortex of the DVS.
The vortex structures move downstream with a nearly constant flow velocity of u=3 m s^{–1}. Because of the flapping motion, the vortices are aligned transversely in the wake and the streamwise wavelength of the DVS is λ_{DVS}=0.22 m=0.57u/f (see Fig. 6). As illustrated in Fig. 7A the DVS and UVS are interlocked via their synchronously rotating stopping and starting vortices, which are located on top of each other. At the transition zones of these vortices an opposite momentum exchange occurs, which causes the vortex tubes to disperse crosswise leading to an expansion of the wake. The DVS moves downwards and the UVS upwards. Therefore, the welldefined vortex structures are more exposed to the free stream and their diffusion is amplified. The connection between the stopping vortex of the UVS and the starting vortex of the DVS diffuses rapidly and has already vanished one wavelength downstream. In contrast to the fast diffusion of the UVS, the DVS is sustained further downstream and starts its decay with the decrease of the stopping vortex.
At medium reduced frequencies (k=0.34, f=8 Hz and u=6 m s^{–1}) the vortex generation mechanism is similar to that of high reduced frequencies during the downstroke but differs during the upstroke. At the beginning of the downstroke, an anticlockwise rotating starting vortex superimposed by small scale vortex structures sheds to the wake behind each wing and generates a vortex sheet (see Fig. 4B, Fig. 5B). The vortex sheet consists of anticlockwise and clockwise rotating small scale structures and is connected to the root vortex. During the downward wing motion, the vortex sheet starts to separate at the wing tip and is fully separated at the end of the downstroke. The starting, tip and root vortices form a vortex tube, whereas the clockwise rotating bound vortex is still attached to the wing and has not been shed. During the upstroke, detaching small scale vortex structures of the bound vortex generate a further vortex sheet at the outer part of the wing. This is connected to the tip vortex and yields a continuous and wellpronounced tip vortex structure. The vortex sheet expands to the wing base and thus the root vortex starts to decrease and moves to the rear of the body. Simultaneously, the vortex sheet starts to separate at the wing tip and remains attached for longer at the inner part of the wing, leading to a partial detachment of the bound vortex (stopping vortex) throughout the upstroke. On passing the midwing position, weak oppositely rotating tip and root vortices develop and are shed to the wake. At the upper reversal point of the wingbeat cycle, the bound vortex (stopping vortex) and its corresponding vortex sheet are completely detached from the wing, generating a closed vortex structure (DVS) with the starting, tip and root vortices of the downstroke. During wing pitching an anticlockwise rotating starting vortex is generated and the downstroke begins again. The DVS is dominated by its large and welldefined starting and stopping vortices which are merged partially with its corresponding vortex sheets and, thus, develop nonhomogeneously. Because the bound vortex remains attached for longer throughout the upward wing motion, the UVS is only observed at the end of the upstroke and is represented weakly by its weak tip and root vortices. These vortices are connected to the starting vortex of the DVS. While moving downstream with a nearly constant flow velocity of u=6 m s^{–1}, the tip vortex is amplified, whereas the root vortex diffuses quickly. Because of the flow velocity, streamwise vortex stretching is induced and the DVS is more elliptical leading to a DVS wavelength of λ_{DVS}=0.63 m=0.84u/f. The stopping vortices are connected spanwise with each other providing heartshaped vortex loops. One wavelength further downstream, the connection between the affected root and stopping vortices vanishes. As with higher reduced frequencies, the vortices are aligned transversely in the wake. However, as a result of the weakly defined UVS, opposite momentum exchange is almost completely inhibited and the crosswise expansion of the wake is weak.
The tip and root vortices increasingly dominate the wake pattern at lower reduced frequencies (k=0.22, f=8 Hz and u=9 m s^{–1}) (see Fig. 4C, Fig. 5C). As with medium reduced frequencies, spanwise vortex structures are generated at the upper reversal point, providing one closed vortex structure (DVS) per wingbeat cycle. At the beginning of the downstroke, an anticlockwise rotating starting vortex is shed to the wake and the corresponding vortex sheet of small vortex structures occurs. Small scale structures detaching continuously from the upper and lower sides of the wing surface rotate anticlockwise and clockwise, respectively. The vortex sheet separates at the wing tip and is present for longer at the base of the wing causing a continuous and wellpronounced root vortex throughout the downstroke. At the end of the downstroke, during the maximum wing pitching, the vortex sheet is detached and small scale structures related to the bound vortex start to shed at the outer part of the wing. These small scale structures form a further vortex sheet which stretches out over the wing during the upstroke and is interlocked with the tip vortex, retaining its rotational direction. Weak oppositely rotating tip and root vortices develop while passing the mid wing position and diffuse quickly as they move downstream. The root vortex starts to separate from the wing before the upper reversal point is reached. At the upper reversal point, the clockwise rotating bound vortex (stopping vortex) starts to shed at the outer part of the wing and connects with the already detached root vortex further downstream, producing a closed vortex structure (DVS). At the wing base the vortex sheet remains attached for longer and forms a connection to the anticlockwise rotating starting vortex of the downstroke which is shed during the wing pitching. As a result of the almost nonexistent UVS, the DVS covers nearly the entire wake with a streamwise wavelength of λ_{DVS}=1.01 m=0.896u/f. The continuous shedding of the small scale vortex structures and the omnipresent bound vortex during the downstroke and upstroke lead to a decrease of starting and stopping vortices and thus the wake is dominated by the tip and root vortices. As with higher reduced frequencies, stopping vortices are connected spanwise with each other and the connection between the induced root and stopping vortices vanishes as the root vortices diffuse. Because of the horizontal alignment of the stopping and starting vortices and the almost nonexistence of the UVS, the crosswise expansion does not occur.
Considering the wake behind the body, secondary vortex structures can be observed which are produced as a result of flow separation around the body. Those structures diffuse quickly further downstream. Compared with the strength of the vortices generated by the wing motion, their effects on the development of the wake configuration can be neglected.
Lift and drag development
The wingbeat motion indicates continuously alternating pressure and friction forces on the wing throughout the wingbeat. Therefore, the drag/thrust occurring and the lift production which depends on the flow velocity and flapping frequency, also change. The drag and lift coefficients^{d} with highlighted maxima and minima and the corresponding time integrated values over one wingbeat cycle are illustrated in Figs 8 and 9. The polar diagrams for the recorded drag/thrust and lift forces of the wings are shown in Fig 10. In addition, characteristic wing positions are marked.
The maximum lift coefficient c_{l,max}=7.136 and minimum lift coefficient c_{l,min}=–5.596 are produced for a flow velocity of u=3 m s^{–1} and a wingbeat frequency of f=12 Hz after passing the midstroke position. The corresponding dimensionless times are t/t_{0}(c_{l,min})=0.048 and t/t_{0}(c_{l,max})=0.540. At these reduced frequencies, a huge vertical momentum force is induced, as a result of the high flapping frequency and, consequently, the aerodynamic performance is dominated by the wing motion and nearly independent of the flow velocity. In this case, lift and downforce are generated to a comparable degree throughout the upstroke and downstroke. For decreasing reduced frequency, the lift coefficient decreases. Because the bound vortices (stopping vortices of the downstroke) remain attached for longer, the upstroke becomes aerodynamically active and lift is generated continuously, whereas the downforce decreases. The smallest maximum lift coefficient c_{l}=0.664 is generated for k=0.22. In each case, the maximum and minimum lift coefficients shift increasingly from the midstroke position to the reversal point (see Fig. 8). The decreasing minimum lift coefficient converges to c_{l}=0 for decreasing reduced frequency and shifts to the upper reversal point.
The maximum negative drag (maximum thrust) coefficient is produced for a flow velocity of u=3 m s^{–1} and a wingbeat frequency of f=12 Hz after passing the middownstroke position. Because of the airtight wings and the synchronous beat motion, negative drag coefficients are also generated during the upstroke for high reduced frequencies of k>0.34 (see Fig. 8). For increasing flow velocity or/and decreasing flapping frequency, the horizontal momentum exchange which is produced by the wing motion decreases and thus the thrust generation suffers. In this case, the maximum negative drag coefficient shifts increasingly to the reversal points of the wingbeat. The maximum drag coefficient correlates with the shedding of the bound vortex which occurs during the upstroke. Because the bound vortices remain attached for longer for decreasing reduced frequency, the maximum drag coefficient shifts to the midupstroke position. Maximum drag is generated for a flow velocity of u=9 m s^{–1} and a wingbeat frequency of f=8 Hz at the lower reversal point.
The development of drag and lift for one wingbeat cycle depends on the reduced frequency. For decreasing beat frequency the timeintegrated drag increases and the lift decreases (see Fig. 9). Compared with real birds, the lift force determined is small, ranging from F_{l,min}=0.055 N (f=8 Hz, u=3 m s^{–1}) to F_{l,max}=0.507 N (f=12 Hz, u=9 m s^{–1}). However, for all reduced frequencies, lift is positive and the corresponding coefficients range between c_{l}=0.290 (f=8 Hz, u=6 m s^{–1}) and c_{l}=0.477 (f=12 Hz, u=3 m s^{–1}). For reduced frequencies of k>0.335, thrust is generated. The timeintegrated drag forces range from F_{d,min}=–0.072 N (f=10 Hz, u=6 m s^{–1}) to F_{d,max}=–0.285 N (f=12 Hz, u=3 m s^{–1}) with c_{d}=–0.106 to c_{d}=–1.675, respectively.
The relationship between lift and drag coefficients can be seen in polar diagrams (Fig. 10). For a flow velocity of u=3 m s^{–1} thrust is produced throughout the wingbeat cycle, and lift and downforce are generated during both the downstroke and upstroke (see Fig. 10A). Because of the reduced influence of the wing motion on the horizontal momentum exchange for increasing flow velocity, the ratio of drag to thrust increases. The value of the corresponding force coefficient decreases. For reduced frequencies of k<0.50, thrust is generated only during the downstroke (see Fig. 10B). The strong development and longer existence of the bound vortices lead to an increase in the ratio of lift to downforce. In the range k=0.22–0.34 downforce is produced close to the upper reversal point.
Fig. 11B shows the magnitude and components of the relative movement of the reference point RPO (see Fig. 11A), over one wingbeat cycle at a reduced frequency of k=0.67 (f=8 Hz, u=3 ms^{–1}). The deformation of the flexible airfoil relative to the rigid airfoil is illustrated in Fig. 11C for four different times. It can be seen that the curvature of the airfoil decreases leading to a reduction of lift. Therefore, the lift obtained for the FSI simulations is smaller than that for the corresponding flow simulations with rigid wings. A direct proportionality between the difference and the reduced frequencies is not inferred from the results. However, for a given flapping frequency, the difference in lift production between the FSI simulation and the flow simulation decreases for increasing flow velocities, and for a given flow velocity, the difference increases for increasing flapping frequencies. The difference in the timeintegrated lift coefficient is in the range 5.0% (f=8 Hz, u=3 m s^{–1}) to 17.4% (f=12 Hz, u=3 m s^{–1}). The drag coefficients also depend on the wing deformation. The difference in the drag coefficient between the FSI simulation and the flow simulation demonstrates no regularity and is in the range 0.2% (f=8 Hz, u=9 m s^{–1}) to 15.4% (f=12 Hz, u=9 m s^{–1}).
In contrast to real flying animals, the deformation of the model wings is dominated by the timedependent inertia forces and the contribution of aeroelastics is small. Therefore, the deformation of the model wing depends on flapping frequency and wing weight, while the influence of the flow velocity can be neglected. The largest relative motion of the trailing edge occurs in the direction of the stroke direction (zdirection).
DISCUSSION
The present results provide a macroscopic qualitative description of the spatial and timedependent generation and development of vortex structures behind flapping wings. The vortex wake behind the avian model at high reduced frequencies (k=0.67) differs from the wake models presented for birds. As a result of the omnipresent aerodynamic upstroke of the model, the wellknown succession of discrete vortex rings in the wake at low flight speeds (the vortex ring gait), which have been postulated (Rayner, 1979a; Rayner, 1979b) and experimentally verified (Spedding et al., 1984; Spedding, 1986), was not found. The combined shedding of starting and stopping vortices at the reversal points is not comparable to that found in the doubleringed models of thrust nightingales (Spedding et al., 2003) and robins (Hedenström, 2006). These differences are based on the simplified wing motion of the avian flight model during the upstroke and its airtight and solid wing design. However, it can be assumed that for birds using an almost inflexible wing motion, for acceleration in the forward direction of flight at medium and fast flight speeds, the present results show a plausible wake model. In this case, the increased stroke velocity of the wings causes an alternating circulation of the bound vortices during the upstroke and downstroke and leads to an enhanced thrust production whereas the lift generation suffers. As described, the corresponding vortex structures shed at the reversal point of the wingbeat cycle generate another type of ladder wake (Pennycuick, 1988; Hedrick, 2002).
The occurrence of root vortices close to the body at the wing base leads to a closed vortex behind each wing and has been described for slowflying bats (Muijres et al., 2008; Hedenström et al., 2009) and an insect model of a hovering hawkmoth (Van den Berg and Ellington, 1997a; Van den Berg and Ellington, 1997b). As with previous aerodynamic studies, the bound vortices of the downstroke remain attached to the wing for longer for decreasing reduced frequency. Consequently, spanwise vortex shedding occurs, leading to diffusion of the starting and stopping vortices and a reduction of the upstroke vortex structure. The transition of the wake configuration forms intermediate vortex structures (Spedding, 1987) which do not belong to one of the two idealized vortex wakes (Rayner and Gordon, 1998). Considering the wake models of bats (Johansson et al., 2008; Hedenström et al., 2009) for flow velocities around 4 m s^{–1}, the wake configuration seems to be consistent with the present results at reduced frequencies of k=0.34. The diffused upstroke vortex structure correlates with the oppositely rotating vortices observed for bats at the end of the upstroke. Furthermore, the strong root vortices which are shed at the wing base accord with the root vortices of bats, but the spanwise ring structures are not the same. For decreasing reduced frequency an almost continuous spanwise vortex shedding occurs. The wake structures are stretched streamwise and the corresponding tip and root vortices are well pronounced, whereas the starting and stopping vortices diffuse. The wake structure transforms more and more to a continuous vortex gait type (Rayner, 1986; Rayner and Gordon, 1998; Spedding, 1987), but with a doubleringed configuration.
Although the kinematics and geometrical entities of real birds and the avian model are not identical, some fundamental aerodynamics can be compared and applied to one another. In this case, a basic question concerning root vortex structures arises: do root vortices occur during flapping flight of birds and form vortex rings behind each wing? On the basis of results of flow visualization tests with flying birds and the present simulation results, the appearance of such root vortices cannot be excluded.
The generation of root vortex structures has not been verified until now. On the one hand, the formation of root vortices is amplified by a continuous pressure exchange at the gap between the wing and the body which does not exist for birds. Furthermore, the huge body and the mounting suspension disturb the flow excessively and represent an additional source of errors. On the other hand, the kωSST turbulence model overpredicts the extent of the separated region (Menter et al., 2003). However, the single vortex ring configuration requires a wing motion synchronized circulation around the body to ensure the consistent spanwise circulation of the shed vortices. Furthermore, small bodies are needed to prevent additional vortex separation. Neither of these conditions has been clearly verified for flying birds. Therefore, the vortex formation across the wingspan can be interrupted and weak root vortices are shed and may lead to the generation of closed vortex rings behind each wing.
Lift and drag forces at high reduced frequencies have been computed using FSI simulation techniques. The lift and drag production correlates with the airtight wing construction of the avian model and its kinematics. Therefore, the coefficient diagrams show the expected curves for one steadystate wingbeat cycle depending on the flapping and pitching angle. It was shown that, for decreasing reduced frequency, the bound vortex of the downstroke remains attached to the wing for longer, and consequently the lift generation is enhanced. As with the lift production of bats (Hedenström et al., 2007), the aerodynamically active upstroke leads to the generation of a downforce at the end of the upstroke when the bound vortex is shed into the wake. In this case, the momentum exchange induced by the wing motion suffers and the thrust decreases.
For all flow velocities and flapping frequencies, lift is produced and thrust appears for reduced frequencies of k>0.34. The ratio of lift to thrust (negative drag) F_{l}/F_{d} ranges from (F_{l}/F_{d})_{min}=0.369 (k=1.00) to (F_{l}/F_{d})_{max}=3.300 (k=0.42) and is small compared with that of real birds, e.g. for a pigeon F_{l}/F_{d}=9 (Dial and Biewener, 1993) and for a swift F_{l}/F_{d}=11 (Henningsson et al., 2008).
Concluding remarks
For the first time threedimensional unsteady FSI simulations have been conducted for an avian flight model at a wide range of flapping frequencies and high Reynolds numbers. The results give an improved insight into the complex vortex configurations behind flapping wings, helping us to understand the threedimensional and timedependent formation mechanism. The wakes consist of interlocked elliptical vortex rings on each wing. For high reduced frequencies, an upstroke and a downstroke vortex structure were produced which contain starting, stopping, root and tip vortices. It was shown that, for decreasing reduced frequency, the upstroke vortex structure diffuses and the tip and root vortices of the downstroke dominate the wake. The vortex configurations occurring are not consistent with those of two idealized wake configurations postulated for birds (Rayner, 1979a; Rayner, 1979b) and differ from those recorded for birds (Spedding et al., 2003; Rosén et al., 2004; Hedenström, 2006). The analysed vortex structures are related to the wake of birds with an inflexible wing motion and may occur during acceleration in the forward direction at medium and fast flight speeds. Furthermore, the results show similarities to those of bats (Muijres et al., 2008; Hedenström et al., 2009) and an insect model (Van den Berg and Ellington, 1997a; Van den Berg and Ellington, 1997b). Although the model does not represent real bird geometry and wingbeat kinematics, flow characteristics that may appear in the aerodynamic performance of birds are described. The nonappearance of the root vortex during bird flight has been discussed, but the question has not yet been answered. Further numerical and experimental studies will be necessary using avian models without a gap and with real bird geometry and structural wing behaviour. In this case, numerical representation of a real bird wing for FSI simulations will be desirable to predict the contribution of the FSI to lift and drag production.
APPENDIX 1
The kωSST turbulence model
In the present paper the RANS equations were used to solve the turbulent flow field. The conservation equations for mass and momentum for the timeaveraged continuity equation and the RANS equations given by: (A1) (A2) where τ_{t}_{,}_{ij} is turbulent Reynolds stress, u_{i} is timeaveraged flow velocity in the idirection, u_{j} is timeaveraged flow velocity in the jdirection, ρ is density, is timeaveraged volume force, is timeaveraged viscous stress and is timeaveraged pressure.
Using the Boussinesque approximation, the Reynolds stresses can be related to the mean velocity gradients: (A3) where u_{i}′ is velocity fluctuation in the idirection, u_{j}′ is velocity fluctuation in the jdirection, μ_{t} is turbulent viscosity, k is turbulent energy and δ_{ij} is the Kronecker delta.
The turbulent viscosity is modelled using the two transport equations of the kωSST model which was developed by Menter (Menter 1993; Menter, 1994; Menter et al., 2003). It combines the robust formulation of the kωmodel of Wilcox (Wilcox, 1988) at near wall regions with a kϵmodel (Launder and Spalding, 1972) at the outer flow, transforming the kϵmodel to the kωformulation. The original kωequations are multiplied by the blending function F_{1} and the transformed kϵmodel is multiplied by (1–F_{1}). The blending function is F_{1}=0 at the outer flow and F_{1}=1 at the wall. The corresponding equations are added together to give the transport equations for the kωSST model: (A4) (A5) The definition of the functions and model constants may be found in the literature (Menter et al., 2003) (FLUENT 6.3 User Guide 2006). In the present paper, the standard settings of FLUENT were used for the constants. Note that the character of the turbulence is not fully resolved and the accuracy suffers when the RANS approach is applied. The accuracy of the kωSST turbulence model is comparable to that of other twoequation turbulence models and is in acceptable agreement with experimental data in terms of force and flow development (Ruck and Tischmacher, 2010; Rival et al., 2007; Rival et al., 2008; Menter et al., 2003; Oertel et al., 2008).
APPENDIX 2
Statistics
By means of Eqns A5–A9, the statistics of the grid sensitivity study were determined (Tables A1 and A2).
Maximum standard deviation SD_{max}: (A6) Mean force coefficient m_{ma}_{x} of SD_{max}: (A7) Maximum standard error of the mean SEM_{max}: (A8) Timeaveraged standard error of the mean SEM_{ta}: (A9) where c_{ij} is the force coefficient for mesh i=1,..., N at time step j=1,..., n and is the mean force coefficient at time step j=1,..., n.
FOOTNOTES

↵a The Reynolds number belongs to the class of dimensionless numbers in fluid mechanics. It is used to characterize the ratio of inertial force to viscous force in the flow. The Reynolds number for wings can be calculated using the mean wing chord c, the undisturbed flow velocity u and the kinematic viscosity ν: Re=uc/ν.

↵b Analogous to the Strouhal number St, the reduced frequency k is a measure of the unsteadiness behind flapping wings. It is used to describe the ratio of local fluid acceleration to inertial force. Using the mean wing chord c and the wing flapping amplitude a, the undisturbed flow velocity u and the flapping frequency f: k=cfπ/u and St=af/u.

↵c See footnote b.

↵d The nondimensional lift and drag coefficients describe the aerodynamic performance of flow profiles and are defined as follows: c_{l}=F_{l}/(qA) and c_{d}=F_{d}/(qA) where F_{l} is lift force, F_{d} is drag force, q is dynamic pressure and A is wing area.

LIST OF SYMBOLS AND ABBREVIATIONS
 a
 flapping amplitude
 A
 wing area
 ALE
 arbitrary Lagrangian–Eulerian
 AR
 aspect ratio
 ASF
 cell average shape factor
 b
 wing span
 bf
 body frontal area
 bg
 body ground area
 bl
 body length
 bt
 body thickness
 c
 mean chord
 c_{d}
 drag coefficient
 c_{l}
 lift coefficient
 d
 deformation
 DVS
 downstroke vortex structure
 E
 Young's modulus
 f
 flapping frequency
 f_{i}
 timeaveraged volume force
 f_{w}
 write frequency
 F_{d}
 drag force
 F_{l}
 lift force
 FSI
 fluid–structure interaction
 k
 reduced frequency
 k
 turbulent energy
 m_{max}
 mean force coefficient
 normalized vector
 p
 pressure
 p_{∞}
 ambient pressure
 PIV
 particle image velocimetry
 q
 dynamic pressure
 RANS
 Reynoldsaveraged Navier–Stokes
 Re
 Reynolds number
 RPL
 reference plane
 RPO
 reference point
 s
 semispan
 S
 boundary surface
 SD_{max}
 maximum standard deviation
 SEM_{max}
 maximum standard error of the mean
 SEM_{ta}
 timeaveraged standard error of the mean
 SST
 shear–stress transport
 St
 Strouhal number
 t
 time
 t/t_{0}
 dimensionless time
 u
 undisturbed flow velocity
 u_{i}′
 velocity fluctuation in the idirection
 u_{j}′
 velocity fluctuation in the jdirection
 U_{i}
 timeaveraged flow velocity in the idirection
 U_{j}
 timeaveraged flow velocity in the jdirection
 UVS
 upstroke vortex structure
 velocity
 local mesh velocity
 x,y,z
 coordinates
 γ
 flapping angle
 δ_{ij}
 Kronecker delta
 Δp
 pressure difference
 θ
 angle of attack
 λ
 wavelength
 λ_{D}
 DVS wavelength
 μ_{i}
 turbulent viscosity
 ν
 kinematic viscosity
 ρ
 density
 ρ_{m}
 material density
 σ
 Cauchy stress tensor
 stress vector
 timeaveraged viscous stress
 τ_{t,ij}
 turbulent Reynolds stress
 ω
 angular frequency
 © 2010.