Generation of the vorticity field by the flapping profile

Birds and insects move due to flapping their wings. In the paper we presented unsteady effects that led to the lift and thrust force generation on the flapping foil. It was assumed that the flow was two-dimensional. We demonstrated that flapping motion with low Reynolds produced well-ordered vorticity around the profile. The well-ordered vorticity atmosphere generates predictable aerodynamic force distribution on the profile. We demonstrated that high enough parameters of flapping motion caused the disordered vorticity distribution around the profiles that led to the loss of the lift force and limited the possibility of a flight.


Introduction
In contrast to classical airplanes, most insects and birds operate in low Reynolds number regime, below, 10 4 [1]. Such a flow is dominated by unsteady effects caused by vorticity dynamics. For low Reynolds number flows around the flapping foil one can observe the ordered fluid structures [2,3]. The ordered fluid structures behind the flapping profile create vortex atmosphere which can be divided into several types of vortex streets [6]. The type of vortex wake is related to aerodynamic force distribution on the profile [7]. Proper vortex wake generation is extremely important, especially for small insects which fly in very low Reynolds number flows (up to 500) [1,5]. In the paper we demonstrated three different flapping kinematics. First, we presented the formation of the vortex street and its transitions (into different types of vortex wake) for a simple plunging profile. We observed the whole family of vortex streets and created diagrams which linked the type of the observed vorticity field with the Reynolds, the Strouhal number and the plunging amplitude. The main purpose of our simulations was to investigate the flow phenomena related to transitions of the vortex wake and how this transition depends on the Reynolds number. The second part of our investigation was related to the simulation of flapping motion which is appropriate for small insects. It corresponds to the plunging and pitching profile [8]. We presented the dependence of aerodynamic forces calculated on the profile on the Reynolds number and different types of vorticity field. In the end, we simulated the free flapping motion. We demonstrated various profile trajectories for different flapping kinematics and different types of vortex wake behind the moving profile. We limited our research to two-dimensional flows. Although a flight in nature is three-dimensional, a two-dimensional simplification is widely used and seems to be appropriate to capture the essence of flapping flight aerodynamics [1,4,5,6]. For the study of a moving object in fluid we chose the vortex particle method -Vortex-in-Cell method. The importance of the VIC method lies in the possibility to analyze more directly the vorticity field since, in computations, particles that carry information about the vorticity field are used. The method was carefully tested for problems with known analytical solutions and well-documented experimental data.

Vortex particle method
The Navier-Stokes equations of the incompressible viscous fluid in primitive variables (u,p) can be transformed to the Helmholtz equations in vorticity-stream function formulation. The Helmholtz equations in the non-inertial reference frame have the form where = ( , ) means vector velocity, stream function and viscosity. For our purposes, it was convenient to change the reference frame. The motion of the profile is made by the boundary condition for the stream function far from the body with specified angular  and translational U 0 velocity.
where  means instantaneous angle of attack of the moving object in fluid and U 0 = |u|. A detailed description of the solution of the Helmholtz equations in the moving reference frame can be found in [9]. Numerical calculations of the moving profile were carried out by the vortex in cell method (VIC). In the calculation we used the numerical grid. To better fit the numerical grid to the solid boundary we used the conformal transformation. Through the conformal mapping the physical non-uniform flow region (x,y) is replaced by the rectangular one (), where the calculations were performed very efficiently. The following conformal transformation was applied ( Fig. 1) (4) In new variables the equations of motion have the form where J denotes the Jacobian of the transformation.
In the vortex in cell method the continuos vorticity field is approximated with the discrete particle distribution. The calculation area is covered with the uniform grid h = and in every grid node the particle with circulation: Γ = ∫ is placed, where A = h 2 and The equations of motion were solved by means of viscous splitting algorithm [10]. First, the inviscid fluid motion equation was solved The equation (8) proves that the vorticity is constant along the trajectory of the fluid particles. According to the Helmholtz theorem, vortex particles move in the same manner as material fluid particles [7,10], partial differential equation is replaced by the set of ordinary differential equations which describe the trajectory of the particles where     denotes the Lagrangian coordinate of fluid particles. The number of particles equals the number of grid nodes. The set of ordinary equations (6) was solved by the four order Runge-Kutta method. The velocity field was obtained through the solution of Poisson's equation for the stream function and vorticity on the numerical grid (5). The velocities of the particles between the grid nodes were calculated by the two-dimensional bilinear interpolation.
In the second stage the viscosity was taken into account by the solution of the diffusion equation where  s denotes the vorticity at the solid boundary,   means the Euler inviscid distribution of vorticity in time step t n . The motion of the body in the fluid and the vorticity in the flow region generate non-zero tangent velocity component at the wall. The no slip condition is realized by putting a proper amount of vorticity at the wall. We used the Briley formula of the fourth order [11]. Nullification of the normal velocity component at the wall was performed by the constant stream function at the wall  const. The diffusion equation was solved by the alternating direction implicit method (ADI). In calculations we used the numerical grid, therefore, it was necessary to transfer the information from the vortex particle to the mesh nodes, Fig. 2. The redistribution of the vorticity from fluid particles to the mesh nodes was performed according to the formula where  h (•) denotes the interpolation kernel. The redistribution process plays a fundamental role in the accuracy of the vortex particle methods. In the current study we used the Z-splines interpolation functions [12], which, when applied to vortex particle methods have very good interpolation properties.

