Xue-ming SHAO （邵雪明）， Xiao-long ZHANG （張曉龍）， Zhao-sheng YU （余釗圣）

State Key Laboratory of Fluid Power Transmission and Control， Department of Mechanics， Zhejiang University，Hangzhou 310027， China， E-mail: mecsxm@zju.edu.cn

?

Numerical studies of the hysteresis in locomotion of a passively pitching foil*

Xue-ming SHAO （邵雪明）， Xiao-long ZHANG （張曉龍）， Zhao-sheng YU （余釗圣）

State Key Laboratory of Fluid Power Transmission and Control， Department of Mechanics， Zhejiang University，Hangzhou 310027， China， E-mail: mecsxm@zju.edu.cn

The direct-forcing fictitious domain method is extended to simulate the locomotion of a passively pitching foil. Our study focuses on the hysteresis phenomenon that the critical frequency for the reverse of the locomotion direction of the wing in case of decreasing frequency is smaller than that in case of increasing frequency. In our simulations， the hysteresis phenomenon is produced by imposing different initial conditions at a same frequency. Our results indicate that the ratio of the heaving amplitudes of two foil edges is crucial to the direction of the foil’s horizontal motion， and the amplitude of the leading edge is generally smaller. The critical frequencies for the reverse of the locomotion direction are increased， when the foil-fluid density ratio is decreased or the spring constant is increased. The critical frequencies in the bi-stability regime also depend on the initial velocity imposed， and the hysteresis loop generally becomes larger if the initial velocities are closer to the terminal locomotion velocities of the foil.

flapping foil， passively pitching， hysteresis， fictitious domain method

Introduction

The flapping locomotion of insects， birds and fishes was extensively studied［1-3］. At the early stage，the focus is on the production of the thrust and the lift on a wing model fixed in the streamwise direction and subject to a prescribed lateral flapping motion［4-6］. In recent years， the self-propelled flight or swimming was studied. Vandenberghe et al.［7，8］experimentally investigated the dynamics of a rigid symmetric foil that flaps vertically and moves freely in the horizontal direction in a fluid. They observed the locomotion of the foil above a critical flapping frequency. The dynamics of the self-propelled body under a given vertical flapping motion were numerically studied by Alben and Shelley［9］， Lu and Liao［10］， Zhang et al.［11］and Hu and Xiao［12］.

Spagnolie et al.［13］experimentally examined the locomotion of a flapping foil which pitched passively. One edge of the foil was driven to heave vertically，and the wing was able to passively pitch with a restoring torque from a torsional spring and to translate in the horizontal direction. They observed a hysteresis phenomenon: the foil reversed its motion direction as the heaving frequency was increased or decreased， but the critical frequency in the case of decreasing frequency was smaller than that in the case of increasing frequency. The dynamics of such a model was numerically studied by Spagnolieet al.［13］， Zhang et al.［14］and Xiao et al.［15］， but the hysteresis phenomenon has not been reproduced and studied. The aim of the present study is to numerically reproduce the hysteresis phenomenon， elucidate its mechanism， and examine the effects of the density ratio and the spring constant.

1. Physical model and numerical method

1.1 Physical model

The flapping foil is modeled as a two-dimensional elliptic body （as shown in Fig.1）. The major and minor axis lengths areand， respectively. The right edge of the elliptic wing is driven to heave vertically， with the position given as

Fig.1 Schematic diagram of a two-dimensional elliptic body as a model of the foil. The right edge point is driven to move with a given vertical displacementis the pitching angle between the elliptic major axis and the horizontal plane

The foil can rotate freely， and a torsional spring is introduced to produce a restorative torque， whose amplitude is linearly proportional toθ，

1.2 Numerical method

Yu and Shao［16］proposed a direct-forcing fictitious domain method for the simulation of the particle motion in a fluid， by modifying the original fictitious domain method proposed by Glowsinki et al.［17］. The key idea of the fictitious domain method is that the interior domain of the solid body is assumed to be filled with the fluid and the inner fictitious fluid is constrained to move as a rigid body by a pseudo body force （i.e.， Lagrange multiplier）［17］. In the present study， we further extend the fictitious domain method to deal with the flapping foil with passive pitching and translational motions as introduced earlier.

