TThe optimal swimming sheet
Thomas D. Montenegro-Johnson and Eric Lauga ∗ Department of Applied Mathematics and Theoretical Physics, Centre for Mathematical Sciences,University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK (Dated: July 23, 2018)Propulsion at microscopic scales is often achieved through propagating traveling waves along hair-like organelles called flagella. Taylor’s two-dimensional swimming sheet model is frequently usedto provide insight into problems of flagellar propulsion. We derive numerically the large-amplitudewaveform of the two-dimensional swimming sheet that yields optimum hydrodynamic efficiency; theratio of the squared swimming speed to the rate-of-working of the sheet against the fluid. Usingthe boundary element method, we show the optimal waveform is a front-back symmetric regularizedcusp that is 25% more efficient than the optimal sine-wave. This optimal two-dimensional shapeis smooth, qualitatively different from the kinked form of Lighthill’s optimal three-dimensionalflagellum, not predicted by small-amplitude theory, and different from the smooth circular-arc-likeshape of active elastic filaments.
I. INTRODUCTION
The microscopic world is teeming with organisms andcells that must self-propel through their fluid environ-ment in order to survive or carry out their functions [1].At these very small scales, viscous forces dominate inertiain fluid flows, and a common method of overcoming thechallenges of viscous propulsion through the fluid envi-ronment is by propagating waves along slender, hair-likeorganelles called flagella or cilia [2, 3]. To examine thefluid mechanical basis for microscopic propulsion, Taylor[4] considered a simplified flagellum model comprising atwo-dimensional sheet exhibiting small amplitude travel-ing waves. This seminal work subsequently sparked thedevelopment of other techniques for examining Newto-nian viscous flows such as slender-body theory [5–7] andresistive-force theory [2, 8], as well as other models fornon-Newtonian swimming based on distribution of forcesingularities [9, 10] .Due to its analytical tractability and agreement withmore involved approaches in the small-amplitude limit,Taylor’s swimming sheet has been used to give insightinto many fundamental problems in microscale propul-sion, such as hydrodynamic synchronization betweenwaving flagella [4, 11–13], swimming in non-Newtonianfluids [14–16] and swimming past deformable membranes[10]. These approaches are typically characterized byasymptotic expansion of the flagellar waveform under thecondition that the amplitude of the waves is small whencompared to the wavelength. Recently, Taylor’s small-amplitude expansion was formally extended to arbitrar-ily high order for a pure sine-wave, a method able toproduce results comparable to full numerical simulationsof large amplitude sine-waves with the boundary elementmethod [17].Motivated by the role of evolutionary pressures onthe shape and kinematics of swimming microorganisms, ∗ Electronic address: [email protected] it is relevant to investigate which flagellar waveformis the most energetically efficient for the cell. Foran infinite flagellum, Lighthill [2] showed that in thelocal drag approximation of resistive-force theory, thehydrodynamically-optimal flagellar waveform has a con-stant tangent angle to the swimming direction. Thisleads to the shape of a smooth helix in three dimen-sions, and a singular triangle wave in two dimensions.Whilst the helical waveform is commonly observed inbacterial flagella [18, 19], unsurprisingly the kinked pla-nar waveform is not. Spagnolie and Lauga [20] showedthat this shape singularity in Lighthill’s flagellum can beregularized by penalizing the swimming efficiency by theelastic energy required to bend a flagellum, which mightprovide one explanation for its absence in nature. Thismodel was then improved upon by Lauga and Eloy [21]by proposing an energetic measure based on the internalmolecular cost necessary to deform the active flagellum.For finite-length flagella, Pironneau and Katz [22] showedthat traveling waves are fundamental to optimal propul-sion. Using resistive force theory, they then analyzed op-timal patterns for model spermatozoon exhibiting smallamplitude planar sinusoidal waves and finite amplitudetriangle waves. The optimal stroke pattern of Purcell’sfinite three-link swimmer was found by Tam and Hosoi[23], who then went on to consider optimal gaits for thegreen alga
Chlamydomonas [24].While all past optimization work has focused on three-dimensional slender filaments, the optimal waveform ofTaylor’s two-dimensional swimming sheet has yet to beconsidered. Although the two-dimensionality of the prob-lem makes it less realistic as a model for swimming cells,the fluid dynamics around a sheet can be computed veryaccurately, allowing us to bypass various hydrodynamicmodeling approximations employed in three dimensions.In the present study, we use the boundary elementmethod to examine a sheet propagating large amplitudewaves of arbitrary shape. We derive computationally thewaveform leading to swimming with maximum hydrody-namic efficiency. We show that the optimal waveformfor the swimming sheet is a regularized cusp wave not a r X i v : . [ phy s i c s . f l u - dyn ] J un Wave propagationSwimming direction ⇡ ⇡/ ⇡/ ⇡ x y FIG. 1: A schematic of the wave propagation of a swimmingsheet, here illustrated by a single sine-wave mode, showing thecomputational domain ( − π to π along the x axis), directionof wave propagation (in the swimming frame) and directionof swimming. predicted by small-amplitude analysis. The optimal isqualitatively different from three-dimensional swimmers,both the kinked triangle of Lighthill’s hydrodynamically-optimal flagellum [2] and the circular arcs of internally-optimal active filaments [21] and indicates a qualitativedifference between two- and three-dimensional swimmingat large amplitude. II. FORMULATION OF THE PROBLEM
Newtonian fluid mechanics at microscopic scales is gov-erned by the incompressible Stokes flow equations, µ ∇ u = ∇ p, ∇ · u = 0 , (1)where u is the fluid velocity field, p the dynamic pres-sure and µ is the dynamic viscosity, hereafter non-dimensionalized to µ = 1.We consider the waving sheet model illustrated inFig. 1, and assume the waveform to be fixed and to travelalong the positive x direction at unit speed. The infinitesheet is periodic over the interval [0 , π ], and swimmingis expected to occur in the − x direction, opposite to thedirection of propagation of the wave [4]. Since the sheetis infinite, there is no extrinsic length-scale to the prob-lem, and as such we hereafter non-dimensionalise lengthsusing the reciprocal of the wavenumber, k − . In order fornet swimming to take place with no rotation, we requirethe wave to be odd about the axis x = 0, i.e. ask that y ([0 , − π ] , t ) = − y ([0 , π ] , t ). Without loss of generality,the waveform may therefore be described as a Fourier-sine series where the shape in the swimming frame isdescribed by y ( x, t ) with y ( x, t ) = P (cid:88) p =1 B p sin[ p ( x − t )] . (2)Even modes for the shape, B q with q any integer, are al-ways obtained by our optimization algorithm to be zero,indicating that optimal waveforms are front-back sym-metric waves. The physical reason underlying this front-back symmetry is unclear. Due to kinematic re-versibility, if the shape was asymmetric then an equallyoptimal waveform would be its front-back mirror image,and thus the optimization procedure would always leadto two symmetric solutions. This is not the case and aunique, front-back symmetric shape is always obtained.We thus consider a general waveform represented by y ( x, t ) = N (cid:88) n =1 B n sin[(2 n − x − t )] , (3)and use our computational approach to derive the opti-mal series of coefficients { B n } N , n ≤ N for increasing N.The lack of an extrinsic length-scale to the infinite sheetmeans that our choice of first mode is in some sense ar-bitrary, and thus we will not consider solutions for which | B | >
0, in order to define a fundamental period to thewave.To derive the hydrodynamically-optimal waveform weuse the standard definition of swimming efficiency intro-duced by Lighthill [2]. We therefore compare a useful rateof swimming, ∼ U , to the rate of working of the sheetagainst the fluid, W = − (cid:82) S u · σ · n d s , where the fluidstress is σ = − p I + ( ∇ u + ( ∇ u ) T ) / S is the surfaceof the swimmer over one wavelength, and n is the unitnormal to the sheet into the fluid. We thus seek the setof coefficients { B n } N that maximize the hydrodynamicefficiency, E , defined as E = U − (cid:82) S u · σ · n d s , (4)and numerically compute the value of the swimmingspeed, U , and surface stress, σ · n .In order to impose velocity conditions on the surfaceof the sheet, we solve the problem in a frame of referencethat moves with the propagating wave. Since the sheet isstationary in this frame, the velocity of material elementsis purely tangential [2]. By subtracting the normalizedwave speed, this allows us to retrieve the boundary con-ditions everywhere along the sheet as u ( x ) = − Q cos θ ( x ) + 1 , v = − Q sin θ ( x ) , (5)where Q denotes the ratio between the arclength of thewaveform in one wavelength to the wavelength measuredalong the x direction and θ ( x ) is the tangent angle of thesheet measured about the x axis. Both the value of Q andthe distribution of θ are functions of the wave geometryonly, and thus of the coefficients { B n } N . III. COMPUTATIONAL APPROACH
In order to compute the flow field generated by thesheet, and the resultant surface stresses, we employ theboundary element method [25] with two-dimensional, pe-riodic Green’s functions as in Pozrikidis [26] and Sauzadeet al. [17]. At any point x along the sheet, the velocityat that point is given by the surface integral u j ( x ) = 12 π (cid:90) S [ S ij ( x − x (cid:48) ) f i ( x (cid:48) ) − T ijk ( x − x (cid:48) ) u i ( x (cid:48) ) n k ( x (cid:48) )] d s ( x (cid:48) ) , (6)where n ( x (cid:48) ) is the unit normal pointing into the fluidat x (cid:48) and f i = σ ij n j . Using the notation r = | ˆ x | , thestokeslet tensor, S ij , and stresslet, T ijk , are given by S ij (ˆ x ) = δ ij ln r − ˆ x i ˆ x j r , T ijk (ˆ x ) = 4 ˆ x i ˆ x j ˆ x k r , (7)and represent the solution to Stokes flow due to a pointforce in two dimensions and the corresponding stress re-spectively. Since we are modeling an infinite, 2 π -periodicsheet, we have velocity contributions at x from an infinitesum of stokeslets and stresslets, S p = ∞ (cid:88) n = −∞ I ln r n − ˆ x n ˆ x n r n , T p = ∞ (cid:88) n = −∞ x n ˆ x n ˆ x n r n , (8)where ˆ x n = ( x − x (cid:48) + 2 πn, y − y (cid:48) ) for singularities posi-tioned at x (cid:48) . These infinite sums may be convenientlyexpressed in a closed form as S pxx = A + ˆ y∂ ˆ y A − ,S pyy = A − ˆ y∂ ˆ y A,S pxy = − ˆ y∂ ˆ x A = S pyx ,T pxxx = 2 ∂ ˆ x (2 A + ˆ y∂ ˆ y A ) , T pxxy = 2 ∂ ˆ y (ˆ y∂ ˆ y A ) ,T pxyy = − y∂ ˆ x ˆ y A,T pyyy = 2( ∂ ˆ y A − ˆ y∂ ˆ y ˆ y A ) ,T pijk = T pkij = T pjki , (9)where A = ln [2 cosh(ˆ y ) − x )]. Note that theseare equivalent to, but differ by a minus sign from, thosefound in Sauzade et al. [17], due to our adoption of thesign convention from Pozrikidis [26].For the computational procedure, the sheet is dis-cretized into 500 straight line segments of constant forceper unit length, i.e. the components f , in equation (6)are constant over each straight line element. This dis-cretization breaks equation (6) into a sum of line integralsof singularities, multiplied by the unknown force per unitlength. Numerical evaluation of each non-singular line in-tegral is performed with four-point Gaussian quadrature,whilst singular integrals are treated analytically. This nu-merical discretization produced a < .
02% relative differ-ence in the calculated swimming velocity and efficiencyfor both kinked and unkinked example sheets when com-pared to simulations with 600 and 800 elements with 8-and 16-point Gaussian quadrature, whilst still allowingcalculation in a reasonable time. By decoupling the nu-merical quadrature from the force discretization, com-parable accuracy is achieved for relatively smaller linearsystems [27]. The computational mesh is refined locallyaround regions of high curvature appearing as a result ofthe optimization in order to resolve potential kinks andsingular shapes. Using these parameters, we find conver-gence of the waveform and optimal efficiency for N = 30.In order to reach this number of coefficients quickly, it is important to provide a good starting guess for the coeffi-cients B n at the beginning of the optimization procedure.Since we are interested in the convergence of our wave-form with an increased number of points, we numericallyoptimize for n = 1 , . . . ,
30 sequentially, and use the con-verged optimal coefficients for the n = N − n = N waveform, with B N initially set to zero. In order to ensure that our solutionrepresents the global maximum of the hydrodynamic effi-ciency (4), multiple initial conditions were tried, all foundto be optimized to the same waveform. Optimization iscarried out using the standard fminunc function in Matlab using the quasi-newton algorithm.
IV. THE SHAPE OF THE OPTIMALSWIMMING SHEET
With this framework established, we are able to com-pute the shape of the optimal swimming sheet. It isfirst instructive to ask what would be predicted from theanalytical small-amplitude approach. In that case, thewaveform is written as y ( x, t ) = (cid:15) ∞ (cid:88) n =1 B n e in ( x − t ) , (10)where (cid:15) is a small dimensionless amplitude. To secondorder in the amplitude [4], the swimming speed, U , andwork done against the fluid, W , are given by U ∼ ∞ (cid:88) n =1 n | B n | , W ∼ ∞ (cid:88) n =1 n | B n | . (11)Because the work done is proportional to n , whereasthe swimming velocity is proportional to n , higher-ordermodes are energetically penalized compared to lowermodes. Therefore for a fixed amount of mechanical powerexpended by the swimmer, it is more efficient to dis-tribute all of that power to the first mode n = 1; un-der the small amplitude approximation, the most efficientwaveform is therefore a single sine-wave of period 2 π .Relaxing the small-amplitude constraint, we show inFig. 2a the optimal waveform obtained numerically for N = 30 odd Fourier modes. The set of coefficients B n that describe this waveform are given in Appendix A anda three-dimensional sheet propagating this wave is fur-ther shown in Fig. 2b. The optimal swimming sheet ap-pears to take the shape of a regularized cusp wave, qual-itatively different from the single-mode sine-wave pre-dicted by the asymptotic analysis. We further plot inFig. 2c the distribution of slopes along the sheet, show-ing that the wave has an almost straight section (con-stant slope), steepening towards smooth wave crest. Theangle of the slope at the point of symmetry x = 0 is ap-proximately 36 . ◦ , close to the optimum value of 40 . ◦ obtained for Lighthill’s three-dimensional flagellum viaresistive-force theory with a drag anisotropy ratio of 1 / Optimal waveform x y Optimal sheetSlope x d y / d x Curvature x c u r v a t u r e Hydrodynamic optima x y Large-amplitude sheetSmall amplitude sheetInfinite flagellum
FIG. 2: The optimal swimming sheet. (a) Optimal waveform obtained computationally, which exhibits an almost straight regionfollowed by a steepening of the wave crest into a rounded cap; (b) Three-dimensional visualization of the optimal swimmingsheet; Distribution of slope (c) and curvature (d) of the optimal waveform showing largely straight regions; (e) Waveform of theoptimal sheet shown for comparison with the pure sine-wave predicted by small amplitude theory [4] and Lighthill’s optimaltriangle wave for a flagellum with drag anisotropy ratio of 1 / [2]. We display in Fig. 2d the distribution of curvaturealong the sheet; while the curvature at the wave crest in-creases, it remains finite. For comparison, our predictedoptimal waveform is plotted against the pure sine-waveof small amplitude theory [4] and the triangle wave pre-dicted by Lighthill [2] in Fig. 2e.Since our optimization procedure finds the optimal so-lution for incremental values of the number of coefficients, N , used to describe the wave, we can investigate conver-gence of all optimal waveforms described by N rangingfrom 1 to 30. The convergence for the swimming effi-ciency is shown in Fig. 3a while the dependence of themaximum curvature on N is plotted in Fig. 3b. For N = 30, the optimal waveform is over 25% more effi-cient than the optimal one-mode sine-wave ( N = 1). Theswimming efficiency appears to reach its asymptote near N = 13, which corresponds to the peak in the maximumcurvature, but thereafter continues to increase slightlybefore reaching its converged value of E ≈ . N ≥ N = 13, it appears that subsequent modes serve tosteepen the waveform as it approaches around the crest.Such steepening is likely hydrodynamically favorable intwo dimensions since fluid cannot pass around the sheetas it would around three-dimensional flagellum. How-ever, steepening results in a region of high curvature atthe wave crest, which induces locally high viscous dissipa-tion in the fluid, and so there appears to be an efficiencytrade-off between wave steepening and minimizing curva-ture. For N ≥
14, the wavelength of the Fourier modesis on the order of the length of the cap on the wave crest.These modes are then able to decrease the maximum cur-vature without decreasing the slope of the wave, yieldingsmall increases in efficiency until the curvature convergesfor N ≥
30. We further display the convergence of theoptimal waveform as a function of the number of coeffi-cients, N , in Fig. 4. Despite the decrease in maximum E cicency0 10 20 300 . . .
11 Number of coe cients E c i e n c y E Curvature0 10 20 300102030 Number of coe cients M a x i m u m c u r v a t u r e FIG. 3: Convergence of the swimming efficiency (a) and max-imum wave curvature (b) as a function of the number of oddFourier modes in the optimisation, N . curvature seen in Fig. 3b for N ≥
14, all waveforms be-tween N = 10 and 30 are virtually indistinguishable byeye.The trade-off between wave steepening and reducingcurvature can be further investigated by examining afamily of waves of the form B n = C ( − n − (2 n − m , n ≥ , (12)where C is a constant. The value of m dictates the de-cay of the Fourier coefficients with the Fourier mode, andwith the choice m = 2 , Eq. (12) leads to the triangularwaveform of Lighthill’s optimal flagellum. The choice ofalternating sign is informed both by the series for thetriangle wave, and by the coefficients of our optimal so-lution (up to N = 14). Cusps are obtained for m < m > N , we retrieve an approximateregularized cusp wave. Fig. 5 shows iso-contours of the . . . x y
30 Coe↵s10 Coe↵s5 Coe↵s3 Coe↵s2 Coe↵s1 Coe↵s
FIG. 4: Convergence of the shape of the optimal waveformfor increasing numbers of odd coefficients. Shapes are shiftedby 0 . π along the y direction for visualization. ⇡ ⇡ x y FIG. 5: Efficiences of optimal truncated ‘cusp-like’ waveforms(12) as a function of number of odd modes in the description N and decay rate m , showing a maximum at N = 9 , m ≈ . efficiency of waves described by Eq. (12) for the optimalvalue of the amplitude, C , as a function of the numberof coefficients used to describe the wave, N , and the de-cay rate, m . The optimal efficiency E = 0 . N = 9, for C = 1 .
237 and m = 1 . N ≥
50) isused to describe the curve, the kink at the wave crestis sufficiently resolved as to no longer be regularized. Inthis case, the optimal jumps to an unkinked profile withwhich C = 0 . m = 2 . . V. DISCUSSION
Taylor’s swimming sheet model is commonly used toaddress a range of phenomena in the biological physicsof small-scale locomotion. A natural question to raiseis the relevance of a two-dimensional geometry to thethree-dimensional locomotion of flagellated cells. Inthis paper, we used the boundary element method tocompute the swimming efficiency of arbitrary wave-forms in two dimensions. By focusing on the ques-tion of optimal waveform for locomotion, we show thatthe optimal two-dimensional waveform is a regularizedcusp, which is about 25% more efficient than a sim-ple sine-wave. This result is different from the three-dimensional hydrodynamically-optimal triangle wave de-rived by Lighthill [2]; the slope of the straight section isshallower, the waveform steepens towards the wave crestand there is no discontinuity in the slope but rather aregularized cusp. The result is also different to the three-dimensional internally-optimal wave, which is composedof circular arcs joined by straight lines [21]. Although it isknow that the dynamics of a swimming sheet can providequalitative insight into the hydrodynamics of small-scale locomotion, differences with three-dimensional results ex-ist therefore at large amplitude.
Acknowledgements
The authors would like to thank Gwynn Elfring foruseful discussions. This work was funded in part by theEuropean Union through a Marie Curie CIG to EL.
Appendix A: Coefficients of the optimal waveform
The Fourier coefficients for the optimal waveform for N = 30 are given by × − : B = 114 . B = − . B = 12 . B = − . B = 4 . B = − . B = 2 . B = − . B = 0 . B = − . B = 0 . B = − . B = 0 . B = − . B = − . B = 0 . B = − . B = 0 . B = − . B = 0 . B = − . B = 0 . B = − . B = 0 . B = − . B = 0 . B = − . B = 0 . B = − . B = 0 . [1] D. Bray, Cell Movements (Garland Publishing, NewYork, NY, 2000).[2] J. Lighthill,
Mathematical Biofluiddynamics (SIAM,1975).[3] S. Childress,
Mechanics of Swimming and Flying (Cam-bridge Univ. Press, 1981).[4] G. I. Taylor, Proc. Roy. Soc. Lond. A , 447 (1951).[5] G. J. Hancock, Proc. R. Soc. Lond. A , 96 (1953).[6] J. B. Keller and S. I. Rubinow, J. Fluid Mech. , 705(1976).[7] R. E. Johnson, J. Fluid Mech. , 411 (1980).[8] J. Gray and G. J. Hancock, J. Exp. Biol. , 802 (1955).[9] T. D. Montenegro-Johnson, A. A. Smith, D. J. Smith,D. Loghin, and J. R. Blake, Euro. Phys. J. E , 1 (2012).[10] J. C. Chrispell, L. J. Fauci, and M. Shelley, Phys. Fluids , 013103 (2013).[11] G. J. Elfring and E. Lauga, Phys. Rev. Lett. , 088101(2009).[12] G. J. Elfring, O. S. Pak, and E. Lauga, J. Fluid Mech. , 505 (2010).[13] G. J. Elfring and E. Lauga, J. Fluid Mech. , 163(2011).[14] E. Lauga, Phys. Fluids , 083104 (2007).[15] J. Teran, L. Fauci, and M. Shelley, Phys. Rev. Lett. , 038101 (2010).[16] T. D. Montenegro-Johnson, D. J. Smith, and D. Loghin,Phys. Fluids , 081903 (2013).[17] M. Sauzade, G. J. Elfring, and E. Lauga, Physica D ,1567 (2011).[18] C. Brennen and H. Winet, Ann. Rev. Fluid Mech. , 339(1977).[19] S. E. Spagnolie and E. Lauga, Phys. Rev. Lett. ,058103 (2011).[20] S. E. Spagnolie and E. Lauga, Phys. Fluids , 031901(2010).[21] E. Lauga and C. Eloy, J. Fluid Mech. , R1 (2013).[22] O. Pironneau and D. F. Katz, J. Fluid Mech. , 391(1974).[23] D. Tam and A. E. Hosoi, Phys. Rev. Lett. , 068105(2007).[24] D. Tam and A. E. Hosoi, Proc. Nat. Acad. Sci. USA(2011).[25] G. K. Youngren and A. Acrivos, J. Fluid Mech. , 377(1975).[26] C. Pozrikidis, J. Fluid Mech. , 495 (1987).[27] D. J. Smith, Proc. Roy. Soc. A465