Figure 2.
Redistribution of the particle masses onto the grid nodes, a) for particles lying inside the computational domain, b) for the particles in the vicinity of the wall.

Plunging simulation
We performed a numerical calculation for the simple plunging profile which moves according to the equation where y(t) denotes the instantaneous position of the profile center, A 0 is the amplitude and f denotes the frequency of the plunging. The profile's vertical velocity was calculated from the formula: v = y'(t). We assumed that velocity of the fluid U 0 is constant far from the body, see figure 3. The flow over the plunging body can be characterized by the following set of non-dimensional parameters: The Reynolds number Re, the Strouhal number St and the plunging amplitude A C , [1,2]: where c is a chord of the profile and  denotes a kinematic coefficient of fluid viscosity. The profile chord was set to c = 2, we assumed the thickness of the profile: e = b/c = 0.2. We performed numerical calculations for the Reynolds number Re = 100, Re = 250 and Re = 500, with the density: ρ = 1. The plunging frequency was fixed to f = 0.5 and the Strouhal number (St) was controlled by changing the free stream velocity U 0 . The calculations were performed for dimensionless time T = ft, in the range of T = (0, 10), which corresponds to ten periods in the equation (15). We used the elliptical grid shown in figure 1, with 256 grid nodes in a radial direction and 256 grid nodes in an azimuthal direction. In figure 4, the visualization of vortex wake around the profile was presented. If the amplitude and the Strouhal number were low, we observed a vortex Karman street (figure 1a).  We can summarize the above results in a phase transition diagram. In figure 5 we presented the relation between the vortex wake type and the plunging parameters. The phase diagram for the Reynolds number Re = 500 can be found on the right. The increase in the plunging parameters caused vanishing of the stable deflected Karman vortex street. For these flapping parameters the deflection of the vortex wake was constantly changing. A similar effect for flapping with high Strouhal number: St ~ 1 was observed in [13]. We showed in [6] that the effect which is responsible for changing the deflection is the boundary layer eruption phenomenon. A description of the discussed phenomenon can be found in [14].
The first part of our research was related to simple plunging motion. We found that for the low Reynolds number flows, flapping kinematics is very important to create the ordered vorticity field around the profile. The increase of the Reynolds number or improper flapping parameters caused the disordered vorticity field around the profile.