1.2.1 Governing equations

The original momentum equation for the fluid field can be written as follows

and the Newton’s equations of motion for the foil are:

In the following， we deduce the fictitious domain formulation based on its key idea introduced earlier. First， we extend the fluid momentum Eq.（3） from the real fluid domainto the entire domainincludingand the foil domain， and a pseudo body forceis introduced to make the fictitious fluid inside the foil satisfy the rigid body motion constraint，

Substituting Eq.（14） into Eq.（4） yields

Substituting Eq.（15） into Eq.（5） yields

Substituting Eq.（13） into Eq.（6） yields

Then， substituting Eq.（17） into Eq.（18）， we obtain

Now we have the new governing equations for the incompressible Newtonian fluid and the foil motion:

The momentum equation Eq.（20） is obtained by substituting the constitutive equation of the Newtonian fluid into Eq.（9）. In Eq.（20），pis the fluid pressure andis the fluid viscosity. Equation （21） is the continuity equation for the incompressible fluid. The following scales are used for the non-dimensionlization:for length，for velocity，for time，for the pressure，for the pseudo body force， andfor the pseudo body force， and then the dimensionless forms of Eqs.（20）-（24） are

1.2.2 Numerical schemes

A fractional-step time scheme is used to decouple the system （25）-（29） into the following two sub-problems. The fluid sub-problem forand

A finite-difference-based projection method on a homogeneous half-staggered grid is used for the solution of the above fluid subsystem Eqs.（30）， （31）. All spatial derivatives are discretized with the secondorder central difference scheme.

The foil-motion sub-problem forand

In the following， we re-formulate the foil-motion subproblem so that it can be solved explicitly. Integrating thecomponent equation of Eq.（32） over the foil domain and adding the resulting equation to Eq.（33），we have

Considering that

and

After the solution of the fluid sub-problem， the fluid velocityis obtained， and then all the right-hand side terms in Eqs.（35） and （38） are known quantities，therefore， we can compute the translational velocity and the angular velocity of the foil from Eqs.（35） and（38） explicitly.

The above reformulation of the foil-motion subproblem not only simplifies the computation by changing a coupled system into a de-coupled system without any sacrifice of the stability， but also circumvents the singularity atin the original subproblem of Eqs.（33）， （34）. Finally， we note that a bilinear function is used to interpolate the fluid velocity from the Eulerian nodes to the Lagrangian nodes， and to distribute the pseudo body force from the Lagrangian nodes to the Eulerian nodes［16］.

The direct-forcing fictitious domain method was previously proposed by Yu and Shao［16］for the particle-laden flow， but the extension to the foil problem with the prescribed motion at one point in the lateral direction is non-trivial， as seen from the above derivation. The extension of the fictitious domain method might be a meritorious point of this study.

In our simulations， the computational domain is， and the grid number is 1 024×512， thus the mesh size isThe dimensionless time step isThe zero-velocity condition is imposed on the top and bottom lateral boundaries and the upstream boundary， and the outflow condition is imposed on the downstream boundary. The streamwise size ofof the domain is definitely not long enough for the foil to reach the steady state， so we move the computational domain with the foil in thedirection so that the foil center is always located in a mesh at the center of the domain. The test result shows that further increase in the lateral size of the domain does not make a significant difference.

1.2.3 Validation

The locomotion of a foil is driven to heave vertically without pitching motion， which is a classic selfpropulsion problem. Therefore， it is chosen to validate the accuracy of our code. For this problem， two Reynolds numbers are defined

Fig.2 Comparison between the present results and previous data for the locomotion Reynolds numberversus the frequency Reynolds numberThe foil is driven to heave vertically as a whole and moves freely in the horizontal direction forand ?

2. Results and discussions

