Abstract
Force generation in avian and aquatic species is of considerable interest for possible engineering applications. The aim of this work is to highlight the theoretical and physical foundations of a new formulation of the unsteady Kutta condition, which postulates a finite pressure difference at the trailing edge of the foil. The condition, necessary to obtain a unique solution and derived from the unsteady Bernoulli equation, implies that the energy supplied for the wing motion generates trailing-edge vortices and their overall effect, which depends on the motion initial parameters, is a jet of fluid that propels the wing. The postulated pressure difference (the value of which should be experimentally obtained) models the trailing-edge velocity difference that generates the thrust-producing jet. Although the average thrust values computed by the proposed method are comparable to those calculated by assuming null pressure difference at the trailing edge, the latter (commonly used) approach is less physically meaningful than the present one, as there is a singularity at the foil trailing edge. Additionally, in biological applications, that is, for autonomous flapping, the differences ought to be more significant, as the corresponding energy requirements should be substantially altered, compared to the studied oscillatory motions.
1. Introduction
Generating force by active flapping of fins and wings has become a focus of research in multiple disciplines from zoology through fluid mechanics and mathematics. Strong interest in the field started back in the 1950s with the landmark studies of Weis-Fogh culminating in the paper “Unusual Mechanisms for the Generation of Lift in Flying Animals” [1]. A relevant theoretical model was derived by Lighthill [2] and other seminal works include those of Ellington [3], on insect flight, and Pedley and Hill [4], on fish swimming. The idea of utilizing the forces generated by flapping wings for the propulsion of man-made vehicles emerged consequently from zoological observations and numerous recent works address the related issue of propulsive efficiency; see, for example, Eloy [5] and references therein. The thorough reviews by Rozhdestvensky and Ryzhov [6] and Triantafyllou et al. [7] focus primarily on experimental studies of flapping-wing propulsors. Key areas of interest, for example, flow unsteadiness (already crucial in [1, 2]), wing flexibility, and three-dimensionality, were identified by the authors who, additionally, highlighted that a lot more is to be learned about how moving wings function as well as about the role of the relevant design parameters. A further in-depth review is that by Platzer et al. [8] who specifically discussed the biomimetic design of micro air vehicles utilizing flapping propulsion; see also Jones et al. [9].
However, experiments on their own are not sufficient to understand the underlying physical mechanisms as they also rely to some extent on assumptions. Over the years mathematical and computational models were developed and dedicated experiments performed, resulting in a steady flow of high quality scientific works, of which we mention here only a small but relevant sample.
Triantafyllou et al. [10] underlined that the large-scale patterns observed in the wake of oscillating wings are related to the production of a jet-like average flow. Jones et al. [11] investigated both experimentally and numerically the thrust produced by a heaving hydrofoil. The performed experiments provided important data about the unsteady wake formed by the moving wing. The resulting vortical structures were compared with the computational results obtained by a modification of the well-known code developed by Basu and Hancock [12]. Qualitative and quantitative comparisons indicated that the evolution of the wake structures is primarily an inviscid phenomenon, although their formation is ultimately due to viscosity. The latter are crucial findings that justify the use of potential methods in the analysis of thrust generation, providing that suitable assumptions are made, and this computational approach is also employed in the present work, as detailed below.
Basu and Hancock [12] were indeed among the first researchers to devise a panel method for the calculation of the forces acting on a two-dimensional airfoil undergoing an unsteady motion in an inviscid and incompressible fluid. An auxiliary condition, known as the Kutta condition and related to assumptions on the flow characteristics at the trailing edge of the foil, was added to obtain a unique solution. This condition, initially developed solely for steady flows, was proposed in order to ensure that the flow passes the trailing edge smoothly, with a finite velocity (see, e.g., [13] and references therein for further details).
In a broad review of various panel methods Hess [14] remarked that the specification of a proper Kutta condition is more important than any other detail of the numerical implementation. In an earlier work [15] he elucidated that this mathematical treatment of the lifting problem is merely a model that simplifies the physical phenomenon to a potential flow while in reality it is more complex and involves viscosity. In other words, a suitable condition ought to account for the viscosity effects, in order to include them implicitly in a computational model that is essentially inviscid. Poling and Telionis [16] investigated experimentally a number of unsteady flow fields and their results indicate that the classical Kutta condition, which states that the velocity at the trailing edge is finite and the pressure difference there is zero, is not valid in certain conditions; that is, it is a function of the initial parameters. Strong viscous-inviscid interactions of the boundary layer can appear in the trailing-edge region of the moving wing and lead to large streamlines curvatures and nonnegligible pressure gradients. The numerical results of Young and Lai [17] showed that flow separation occurs at the trailing edge of heaving foils meaning that the flow streamlines form a time-dependent trailing-edge vortex rather than smoothly enveloping it on both sides. The edge flow mechanism was independently studied by Liebe [18] who proposed to replace the classical (steady) Kutta condition with a more general condition based on the formation and periodic shedding of trailing-edge vortices.
In summary, there is no experimental evidence that conclusively supports the notion of the pressure difference at the trailing edge being equal to zero for unsteady motions but, despite the lack of such evidence, it has been commonly assumed in the literature (e.g., see again [13] and references therein). Furthermore, the application of linear (potential) theories to flapping flight is backed by the fact that experimental data show clear linear correlations, as reported, for example, by Ol et al. [19] and, more recently, by Mackowski and Williamson [20]. The latter authors stated that the vortex wake behind an oscillating airfoil seems to be a remnant of force generation, rather than a driver, and that viscous effects may have a noticeable influence on the magnitude of the produced forces but do not seem to significantly affect the unsteady vortex dynamics.
The first step in the formulation of a comprehensive unsteady Kutta condition appears consequently to be the relaxation of the postulated null pressure difference at the trailing edge. This allows considering the variation in direction and magnitude of the velocities in the vicinity of the trailing edge, that is, the formation and shedding of trailing-edge vortices, resulting in the generation of the experimentally observed thrust-producing jet. Such a modified condition, postulating a finite pressure difference at the trailing edge of the oscillating wing, was already implemented by us in two-dimensional and three-dimensional panel methods, presented in La Mantia and Dabnichki [21, 22], which were found to produce accurate results in terms of experimental verification and numerical reliability (see also [23] and references therein).
The current work extends the numerical study reported in our previous publication on the influence of the wake model on thrust generation [13]. It is specifically shown below how the choice of the wake initial and boundary conditions affects the calculated thrust, in exemplary two-dimensional cases, different from those discussed previously. The aim of this paper is therefore to reinforce the conclusion of La Mantia and Dabnichki [13], that is, the idea that the wake parameters, necessary to obtain a solution, might be chosen, in potential methods, on a physical basis, which is here represented by the proposed formulation of the unsteady Kutta condition, rather than being just free parameters.
2. Theoretical and Numerical Approach
Theoretically it has been firmly established that the potential flow approach is well suited for the prediction of gross hydrodynamic forces acting on moving bodies (see, e.g., [23] and references therein). We discuss here its particular implementation in the form of an unsteady Boundary Element Method (BEM) computer program. The latter was developed to estimate the hydrodynamic forces generated by oscillating wings and corresponding comparisons with published experimental data showed good agreement with the computational results [21, 22].
Following Katz and Plotkin [24], the fluid is assumed to be incompressible and the flow irrotational. In the (exemplary) two-dimensional formulation discussed here the wing section is represented by a finite number of linear panels and its outer shape is that of a NACA 0012 airfoil (often used in the literature as a sample airfoil shape, see Figure 1). Constant strength distributions of a source and a doublet are situated on each panel, the midpoint of which is called collocation point. The flow potential function at each collocation point is defined as the sum of a local (perturbation) potential , related to the unknown doublet strength, and a free-stream potential , linked to the fluid kinematic velocity. An internal Dirichlet boundary condition is imposed; that is, an inner potential function is specified on the internal foil surface in order to meet the nonpenetration condition. At each collocation point the source strength is known, , where is a unit vector normal to the foil surface pointing into the body and the fluid kinematic velocity due to the imposed foil motion. The governing integral equation is derived by using the Laplace equation and the Green third identity. In the body-fixed coordinate system, at time and for each collocation point, at position , it can be written aswhere and indicate the foil and wake surface (length), respectively, and is the strength of the wake doublet distribution. To define uniquely the problem has to be known or related to the unknown doublets on by means of a suitable condition, that is, the Kutta condition. Besides, changes with time; that is, new portions of the wake surface are added as the time advances, and the wake shape has to be properly modelled. The latter can be seen as an additional assumption that in effect represents initial conditions for the dynamic problem.
The discretized form of (1) can be written aswhere is the number of wake panels at time , that is, at the first time step . Consequently, at each time step, a wake panel is added and its contribution evaluated (each linear wake panel has an assigned chordwise length and a constant strength doublet distribution on it). Moreover,indicates the appropriate two-dimensional doublet influence coefficient of panel at the considered collocation point, anddenotes the corresponding source influence coefficient. They are only dependent on the foil geometry, where is the distance between the panel and the respective collocation point and the length of panel . is defined as ; that is, is the distance between the wake panel and the respective foil collocation point and the length of the wake panel .
As previously stated, the wake has to be modelled. It was indeed chosen to study how the generated forces are affected by the use of two different (typical) modelling approaches (see below the corresponding results). More specifically, the positions of the shed wake panels were obtained either (i) by estimating the influence on each wake panel of the source and doublet distributions of the foil and wake (“Estimated Wake” approach, following [22]) or (ii) by assuming that the wake panels remain static in the inertia coordinate system; that is, the wake follows the foil path (“Assumed Wake” approach, following [21]; also used in the mentioned related study [13]).
Additionally, since the foil is moving, the positions of the wake collocation points in the body-fixed coordinate system have to be calculated at each time step starting from their positions in the inertia frame of reference. The body-fixed position of the wake panels closest to the trailing edge, that is, those added at each time step, was set parallel to the chord.
Besides, in our studies the wake panel length was either (iii) fixed or (iv) set proportional to the time step length (in [13], the latter modelling approach was employed).
To summarize, the numerical results discussed below are obtained, in exemplary two-dimensional cases, as specified, by using options (i) or (ii) for the wake shape and (iii) for the wake panel length. In the earlier study [13] results were obtained using options (ii) and (iv).
Equation (2) represents an algebraic system of equations but the unknowns are since, at time , the doublet strengths of the previously shed wake panels are already derived. This means that, at each time step, an unsteady Kutta condition at the foil trailing edge is needed to solve the system of equations.
The trailing-edge condition, necessary to obtain a unique solution, is derived from the unsteady Bernoulli equation, that is, the conservation of momentum equation for incompressible fluid and irrotational flow. It implies that the second derivative of every involved function, such as the fluid velocity and pressure, exists and is continuous ( space). This is not the case everywhere in the considered domain as the airfoil is a profile with one singularity at the trailing edge. Following Cebeci et al. [25] the unsteady Bernoulli equation can be written for a generic point on the foil aswhere is the fluid pressure far from the oscillating foil, the pressure, the velocity, the kinematic velocity due to the motion of the foil, and the potential function. The fluid velocity is the sum of and the perturbation velocity, which is estimated by means of the derivative of the perturbation potential over the foil. Besides, kg/m3 is the fluid (water) density, which is assumed to be constant and uniform. At the trailing edge the unsteady Bernoulli equation can be written aswhere the subscripts and indicate the collocation points of the lower and upper panels that meet at the trailing edge, respectively. It can be rearranged asto highlight the pressure difference in the vicinity of the trailing edge.
It is important to underline that as the trailing edge is the domain singularity point, the velocity and pressure do not have there a unique value; that is, a pressure difference can be present there. Besides, (7) shows that the unsteady Bernoulli equation on its own is not sufficient to obtain the solution. It is necessary to specify physically relevant boundary conditions, for example, velocity or pressure values. The jump in the potential at the trailing edge, whose time derivative is the third term on the right hand side of (7), is usually assumed to be equivalent to the bound circulation , that is, the line integral of the fluid velocity around the foil. The latter is also represented by , that is, the constant strength of the wake doublet distribution shed at the considered time step [24]. The existence of a singularity at the foil trailing edge ( space) allows such a result. In other words, the jump in the potential at the trailing edge can be generally different from zero and a component of the fluid velocity perpendicular to the wake panel at the trailing edge can then exist.
As already anticipated, it is here postulated that the pressure difference at the trailing edge could be finite rather than zero, for unsteady motion. Following La Mantia and Dabnichki [22], the flow velocity and pressure difference are assumed to be finite at the trailing edge by imposing that the third term on the right hand side of (7) is also finite. This allows considering the variation of the velocities in the vicinity of the trailing edge.
More generally, such an assumption implies that the energy supplied for the wing motion generates time-dependent trailing-edge vortices. Their overall effect, which depends on the motion initial parameters, is a jet of fluid that propels the wing. As the kinetic energy is transferred from the jet back to the wing, the vortices would disappear and it is consequently not necessary to assume the fluid to be viscous.
The postulated pressure difference at the trailing edge is then fundamental for such a model as it can justify the velocity difference that generates the thrust-producing jet. It has to be noted that the trailing-edge vortices could also rotate in the opposite direction; that is, depending on the initial conditions the generated force can propel the wing or oppose its forward motion.
The mentioned pressure gradient (perpendicular to the wake panel at the trailing edge) can be viewed mathematically as a function applied at the domain singularity point, that is, the trailing edge. It could be linked to the energy supply by means of a suitable weight factor that depends on the motion initial parameters. The integral of such a function over the domain would then represent the pressure difference at the trailing edge.
The thrust production can consequently be modelled by taking into account the singularity at the foil trailing edge (which exists as the domain is a space) and the Boundary Element Method allows such an approach. Moreover, the unsteady motion of the foil causes the pressure difference at the trailing edge and the resulting thrust-producing jet. In other words, as discussed below, imposing a null pressure difference at the domain singularity point might lead to computational results that agree well with experimental data but do not actually account for the physics of the flapping-wing problem.
Following our previous work [13] here we investigate, in exemplary two-dimensional cases, different from those reported earlier, how the numerical treatment of the problem affects the magnitude of the obtained forces. Our aim is to stress, once more, that panel methods can accurately model the physics of flapping flight, once adequate assumptions are made. Additionally, these methods are much less time consuming than other numerical approaches, such as RANS and DNS methods, and can therefore provide meaningful insight into the underlying physics at a much lower computational cost.
3. Results and Discussion
The wake panel length can be set proportional to the time step length; see La Mantia and Dabnichki [13, 22] and option (iv) above, that is, inversely proportional to the Strouhal number St. The latter nondimensional number, which is often employed as a measure of the strength of oscillating motions (see, e.g., Eloy [5]), is mathematically expressed by the oscillation frequency (in Hz) times a suitably chosen motion amplitude, divided by the flow speed . This choice of wake panel length led to trailing-edge pressure differences that are, in most cases, lower than 10% of the maximum pressure coefficient, the pressure coefficient being defined aswhere is the fluid pressure, at the considered collocation point, obtained from the calculated fluid velocities, and indicates, as just mentioned, the steady fluid velocity far away from the oscillating foil. Additionally, the foil motion is imposed; that is, the NACA 0012 section, of chord and span , is moving forward at the velocity and oscillates harmonically with a heaving motion , perpendicular to , and a pitching motion . The oscillation amplitude used to define St is set to be two times the heaving amplitude , as it can be seen as a measure of the wake width in the proximity of the foil. The main kinematic parameters employed in the study are m, m, pitching axis position from the foil leading edge , m/s, , and maximum angle of attack deg. (e.g., for St = 0.3, the pitch amplitude deg. and Hz). The instantaneous angle of attack is not explicitly related to the initial conditions and the flow is considered attached on the oscillating foil. See our previous publications ([23] and references therein) for further details on the just sketched set of kinematic parameters.
In this section we initially show how a different choice of affects the computational results, option (iii) above. First of all, a wake panel length was defined to ensure a pressure difference close to zero (i.e., less than or equal to 10−3) at the trailing edge of the oscillating foil, at the first time step. A Newton-Raphson iterative procedure was developed to fulfil the constraint. Following (6) and (7) the pressure difference at the trailing edge is defined, at the first time step, aswhere the subscripts and indicate the lower and upper surfaces of the trailing edge, respectively. At the first iteration, a wake panel length is imposed and the computer program evaluates the forces acting on the oscillating foil. Then a second iteration starts and the new wake panels length is set equal to times , where . It is then possible to write the following equations for the third iteration, that is, to estimate :where and are the pressure differences at the trailing edge, at the first time step, and at the first and second iteration, respectively. If (estimated by using as the wake panel length) is less than or equal to 10−3, the iterative procedure stops. If is larger than 10−3, the iterative procedure continues similarly till the imposed condition is fulfilled. The corresponding convergence was reached within a few cycles, for values of up to ca. 0.2.
In Figure 2 the pressure coefficients in the trailing-edge region of the foil, at the first time step, are plotted for two different wake panels lengths (St = 0.3 and ). One hundred and eight time steps were used for each calculation, that is, three complete oscillations and thirty-six time steps for each foil cycle. It can be easily seen that is the wake panels length that ensures a null pressure difference at the trailing edge, at the first time step, for the considered set of parameters. It was calculated by using the Newton-Raphson iterative procedure explained above.
Figure 3 displays the thrust coefficient as a function of , for St = 0.3, being defined aswhere , which is the time-averaged force in the direction of the inertia frame of reference, which parallels the foil forward motion, is evaluated as discussed in La Mantia and Dabnichki [22], that is, by using the “Estimated Wake” assumption, option (i) above (such a force is negative if it is directed as the foil forward motion, i.e., parallel but opposite to the flow mean speed).
As the length of the wake panels increases, the thrust coefficient magnitude becomes larger; that is, if a wake panel is larger, its influence on the unsteady motion results in being larger. It can also be noted that the relation between the thrust coefficient and wake panel length is approximately linear, as it is expected from the equations governing the problem, and that a 50% increase of corresponds roughly to a 20% magnitude increase of the thrust coefficient. Additionally, as reported by La Mantia and Dabnichki [22], the thrust coefficient in the considered (experimentally relevant) conditions is indeed expected to be of a magnitude similar to that calculated here.
For the sake of completeness, in Figure 4 the influence on the thrust coefficient of the number of linear panels modelling the foil shape is displayed for St = 0.3 and . remains constant for larger than 200 (other relevant convergence studies can be found in [13, 22]).
The size of the linear wake panels has therefore a noticeable influence on the computational results. In some previous studies (see especially La Mantia and Dabnichki [13, 22]) a wake panel length proportional to the time step length was chosen, following Katz and Plotkin [24], in order to make its choice explicitly independent of the oscillation frequency and consequently allow positive comparisons with corresponding experimental results, at various St. As already mentioned, under this assumption, the trailing-edge pressure differences resulted, in most cases, in being lower than 10% of the maximum pressure coefficient (which is usually at the foil leading edge). In other words, an assumption about the wake panel size was needed and the choice was done on physical basis, that is, to accomplish the new formulation of the unsteady Kutta condition, although, as shown, for example, in Figure 3, a null trailing-edge pressure difference does not affect largely the value of the average thrust.
The latter outcome is also evident in Figure 5, where the thrust coefficient is plotted as a function of the Strouhal number (see again [22], for further details and relevant experimental values). If the pressure difference at the trailing edge is assumed to be null, the average values of thrust are comparable to those obtained by postulating the finite pressure difference enabling the existence of the thrust-producing jet. It can therefore be said that what is needed to obtain, from standard panel methods, thrust values that agree with the experimental ones is a suitable choice of the wake panel length, which can then be interpreted as one of the free parameters of the method. We propose here to make this choice on a physical basis, in order to account for the existence of the experimentally observed thrust-producing jet, which appears difficult to justify, in a potential flow, if the pressure at the trailing edge is zero.
Similar motion parameters were used to analyse the influence of the wake shape on the numerical results for a sample case (St = 0.3, , and ). One hundred and eight time steps were also employed for these calculations, that is, three complete oscillations and thirty-six time steps for each foil cycle.
If the wake shape is estimated by using the potential formulas, the evaluated thrust coefficient is −0.45 (“Estimated Wake,” option (i) above). If it is imposed that the wake follows the foil path (i.e., the wake panels remain static in the inertia coordinate system), (“Assumed Wake,” option (ii) above). In the latter situation the wake panels are generally closer to the foil than in the former case; that is, their influence on the foil is larger because they are closer.
Figure 6 shows the instantaneous force in the direction of the inertia frame of reference (i.e., the instantaneous thrust) as a function of the nondimensional time for the two different wake shapes ( indicates the beginning of the second cycle of the oscillating motion, while is its end). Note the slightly larger maximum thrust values achieved when the wake shape follows the foil path.
It should be emphasised once more that in the computer program the wake shape, that is, the positions of the wake panels in the global or local frame of reference, has to be set or estimated on the basis of the chosen assumptions, for example, the potential flow theory. However, as it can be seen in Figure 6 for a sample case, the wake shape exhibits relatively small influence on the obtained results; see also Mackowski and Williamson [20].
As mentioned above, we believe that the numerical implementation of the flapping-wing problem should be defined on a physical basis. This means that all the assumptions needed in the code development were related to experimental evidence or to theoretical (physically based) models. Those that could not be experimentally verified or theoretically justified were chosen for the sake of simplicity, for example, to have faster computations.
There is then a need for considering the differences between the actual flow conditions and potential assumptions. Such a consideration would provide an appropriate option to properly choose the mathematical model able to better represent the physical problem, once relevant experimental data on the pressure and velocity distribution over flapping wings are available.
It should be noted that the employed numerical procedure affects specifically the local pressure coefficients while the average thrust values seem to be less dependent on the chosen parameters; see also Ol et al. [19]. Physically this is quite an interesting phenomenon suggesting that thrust estimates in the present method are related to the global fluid conditions rather than to the local pressure distributions (on the foil boundary) and, in essence, it represents a further argument in favour of the potential methods ahead of highly computationally expensive methods using Navier-Stokes solvers with moving boundaries that are strongly dependent on local boundary conditions. The argument is further supported by the fact that potential methods are superior in dealing with one of the most important unsteady effects, that is, the added mass effect [26, 27], which has profound influence on the generated forces and the corresponding structural requirements [23]. In effect the method, with these additional features, is very robust in estimating the unsteady forces; hence adding reliability in the wake modelling is an essential step in order to ensure that it is also robust in assessing the pressure distribution, which should affect the stability of the autonomic motion of flapping wings. Virtually all up to date models (including the present one) do not allow autonomous motion but impose kinematic motion patterns. This partially explains why force generation remains independent of the accuracy of the pressure estimates along the foil. Current results point in the direction that an autonomous moving wing model would be dependent on a robust formulation of the unsteady Kutta condition.
4. Conclusions
The improved understanding of force generation and related energy requirements in flapping flight needs accurate solutions for the corresponding physical problems. Potential methods yield results very close to the experimental ones but often utilise simplified descriptions of the physical phenomena. An appropriate Kutta condition offers the opportunity to more accurately describe existing physical effects by incorporating relevant experimental results. In this work we discussed a novel approach to bridge this gap in the case of unsteady flows around an oscillating hydrofoil. The proposed Kutta condition takes into account an existing domain singularity at the foil trailing edge (generally neglected in the literature) to model closely the physics of the problem, that is, the trailing-edge pressure difference and the resulting thrust-producing jet. The computational results obtained by methods that do not account for such physical features agree well with experimental data but do not explain the thrust-generation mechanisms that are ultimately due to the trailing-edge pressure difference. Further research is needed in order to define theoretically how the pressure in the trailing-edge region should be defined and how this integral measure (in effect a Heaviside step function) could be represented as a function of the foil edge geometry, fluid properties, and flapping kinematics. Such a research should also allow solving the dynamics problem of autonomous flapping-wing movement, as a result of thrust generation, and should be applied to three-dimensional flows, too, on the basis, for example, of the related approach presented in La Mantia and Dabnichki [21].
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.