Fluid–structure interaction simulation of an avian flight model

The biomechanics of flying birds and insects plays a crucial role in the domain of aeronautics. In view of the complexity of wing-beat 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 slow-flying birds and was the basis for a well-known 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 slow-flying 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, double-ringed 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). Two-dimensional velocity distributions were used to reconstruct three-dimensional 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.

A three-dimensional 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², 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 half-wing 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^bk, 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 off-centre axis parallel to the body axis provides the flapping motion of the model. The angle of attack was controlled by phase-delayed 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.

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 time-averaged 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₁ was placed 0.3 m downstream from the trailing edge of the wing and point p₂ 0.6 m downstream at the mid-wing 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.

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.

For the structural stiffness and the ratio of fluid density to structural density, ‘added-mass’ effects causing instabilities (Oertel, 2008) can be neglected for aeroelastic problems. Therefore, an explicit serial coupling scheme was applied.

The numerical simulations of the three-dimensional flow were performed by the finite volume method in the ALE formulation as implemented in the FLUENT solver. Using the Reynolds-averaged 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. Non-linear 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 TR-PIV 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.

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 post-processing one steady-state 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.

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 λ₂ (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 mid-upstroke and mid-downstroke 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 well-known vortex-forming 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 mid-position 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 vortex-forming 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 semi-circular 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 heart-shaped vortex loops while the connection to the root vortices vanishes. The well-pronounced 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 well-defined 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 well-pronounced 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 mid-wing 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 well-defined starting and stopping vortices which are merged partially with its corresponding vortex sheets and, thus, develop non-homogeneously. 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 heart-shaped 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 well-pronounced 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 non-existent 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 non-existence 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.

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 mid-stroke position. The corresponding dimensionless times are t/t₀(c_l,min)=0.048 and t/t₀(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 mid-stroke 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 mid-downstroke 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 mid-upstroke 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 time-integrated 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 time-integrated 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 time-integrated 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 time-dependent 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 (z-direction).

The present results provide a macroscopic qualitative description of the spatial and time-dependent 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 well-known 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 double-ringed 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 slow-flying 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 double-ringed 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 steady-state 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).

For the first time three-dimensional 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 three-dimensional and time-dependent 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 non-appearance 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.

By means of Eqns A5–A9, the statistics of the grid sensitivity study were determined (Tables A1 and A2).