For the problem of interest， the dimensionless physical parameters include the Reynolds numberthe density ratio， the spring constantand the amplitudeThe effect of the frequency is implicitly reflected in the values of the Reynolds number and the dimensionless spring constant. Nevertheless the hysteresis phenomenon is related to the change of the frequency， therefore， it is better to directly examine the effect of the frequency. Thus， following Spagnolie et al.［13］， we introduce a dimensionless frequencywhich is the relative frequency with respect towithbeing determined fromSimilarly，is a spring constant normalized with. With the above definitions， we have the relationships:and. The effect of theamplitudeis not considered. Throughout the study， we fix the following parameters:We first study the hysteresis phenomenon forandin Sections 2.1， and then examine the effects ofandin Section 2.2.

2.1 Hysteresis phenomenon

In the experiment of Spagnolie et al.［13］， the hysteresis phenomenon was produced by continuously increasing and decreasing the heaving frequency，respectively. The frequency in the hysteresis loop sees two stable states with different foil locomotion directions. Essentially， the bi-stable states are caused by the different initial conditions for the foil motion. Therefore， in the present study， at the initial stage，instead of releasing the foil immediately， we enforce the foil to move at a constant horizontal velocitytill the system reaches the stable state （we actually set the simulation time of this initial stage to be 40， which we found is sufficiently long） and then release it. We taketo simulate the experimental condition of increasing frequency andfor decreasing frequency， since outside the hysteresis loop， the foil moves forward at a low frequency and backward at a high frequency.

Fig.3 The foil’s mean horizontal velocityas a function of heaving frequency under two different initial conditions atand

The results of the mean horizontal velocity of the foil at the stable （quasi-steady） stateas a function of the heaving frequency for two different initial conditions （i.e.and ˉ1） andandare shown in Fig.3. In the case， the foil reverses its locomotion velocity direction at the critical frequency， whereas in the casethe foil reverses its velocity direction at the critical frequency. There exists a bi-stable regime in the frequency range of. Thus， we reproduce qualitatively the hysteresis phenomenon observed in the experiment of Spagnolie et al.［13］. Outside the hysteresis loop， the computed velocities under two different initial conditions are not exactly the same， which is probably a numerical artifact. Nevertheless， the difference is not significant.

Fig.4 Time-development forzndThe foil is release at

In the bi-stable regime aforementioned， the foil’s locomotion direction at the stable state depends on its initial condition. The evolutions of the lateral displacement of the foil left edge （or point） and the pitching angle in the caseare shown in Fig.4. The foil reaches the stable state after the initial transient stage of a few oscillating periods when its locomotion speed is constrained to be unity. After the foil’s horizontal motion is released at， the foil’ motion evolves into another stable state with the motion amplitudes not differing much from those at the constrained stage. Similar behavior is observed in the case ofWe also find that if the foil’s horizontal motion is released with the initial horizontal velocitywithout being driven for a while， the hysteresis phenomenon cannot be observed. This indicates that the well-developed flow fields are important for the occurrence of the bi-stable states.When the flow fields are well developed， the foil can extract the energy from the flow and sustain its original motion state. This is the mechanism for the bi-stability in the present problem， as also for other bi-stability phenomena associated with flow-structure interactions such as the flow over a flag or filament［18-20］. In the previous numerical studies of Spagnolieet al.［13］and Zhang et al.［14］， the same initial condition （presumably zero velocity） was given for each frequency and therefore the bi-stability was not observed.

Fig.5 Vorticity contours of one period at the quasi-steady state

Figure 5 shows the vorticity contours atandof one period in the quasisteady state in the casesandrespectively. The reverse Kármán vortex street in the wake of the foil and the large leading edge vortex are observed in both cases， and there are no essential differences in the flow structures between the two cases except that the motion direction is opposite and the leading vortex in the forward-moving case of1 appears a little stronger.

Fig.6 Amplitude ratiosat different frequencies forandThe amplitude at the right edgeis fixed as 0.5， i.e.

