OOn the forced flow around a flapping foil.
F. Mandujano and C. M´alaga.October 8, 2018
Abstract
The two dimensional incompressible viscous flow past a flapping foilimmersed in a uniform stream is studied numerically. Numerical simu-lations were performed using a Lattice-Boltzmann model for moderateReynolds numbers. The computation of the hydrodynamic force on thefoil is related to the the wake structure. In particular, when the foil’s cen-tre of mass is fixed in space, numerical results suggest a relation betweendrag coefficient behaviour and the flapping frequency which determinesthe transition from the von K´arm´an (vKm) to the inverted von K´arm´anwake. Beyond the inverted vKm transition the foil was released. Upstreamswimming was observed at high enough flapping frequencies. Computedhydrodynamic forces suggest the propulsion mechanism for the swimmingfoil.
The flow around a flapping foil has been studied in a variety of conditions. Foilson a stream, flapping as a result of the hydrodynamic forces acting on them,have been studied in connection to energy extraction processes (Wu et al.,2015; Wu, 1972). The wake of foils with an imposed flapping and translationalmotion, and its relation with the reacting forces on them has deserved attentiondue to its relation to flying and swimming (Anderson et al., 1998; Dong et al.,2006; Triantafyllou et al., 2000; Taylor et al., 2003). Little has been done inthe case of an imposed flapping on a foil free to swim, although work can befound related to other type of swimmers at moderate Reynolds numbers (Kernand Koumoutsakos, 2006; Gemmell et al., 2015; Chisholm et al., 2016). Herewe study numerically a two-dimensional foil with an imposed flapping motion,immersed in a stream, that can be released and become free of translationalmotion. Both phenomena were studied and compared for moderate Reynoldsnumbers, where inertia and viscous forces are important.We used the lattice Boltzmann model (LBM) of He et al. (1998), with aprocedure to include immersed moving boundaries of arbitrary shape subject tohydrodynamic forces (Mandujano and Rechtman, 2008). The LBM is an effi-cient algorithm to approximate solutions to the Navier- Stokes equations whenimplemented on massively parallel architectures, where it can handle simula-tions of high resolution and large computational domains (Clausen et al., 2010).For the problem in hand the LBM seemed a good option, since wakes behindobjects may show interesting behaviour far away behind the object (Vorovieff1 a r X i v : . [ phy s i c s . f l u - dyn ] A ug D h
Figure 1: Schematic illustration of the flow around a flapping foil.et al., 2002), and so a large domain relative to the size of the foil was desirable,along with enough resolution to observe vortex formation and evolution.The article is divided as follows, in section 2 the general problem is presentedand the numerical method used is described briefly in section 3. In section 4 theresults are presented and discussed. Conclusions are summarised in in section5.
Consider an unbounded two dimensional and incompressible flow of a viscousfluid around a rigid flapping foil made by a semicircle, with diameter D , inter-sected with an isosceles triangle of high h , as shown in figure 1. Two cases areconsidered, a foil with a point fixed in space and a foil free to move. In bothcases the flapping motion consist of an oscillation with respect to the centreof the semicircle, which is fixed in space in one case and free in the other. Tosimplify the calculations, it is assumed that the foil mass density distribution, ρ s , is such that its centre of mass is at the centre of the semicircle.The governing equations for the fluid, with mass density ρ f and kinematicviscosity ν , are the Navier-Stokes equations given by ∇ · u = 0 , (1) ∂ u ∂t + u · ∇ u = − ρ f ∇ P + ν ∇ u , (2)where u ( r , t ) and P ( r , t ) are the velocity and pressure fields, respectively. Con-sider no-slip condition at the foil’s surface. Constant pressure and a uniformand horizontal velocity field U are imposed far from the foil (see figure 1).The rigid foil is forced to rotate around the centre of the semicircle (also itscentre of mass) so that the angle of its axis of symmetry with respect to thehorizontal coordinate evolve with θ ( t ) = θ sin ωt (see Figure 1). The velocityof a point at the foil’s surface is v s ( r s , t ) = V ( t ) + ˙ θ ( t )( r s − R ), where R is theposition of the centre of mass, V its velocity and r s is the position of a pointat the foil’s surface. Therefore, the no-slip boundary condition takes the form u ( r s , t ) = v s ( r s , t ). In the case of the free flapping foil, the motion of its centreof mass is obtained by solving md V /dt = F , where m is the mass of the foil.The hydrodynamical force F acting on the foil, is obtained directly from thecomputed flow fields. 2o work with non-dimensional variables we choose D as the characteristiclength, U as the characteristic velocity and the fluid density ρ f . With thischoice of scaling the flow regime is determined by five dimensionless parameters:the Reynolds number Re = U D/ν , the Strouhal number St = ωD/ πU thedimensionless amplitude A = 2 θ h/D , the ratio between the chord of the foil andthe circle’s diameter C = 1 / h/D , and the fluid-foil’s mass density ratio µ = ρ f /ρ s which is set to unity to work with a non-buoyant foil. The hydrodynamicforce is scaled with ρ f U D/
2. We define the drag C D and lift C L coefficients asthe horizontal and vertical components of 2 F /ρ f U D , respectively. Notice thatnegative values of C D correspond to a thrust dominated hydrodynamic force.Solutions to equations (2) and (1) for the pressure and velocity fields areapproximated using a lattice-Boltzmann model. The numerical procedures aredescribed in the following section. To compute the flow around the flapping foil a two dimensional, nine neighbours(D2Q9) lattice-Boltzmann model was used (He et al., 1998). The proposed algo-rithm, that includes moving immersed boundaries, was validated in a previouswork concerning a problem of a moving cylinder in a convective flow (Mandujanoand Rechtman, 2008).In this method space is discretised using a square lattice. Lattice spacing aswell as time steps can be conveniently set to unity. The state of the fluid at thenode with vector position r at time t , is described by the particle distributionfunction f k ( r , t ) that evolves in time and space according to f k ( r + e k , t + 1) = f k ( r , t ) − τ (cid:104) f k ( r , t ) − f ( eq ) k ( r , t ) (cid:105) , (3)where τ is a relaxation time related to the fluid kinematic viscosity ν = ( τ − / /
3. The distribution function f ( eq ) k is given by f ( eq ) k ( r , t ) = w k ρ (cid:18) e k · u + 92 ( e k · u ) − u (cid:19) , (4)which corresponds to a discrete Maxwell distribution function for thermal equi-librium. In the above expressions, the macroscopic density and velocity fieldsare computed using ρ ( r , t ) = (cid:88) k f k ( r , t ) and ρ u ( r , t ) = (cid:88) k e k f k ( r , t ) . The microscopic set of velocities e k is given by e = (0 , , e k = (cos( π ( k − / , sin( π ( k − / k = 1 , . . . , , e k = √ π ( k − / / , sin( π ( k − / / k = 5 , . . . , , where w = 4 / w k = 1 / k = 1 , . . . , w k = 1 /
36 for k = 5 , . . . , e k , expression (3) is al-ways evaluated at lattice points. It is well known that the above procedure ap-proximates solutions to the Navier-Stokes equations in the limit of small Machnumbers (He et al., 1998). 3quation (3) provide an explicit algorithm for updating all the distributionfunctions f k at a given node in the lattice, as long as its 8 nearest neighbouringnodes are inside the fluid domain. For nodes adjacent to a solid wall, thedistribution functions coming from neighbouring nodes outside the fluid domainmust be provided as a boundary condition for the method. We choose to adoptthe set of boundary conditions proposed by Guo and Zheng (2002) for curvedrigid walls.The force and torque acting on the body are computed using the momentum-exchange method of Mei et al. (2002). For the free flapping foil, as rotationalmotion is imposed, only the force is used to compute the foils translationalmotion using a forward Euler integration in time at each time step. This is aparticular case of the scheme described in Mandujano and Rechtman (2008) forthe motion of a free particle in a convective flow using the lattice-Boltzmannmodel.The computational domain was a rectangular lattice of 10000 × × h of 280 to 600 nodes, a width D of 80 nodes (about a 1 /
30 of the width andlength of the domain), and was placed at 3000 nodes from the left boundary ofthe domain (5 foil lengths). To simulate conditions far from the foil, the velocitywas set to U at the left boundary following the procedure of Guo and Zheng(2002). On the rest of the boundaries, the normal components of the velocitygradients were set to zero setting the unknown velocity at the boundary equalto that of the adjacent node normal to the wall.The numerical scheme was implemented to run in parallel in GPU’s due tothe large number of nodes involved in the simulations. Typically, two days areneeded to obtain a periodic flow (3 × time steps) of a free foil swimmingupstream running on an Nvidia R (cid:13) Tesla K40 processor.
The foil profile geometry was chosen to compare with experiments performedby Schnipper et al. (2009) in a soap films and by Godoy-Diana et al. (2008) ina wind tunnel. Numerical results confirm the relation between the drag-trusttransition and the behaviour of the von Karman (vKm) vortex street.
Simulations start with the foil at rest facing a uniform flow at the left boundaryto allow for a wake to form before flapping. After 5 × time steps, the foilstarts to flap increasing its amplitude exponentially in time until it reaches A . The simulations were performed for a series of values of St ∈ [0 . , . A ∈ [0 . , Re ∈ [100 , C took values 4 and 6.Following Schnipper et al. (2009), in the first set of simulations Re = 220 and C = 6. The simulations performed show the formation of patterns in approxi-mately the same parameter region as reported there (see Figure 2); transitionsbetween vKm and inverted vKm wake, formation of 2P, 4P, 2P+2S, 4P+2Swakes (following the nomenclature of Williamson and Govardhan (2004)) andthe transitions between them. Quantitative comparison can not be expected asfilm flow is not a two dimensional phenomena.4igure 2: Vorticity distribution for a fixed flapping foil for Re = 220. Coloursrepresent positive (red) and negative (blue) vorticity. a) vKm wake for St = 0 . A = 0 .
98. b) inverted vKm wake for St = 0 .
12 and A = 2. c) 2P wakefor St = 0 .
085 and A = 1 .
4. d) 2P+2S wake for St = 0 .
05 and A = 1 .
4. e) 4Pwake for St = 0 .
035 and A = 1 .
34. f) 4P+2S wake for St = 0 .
025 and A = 1 . C L oscillates with the flapping frequency, corre-sponding to a Fourier mode n = 1. The drag coefficient C D oscillates withtwice the flapping frequency and has modes n = 0 and n = 2. Notice that C D minima (maximum thrust) is obtained after extremal values of the angle ofattack θ are reached.The flapping frequency sets the hydrodynamic force frequencies and higherharmonics. For patterns with combinations of 2P, 4P and 2S, Fourier coefficientscorresponding to higher harmonics become important. In Figure 2-(d)-(f), thelift and torque coefficients are composed by the first few odd modes while C D modes are even (see Figure 3-(a)). In the combination 2P+2S, C D have n = 2and 4 wile C L is composed by n = 1 and 3 modes. The 4P wake includesmodes n = 6 and n = 7 and the 4P+2S wake has n = 8 and n = 9 modes. C D always showed even modes while C L modes were all odd. There is a stronghorizontal push with every stroke and higher harmonics are produced at lowflapping frequencies where many vortices are shed during each stroke.When St = 0 .
12 and A ∈ [0 . ,
3] there is a transition between the vKmand the inverted vKm wake, known to be relate with the drag-thrust transi-tion (K´arm´an and Burgers, 1935; Schnipper et al., 2009; Godoy-Diana et al.,2008). Our LBM simulations confirm that there is an inversion of the vKmwake in this interval accompanied by a mostly sinusoidal C D with a minimumvalue that crosses zero at the transition, while retaining a positive mean value(see Figure 3-(b)).In Figure 4 it is shown the mean drag coefficient (cid:104) C D (cid:105) and the differencebetween (cid:104) C D (cid:105) and the norm of its Fourier decomposition C DA , as a function ofthe product StA . For a sinusoidal C D , (cid:104) C D (cid:105) − C DA is the minimum value of C D . 5 (2 π St) t C D a) (2 π St) t C D b) θ (t) Figure 3: Drag coefficient as function of time for some of the vorticity distri-butions of Figure 2. a) 4P wake for St = 0 .
035 and A = 1 .
34 and 4P+2S wakefor St = 0 .
025 and A = 1 . St = 0 .
12 and A = 0 .
98 andinverted vKm wake for St = 0 .
12 and A = 2. The angle of attack θ ( t ) of thefoil is shown in (b) for the inverted vKm case. A St 〈 C D 〉 a) A=1.4, Re=100A=1.7, Re=100A=2.0, Re=100A=1.4, Re=125A=1.4, Re=150A=1.4, Re=200St=0.12, Re=220St=0.22, Re=255
A St 〈 C D 〉 − C D A b) A=1.4, Re=100A=1.7, Re=100A=2.0, Re=100A=1.4, Re=125A=1.4. Re=150A=1.4, Re=200St=0.12, Re=220St=0.22, Re=255Fit
Figure 4: a) Mean drag coefficient (cid:104) C D (cid:105) and b) (cid:104) C D (cid:105) − C DA as functions of the StA . The lines connecting points mean that either A or St were kept constantas marked. The last two experiments correspond to the parameters used inSchnipper et al. (2009) and Godoy-Diana et al. (2008), respectively. C = 4when Re = 255 and C = 6 in the rest of the cases.6igure 5: Evolution of the wake for Re = 100 and C = 6. a) StA = 0 .
21, b)
StA = 0 .
42, c)
StA = 0 .
49, d)
StA = 0 .
56, e)
StA = 0 .
21, f)
StA = 0 .
42, g)
StA = 0 .
51 and h)
StA = 0 .
6. In the left column the dimensionless amplitudewas keep fixed at A = 1 . S = 0 .
3. Figures show the region of interest within thecomputational domain. Colours represent positive (red) and negative (blue)vorticity.Figure 4(a) shows that the mean drag decreases monotonically and becomesnegative, as
StA increases for the cases where
Re < Re show that (cid:104) C D (cid:105) decreases and can be negative but starts grow-ing beyond certain value of StA . These observations suggest that for a givenimposed flow U , to which Re is related, there is a critical flapping frequencybeyond which the mean drag starts growing. We could not confirm this ideafor the cases where Re <
150 because the numerical model became unstable forhigh flapping frequencies. Figure 4(b) shows that the minimum of C D crosseszero at StA ∼ .
16 regardless of the values of Re and C used. The numericalexperiments show that below this value of StA vortices at the wake correspondsto the vKm wake type, but can be combinations of 2S and 2P patterns when St is small. The inverted vKm is found when StA > .
16 (see figures 2(b), 5(a)and 5(e)). Figure 4(b) shows that all cases collapse approximately on a singlecurve given by the fit (cid:104) C D (cid:105) − C DA = − . StA ) . + 0 . . (5)For Re = 220 and C = 6, the wake remains symmetric with a positive valueof (cid:104) C D (cid:105) for all values of StA explored by Schnipper et al. (2009). The evolution7f the wake as function of A and St , starting form a vKm wake, is shown infigure 5 for Re = 100. Beyond the inverted vKm transition, for StA > .
28, thewake is deflected from the horizontal centre line and (cid:104) C D (cid:105) becomes negative.The results shown that the distance between the tail of the foil and the point ofdeflection seems to decrease as StA increases, in agreement with the numericalsimulations made by Deng and Caulfield (2015) (see figures 5(b), (c), (f) and(g)). For Re = 255, the performed simulations shown that the wake deflectswhen StA > .
23, in good agreement with the 3D experimental observationsof Godoy-Diana et al. (2008) and the numerical simulations of He et al. (2012)and Deng and Caulfield (2015).The direction of deflection is related with the initial condition, and a down-wards deflection was found when changing the sign of the angle θ → − θ .With the set of parameters used by Godoy-Diana et al. (2008), the simulationsshowed that the position and angle of deflection can change with time after sev-eral vortex shedding periods of observation, this long time behaviour has beenobserved by Lewin and Haj-Hariri (2003).The results with Re = 100 at values of StA ∼ . C L showing only even modes while C D showing contribution from both even and odd modes. The simulations for the free foil were made using similar initial conditions as inthe fixed case. After 5 × iterations, where the foil was fixed without flapping,the foils is released and starts flapping with the same exponential increase inamplitude as the fixed case. The values of StA ∈ [0 . , .
5] were chosen incorrespondence to patterns beyond the transition to the inverted vKm wake forthe fixed case.A stationary symmetric flow without vortex shedding was observed beforethe foil was released in all cases. After released, a transient time was observedwhere the hydrodynamic force changed from a constant value to a periodicfunction of time. During this time interval, the centre of mass accelerateduntil it reached an almost uniform motion in the horizontal direction with smalltransversal oscillations of the order of D (see figures 6(a) and (d)). The centre ofmass deviates vertically a small distance from the horizontal centre line duringthe transient, probably related to the deflection of the wake.The centre of mass velocity is a periodic function of time, the horizontalcomponent oscillates around a constant value that is a fraction of the velocityof the free stream ( U = 1 in dimensionless units) with an amplitude that is neg-ligible, hence the swimming velocity is practically constant as shown in figures6 (a) and (b). The observed amplitude of oscillation of C L is ten times biggerthan that of C D , which produces a more appreciable oscillation in the verticaldirection (see Figures 6 (d)-(f)).As StA is increased, V x diminishes until it becomes negative and the foilswims upstream. The mean drag becomes very small, (cid:104) C D (cid:105) ∼ − in all casesobserved, which is consistent with an almost uniform motion observed at longtimes in other swimmers (Chisholm et al., 2016; Kern and Koumoutsakos, 2006).Figure 6 (c) shows that minimum values of C D (maximum thrust) are reachedlittle after extremal values of the angle of attack of the foil, suggesting that8 (2 π St) t C D
40 40.5 41 41.5 42−505 (2 π St) t C L (2 π St) t V x (2 π St) t V y (2 π St) t X (2 π St) t Y a)b)c) d)e)f) Figure 6: Dynamics of the centre of mass of the foil for A = 1 .
7, a) horizontalposition X ( t ), b) horizontal velocity V x ( t ), c) Drag coefficient C D ( t ), d) verticalposition Y ( t ) , e) vertical velocity V y ( t ) and f) Lift coefficient C L ( t ). (Blue) St = 0 .
2, (red) St = 0 . St = 0 .
45 . The black line in (c) shows theangle of attack θ ( t ).swimming is not produced by pushing fluid but by a suction mechanism resem-blant of that observed by Gemmell et al. (2015) in efficient animal swimming.The wake behind the foil in this region corresponds to a inverted asymmetricvKm wake. When V x is positive (figure 7-(a)) the angle of deflection of thewake seems to be smaller than those observed in the fixed foil cases of similar StA values (see figures 2 (b) and (f)). As
StA is increased and V x becomesnegative (figure 7-(b)), the deflection of the wake is swept away and becomessymmetric and much longer than the wakes of foils unable to swim upstream.The parameter values for the free foil swimming upstream correspond to theunstable wake for the fixed case (see figure 5-(h)), a P+S mode as seen in figure8(a). Figure 8 shows how vortical structures at the foil’s surface produce by fixedfoil of unsteady wake are larger than those produced by the upstream swimmingfoil. Probably, the momentum produced by the foil, enough to swim upstreamwhen released, competes with that of the stream to produce the instability. In this work numerical results on the two-dimensional incompressible viscousflow around a flapping foil are presented. Both a fixed and a free flapping foilwere considered. The flow field was computed using a lattice Boltzmann model9 ) b)
Figure 7: Vorticity distribution for a moving flapping foil for C = 6, A = 2 . Re = 100. Colours represent positive (red) and negative (blue) vorticity. a) St = 0 . V x >
0. b) St = 0 . V x <
0, the foil swims upstream. Figuresshow the complete computational domain. a) b)
Figure 8: Close up view of the vorticity generated by the foil for Re = 100, St = 0 .
3, and A = 2. (a) Fixed foil of figure 5-(h) that shows a P+S wake. (b)Free foil of figure 7-(b) that shows an inverted vKm wake.10hat includes a procedure for moving internal boundaries with arbitrary geome-try. Thus allowing the computation of the translational motion of a free foil dueto hydrodynamic forces. Numerical results were compared with experiments onsoap films (Schnipper et al., 2009) and on wind tunnels (Godoy-Diana et al.,2008).The numerical simulations show that the transition between the vKm wakeand the inverted vKm wake on a fixed flapping foil is found when the flappingStrouhal number StA ∼ .
16 and coincides with the minimum value of the hor-izontal hydrodynamic force on the foil crossing zero. When St (cid:46) .
16 the wakecan be a 2S (vKm) or 2P. Above this value the wake corresponds to an invertedvKm that deflects as
StA increases and the mean horizontal hydrodynamic forcebecomes negative.Whenever (cid:104) C D (cid:105) is a decreasing function of StA , there is a universal be-haviour for all Re and C explored, with a power law fit given by equation (5).This result can be related with the narrow range of ASt where many mam-mals swims Taylor et al. (2003). Results suggest that above a certain valueof the flapping frequency the drag exerted on the foil starts growing probablyproducing a flapping foil unable to swim upstream if released.The cases studied for a flapping foil free of translational motion alwaysreached an almost uniform horizontal speed with a drag coefficient value prac-tically zero. As flapping frequency increases the horizontal velocity becomesnegative (upstream motion) and the deflected portion of the wake is sweptaway. Upstream swimming produces long, symmetric and inverted vKm wakes,and is observed beyond the vKm wake transition, when the fixed foil’s wakebecomes unstable. Maximum thrust is observed close to maximal values of thefoil’s angle of attack, suggesting swimming through a suction mechanism.Partial support from project UNAM-PAPIIT-IN115216 and IN115316 is ac-knowledged. Authors thank Dr. Eduardo Ramos and Dr. Ra´ul Rechtman fortheir support.
References
J. M. Anderson, K. Streitlien, D. S. Barrett, and M. S. Triantafyllou. Oscillatingfoils of high propulsive efficiency.
J. Fluid Mech. , 360:41–72, 1998.N. G. Chisholm, D. Legendre, E. Lauga, and A. S. Khair. A squirmer acrossReynolds numbers.
J. Fluid Mech. , 796:233–256, 2016.J. R. Clausen, D. A. Reasor, and C. K. Aidun. Parallel performance of a lattice-Boltzmann/finite element cellular blood flow solver on the IBM Blue Gene/Parchitecture.
Comp. Phys. Comm. , 181(6):1013–1020, 2010.J. Deng and C. P. Caulfield. Three-dimensional transition after deflection behinda flapping foil.
Phys. Rev. E , 91:043017, 2015.H. Dong, R. Mitall, and F. M. Najjar. Wake topology and hydrodynamic perfor-mance of low-aspect-ratio flapping foils.
J. Fluid Mech. , 566:309–343, 2006.11. J. Gemmell, S. P. Colin, J. H. Costello, and J. O. Dabiri. Suction-basedpropulsion as a basis for efficient animal swimming.
Nature Comm. , 6:8790,2015.R. Godoy-Diana, J.-L. Aider, and J. E. Wesfreid. Transitions in the wake of aflapping foil.
Phys. Rev. E , 77:016308, 2008.Z. Guo and C. Zheng. An extrapolation method for boundary conditions inlattice Boltzmann method.
Phys. Fluids , 14(6):2007–2010, 2002.G.-Y. He, Q. Wang, X. Zhang, and S.-G. Zhang. Numerical analysis on transi-tions and symmetry-breaking in the wake of a flapping foil.
Acta Mech. Sin. ,28(6):1551–1556, 2012.X. He, S. Chen, and G. D. Doolen. A Novel Thermal Model for the Lattice Boltz-mann Method in Incompressible Limit.
Journal of Computational Physics ,146:282–300, 1998.T. v. K´arm´an and J. M. Burgers.
General aerodynamic theory – perfect fluids .Dover Publications, 1963, 1935.S. Kern and P. Koumoutsakos. Simulations of optimized anguilliform swimming.
J. Exp. Biol. , 209:4841–4857, 2006.G. C. Lewin and H. Haj-Hariri. Modelling thrust generation of a two-dimensional heaving airfoil in a viscous flow.
J. Fluid Mech. , 492:339–362,2003.F. Mandujano and R. Rechtman. Termal levitation.
J. Fluid Mech. , 606:105–114, 2008.R. Mei, D. Yu, W. Shyy, and L. Luo. Force evaluation in the lattice Boltzmannmethod involving curved geometry.
Phys. Rev. E , 64:041203, 2002.T. Schnipper, A. Andersen, and T. Bohr. Vortex wakes of a flappig foil.
J.FluidMech. , 633:411–423, 2009.G. K. Taylor, R. L. Nudds, and A. L. R. Thomas. Flying and swimming animalscruise at a Strouhal number tuned for high power efficiency.
Nature , 425:707–711, 2003.M. S. Triantafyllou, G. S. Triantafyllou, and D. K. P. Yue. Hydrodynamics offishlike swimming.
Annu. Rev Fluid Mech. , 32:33–53, 2000.P. Vorovieff, D. Georgiev, and M. S. Ingber. Onset of the second wake: depen-dence on the Reynols number.
Phys. Fluids , 14(7):53–57, 2002.C. H. K. Williamson and R. Govardhan. Vortex-induced vibrations.
Annu. Rev.Fluid Mech , 36:413–455, 2004.J. Wu, Y. L. Chen, and N. Zhao. Role of induced vortex interaction in a semi-active flapping foil based energy harvester.
Phys. Fluids , 27:093601, 2015.T. Y. Wu. Extraction of flow energy by a wing oscillating in waves.