Pitching and plunging simulation
One of the most popular two-dimensional flapping kinematics is known as the incline hovering [1,5,15]. It may be described by the equations: The visualization of one stroke is presented in figure 6. After the first four strokes the vorticity field was ordered and generated periodically. It can be noticed that the role of the flapping kinematics is to create a vortex pair structure and get rid of it [15]. The vortex pair structure moves due to selfinduction. In the range of the Reynolds number below 125, it causes the removal of vorticity generated in the previous stroke far from the profile. Further, such a vortex pair creates the vortex wake. To check the generation of the vortex pair structures we performed numerical calculations in the wide range of the Reynolds numbers. In figure 7 the comparison of the vorticity field generated around the profile for different Reynolds numbers was presented. We noticed that the creation of the vortex pair did not occur for the Reynolds number below Re = 40 (in the flow over the cylinder the Karman vortex wake created from a vortex  To check the influence of the vorticity field around the profile for aerodynamic forces we calculated the lift force on the profile. In figure 8 we presented the lift force coefficient (22) for different Reynolds numbers. One can observe that for the Reynolds number: Re = 125 the generation of the aerodynamic forces was periodic. If the Reynolds number increased the characteristics of the lift force coefficient were disordered. It is difficult to estimate the difference in aerodynamic force creation from the lift force characteristic. To be able to compare the lift force generation for different Reynolds numbers we calculated the average lift force coefficient The average lift force characteristic is presented in figure 9. Our calculations prove that the highest growth of the average lift force coefficient was in the range of the Reynolds number: 5 < Re <125.
The increase of the Reynolds number from 125 up to Re = 1000 caused only 18% gain of the average lift force coefficient. We believe this is related to the disordered vorticity around the profile with higher Reynolds number flows.

Free flapping flight
To check in what way the disordered vorticity around the profile can affect the trajectory of the real flight we simulated a free flapping flight. The profile was replaced by the mass point which can move according to Newton's second law where m denotes the mass of the profile, F fluid was calculated from the equations (20) and (21), and F denotes the driving force which is necessary to move the profile. The driving force was calculated on the basis of the equations (17), (18). We calculated the mass of the profile assuming the density of the profile: = 18kg/m 3 . The equation (20) was solved by the second order Runge-Kutta method. Parameters of the flapping kinematics were the same as in section 4, however, the frequency of the motion was modified. We noticed that we can simply control the flight by a small modification of the angle  0 in the equation (19).  In figures 10 and 11 the vorticity field generated by the free flapping profile and the instantaneous velocity of the profile were presented. It is well known that all birds and insects have some minimal frequency of flapping motion which is necessary to maintain in the air [16]. If this frequency is too low, we observed the disordered vorticity field and noticed that in such conditions a stable flight is not possible (see figure 10). For proper parameters of flapping motion we can observe that after a few strokes the vorticity field is generated periodically ( figure 11). In figure 12 the trajectory of the flapping profile was presented. We can check how the angle of climb changes with a parameter  0 . Moreover, one can notice that the direction of the vortex wake deflection altered with the increase of the angle of climb (see figure 12a and b). It is worth noting that, although we set the plunging amplitude A 0 = 5, the flapping amplitude was lower due to the resistance of the motion, see figure 12.

Conclusion
In the paper the simulations of the flapping motion for low Reynolds number were presented. Due to a large number of fluid mechanics phenomena related to moving objects in fluid, we limited our research to two-dimensional flows. Despite the two-dimensional simplification the aerodynamic of flapping motion is very rich and full of vortex structures creation, and their interaction. All fluid mechanics phenomena depend strongly on the dynamic of vorticity. We showed how the vortex pair structures generated by the flapping profile form the vortex wake around the profile. For plunging motion of the profile we observed the whole family of different types of vortex wake. We built transition diagrams parameterised with the Strouhal number and the dimensionless plunging amplitude. For appropriate parameters of plunging motion we observed the ordered vorticity field around the profile. We noticed that for higher Reynolds number flows the generation of the stable lift force might be difficult. This being caused by the disordered vorticity field. Additionally, we discovered that the disordered vorticity field around the flapping profile led to the disordered, unpredictable aerodynamic forces. Through the simulation of a free flapping flight we proved that in such circumstances the stable flight may not be possible.