Forces and energetics of intermittent swimming
AActa Mech. Sin.DOI
RESEARCH PAPER
Forces and energetics of intermittent swimming
Daniel Floryan ···
Tyler Van Buren ···
Alexander J. Smits
Received: April 5, 2017 / Revised: / Accepted: © The Authors
Abstract
Experiments are reported on intermittent swim-ming motions. Water tunnel experiments on a nominallytwo-dimensional pitching foil show that the mean thrust andpower scale linearly with the duty cycle, from a value of 0.2all the way up to continuous motions, indicating that individ-ual bursts of activity in intermittent motions are independentof each other. This conclusion is corroborated by PIV flowvisualizations, which show that the main vortical structuresin the wake do not change with duty cycle. The experimentaldata also demonstrate that intermittent motions are generallyenergetically advantageous over continuous motions. Whenmetabolic energy losses are taken into account, this conclu-sion is maintained for metabolic power fractions less than1.
Keywords unsteady propulsion · burst and coast · bio-inspired Many aquatic animals, such as large sharks and seals [1]to small schooling fish [2], exhibit an intermittent swim-ming behavior, sometimes called burst-and-coast swimming.Fish practice intermittent swimming while hunting, fleeing apredator, pursuing a mate, or while starving, and exhibit awide range of ratios of burst to coast times [3]. Our primaryinterest here is to examine the potential energy benefit of in-termittent swimming in cruising conditions.In this respect, Weihs [4] developed a simple the-oretical model and showed that fish could achieveup to 50% savings in energy through intermittent
The project was supported by the US O ffi ce of Naval Research, un-der N00014-14-1-0533.D. Floryan ( (cid:12) ) · T. Van Buren · A. J. SmitsDepartment of Mechanical and Aerospace Engineering, PrincetonUniversity, Princeton, NJ 08544, USAe-mail: dfl[email protected] swimming, which was supported by observations on theswimming range of salmon. Videler and Weihs [5] thenfound that cod and saithe choose minimum and maximumswimming velocities that match predicted theoretical optima.The key parameter that di ff erentiates intermittent fromcontinuous swimming in Weihs’ model is α , the ratio of dragduring active swimming to drag during coasting. The modelshowed that intermittent swimming is potentially energeti-cally advantageous only for values of α >
1. Lighthill [6]and Webb [7] noted that the total drag on fish is 3-5 timeshigher when they are actively swimming than when they donot actuate their bodies, which corresponds to a range of α where intermittent swimming is expected to be energeticallyadvantageous. Lighthill suggested that the drag increasesduring active swimming due to the thinning of the bound-ary layer, and so boundary layer thinning may be hypothe-sized as the mechanism responsible for the energetic bene-fits of intermittent swimming; we emphasize that this is aviscous mechanism. As the author himself noted, however,this hypothesis was tentative and untested, and more recentwork suggests that boundary layer thinning during unsteadyswimming motions could lead to only about a 90% increasein drag [8].Similar simplified hydromechanic approaches were fol-lowed by Blake [9] and Chung [10], but the main focus oftheir work was on the impact of the shape of the body of thefish and how to minimize drag, factors that directly impactthe energy expenditure of intermittent motions. Numericalsimulations by Wu et al. [11] show similar energetic benefitsas the theory.To understand the consequences of intermittent swim-ming, we have conducted an experimental investigation ofintermittent propulsion by a rigid pitching foil. The perfor-mance of pitching and heaving foils in continuous motion isrelatively well understood, and Floryan et al. [12] recentlypresented a comprehensive scaling argument that collapsesthe available data on thrust and power well. In that context,our aim is to understand how intermittent swimming motionsa ff ect the forces produced, and the energy expended, and todevelop a scaling analysis to describe these changes. Weseek to answer three specific questions: (1) do intermittent a r X i v : . [ phy s i c s . f l u - dyn ] A p r D. Floryan, T. Van Buren, and A. J. Smits
RC motor+ encoderLoad cellFlow
Pitch
Laser sheetPIV camera90 o mirror Figure 1: Experimental setup. motions reduce the energy needed to traverse a given dis-tance if metabolic losses are ignored, (2) how does this con-clusion change when metabolic energy losses are taken intoaccount, and (3) do intermittent motions reduce the energyneeded to traverse a given distance when the average speedis fixed?
Experiments were conducted using a rigid pitching foil ina recirculating free-surface water channel, as shown in fig-ure 1. The water channel test section is 0.46 m wide, 0.3 mdeep, and 2.44 m long. Ba ffl es upstream and downstream ofthe foil minimized surface waves. The free-stream velocity, U ∞ , was fixed at 100 mm / s for performance testing and 60mm / s for wake measurements.The propulsor was a two-dimensional teardrop foil witha chord of c =
80 mm, maximum thickness of 8 mm, andspan of s =
279 mm. The performance was measured by asix component force / torque sensor (ATI Mini40), with forceand torque resolutions of 5 × − N and 1 . × − N · m,respectively. The pitching motions were generated by an RCmotor (Hitec HS-8370TH) monitored via a separate encoder.Pitching motions were sinusoidal, and the amplitudes werevaried from θ = ◦ to 15 ◦ every 5 ◦ . For intermittent mo-tions, duty cycles of ∆ = . ∆ = θ = f = S t = f c sin ( θ ) / U ∞ , range from 0 .
05 to 0.4. Each trialconsisted of 30 cycles, and the data were averaged over themiddle 20 cycles. Each trial was repeated a minimum of 3times to ensure repeatability and reduce uncertainty.The wake velocity measurements were taken at thecenter-span of the foil with particle image velocimetry (PIV).Silver coated hollow ceramic spheres (Potter Industries Inc.Conduct-O-Fil AGSL150 TRD) were used to seed the flow,illuminated by a CW argon-ion laser (Spectra Physics 2020).An 8-bit monochrome CCD camera (MotionXtra HG-LE)with 1128 ×
752 resolution was used to acquire images at 25Hz. Images were processed sequentially with commercialDaVis software using spatial correlation interrogation win-dow sizes of 64 ×
64 and twice at 32 ×
32 with 50% overlap.The (cropped) measurement region covered 86 mm in thestreamwise direction and 84 mm across, with a resolution of7 vectors per 10 mm. Average and instantaneous velocityerrors are estimated to be 2.7% and 1-5%, respectively [13].
We begin by considering experiments in which the foil isfixed in the streamwise direction and pitched sinusoidally orces and energetics of intermittent swimming 3 either continuously or intermittently. The propulsive per-formance is described using the conventional definitions ofthrust and power coe ffi cients, where C T = F x ρ U ∞ sc , C P = M z ˙ θ ρ U ∞ sc . (1)Here, F x is the streamwise force, M z is the spanwise torque,˙ θ is the angular velocity of pitching, and ρ is the fluid density.In switching from swimming to non-swimming, theform of the transition is important since it has an impact onboth the thrust and power, as shown in figure 2. We thereforedefine a smoothing parameter, ξ , where ξ = ξ is the width of the Gaussian kernelrelative to the width of the active portion of the wave. As thewaveform becomes smoother, the thrust and power decreasebecause the starting and stopping accelerations decrease. Tomake comparisons meaningful, we apply a smoothing pa-rameter of ξ = . ffi cients are shown in figure 3for a range of duty cycles and pitch amplitudes. Both thethrust and power increase nonlinearly with Strouhal numberand duty cycle in all cases. In figure 4 we show the sameresults scaled with the duty cycle, which collapses all of thedata points onto a single curve. For intermittent motions,it appears to be su ffi cient to correct the time-averaged forcesuch that the averaging is only done on the active portion ofthe cycle. This procedure is similar to that used by Akoz &Moored [14], who scale the thrust for intermittent motionswith the frequency of motion, which e ff ectively accounts forduty cycle. This result has the important implication that thethrust and power generated in each actuation cycle can betreated independently. The whole is, indeed, the sum of itsparts. Figure 2: Time-averaged thrust (circular symbols) andpower coe ffi cients (square symbols) as they vary withsmoothing parameter ξ ( θ = ◦ , f = ∆ = . Figure 5 shows the instantaneous wake vortex structure atthe point in the generation cycle where the primary positivevortex is generated (labeled 1 in each figure). The contin-uous sinusoidal motion produces a typical thrust-producing“reverse von K´arm´an vortex street,” generating two opposite-sense vortices per actuation cycle. The intermittent motion,however, produces two primary vortices (1 & 2) and twosmaller secondary vortices (1’ & 2’). These secondary vor-tices are due to the rapid start and stop of the intermittentmotion. Similar vortex formations have been shown numer-ically for intermittent motions [11] and secondary vorticesare also found in rapid starting and stopping square wavemotions [15].Changing the duty cycle appears to have little e ff ect onthe location and strength of the vortices produced per cycle,suggesting that there is minimal interaction between separatecycles. This finding is consistent with the observed scalingof the thrust and power coe ffi cients with the duty cycle. We now consider a free swimmer, no longer constrained tomove at a fixed speed. Its motion is governed by Newton’sSecond Law m ˙ u = T − D , (2)where m is the mass, u is the speed, T is the thrust, and D is the drag. We make the assumptions that the swimmer canbe e ff ectively split into a drag-producing part (the body) anda thrust-producing part (the propulsor), as illustrated in fig-ure 6, and that these parts are independent of each other.We further assume a quadratic drag law [16], such that D = ρ u A w C D , (3)where ρ is the fluid density, A w is the wetted area of the body,and C D is the drag coe ffi cient. We know from [12] that thethrust produced by a pitching foil is independent of the free-stream velocity, and so we can assume that the thrust mea-sured from the fixed velocity experiments is the same as whatwould be measured if the foil were allowed to move. Fur-thermore, the appropriate velocity scale for a pitching foil is f A , where f is the frequency of the motion, and A is the am-plitude of the motion. The non-dimensionalized governingequation is then2 m ∗ A ∗ du ∗ dt ∗ = C ∗ T − C D A ∗ w u ∗ , (4)where m ∗ = m /ρ c s is the ratio of the body mass to thepropulsor added mass, A ∗ = A / c is the ratio of the ampli-tude to the chord, u ∗ = u / f A is the non-dimensional speed, t ∗ = t f is the non-dimensional time, C ∗ T = T / ρ ( f A ) sc isthe thrust coe ffi cient using the new velocity scale, C D is thedrag coe ffi cient as in (3), and A ∗ w = A w / sc is the ratio of the D. Floryan, T. Van Buren, and A. J. Smits (a) (b)
Figure 3: Time-averaged (a) thrust and (b) power coe ffi cients as functions of Strouhal number. Dark to light sym-bols represent increasing duty cycles, ranging from ∆ = . θ = ◦ (circle), 10 ◦ (square), and 15 ◦ (triangle). (a) (b) Figure 4: Time-averaged (a) thrust and (b) power coe ffi cients normalized by duty cycle as functions of Strouhalnumber. Symbols and tones as in figure 3. orces and energetics of intermittent swimming 5 Figure 5: Instantaneous wake vorticity with in-plane velocity arrows overlaid for (a) continuous motion, (b) dutycycle of ∆ = .
25, (c) ∆ = .
5, and (d) ∆ = .
75 and a Strouhal number of
S t = .
4. Small graphics in the bottomright show the pitch angle over time and the generation points of each vortex in the the actuation cycle. body’s wetted area to the propulsor’s wetted area. The pref-actor m ∗ / A ∗ is the ratio of the body’s mass to the mass of thefluid displaced by the propulsor.We find the mean speed for the free swimmer by in-putting the phase-averaged thrust measured in our experi-ments into (4) and solving numerically. We solve usingMATLAB’s built-in integrator ODE45 [17, 18], evolving thesystem forward until it settles on the limit cycle. The dragcoe ffi cient is kept constant at 0.01 [19], A ∗ varies with thedata, and we vary the mass and area ratios over some ordersof magnitude. body propulsoradded mass Figure 6: A swimmer can be simply represented asa combination of a drag-producing body and a thrust-producing propulsor. U ∗ mean to the“steady-state” speed U ∗ steady , that is, the speed found by set-ting du ∗ / dt ∗ = m ∗ and large values of the area ratio A ∗ w , which representunphysical realizations. Points with mean speed greater thanthe steady speed seem to be spurious.The relation between the calculated mean and steadyspeeds can be understood by considering the problem in thefrequency domain. Let C ∗ T = X n C n e i ω nt , u ∗ = X n u n e i ω nt (5)be the Fourier series representing the thrust and speed, re-spectively. Since the thrust and speed are real, we requirethat C − n = C n and u − n = u n , where the overbar denotesthe complex conjugate. The equation for the zeroth Fouriermode of speed is then u = C C D A ∗ w − X n , u n u n . (6)We see that the higher harmonics will act to decrease themean speed. The e ff ects of the mean thrust, drag coe ffi cient,and area ratio on the mean speed are as expected.To understand the e ff ect of the mass ratio, we lin-earize (4) about the steady speed and take the Laplacetransform. The system has a single left-half-plane pole at − bC D A ∗ w A ∗ / m ∗ , where b is a constant arising from lineariza-tion. A greater mass ratio moves the pole closer to the ori-gin, thus attenuating high frequencies. Intuitively, a greatermass ratio corresponds to a swimmer with greater inertia; thegreater inertia will cause the speed to fluctuate less about itsmean. According to (6), the attenuation of high frequencieswith greater mass ratio brings the mean speed closer to thesteady speed, explaining the observed behavior.The e ff ect of the area ratio on the mean speed can beunderstood in the same way. A smaller area ratio moves thepole closer to the origin, attenuating high frequencies. The D. Floryan, T. Van Buren, and A. J. Smits mean speed should thus be closer to the steady speed as thearea ratio is decreased.We end by noting that the mean speed only deviatesfrom the steady speed for extreme values of the parameters,unlikely to be encountered in nature, and for low values ofthe mean speed where experimental errors in our force mea-surements are relatively larger. For a wide range of the pa-rameters, estimating the mean speed to be the steady speedappears to be a reasonable approximation.
Figure 7: Mean speed versus steady speed for all fre-quencies, amplitudes, and duty cycles. Mass ratios are m ∗ = . , . , ,
10 (light to dark symbols), and arearatios are A ∗ w = φ = E / E , the ratio of the energy expendedby intermittent motion to the energy expended by continu-ous motion with the same actuation. This energy ratio canbe otherwise expressed as φ = EE = PtP t = Pd / UP d / U = PU P U , where P is the mean power, and t is the time taken to travela distance d at a mean speed U . For the problem consideredhere, d = d . Values of φ < φ > m ∗ = A ∗ w =
10, as trends should not vary withthe mass and area ratios for reasonable values, as found ear-lier. We see that intermittent motions are almost always en-ergetically favorable compared to continuous motions, with greater energy savings with decreasing duty cycle. We knowfrom section 3 that T ∼ T ∆ ( ∆ is the duty cycle), where T is the mean thrust, and the symbol ∼ means “proportionalto.” Similarly, P ∼ P ∆ . Furthermore, figure 7 shows that U ∼ √ T for a wide range of parameters. The energy ratioshould then behave according to φ = PP U U ∼ ∆ √ ∆ = √ ∆ . Figure 8 indicates that the data follow this trend, which isnot unexpected: a swimmer in this scenario does not expendany energy during the inactive portion of intermittent mo-tions, but still coasts forward, so that intermittent motionsare always energetically favorable.
Figure 8: Ratio of energy expended by intermittent mo-tions to energy expended by continuous motions as afunction of duty cycle for θ = ◦ , all frequencies, m ∗ =
1, and A ∗ w = Apart from energy spent on swimming, aquatic animalsalso expend energy on metabolic processes. To capture this,we consider the metabolic energy ratio ψ = E + E m E + E m , where E m is the energy spent due to metabolic processes.We will assume that the mean power spent on metabolic pro-cesses is the same for continuous and intermittent motion,and that it is a constant fraction c m of the power lost in con-tinuous swimming, P . Hence, ψ = E + E m E + E m , = Pt + c m P tP t + c m P t = ( P + c m P ) d / U (1 + c m ) P d / U = P + c m P (1 + c m ) P U U , since d = d . Values of ψ < c m from 0 to 2, typical in biology (Di Santo and Lauder, private orces and energetics of intermittent swimming 7 communication). With c m = c m (approximately bounded by c m > ψ = P + c m P (1 + c m ) P U U ∼ ( ∆ + c m ) P (1 + c m ) P √ ∆ = ∆ + c m (1 + c m ) √ ∆ . From figure 9 we see that the data and the scaling tell thesame story: even though intermittent motions expend lessenergy on swimming, the metabolic losses can play a sig-nificant role because intermittent motions take more time totraverse a given distance than continuous motions. The extratime to travel will increase the metabolic energy losses, andthis e ff ect may dominate the benefits gained in swimmingenergy losses. Figure 9: Ratio of energy expended by intermittent mo-tions to energy expended by continuous motions, in-cluding metabolic energy losses, as a function of dutycycle for θ = ◦ , all frequencies, m ∗ =
1, and A ∗ w =
10. Each point is an average over all frequencies.The color denotes the value of the metabolic powerfraction c m , 0 to 2 in intervals of 0.25 (dark to light). Another comparison we can make is to consider the en-ergy ratio φ , but restrict ourselves to motions which producethe same mean speed, that is, φ | U mean . For example, it may bethat a continuous motion produces the same mean speed asan intermittent motion with a duty cycle of 0.5 actuating attwice the frequency, but which motion expends more energy?The answer to this question will reveal which gait is bestif a swimmer wants to traverse a given distance in a givenamount of time. The results are plotted in figure 10. Inter-estingly, it appears that intermittent motions continue to be energetically favorable with the added time restriction. En-ergetically optimal duty cycles exist, and savings are greaterfor lower speeds, at least for the data considered here. De-spite having to increase the frequency of actuation in orderto match the mean speed of continuous motions, intermit-tent motions are nevertheless energetically favorable in thiscontext. Figure 10: Ratio of energy expended by intermittentmotions to energy expended by continuous motions, re-stricted to equal mean speeds, as a function of duty cy-cle for θ = ◦ . Dashed lines correspond to equation 7,with m ∗ = A ∗ w =
10. The symbol grey scale cor-responds to three values of (dimensional) mean speedchosen, U mean = . , . , . We can understand this behavior by using the scaling re-lations obtained here and in previous work [12]. The thrustcoe ffi cient approximately follows C T ∆ ∼ c S t − C D , u , where C D , u is the o ff set of the thrust curve. Similarly, thepower coe ffi cient approximately follows C P ∆ ∼ c S t . Note that this relationship is di ff erent than given in [12], butit fits our data almost as well and is more convenient for ourpurposes here. Solving for the power gives P ∼ ∆ · c f A ρ sc = ∆ · c T ∆ · c ρ sc + C D , u U c ! / ρ sc , As shown earlier, the mean speed can be estimated well byequating the mean thrust of the propulsor to the mean dragof the body, so that T = D = ρ U A w C D . D. Floryan, T. Van Buren, and A. J. Smits
The power then becomes P ∼ ∆ c U c / ∆ A w sc C D + C D , u ! / ρ sc . In this case, the energy ratio between intermittent and con-tinuous motions is φ | U mean = EE = PtP t = PP , since t = t . Hence, φ | U mean ∼ ∆ · (cid:16) ∆ A ∗ w C D + C D , u (cid:17) / (cid:0) A ∗ w C D + C D , u (cid:1) / . (7)For the cases analyzed here, we chose A ∗ w =
10 and C D = .
01. This expression is plotted in figure 10 for C D , u = { . , . , . } (for our data, C D , u = . C D , u . The ex-pression is in qualitative agreement with the data, having aU-shape and indicating that there may be a particular dutycycle which is optimal. The forces and energetics of intermittent swimming motions,characterized by the duty cycle, were analyzed and com-pared to continuous swimming. Water tunnel experimentson a nominally two-dimensional rigid foil pitching about itsleading edge showed that mean thrust and power increasedwith increasing duty cycle, all the way up to continuous mo-tion. The mean thrust and power data for di ff erent duty cy-cles were collapsed by dividing the mean thrust and power bythe duty cycle, indicating that thrust and power in one burstcycle are not a ff ected by the previous burst cycle. PIV mea-surements of the wake showed that the dominant structuresshed into the wake were always counter-rotating pairs of vor-tices, regardless of the duty cycle. The PIV measurementscorroborated the assertion that individual cycles of activityare una ff ected by previous cycles.Free swimming speed and energetics were analyzed bynumerically integrating the measured experimental data. Al-though the experimental data was acquired for a foil movingat a constant speed, previous work showed that the thrustproduced by a pitching foil is independent of speed so weexpect that the thrust measured in stationary experimentswould be the same as what would be measured in free swim-ming experiments. For a large range of area and mass ra-tios, the mean speed was then found to be the same as whatwould be calculated by assuming a constant speed and equat-ing mean thrust with mean drag. The mean speed was lowerthan this steady speed only for parameters unlikely to be en-countered in nature, indicating that the constant steady speedis a good approximation to the mean speed.The energetics of intermittent motions were then com-pared to continuous motions according to three criteria: (i) energy expended in traversing a given distance; (ii) energyexpended in traversing a given distance, including metabolicenergy; and (iii) energy expended in traversing a given dis-tance in a given time. Intermittent motions were generallyenergetically favorable as no energy is spent during the in-active portion of the motion where the swimmer still coastsforward. When metabolic energy losses were added, theycould be high enough to make continuous swimming ener-getically advantageous.The assertion that forces are una ff ected by speed forpitching motions, shown in a previous study, was instrumen-tal in being able to use the stationary experiments in the freeswimming analysis. Although the forces produced by gen-eral motions (for example, combinations of pitch and heave)will depend on the speed, this work nonetheless highlightsthe potential benefits and pitfalls of intermittent motions as aswimming protocol.This work was supported by ONR Grant N00014-14-1-0533 (Program Manager Robert Brizzolara). We would alsolike to thank Dr. Keith Moored for stimulating our interestsin intermittent swimming. References [1] A. C. Gleiss, S. J. Jorgensen, N. Liebsch, J. E. Sala,B. Norman, G. C. Hays, F. Quintana, E. Grundy,C. Campagna, A. W. Trites, B. A. Block, and R. P.Wilson. Convergent evolution in locomotory patternsof flying and swimming animals.
Nature Communica-tions , 2:352, 2011.[2] F. E. Fish, J. F. Fegely, and C. J. Xanthopoulos. Burst-and-coast swimming in schooling fish (notemigonuscrysoleucas) with implications for energy economy.
Comparative Biochemistry and Physiology Part A:Physiology , 100(3):633–637, 1991.[3] D. L. Kramer and R. L. McLaughlin. The behavioralecology of intermittent locomotion 1.
American Zool-ogist , 41(2):137–153, 2001.[4] D. Weihs. Energetic advantages of burst swimming offish.
Journal of Theoretical Biology , 48(1):215–229,1974.[5] J. J. Videler and D. Weihs. Energetic advantages ofburst-and-coast swimming of fish at high speeds.
Jour-nal of Experimental Biology , 97(1):169–178, 1982.[6] M. J. Lighthill. Large-amplitude elongated-body the-ory of fish locomotion.
Proceedings of the Royal Soci-ety of London B: Biological Sciences , 179(1055):125–138, 1971.[7] P. W. Webb. Hydrodynamics and energetics of fishpropulsion. 1975. orces and energetics of intermittent swimming 9 [8] E. J. Anderson, W. R. Mcgillis, and M. A. Grosen-baugh. The boundary layer of swimming fish.
Journalof Experimental Biology , 204(1):81–102, 2001.[9] R. W. Blake. Functional design and burst-and-coastswimming in fishes.
Canadian Journal of Zoology ,61(11):2491–2494, 1983.[10] M. H. Chung. On burst-and-coast swimming per-formance in fish-like locomotion.
Bioinspiration & Biomimetics , 4(3):036001, 2009.[11] G. Wu, Y. Yang, and L. Zeng. Kinematics, hydro-dynamics and energetic advantages of burst-and-coastswimming of koi carps (cyprinus carpio koi).
Journalof Experimental Biology , 210(12):2181–2191, 2007.[12] D. Floryan, T. Van Buren, C. W. Rowley, and A. J.Smits. Scaling the propulsive performance of heavingand pitching foils.
Journal of Fluid Mechanics , underreview, 2017.[13] A. Sciacchitano, B. Wieneke, and F. Scarano. PIV un-certainty quantification by image matching.
Measure-ment Science and Technology , 24(4):045302, 2013. [14] E. Akoz and K. W. Moored. Unsteady propulsion by anintermittent swimming gait.
ArXiv:1703.06185 , 2017.[15] T. Van Buren, D. Floryan, D. Quinn, and A. J. Smits.Non-sinusoidal gaits for unsteady propulsion.
PhysicalReview Fluids , under review, 2017.[16] F. M. White. Fluid mechanics, 2011.[17] J. R. Dormand and P. J. Prince. A family of embeddedRunge-Kutta formulae.
Journal of Computational andApplied Mathematics , 6(1):19–26, 1980.[18] L. F. Shampine and M. W. Reichelt. The MATLABODE suite.
SIAM Journal on Scientific Computing ,18(1):1–22, 1997.[19] D. S. Barrett, M. S. Triantafyllou, D. K. P. Yue, M. A.Grosenbaugh, and M. J. Wolfgang. Drag reductionin fish-like locomotion.