PHYSICS OF FLUIDS 17, 113601 !2005" Numerical experiments on flapping foils mimicking fish-like locomotion Paolo Blondeaux, Francesco Fornarelli, and Laura Guglielminia! Department of Environmental Engineering, University of Genova, Via Montallegro 1, 16145 Genova, Italy Michael S. Triantafyllou Department of Ocean Engineering, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, Massachusetts 02139 Roberto Verzicco DIMeG and CEMeC, Polytechnic of Bari, Via Re David 200, 70123 Bari, Italy !Received 9 February 2005; accepted 26 September 2005; published online 16 November 2005" The results of numerical experiments aimed at investigating the topology of the vortex structures shed by an oscillating foil of finite span are described. The motion of the foil and its geometry are chosen to mimic the tail of a fish using the carangiform swimming. The numerical results have been compared with the flow visualizations of Freymuth #J. Fluids Eng. 111, 217 !1989"$ and those of von Ellenrieder et al. #J. Fluid Mech. 490, 129 !2003"$. The results show that a vortex ring is shed by the oscillating foil every half a cycle. The dynamics of the vortex rings depends on the Strouhal number St. For relatively small values of St, the interaction between adjacent rings is weak and they are mainly convected downstream by the free stream. On the other hand, for relatively large values of St, a strong interaction among adjacent rings takes place and the present results suggest the existence of reconnection phenomena, which create pairs of longitudinal counter-rotating vortices. © 2005 American Institute of Physics. #DOI: 10.1063/1.2131923$ I. INTRODUCTION Fish usually propel themselves in water by means of a transverse wave that moves backwards along the body from the head to the tail. This propulsive mode can be broadly divided into two classes: the anguilliform mode and the carangiform mode. Fish using the anguilliform mode are flexible over the whole body and the propulsive wave has an amplitude that is almost constant along the body. In the carangiform mode, the amplitude of the wave becomes significant only in the posterior part of the fish, the remainder of the body being relatively rigid. Hence, thrust generation is confined almost exclusively to the caudal fin, which can be considered to oscillate in an otherwise undisturbed free stream. An analysis of the tail kinematics shows that, in order to investigate the propulsion mechanism, the tail can be approximated by a foil that moves with a constant forward speed, oscillates in the transverse direction, and simultaneously rotates, the two periodic motions being characterized by a 90° phase shift. To understand fish propulsion, a great deal of work has been devoted to the study of the dynamics of oscillating foils. The papers of Triantafyllou et al.1 and Mittal2 provide an excellent review of the experimental and numerical investigations, respectively. On the basis of two-dimensional experiments, Anderson et al.3 noted that the production of a significant leading-edge vortex, which travels downstream and interacts with the developing trailing edge vortex, is important for high efficiency. This conclusion has been supported by the recent two-dimensional numerical simulations a" Electronic mail: [email protected] 1070-6631/2005/17"11!/113601/12/$22.50 of Wang,4 Lewin and Haj-Hariri,5 Pedro et al.,6 and Guglielmini and Blondeaux.7 In particular, Guglielmini and Blondeaux7 obtained results for a wide range of parameters #in particular, the Strouhal number St ranged between 0.1 and 0.6, the ratio between the amplitude of the heaving oscillations and the foil chord ranged between 0.1 and 1.1, the phase between the pitching and heaving oscillations was varied between 65° and 105°, and the distance of the pitching axis from the leading edge of the foil made dimensionless with the foil chord fell in the range !−0.25, 1"$ and identified those giving rise to optimal propulsive performances, even though only moderate values of the Reynolds number were considered #in particular, the Reynolds number fell in the range !1000, 3000"$. In the two-dimensional approximation, only the crossstream component of the vorticity is taken into account. Therefore, the propulsive efficiency is overestimated, due to the underestimation of the wake energy loss. A threedimensional analysis is clearly required for more accurate and realistic estimates of the propulsive efficiency. Moreover, a three-dimensional approach able to handle different foil geometries should clarify the apparent evolutionary convergence on the fish tail toward a lunate form. A recent contribution to the understanding of the flow structure behind an oscillating foil of finite span is the work of von Ellenrieder et al.,8 where the overall structure of the flow is examined in detail. von Ellenrieder et al.8 used flow dye visualizations and tried to identify the vortex structures shed in the flow by looking at the time development of streak lines. Even though the approach is certainly valuable to determine the gross features of the flow, it has a few limitations. As discussed by von Ellenrieder et al.,8 dye filaments 17, 113601-1 © 2005 American Institute of Physics Downloaded 14 Dec 2005 to 18.172.7.243. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 113601-2 Phys. Fluids 17, 113601 "2005! Blondeaux et al. contain information that starts at the point of introduction and is integrated along the path of the dye lines. Hence, it is difficult to ascertain whether the shape of a streak line is the result of fluid distortion or flow memory, and a great deal of caution must be exercised in interpreting its behavior.9 For example, Hama10 analyzed the case of a shear flow perturbed by an unamplified traveling sinusoidal wave and showed that the streak line that emerges from the critical layer always appears to roll up regardless of the true nature of the superposed perturbation. Moreover, von Ellenrieder et al.8 pointed out that the Schmidt number of the dye was O!103" and, hence, it is likely that the dye followed the vorticity only during the initial stage of the flow development. Quantitative information on the wake structure behind a thrust-producing flapping foil of finite span began to be collected with use of the particle image velocimetry !PIV" technique,11,12 which allows the measurement of the velocity field to be made along with the evaluation of its phaseaveraged and random components. Moreover, the strength of the vortices shed by the foil can be estimated too. Presently, we investigate the flow that is generated by an oscillating finite-span foil by means of a numerical approach. One popular method that is available to compute the unsteady flow generated by arbitrary bodies moving at high Reynolds numbers is the potential flow panel method. A review of the various panel methods is given by Katz and Plotkin.13 Recently, Vest and Katz14 determined the shape of the trailing wake, by assuming that all vorticity is confined into the boundary layer on the surface of the body and in a thin vortex wake behind the body and by using a timestepping procedure. Smith et al.15 also developed an unsteady aerodynamic three-dimensional panel method that can realistically describe the kinematics of a flapping wing and can also account for wing flexibility. However, these panel methods assume an infinite value of the Reynolds number and, by ignoring the viscosity of the fluid, they fail to predict the details of the near-field flow, e.g., the large attached leading-edge vortex that has been recognized to play a key role in the propulsive efficiency3 and in the lift generation.16 In the present work, the solution of the full NavierStokes equations for the three-dimensional, unsteady flow generated by a flapping foil of finite span is determined. Numerical simulations have some advantages with respect to laboratory experiments. First, via numerical simulations one has access to velocity and pressure fields in the threedimensional space and in time. Hence, it is possible to make simultaneous velocity and pressure measurements in different points of the domain, while this would be prohibitively expensive and would be subject to unacceptable errors due to probe interference in laboratory experiments. Moreover, using a numerical approach, it is possible to “measure” quantities such as vorticity, strain rate, dissipation, etc., which are impossible to measure accurately in the laboratory. Even though computational approaches are limited in their reproduction of the actual flow by the relatively small values of the Reynolds number that can be efficiently computed, by the accuracy of the numerical methods !which rapidly deteriorates as the grid size is increased", and by the realism of the initial and boundary conditions, the use of numerical ap- ˜ , ˜Y , ˜Z" moving with the forward velocity −U FIG. 1. The reference frame !X 0 of the foil and the reference frame !X , Y , Z" moving with the foil. proaches is certainly valuable to investigate the details of the flow structure. In this paper, our attention is mainly focused on the vortices that are generated by flow separation both at the leading and trailing edges of the foil and on the vorticity shed at its tips. As in the paper by Dong et al.,17 a study of the topology and dynamics of the vortex structures is carried out to provide some physical insight into the mechanisms that control momentum transfer and the generation of lift and propulsive thrust. The procedure used in the rest of the paper is as follows. In the next section, the problem is formulated and the dimensionless parameters that control the results are listed. In Sec. III, the numerical code used to perform the numerical simulations is briefly described. In Sec. IV, the code is validated by comparing the results with the experimental visualizations of Freymuth.18 Then, the phenomenon is investigated for values of the parameters close to those used in the laboratory experiments of von Ellenrieder et al..8 The numerical results are compared against dye visualizations and discussed to point out the main features of the vortical flow. Some conclusions are drawn in Sec. V, where both ongoing and future research are also described. II. THE PROBLEM A foil, characterized by a finite aspect ratio, moves with a constant forward speed −U0 and oscillates with a combination of harmonic vertical and angular oscillations !see Fig. 1". The foil has a symmetrical cross section and a rectangular planform. Its chord and thickness are denoted with C and S, respectively, while its span is B. If a right-handed reference ˜ , ˜Y , ˜Z" moving with velocity −U is introduced, such frame !X 0 that the ˜X coordinate is aligned with U0 and the ˜Z coordinate is in the spanwise direction, the motion of the foil is described as function of the time T by providing the ˜Y coordinate A!T" of the point O around which the foil pitches and the angle !!T" that the foil chord forms with the ˜X axis, ˜Y = A!T" = Aˆ cos!2"FT", O !1" Downloaded 14 Dec 2005 to 18.172.7.243. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 113601-3 Phys. Fluids 17, 113601 "2005! Numerical experiments on flapping foils !!T" = ¯! + !ˆ cos!2"FT + #". !2" Here Aˆ and !ˆ denote the amplitudes of the heaving and pitching oscillations, respectively, F is the frequency of the oscillations, and # is the phase lead of the angular motion with respect to the lateral motion. Moreover, ¯! denotes the mean angle of attack. The flow generated by the fluid flowing around the oscillating foil can be determined by solving momentum !Navier-Stokes" and continuity equations. By choosing the foil chord C as length scale, the period 1 / F of the foil oscillations as time scale, and U0 as velocity scale, it can be easily verified that the dimensionless problem is characterized by eight dimensionless parameters beside the geometry !thickness and shape" of the foil section: !i" the ratio b = B / C between the span B of the foil and its chord C, !ii" the ratio a = Aˆ / C between the amplitude of the heaving oscillations and C, !iii" the amplitude !ˆ of the pitching oscillation, !iv" the phase # between the pitching and heaving oscillations, !v" the mean angle of attack ¯!, !vi" the ratio d = D / C, where D is the distance of the pitching axis from the leading edge of the foil, !vii" the Reynolds number of the flow Re= U0C / $, defined with the forward velocity of the foil, the chord of the foil, and the kinematic viscosity $ of the fluid, and !viii" the Strouhal number St= 2FAˆ / U0, defined using the frequency F of the foil oscillations. The numerical approach solves the problem in the Cartesian coordinate system !X , Y , Z", which moves with the foil !see Fig. 1", using primitive variables. The details of the problem are given in the Appendix. III. THE NUMERICAL APPROACH In the present work, the problem has been solved numerically looking for the values of the velocity and pressure, at a finite number of grid points in the fluid domain, which approximate the continuous distribution of the actual velocity and pressure fields. The numerical method uses a finite-difference approach. Standard centered second-order finite-difference approximations of the spatial derivatives are used. The time advancement of the Navier-Stokes equations, written in the noninertial reference system !X , Y , Z", employs the fractional-step method extensively described in Refs. 19 and 20. The nonsolenoidal intermediate velocity field is evaluated by means of a third-order Runge-Kutta scheme to discretize convective terms, together with a Crank-Nicholson scheme for the diffusive terms. The implicit treatment of the viscous terms would require the inversion of large sparse matrices, which are reduced to three tridiagonal matrices by a factorization procedure with an error of order !%T"3.21 Then, by forcing a continuity equation, a Poisson equation for the pressure field is obtained, which is solved by standard methods. The complex geometry that characterizes the foil section is handled by means of the immersed boundary technique,22,23 which is based on the introduction of fictitious body forces in the fluid domain along the contour of the foil. The strength of these fictitious forces is chosen in such a way that the no-slip con- dition at the foil surface is forced without the need to introduce any special treatment of the problem at the surface of the foil. A nonreflective boundary condition at the outflow boundary allows the use of a short computational domain in the X direction. Different computational boxes are used, the size of which is given in the following. The mesh is uniform in the spanwise direction, while in the streamwise and vertical ones a nonuniform mesh has been used to cluster the grid points in the vicinity of the foil where gradients are expected to be stronger. Visualizations of the trajectories of passive tracers have been used to obtain preliminary information on the velocity field. However, the results are not discussed hereafter because the analysis of the streak lines cannot be used to surely identify vortex structures. After having gained a qualitative picture of the flow structure, a quantitative investigation is made of the coherent vortices, having access to velocity and vorticity fields in the three-dimensional space and time. In addition, since the analysis of the isovorticity surfaces may be considered not fully adequate to detect vortices, we also use the definitions of vortex structures proposed by Jeong and Hussain24 and Zhou et al.25 If the unsteady irrotational straining and the viscous effects are neglected, the pressure Hessian is proportional to D2 + R2, where D and R are the symmetric and antisymmetric part of the velocity gradient tensor, respectively. The above fact can be used to determine the existence of a local pressure minimum due to vortical motion. The vortex core is so defined as a connected region with two negative eigenvalues of D2 + R2.24 Since D2 + R2 is a symmetric tensor, all eigenvalues are real. If the eigenvalues are ordered in accordance with their decreasing values !&1 ' &2 ' &3", the new definition is equivalent to the requirement &2 ( 0 within the vortex core. As discussed by Jeong and Hussain,24 this definition captures the pressure minimum in a plane perpendicular to the vortex axis at high Reynolds numbers and also accurately defines vortex cores at low Reynolds numbers, unlike a pressure minimum criterion, which sometimes is used as an intuitive indicator of a vortex. Moreover, as suggested by Zhou et al.,25 we visualize also the imaginary part &ci of the complex eigenvalues of the velocity gradient tensor. Indeed, if &r is the real eigenvalue of the velocity gradient tensor, with a corresponding eigenvector Vr, and &cr ± ı&ci are the conjugate pair of the complex eigenvalues, with complex eigenvector Vcr ± ıVci, the local flow is either stretched or compressed along the axis Vr, while on the plane spanned by the vectors Vcr and Vci the flow is swirling. Moreover, the strength of the local swirling motion is quantified by &ci and hence the imaginary part of the complex eigenvalue pair is referred to as the local “swirling strength” of a vortex. Numerical experiments suggest the values of the modulus %!%, the values of the second eigenvalue &2 of the matrix D2 + R2 and of &ci that should be used to visualize the vortex structures generated by the foil oscillations. While %!% and &ci are of order U0 / C, the second eigenvalue of D2 + R2 is of order !U0 / C"2. Of course, the larger the values are, the smaller are the regions visualized by the approaches. Moreover, it is worth pointing out that the approaches do not visualize the vortex tubes. Hence, even Downloaded 14 Dec 2005 to 18.172.7.243. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 113601-4 Blondeaux et al. Phys. Fluids 17, 113601 "2005! FIG. 3. Square wing in rapid pitch from 0° to 30° in 0.5 s. On the left, the vortical structures observed by Freymuth !Ref. 18" are sketched on the basis of his Fig. 3. On the right, the isosurface of %!% = 10U0 / C obtained by means of the present numerical simulation is plotted. %T = 1 / 8 s. results it provides with the results obtained with the code N)*+!r,26,27 developed at Brown University and based on a finite-element approach. The interested reader is referred to Blondeaux et al.,28 where a detailed quantitative comparison is described. A. Simulation of Freymuth’s experiments FIG. 2. Rectangular wing in pitch up-down maneuver from 0° to 30° to 0° within 0.5 s. On the left, the vortical structures observed by Freymuth !Ref. 18" are sketched on the basis of his Fig. 4. On the right, the isosurface %!% = 12U0 / C obtained by means of the present numerical simulation is plotted. %T = 1 / 8 s. though sometimes the vortex structures described in the following seem to begin and end in the flow, the Helmholtz vortex theorems do not break down. IV. THE RESULTS The reliability and the accuracy of the numerical code presently used is supported by a comparison between the To further support the numerical results, some of the flow visualizations of Freymuth18 are reproduced. In a first run, b is set equal to 2 and the foil profile is chosen similar to a NACA0012 profile. The computational box in the !X , Y , Z" coordinate system is 5.5 foil chords long, 3 chords wide, and 4 chords high. Moreover, 240, 151, 220 grid points have been used in the streamwise, spanwise, and cross-stream directions, respectively. The foil is described using 60, 50 , 13 points, respectively. In order to limit the computational costs, the Reynolds number is fixed equal to 1000, rather than 5100 as in Ref. 18. For both Re= 1000 and Re= 5100, the size of the vortex structures shed by the foil is much larger than the thickness of the boundary layer that develops along the foil surface and the thickness of the free shear layers shed by the foil. Hence, no relevant difference is expected between the dynamics of the vortex structures presently computed and that observed by Freymuth.18 Only, the viscous decay of the vortices, taking place far from the foil, is envisaged to be different. Since Freymuth18 visualized the flow generated by an aperiodic motion of the foil, relationships !1" and !2" are modified in accordance with the experi- Downloaded 14 Dec 2005 to 18.172.7.243. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 113601-5 Numerical experiments on flapping foils Phys. Fluids 17, 113601 "2005! FIG. 4. Isocontours of the streamwise vorticity -x at . = −" / 2, when the foil is at its middle heave position, for a = 0.5, !ˆ = 5 ° , St= 0.35, # = 90° , Re ˜ = 0. Second row: portion of the plane ˜Z = 1.55C, located at distance equal to C / 20 from the foil tip. The results on the = 164. First row: portion of the plane X left have been obtained with 311, 188, 325 grid points, the results on the right have been obtained with 311, 188, 257 grid points. In the figure, the spatial coordinates are made dimensionless using the foil chord C. mental data. In particular, the transverse position of the foil is kept fixed and equal to zero. Starting from ! = 0 at T = 0, the foil is rotated up to a value of ! = 30° at T = 0.25 s and then rotated back to 0° at T = 0.5 s, where it is again kept fixed. Figure 2 shows the numerical predictions of the vorticity field at different times, along with the streak lines visualized by Freymuth.18 A vortex ring is formed by the two counter-rotating trailing edge vortices that develop during the pitch-up and pitch-down phases and by the vorticity shed by the tips of the foil. The remaining parts of the edge vortex loops, which were previously connected to the leading-edge corners, annihilate. Also, the leading-edge dynamic stall vortex filaments seem to form a ring. However, in the experiments the dye pattern becomes so entangled that the leadingedge ring has the appearance of a puff. In the numerical simulation, the leading-edge vortices give rise to a vortex structure that moves downstream along the lee side of the foil and then is shed in the free stream. The qualitative agreement between the numerical simulations and the flow visualizations of Freymuth18 is good and supports the reliability of the numerical results, even though the vortex structures in the laboratory experiments are slightly farther away from the foil than the numerically simulated vortices. Some of the discrepancy might be due to the deficiencies of the numerical approach and in particular to the difference between the experimental and numerical Reynolds numbers. However, part of the quantitative discrepancy can be certainly ascribed to the possible differences between the motion of the foil in the laboratory experiment and that in the numerical simulation. Indeed, the details of the foil kinematics are not given in Ref. 18 so that a sinusoidal time development has been presently chosen for the foil angular velocity. Figure 3 shows a comparison between the computed and the visualized vorticity fields for a foil that, starting from a vanishing value of the angle of attack, is rotated up to 30° in 0.5 s and then kept fixed. In this case, the dimensionless span of the foil is equal to 1 !b = B / C = 1". The large vortex loop that develops on the right of the foil is formed by the vorticity generated along the foil and separated from the trailing edge and from the wing tips during the pitch-up. The loop connects to the front corners of the wing. Simultaneously, because of the high angle of attack, vorticity also separates from the leading edge and a dynamic stall vortex appears. This connects to the front corners of the wing, thus completing a closed vortex system in accordance with Helmholtz’s law. Subsequently, the anchoring of the combined leadingand trailing-edge vortex system to the front corners gets Downloaded 14 Dec 2005 to 18.172.7.243. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 113601-6 Blondeaux et al. Phys. Fluids 17, 113601 "2005! FIG. 5. Isosurface of %!% = 1.5U0 / C for a = 0.5, !ˆ = 5 ° , ¯! = 0 , # = 90° , Re= 164, St= 0.35. Side and plan views at !a" . = −" / 2, !b" . = −" / 4, !c" . = 0, !d" . = " / 4, after about four cycles. In the figure the spatial coordinates are made dimensionless using the foil chord C. somewhat lost, presumably because of viscous and turbulent annihilation of the nearly counter-rotating leading- and trailing-edge tip vortices. It should be mentioned that in the experimental procedure of Freymuth,18 some smoke was inadvertently introduced in regions where no vorticity exists and contributed to the “optical” noise of the photograph. Also in this case the qualitative agreement between the flow visualizations and the numerical results is good, even though quantitative differences are present. As previously pointed out, these differences might be due to deficiencies of the numerical simulation and to the possible differences between the kinematics of the foil in the laboratory and numerical experiments. Indeed, after 1 / 8 s, the trailing-edge vortex is well developed and it has moved far from the foil in the laboratory experiment, while it is weaker and still close to the foil in the numerical simulation, thus indicating that the actual rotation of the foil is different. B. Simulation of von Ellenrieder et al.’s experiments As pointed out in the formulation of the problem, the results of the numerical code depend on eight dimensionless parameters beside the geometry !thickness and shape" of the foil. Hence, the high computational costs do not allow an exhaustive investigation of the phenomenon in the parameter space. In this work, the phenomenon has been simulated for values of the parameters equal to those considered in the experiments carried out by von Ellenrieder et al.8 A detailed investigation of the velocity and vorticity fields has been carried out to identify possible differences between the streak-line behavior described by von Ellenrieder et al.8 and the actual vorticity field and to support or modify the isometric sketch of the flow structure proposed by the above authors. A foil characterized by a value of b equal to 3 has been Downloaded 14 Dec 2005 to 18.172.7.243. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 113601-7 Phys. Fluids 17, 113601 "2005! Numerical experiments on flapping foils FIG. 6. Isosurface of %!% = 2.0U0 / C for the region specified in Fig. 5!c". Side and plan views. In the figure, the spatial coordinates are made dimensionless using the foil chord C. considered and the section of the foil has been taken similar to that of a NACA0030 profile. The dimensionless distance d of the pitching axis from the leading edge of the foil is equal to 1 / 4 and the ratio a between the amplitude of the heaving oscillations and the foil chord is fixed equal to 0.5. Moreover, the Reynolds number is about 164, as the experimental value. The pitching amplitude !ˆ and the phase # between the two oscillatory motions have been fixed equal to 5° and 90°, respectively, while a vanishing value of ¯! is considered. Simulations have been performed for two different values of the Strouhal number, namely St= 0.35 and St= 0.175. The computational box in the !X , Y , Z" coordinate system is 11.75 foil chords long, 8 chords high, and 5.6 chords wide. Moreover, 311, 188, 325 grid points have been used in the streamwise, cross-stream, and spanwise directions, respectively. The foil is described using 118, 36, 177 points, respectively. It takes about 24 h to simulate one oscillation cycle using a PC AMD Athlon XP 3000. Some preliminary tests have been made to ascertain that the size of the computational grid is small enough to assure a fair accuracy of the results. For example, Fig. 4 shows, at a certain phase of the cycle, a comparison between the distribution of the streamwise vorticity component obtained with the number of grid points previously given and that obtained decreasing the number of grid points in the spanwise direction !257 grid points have been used". The plots in the first row show a portion of the plane ˜X = 0, while the plots in the second row show a portion of the plane ˜Z = 1.55C, which is located at a distance equal to C / 20 from the foil tip. The agreement between the two vorticity fields supports the accuracy of the obtained results. The streamwise vorticity component and the region around the foil tips are shown because strong velocity gradients are expected along the spanwise direction at the tips of the foil and the code employs a uniform grid spacing in the Z direction. A run has been also made with a shorter computational box !7.5, 8 , 5.6 foil chords" and using the same number of points, so that it has been verified that the flow close to the foil and in the far wake is not affected by the length of the computational box. Moreover, at the beginning of the calculations, both the fluid and the foil, which are initially at rest, are impulsively set into motion. All the simulations have been carried out for a number of cycles large enough to assure that the flow close to the foil no longer feels the influence of the initial conditions. The time development of the flow over half a heave/ pitch cycle is shown for St= 0.35 in Fig. 5, where an isosurface of the modulus of vorticity !%!% = 1.5U0 / C" is plotted at different phases of the cycle, after four cycles of oscillation. As in von Ellenrieder et al.,8 the figure displays both the plan and side views in order to give a sense of the threedimensional structure of the flow. Motion parameters and cycle phases are the same of von Ellenrieder et al.’s8 Figs. 3!c"–3!f". By looking at Fig. 6, which is an enlargement of Fig. 5!c", it is possible to recognize the following: !i" !ii" !iii" !iv" !v" The leading-edge vortex !side view" moving along the lower surface of the foil. The trailing-edge vortex rolling close to the foil in the span-wise direction !plan view". The leading-edge vortex, shed during the previous half-cycle, which merged with trailing-edge vorticity and constitutes the head of the developing vortex ring !side view". The tip vorticity !plan view". The vortex rings shed at the previous half-cycles. These findings are similar to those of Ref. 8 even though the ring vortex structures identified in the laboratory experiments are slightly longer and larger. Since the locations of dye emission, which strongly influence streak-line patterns, are not specified in Ref. 8 and taking into account that it is quite difficult to reproduce the details of dye injection !geometry of the injectors, initial velocity of dye, etc.", no attempt has been made to make a quantitative comparison between the computed streak lines and the visualizations of von Ellenrieder et al.8 Since the key aspect of the work is the identification of vortical structures, we have tried some of the various techniques proposed in the literature. A first, simple, and quite direct method to identify vortex structures is represented by the isosurfaces of the strength of vorticity !see Fig. 5", even if the largest values of the vorticity occur at the foil surface or close to it, where the flow may exhibit shear but not necessary a swirling motion.24 Moreover, in the experiments, dissipative effects are rather strong and the vorticity decays along the wake quite rapidly, so that small values of %!% should be used if the isovorticity surface should identify the vortex rings in their completeness. Hence, we cannot avoid losing the details of the vortex structures in those regions where vorticity is large. Another simple and quite direct method to identify vortex structures is represented by the isosurfaces of the pressure distribution !see Blondeaux et al.28". Indeed, local pressure minima are related to the presence of vortex structures. However, there is an inherent scale difference between a vor- Downloaded 14 Dec 2005 to 18.172.7.243. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 113601-8 Blondeaux et al. Phys. Fluids 17, 113601 "2005! FIG. 7. Isosurface of &ci2 #&ci2 = 0.4!U0 / C"2$ for a = 0.5, !ˆ = 5 ° , ¯! = 0 , # = 90° , Re= 164, St= 0.35, in a perspective view that looks at the foil from the bottom. !a" . = −" / 2, the foil is drawn in black gray; !b" . = −" / 4, the isosurface is transparent and superimposed to &ci2 contour levels in the plane ˜Z = 0; !c" . = 0, the isosurface is transparent and superimposed to &ci2 contour levels in the planes ˜Z = 0 and ˜Y = 0.2C; !d" . = " / 4, the isosurface is transparent and superimposed to /ci2 contour levels in the planes ˜Z = 0 and ˜Y = 0. In !b", !c", and !d", the foil is not displayed. tex core and the associated low-pressure region, and this makes the demarcation of vortices using an isosurface of pressure problematic.24 The so-called “&ci” criterion seems to be more appropriate to the description of the vortex dynamics both in the region close to the foil and along the wake. In Fig. 7, a surface characterized by a constant value 2 is shown at four phases of the cycle. Now the structure of &ci of the vorticity field close to the foil becomes more understandable, as does the structure of the vortex rings that leave the foil every half a cycle. Looking at dye filaments, von Ellenrieder et al.8 distinguished between the vortex structures shed by the leading edge of the foil, which they labeled Li !the index i indicating the half-cycle", and the vortices shed by the trailing edge, which they labeled Ti. However, if the vorticity time development is analyzed in the plane of symmetry of the foil !see Fig. 8", the present results show that the vorticity shed by the leading edge merges with that generated all along the foil surface and that shed by the trailing edge, to produce a unique vortex structure every half-cycle. An analysis of the sequential perspective views of the regions characterized by a constant value of &2ci shows that these unique vortex structures are shed in alternating order near to the extrema of the heaving motion. These vortices tend to be concave and form vortex rings being connected to the tip vorticity and the spanwise vorticity of opposite sign shed during the following half-cycle. When the vortex rings are convected downstream, the vortex structures are stretched and tend to assume a hook shape in the vertical plane. As the cycle continues, the vortex filaments move toward the midspan plane and simultaneously wrap around and combine with the head of the following vortex. Comparing the results displayed in Figs. 5 and 7 with those of Fig. 3 of von Ellenrieder et al.,8 it can be ˜ , ˜Z" seen that the vortex wake of the foil is smaller in the !X ˜ , ˜Y " plane. plane and larger in the !X Downloaded 14 Dec 2005 to 18.172.7.243. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 113601-9 Numerical experiments on flapping foils Phys. Fluids 17, 113601 "2005! ˜ = 0" for a = 0.5, !ˆ = 5 ° , ¯! = 0 , # = 90° , Re= 164, St= 0.35, at !a" . = −" / 2, !b" FIG. 8. Contour levels of -z in the plane of symmetry of the foil !Z . = −" / 4, !c" . = 0, !d" . = " / 4 , %-z = 1.85U0 / C. In the figure, the spatial coordinates are made dimensionless using the foil chord C. If we look at the contour plot of the spanwise vorticity -z in the plane of symmetry of the foil ˜Z = 0 !see Fig. 8", we find that the vorticity generation along the surface of the foil is not much different from that characterizing the twodimensional flow around a foil of infinite length !see Fig. 9". Plots refer to four phases during half-cycle. During the upward motion of the foil, flow separation from the leading edge is observed, associated with a weak back flow along the lower surface of the foil. The free shear layer rolls up and gives rise to a negative vortex structure that is convected downstream. When the foil reverses its motion, new negative vorticity is generated around the trailing edge of the foil because of its downward motion. The evolution of these two vortices can be followed looking at the plot for . = −" / 2 #Fig. 8!a"$ and considering the mirror image of the positive vorticity shed during the previous half-cycle at the upper FIG. 9. Contour levels of -z for the 2D simulation of a high aspect ratio foil similar to a NACA0030, with a = 0.5, !ˆ = 5°, ¯! = 0, # = 90°, Re= 164, St= 0.35, at !a" . = −" / 2, !b" . = −" / 4, !c" . = 0, !d" . = " / 4 , %-z & 0.46U0 / C. In the figure, the spatial coordinates are made dimensionless using the foil chord C. Downloaded 14 Dec 2005 to 18.172.7.243. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 113601-10 Phys. Fluids 17, 113601 "2005! Blondeaux et al. FIG. 11. Isosurface of &2 = −0.7!U0 / C"2 for !ˆ = 5 ° , a = 0.5, # = 90° , St = 0.175 at . = 3" / 4. FIG. 10. Isosurface of &2#&2 = −0.3!U0 / C"2$ for !ˆ = 5 ° , a = 0.5, # = 90° , St= 0.35, at . = 3" / 4 after about four cycles. The isosurface is translucent and superimposed to contour levels of the streamwise vorticity -x at planes ˜X = −2 . 0C , − 3.3C , − 4.65C !negative values are dashed". Arrows indicate the approximate direction and versus of vorticity. surface of the foil. Far from the foil, some differences appear between the 2D and 3D cases. Indeed, in the 2D case, the vorticity strength decreases less rapidly as the vortex structures are convected downstream. Moreover, looking at the isosurface of &ci at the same phases of the cycle in the threedimensional case !see Fig. 7, where the foil is seen in a bottom view", the same vortices with positive -z generated at the upper surface of the foil, namely the leading-edge vortex, which is going to constitute the head of the new vortex ring, and the trailing-edge vortex, can be recognized. However, these vortices are connected by the vorticity generated at the tips of the foil: the two arms of the tip vorticity elongate, lose connection with the tips, and a ring is formed and shed when the foil is leaving its extreme transversal position. The leading-edge vortex with negative -z shed at the lower surface of the foil contributes to the formation of the ring connecting the tip vorticity at the front. Looking at the spanwise velocity component along a different vertical plane crossing the foil at different values of the coordinate ˜X, we can observe a spanwise symmetrical transport of momentum toward the foil extrema to feed the strong tip vortices. These results agree with Maxworthy’s29 finding, which shows a strong spanwise transport of vorticity toward the wing tips of a three-dimensional foil. The arrows in Fig. 10, lying along the vortex structures identified by a constant negative value of &2, indicate the direction of circulation and have been obtained plotting the different vorticity components provided by the calculation. Indeed, even though the &ci and &2 criteria are reliable approaches to visualize the regions where vortex structures are present, they are unable to provide the sense of rotation of the vortex structures. In accordance with the rule that no vortex can simply begin or end in the flow, the vortices form a continuous system of interconnected loops that extends downstream of the foil. To give another example, further results are shown in Fig. 11, where an isosurface of &2 is represented for a = 0.5, !ˆ = 5 ° , ¯! = 0 , # = 90° , Re= 164, St= 0.175, i.e., for a smaller value of the Strouhal number. In this case, the shed vorticity is weaker, the vortex filaments turn out to be less curved, and the characteristic length of the vortex structures is approximately twice that for St= 0.35. Moreover, the wake tendency to become narrower in the ˜Z direction and to open and diverge in the ˜Y direction is less evident. V. CONCLUSIONS In Lighthill,30 it is first explained how vortex lines shed from the central part of the tail/foil must bend downstream, then turn inward and close up when they come close to the vorticity of opposite sign shed by the tail/foil half a period later. In accordance with Lighthill,30 if the foil is supposed to flap in the vertical plane moving alternatively up and down, the wake turns out to consist of a series of horizontal vortex rings. The experimental works of Drunker and Lauder31,32 and Lauder,33 as well as the recent flow visualizations of von Ellenrieder et al.,8 have confirmed the formation of ringlike vortical structures by oscillating tails/foils. The present numerical simulations also show the generation of vortex rings by an oscillating foil. Indeed, the results show that, for the chosen parameters, the leading vortices, shed by the foil for high angle of attack, merge with the trailing-edge vorticity, giving rise to large vortex structures. These originate ring-like structures being connected to alternating sign tip vortices and to the spanwise vorticity of opposite sign shed during the following half-cycle. Numerical simulations show that sometimes the rings have a complex geometry and are characterized by a complex dynamics. In particular, on the basis of the present results, it is possible to guess the existence of reconnection phenomena at relatively high Strouhal numbers. As described in Kida and Takaoka,34 first the vorticity tubes tend to approach each other advected by the mutual- and selfinduction velocity. As the two vortex tubes get closer, the shape of the vortex core is deformed.35 This can be seen in Fig. 12, which shows a detail of the vorticity field plotted in Fig. 7. Then, the generation of vorticity in the direction per- Downloaded 14 Dec 2005 to 18.172.7.243. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 113601-11 Phys. Fluids 17, 113601 "2005! Numerical experiments on flapping foils FIG. 12. Isosurface of &2 #&2 = −0.3!U0 / C"2$ for !ˆ = 5 ° , a = 0.5, # = 90° , St= 0.35 at . = 3" / 4 after about four cycles. The region that is emphasized in the top figure is enlarged in the two bottom plots. In the bottom left figure, the isosurface of &2 is colored with the contour levels of -x. In the right bottom figure, the isosurface of &2 is colored with the contour levels of -z. ' ( pendicular to the original vortex tubes, which can be observed in Fig. 12, gives rise to vorticity reconnection. However, the viscous cancellation of opposite signed vorticity is not complete in Fig. 12 and the original rings are not definitively broken. Recent detailed data on three-dimensional foils by Martin et al.36 show little degradation of propulsive performances for moderate aspect ratio foils, compared with the results of two-dimensional foils. The evaluation of the forces acting on the foil should clarify this point and the simulation of different foil geometries should clarify the apparent evolutionary convergence of the fish tail toward a lunate form and/or suggest another form to obtain optimal efficiency. # V d 2! d! − 2X+ # T dt dt +V 2 Y− #V d 2A d! U+U 2 cos ! − 2 #X dt dt ) * #V #V 1 #P # 2V # 2V # 2V +W =− +$ + + , #Y #Z 0 #Y # X2 # Y 2 # Z2 !A3" #W #W #W #W +V +W +U #X #Y #Z #T =− ) * 1 #P # 2W # 2W # 2W +$ + + . 0 #Z # X2 # Y 2 # Z2 !A4" ACKNOWLEDGMENTS Thanks are due to the Office of Naval Research for funding the research under Contract No. 000140310193 monitored by Dr. Thomas F. Swean. APPENDIX Continuity and momentum equations in the moving reference frame !X , Y , Z": #U #V #W + + = 0, #X #Y #Z ' ( # U d 2! d! + Y− # T dt2 dt +V !A1" 2 X− #U d 2A d! sin ! + 2 V + U #X dt2 dt ) * #U #U 1 #P # 2U # 2U # 2U +W =− +$ + + , #Y #Z 0 #X # X2 # Y 2 # Z2 !A2" 1 M. S. Triantafyllou, A. H. Techet, and F. S. Hover, “Review of experimental work in biomimetic foils,” IEEE J. Ocean. Eng. 29, 585 !2004". 2 R. Mittal, “Computational modeling in biohydrodynamics: Trends, challenges, and recent advances,” IEEE J. Ocean. Eng. 29, 595 !2004". 3 J. M. Anderson, K. Streitlien, D. S. Barret, and M. S. Triantafyllou, “Oscillating foils of high propulsive efficiency,” J. Fluid Mech. 360, 41 !1998". 4 Z. J. Wang, “Vortex shedding and frequency selection in flapping flight,” J. Fluid Mech. 410, 323 !2000". 5 G. C. Lewin and H. Haj-Hariri, “Modelling thrust generation of a twodimensional heaving airfoil in a viscous flow,” J. Fluid Mech. 492, 339 !2003". 6 G. Pedro, A. Suleman, and N. Djilali, “A numerical study of the propulsive efficiency of a flapping hydrofoil,” Int. J. Numer. Methods Fluids 42, 493 !2003". 7 L. Guglielmini and P. Blondeaux, “Propulsive efficiency of oscillating foils,” Eur. J. Mech. B/Fluids 23, 255 !2004". 8 K. D. von Ellenrieder, K. Parker, and J. Soria, “Flow structures behind a heaving and pitching finite-span wing,” J. Fluid Mech. 490, 129 !2003". 9 D. I. Pullin and A. E. Perry, “Some flow visualization experiments of the starting vortex,” J. Fluid Mech. 97, 239 !1980". 10 F. R. Hama, “Streaklines in a perturbed shear flow,” Phys. Fluids 5, 644 !1962". Downloaded 14 Dec 2005 to 18.172.7.243. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp 113601-12 11 Blondeaux et al. K. Parker, K. Ellenrieder, and J. Soria, “Using stereo multigrid DPIV !SMDPIV" measurements to investigate the vortical skeleton behind a finite-span flapping wing,” Exp. Fluids 39, 281 !2005". 12 K. Parker, K. von Ellenrieder, and J. Soria, “The flow of a threedimensional thrust producing wing,” in Advances in Turbulence, 10th European Turbulence Conference, Fluid Mechanics and Its Applications, edited by P. A. Krogstad and H. I. Andersson !CIMNE, Trondheim, Norway, 2004", Vol. 49, pp. 11–14. 13 J. Katz and A. Plotkin, Low-Speed Aerodynamics from Wing Theory to Panel Methods !McGraw-Hill, New York, 1991". 14 M. S. Vest and J. Katz, “Unsteady aerodynamic model of flapping wing,” AIAA J. 34, 1435 !1996". 15 M. Smith, P. Wilkin, and M. Williams, “The advantages of an unsteady panel method in modelling the aerodynamic forces on rigid flapping wings,” J. Exp. Biol. 199, 1073 !1996". 16 C. P. Ellington, C. van der Berg, A. P. Willmott, and A. L. R. Thomas, “Leading-edge vortices in insect flight,” Nature 384, 626 !1996". 17 H. Dong, M. Mittal, M. Bozkurttas, and F. Najjar, “Wake structure and performance of finite aspect ratio flapping foils,” AIAA Technical Papers 2005-81, 43rd AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV !2005". 18 P. Freymuth, “Visualizing the connectivity of vortex systems for pitching wings,” J. Fluids Eng. 111, 217 !1989". 19 J. Kim and P. Moin, “Application of fractional step method to incompressible Navier-Stokes equations,” J. Comput. Phys. 59, 308 !1985". 20 M. M. Rai and P. Moin, “Direct simulations of turbulent flow using finitedifference schemes,” J. Comput. Phys. 96, 15 !1991". 21 R. M. Beam and R. F. Warming, “An implicit finite-difference algorithm for hyperbolic systems in conservation law form,” J. Comput. Phys. 22, 87 !1976". 22 J. Mohd-Yusof, “Combined immersed boundaries/B-splines methods for simulations of flows in complex geometries,” Technical Report, CTR Annual Research Briefs, NASA Ames, Stanford University, 1997. 23 E. A. Fadlun, R. Verzicco, P. Orlandi, and J. Mohd-Yusof, “Combined Phys. Fluids 17, 113601 "2005! immersed-boundary finite-difference methods for three-dimensional complex flow simulations,” J. Comput. Phys. 161, 35 !2000". 24 J. Jeong and F. Hussain, “On the identification of a vortex,” J. Fluid Mech. 285, 69 !1995". 25 J. Zhou, R. J. Adrian, S. Balachandar, and T. M. Kendall, “Mechanisms for generating coherent packets of hairpin vortices in channel flow,” J. Fluid Mech. 387, 353 !1999". 26 S. J. Sherwin, Ph.D. thesis, Princeton University !1995". 27 T. C. Warburton, Ph.D. thesis, Brown University, Division of Applied Mathematics !1999". 28 P. Blondeaux, F. Fornarelli, L. Guglielmini, M. S. Triantafyllou, and R. Verzicco, “Vortex structures generated by a finite-span oscillating foil,” AIAA Technical Papers 2005-84, 43rd AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV !2005". 29 T. Maxworthy, “Experiments on the Weis-Fogh mechanism of lift generation by insects in hovering flight. Part 1. Dynamics of the ‘fling’,” J. Fluid Mech. 93, 47 !1979". 30 M. J. Lighthill, “Hydromechanics of aquatic animal propulsion,” Annu. Rev. Fluid Mech. 1, 413 !1969". 31 E. G. Drucker and G. V. Lauder, “Locomotor forces on a swimming fish: Three-dimensional vortex wake dynamics quantified using digital particle image velocimetry,” J. Exp. Biol. 202, 2393 !1999". 32 E. G. Drucker and G. V. Lauder, “A hydrodynamic analysis of fish swimming speed: Wake structure and locomotor force in slow and fast labriform swimmers,” J. Exp. Biol. 203, 2379 !2000". 33 G. V. Lauder, “Function of the caudal fin during locomotion in fishes: Kinematics, flow visualization, and evolutionary patterns,” Am. Zool. 40, 101 !2000". 34 S. Kida and M. Takaoka, “Vortex reconnection,” Annu. Rev. Fluid Mech. 26, 169 !1994". 35 S. Kida, M. Takaoka, and F. Hussain, “Formation of head tail structure in a 2-dimensional uniform straining flow,” Phys. Fluids A 3, 2688 !1991". 36 C. B. Martin, F. S. Hover, and M. S. Triantafyllou, “Propulsive performances of a rolling and pitching wing,” Proceedings of the UUST, International Symposium on Unmanned Untethered Submersible Technology, Durham, NH !2001". Downloaded 14 Dec 2005 to 18.172.7.243. Redistribution subject to AIP license or copyright, see http://pof.aip.org/pof/copyright.jsp
© Copyright 2024