Figure 6 shows the ratio of the amplitudes of the vertical motions of the left and right edgesas a function of frequency. Note that the amplitude at the right edgeis fixed as 0.5， i.e.There also exists a hysteresis loop for the amplitude ratio. Outside the hysteresis loop， the amplitude of the left edge is larger than that of the right edge at the frequencies below the lower critical frequency， corresponding to the forward motion， whereas the amplitude of the left edge is smaller than that of the right edge at the frequencies above the upper critical frequency，corresponding to the backward motion. In other words，outside the bi-stable regime， the edge with smaller heaving amplitude plays a role as the leading edge. In the bi-stable frequency regime， in the forward motion，the amplitude of the left edge can be smaller than that of the right edge， but is clearly larger than that of the left edge in case of backward motion. Therefore， we conclude that the heaving amplitude of the edge is crucial to the motion direction of the foil， and the edge with smaller amplitude tends to serve as the leading edge. Of course， the amplitude of the left edge is a passive outcome， and the motion direction of the foil should be determined by the control parameters and the initial conditions. In the following， we study the effects of the density ratio and the spring constant.

2.2 Effects of the density ratio and the spring constant

We first examine the effects of the density ratio on the hysteresis. The results in the bi-stable frequency regimes and the amplitude ratioas a function of frequency in the caseandare shown in Figs.7 and 8， respectively. Figure 7 indicates that both lower and upper critical frequencies increase with the decrease of the density ratio. From Fig.8， our earlier observations on the heaving amplitude of the left foil edge hold true for all density ratios: the amplitude of the left foil edge is always smaller for the backward motion， and larger for the forward motion below the lower critical frequency， compared to that of the right edge.

Fig.7 Bi-stable regimes for different density ratios andThe regime is marked by the line with the critical frequencies indicated

Fig.8 Amplitude ratiosas a function of frequency for different density ratiosand

Fig.9 Bi-stable regimes for different spring constants， and. The regime is marked by the line with the critical frequencies indicated

The results in the bi-stable regimes forare shown in Fig.9， which shows that both lower and upper critical frequencies increase with the decrease of the spring constant. In addition， from Figs.7 and 9， when the critical frequencies become larger， the frequency interval in the bi-stable regime generally becomes larger.

Zhang et al.［14］investigated numerically the locomotion of a passively flapping flat plate， and observed that the locomotion direction of the plate depended solely on the ratio of the natural frequency of the systemand the heaving frequency， The plate moved forward whenand backward when. The natural frequency is defined as follows

Fig.10 Amplitude ratiosas a function of frequency for different velocities imposed at initial stage.and

Finally we explain why the frequency interval for the bi-stable regime is expanded for a lower density ratio and particularly for a higher spring coefficient， as shown in Figs.7 and 9. The reason may be related to the effect of the initial velocity imposed. Figure 10 shows that the hysteresis loop is shrinked when the initial velocitiesare changed fromtoin the case ofandwhich could be explained by the fact that the initial velocities ofare closer to the terminal locomotion velocities of the foil （being aroundin the bi-stability frequency regime as seen from Fig.3 or Fig.11）. From Fig.11， the locomotion velocity of the foil whenapproaches unity （i.e.， the value of the initial velocity）， as the spring coefficient is increased from 500 to 1 500，which is a plausible reason for the larger bi-stability frequency interval at a larger

Fig.11 The foil’s mean horizontal velocitiesas a function of heaving frequency for different spring coefficients atand

3. Conclusion

We have extended the direct-forcing fictitious domain method to simulate the locomotion of a passively pitching foil. The hysteresis phenomenon is produced by imposing different initial conditions at the same frequency. Our results indicate that the ratio of the amplitudes of two foil edges is closely related to the locomotion direction of the foil， and the edge with smaller heaving amplitude tends to play a role as the leading edge. In addition， the effects of the density ratio and the spring constant on the hysteresis are investigated， and our results show that the critical frequencies for the reverse of the locomotion direction are increased， as the foil-fluid density ratio is decreased or the spring constant is increased， as is consistent with the observations of Zhang et al.［14］for a passively pitching flat plate. The critical frequencies in the bistability regime also depend on the initial velocity imposed， and the hysteresis loop generally becomes larger if the initial velocities are closer to the terminal locomotion velocities of the foil.

References

［1］ LEHMANN F. O. The mechanisms of lift enhancement ininsect flight［J］. Naturwissenschaften， 2004， 91（3）: 101-122.

