Inertialess Swimming and Propulsion of Slender BodiesbyZhiwei PengB.Sc., Beijing University of Aeronautics & Astronautics, 2014A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Mechanical Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)April 2016c© Zhiwei Peng 2016AbstractIn this thesis, two problems relevant to the biological locomotion in inertialess environmentsare studied, one is the characteristics of undulatory locomotion in granular media, the other isthe optimal flexibility of a driven microfilament in a viscous fluid.Undulatory locomotion is ubiquitous in nature and observed in different media, from theswimming of flagellated microorganisms in biological fluids, to the slithering of snakes on land,or the locomotion of sandfish lizards in sand. Despite the similarity in the undulating pattern,the swimming characteristics depend on the rheological properties of different media. Analysisof locomotion in granular materials is relatively less developed compared with fluids partiallydue to a lack of validated force models but recently a resistive force theory in granular media hasbeen proposed and shown useful in studying the locomotion of a sand-swimming lizard. In thiswork, we employ the proposed model to investigate the swimming characteristics of a slenderfilament, of both finite and infinite length, undulating in a granular medium and compare theresults with swimming in viscous fluids. In particular, we characterize the effects of driftingand pitching in terms of propulsion speed and efficiency for a finite sinusoidal swimmer. Wealso find that, similar to Lighthill’s results using resistive force theory in viscous fluids, thesawtooth swimmer is the optimal waveform for propulsion speed at a given power consumptionin granular media.Though it is understood that flexibility can improve the propulsive performance of a fila-ment in a viscous fluid, the flexibility distribution that generates optimal propulsion remainslargely unexplored. In this work, we employ the resistive force theory combined with theEuler-Bernoulli beam model to examine the optimal flexibility of a boundary driven filamentin the small oscillation amplitude limit. We show that the optimality qualitatively dependson the boundary actuation. For large amplitude actuation, our numerics show that complexasymmetry in the waveforms emerge. The results complement our understanding of inertia-less locomotion and provide insights into the effective design of locomotive systems in variousenvironments.iiPrefaceThe research presented in this thesis is original work of the thesis author under the supervi-sion of Professor Gwynn J. Elfring in collaboration with Professor On Shun Pak (Santa ClaraUniversity, USA).A version of Chapter 2 was presented at the American Physical Society 68th Annual Di-vision of Fluid Dynamics Meeting (Boston, MA, November 22-24, 2015) by the thesis authorand published in Physics of Fluids. Z. Peng, O.S. Pak, G.J. Elfring (2016) Characteristics ofundulatory locomotion in granular media, Phys. Fluids, 28, 031901.Chapters 3 and 4 are original, unpublished work of the thesis author.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Swimming at low Reynolds number . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Slender body hydrodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.3 Thesis outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 Undulatory locomotion in granular media . . . . . . . . . . . . . . . . . . . . . 42.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.2 Mathematical Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.2.1 Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.2.2 Resistive force theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.2.3 Swimming efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.2.4 Waveforms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.3 Bodies of infinite length . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.3.1 Optimal shape: numerical results . . . . . . . . . . . . . . . . . . . . . . 112.3.2 Sawtooth and sinusoid . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.4 Bodies of finite length . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.4.1 Geometries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.4.2 Pitching, drifting and reorientation . . . . . . . . . . . . . . . . . . . . . 172.4.3 Swimming performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22ivTable of Contents3 Optimal flexibility of a driven microfilament . . . . . . . . . . . . . . . . . . . 243.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 243.2 Mathematical formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253.2.1 Enthalpy functional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.2.2 Elastohydrodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 273.2.3 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 283.2.4 Non-dimensionlization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 293.2.5 Propulsive thrust . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 293.3 Two-segment filaments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303.3.1 Elastic-elastic case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 323.3.2 Rigid-elastic case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 343.3.3 Elastic-rigid case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 353.4 Torsional spring . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 373.5 Continuous optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 383.6 Effect of boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 403.6.1 Torque-free displacement oscillation . . . . . . . . . . . . . . . . . . . . . 403.6.2 Angle oscillation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 413.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 434 Nonlinear dynamics of a driven microfilament . . . . . . . . . . . . . . . . . . 444.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 444.2 Mathematical formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 444.3 Numerical implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 454.4 Uniform stiffness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49AppendicesA Numerical implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54A.1 Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54A.2 Numerical solution for finite swimmers . . . . . . . . . . . . . . . . . . . . . . . 54B Solution of an elastic-elastic filament . . . . . . . . . . . . . . . . . . . . . . . . 56B.1 Equations of motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56B.2 Propulsive force . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57C Boundary conditions for angle oscillation . . . . . . . . . . . . . . . . . . . . . . 59vList of Tables2.1 The parameters of the resistive force model in LP and CP GM as obtained byMaladen et al. [20] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9viList of Figures2.1 Illustration of an undulating slender filament and the resistive force theory ingranular media. The body propagates a prescribed waveform to propel itself.Each element ds experiences a drag force dF = f ds. The basis vectors {ex, ey}and the position vectors of its head x(0, t) and a material point x(s, t) on thebody in the lab frame are shown (ez = ex × ey). The angle between the localvelocity u and unit tangent vector t is ψ(s, t). . . . . . . . . . . . . . . . . . . . . 62.2 Undulating filaments with a single wave (N = 1). (a): sinusoid, kX0 = 0; (b)sawtooth, kX0 = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.3 Optimal shapes in terms of swimming efficiency for an infinite filament in agranular substrate (LP, CP) and Newtonian fluid. The spatial coordinates arescaled to the same wave length. For loosely packed granular material, the optimalshape is almost the same as the analytical result of Lighthill’s in Newtonian fluid. 122.4 (a): Swimming speed of infinite sawtooth waveforms as a function of amplitude (or bending angle β) in granular material and Newtonian fluids. The dashedlines indicate the small and large amplitude asymptotic solutions. (b): Efficiencyof infinite sawtooth waveforms as a function of amplitude (or bending angle β)in granular material and Newtonian fluids. . . . . . . . . . . . . . . . . . . . . . . 152.5 A comparison of the swimming speed (a) and efficiency (b), as a function of waveamplitude for sawtooth and sinusoidal waveforms in granular substrates andNewtonian fluids. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.6 Shapes of swimming finite length single wave (N = 1) sinusoidal filaments fordifferent wave amplitude . (a): the odd sine configuration, with kX0 = 0, (b):the even cosine configuration, with kX0 = pi/2. The waveforms are rescaled tothe same wave length for better comparison. . . . . . . . . . . . . . . . . . . . . . 162.7 Trajectory of the head x(0, t) (black solid lines) and trajectory of the swimmercentroid (dotted lines) for swimming finite sinusoidal filaments with N = 1 and = 1 that possess odd/even symmetry at t = 0 in loosely packed GM. Thefilament swims towards the left when the wave propagates to the right. If theconfiguration possesses even symmetry it does not undergo a net reorientation. . 17viiList of Figures2.8 Parametric plots for the magnitude of the reorientation angle |〈θ〉 − θ0| for asingle wave (N = 1) sinusoid in (a) loosely packed GM, (b) closely packed GMand (c) Newtonian fluids. (d): Plots of |〈θ〉 − θ0| against the wave amplitude for the odd sine configuration in GM and Newtonian fluids. |〈θ〉− θ0| is periodicwith a period of pi. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.9 Maximum instantaneous pitching angle θmp as a function of the wave amplitude for single wave (N = 1) sinusoidal swimmers in GM and Newtonian fluids. . . 192.10 Swimming speed U/V as a function of the dimensionless amplitude for differentnumber of waves N in (a) loosely packed GM and (b) closely packed GM. Thesolid lines denote the swimming speed of an infinite sinusoid. . . . . . . . . . . . 202.11 Swimming efficiency η as a function of the dimensionless amplitude for differentnumber of waves N in (a) loosely packed GM and (b) closely packed GM. Theshaded regions represent the observed values of for lizards reported in theliterature [20, 21]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212.12 (a) Swimming speed as a function of the number of waves in GM. (b) Swim-ming efficiency as a function of the number of waves in GM. The dimensionlessamplitude is fixed ( = 1). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212.13 Maximum instantaneous pitching angle as a function of the number of waves inGM. The dimensionless amplitude is fixed ( = 1). . . . . . . . . . . . . . . . . . 223.1 Propulsive thrust generated by a cantilevered filament of uniform bending stiff-ness under a displacement actuation at one extremity. The two blue dots indicatetwo filaments with different bending stiffnesses generating the same propulsiveforce. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303.2 Shapes of the displace driven cantilevered filament in different times over a pe-riod. From top to bottom, the sperm number increases (top:Sp = 0.5, middle:Sp =1.89, bottom:Sp = 8), namely, the filament becomes more flexible. . . . . . . . . 313.3 Schematic of a serially connected filament. The segment at the actuation end hasa bending flexibility A1 and the segment at the free end has a bending flexibilityA2, β = A2/A1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 323.4 Propulsive thrust generated by two-segment filaments as a function of α at dif-ferent values of stiffness ratio β. For β = 0.144 (A2 < A1, more flexible materialsat the free end), the maximum propulsive force generated is greater than themaximum achievable thrust of a filament with any uniform stiffness. . . . . . . . 333.5 Propulsive thrust generated by a cantilevered rigid-elastic filament under dis-placement oscillation as a function of α for different sperm numbers Sp2, whichare indicated by the numeric values around each line. The horizontal dash-dottedline denotes the maximum propulsion of a filament with uniform stiffness. . . . . 34viiiList of Figures3.6 Propulsive thrust generated by an elastic-rigid filament (A2 =∞) as a functionof α for different sperm numbers of the flexible part. For small sperm num-bers, a monotonic increasing of the propulsive thrust is observed as α increases.For larger sperm numbers, however, the variation of propulsive thrust is non-monotonic. Note that the sperm number of the rigid part would be zero. Thehorizontal dash-dotted line denotes the maximum propulsive force achievable bya filament with uniform stiffness. . . . . . . . . . . . . . . . . . . . . . . . . . . . 363.7 Propulsive thrust generated by a torsional spring arrangement at the actuationend connected to a rigid filament as a function of the spring constant. Theoptimum propulsive force is F(2)p = 0.1875. . . . . . . . . . . . . . . . . . . . . . . 383.8 (a) Optimal linear stiffness distribution, Sp0 ≈ 1.77. (b) Optimal quadraticstiffness distribution denoted by the dashed line (Sp0 ≈ 1.51, F (2)p ≈ 0.2224) . . . 393.9 (a) Propulsive thrust generated by a filament of uniform bending stiffness ac-tuated at one end [4]. The two blue dots indicate two filaments with differentbending stiffnesses generating the same propulsive force. Propulsive thrust gen-erated by two-segment filaments as a function of α at different values of stiffnessratio β. For β = 4.17 (A2 > A1, more flexible materials at the actuation end),the maximum propulsive force generated is greater than the maximum achievablethrust of a filament with any uniform stiffness. . . . . . . . . . . . . . . . . . . . 413.10 (a) Propulsive thrust generated by a filament of uniform bending stiffness withone end under angle oscillation [57]. (b) Propulsive thrust generated by a rigid-elastic filament (A1 = ∞) as a function of α for different sperm numbers (Sp2)of the flexible part. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 424.1 Shape (not to scale) of a deforming filament with uniform stiffness Sp = 4 at = 0.1 at different times npi/4, where n = 1, 2, ..., 8, the intensity of colordecreases as time increases. The numerical solution is calculated with N = 160and ∆t = 5 × 10−4. The red dashed lines denote the shapes obtained from thesmall amplitude linear theory [4]. The shape from numerical simulation matchesvery well with those from the small amplitude asymptotic expansion. . . . . . . . 474.2 Shape of a deforming filament with uniform stiffness Sp = 4 at = 1 at differenttimes npi/4, where n = 0, 1, 2, ..., 8, the intensity of color decreases as time in-creases. The solution is calculated with N = 160 and ∆t = 5×10−4. Pronouncedbuckling is observed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48ixAcknowledgementsFirst of all, I would like to thank my thesis advisor, Gwynn Elfring and my collaboratorand mentor On Shun Pak, for being awesome and inspiring. I feel very lucky to have the chanceto work with them. I thank Gwynn and On Shun for guidance and encouragement throughoutthe journey of my research at UBC.I have been very fortunate to be surrounded by wonderful and talented colleagues, fromwhom I learned a lot both academically and personally. I would like to thank my lab mates,Gaurav, Babak, Charu and Silvia for being brilliant and supportive. Initially, Gaurav and Iwere sharing an office in the Rusty Hut, which indeed is a rusty box. Luckily enough, we movedout to a new lab where all the group members have been creating glorious work in the ElfringAlley.Finally, I would like to thank my friends and family without whose support I could not havemade this far.xDedicationTo my parents, for their love and support.xiChapter 1Introduction1.1 Swimming at low Reynolds numberBiological locomotion in fluids is ubiquitous in nature and spans a wide range of lengthscales from the undulatory locomotion of flagellated bacteria in a microscopic world to theswimming of sharks in the vast ocean. It plays a crucial role in predation, avoiding predators orreproduction throughout the whole lifespan of a swimming animal [1]. Biological locomotion influids have received substantial attention from biologists, engineers and mathematicians alikein recent decades, among which the study of locomotion at small scales relevant to bacteria orflagella remains more recent [2, 3].The physics governing locomotion in fluids at small scales are qualitatively different fromthose of the mesoscale or macroscale locomotion. For swimming bacteria or flagella, inertiaplays a negligible role while viscous forces dominate. In general, the motion of a Newtonianfluid is governed by the Navier-Stokes equation,ρ(∂u∂t+ u · ∇u)= −∇p+ µ∇2u (1.1)where u(x, t) is the velocity of the fluid at x in space and time t, µ is the viscosity of thefluid, ρ is the density of the fluid and p is the pressure. If we non-dimensionalize the Navier-Stokes equations with respect to a characteristic length scale of the microorganism L and acharacteristic velocity scale U0, we have the dimensionless equation given byRe(∂u∂t+ u · ∇u)= −∇p+∇2u (1.2)where the same variables as the dimensional ones are used andRe =ρU0Lµ(1.3)is the Reynolds number which compares the relative importance of inertia effects to viscousforces. For a typical microorganism such as E. Coli (U0 ≈ 30µm/s, L ≈ 10µm) swimmingin water, the Reynolds number is on the order of 3 × 10−4 and thus negligible [4]. In themathematical limit of zero Reynolds number, we have the Stokes equations governing the motion11.2. Slender body hydrodynamicsof a fluid,−∇p+∇2u = 0, ∇ · u = 0. (1.4)The linearity and time independence of the Stokes equations leads to kinematic reversibility,which is well-described by Purcell’s famous “scallop theorem” states that a reciprocal motion (adeformation that exhibits time-reversal symmetry) cannot generate any net propulsive thrust[5]. In order to break the constraint of time-reversibility, many microorganisms including flag-ellated bacteria and spermatozoa achieve self-propulsion by passing deformation waves alongtheir slender flexible bodies [3, 6, 7].1.2 Slender body hydrodynamicsIf the length of a swimming cylindrical body L is much larger than its radius r, r/L 1.In this case, instead of solving the Stokes equations in the fluid domain, we can obtain a localdrag law, which is the so-called Resistive Force Theory (RFT)[6, 8, 9]. RFT states that theviscous force per unit length on the body at a point is related linearly to the local filamentvelocity, namely,fvis = −(ξ⊥nn + ξ‖tt) · xt, (1.5)where the hydrodynamics at this order is characterized by the tangential ξ‖ and normal ξ⊥resistive coefficients, the subscript t indicates differentiation with respect to time and n and tare the local unit normal and tangent vectors along the body. For a slender, cylindrical rod,the resistance ratio γ = ξ⊥/ξ‖ → 2 as L/r → ∞. This drag anisotropy (γ 6= 1) is crucial forlocomotion at low Reynolds number [3, 4].1.3 Thesis outlineThis thesis studies the characteristics and optimality in swimming and propulsion of aslender filament in granular media and viscous fluids.Chapter 2 examines the characteristics of undulatory locomotion of a slender swimmerin granular media and compare the results with those for swimming in a viscous fluid. Wecharacterize the complex kinematics and optimal swimming and discuss the similarities toswimming in viscous fluids.Chapter 3 presents a mathematical modeling of the propulsion of a filament with nonuniformflexibility along the body under a boundary actuation. We explore the optimal flexibilitydistribution that maximizes the propulsive force of a driven filament for small amplitude. Wenote that the optimality qualitatively depends on the boundary actuation, so that one may notextend the results of one case to other cases of propulsion where the actuation might differ.21.3. Thesis outlineChapter 4 investigates the fully nonlinear dynamics of a boundary driven filament using anumerical approach. For small amplitude, we compare the results from numerics to the smallamplitude asymptotic solution and observe good agreement between them. For large amplitude,we show that the linear theory breaks down and complex asymmetry in the waveforms emerge.The results complement our understanding of inertialess locomotion and provide insights intothe effective design of locomotive systems in various environments.3Chapter 2Characteristics of undulatorylocomotion in granular media12.1 IntroductionUndulatory locomotion, the self-propulsion of an organism via the passage of deformationwaves along its body, is ubiquitous in nature [10, 11]. Flagellated microorganisms swim in fluids[5, 8, 9, 12–14], snakes slither on land [15–18]and sandfish lizards (Scincus scincus) undulate ingranular substrates [19–21]. Yet the underlying physics differ: from viscous forces [3] in fluidsto frictional forces [20] in terrestrial media. The investigation of these undulatory mechanismsin different environments advances our understanding of various biological processes [2, 11] andprovides insights into the effective design of biomimetic robots [22, 23].The swimming of microorganisms in Newtonian fluids, where viscous forces dominate in-ertial effects, is governed by the Stokes equations [3]. Despite the linearity of the governingequation, locomotion problems typically introduce geometric nonlinearity, making the problemless tractable [24]. For slender bodies such as flagella and cilia, Gray and Hancock [8] exploitedtheir slenderness to develop a local drag model, called resistive force theory (RFT), which hasbeen shown useful in modeling flagellar locomotion and the design of synthetic micro-swimmers[3, 4]. In this local theory, hydrodynamic interactions between different parts of the body areneglected and the viscous force acting on a part of the body depends only on the local velocityrelative to the fluid. Using RFT, Lighthill showed that, for an undulating filament of infi-nite length, the sawtooth waveform is the optimal beating pattern maximizing hydrodynamicefficiency [9].Locomotion in granular media (GM) is relatively less well understood due to their complexrheological features [25, 26]. The frictional nature of the particles generates a yield stress, athreshold above which the grains flow in response to external forcing [26]. Different from viscousfluids, the resistance experienced by a moving intruder originates from the inhomogeneous andanisotropic response of the granular force chains, which are narrow areas of strained grainssurrounded by the unstrained bulk of medium [27]. At low locomotion speed, where the granularmatter is in a quasi-static regime, the effect of inertia is negligible compared to frictional andgravitational forces from granular media [21], which is similar to that of a low Reynolds-number1A version of Chapter 2 has been published in Physics of Fluids. Peng, Zhiwei, Pak, On Shun and Elfring,Gwynn J., Characteristics of undulatory locomotion in granular media, Phys. Fluids, 28, 031901 (2016)42.2. Mathematical Formulationfluid. In this regime, studies measuring the drag force of an intruder moving through a GMreveal that the drag force is independent of the speed of the intruder, but it increases with thedepth of GM and proportional to the size of the intruder [27–31].Recently, Maladen et al. [20] studied the subsurface locomotion of sandfish in dry granularsubstrates. While the crawling and burying motion of a sandfish is driven by its limbs, anundulatory gait is employed for subsurface locomotion without use of limbs. Using high speedx-ray imaging, the subsurface undulating pattern of the sandfish body was found to be well de-scribed by a sinusoidal waveform. A major challenge in the quantitative analysis of locomotionin granular materials is a lack of validated force models like the Stokes equation in viscous fluids[25, 26]. But inspired by the success of RFT for locomotion in viscous fluids, Maladen et al. [20]developed an empirical RFT in dry granular substrates for slender bodies (Sec. 2.2.2), whichwas shown effective in modeling the undulatory subsurface locomotion of sandfish [20]. Theproposed force model thus enables theoretical studies to address some fundamental questionson locomotion in granular media. In this paper we employ the proposed RFT to investigatethe swimming characteristics of a slender filament of finite and infinite length undulating in agranular medium and compare the results with those in viscous fluids. In particular, previousanalysis using the granular RFT considered only force balance in one direction [20] and hencea swimmer can only follow a straight swimming trajectory in this simplified scenario. Herewe extend the results by considering a full three-dimensional force and torque balances, result-ing in more complex kinematics such as pitching, drifting and reorientation. The swimmingperformance in relation to these complex kinematics is also discussed.This chapter is organized as follows. We formulate the problem and review the recentlyproposed RFT in granular media in Sec. 2.2. Swimmers of infinite length are first considered(Sec. 2.3): we determine that the optimal waveform maximizing swimming efficiency, similar toresults in viscous fluids, is a sawtooth (Sec. 2.3.1); we then study the swimming characteristics ofsawtooth and sinusoidal swimmers in granular media and compare the results with swimming inviscous fluids (Sec. 2.3.2). Next we consider swimmers of finite length (Sec. 2.4) and characterizethe effects of drifting and pitching in terms of propulsion speed and efficiency, before concludingwith remarks in Sec. 2.5.2.2 Mathematical Formulation2.2.1 KinematicsWe consider an inextensible cylindrical filament of length L and radius r such that r L, and assume that it passes a periodic waveform down along the body to propel itself ingranular substrates. Following Spagnolie and Lauga [32], the waveform is defined as X(s) =[X(s), Y (s), 0]T, where s ∈ [0, L] is the arc length from the tip. The periodicity of the waveform52.2. Mathematical Formulationcan then be described asX(s+ Λ) = X(s) + λ, Y (s+ Λ) = Y (s), (2.1)where λ is the wave length and Λ the corresponding arc length along the body. N is the numberof waves passed along the filament. Note that L = NΛ and λ = αΛ, where 0 < α < 1 is due tothe bending of the body [32].✓exeyx(0, t) x(s, t)t'(0, t)sexeyx(s, t)x(0, t)✓'(0, t)stfk = Ck cos f? = C?( ) sin fu? uk ufontsize:8Figure 2.1: Illustration of an undulating slender filament and the resistive force theory ingranular media. The body propagates a prescribed waveform to propel itself. Each element dsexperiences a drag force dF = f ds. The basis vectors {ex, ey} and the position vectors of itshead x(0, t) and a material point x(s, t) on the body in the lab frame are shown (ez = ex×ey).The angle between the local velocity u and unit tangent vector t is ψ(s, t).Initially, the filament is oriented along the x-axis of the lab frame with its head at x0. Attime t, the filament is passing the waveform at a phase velocity V (with constant phase speedV ) along the waveform’s centerline, which is oriented at an angle θ(t) to the x-axis (Fig. 2.1).In a reference frame moving with the wave phase velocity V, a material point on the filament ismoving tangentially along the body with speed c = V/α, and hence the period of the waveformis T = λ/V = Λ/c. By defining the position vector of a material point at location s and time tin the lab frame as x(s, t), we obtainx(s, t)− x(0, t) = Θ(t) ·R(s, t), (2.2)62.2. Mathematical FormulationwhereΘ(t) =cos θ(t) − sin θ(t) 0sin θ(t) cos θ(t) 00 0 1 (2.3)is the rotation matrix, and R(s, t) = X(s, t)−X(0, t), and note that X(s, t) = X(s− ct). Then,the velocity of each material point in the lab frame would beu(s, t) = x˙(0, t) + θ˙Θ ·R⊥ + Θ · R˙, (2.4)where R⊥ = ez ×R, and dot denotes time derivative. The unit tangent vector in the directionof increasing s ist = xs = Θ ·Xs(s, t), (2.5)where the subscript s denotes the derivative with respect to s. The angle between the localvelocity vector u and the local unit tangent vector t is ψ:cosψ = uˆ · t, uˆ = u‖u‖· (2.6)Now, to define the waveform we specify the tangent angle made with the centerline of thewaveformϕ(s, t) = arctanYsXs, (2.7)orXs = [cosϕ, sinϕ, 0]T. (2.8)Note that we have the following geometric relations:R =∫ s0Xs ds, R˙ =∫ s0ϕ˙X⊥s ds, (2.9)t = Θ ·Rs = Θ ·Xs, (2.10)where X⊥s = ez ×Xs, andα =λΛ=1Λ∫ Λ0cosϕds. (2.11)The inextensibility assumption requires that ∂[xs ·xs]/∂t = 0, and the arc-length parameteriza-tion of the swimming filament naturally satisfies this constraint. The tangent angle is specified72.2. Mathematical Formulationas a composition of different Fourier modes2:ϕ(s, t) =n∗∑n=1{an cos[2pinΛ(s− ct)]+ bn sin[2pinΛ(s− ct)]}, (2.12)wherean =2Λ∫ Λ0ϕ(s, 0) cos[2pinsΛ]ds, (2.13)bn =2Λ∫ Λ0ϕ(s, 0) sin[2pinsΛ]ds, n = 1, 2, 3, ... (2.14)2.2.2 Resistive force theoryIn low Reynolds number swimming of a slender filament in a Newtonian fluids, the resistiveforces are linearly dependent on the local velocity. The force per unit length exerted by thefluid on the swimmer body at location s and time t is given byf(s, t) = −KTu · tt−KN (u− u · tt), (2.15)where KN and KT are, respectively, the normal and tangential resistive coefficients. The self-propulsion of elongated filaments is possible because of drag anisotropy (KN 6= KT ). A detaileddiscussion on this property can be found in the review paper by Lauga and Powers [3]. Recentexperimental studies of direct force and motion measurements on undulatory microswimmers inviscous fluids find excellent agreement with RFT predictions [33, 34]. The ratio rK = KN/KTvaries with the slenderness (L/r) of the body. In the limit of an infinitely slender body, L/r →∞, rK → 2, which is the value adopted in this study.For undulatory locomotion in dry granular media, we only consider the slow motion regimewhere grain-grain and grain-swimmer frictional forces dominate material inertial forces [20].The motion of the swimmer is confined to the horizontal plane such that the change of resistancedue to depth is irrelevant. In this regime the granular particles behave like a dense frictionalfluid where the material is constantly stirred by the moving swimmer [25]. The frictional forceacting tangentially everywhere on the surface of a small cylindrical element is characterizedby CF , which is refered to as the flow resistance coefficient [20]. The other contribution tothe resistive forces is the in-plane drag-induced normal force, which is characterized by CS .Note that CS is a constant because the drag is independent of the velocity magnitude. Thenormal resistive coefficient C⊥ depends on the orientation (ψ) of the element with respect tothe direction of motion (Fig. 2.1). In other words, the resistive force exerted by the granular2For numerical computations, the number of modes n∗ has to be finite. In this chapter, n∗ = 100 as detailedin Appendix A82.2. Mathematical Formulationmaterial on the swimmer per unit lengthf(s, t) = −C‖uˆ · tt− C⊥(uˆ− uˆ · tt), (2.16)whereC‖ = 2rCF , (2.17)C⊥(ψ) = 2rCF +2rCS sinβ0sinψ= C‖(1 +CS sinβ0CF sinψ), (2.18)tanβ0 = cot γ0 sinψ and γ0 is a constant related to the internal slip angle of the granularmedia[20]. Although a complete physical picture of the dependence of C⊥ on the orientationψ remains elusive, the application of the granular RFT proves to be effective. Several studieshave applied the granular RFT to study the locomotion of sand-swimming animals and artificialswimmers and found good agreement with experiments and numerical simulations [23, 25]. Adetailed discussion about the effectiveness of granular RFT on modelling sand-swimming canbe found in a review article by Zhang and Goldman [25].An important parameter characterizing the response of dry GM to intrusion is the volumefraction φ, which is defined as the ratio of the total volume of the particles divided by theoccupied volume. The level of compaction affects drag response as closely packed (high φ)GM expands to flow while loosely packed (low φ) material would consolidate [20]. The dragparameters CS , CF and γ0 depend on the volume fraction of the GM. In our study, we refer tothe GM with φ = 0.58 as loosely packed (LP) whereas φ = 0.62 as closely packed (CP). Thenumerical values of the drag parameters are adopted from the paper by Maladen et al. [20],where the forces at a fixed depth of 7.62 cm were measured by towing a cylinder of stainlesssteel. The drag parameters are presented in Table 2.1.Packing φ CS ,N/m2 × 10−4 CF ,N/m2 × 10−4 γ0,degreeLP 0.58 0.51 0.28 13.84CP 0.62 0.77 0.59 12.21Table 2.1: The parameters of the resistive force model in LP and CP GM as obtained byMaladen et al. [20]Without external forcing, the self-propelled filament satisfies force-free and torque-free con-ditions:F =∫ L0f(s, t) ds = 0, (2.19)T =∫ L0[x(s, t)− x(0, t)]× f(s, t) ds = 0. (2.20)The granular RFT exhibits the symmetry property that u→ −u results in f → −f . Com-92.2. Mathematical Formulationbining this symmetry with the kinematics of the undulatory locomotion (see Sec. 2.2.1), onecan show that the velocities −x˙(0, t) and −θ˙ are solutions to the instantaneous motion undera reversal of the actuation direction (c → −c) provided that x˙(0, t) and θ˙ are solutions to theoriginal problem (without reversal of the actuation). This symmetry is of course present inviscous RFT and this commonality, as we shall show, leads to qualitatively similar swimmingbehaviors.2.2.3 Swimming efficiencyThe instantaneous swimming speed of the filament is given by x˙(0, t), and the mean swim-ming velocity is defined as U = 〈x˙(0, t)〉 = Uxex + Uyey with the magnitude U = ‖U‖. Theangle brackets 〈...〉 denote time-averaging over one period T . The efficiency of the undulatorylocomotion for a given deformation wave is defined by the ratio of the power required to dragthe straightened filament through the surrounding substance to the power spent to propel theundulating body at the same velocity [35]. Hence, the efficiency for undulatory swimming ofslender filaments in viscous fluid (ηf ) and granular substance (ηg), respectively, areηf =KTLU2P, ηg =C‖LUP, (2.21)whereP =〈∫ L0f(s, t) · u(s, t) ds〉· (2.22)The optimal swimming can then be interpreted as either swimming with the maximum speedat a given power or swimming with the minimum power at a given speed.2.2.4 WaveformsWe consider two typical planar waveforms that have been well studied in Newtonian swim-ming: the sinusoidal waveform, and the sawtooth waveform (Fig. 2.2). The sinusoidal waveformcan be described by its Cartesian coordinates3:Y = b sin k(X +X0), (2.23)where k = 2pi/λ is the wave number, kX0 is the initial phase angle of the waveform, and b thewave amplitude. The dimensionless wave amplitude is defined as = kb.The sawtooth waveform, which consists of straight links with a bending angle β (ϕ = ±β/2),3The local tangent ϕ(s, t) of the Cartesian sine can be obtained using the relation ds =√1 + Y 2X dX combinedwith Eq. (2.7). Then we obtain the Fourier coefficients of the Cartesian sine using fsolve because the functionsare nonlinear.102.3. Bodies of infinite lengthcan be described asY =2bpiarcsin[sin k(X +X0)], (2.24)The dimensionless amplitude = kb = (pi/2) tan(β/2).(a) (b)b bsmallFigure 2.2: Undulating filaments with a single wave (N = 1). (a): sinusoid, kX0 = 0; (b)sawtooth, kX0 = 02.3 Bodies of infinite lengthFor bodies of infinite length (L → ∞), the swimming motion is steady and unidirectional,and hence ˙θ(t) = 0. Without loss of generality, we assume the filament propagates the defor-mation wave in the positive x-direction. Then the velocity of a material point on the body canbe written asu = −Uex + V ex − ct, (2.25)where U is the swimming speed [35]. For an infinite swimmer, the unidirectional swimmingvelocity for a given waveform can be obtained from only the force balance in the x-direction,F · ex = 0, over a single wavelength,∫ Λ0(CS sinβ0sinψ+ CF)uˆ · ex ds−∫ Λ0CS sinβ0sinψ(uˆ · t)t · ex ds = 0. (2.26)The above integral equation can be solved for U numerically for a given waveform in generalbut is analytically tractable in certain asymptotic regimes, which we discuss below.2.3.1 Optimal shape: numerical resultsA natural question for swimming organisms is how their swimming gaits evolve under thepressure of natural selection [36], since being able to swim does not necessarily mean one does itefficiently. The understanding of optimal swimming may reveal nature’s design principles andguide the engineering of robots capable of efficient self-propulsion. As a response, the optimalstrategies of several Newtonian swimming configurations have been studied. Becker et al.[37]determined the optimal strategy of Purcell’s three-link swimmer under constant forcing and112.3. Bodies of infinite lengthminimum mechanical work. Tam and Hosoi [38] improved the swimming speed and efficiencyof the optimal strategy of Purcell’s three-link swimmer by allowing simultaneous rather thansequential movement of both hinges (kinematic optimization). Using viscous RFT, Lighthillshowed that the optimal flagellar shape has constant angle between the local tangent to theflagellum and the swimming direction [9]. In 2D, the sawtooth profile with a tangent angleϕ ≈ ±40◦ (bending angle β ≈ 80◦) was found to optimize the swimming efficiency of an infinitelength swimming filament. Alternatively, this solution can be obtained through a variationalapproach [32]. In 3D, Lighthill’s solution leads to an optimal shape of a rotating helix. Morerecently, Spagnolie and Lauga studied the optimal shapes for both finite and infinite elasticflagellum by incorporating physical constraints such as bending and sliding costs [32]. Inspiredby the investigations of optimal strategies for Newtonian swimming, we study the optimal shapefor infinite swimmers in granular substrates using resistive force theory.For bodies of infinite length, the optimal shape is time, scale and phase invariant [32].Therefore, we take Λ = L = 1 and consider the optimization for t = 0. In other words, thelocal tangent angle for the optimization problem would beϕ(s, t = 0) =n∗∑n=1an cos(2pins). (2.27)We consider the optimal filament shape by maximizing the swimming efficiency η defined inSec. 2.2.3. Once the local tangent angle is obtained, the shape itself can be recovered byintegration. The numerical methods used in this optimization can be found in Appendix A.1.0 pi/2 pi 3pi/2 2pikx−101kyLighthillLPCPFigure 2.3: Optimal shapes in terms of swimming efficiency for an infinite filament in a granularsubstrate (LP, CP) and Newtonian fluid. The spatial coordinates are scaled to the same wavelength. For loosely packed granular material, the optimal shape is almost the same as theanalytical result of Lighthill’s in Newtonian fluid.The optimal shapes found by maximizing the swimming efficiency are presented in Fig. 2.3for a LP granular substrate (red dashed line), a CP granular substrate (blue dash-dot line), and122.3. Bodies of infinite lengtha viscous Newtonian fluid (black solid line) as a comparison. First, it is interesting that theoptimal shape stays as sawtooth despite the nonlinearity in the resistive force model of granularsubstrates. The optimal bending angles for LP and CP granular media are, respectively, β ≈ 80◦and β ≈ 87◦. The associated efficiencies of the optimal shapes are around 0.56 for LP and 0.51for CP granular substrates, which are much greater than that of Newtonian swimming. In spiteof the difference in the surrounding media, the optimal bending angle for granular substratesand viscous Newtonian fluids lie within the same range; in particular, the optimal sawtooth inLP closely resembles that in Newtonian fluids.We argue that it is not surprising that the sawtooth waveform is optimal in both the viscousRFT and the nonlinear granular RFT. Given an angle that maximizes the efficiency of a localelement. Without any penalty, the globally optimal shape would be the one that is locallyoptimal everywhere along the body. As a result, a local resistive force model should exhibitan optimal shape of a certain sawtooth waveform. Using this argument, we can simply dropthe integration (or assume it is a sawtooth) in Eq. (2.26) and consider the local optimality.The local optimal angle obtained is indeed the same as that found using numerical globaloptimization (see Sec. 2.3.2).The existence of a locally optimal tangent angle ϕ originates from the physical picture in-troduced by the drag-based propulsion model [3] (Fig. 2.1). Let ud = udey be the transversedeformation velocity of an infinite swimming filament. Then a propulsive force, which is per-pendicular to the direction of the deformation velocity, generated by this deformation can begiven byfprop = −(C⊥(ψ)− C‖) sinϕ cosϕex. (2.28)Therefore, the propulsive force arising from a local deformation of the filament scales with itsorientation assinϕ cosϕ√tan2 γ0 + cos2 ϕ, (2.29)the maximum of which is achieved when ϕ ≈ 64◦. However, as the tangent angle increases,the power consumption of the swimming filament increases. As a result, the swimmer tends toreduce the tangent angle to decrease the energy expenditure while maintaining a relatively highpropulsive force. It is the interplay of these two factors that determines the optimal tangentangle.2.3.2 Sawtooth and sinusoidThe swimming speed of an infinite sawtooth in viscous fluids can be expressed asUV=1− cosβ3− cosβ · (2.30)132.3. Bodies of infinite lengthFor a sawtooth profile in granular substrates, although an explicit analytical solution cannotbe extracted, an implicit algebraic equation for the swimming speed U can be obtained sincethe local resistive forces do not vary along the body:(CS sinβ0sinψ+ CF)uˆ · ex − CS sinβ0sinψ(uˆ · t)t · ex = 0, (2.31)where t · ex = cos(β/2). We then solve Eq. (2.31) numerically (see Appendix) with the sameconvergence criterion as in the optimization (Sec. 2.3.1). For a sinusoidal wave in granularmedia, a simplification like Eq. 2.31 is not available and we therefore directly solve Eq. (2.26)with the numerical method outlined in the Appendix.For small amplitude sawtooth waveforms ( 1), or small bending angle β, we obtain anasymptotic solution of the swimming speed U . Note that the swimming speed is invariant undera phase shift of pi, which is equivalent to a sign change in the amplitude: → −. Assuming aregular expansion in , this symmetry argument leads to a quadratic scaling of the swimmingspeed in the wave amplitude[4]UV∼ 4 cos γ0CSpi2CF2· (2.32)When the bending angle is large, another asymptotic limit can be obtained. The swimmingspeed U/V approaches a constant as β → pi and analytically we find thatUV∼ CSCS + CF tan γ0· (2.33)One can also show that this large amplitude asymptotic limit for a sawtooth equals that of asinusoidal wave. For small amplitude sinusoidal waveforms, however, the nonlinearity of theshape and the resistive forces results in a non-uniform integral and a slowly converging asymp-totic series. To leading order, the swimming speed U/V scales as 2/ ln(1/||), which does notagree well with the numerical results even for < 0.1 as the higher order terms being truncatedare not significantly smaller. We present the small and large amplitude asymptotic solutionsfor the granular swimming of a sawtooth profile in Fig. 2.4(a). The asymptotic solutions agreewell with the numerical solutions even for wave amplitudes close to one. Fig. 2.4(b) showsthe efficiency of swimming as a function of the bending angle for an infinite sawtooth in bothgranular media and viscous fluids. For swimming efficiency, a global maximum in bending angleexists for both viscous and granular swimming. Note that the optimal angles obtained here areequal to those obtained via the global optimization (Sec. 2.3.1).In Fig. 2.5, we compare the swimming speed and efficiency of sawtooth and sinusoidalwaveforms in both GM and Newtonian fluids as a function of the wave amplitude . In bothGM and Newtonian fluids, the swimming speed of a sawtooth is only slightly different from thatof a sinusoid with the same dimensionless amplitude. This small difference indicates that theeffects of the local curvature variations are not significant in both the granular and viscous RFT.142.3. Bodies of infinite length0 0.5 1 2 3 4 ∞ǫ00.30.60.9U/VLPCPNewtonian(a)Asymptotic0 π/4 π/2 3π/4 πβ0 0.5 1 2 3 4 ∞ǫ00.20.40.6ηLPCPNewtonian(b)0 π/4 π/2 3π/4 πβFigure 2.4: (a): Swimming speed of infinite sawtooth waveforms as a function of amplitude (or bending angle β) in granular material and Newtonian fluids. The dashed lines indicate thesmall and large amplitude asymptotic solutions. (b): Efficiency of infinite sawtooth waveformsas a function of amplitude (or bending angle β) in granular material and Newtonian fluids.0 1 2 3 4ǫ00.30.60.9U/VCPLPNewtonian(a)sawtoothsinusoid0 1 2 3 4ǫ00.10.20.30.40.50.6η CPLPNewtonian(b)sawtoothsinusoidFigure 2.5: A comparison of the swimming speed (a) and efficiency (b), as a function of waveamplitude for sawtooth and sinusoidal waveforms in granular substrates and Newtonian fluids.Although the sawtooth is found to be the mathematically optimal shape, the undulatory gaitof a sandfish resembles a smooth sinusoidal waveform [20]. The slight difference in swimmingperformance between the two waveforms presented in this section might justify the adoption of asinusoidal waveform instead of the mathematically optimal sawtooth waveform, since the kinksin the sawtooth may involve other energetic costs associated with bending and the deformationof the internal structure of the body [32].152.4. Bodies of finite length2.4 Bodies of finite lengthThe infinite swimmer model only enforces a force balance in one direction and hence aswimmer is confined to swim only unidirectionally without any rotation. In reality, however, aswimmer has a finite size and more complex swimming kinematics, including transverse motionrelative to the wave propagation direction and rotation. Previous studies employed slender bodytheory to investigate the swimming motion of finite filaments in a viscous Newtonian fluidand their swimming performance in relation to number of wavelengths and filament length[14, 32, 39, 40]. In this section, we investigate the swimming characteristics of finite-lengthsinusoidal swimmers in a granular medium and compare with their Newtonian counterparts.The numerical methods implemented to solve the equations of motion of a finite length swimmerare given in Appendix A.2.2.4.1 Geometries(a) (b)✏ = 0.2✏ = 1✏ = 4Figure 2.6: Shapes of swimming finite length single wave (N = 1) sinusoidal filaments fordifferent wave amplitude . (a): the odd sine configuration, with kX0 = 0, (b): the even cosineconfiguration, with kX0 = pi/2. The waveforms are rescaled to the same wave length for bettercomparison.For an undulating sinusoidal filament, the initial shape of the swimmer is determined bythe number of waves N , the wave amplitude , and the initial phase angle kX0 (Eq. (2.23)).The two specific categories of shapes that possess odd or even symmetry for a single wavesinusoidal swimmer are shown in Fig. 2.6. A swimmer in an odd configuration is the one thathas point symmetry about the midpoint of the filament as seen in Fig. 2.6(a), while an evenconfiguration is the one that possesses mirror symmetry about the vertical line through themidpoint as in Fig. 2.6(b). In our paper, the shapes shown in Fig. 2.6(a) are referred to asodd sine swimmers, while even cosine swimmers are those shown in Fig. 2.6(b). Note that aneven sine swimmer would be the one that has the number of waves N ∈ {1/2, 3/2, 5/2, ...} anda phase angle kX0 ∈ {0,±pi,±2pi, ...}; an even cosine swimmer is the one that has the numberof waves N ∈ {1, 2, 3, ...} and a phase angle kX0 ∈ {±pi/2,±3pi/2, ...}.162.4. Bodies of finite length2.4.2 Pitching, drifting and reorientationUnlike the swimming of an infinite length undulatory swimmer whose motion is steadyand unidirectional, the locomotion of a finite filament may also experience net motion normalto the initial direction wave propagation direction, also referred to as drifting, and unsteadyrotational motion, known as pitching. Here we characterize in GM the re-orientation of a finiteswimmer that results in drifting, and the dependence of swimming performance on pitchingmotion, previously reported to diminish performance in viscous Newtonian media [32, 40].−4 −3 −2 −1 0 1x/L−0.500.511.5y/Ldirection of motionFigure 2.7: Trajectory of the head x(0, t) (black solid lines) and trajectory of the swimmercentroid (dotted lines) for swimming finite sinusoidal filaments with N = 1 and = 1 thatpossess odd/even symmetry at t = 0 in loosely packed GM. The filament swims towards theleft when the wave propagates to the right. If the configuration possesses even symmetry itdoes not undergo a net reorientation.For an even symmetry filament in viscous fluids, Koehler et al. [40] showed that the velocityof the center of mass is along the centerline of the waveform, hence the net drifting is zero. Thisargument relies on the kinematic reversibility of Stokes flow: reflection about the vertical line isequivalent to a time reversal (or reversing the direction of the actuation), so the instantaneousswimming is identical to the mirror reflection of its time-reversal, and the linearity requires thereverse of velocity due to time-reversal, thus one can show that the transverse component ofthe velocity is zero. As a result, the net displacement in one period for a filament starts withthe even configuration is along the initial waveform centerline.Although the granular RFT is nonlinear, the aforementioned symmetry property (u →−u⇒ f → −f , see Sec. 2.2.2) means that the same argument for an even symmetry swimmercan be made in GM. Therefore, zero net transverse motion is achieved if the swimmer starts withan even symmetry, which is also corroborated by the numerical simulation. Fig. 2.7 shows thehead trajectories of two swimming sinusoidal filaments with the same wave amplitude ( = 1),one starts with even symmetry while the other starts with odd symmetry. The net displacementof the even cosine swimmer is in the negative x-direction, which is the opposite direction of thewave propagation at t = 0. The odd sine swimmer, however, appears to be drifting upwards to172.4. Bodies of finite lengththe positive y-direction through time.0 1 2 3 4ǫ−π/2−π/40π/4π/2kX0(a) |〈θ〉 − θ0|LP05◦10◦15◦20◦0 1 2 3 4ǫ−π/2−π/40π/4π/2kX0(b) |〈θ〉 − θ0|CP05◦10◦15◦20◦0 1 2 3 4ǫ−π/2−π/40π/4π/2kX0(c) |〈θ〉 − θ0|Newtonian05◦10◦15◦20◦0 1 2 3 4ǫ05◦10◦15◦20◦|〈θ〉−θ 0|NewtonianCPLPkX0 = 0, π, 2π, ...(d)Figure 2.8: Parametric plots for the magnitude of the reorientation angle |〈θ〉 − θ0| for a singlewave (N = 1) sinusoid in (a) loosely packed GM, (b) closely packed GM and (c) Newtonianfluids. (d): Plots of |〈θ〉 − θ0| against the wave amplitude for the odd sine configuration inGM and Newtonian fluids. |〈θ〉 − θ0| is periodic with a period of pi.The swimming behavior presented in Fig. 2.7 can be understood by examining the periodicinstantaneous motion of the swimmer. In the moving frame, or the Lagrangian frame, theinstantaneous motion of the swimmer can be viewed as being pulled through a waveform-shapedtube [40]. This motion, in turn, causes rotation and translation of the Lagrangian frame. Theinstantaneous rotation of the Lagrangian frame is described by θ(t), which is periodic dueto the periodicity of the wave propagation. The average of θ(t) over one period, denoted as〈θ〉, describes the average swimming direction. This angle 〈θ〉 is the same in every periodwhich results in a straight line trajectory on average. If a filament, starts with an odd (even)configuration at t = 0 (if aligned with the x-axis then θ0 = 0), it would possess even (odd)symmetry at t = T/4. Thus the filament alternates between even symmetry and odd symmetryafter successive time steps of T/4. In this viewpoint, 〈θ〉 − θ0 characterizes the amount of timet1 required for the filament to reorient itself such that it reaches an even symmetry. After that,the swimmer would move in the direction of the waveform centerline at t = t1. For a fixednumber of waves N and amplitude , the odd configuration requires the largest amount of time(T/4) to reach an even symmetry, therefore has the largest angle of reorientation. Note thatthe angle of reorientation should be distinguished from pitching of the swimmer, which is theinstantaneous rotation of the swimmer about its waveform centerline.182.4. Bodies of finite length0 1 2 3 4ǫ05◦10◦15◦20◦θ mpNewtonianLPCPFigure 2.9: Maximum instantaneous pitching angle θmp as a function of the wave amplitude for single wave (N = 1) sinusoidal swimmers in GM and Newtonian fluids.In Fig. 2.8, we present parametric plots of absolute value of the angle of reorientation|〈θ〉 − θ0| by varying the wave phase angle kX0 and the amplitude in both GM and viscousfluids. The number of waves is fixed as N = 1, which approximates the shape of an undulatingsandfish body [20]. Note that a phase shift of pi would result in a reversal of the directionof the transverse motion, hence the sign of 〈θ〉 − θ0. In both GM and Newtonian fluids, themaximum in |〈θ〉 − θ0| is obtained when the filament possesses an odd symmetry at t = 0, i.e.,kX0 ∈ {0, pi, 2pi, ...}. For shapes that possess even symmetry, namely, kX0 ∈ {pi/2, 3pi/2, ...},zero transverse motion is observed. Within our parameter range, a maximum in |〈θ〉 − θ0| isachieved around an intermediate value of the amplitude for a given phase angle. As an example,the variation of |〈θ〉−θ0| with the amplitude for the odd configuration is shown in Fig. 2.8(d).The largest amount of reorientation of an odd swimmer is achieved when ≈ 1 − 1.2 in GMwhile ≈ 2.2 in viscous fluids. We also note that the angle of reorientation decreases with theincreasing of wave amplitude in the large amplitude region ( > 2).Although the transverse motion of the even configuration is minimal, the instantaneouspitching, θ(t)−〈θ〉, which generally diminishes performance, can be significant. Multiple metricshave been used to characterize pitching of a swimmer [32, 40], here we use the maximal amountof instantaneous pitching a swimmer can experience in one cycle of its motion θmp = |θ(t) −〈θ〉|max. Fig. 2.9 shows the maximal instantaneous pitching angle θmp for single wave sinusoidalswimmers in GM and Newtonian fluids. The maximal instantaneous pitching angle of a singlewave sinusoid goes up to about 15◦ in loosely packed GM while around 19◦ in closely packedGM.The instantaneous pitching of the swimmer results in a tortuous motion with a net swimmingspeed smaller than that of an infinite sinusoid. For a fixed number of waves and wave amplitude,a phase shift only leads to a variation in the direction of swimming. In other words, the velocitymagnitude U is independent of kX0 but the x and y components vary. From a control point of192.4. Bodies of finite lengthview, one can change the phase angle of an artificial sinusoidal swimmer to obtain the desireddirection of swimming.2.4.3 Swimming performanceThe two typical metrics for swimming performance used in the literature are the dimen-sionless swimming speed U/V and the swimming efficiency η, see Eq. (2.21). For a sinusoidalswimmer, the performance depends on the dimensionless amplitude and the number of wavesN . Note that the initial phase angle kX0 does not affect the two performance metrics. Thedesired motion of a finite swimmer is its translation, therefore the optimization of a finitesinusoidal filament requires minimizing pitching.0 1 2 3 4ǫ00.30.60.9U/VLP(a)N = 1N = 2N = 3N =∞0 1 2 3 4ǫ00.30.60.9U/VCP(b)N = 1N = 2N = 3N =∞Figure 2.10: Swimming speed U/V as a function of the dimensionless amplitude for differentnumber of waves N in (a) loosely packed GM and (b) closely packed GM. The solid lines denotethe swimming speed of an infinite sinusoid.For an undulatory finite filament in viscous fluids, several studies have characterized theswimming performance and optimal strategies. Spagnolie and Lauga reported that the lo-cal maxima in swimming efficiency occur for around half-integer number of waves (N ≈3/2, 5/2, ...,) when the bending cost is small [32]. Later studies by Koehler et al. [40] andBerman et al. [41] also showed that, for a sinusoidal swimmer, local maxima in performanceare achieved for close to half-integer number of waves where pitching is small.We first verify that the swimming velocity (Fig. 2.10) and efficiency (Fig. 2.11) of a finitesinusoidal swimmer in GM both converge to that of an infinite sinusoidal swimmer as thenumber of waves N increases. For a single wave sinusoid (N = 1) in loosely packed GM, theoptimal dimensionless amplitude that maximizes the efficiency is ≈ 1.68. As the number ofwaves increases, the optimal dimensionless amplitude approaches that of an infinite sinusoid( ≈ 1.33). Similar observations can be made for closely packed GM. We also observe that for agiven dimensionless amplitude , the difference in the swimming velocity (or efficiency) betweena short swimmer (N = 1) and an infinite swimmer can be associated with the pitching motion:202.4. Bodies of finite length0 1 2 3 4ǫ00.10.20.30.40.50.6ηLP(a)N = 1N = 2N = 3N =∞0 1 2 3 4ǫ00.10.20.30.40.50.6ηCP(b)N = 1N = 2N = 3N =∞Figure 2.11: Swimming efficiency η as a function of the dimensionless amplitude for differentnumber of waves N in (a) loosely packed GM and (b) closely packed GM. The shaded regionsrepresent the observed values of for lizards reported in the literature [20, 21].1 1.5 2 2.5 3N0.20.30.40.5η(b)LPCP1 1.5 2 2.5 3N0.20.30.40.50.6U/V(a)LPCPFigure 2.12: (a) Swimming speed as a function of the number of waves in GM. (b) Swimmingefficiency as a function of the number of waves in GM. The dimensionless amplitude is fixed( = 1).the largest difference in swimming speed (or efficiency) between the N = 1 and N = ∞swimmers occurs in the region ≈ 1 in Figs. 2.10 and 2.11, which is also the region wherepitching is the most significant (Fig. 2.9).For a given waveform, the amount of pitching can be altered by changing the number ofwaves N . We investigate in Fig. 2.12 the dependence of the performance metrics on the numberof waves for a finite sinusoidal swimmer, keeping dimensionless amplitude fixed at = 1. Ratherthan approaching the swimming velocity (or efficiency) of the corresponding infinite sinusoidmonotonically with increasing number of waves, the swimming speed and efficiency exhibitlocal maxima and minima. Similar to the Newtonian case, the local maxima in efficiency andswimming speed occur for the number of waves close to (but not equal) half-integers. The212.5. Conclusion1 1.5 2 2.5 3N05◦10◦15◦20◦θ mpLPCPFigure 2.13: Maximum instantaneous pitching angle as a function of the number of waves inGM. The dimensionless amplitude is fixed ( = 1).volume fraction of the GM has no significant influence on the number of waves where localmaxima in swimming performance occur. As shown in Fig. 2.12, the first local maximum inswimming performance for the number of waves greater than one occurs around N ≈ 1.4. Themaxima in swimming performance are associated with minimal pitching as shown in Fig. 2.13.Finally we note that although both the first two local maxima have minimal pitching (Fig. 2.13),the swimmer with more number of waves (N ≈ 2.5) still displays better swimming performance,which can be attributed to a smaller bobbing motion [40] (the relative motion of the center ofmass of the swimmer to the net swimming direction) for the swimmer with more number ofwaves.Finally, we relate our findings to biological observation; we show, in the shaded regionsof Fig. 2.11, the observed dimensionless amplitude (amplitude-to-wavelength ratio) for lizardsreported in the literature ( = 1.20−1.38) [20, 21]. We see in the case of both loosely-packed andclosely-packed granular media, that the biologically observed range of wave amplitudes samplehigh efficiencies not far from optimal ( ≈ 1.69 for LP and ≈ 1.95 for CP for N = 1). Sincethe efficiency peak is broad, a swimmer may adopt a close-to-optimal shape at the expense ofonly a modest drop in swimming efficiency to address other constraints (such as bending costsor internal dissipation).2.5 ConclusionIn this chapter, we have investigated locomotion of slender filaments in granular media usinga resistive force theory proposed by Maladen et al. [20]. While previous work focused on infiniteswimmers (or 1-D swimming) in reality a swimmer has a finite size, which leads to more complexswimming motion. By taking into account full force and torque balances, a finite swimmer isno longer only confined to swim in a straight trajectory. The orientation of the swimmer can be222.5. Conclusioncontrolled by adjusting the features of the waveform such as the amplitude, phase, and numberof wavelengths, allowing a swimmer to move from an initial position to a final destination via amore complex, designated trajectory. These degrees of freedom enable the control of swimmerswithout the use of any external fields to actively steer the swimmer. Our studies characterizethis complex swimming motion in granular media, which may be useful for the developmentof programmable and efficient autonomous locomotive systems in such environments, but alsosuggest that swimmers in nature are themselves closely tuned for optimality.We also find that undulatory locomotion of filaments in granular media is distinctly similarto that in viscous fluids. We compared a number of observations made for swimming in viscousfluids with RFT both for finite and infinite swimmers and found qualitatively similar behaviorusing granular resistive force theory despite the nonlinearity of the force law. The reason islargely down to two distinct similarities. The first, is that both laws are still local and thusignore interactions of distinct parts of the body through the medium in which they swim.Ultimately this leads to finding that a sawtooth profile optimizes locomotion in both viscousfluids and granular media. The second, is that both force laws display the symmetry thatu → −u results in f → −f . This leads to a kinematic reversibility in both cases, where areversal of the wave speed leads to a reversal of the translational and rotational motion of theswimmer, and hence a myriad of qualitatively similar behaviors that we have explored andquantified in this chapter.23Chapter 3Optimal flexibility of a drivenmicrofilament in a viscous fluid3.1 IntroductionIn many situations of biological locomotion in fluids, such as flapping birds, swimmingfish and beating flagella, propulsive thrust is achieved by periodic motion of a flexible bodyor appendage [3, 36, 42]. Elastic deformations of flexible parts on a swimmer due to thefluid-structure interaction can lead to significantly improved propulsive performance acrossReynolds numbers [4, 43, 44]. However, the physics of swimming at small scales relevant tomicroorganisms are fundamentally different from that of the mesoscopic swimming of flexiblebodies.For swimming microorganisms such as flagella and bacteria, the dominance of viscous forcesover inertial effects leads to the time-reversible Stokes equations governing the fluid motion. Inthis low Reynolds number regime, Purcell’s famous “scallop theorem” states that a reciprocalmotion (a deformation that exhibits time-reversal symmetry) cannot generate any net propul-sive thrust [5]. In order to break the constraint of time-reversibility, many microorganismsincluding flagellated bacteria and spermatozoa achieve self-propulsion by passing deformationwaves along their flexible bodies [3, 6, 7]. Advances in fabrication technologies at small scalesallow the recent rapid development of synthetic micro-propellers capable of swimming at speedscomparable with microorganisms. In particular, slender flexible filaments have been employedto enable locomotion at small scales [45–48].A rigid filament driven at one end cannot propel itself because the motion is reciprocal. Byintroducing flexibility in the filament, the coupling between the viscous and elastic forces pro-duces deformation along the filament that can lead to propulsion [49, 50]. For a given actuationfrequency and filament length, an optimal bending stiffness of the filament can be determinedto produce the largest propulsive force [4, 49, 51]. However, the possibility of further improvingthe propulsion by allowing variable flexibility along the filament remains largely unexplored.Flying animals such as hoverflies and hummingbirds exhibit non-uniform flexibility distributionalong their wings, which can potentially enhance the propulsive performance [52, 53]. For aflapping wing in this high Reynolds number regime, Shoele and Zhu [54] have compared theperformance of several cases of nonuniform flexibility distributions and a recent work by Moore[55] has shown that optimal propulsion can be achieved by a highly localized flexibility arrange-243.2. Mathematical formulationment at the front of the wing using a torsional spring. At low Reynolds numbers, Maier et.al. [48] have constructed a flagellar bundle attached to a magnetic head using DNA tile-tubeassembly where they achieved an exponentially decreasing stiffness profile. They found thatfor several starting stiffnesses at the basal end, the swimmer with an exponentially decreasingstiffness profile outperforms the one that has uniform stiffness down along the flagellum un-der a rotational actuation using an external magnetic field. However, they did not attempt asystematic parameter study or an optimization.In this work, we consider the propulsive force generated by a boundary displacement-drivenpassive cantilever filament at low Reynolds number. The mathematical formulation allows vari-able bending flexibility, and we derive analytically the expressions for the propulsive force ofa filament with two segments of different bending flexibilities connected serially together, atorsional spring connected to a rigid rod and an arbitrarily continuous distribution of bendingflexibility using asymptotic analysis for small actuation amplitude. We show that, differentfrom the high Reynolds number case, the torsional spring arrangement does not optimize thepropulsion. From a numerical optimization, we show the optimal linear and quadratic distri-bution of flexibility. By considering two other different boundary conditions, we show that theoptimal flexibility arrangement can be qualitatively modified so that one can not simply extendthe optimality for one case to other cases where the boundary actuation mechanisms differ.This Chapter is organized as follows. We formulate the problem and review the classicalresults of a boundary-driven passive filament with uniform bending stiffness in Sec. 3.2. A fila-ment with two segments that have different bending stiffnesses is first considered (Sec. 3.3): weshow that this arrangement can achieve a higher propulsive thrust than the maximum of a uni-form case. In Sec. 3.4, we calculate the propulsion generated by a torsional spring arrangement.Then we consider a numerical optimization over several cases of continuous stiffness distribu-tion in Sec. 3.5. Next, we discuss the effect of boundary conditions by considering two differentmechanisms of actuation in Sec. 3.6, before concluding the work with remarks in Sec. 3.7.3.2 Mathematical formulationWe consider a slender cylindrical filament of length L and uniform radius a such that a Land assume the filament is elastic and inextensible. The deformation of the filament is assumedto be confined in the x-y plane. We define the position vector of a material point on thefilament neutral line relative to the laboratory frame as x(s, t), where s is the arclength alongthe filament with s ∈ [0, L] at time t. It is convenient to describe the shape of the filament bythe local tangent angle made with the x-axis as ψ(s, t) such thatxs = [cosψ, sinψ]T, (3.1)where the subscript s denotes differentiation with respect to s, namely xs ≡ ∂sx. The local unittangent and normal vectors at location s along the filament are defined as t and n respectively,253.2. Mathematical formulationwith t = xs. The local geometry is thus characterized by the Frenet-Serret formulas:ts = xss = κn, ns = −κt, (3.2)where κ = ‖xss‖ = ψs is the local curvature.3.2.1 Enthalpy functionalIn order to model the dynamics of the driven filament, we start from an enthalpy functionalE = 12∫ L0Ax2ss ds+12∫ L0σ(x2s − 1)ds, (3.3)where A = A(s) = EI is the bending stiffness which we allow to vary along the filament withE the Young’s modulus and I the second moment of inertia of the cross-section. The localinextensibility condition xs · xs = 1 is enforced by introducing the Lagrange multiplier σ(s, t).The internal elastic force density is determined by a variation δE with respect to a variation δxof the shape x. Noting thatδκ = δ‖xss‖ = δ(√x2ss)= n · δxss, (3.4)we haveδE =∫ L0(Aκn · δxss + σxs · δxs) ds. (3.5)Upon integration by parts, we obtainδE =[Aκn · δxs]s=Ls=0−[∂s (Aκn) · δx]s=Ls=0+[σxs · δx]s=Ls=0+∫ L0∂s[∂s(Aκ)n− (Aκ2 + σ)t] · δx ds. (3.6)So the internal elastic force density is given byfelastic = −δEδx= −∂s [∂s(Aκ)n− τt] , (3.7)where we have defined τ = σ + Aκ2. The boundary terms in Eq. (3.6) can be interpreted asexternal forces and toques applied at the two ends [56]. In other words,Text = Aκ, Fext = τt− ∂s(Aκ)n at s = L,Text = −Aκ, Fext = −τt + ∂s(Aκ)n at s = 0. (3.8)263.2. Mathematical formulation3.2.2 ElastohydrodynamicsWe describe the hydrodynamics of the viscous fluid by the resistive force theory, which is aleading-order approximation in the small filament aspect ratio a/L. The theory states that theviscous force per unit length on the filament is given byfvis = −(ξ⊥nn + ξ‖tt) · xt, (3.9)where the subscript t denotes differentiation with respect to time and ξ⊥ and ξ‖ are the normaland tangential resistive coefficients respectively. The local balance between the viscous andelastic forcesfvis + felastic = 0 (3.10)can be written asxt =(1ξ⊥nn +1ξ‖tt)· felastic. (3.11)From the variational formulation in Sec. 3.2.1 it follows that∫ Lsfvis ds =∫ LsδEδxds = [∂s(Aκ)n− τt] (s = L)− [∂s(Aκ)n− τt] , (3.12)orτt− ∂s(Aκ)n =∫ Lsfvis + Fext(L). (3.13)In other words,τ = t ·(∫ Lsfvis ds+ Fext(L)),∂s(Aκ) = −n ·(∫ Lsfvis ds+ Fext(L)). (3.14)Eq. (3.14) shows that τ acts as the physical tension along the filament. Following Eq. (3.11),we obtainxt =1ξ⊥n (−∂ss(Aψs) + ψsτ) + 1ξ‖t (ψs∂s(Aψs) + τs) . (3.15)Noting thattt = (− sinψ, cosψ)T ψt = ψtn, ns = −ψst, ts = xss = ψsn, (3.16)273.2. Mathematical formulationwe differentiate Eq. (3.15) with respect to the arclength and obtainψtn =1ξ⊥(−ψst) (−∂ss(Aψs) + ψsτ) + 1ξ‖t∂s (ψs∂s(Aψs) + τs)+1ξ⊥n∂s (−∂ss(Aψs) + ψsτ) + 1ξ‖ψsn (ψs∂s(Aψs) + τs) . (3.17)Finally, the equation of motion for the tangent angle ψ(s, t) readsψt =1ξ⊥(−∂sss(Aψs) + ∂s(ψsτ)) + 1ξ‖ψs (ψs∂s(Aψs) + τs) . (3.18)Another equation obtained is a partial differential equation (PDE) for tension τ(s, t), or equiv-alently from the inextensibility condition t · tt = 0:τss −ξ‖ξ⊥ψ2sτ = −∂s(ψs∂s(Aψs))−ξ‖ξ⊥ψs∂ss(Aψs). (3.19)Eqs. (3.18) and (3.19) determine the filament dynamics. Note that (−ψ, τ) satisfies Eqs.(3.18) and (3.19) if (ψ, τ) is a solution. Once the tangent angle ψ is solved, the filament shapecan be recovered by integrationx(s, t) = x(0, t) +∫ s0(cosψ, sinψ)T ds′, (3.20)where x(0, t) can be obtained from Eq. (3.15) by integration with respect to time evaluated ats = 0.3.2.3 Boundary conditionsThe filament dynamics governed by Eqs. (3.18) and (3.19) depends on the prescribedboundary conditions. In this work, we consider a boundary driven passive filament where oneend is oscillated in a controlled manner while the other end is free to move in the surroundingfluid. At the free end (s = L), we have force-free and torque-free boundary conditions,Fext(L) = [τt− ∂s(Aκ)n]s=L = 0, Text(L) = [Aκ]s=L = 0. (3.21)At the actuation end (s=0), we consider a harmonic oscillation with frequency ω of the trans-verse position of a cantilevered filament i.e.,y(0, t) = y0 sinωt, x(0, t) = 0, ψ(0, t) = 0. (3.22)283.2. Mathematical formulation3.2.4 Non-dimensionlizationWe non-dimensionalize the governing equations with respect to a length scale L, time scaleω−1, velocity scale Lω and a force scale4 A/L2. The resistance ratio γ, sperm number Sp andthe dimensionless oscillation amplitude , respectively, are defined asγ =ξ⊥ξ‖, Sp = L(ξ⊥ωA)1/4, =y0L. (3.23)The sperm number Sp compares the magnitude of viscous and elastic forces. For a givenoscillation frequency ω and length of the filament L, a larger sperm number indicates a moreflexible material. For a rigid filament, the sperm number Sp→ 0. To make analytical progress,we perform asymptotic analysis in the small amplitude oscillation limit, 1, to determinethe filament shape order by order.3.2.5 Propulsive thrustIn our study, we are considering the effect of flexibility distribution on the propulsive per-formance, hence we non-dimensionalize the thrust by L2ξ⊥ω which is hold constant. Theinstantaneous propulsive thrust is defined as the total hydrodynamic drag force in the directionof swimming:Fx = −ex ·∫ 10fvis ds (3.24)where ex is the unit vector in the swimming direction (x-direction in our formulation). Thenet propulsive thrust is the time-average of the instantaneous thrust over one period of theactuation, Fp = 〈Fx〉 with 〈...〉 =∫ 2pi0 (...) dt/2pi.For small amplitude oscillation, we write the net propulsive force as a regular series expan-sion in ,Fp = 2F (2)p +O(4). (3.25)For a continuous stiffness profile A(s) (dimensional), we haveF (2)p =γ − 1γSp40〈∫ 10ψ(1)∂2s (A∗∂sψ(1)) ds〉, (3.26)where Sp0 = L[ξ⊥ω/A(0)]1/4 is the sperm number evaluated using the stiffness at the basal end,A∗ = A(s)/A(0) and ψ = ψ(1) +O(2). Note that an order boundary actuation generates apropulsive force of O(2) to leading order. For A∗ ≡ 1, the propulsive force given by Eq. (3.26)reduces to the case of a uniform stiffness profile as given by previous work [4, 57].4Since the bending stiffness A can be varying along the filament, we use the value of A at a certain materialpoint on the body.293.3. Two-segment filaments0 2 4 6 8 10Sp0.000.050.100.150.20F(2)pFigure 3.1: Propulsive thrust generated by a cantilevered filament of uniform bending stiffnessunder a displacement actuation at one extremity. The two blue dots indicate two filamentswith different bending stiffnesses generating the same propulsive force.For the classical case of uniform bending stiffness (A = constant) along the filament, thevariation of propulsive force as a function of the sperm number is shown in Fig. 3.1. In thelow sperm number limit (Sp 1), the filament becomes a rigid rod undergoing a reciprocalmotion, which produces no net propulsion. For finite values of sperm number, flexibility enablesthe propagation of deformation waves along the body, which breaks the constraint of kinematicreversibility and generates net propulsive force. In the large sperm number limit, the dominanceof viscous force over elastic force suboptimally localizes the deformation of the filament aroundthe actuation end. Therefore, the propulsive force exhibits a maximum around an intermediatesperm number Sp ≈ 1.89 and asymptotes to a limiting value for large sperm numbers. Thedeforming shapes of a filament with uniform stiffness are shown in Fig. 3.2 for several differentsperm numbers.As indicated by the two blue dots in the proximity of the maximum propulsive force inFig. 3.4, one can have deforming filaments of two different sperm numbers generating the samepropulsive force. In other words, for a given actuation frequency and filament length one cangenerate the same propulsive thrust with filaments of two different bending rigidities (A1 andA2). Then the question is, can we still obtain the same propulsive force by serially connectingtwo segments that each individually generates the same propulsion, or can we enhance thepropulsion?3.3 Two-segment filamentsWe probe the potential advantages of non-uniform stiffness by connecting the segments withbending flexibilities A1 (at the actuation end) and A2 serially as diagrammed in Fig. 3.3. Therelative proportion of the segment with stiffness A1 is denoted as α, 0 ≤ α ≤ 1. We define the303.3. Two-segment filaments−0.10.00.1−0.10.00.1y0 0.2 0.4 0.6 0.8 1x−0.10.00.1Figure 3.2: Shapes of the displace driven cantilevered filament in different times over a period.From top to bottom, the sperm number increases (top:Sp = 0.5, middle:Sp = 1.89, bottom:Sp =8), namely, the filament becomes more flexible.ratio of stiffness β = A2/A1. In this case, the local tangent ψ(s, t) and tensile force τ(s, t) aresplit into two functions for the two segments, with ψnk denoting the tangent angle of the k-thsegment (k = 1, 2) at O(n). Noting that ψ ∼ , from Eqs. (3.18) and (3.19) one can show thatτ ∼ 2 [4]. This two-segment arrangement considered here is arguably one of the simplest casesof non-uniform flexibility distribution, and we shall show that the asymptotic analysis revealsthe advantage of this simple non-uniform flexibility arrangement over the uniform case.In the small amplitude limit, we seek a perturbative solution by using regular series expan-sions in ,ψ1(s, t) = ψ(1)1 + 2ψ(2)1 +O(3),τ1(s, t) = τ(0)1 + τ(1)1 + 2τ(2)1 +O(3),ψ2(s, t) = ψ(1)2 + 2ψ(2)2 +O(3),τ2(s, t) = τ(0)2 + τ(1)2 + 2τ(2)2 +O(3). (3.27)313.3. Two-segment filaments0 0.2 0.4 0.6 0.8 1α0.0800.0850.0900.0950.1000.1050.1100.115F2β = 0.239β = 4.186A1A2A1 A2↵(a) (b)↵Ls = 0Figure 3.3: Schematic of a serially connected filament. The segment at the actuation end has abending flexibility A1 and the segment at the free end has a bending flexibility A2, β = A2/A1.3.3.1 Elastic-elastic caseIn this section, we consider the case where both A1 and A2 are finite. The leading orderequations of motion are written asSp41∂tψ(1)1 + ∂4sψ(1)1 = 0,Sp41∂tψ(1)2 + β∂4sψ(1)2 = 0, (3.28)where Sp1 = L (ξ⊥ω/A1)1/4 is the sperm number of the material at the driven end.The boundary conditions at the two ends are given byψ(1)1 (0, t) = 0, ∂3sψ(1)1 (0, t) = −Sp41 cos t at s = 0,∂sψ(1)2 (1, t) = 0, ∂2sψ(1)2 (1, t) = 0 at s = 1. (3.29)In order to solve Eq. (3.28), 4 more boundary conditions are required. Note that thetangent angle ψ should be continuous, thus we have ψ(1)1 (α, t) = ψ(1)2 (α, t). Meanwhile, theinternal force and moments are expected to be continuous across the connecting point, i.e.∂sψ(1)1 (α, t) = β∂sψ(1)2 (α, t),∂ssψ(1)1 (α, t) = β∂ssψ(1)2 (α, t). (3.30)The last boundary condition is obtained from the continuity of the viscous force density, whichstates that ∂sssψ(1)1 (α, t) = β∂sssψ(1)2 (α, t).Since we are interested in the steady state solution, we write the solution asψ(1)1 (s, t) = Re{eith1(s)},ψ(1)2 (s, t) = Re{eith2(s)}, (3.31)323.3. Two-segment filamentswhere i =√−1 is the imaginary unit. Then Eq. (3.28) reduces to two ordinary differentialequations (ODEs), which can be solved analytically.The leading order propulsive force in this case is given byF (2)p =γ − 12γSp41〈(∂sψ(1)1 (0, t))2 −[(∂sψ(1)1 (α, t))2 − β(∂sψ(1)2 (α, t))2]〉. (3.32)Note that if β = 1, α = 0 or α = 1, the expression for propulsive force reduces to the result fora filament with uniform stiffness.0 0.2 0.4 0.6 0.8 1α0.130.150.170.190.21F(2)pβ = 0.144β = 6.959Figure 3.4: Propulsive thrust generated by two-segment filaments as a function of α at differentvalues of stiffness ratio β. For β = 0.144 (A2 < A1, more flexible materials at the free end),the maximum propulsive force generated is greater than the maximum achievable thrust of afilament with any uniform stiffness.For given sperm numbers of two materials indicated by the blue dots in Fig. 3.1, by varyingthe relative proportion of the two segments, a non-monotonic variation of the propulsive forceemerges (Fig. 3.4). The limits α = 0, 1 reduce to the case of uniform stiffness with a propulsiveforce F(2)p = 0.1358 as indicated by the horizontal dash-dotted line in Fig. 3.4. The solid linedenotes putting the more flexible material at the free end while the dashed line indicates puttingthe more flexible material at the driven end. It is interesting to note that the propulsive forcegenerated by the two-segment filament can be greater than the original uniform case (α = 0, 1)for both arrangements. For β = 6.959 , the filament outperforms the original uniform case inthe range α . 0.4. For β = 0.144, the two-segment filament outperforms the original uniformcase for almost all α.Furthermore, by putting a more flexible material at the free end i.e. β < 1 (e.g. β =0.144, A2 < A1, solid line, Fig. 3.4), the propulsive force generated by such an arrangement canbe greater than the maximum possible propulsive thrust (indicated by the horizontal dottedline in Fig. 3.4) achievable by a filament with any uniform bending stiffness. This enhancementin propulsion indicates that the filament prefers to have a more flexible part at the free end333.3. Two-segment filamentswith an optimal choice of α.3.3.2 Rigid-elastic caseThe limiting case from the elastic-elastic arrangement would be to consider an rigid-elasticfilament with A1 = ∞. In other words, the segment at the driven end is rigid. For the rigidpart, the local tangent angle is independent of s, i.e.0 = ψ1 = ψ1(t) = ψ2(α, t). (3.33)As a result, one only needs to solve the PDE for ψ2. Following similar procedure to the elastic-elastic case, we solve the governing equations to leading order. The two boundary conditionsat the free end is the same as the elastic-elastic case. The boundary conditions at α are givenby ψ(1)2 (α, t) = 0 and ∂3sψ(1)2 (α, t) = −Sp42 cos t. Note that Sp2 is the sperm number evaluatedusing the stiffness (A2, finite) of the flexible part. The leading order propulsive force in thiscase readsF (2)p =γ − 12γSp42〈(∂sψ(1)2 (α, t))2〉. (3.34)0 0.2 0.4 0.6 0.8 1α0.000.050.100.150.20F(2)p124 6Figure 3.5: Propulsive thrust generated by a cantilevered rigid-elastic filament under displace-ment oscillation as a function of α for different sperm numbers Sp2, which are indicated bythe numeric values around each line. The horizontal dash-dotted line denotes the maximumpropulsion of a filament with uniform stiffness.In Fig. 3.5, we present the leading order propulsive force as a function of the proportionof the rigid part for different sperm numbers of the flexible part. For small sperm numbers,the propulsive force decreases monotonically as α increases because the rigid part does notcontribute. For larger sperm numbers (Sp2 > 2), we observe a non-monotonic variation ofpropulsion as α increases from zero to one. It is interesting to note that the optimal propulsion343.3. Two-segment filamentsa rigid-elastic filament can achieve is around the maximum propulsion of a filament with uniformstiffness as indicated by the horizontal dash-dotted line in Fig. 3.5. In other words, a rigid-elastic filament cannot outperform the optimal uniform case, which means that an optimaltwo-segment arrangement should be an elastic-elastic one, i.e., both A1 and A2 remain finite.3.3.3 Elastic-rigid caseBefore discussing a filament with a continuously varying stiffness profile, it is interesting toinvestigate an elastic-rigid filament with A2 =∞. In other words, the segment at the free endis rigid. For the rigid part, the local tangent angle is independent of s, i.e.ψ2 = ψ2(t) = ψ1(α, t). (3.35)As a result, one only needs to solve the PDE for ψ1. Following similar procedure to the elastic-elastic case, we solve the governing equations to leading order. The two boundary conditions atthe driven end is the same as the elastic-elastic case. A force balance on the rigid part requiresthat∂2sψ(1)1 (α, t) = (α− 1)∂3sψ(1)1 (α, t)− ∂4sψ(1)1 (α, t)(12− α+ α22). (3.36)The fourth boundary condition is obtained via a torque balance on the rigid filament,−∂sψ(1)1 (α, t) +12(1− α)2∂3sψ(1)1 (α, t) +13(1− α)3∂4sψ(1)1 (α, t) = 0. (3.37)The leading order propulsive force is expressed asF (2)p =γ − 1γSp41〈ψ(1)1 (α, t)∂2sψ(1)1 (α, t)−12(∂sψ(1)1 (α, t))2+12(∂sψ(1)1 (0, t))2 + (1− α)ψ(1)1 (α, t)∂3sψ(1)1 (α, t)〉. (3.38)In Fig. 3.6, we present the variation of propulsive force generated by an elastic-rigid filamentas a function of α for different sperm numbers of the flexible part. We note that the optimalelastic-rigid filament cannot outperform the maximum of a uniform case. For a relatively smallsperm number (e.g., Sp1 < 2), we observe in general a monotonic increasing of the propulsiveforce as the proportion of the flexible part increases from 0 to 1. When α = 0, the filamentbecomes a rigid rod so that the propulsion generated is zero due to the time-reversibility of thenon-deforming motion. As α increases, the proportion that can deform due to the interactionbetween the viscous force and elastic force increases. As a result, one can observe an increasein the propulsive performance. However, for a larger sperm number of the flexible part, anon-monotonic behavior of the propulsive force is observed as α increases. Take Sp1 = 3 as anexample, the maximum propulsive force F(2)p ≈ 0.1875, which is achieved when α ≈ 0.039. As353.3. Two-segment filaments0 0.2 0.4 0.6 0.8 1α0.000.050.100.150.20F(2)pSp1 = 1Sp1 = 2Sp1 = 3Figure 3.6: Propulsive thrust generated by an elastic-rigid filament (A2 = ∞) as a functionof α for different sperm numbers of the flexible part. For small sperm numbers, a monotonicincreasing of the propulsive thrust is observed as α increases. For larger sperm numbers,however, the variation of propulsive thrust is non-monotonic. Note that the sperm number ofthe rigid part would be zero. The horizontal dash-dotted line denotes the maximum propulsiveforce achievable by a filament with uniform stiffness.the sperm number Sp1 increases, namely the flexible part becomes more flexible, the optimumpropulsive force is achieved at a smaller α. For sperm numbers larger than 2, the optimal elastic-rigid arrangement exhibits a localized flexibility around the actuation end while the propulsionis predominantly generated by the rigid part at the free end. In this scenario, the tiny flexiblebit at the driven end is effectively setting an optimal angle for the rigid part, which maximizesthe propulsion. From fundamental beam theory, we may scale the deflection δ (dimensional) ofthe connecting point asδ ∼ Fl3A1, (3.39)where l = αL is the dimensional length of the flexible bit and F the effective total hydrodynamicforce at this point. Noting that F ∼ L2ξ⊥ω, and the average angle at the connecting pointψ ∼ δ/l, we recover the physical scalingψ ∼ l2L2ξ⊥ωA1= Sp41α2. (3.40)Note that Sp41α2 is the effective sperm number of the flexible part upon a proper scaling of theelastic force, namely, we haveSpproper =(L2ξ⊥ωA1/l2)1/4. (3.41)363.4. Torsional springThese observations presented above lead us to considering the limiting case of an entirely rigidfilament with a torsional spring at the actuation end, which we shall show would be the optimalelastic-rigid arrangement upon an optimal choice of the spring constant.3.4 Torsional springWe consider a rigid filament of length L connected by a torsional spring with a springconstant C at the actuation end. In this scenario, the local tangent on the entire filament isindependent of s, i.e., ψ(t) = ψ(1) +O(2).The torque balance of the filament readsez ·∫ 10[x(s, t)− x(0, t)]× fvis ds−Kψ = 0, (3.42)where K = C/L3ξ⊥ω is the dimensionless spring constant. At leading order, we have13ψ˙(1) +12cos t+Kψ(1) = 0, (3.43)with ψ(1)(0) = 0. The long time solution is given byψ(1)(t) = − 32 + 18K2(3K cos t+ sin t) . (3.44)Finally, we obtain a simple analytical expression for the dimensionaless propulsive forceFp = 2F(2)p +O(4):F (2)p =γ − 1γ9K4(1 + 9K2)(3.45)We plot the variation of the propulsive force generated by a torsional spring arrangement asa function of the dimensionless spring constant in Fig. 3.7. As a limiting case of the two-segmentfilament, the torsional spring arrangement represents a highly localized flexibility focused at theactuation end with the rest of the filament being rigid, thus a straight line shape is maintainedduring the actuation with the tangent angle ψ(t) adjusted by the torque balance. In otherwords, torsional spring is the optimal elastic-rigid arrangement. The leading order velocity of alocal material point along the rigid filament is in the ±ey direction. If we assume a local elementds at s has a velocity u = uey, the resulting propulsive force density at this point scales withthe geometry as fprop ∼ cosψ sinψ [3], the extremum of which is obtained for ψ = ±45◦. Thisphysical scaling means that a filament can achieve higher propulsion by spending more timearound the optimal tangent angle in one period (i.e.,〈|ψ(1)|〉 should be close to 45◦). For thetorsional spring arrangement,〈|ψ(1)|〉 = 9K/(pi(1 + 9K2)) has a maximum of 3/(2pi) (≈ 27◦)which is obtained exactly at K = 1/3 where the optimal propulsion is also achieved. As Kdeviates from this optimal value,〈|ψ(1)|〉 decreases so that the propulsion becomes suboptimal.373.5. Continuous optimization0 2 4 6 8 10K00.050.100.150.20Fp2Figure 3.7: Propulsive thrust generated by a torsional spring arrangement at the actuation endconnected to a rigid filament as a function of the spring constant. The optimum propulsiveforce is F(2)p = 0.1875.Up to this section, via a discrete approach we have shown that both elastic-rigid and rigid-elastic cannot outperform the optimal uniform case but an elastic-elastic distribution can ac-tually enhance the propulsion as compared to the uniform filament. This begs the question:If we allow the stiffness to vary continuously along the filament, what is the optimum stiffnessprofile?3.5 Continuous optimizationWith the cases of two-segment and localized flexibility arrangements understood, we nowallow the bending stiffness to vary continuously along the entire filament. Using Eqs. (3.18)and (3.19), the asymptotic analysis shows that at order we haveSp40∂tψ(1) + ∂3s(A∗∂sψ(1))= 0, (3.46)where Sp0 and A∗ > 0 are defined in Sec. 3.2. By writing ψ(1) = Re{eith(s)}, we solve theresulting ODE for h(s) using Matlab’s build-in bvp4c solver for a given stiffness profile A∗. Asa validation of the algorithm, we solved h(s) using a regularized tanh stiffness profile to approx-imate an elastic-elastic filament as discussed in Sec. 3.3.1, and the numerical results matchedvery well with the analytical solutions. Once the evolution of the filament shape is obtained,the propulsive force can be evaluated from Eq. (3.26) where Gauss-Legendre quadrature is usedto calculate the numerical integration in space and time.Now, we want to understand which continuous distribution of stiffness optimizes the propul-sive performance of a driven filament. Following the procedure obtained by Moore [55] in thehigh Reynolds number case, we perform numerical optimization of the propulsive force over a383.5. Continuous optimizationlinear, quadratic and cubic distribution of bending stiffness, namely A∗ = as3 + bs2 + cs + 1.As the degree of the polynomial distribution increases, one would expect a better approxima-tion towards an arbitrary continuous stiffness profile. The optimization is performed using thefmincon solver in Matlab in the parameter space (Sp0, a, b, c) with the constraint of positivity:Sp0 > 0, A∗(s) > 0 for any s ∈ [0, 1].We introduce a theorem [58] of positivity of cubic polynomials on a given interval whichstates that the polynomial as3 + bs2 + cs + 1 is nonnegative for all s ∈ [0, 1] if and only if(m,n, p, q) ∈M ∪N , whereM = {(m,n, p, q) : m ≥ 0, n ≥ 0, p ≥ 0, q ≥ 0} ,N = {(m,n, p, q) : m ≥ 0, q ≥ 0, 4mp3 + 4qn3 + 27m2q2 − 18mnpq − n2p2 ≥ 0} .and m = a+ b+ c+ 1, n = b+ 2c+ 3, p = c+ 3, q = 1. Now, we have the constraints expressedonly in terms of the optimization parameters (a, b, c).0 0.5 1s0.00.20.40.60.81.0A∗0 0.5 1s0.00.20.40.60.81.0A∗Figure 3.8: (a) Optimal linear stiffness distribution, Sp0 ≈ 1.77. (b) Optimal quadratic stiffnessdistribution denoted by the dashed line (Sp0 ≈ 1.51, F (2)p ≈ 0.2224)Fig. 3.8 shows the optimal linear and quadratic stiffness distributions. Both cases reveal aless flexible material at the actuation end (s = 0). For example, the optimal linear profile has abasal sperm number Sp0 ≈ 1.77 and the scaled stiffness A∗ decreases as s increases. An optimallinear profile improves the propulsive performance by 3% compared with the maximum of theuniform case. The optimal quadratic profile improves the performance by 13% compared withthe optimal uniform case. The optimal cubic distribution (not shown in Fig. 3.8)only slightlyimproves upon the optimal quadratic case.We note that from a numerical optimization over an exponentially increasing profile, i.e.,A∗ = exp(as), a ≥ 0, the result would be that of the optimal uniform stiffness with a ≈ 0,Sp0 ≈1.89. This reduction to a uniform distribution indicates that if you start with an optimaluniform stiffness, increasing it down along the body would diminish propulsive performance.On the other hand, the optimal exponentially decreasing stiffness profile (Sp0 ≈ 0.1, a ≈ −20)improves the propulsion by 13.5% as compared to the maximum of a filament with uniform393.6. Effect of boundary conditionsstiffness, which is slightly better than the optimal quadratic distribution.3.6 Effect of boundary conditionsWe have shown that, different from the high Reynolds number case, the torsional springdoes not outperform the maximum propulsion achievable by a filament with uniform stiffnessunder a displacement-driven actuation at the cantilevered extremity of the filament. In general,we observe a trend where a decreasing in stiffness can potentially enhance the propulsive per-formance of a filament. This result seems to be consistent with the experimental observationsobtained by Maier et. al. [48] where they found that an exponentially decaying stiffness profilecan outperform a filament of uniform stiffness under a rotational actuation. In this section,by considering a torque-free displacement actuation and an angle-driven actuation, we shallshow that boundary conditions qualitatively affect the optimality of flexibility. As a result, onecannot simply extend the results by Maier et. al. [48] to other situations of propulsion wherethe actuation techniques may differ.3.6.1 Torque-free displacement oscillationWe consider a boundary driven passive filament where one end is under torque-free harmonicoscillation while the other end is free [51]. At the actuation end (s = 0), we havey(0, t) = y0 sinωt, Text(0) = −[Aκ]s=0 = 0, (3.47)After non-dimensionalization, we have y(0, t) = sin t and ψs(0, t) = 0 where = y0/L is theoscillation amplitude.Following the cantilevered case, we investigate an elastic-elastic filament at small oscillationamplitude. The governing equations for ψ(1)1 and ψ(1)2 are the same as those obtained for acantilevered filament, but with different boundary conditions. Once the deforming shape isobtained, the leading order propulsive force is then given byF (2)p =1− γγSp41〈ψ(1)1 (0, t)∂2sψ(1)1 (0, t) +12[(∂sψ(1)1 (α, t))2 − β(∂sψ(1)2 (α, t))2]〉. (3.48)For the classical case of uniform bending stiffness along the filament, the variation of propulsiveforce as a function of the sperm number is shown in Fig. 3.9(a) [4]. For the two values ofstiffness chosen according to the blue dots in Fig. 3.9(a), varying the relative proportion ofthe two segments leads to a non-monotonic variation of the propulsive force (Fig. 3.9(b)). Thelimits α = 0, 1 reduce to the case of uniform stiffness with a propulsive force F(2)p = 0.09 asindicated by the horizontal dash-dotted line in Fig. 3.9(b). By putting a more flexible materialat the actuation end i.e. β > 1 (e.g. β = 4.17, dashed line, Fig. 3.9(b)), the propulsiveforce generated by such an arrangement can be greater than the maximum possible propulsive403.6. Effect of boundary conditionsthrust (indicated by the horizontal dotted line in Fig. 3.9(b)) achievable by a filament with anyuniform bending stiffness. We remark that the proper choice of α is essential because one canobtain a lower propulsion by putting a more flexible material at the driven end. As an example,for β = 4.17, α & 0.55, the propulsion achieved using a two-segment filament is worse thanthe original uniform profile. The aforementioned enhancement in propulsion indicates that thefilament prefers to have a more flexible part at the actuation end with an optimal choice of α.It is interesting to note that these behaviors are qualitatively different from those observedfor a displacement-driven cantilevered filament (recall that more flexible material at the drivenend is suboptimal in the cantilevered case).0 2 4 6 8 10Sp00.030.060.090.12F(2)p(a)0 0.2 0.4 0.6 0.8 1α0.080.090.100.110.12F(2)p(b)β = 0.24β = 4.17Figure 3.9: (a) Propulsive thrust generated by a filament of uniform bending stiffness actuatedat one end [4]. The two blue dots indicate two filaments with different bending stiffnessesgenerating the same propulsive force. Propulsive thrust generated by two-segment filaments asa function of α at different values of stiffness ratio β. For β = 4.17 (A2 > A1, more flexiblematerials at the actuation end), the maximum propulsive force generated is greater than themaximum achievable thrust of a filament with any uniform stiffness.3.6.2 Angle oscillationFor an angle oscillation, we may write the dimensionless boundary actuation as ψ(0, t) = cos t, x(0, t) = 0 [51, 57]. The method of solutions are similar to those presented in Sec. 3.3with different boundary conditions which are detailed in Appendix C. For a filament withuniform stiffness, the optimal propulsive performance can be achieved around Sp ≈ 1.88, seeFig. 3.10(a). Similar to the displacement actuation of a cantilevered filament, an elastic-rigidfilament under angle actuation cannot outperform the optimum propulsion achievable by a fila-ment with uniform stiffness. In fact, by doing a two parameter (Sp0, a) numerical optimizationover an exponentially increasing stiffness profile, namely A∗ = exp(as), a ≥ 0, we found thatthe optimum converges to the maximum of a uniform stiffness with Sp0 ≈ 1.88, a ≈ 0. Thisindicates that if you starts with an optimum uniform stiffness, increasing the stiffness along thebody diminishes your performance.413.6. Effect of boundary conditions0 2 4 6 8 10Sp00.0070.0140.0210.0280.035Fp2(a)0 0.2 0.4 0.6 0.8 1α00.040.080.120.16Fp2(b)Sp2 = 1Sp2 = 2Sp2 = 4Sp2 = 10Figure 3.10: (a) Propulsive thrust generated by a filament of uniform bending stiffness withone end under angle oscillation [57]. (b) Propulsive thrust generated by a rigid-elastic filament(A1 =∞) as a function of α for different sperm numbers (Sp2) of the flexible part.On the other hand, the rigid-elastic arrangement can significantly improve the propulsioncompared with the uniform case as shown in Fig. 3.10(b) where the leading order propulsiveforce is plotted as a function of the relative proportion of the rigid segment for several spermnumbers of the flexible part. Take Sp2 = 10 as an example, the optimum propulsion thisarrangement can achieve is one magnitude larger than the maximum of the uniform case. Fora rigid-elastic filament, the flexible part is effectively driven by a displacement (rather than anangle oscillation) of the connecting point (s = α) though the rigid segment does not contributeto the propulsion. As the length of the flexible part, which contributes to propulsion, decreases(α increases), the displacement of the connecting point increases which tends to enhance theperformance. A competition between these two factors yields an optimal α for a given spermnumber of the flexible part.In summary, we’ve shown that though both elastic-rigid and rigid-elastic cannot outperformthe maximum of a uniform case for a displacement-driven cantilevered filament, the rigid-elastic423.7. Conclusionactually can enhance the propulsive performance for an angle-driven filament.3.7 ConclusionIn this Chapter, we have presented an analytical and numerical treatment of the optimalpropulsive performance of a boundary-driven passive cantilevered filament with variable stiffnessalong the body in the small amplitude oscillation limit. This work differs from and improvesupon previous studies by allowing a varying stiffness profile along the filament. In particular, forthe displacement-driven cantilevered filament, we show that different from the high Reynoldsnumber case, the torsional spring arrangement does not enhance the propulsion as compared toan optimal filament with uniform stiffness. However, an exponentially decreasing stiffness profilecan improve the propulsion. The rich behavior, particularly the enhancement of propulsion,revealed by our study may be useful for the development of efficient synthetic micro-propellersthat can swim as fast as, if not outperform microorganisms nature can offer.To conclude, we have shown that boundary conditions can qualitatively modify the optimalflexibility arrangement of a driven filament. Therefore, the optimal flexibility distribution mayvary from case to case and one may not extend the results obtained from one case to otherproblems of propulsion where the method of actuation might differ.43Chapter 4Nonlinear dynamics of a drivenmicrofilament in a viscous fluid4.1 IntroductionIn Chapter 3, we studied the dynamics and propulsion of a boundary driven microfilament ina viscous fluid in the small oscillation amplitude regime ( 1). For small amplitude, the fullynonlinear governing equations are linearized thus makes the dynamics analytically tractable.In this Chapter, we explore the fully nonlinear dynamics by solving the governing equationsnumerically for amplitudes that are not necessarily small.4.2 Mathematical formulationThe governing equation for x(s, t) as derived in Chapter 3 can be given byξ⊥xt = −(Axss)ss + (γ − 1)xs (2Asxss · xss + 3Axsss · xss) + γσsxs + σxss. (4.1)Differentiating Eq. (4.1) with respect to s and taking a dot product with xs, we obtain anauxiliary PDE for the Lagrange multiplier σ(s, t):γσss + (xs · xsss)σ = xs · (Axss)sss − (γ − 1)(2Asxss · xss + 3Axsss · xss)s. (4.2)In order to simplify the above Equation, we successively differentiate the identity xs · xs = 1and obtainxs · xss = 0, (4.3)xsss · xs + xss · xss = 0, (4.4)xssss · xs + 3xsss · xss = 0, (4.5)xsssss · xs + 4xssss · xss + 3xsss · xsss = 0. (4.6)444.3. Numerical implementationUsing the above identities, we obtainγσss − (xss · xss)σ =− (1 + 2γ)Assxss · xss − (2 + 7γ)Asxss · xsss− 3γAxsss · xsss − (1 + 3γ)Axss · xssss. (4.7)We non-dimensionalize the governing equations with respect to a time scale ω−1, lengthscale L and a force scale L2ξ⊥ω. The resulting dimensionless equations are given bySp4xt =− (A∗xss)ss + (γ − 1)xs (2A∗sxss · xss −A∗xs · xssss)+ Sp4γσsxs + Sp4σxss, (4.8)Sp4 [γσss − (xss · xss)σ] =− (1 + 2γ)A∗ssxss · xss − (2 + 7γ)A∗sxss · xsss− 3γA∗xsss · xsss − (1 + 3γ)A∗xss · xssss, (4.9)where the same symbols as the dimensional ones are used and A∗(s) = A(s)/A0 and Sp =L(ξ⊥ω/A0)1/4 is the sperm number evaluated using the stiffness at the basal end s = 0.The dimensionless boundary conditions at the free end are given byσxs − Sp−4(A∗xss)s = 0, A∗xss = 0 (4.10)The boundary conditions at the driven end s = 0 depends on the mechanisms of actuation. Fora torque free displacement oscillation, we may writex(0, t) = 0, y(0, t) = sin t, xss(0, t) = 0, (4.11)where = y0/L is the oscillation amplitude which is not necessarily small. For an angleoscillation, we can express the boundary conditions asψ(0, t) = cos t, x(0, t) = 0. (4.12)We note that for A∗(s) ≡ 1, the governing equations reduce to the case of a filament withuniform stiffness as given by previous work [59].4.3 Numerical implementationTo solve the equations of motion numerically, we implement a finite difference formulationthat combines two different methods suggested by Tornberg and Shelley [60] and Montenegro-Johnson [61]. In our approach, we treat Eq. (4.8) with a multi-step finite difference in time. In454.4. Uniform stiffnessorder to derive the numerical method, we write Eq. (4.8) symbolically asxt = F (xs,xss,xsss,xssss;σ, σs) (4.13)Now, the finite difference equation is given by3xn+1 − 4xn + xn−12∆t= Fn+1, (4.14)where ∆t is the time step, and tn = n∆t. Note that both linear and nonlinear terms in Fn+1are treated implicitly. In order to solve the nonlinear finite difference equation at tn+1, we usean iterative approach where at each iteration the implicit terms are linearized as−Sp−4(A∗xss)ss + Sp−4(γ − 1)x˜s (2A∗sx˜ss · xss −A∗x˜s · xssss) + γσsx˜s + σx˜ss, (4.15)where variables with tildes indicate values taken from the previous iteration. At the firstiteration of each time step, an extrapolation of the values from previous two time steps (2xn−xn−1) are used as initial guess. Similarly for the auxiliary equation for tension at tn+1, we haveSp4 [γσss − (x˜ss · x˜ss)σ] =− (1 + 2γ)A∗ssx˜ss · xss − (2 + 7γ)A∗sx˜ss · xsss− 3γA∗x˜sss · xsss − (1 + 3γ)A∗x˜ss · xssss+ λ(1− x˜s · xs), (4.16)where the term λ(1−x˜s ·xs) is introduced to penalize any length error that might arise from thenumerical approach with λ being the penalty factor. The filament is uniformly discretized intoN segments (h = 1/N) with the nodes (starting from s = 0) denoted as sj = jh, j = 0, 1, ..., N .Second order accurate centered finite differences as given by Tornberg and Shelley [60] are usedfor spatial derivatives in s.At each iteration in a time step, we solve a linear system of 3(N + 1) equations to obtainthe discrete values of the shape x(sj , tn+1) and Lagrange multiplier σ(sj , tn+1) simultaneously.In general, we terminate the iteration when the relative error of x in two successive iterationsis within 0.5%. To implement the boundary conditions, one-sided finite differences are used.4.4 Uniform stiffnessFor the simple case of a filament with uniform stiffness A∗ ≡ 1, all the terms that havederivatives of A∗(s) vanish. In Fig. 4.1, we plot the shapes of the deforming filament at differ-ent times in one period obtained from linear theory (as given in Chapter 3) along with numericalsimulations from the fully nonlinear theory for amplitude = 0.1 under a torque free displace-ment oscillation. As can be seen from Fig. 4.1, the numerical solution matches very well withthose from the small amplitude asymptotic expansion.464.4. Uniform stiffness0 0.2 0.4 0.6 0.8 1x−0.10−0.050.000.050.10yFigure 4.1: Shape (not to scale) of a deforming filament with uniform stiffness Sp = 4 at = 0.1 at different times npi/4, where n = 1, 2, ..., 8, the intensity of color decreases as timeincreases. The numerical solution is calculated with N = 160 and ∆t = 5 × 10−4. The reddashed lines denote the shapes obtained from the small amplitude linear theory [4]. The shapefrom numerical simulation matches very well with those from the small amplitude asymptoticexpansion.In Fig. 4.2, we show the shapes of a deforming filament at different times over one period ofthe actuation for = 1 and Sp = 4. It is interesting to note that, different from the linear theorywhere the waveforms are always symmetric (see Fig. 4.1), complex beating patterns emerge fromthe intrinsically symmetric model for large amplitude. These asymmetries in waveforms are dueto buckling instability followed by complex shape perturbations [59, 62].From our preliminary results, we’ve shown that the nonlinear model predicts that the com-pression of the filament due to the internal forces leads to a transition to symmetry breaking,as well as the breakdown of the small amplitude linear theory. In future study, we want tolook at the effect of non-uniform stiffness on the beating pattern for large amplitude oscillationand also the swimming characteristics by allowing the filament to swim freely under a properactuation.474.4. Uniform stiffness0 0.2 0.4 0.6 0.8 1x−1.0−0.50.00.51.0yt = 0t = pi/4t = 2pi/4t = 3pi/4t = 4pi/4t = 5pi/4t = 6pi/4t = 7pi/4t = 2piFigure 4.2: Shape of a deforming filament with uniform stiffness Sp = 4 at = 1 at differenttimes npi/4, where n = 0, 1, 2, ..., 8, the intensity of color decreases as time increases. Thesolution is calculated with N = 160 and ∆t = 5× 10−4. Pronounced buckling is observed.48Bibliography[1] Andrew A Biewener. Animal locomotion. Oxford University Press, 2003.[2] Lisa J Fauci and Robert Dillon. Biofluidmechanics of reproduction. Annu. Rev. FluidMech., 38:371–394, 2006.[3] Eric Lauga and Thomas R Powers. The hydrodynamics of swimming microorganisms. Rep.Prog. Phys., 72(9):096601, 2009.[4] On Shun Pak and Eric Lauga. Theoretical models in low-Reynolds-number locomotion.In Camille Duprat and Howard A. Stone, editors, Fluid-Structure Interactions in Low-Reynolds-Number Flows. Royal Society of Chemistry, 2014.[5] Edward M Purcell. Life at low Reynolds number. Am. J. Phys, 45(1):3–11, 1977.[6] Christopher Brennen and Howard Winet. Fluid mechanics of propulsion by cilia andflagella. Annu. Rev. Fluid Mech., 9(1):339–398, 1977.[7] Howard C Berg. E. coli in Motion. Springer Science & Business Media, 2008.[8] J Gray and GJ Hancock. The propulsion of sea-urchin spermatozoa. J. Exp. Biol., 32(4):802–814, 1955.[9] James Lighthill. Flagellar hydrodynamics. SIAM Rev., 18(2):161–230, 1976.[10] J Gray. Undulatory propulsion. Q. J. Microsc. Sci., 3(28):551–578, 1953.[11] Netta Cohen and Jordan H Boyle. Swimming at low Reynolds number: a beginners guideto undulatory locomotion. Contemp. Phys., 51(2):103–123, 2010.[12] AT Chwang and TY Wu. A note on the helical movement of micro-organisms. Proc. R.Soc. B, 178(1052):327–346, 1971.[13] Joseph B Keller and SI Rubinow. Swimming of flagellated microorganisms. Biophys. J.,16(2 Pt 1):151, 1976.[14] JJL Higdon. A hydrodynamic analysis of flagellar propulsion. J. Fluid Mech., 90(04):685–711, 1979.[15] J Gray. The mechanism of locomotion in snakes. J. Exp. Biol., 23(2):101–120, 1946.49Bibliography[16] ZV Guo and L Mahadevan. Limbless undulatory propulsion on land. Proc. Natl. Acad.Sci. USA, 105(9):3179–3184, 2008.[17] David L. Hu, Jasmine Nirody, Terri Scott, and Michael J. Shelley. The mechanics ofslithering locomotion. Proc. Natl. Acad. Sci. USA, 106(25):10081–10085, 2009. doi: 10.1073/pnas.0812533106.[18] Silas Alben. Optimizing snake locomotion in the plane. Proc. R. Soc. A, 469(2159), 2013.ISSN 1364-5021.[19] Werner Baumgartner, Florian Fidler, Agnes Weth, Martin Habbecke, Peter Jakob,Christoph Butenweg, and Wolfgang Bo¨hme. Investigating the locomotion of the sand-fish in desert sand using nmr-imaging. PloS One, 3(10):e3309, 2008.[20] Ryan D Maladen, Yang Ding, Chen Li, and Daniel I Goldman. Undulatory swimming insand: subsurface locomotion of the sandfish lizard. Science, 325(5938):314–318, 2009.[21] Yang Ding, Sarah S Sharpe, Andrew Masse, and Daniel I Goldman. Mechanics of undula-tory swimming in a frictional fluid. PLoS Comput. Biol., 8(12):e1002810, 2012.[22] Brian J Williams, Sandeep V Anand, Jagannathan Rajagopalan, and M Taher A Saif. Aself-propelled biohybrid swimmer at low Reynolds number. Nat. Commun., 5, 2014.[23] Ryan D Maladen, Yang Ding, Paul B Umbanhowar, and Daniel I Goldman. Undulatoryswimming in sand: experimental and simulation studies of a robotic sandfish. Int. J. Robot.Res., 30(7):793–805, 2011.[24] M. Sauzade, G. J. Elfring, and E. Lauga. Taylor’s swimming sheet: Analysis and improve-ment of the perturbation series. Phys. D, 240:1567 – 1573, 2011.[25] Tingnan Zhang and Daniel I. Goldman. The effectiveness of resistive force theory ingranular locomotiona). Phys. Fluids, 26(10):101308, 2014.[26] Daniel I. Goldman. Colloquium : Biophysical principles of undulatory self-propulsion ingranular media. Rev. Mod. Phys., 86:943–958, Jul 2014.[27] R Albert, MA Pfeifer, A-L Baraba´si, and P Schiffer. Slow drag in a granular medium.Phys. Rev. Lett., 82(1):205, 1999.[28] Glen Hill, Susan Yeung, and Stephan A Koehler. Scaling vertical drag forces in granularmedia. Europhys. Lett., 72(1):137, 2005.[29] Matthias Schro¨ter, Sibylle Na¨gle, Charles Radin, and Harry L Swinney. Phase transitionin a static granular system. Europhys. Lett., 78(4):44004, 2007.50Bibliography[30] Fuping Zhou, Suresh G. Advani, and Eric D. Wetzel. Simulation of slowly dragging acylinder through a confined pressurized bed of granular materials using the discrete elementmethod. Phys. Fluids, 19(1):013301, 2007.[31] A. Seguin, Y. Bertho, P. Gondret, and J. Crassous. Dense granular flow around a pene-trating object: Experiment and hydrodynamic model. Phys. Rev. Lett., 107:048001, Jul2011.[32] Saverio E Spagnolie and Eric Lauga. The optimal elastic flagellum. Phys. Fluids, 22(3):031901, 2010.[33] BM Friedrich, IH Riedel-Kruse, J Howard, and F Ju¨licher. High-precision tracking ofsperm swimming fine structure provides strong test of resistive force theory. J. Exp. Biol.,213(8):1226–1234, 2010.[34] Rafael D. Schulman, Matilda Backholm, William S. Ryu, and Kari Dalnoki-Veress. Dy-namic force patterns of an undulatory microswimmer. Phys. Rev. E, 89:050701, May 2014.[35] James Lighthill. Mathematical biofluiddynamics. Society for Industrial & Applied Mathe-matics, US, 1975.[36] Stephen Childress. Mechanics of swimming and flying, volume 2. Cambridge UniversityPress, 1981.[37] Leif E Becker, Stephan A Koehler, and Howard A Stone. On self-propulsion of micro-machines at low Reynolds number: Purcell’s three-link swimmer. J. Fluid Mech., 490:15–35, 2003.[38] Daniel Tam and A. E. Hosoi. Optimal stroke patterns for purcell’s three-link swimmer.Phys. Rev. Lett., 98:068105, Feb 2007.[39] O Pironneau and DF Katz. Optimal swimming of flagellated micro-organisms. J. FluidMech., 66(02):391–415, 1974.[40] Stephan Koehler, Tristan Spoor, and BS Tilley. Pitching, bobbing, and performance met-rics for undulating finite-length swimming filaments. Phys. Fluids, 24(9):091901, 2012.[41] RS Berman, O Kenneth, J Sznitman, and AM Leshansky. Undulatory locomotion of finitefilaments: lessons from caenorhabditis elegans. New J. Phys., 15(7):075022, 2013.[42] MJ Lighthill. Hydromechanics of aquatic animal propulsion. Annu. Rev. Fluid Mech., 1(1):413–446, 1969.[43] John Young, Simon M. Walker, Richard J. Bomphrey, Graham K. Taylor, and Adrian L. R.Thomas. Details of insect wing design and deformation enhance aerodynamic function andflight efficiency. Science, 325(5947):1549–1552, 2009. doi: 10.1126/science.1175928.51Bibliography[44] Theodore Yaotsu Wu. Fish swimming and bird/insect flight. Annu. Rev. Fluid Mech., 43(1):25–58, 2011.[45] Re´mi Dreyfus, Jean Baudry, Marcus L Roper, Marc Fermigier, Howard A Stone, andJe´roˆme Bibette. Microscopic artificial swimmers. Nature, 437(7060):862–865, 2005.[46] On Shun Pak, Wei Gao, Joseph Wang, and Eric Lauga. High-speed propulsion of flexiblenanowire motors: Theory and experiments. Soft Matter, 7:8169–8181, 2011.[47] Brian J Williams, Sandeep V Anand, Jagannathan Rajagopalan, and M Taher A Saif. Aself-propelled biohybrid swimmer at low reynolds number. Nat. Commun., 5, 2014.[48] Alexander M. Maier, Cornelius Weig, Peter Oswald, Erwin Frey, Peer Fischer, and TimLiedl. Magnetic propulsion of microswimmers with dna-based flagellar bundles. NanoLett., 16(2):906–910, 2016.[49] Chris H. Wiggins, D. Riveline, A. Ott, and Raymond E. Goldstein. Trapping and wiggling:Elastohydrodynamics of driven microfilaments. Biophys. J., 74(2):1043 – 1060, 1998. ISSN0006-3495.[50] Chris H. Wiggins and Raymond E. Goldstein. Flexive and propulsive dynamics of elasticaat low reynolds number. Phys. Rev. Lett., 80:3879–3882, Apr 1998.[51] Arthur A. Evans and Eric Lauga. Propulsion by passive filaments and active flagella nearboundaries. Phys. Rev. E, 82:041915, Oct 2010.[52] Hiroto Tanaka, John P Whitney, and Robert J Wood. Effect of flexural and torsional wingflexibility on lift generation in hoverfly flight. Integr. Comp. Biol., 51:142, 2011.[53] Jialei Song, Haoxiang Luo, and Tyson L Hedrick. Wing-pitching mechanism of hoveringruby-throated hummingbirds. Bioinspir. Biomim., 10(1):016007, 2015.[54] Kourosh Shoele and Qiang Zhu. Performance of a wing with nonuniform flexibility inhovering flight. Phys. Fluids, 25(4):041901, 2013.[55] M Nicholas J Moore. Torsional spring is the optimal flexibility arrangement for thrustproduction of a flapping wing. Phys. Fluids, 27(9):091701, 2015.[56] Se´bastien Camalet and Frank Ju¨licher. Generic aspects of axonemal beating. New J. Phys.,2(1):24, 2000.[57] Tony S. Yu, Eric Lauga, and A. E. Hosoi. Experimental investigations of elastic tailpropulsion at low reynolds number. Phys. Fluids, 18(9):091701, 2006.[58] Jochen W Schmidt and Walter Hess. Positivity of cubic polynomials on intervals andpositive spline interpolation. BIT Numerical Mathematics, 28(2):340–352, 1988.52[59] H Gadeˆlha, EA Gaffney, DJ Smith, and JC Kirkman-Brown. Nonlinear instability in flag-ellar dynamics: a novel modulation mechanism in sperm migration? J. R. Soc. Interface,page rsif20100136, 2010.[60] Anna-Karin Tornberg and Michael J Shelley. Simulating the dynamics and interactions offlexible fibers in stokes flows. Journal of Computational Physics, 196(1):8–40, 2004.[61] TD Montenegro-Johnson, Hermes Gadeˆlha, and David J Smith. Spermatozoa scatteringby a microchannel feature: an elastohydrodynamic model. R. Soc. open Sci., 2(3):140475,2015.[62] MK Jawed, NK Khouri, F Da, E Grinspun, and PM Reis. Propulsion and instability of aflexible helical rod rotating in a viscous fluid. Phys. Rev. Lett., 115(16):168101, 2015.53Appendix ANumerical implementationIn this appendix, we present the numerical methods implemented in the optimization ofinfinite filaments and the solution to the equations of motion of finite length filaments.A.1 OptimizationThe numerical optimization (see Sec. 2.3.1) is performed using MATLAB’s built-in fmin-search function, which implements the Nelder-Mead simplex algorithm. We truncate the Fourierseries by taking n∗ = 100 to have a sufficient spectral accuracy and use m = 1000 points forthe Gauss-Legendre integration scheme. Further increase in spectral and spatial resolutionhas a negligible effect on the optimization. The optimization search routine iterates until thealgorithm detects a local solution gradient with a relative error tolerance of 10−14. For eachiteration, the swimming speed U is obtained by solving the force balance in the swimmingdirection using MATLAB’s fzero function, which runs until a relative error of 10−16 is reached.A variety of shapes are provided as the initial guess for the starting of the optimization. Theoptimization calculation is iterated by taking the converged shape of the previous calculationas the initial guess until the shapes acquired in two successive calculations are consistent. Theoptimal shape obtained does not vary with the initial guess.We validate our approach by solving the optimal shape for the Newtonian case. For New-tonian swimming, the swimming speed U can be obtained by a simple matrix inversion due tolinearity. The optimal shape obtained from our numerical approach agrees with the analyticalsolution of Lighthill[35].A.2 Numerical solution for finite swimmersThe study of swimming characteristics requires solving the force and torque balance ofthe finite swimmer as formulated in Sec. 2.2. The force- and torque-free conditions posed inEq. (2.19) and (2.20) provide a system of non-linear ordinary differential equations (ODEs)for the swimmer’s linear and angular velocities in terms of its instantaneous location and ori-entation. The instantaneous velocities in turn, once obtained, can be integrated over time todetermine the trajectory, location and orientation of the swimmer.Having assumed that the centerline of the waveform is initially aligned with the x−axisof the lab frame, we solve the swimming problem numerically. Starting at t = 0 with a time54A.2. Numerical solution for finite swimmersstep ∆t, we denote ti = i∆t. With this notation, we employ a second order multi-step finitedifference method to discretize the ODEs such thatxi+1 =43xi − 13xi−1 +2∆t3(2x˙i − x˙i−1). (A.1)We do similarly for θi+1, and then Θi+1 can be computed. To initialize this numerical scheme,we need both [x0, θ0] and [x1, θ1]. At the first time step, [x1, θ1] is computed using the Runge-Kutta fourth order method. At each time step (i ≥ 1), we obtain [x˙i, θ˙i] by solving the integralequations for force and torque balance using Gauss-Legendre quadrature integration coupledwith MATLAB’s fsolve routine. We use m points along the filament for the Gauss-Legendreintegration method. The fsolve routine in Matlab Optimization Toolbox attempts to solvea system of equations by minimizing the sum of squares of all the components. We set thetermination tolerance on both the function value and independent variables to 10−14. Wegenerally use m = 1000 points along the filament and Tm = 500 time steps for one periodT of the motion. The number of Fourier modes is taken as n∗ = 100. Further increasing ofthe number of spatial points or time steps have no significant influence on the accuracy of theresults. All the numerical simulations are performed using MATLAB.55Appendix BSolution of an elastic-elastic filamentIn this section, we derive the solution of an elastic-elastic filament.B.1 Equations of motionThe dimensionless governing equations are given bySp41ψ1t = −∂4sψ1 + ∂s(ψ1sτ1) + γψ1s(ψ1s∂2sψ1 + τ1s),τ1ss − 1γψ21sτ1 = −∂s(ψ1s∂2sψ1)− 1γψ1s∂3sψ1,Sp41ψ2t = −β∂4sψ2 + ∂s(ψ2sτ2) + γψ2s(βψ2s∂2sψ2 + τ2s),τ2ss − 1γψ22sτ2 = −β∂s(ψ2s∂2sψ2)− β 1γψ2s∂3sψ2. (B.1)The force balance on the entire filament is given byFext(0, t) +∫ 10fvis ds = 0, (B.2)or−τ1t + ∂ssψ1n = −∫ 10fvis ds. (B.3)In other words, we haveτ1(0, t) = t(0, t) ·∫ 10fvis ds, (B.4)∂2sψ1(0, t) = −n(0, t) ·∫ 10fvis ds, (B.5)Since we prescribed the displacement of the point of actuation (s = 0), we may writex(0, t) = [0, sin t]T= x(0, 0) + Sp−41∫ t0dt[n(−∂3sψ1 + ψ1sτ1)+ γt (ψ1s∂2sψ1 + τ1s)]s=0 (B.6)56B.2. Propulsive forceDifferentiating with respect to time, we obtain0 = − sinψ (−∂3sψ1 + ψ1sτ1)+ γ cosψ (ψ1s∂2sψ1 + τ1s) ,Sp41 cos t = cosψ(−∂3sψ1 + ψ1sτ1)+ γ sinψ (ψ1s∂2sψ1 + τ1s) . (B.7)The torque-free condition reads ψ1s(0, t) = 0. At the free end, we have τ2(1, t) = 0, ∂ssψ2(1, t) =0, ψ2s(1, t) = 0. At the connecting point, we have ψ1(α, t) = ψ2(α, t), τ1(α, t) = τ2(α, t), ψ1s(α, t) =βψ2s(α, t), ∂2sψ1(α, t) = β∂2sψ2(α, t). The last boundary condition is given by the continuity ofthe viscous force at s = α.At zeroth order, we have ∂ssτ1,0 = 0, ∂ssτ2,0 = 0, which have trivial solutions τ1,0(s, t) =τ2,0(s, t) = 0. At the first order, we also have ∂ssτ1,1 = 0, ∂ssτ2,1 = 0 with zero solutions. Thisindicates that tension along the filament τ ∼ O(2), which is consistent with previous studiesof uniform stiffness.Now, expanding the equations for ψ1 and ψ2, we haveSp41∂tψ(1)1 + ∂4sψ(1)1 = 0, Sp41∂tψ(1)2 + β∂4sψ(1)2 = 0. (B.8)The corresponding B.C.s at s = 1 are ∂sψ(1)2 (1, t) = 0, ∂2sψ(1)2 (1, t) = 0. By expandingEq. (B.7), one can show that ∂3sψ(1)1 (0, t) = −Sp41 cos t. The other corresponding B.C.s aregiven by ∂sψ(1)1 (0, t) = 0, ψ(1)1 (α, t) = ψ(1)2 (α, t), ∂sψ(1)1 (α, t) = β∂sψ(1)2 (α, t), ∂2sψ(1)1 (α, t) =β∂2sψ(1)2 (α, t), ∂3sψ(1)1 (α, t) = β∂3sψ(1)2 (α, t).The PDEs for the local tangent can be reduced to two ODEs,Sp41h1(s) + ∂4sh1 = 0, Sp41h1(s) + β∂4sh1 = 0, (B.9)with the translated B.C.s given byh1s(0) = 0, h1sss(0) = −sp41, h2s(1) = 0, h2ss(1) = 0, h1(α) = h2(α),h1s(α) = βh2s(α), h1ss(α) = βh2ss(α), h1sss(α) = βh2sss(α). (B.10)These two ODEs fall into the type so-called hyper-diffusion equations. By assuming a solutionof the form h1 = ceks, one can obtain its solution.B.2 Propulsive forceTo leading order, the local tangent vectort1 ∼[0, ψ(1)1 (s, t)]T, t2 ∼[0, ψ(1)2 (s, t)]T. (B.11)57B.2. Propulsive forceBy differentiate x(s, t) with respect to time in Eq. 3.20, we obtain for the first segmentu1(s, t) ∼[0, cos t+ ∫ s0∂tψ(1)1 ds]T. (B.12)Similarly for the second segment, we haveu2(s, t) ∼[0, cos t+ ∫ α0∂tψ(1)1 ds+ ∫ sα∂tψ(1)2 ds]T. (B.13)Then we havefx1 = ex · fvis,1 ∼ −2Sp41(1γ− 1)ψ(1)1(∫ s0∂tψ(1)1 ds+ cos t). (B.14)Noting that ∂tψ(1)1 = −Sp−41 ∂4sψ(1)1 , we can calculate the integral and obtainfx1 =(1γ− 1)2ψ(1)1 (s, t)∂3sψ(1)1 (s, t), (B.15)where the boundary term ∂3sψ(1)1 (0, t) = −Sp41 cos t are used.Integrating fx1 from 0 to α, we haveFx1 =∫ α0fx1 ds ∼(1γ− 1)2(−12∂s(ψ(1)1 )2(α, t)+ ψ(1)1 (α, t)∂2sψ(1)1 (α, t)− ψ(1)1 (0, t)∂2sψ(1)1 (0, t))(B.16)We do similarly for Fx2, and the propulsive force at a given time t is then Fx ∼ Fx1 + Fx2.Taking the time average of Fx over a period 2pi, we obtain the second order propulsive force:F (2)p =1− γγSp41〈ψ(1)1 (0, t)∂2sψ(1)1 (0, t) +12[(∂sψ(1)1 (α, t))2 − β(∂sψ(1)2 (α, t))2]〉, (B.17)where 〈...〉 = (∫ 2pi0 (...) dt)/2pi is the average over one period. We note that terms equate to zeroafter integration are neglected along the way.58Appendix CBoundary conditions for angleoscillationThe method of solutions for a two-segment filament under angle oscillation is the same asthose in Appendix B, but with different boundary conditions.For an elastic-elastic filament, the boundary conditions at the connecting point s = α andthe free end s = 1 are the same as those in Appendix B. However, the boundary conditions atthe actuation end is given byψ(1)1 (0, t) = cos t, ∂3sψ(1)1 (0, t) = 0. (C.1)For a rigid-elastic filament, we have ψ1 = ψ1(t). The two boundary conditions at the freeend are the same as the elastic-elastic case. The boundary conditions at the connecting pointare given byψ(1)2 (α, t) = cos t, ∂3sψ(1)2 (α, t) = α∂4sψ(1)2 (α, t). (C.2)59
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Inertialess swimming and propulsion of slender bodies
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Inertialess swimming and propulsion of slender bodies Peng, Zhiwei 2016
pdf
Page Metadata
Item Metadata
Title | Inertialess swimming and propulsion of slender bodies |
Creator |
Peng, Zhiwei |
Publisher | University of British Columbia |
Date Issued | 2016 |
Description | In this thesis, two problems relevant to the biological locomotion in inertialess environments are studied, one is the characteristics of undulatory locomotion in granular media, the other is the optimal flexibility of a driven microfilament in a viscous fluid. Undulatory locomotion is ubiquitous in nature and observed in different media, from the swimming of flagellated microorganisms in biological fluids, to the slithering of snakes on land, or the locomotion of sandfish lizards in sand. Despite the similarity in the undulating pattern, the swimming characteristics depend on the rheological properties of different media. Analysis of locomotion in granular materials is relatively less developed compared with fluids partially due to a lack of validated force models but recently a resistive force theory in granular media has been proposed and shown useful in studying the locomotion of a sand-swimming lizard. In this work, we employ the proposed model to investigate the swimming characteristics of a slender filament, of both finite and infinite length, undulating in a granular medium and compare the results with swimming in viscous fluids. In particular, we characterize the effects of drifting and pitching in terms of propulsion speed and efficiency for a finite sinusoidal swimmer. We also find that, similar to Lighthill's results using resistive force theory in viscous fluids, the sawtooth swimmer is the optimal waveform for propulsion speed at a given power consumption in granular media. Though it is understood that flexibility can improve the propulsive performance of a filament in a viscous fluid, the flexibility distribution that generates optimal propulsion remains largely unexplored. In this work, we employ the resistive force theory combined with the Euler-Bernoulli beam model to examine the optimal flexibility of a boundary driven filament in the small oscillation amplitude limit. We show that the optimality qualitatively depends on the boundary actuation. For large amplitude actuation, our numerics show that complex asymmetry in the waveforms emerge. The results complement our understanding of inertialess locomotion and provide insights into the effective design of locomotive systems in various environments. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2016-04-21 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0300048 |
URI | http://hdl.handle.net/2429/57718 |
Degree |
Master of Applied Science - MASc |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2016-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2016_may_peng_zhiwei.pdf [ 924.98kB ]
- Metadata
- JSON: 24-1.0300048.json
- JSON-LD: 24-1.0300048-ld.json
- RDF/XML (Pretty): 24-1.0300048-rdf.xml
- RDF/JSON: 24-1.0300048-rdf.json
- Turtle: 24-1.0300048-turtle.txt
- N-Triples: 24-1.0300048-rdf-ntriples.txt
- Original Record: 24-1.0300048-source.json
- Full Text
- 24-1.0300048-fulltext.txt
- Citation
- 24-1.0300048.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0300048/manifest