［2］ SANE S. P. The aerodynamics of insect flight［J］. Journal of Experimental Biology， 2003， 206（23）: 4191-4208.

［3］ WANG Z. J. Dissecting insect flight［J］. Annual Review of Fluid Mechanics， 2005， 37（1）: 183-210.

［4］ TRIANTAFYLLOU M. S.， TRIANTAFYLLOU G. S. and YUE D. K. P. Hydrodynamics of fishlike swimming［J］. Annual Review of Fluid Mechanics， 2000， 32: 33-53.

［5］ TRIANTAFYLLOU M. S.， TECHET A. H. and HOVER F. S. Review of experimental work in biomimetic foils［C］. IEEE Journal of Oceanic Engineering. Durham， USA，2004， 29: 585-594.

［6］ LIU H. Simulation-based biological fluid dynamics in animal locomotion［J］. Applied Mechanics Reviews， 2005，58（4）: 269-282.

［7］ VANDENBERGHE N.， CHILDRESS S. and ZHANG J. On unidirectional flight of a free flapping wing［J］. Physics of Fluids， 2006， 18（1）: 014102.

［8］ VANDENBERGHE N.， ZHANG J. and CHILDRESS S. Symmetry breaking leads to forward flapping flight［J］. Journal of Fluid Mechanics， 2004， 506: 147-155.

［9］ ALBEN S.， SHELLEY M. Coherent locomotion as an attracting state for a free flapping body［J］. Proceedings of the National Academy of Sciences of the United States of America， 2005， 102（32）: 11163-11166.

［10］ LU X. Y.， LIAO Q. Dynamic responses of a two-dimensional flapping foil motion［J］. Physics of Fluids， 2006，18（9）: 098104.

［11］ ZHANG X.， NI S. and WANG S. et al. Effects of geometric shape on the hydrodynamics of a self-propelled flapping foil［J］. Physics of Fluids， 2009， 21（10）: 103302.

［12］ HU J. X.， XIAO Q. Three dimensional effects on the translational locomotion of a passive heaving wing［J］. Journal of Fluid and Structures， 2013， 46: 77-88.

［13］ SPAGNOLIE S.， MORET L. and SHELLEY M. J. et al. Surprising behaviors in flapping locomotion with passive pitching［J］. Physics of Fluids， 2010， 22（4）: 041903.

［14］ ZHANG J.， LIU S. N. and LU X. Y. Locomotion of a passively flapping flat plate［J］. Journal of Fluid Mechanics， 2010， 659: 43-68.

［15］ XIAO Q.， HU J. X. and LIU H. Effect of torsional stiffness and inertia on the dynamics of low aspect ratio flapping wings［J］. Bioinspiration and Biomimetics， 2014，9（1）: 016008.

［16］ YU Z. S.， SHAO X. M. A direct-forcing fictitious domain method for particulate flows［J］. Journal of Computational Physics， 2007， 227（1）: 292-314.

［17］ GLOWINSKI R.， PAN T. W. and HESLA T. I. et al. A distributed Lagrange multiplier/fictitious domain method for particulate flows［J］. International Journal of Multiphase Flow， 1999， 25（5）: 755-794.

［18］ ZHANG J.， CHILDRESS S. and LIBCHABER A. et al. Flexible filaments in a flowing soap film as a form for one-dimensional flags in a two-dimensional wind［J］. Nature， 2000， 408: 835-839.

［19］ YU Z. S.， WANG Y. and SHAO X. M. Numerical simulations of the flapping of a three-dimensional flexible plate in uniform flow［J］. Journal of Sound and Vibration，2012， 331（20）: 4448-4463.

［20］ TANG Chao， LU Xi-yun. Self-propulsion of a threedimensional flapping flexible plate［J］. Journal of Hydrodynamics， 2016， 28（1）: 1-9.

December 17， 2014， Revised June 15， 2015）

* Project supported by the National Natural Science Foundation of China （Grant No. 11372275）， the Program for New Century Excellent Talents in University.

Biography: Xue-ming SHAO （1972-）， Male， Ph. D.， Professor

Zhao-sheng YU，

E-mail: yuzhaosheng@zju.edu.cn