UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Swimming in slime Pachmann, Sydney 2008

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
[if-you-see-this-DO-NOT-CLICK]
[if-you-see-this-DO-NOT-CLICK]
ubc_2008_fall_pachmann_sydney.pdf [ 796.37kB ]
Metadata
JSON: 1.0066549.json
JSON-LD: 1.0066549+ld.json
RDF/XML (Pretty): 1.0066549.xml
RDF/JSON: 1.0066549+rdf.json
Turtle: 1.0066549+rdf-turtle.txt
N-Triples: 1.0066549+rdf-ntriples.txt
Original Record: 1.0066549 +original-record.json
Full Text
1.0066549.txt
Citation
1.0066549.ris

Full Text

Swimming in SlimebySydney PachmannB.Sc. First Class Honours, The University of Calgary, 2006A THESIS SUBMITTED IN PARTIAL FULFILMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMaster of ScienceinThe Faculty of Graduate Studies(Mathematics)The University Of British Columbia(Vancouver)August, 2008c Sydney Pachmann 2008AbstractThe purpose of this thesis is to study the problem of a low Reynolds numberswimmer that is in very close proximity to a wall or solid boundary in a non-Newtonian fluid. We assume that it moves by propagating waves down its lengthin one direction, creating a thrust and therefore propelling it in the oppositedirection. We model the swimmer as an infinite, inextensible waving sheet.We consider two main cases of this swimming sheet problem. In the firstcase, the type of wave being propagated down the length of the swimmer isspecified. We compare the swimming speeds of viscoelastic shear thinning,shear thickening and Newtonian fluids for a fixed propagating wave speed. Wethen compare the swimming speeds of these same fluids for a fixed rate of workper wavelength. In the latter situation, we find that a shear thinning fluidalways yields the fastest swimming speed regardless of the amplitude of thepropagating waves. We conclude that a shear thinning fluid is optimal for theswimmer. Analytical results are obtained for various limiting cases. Next, weconsider the problem with a Bingham fluid. Yield surfaces and flow profiles areobtained.In the second case, the forcing along the length of the swimmer is specified,but the shape of the swimmer is unknown. First, we solve this problem for aNewtonian fluid. Large amplitude forcing yields a swimmer shape that has aplateau region following by a large spike region. It is found that there existsan optimal forcing that will yield a maximum swimming speed. Next, we solvethe problem for moderate forcing amplitudes for viscoelastic shear thickeningand shear thinning fluids. For a given forcing, it is found that a shear thinningfluid yields the fastest swimming speed when compared to a shear thickeningfluid and a Newtonian fluid. The difference in swimming speeds decreases asthe bending stiffness of the swimmer increases.iiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2 Overview of Non-Newtonian Fluids . . . . . . . . . . . . . . . . 51.2.1 Constitutive Models . . . . . . . . . . . . . . . . . . . . . 81.3 Our Problem: Main Assumptions . . . . . . . . . . . . . . . . . 111.3.1 Zero Force on Swimmer Condition . . . . . . . . . . . . . 151.4 Preliminary Problem . . . . . . . . . . . . . . . . . . . . . . . . 171.4.1 Mathematical Formulation . . . . . . . . . . . . . . . . . 171.4.2 Stream Function and Asymptotics . . . . . . . . . . . . . 191.4.3 Leading Order Equations . . . . . . . . . . . . . . . . . . 201.4.4 Second Order Equations . . . . . . . . . . . . . . . . . . 211.4.5 Limiting Cases . . . . . . . . . . . . . . . . . . . . . . . . 231.5 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . 252 Swimming with Lubrication Theory . . . . . . . . . . . . . . . . 272.1 Mathematical Formulation . . . . . . . . . . . . . . . . . . . . . 272.2 Newtonian Problem . . . . . . . . . . . . . . . . . . . . . . . . . 312.3 Viscoelastic Problem . . . . . . . . . . . . . . . . . . . . . . . . 332.3.1 Effective Viscosity . . . . . . . . . . . . . . . . . . . . . . 34iiiTable of Contents2.3.2 Numerical Formulation . . . . . . . . . . . . . . . . . . . 362.3.3 Fixed Dimensional Work . . . . . . . . . . . . . . . . . . 392.3.4 Limiting Cases . . . . . . . . . . . . . . . . . . . . . . . . 412.4 Bingham Problem . . . . . . . . . . . . . . . . . . . . . . . . . . 432.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 533 Elastic Swimming with Lubrication Theory . . . . . . . . . . . 553.1 Mathematical Formulation . . . . . . . . . . . . . . . . . . . . . 563.2 Newtonian Problem . . . . . . . . . . . . . . . . . . . . . . . . . 603.2.1 Numerical Formulation . . . . . . . . . . . . . . . . . . . 613.2.2 Low Amplitude Solution . . . . . . . . . . . . . . . . . . 643.3 Viscoelastic Results . . . . . . . . . . . . . . . . . . . . . . . . . 663.3.1 Low Amplitude Solution . . . . . . . . . . . . . . . . . . 693.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 704 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.1 Summary of Results . . . . . . . . . . . . . . . . . . . . . . . . . 744.2 Further Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77AppendicesA Full Analysis of Parent Problem . . . . . . . . . . . . . . . . . . 81A.1 Mathematical Formulation . . . . . . . . . . . . . . . . . . . . . 81A.2 Stream Function and Asymptotics . . . . . . . . . . . . . . . . . 83A.3 Leading Order Equations . . . . . . . . . . . . . . . . . . . . . . 85A.4 Second Order Equations . . . . . . . . . . . . . . . . . . . . . . 89A.5 Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93A.6 Limiting Cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94A.6.1 Distance Between Wall and Swimmer Approaches Infinity 94A.6.2 Swimmer Approaches the Wall . . . . . . . . . . . . . . . 95B Details of Bingham Analysis . . . . . . . . . . . . . . . . . . . . . 97B.1 Case I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97B.2 Case II . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98B.3 Case III . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99B.4 Case IV . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100ivTable of ContentsB.5 Case V . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101C Code For Chapter Two . . . . . . . . . . . . . . . . . . . . . . . . 102C.1 Viscoelastic Problem Code . . . . . . . . . . . . . . . . . . . . . 102C.2 Bingham Problem Code . . . . . . . . . . . . . . . . . . . . . . . 108D Code For Chapter Three . . . . . . . . . . . . . . . . . . . . . . . 112D.1 Newtonian Problem Code . . . . . . . . . . . . . . . . . . . . . . 112D.2 Viscoelastic Problem Code . . . . . . . . . . . . . . . . . . . . . 115vList of Tables2.1 Fluid types and the corresponding parameter values. . . . . . . . 352.2 The specific parameter values chosen for the various fluid types. . 362.3 Bingham Fluid Flow Profile Regions. In the outer regions, wesee parabolic flow profiles typically characterized by Newtonianflows. In the inner region, we see what we refer to as pseudo-plugflow. In this region the velocity profile is constant. . . . . . . . . 452.4 Location of Yield Surfaces. . . . . . . . . . . . . . . . . . . . . . 47viList of Figures1.1 Geometry of Problem . . . . . . . . . . . . . . . . . . . . . . . . 121.2 Geometry of Problem in Moving Frame of Reference . . . . . . . 142.1 Lubrication geometry . . . . . . . . . . . . . . . . . . . . . . . . 272.2 Effective Viscosities. Red: Newtonian, Solid Blue: Shear Thin-ning, Dashed Blue: Unrealistic Shear Thinning, Solid Green:Shear Thickening, Dashed Green: Unrealistic Shear Thickening . 352.3 Non-Dimensional Swimming Speeds. Red: Newtonian. Blue:Shearing Thinning. Green: Shear Thickening. . . . . . . . . . . . 392.4 Dimensional Swimming Speeds. Red: Newtonian. Blue: Shear-ing Thinning. Green: Shear Thickening. . . . . . . . . . . . . . . 412.5 Newtonian Flow Profiles for ah = 0.5. The speed of the swimmer,U, is depicted by the vertical red line. . . . . . . . . . . . . . . . 462.6 Pressure Gradients and Yield Surfaces for B = 1 and ah = 1/2.The Newtonian pressure gradient is blue, the Bingham red. Inthe lower figure, Y+ is blue, Y- is red, and the shaded greenrepresents the pseudo-plug. The surface of the swimmer is repre-sented by the black curve and the pink vertical line shows wherethe pressure gradient changes sign. . . . . . . . . . . . . . . . . . 502.7 Characteristic Bingham Velocity Profiles for B = 1 and ah = 1/2.Blue: Yielded Region. Red: Pseudo-plug Region. . . . . . . . . . 512.8 Pressure Gradients and Yield Surfaces for large and small B,ah = 1/2. In the left figures, the Newtonian pressure gradient isblue, the Bingham red. In the right figures, Y+ is blue, Y- is red,and the shaded green represents the pseudo-plugs. The surfaceof the swimmer is represented by the black curve and the pinkvertical line shows where the pressure gradient changes sign. . . . 522.9 Swimming speed and mass flux for varying Bingham number.The corresponding Newtonian values are shown in red. . . . . . . 53viiList of Figures2.10 Example of alternate wave shape and resulting swimming speeds.Red: Newtonian. Blue: Shearing Thinning. Green: Shear Thick-ening. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 543.1 Geometry of Lubrication Problem. This time, the swimmer shapein unknown as well as the swimming speed . . . . . . . . . . . . 563.2 Newtonian Solution of Elastic Lubrication Problem. In the bot-tom left hand corner of each figure we see the swimming speed(red) and the mass flux (green) versus the forcing amplitude.The magenta horizontal lines are the calculated plateau valuesfor each case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 633.3 Phase portraits for Newtonian Solution of Elastic LubricationProblem with large amplitude forcing. . . . . . . . . . . . . . . . 633.4 Solution for viscoelastic elastic lubrication problem for variousfluid types with D=0.2. Red: Newtonian. Blue: Shearing Thin-ning. Green: Shear Thickening. . . . . . . . . . . . . . . . . . . . 673.5 Solution for viscoelastic elastic lubrication problem for variousfluid types with D=7. Red: Newtonian. Blue: Shearing Thin-ning. Green: Shear Thickening. . . . . . . . . . . . . . . . . . . . 683.6 Solution for viscoelastic elastic lubrication problem for variousfluid types with D=27. Red: Newtonian. Blue: Shearing Thin-ning. Green: Shear Thickening. . . . . . . . . . . . . . . . . . . . 693.7 Comparison between analytical and numerical results for low am-plitude swimming problem. Numerical results are solid lines, an-alytical results are shown by *. Red: Newtonian. Blue: ShearingThinning. Green: Shear Thickening. . . . . . . . . . . . . . . . . 703.8 Optimal swimmer shape for D=0.2. . . . . . . . . . . . . . . . . . 713.9 Strain rates on the swimmer (cyan) and on the wall (black) fora=100 and D=1.5. Note that the largest strain rates occur in theregion where the plateau ends and the spike begins. . . . . . . . 73viiiAcknowledgementsI would like to thank my supervisors, Neil Balmforth and Daniel Coombs, forintroducing me to the field of fluid dynamics and allowing me to work on thisproject. I am grateful for all the help they have provided.I would like to thank the other mathematical faculty at UBC from whom Itook courses from. I feel I have learned a lot from you all. I would also like tothank another UBC graduate student, Samara Pillay, for her support and helpthroughout my time here at UBC. I would also like to make a special thanks toall the students in the IAM for just being great and making my experience heresomething truly fantastic. Finally, I express my sincere appreciation to all thehelpful staff of the Mathematics Department and the IAM.This research was supported by an NSERC PGS M.ixDedicationTo my parents, for always being there for me, my sister and brother, who tryvery hard find what I do interesting, and to Tyler, for putting up with methroughout.xChapter 1IntroductionYears ago, the question of how microscopic organisms propel themselves throughfluid was asked. The standard mechanism of propulsion for larger bodies wasunderstood to rely on the shedding of vortices in an inviscid fluid. It wassomething entirely new to consider an organism propelling itself through fluidwith a very small Reynolds number [37]. The mechanism of vortex shedding nolonger was valid, and so the field of low Reynolds number swimming was born.Typically, the microorganism under specific analysis in this field is a sperma-tozoon. This is due to the biological importance of spermatozoa locomotion aswell as the fact that the mechanisms of spermatozoa locomotion are significantlymore amenable to mathematical modelling than other microorganisms. Manyflagellated microorganisms are propelled by a rather complicated spiralling orcorkscrew motion of their flagella. However, the flagellar motion of spermatozoatends to be more planar. In certain situations it can also be low amplitude [21].Thus, it is clear why spermatozoa are desirable to model from a mathemati-cian?s perspective. Contributions to this field have been made by people such asTaylor, Reynolds, Gray and Hancock, with continual contributions from leadingresearchers to this day [37, 33, 14, 15, 22, 25].The purpose of this thesis is to study the problem of a low Reynolds numberswimmer that is in very close proximity to a wall or solid boundary in a non-Newtonian fluid. We assume that it moves by propagating waves down itslength in one direction, creating a thrust and therefore propelling it in theopposite direction. We model the swimmer as an infinite, inextensible wavingsheet, a standard approach in this field [37, 33, 22, 25]. Of course, realisticmicroorganisms have a body or head and flagellum of finite length and width,but this simplification to a swimming sheet provides useful insight into thegeneral problem of swimming. This formulation of the problem also proves tohave a wider range of applications [8, 26, 39, 2, 40, 18].We consider two main cases of the swimming sheet problem. In the firstproblem, we suppose that the type of wave being propagated down the length1Chapter 1. Introductionof the swimmer is specified. In the second problem, the forcing along the lengthof the swimmer is specified, but the shape of the swimmer is unknown.As described above, these questions are related to many other similar prob-lems that have previously been studied. We take a unique twist by consideringthe lubrication approximation, introducing a wall and a non-Newtonian fluidand eventually an elastic swimmer. In this chapter we will outline some previ-ous work done on similar problems, setting the backdrop for the current topic.As well, we will briefly introduce the reader to viscoelastic fluids in preparationfor the following discussions and mathematical analysis. We will then presentthe general problem under consideration, emphasizing the assumptions that aremade and used in all following chapters, and the physical rationalizations forthese assumptions. Finally, we will detail the problem of a two dimensionalswimmer in the presence of a wall without the lubrication approximation tomotivate the use of the thin gap theory.1.1 Literature ReviewIn the literature, there are several trains of thought about how to approach thelow Reynolds number swimming problem. In all considerations of the problem,it is the swimming speed of the organism that is of primary interest. One ap-proach is to simplify the problem to the two dimensional problem of a swimmingsheet, as we will be doing in our discussion of this topic. Typically, an asymp-totic approach is used where one assumes the form of the swimmer to be thatof a low amplitude wave. The solution is expanded as a Taylor series in powersof the non-dimensional wave amplitude [37, 33, 9, 22, 25]. A common result isthat the leading order swimming speed is proportional to the amplitude of theswimmer squared.G. I. Taylor, in [37], approached the swimming problem in this way. He con-sidered an infinite, inextensible, two dimensional waving sheet in a Newtonianfluid. The sheet was fixed to be a sinusoid, which he gave dimensionally byy0 = bsin(kx -sigmat), (1.1)where b is the amplitude, k the wave number, and sigma is the frequency. Nodynamical considerations were made for the sheet. He assumed that the waveswere low amplitude and found that at first order, the swimming speed waszero. Continuing the calculation to second order, he found a non-zero swimming2Chapter 1. Introductionspeed, which is henceforth referred to as the leading order swimming speed. Thisswimming speed, dimensionally, was calculated to beVc =12b2k2, (1.2)where V is the non-dimensional swimming speed, c is the speed of the propa-gating waves, b is the amplitude of the waves, and k is the wave number. Asis evident, this leading order swimming speed is proportional to the amplitudeof the swimmer squared. In this paper, it was observed that real microscopicorganisms are in contact with fluid on both sides. When adding a fluid to bothsides of the swimmer, only the rate of dissipation of energy changed, from Wto 2W. The speed of the swimmer did not. When referring to Taylor?s problembelow, we are referring to the problem of an infinite, inextensible waving sheetwith fluid on one side.Another approach to the problem, in three dimensions, is to model the swim-mer as an infinite cylinder. This method involves the use of what is commonlyreferred to as resistive force theory. Here, a distribution of stokeslets and dou-blets is placed along the center of the cylinder [15, 14, 7, 27]. Similar to resistiveforce theory is slender-body theory, which was introduced when incorporatingthe head of a microorganism into the analysis. It provides more accuracy in thisinstance than resistive force theory [16, 19]. In this swimming cylinder approach,the assumption of low amplitude waves is not necessarily made. However, whenthe assumption is made, the results agree qualitatively with those of the swim-ming sheet approach, that the swimming speed is proportional to the amplitudeof the swimmer squared.Of these above approaches, several include solid boundaries [33, 6, 22]. Katz,in [22], considered a swimmer in a channel that moves by propagating a pre-scribed wave down its length. This wave form is given by, letting (x0,y0) be apoint on the sheet surface and (xm,0) be the mean position of (x0,y0),x0 = xm + acos k[xm + (c -V )t] + dsin k[xm + (c -V )t],y0 = bsin k[xm + (c -V )t,] (1.3)where c is the speed of the propagating wave and V the swimming speed. Thiswave is slightly more general than that considered by G. I. Taylor, in [37]. In [22],Katz found that the leading order velocity of the organism was proportional tothe amplitude of the swimmer squared. Therefore, introducing this more general3Chapter 1. Introductionwave form did not change the result qualitatively from those that assumed nowalls. In the limit as one wall of the channel goes to infinity, Katz?s resultsreturn the results of Blake, in [6], who considered a swimmer in the presenceof one wall with a general wave form. One thing to note about this generalwave form is that if the horizontal wave motion is large enough, the sheet canactually be propelled in the opposite direction.In [22], Katz also considered the problem using lubrication theory. Lubri-cation theory is appropriate for the problem where the swimmer is very closeto the wall. It was found that the leading order swimming speed, V0, is largerthan in the case with the wall being further away. The dimensional result wasV0c =32 + parenleftbighb parenrightbig2, (1.4)where h is the gap width. This implies that swimming speed is bounded by thespeed of the propagating waves. For this lubrication problem, the velocities dueto horizontal motion of the sheet are of a higher order and do not contributeto the leading order solution. Therefore, the sheet cannot reverse directions inthis problem, as it could in the problem without the lubrication approximationapplied.The problem of swimmers in a channel has also been considered numerically.In [12], the authors look at the problem of one and multiple swimmers in aNewtonian fluid with walls that are both rigid and elastic. Here, the resultswere very comparable to real observations of sperm movement.The problem of a swimmer in a viscoelastic fluid has also been discussedin various contexts. Chaudhury, in [9], modified Taylor?s problem to include asecond-order viscoelastic fluid. He found that viscoelasticity augments propul-sion in low Reynolds numbers, but as Reynolds number increases, there is arange in which the viscoelasticity hampers it. More recently, a paper written byFu et al. [13] considers the swimming infinite cylinder problem using an upperconvected Maxwell model as the constitutive equation for the viscoelastic fluid.They found that viscoelasticity tends to decrease the swimming speed and thatswimming speed can reverse direction. This negative effect of viscoelasticity isin direct contrast with the work done by Chaudhury in the two dimensionalproblem. Lauga, in [25], considered Taylor?s problem in numerous different vis-coelastic fluids. For all cases, he found that the swimming speed was dependenton the viscoelasticity of the fluid and that it was always less that the correspond-4Chapter 1. Introductioning Newtonian swimming speed. He also found that the result was unchangedif more general wave shapes, with both normal and tangential motion, wereconsidered. In all the above viscoelastic problems, the assumption that theamplitude of the swimmer is small compared the wavelength of the swimmerwas employed and so swimming speed is technically referring to leading orderswimming speed.In this thesis, we also consider the problem of an elastic swimmer. Here,both the the unknown swimming speed and swimmer shape are of primaryimportance. Some related work has been done in this area. In [10], the authorsconsidered a single elastic swimmer in a channel with a Newtonian fluid. Theysolved the problem using Peskin?s immersed boundary method. Their resultswere comparable to observed flagella and ciliary motion. There also exists somework done with elastic swimmers in slightly different situations. Lauga in [24]considered a lone, finite elastic swimmer in a Newtonian fluid. This problem ofa low amplitude swimmer was solved using resistive force theory and analyticalformulae were obtained. Work on finite elastic sheets in the presence of wallshas been done as well, namely in the context of settling elastic sheets, or elasticsheets with one end clamped and oscillating vertically [2], [18]. Aside from theformulation of the force balances in these papers, particularly in [2], the resultsare relatively unrelated to our current discussions as they are for sheets andswimmers of finite length.In this thesis, we use analytic methods, rather than numerics so that wemay gain more understanding of the problem, elucidate the dependence onthe physical parameters, and determine the importance of the non-Newtonianmedium.1.2 Overview of Non-Newtonian FluidsWhen modelling fluids, we arm ourselves with equations that describe how thefluid flows. The Navier-Stokes equations governing Newtonian flow in vectorform,rhoparenleftbiggpartialdiffupartialdifft + u ? triangleinvuparenrightbigg= -triangleinvp + ? triangleinv2 u + F. (1.5)They represent the conservation of momentum.1 Here rho is the density of thefluid, u = (u,v) is the velocity field, p is the pressure within the fluid, ? is the1Note that we are considering the two dimensional equations here, not the full three di-mensional ones.5Chapter 1. Introductionviscosity of the fluid and F represents any body forces acting. These equationsmust be supplemented with continuity or conservation of mass,triangleinv? u = 0. (1.6)This is equivalent to assuming the fluid is incompressible or that density isconstant. The question is, how do we describe the flow fluids that have behaviourthat deviates from that of Newtonian fluid?To begin with, we must consider the equations of motion in a more generalform. Still assuming an incompressible fluid, the conservation of mass equationremains unchanged. However, let us write the conservation of momentum equa-tions in a more general manner, referred to as Cauchy?s equations of motion:rhoparenleftbiggpartialdiffupartialdifft + u ? triangleinvuparenrightbigg= -triangleinvp + triangleinv? tau + F. (1.7)Here tau is the deviatoric stress tensor. In the above Navier-Stokes equations fora Newtonian fluid, the relationship between the deviatoric stress and rate ofstrain is assumed to betau = ? ? , (1.8)thus, giving the ? triangleinv2 u term in equation (1.5). This constitutive relationshipbetween stress and rate of strain is unable to describe the properties of non-Newtonian fluids. Let us put for now the constitutive relation between stressand rate of strain in a general form:f(tau, ?) = 0. (1.9)To get an idea of how a non-Newtonian fluid can differ from a Newtonianfluid, let us look at some standard examples of non-Newtonian fluids and relatedbehaviours exhibited by these fluids.To begin with, let us consider a cornstarch solution. This is an interestingexample of what is called a shear thickening fluid.2 Anyone that has ever madeup their own home solution of cornstarch and water knows that when you playwith it, in particular if you hit it or smack it with your palm, it immediatelygets hard like a solid. But, when you move your hand through it gently, itbehaves again like a liquid. Shear thickening means that the viscosity of the2It can be shown that cornstarch solutions have behaviour that is more complicated thatthat which would be expected of a simple shear thickening fluid [3]6Chapter 1. Introductionfluid increases with the shear rate, or deformation rate. Therefore, when yousmack the surface of your cornstarch solution, the rate of deformation increases,and so the viscosity increases, making the solution act momentarily like a solid,or a fluid of higher viscosity. Shear thickening behaviour is often observed infairly concentrated suspensions of small particles.Shear thickening is not the only interesting behaviour that non-Newtonianfluids can exhibit. The opposite effect has also been observed, that of shearthinning. In this instance, the viscosity decreases as the strain rate increases.Examples of fluids that exhibit shear thinning are molten polyethylene andpolypropylene, as well as solutions of caryboxymethylcellulose in water andmany others.There also exists what are called viscoplastic fluids. These fluids exhibit ayield stress. If the yield stress of the fluid is exceeded, then the fluid flows.If the yield stress is not exceeded, then the material acts like a solid. Classic,everyday examples of viscoplastic fluids are mayonnaise and toothpaste.The above behaviours are certainly not the only ones observed in non-Newtonian fluids that make them unique. An important phenomenon exhibitedby visocelastic fluids is recoil, also referred to as fading memory or viscoelas-ticity, as the name would suggest. A good example to illustrate it is flow in acylinder. Consider a viscoelastic fluid that starts at rest and is then subject toa pressure gradient, causing it to flow. This pressure gradient is then removed.What is observed is that the fluid actually recoils or moves backwards towardsits original state. The amount of recoil in this situation is dependent on the spe-cific viscoelastic fluid that is used. Some viscoelastic fluids exhibit very strongrecoil, while with others it is less prominent. It is easy to see why this is oftenreferred to as fading memory. A fully elastic material would return fully to itsoriginal state.Two common examples of phenomena observed in viscoelastic fluids are jetswell and rod climbing. After exiting an orifice in a downward direction, aNewtonian fluid has a stream of a set diameter that is the same as that of thejet. In comparison, viscoelastic fluids have a diameter that increases twofoldfrom the orifice diameter. This is also called the Barus or Merrington effect.Rod climbing, or the Weissenburg effect, refers to the situation where a thinsolid cylinder is rotating in a viscoelastic fluid. It is observed that the fluidclimbs up the rod, moving in the opposite direction to a Newtonian fluid inthe same situation. A Newtonian fluid in this situation, besides moving in adifferent direction, does not climb the rod. In fact, a small indentation at the7Chapter 1. Introductionbase of the rod is observed.Due to the composition of cervical mucous, namely that it is composed ofvarious macromolecules, it is important to mention polymeric fluids.3 Poly-meric fluids are a specific type of viscoelastic fluid. Unlike Newtonian, fluids(and many non-Newtonian fluids), that are composed of micromolecules, poly-meric fluids are composed of marcomolecules.4 The rheological properties of apolymeric fluid are very dependent on the molecular weight distribution of themacromolecules that comprise it. That is, the rheological properties of a poly-meric fluid depend not only on the number of molecules composing it, but thesizes and locations of those molecules as well. Polymeric fluids can exhibit manyof the behaviours discussed above, including shear thinning and thickening, jetswell and rod climbing, and, obviously, viscoelasticity.51.2.1 Constitutive ModelsThere are many different ways to model non-Newtonian fluids ranging fromrelatively simple, to immensely complex, depending on the properties of thespecific fluid one is interested in. On the simpler end of the scale there aregeneralized Newtonian fluid models and, for viscoelastic fluids, general linearviscoelastic fluid models. On the complex end of things there are corotationaland codeformational models.To begin with, let us look at generalized Newtonian fluids. This is thesimple case when you take the constitutive model for a Newtonian fluid, as seenin equation (1.8), and relax the assumption that the viscosity, ?, is constant.The constant viscosity is replaced with a viscosity that varies with the shearrate, ?, or even sometimes with |tauxy|. A standard of this type of constitutivelaw is the power law. This law states that? = m?n-1, (1.10)where m and n are determined by the specific fluid under consideration.Viscoplastic fluids exhibit a yield stress. A common example of a viscoplastic3Cervical mucous is a gel that consists of a network of glycoprotein macromolecules and aplasma component that contains soluble proteins, inorganic salts and other chemical species[21].4Undiluted polymer solutions are typically referred to as melts. In any concentrated poly-mer solutions, the viscosity depends on concentration in a highly nonlinear fashion. In dilutesolutions, the macromolecules are far apart from each other and therefore have negligibleinfluence on each other.5For a full, in depth discussion on polymeric fluids and the associated phenomena, see [5]8Chapter 1. Introductionfluid is a Bingham fluid, which is modelled by? = infinity, tau <=< tauy,? = ?0 + tauy?gamma , otherwise. (1.11)Another example of a viscoplastic fluid is a Herschel-Bulkley fluid, which alsoincorporates shear thinning, given by? = infinity, tau <=< tauy,? = (tauy + ?| ?|n)sgn(?), otherwise. (1.12)Generalized linear viscoelastic fluid models are simple models of viscoelasticfluids that display some of the time dependent properties (i.e. recoil or viscoelas-ticity). These, however, are only applicable to situations that have low strainsand shear rates. The most standard example of a generalized linear viscoelasticmodel is Maxwell?s modeltau + lambda partialdiffpartialdiffttau = ??, (1.13)where lambda is a time constant and ? has units of viscosity. Another generalizedlinear viscoelastic model is the Jeffreys model,tau + lambda1 partialdiffpartialdiffttau = ?parenleftbigg? + partialdiffpartialdifft ?parenrightbigg. (1.14)As you can see, this is just Maxwell?s model with an extra term added in on theright-hand side, a time derivative of the strain rate and another time constant,lambda2. The Jeffreys model is important because it has been the starting point forother, more complicated rheological equations of state.Both of these types of models are limited in scope. Ideally, we want a modelthat describes all flow phenomena, at least qualitatively. Two types of modelsthat do this are corotational models and codeformational models. These arebasically just versions of the linear viscoelastic models revised so that they arevalid in situations when the fluid gets largely deformed from its initial position.Corotational models are formulated in a reference frame that corotates with thefluid. That is, the coordinates, call them thetaj, have a time rate of change givenbypartialdiffpartialdifftthetaj = -12omega ? thetaj, (1.15)where omega = triangleinv? u is the fluid vorticity. Observers in this frame will not detect9Chapter 1. Introductionrigid rotation. Codeformational models are formulated in codeforming coordi-nates, call them phij, with time rate of changepartialdiffpartialdifftphij = -phij ? triangleinvu. (1.16)The models being used in this thesis are codeformational. For a more in depthdiscussion and analysis of these models, the reader is referred to [5] and [31].Of specific interest are the Oldroyd models. The Oldroyd models can be bothcodeformational and corotational and are based on the Jeffreys model. We willrestrict our discussion here to the codeformational type, since that is the typethat will be used in this thesis. To begin with, let us look at the Oldroyd-Bmodel. This model is in fact the Jeffreys model generalized in codeformationalcoordinates.6 It is given by,tau + lambda1tautriangleinv = ?parenleftBig? + lambda2 ?triangleinvparenrightBig. (1.17)Here, lambda1 is the polymer relaxation time, and lambda2, another relaxation time, rep-resents the rate of decay of the residual rate of strain when the fluid is stressfree, or the retardation time scale. The upper convected derivative of a tensorA, is defined to beAtriangleinv = partialdiffApartialdifft + u ? triangleinvA -parenleftbig(triangleinvu)T ? A + A ? triangleinvuparenrightbig . (1.18)It is possible to eliminate one relaxation time and write it in terms of the other bythinking of the fluid as being composed of a Newtonian solvent and a polymericsolute. This allows one to write the viscosities and deviatoric stresses as sums.We put ? = ?s + ?p and tau = taus + taup. In this situation, it is generally assumedthat ? = lambda1lambda2 ?s. This assumption allows us to remove one relaxation numberand also allows us to get rid of the convected strain rate term. Applying theseassumptions to the above Oldroyd-B constitutive equation simplifies it to theupper convected Maxwell model for the polymeric component,taup + lambda1tautriangleinvp = ?p ? . (1.19)6Technically, there are two possible generalizations of the Jeffreys model in codeformationalcoordinates, the Oldroyd-A model and the Oldroyd-B model. Oldroyd-A is very similar toOldroyd-B, with the difference being in the convected derivative. Oldroyd-A is given bytau + lambda1parenleftBigDtauDt + tau ? (triangleinvu)T + triangleinvu ? tauparenrightBig= ?parenleftBig? + lambda2parenleftBig D?gammaDt + ?gamma ? (triangleinvu)T + triangleinvu ? ?gammaparenrightBigparenrightBig.10Chapter 1. IntroductionAnother Oldroyd model is the Oldroyd-8 model, aptly named as there areeight parameters in it. This model is also a generalization of the Jeffreys model,in which all possible quadratically nonlinear terms involving the deformationrates allowed by the symmetries of the problem are included. This model isgiven by,tau + lambda1tautriangleinv + 12?0(trtau)? - 12?1 (tau ? ? + ? ? tau) + 12nu1(tau : ?)delta =?parenleftbigg? + lambda2 ?triangleinv -?2( ? ? ?) + 12(? : ?)deltaparenrightbigg. (1.20)The eight constants and extra terms in this model allow for a much largervariety in rheological response compared to the Jeffreys model and the Oldroyd-B model [5]. This model is also convenient as it contains several other modelswith in it that can be obtained by merely setting some of the constants equal tozero, or multiples of each other. For example, it contains the Oldroyd-B modelby setting all constants to zero except lambda1 and lambda2. It also contains the upperconvected Maxwell model, the second-order fluid model, and the Oldroyd fourand six constant models, to name a few.Codeformational constitutive models can get significantly more complicatedthan those discussed above. However, a full discussion these rheological equa-tions of state is not necessary here, as the only codeformational models we willbe using are Oldroyd-B and Oldroyd-8 constitutive models.1.3 Our Problem: Main AssumptionsNow that we are somewhat comfortable with viscoelastic fluids, we are in aposition to set up our main problem mathematically.Consider a two-dimensional, infinite periodic swimmer in the presence of awall. We would like to consider the case when the fluid between the wall andthe swimmer is viscoelastic. This is because the physical properties of cervicalmucous are viscoelastic. We are assuming that the swimmer propagates wavesdown its length therefore propelling itself in the opposite direction. The shapeof the swimmer is referred to as Y and is centered about y = 0. The wall isfixed at y = h, which is also the mean distance from the swimmer to the wall.The fact that the swimmer is periodic tells us that angbracketleftYxangbracketright = 0 and we assume that11Chapter 1. IntroductionangbracketleftY angbracketright = 0, so that h is the mean gap width. Hereangbracketleft...angbracketright = 12piintegraldisplay 2pi/k0(...)dx. (1.21)Figure 1.1: Geometry of ProblemThe equations governing such a flow aretriangleinv?u = 0, (1.22)rhoparenleftbiggpartialdiffupartialdifft + u ? triangleinvuparenrightbigg= -triangleinvp + triangleinv? tau + F, (1.23)f(tau, ?) = 0. (1.24)Recall that the constitutive relation between stress and rate of strain in a generalform. It is assumed that the swimmer propagates waves down its length in thepositive x-direction, therefore propelling itself in the negative x-direction. Theswimming speed is referred to as U, thus its velocity is -U. We assume thatU > 0. Hence, we have boundary conditions as follows. For the horizontalvelocity component,u(x,h) = 0, (1.25)u(x,Y ) = -U, (1.26)and for the vertical velocity component, assuming simple vertical motion rather12Chapter 1. Introductionthan complex elliptical trajectories of a material point on the swimmer we have,v(x,h) = -UYx + Yt, (1.27)v(x,Y ) = Yt. (1.28)We will, however, not be solving this system in its entirety. We will makeseveral simplifying assumptions to help ourselves along.Due to the physical situation under consideration, a low Reynolds numberregime is considered0 < Re lessmuch1, (1.29)allowing the inertial terms or the left hand side of Cauchy?s equation of motionto be neglected. It will be assumed that there are no external forces acting inthe problemF = 0, (1.30)that the unknown swimming speed of the sheet is constant and positive,U = Const > 0. (1.31)It will also be assumed that the amplitude of the propagating waves is smallcompared to the wavelength.7 It has been observed that in low viscosity medi-ums, activated (rather than hyperactivated) sperm generate relatively low am-plitude, symmetrical waves.8 This low amplitude assumption can be formulatedmathematically by stating that0 < parenleftbigmaxxelement[0,2pi]|Y |parenrightbig k = epsilon lessmuch1, (1.32)where k is the wavenumber associated with Y . It is also assumed that theproblem in its entirety is periodic in (kx - (omega - Uk)t). Of specific use will bethe periodicity of the pressure, angbracketleftpxangbracketright = 0.To simplify things, we will consider the frame moving with the swimmer.7Note that this is not the same as saying that the amplitude of the propagating wave issmall compared to the mean distance from the swimmer to the wall, or that a/h is small.8Generally, in oviduct, sperm are hyperactivated, thus producing larger amplitude, moreasymmetrical waves [36].13Chapter 1. IntroductionThat is, putx = ? -Ut, (1.33)dxdt =d?dt -U. (1.34)This implies that the solutions will now be periodic in (k? - omegat). For theremainder of this thesis, we will assume that it is understood we are in themoving frame, and so will drop the tildes on the above changed variables.Figure 1.2: Geometry of Problem in Moving Frame of ReferenceThis simplifies the boundary conditions. In the new frame of reference, forthe horizontal velocity component, we haveu(x,h) = U, (1.35)u(x,Y ) = 0. (1.36)For the vertical velocity component,v(x,h) = 0, (1.37)v(x,Y ) = Yt. (1.38)So far, the equations of motion, most assumptions of the problem and thenew frame of reference have been specified.14Chapter 1. Introduction1.3.1 Zero Force on Swimmer ConditionIn addition to the above kinematical assumptions, we will also introduce a dy-namical one. We assume that the average force per wavelength on the swimmeris zero. This is a reasonable assumption due to the fact that it has already beenassumed that the speed of the swimmer is constant [22].To formulate this condition precisely, we must establish what the force bal-ances are in the problem.We know that the force on the swimmer by the fluid must be equal to minusthe total force on the fluid by the swimmer (per unit length). Thus, we canformulate this condition by considering the total force on the fluid.It known that the total stress tensor, sigma, is given bysigma = -pI + tau =parenleftBigg-p + tauxx tauxytauyx -p + tauyyparenrightBigg,and it is found that the unit outward normal to the swimmer is1radicalbig1 + Y 2xparenleftBiggYx-1parenrightBigg.We can write the total force on the fluid per unit length asintegraldisplay 2pi0parenleftBigg-p + tauxx tauxytauyx -p + tauyyparenrightBiggparenleftBiggYx-1parenrightBiggdx.Do to the force balance, we know that this is equal to minus the force on theswimmer, which is assumed to be zero. Thus, by writing things in componentform, one can see that (recall the stress tensor is symmetric)angbracketleft(-p + tauxx)Yx -tauxyangbracketright = 0,angbracketlefttauxyYx -(-p + tauyy)angbracketright = 0,where the top equation is the tangential force, the bottom the normal force,recalling that angbracketleft...angbracketright denotes the x-average.Now, consider the conservation of momentum equations in component form,slightly rearranged(p -tauxx)x = tauxy,y,(p + tauyy)y = tauxy,x.15Chapter 1. IntroductionConcerning ourselves with only the first equation 9 and by applying the averagein x while maintaining the assumed periodicity of the problem, we discoverangbracketlefttauxyangbracketrighty = 0,and by integratingangbracketlefttauxyangbracketright = Const.Now, consider the same momentum equation. Rather than applying an averagein x, integrate in y. This givesintegraldisplay hkY(p -tauxx)xdy = tauxy(hk) -tauxy(Y ).Pulling the x-derivative out of the integral,partialdiffxintegraldisplay hkY(p -tauxx)dy + (p -tauxx)vextendsinglevextendsingley=Y Yx = tauxy(hk) -tauxy(Y ),and again average in x to obtain,angbracketlefttauxy(hk)angbracketright = angbracketlefttauxy(Y )angbracketright + angbracketleft(p -tauxx)vextendsinglevextendsingley=Y Yxangbracketright.But, notice that the right-hand side of this above equation is merely the tan-gential force balance equation evaluated at y = Y . Knowing then that this forceis equal to zero leavesangbracketlefttauxy(hk)angbracketright = 0,and since it was calculated above that angbracketlefttauxyangbracketright = Const, we finally haveangbracketlefttauxyangbracketright = 0. (1.39)This is a constraint on the system that will be used in both Chapters Two andThree.9At this time we will ignore the vertical force balance due to the fact that either theswimmer is assumed to be fixed vertically by some external force or a more complex verticalforce balance is required utilizing elasticity theory, as will be seen in Chapter Three.16Chapter 1. Introduction1.4 Preliminary ProblemIn this section, we discuss the fundamental problem of a swimmer in the presenceof a wall and a viscoelastic fluid. We do not apply the lubrication approximationhere. We will assume that the amplitude of the waves is small in comparison tothe mean distance to the wall and that the fluid between the swimmer and thewall is an Oldroyd-B viscoelastic fluid. This is a nice warm up problem as it isa combination of Katz?s biharmonic problem and Lauga?s viscoelastic version ofthe Taylor problem [22, 25]. The method of solution applied here follows closelywith that used in [25]. See figure (1.2) for geometry.1.4.1 Mathematical FormulationWe assume that the sheet has prescribed motion, that of waves propagating tothe right, see figure (1.2). In the frame moving with the swimmer, these waveslook likey = asin(kx -omegat). (1.40)The equations of motion for this problem aretriangleinv?u = 0, (1.41)triangleinvp = triangleinv? tau, (1.42)tau + lambda1tautriangleinv = ?parenleftBig? + lambda2 ?triangleinvparenrightBig, (1.43)which are the conservation of mass, momentum and the Oldroyd-B constitutiveequation respectively. 10 Recall that the upper convected derivative of a tensorA, is defined to beAtriangleinv = partialdiffApartialdifft + u ? triangleinvA -parenleftbig(triangleinvu)T ? A + A ? triangleinvuparenrightbig . (1.44)We also have the zero force condition11angbracketlefttauxyangbracketright = 0. (1.45)10Recall that lambda1 is the polymer relaxation time, and lambda2 represents the rate of decay of theresidual rate of strain when the fluid is stress free, or the retardation time scale.11This equation still holds even though there is no lubrication approximation here due tothe low amplitude compared to mean height approximation.17Chapter 1. IntroductionNon-dimensionalizing the equations as followst = 1omega ?, (x,y) = 1k(?, ?),u = c?, ? = omega?,(p,tau) = ?omega(?, ?), ? = Ucwhere, c = omegak , is the speed of the travelling wave, and dropping the hat notation,yieldstriangleinvp = triangleinv? tau, (1.46)triangleinv? u = 0, (1.47)which are unchanged, andtau + De1tautriangleinv = ? + De2 ? triangleinv, (1.48)the constitutive equation, where De1 and De2 are Deborah numbers. TheDeborah number is the ratio of the characteristic time for the fluid over thecharacteristic time for the flow system. Here, we have two characteristic timesfor the fluid, and hence two Deborah numbers.It remains to non-dimensionalize the boundary conditions. Using the lowamplitude approximation, we putepsilon = ak lessmuch1, (1.49)and find that (dropping hats):y = epsilonsin(x -t). (1.50)The boundary conditions become, for the horizontal velocity componentu(x,hk) = U, (1.51)u(x,Y ) = 0, (1.52)and for the vertical velocity componentv(x,hk) = 0, (1.53)v(x,Y ) = -epsiloncos(x -t). (1.54)18Chapter 1. Introduction1.4.2 Stream Function and AsymptoticsSince the problem at hand is an incompressible, two-dimensional problem, it ismost useful to consider a stream function psi such that 12u = partialdiffpsipartialdiffy , v = -partialdiffpsipartialdiffx . (1.55)To solve this system, consider low amplitude solutions of the formU = epsilonU1 + epsilon2U2 + ... (1.56)psi = epsilonpsi1 + epsilon2psi2 + ... (1.57)tau = epsilontau1 + epsilon2tau2 + ... (1.58)That is, regular perturbation expansions in epsilon for the swimming speed, the streamfunction and the stress.Since it is assumed that the speed of the swimmer is constant, it must bethe case that each Uk in the expansion of U is also constant. Notice that, sincethe velocity field components can be written in terms of the stream functionand ? can be written in terms of the velocity field components, we can put? = epsilon ?1 + epsilon2 ?2 + ... (1.59)u = epsilonu1 + epsilon2u2 + ... (1.60)Removing the pressure term from the conservation of momentum equations andthen summarizing the equations of motion givestriangleinv?triangleinv?tau = 0, (1.61)tau + De1tautriangleinv = ? + De2 ? triangleinv. (1.62)The next step is to then substitute in the perturbation expansions and equateorders of epsilon.12Conservation of mass holds automatically sincetriangleinv? u = ux + vy = partialdiffpartialdiffxparenleftBig partialdiffpsipartialdiffyparenrightBig+ partialdiffpartialdiffyparenleftBig-partialdiffpsipartialdiffxparenrightBig= 0.19Chapter 1. Introduction1.4.3 Leading Order EquationsAt O(epsilon) we have:triangleinv?triangleinv?tau1 = 0, (1.63)tau1 + De1 partialdiffpartialdiffttau1 = ?1 + De2 partialdiffpartialdifft ?1. (1.64)We find the boundary conditions at this order by substituting in the regularperturbation expansion for psi and expanding about y = 0. This givespsi1xvextendsinglevextendsingle(x,0) = 0, psi1yvextendsinglevextendsingle(x,0) = 0,psi1xvextendsinglevextendsingle(x,hk) = cos(x -t), psi1yvextendsinglevextendsingle(x,hk) = U1. (1.65)The zero force constraint angbracketlefttauxyangbracketright = 0 yields thatangbracketlefttau1xyangbracketright = 0. (1.66)If one stares at these leading order equations for a moment, it becomesapparent that the easiest thing to do is to take the divergence and then the curlof equation (1.64),thus eliminating leading order stress. By doing this we obtainparenleftbigg1 + De2 partialdiffpartialdifftparenrightbiggtriangleinv4 psi1 = 0, (1.67)where triangleinv4psi = psixxxx + 2psixxyy + psiyyyy.To solve this equation, we simply use separation of variables with hyperbolicsines and cosines. Assume that the form of the solution ispsi1(y) = [(A + BH)cosh(H) + (C + DH)sinh(H)]sin(x -t),where H = hk -y. We find then thatU1 = 0, (1.68)andpsi1 = W(y)W(0) sin(x -t), (1.69)20Chapter 1. IntroductionwherePsi1(y) = (hk -y)sinh(hk -y), (1.70)Psi2(y) = sinh(hk -y) -(hk -y)cosh(hk -y), (1.71)andW(y) = Psiprime1(0)Psi2(y) -Psiprime2(0)Psi1(y). (1.72)So, the leading order swimming speed is zero, a result that is precedented. Thus,we must look to second order to find the speed. To do this, we must calculatethe leading order stress.To solve for the leading order stress, we employ the method of ?complexifica-tion?. That is, we think of the stress and rate of strain as being the real portionof a larger, complex problem that satisfies the same equations.13 We puttau1 = Rfractur[ ?1ei(x-t)], ?1 = Rfractur[ ?1ei(x-t)]. (1.73)By employing this method of solution, we can easily obtain?1 = 1 -iDe21 -iDe1??gamma1, (1.74)and upon substituting in ?1, we gettau1 = RfracturbraceleftBigg(1 -iDe2)ei(x-t)(1 -iDe1)W(0)parenleftBigg2Wprime -i(Wprimeprime + W)-i(Wprimeprime + W) -2WprimeparenrightBiggbracerightBigg. (1.75)Notice that the condition,angbracketlefttau1xyangbracketright = 0,automatically holds due to the fact that tau1xy proportional ?1xy proportional sin(x -t).1.4.4 Second Order EquationsThe leading order problem has been solved and it has been shown that in orderto maintain the assumption that the speed of the swimmer is constant, at leadingorder it must be equal to zero. Thus, in order to calculate the swimming speedof the organism, we must move on to the second order problem. This problem is13Note the periodic dependence explicitly given due to the assumed periodicity of the prob-lem21Chapter 1. Introductionnot quite as straightforward as the leading order problem due to the nonlinearterms of the upper convected derivative. At O(epsilon2) we havetriangleinv?triangleinv?tau2 = 0, (1.76)tau2 + De1 partialdiffpartialdiffttau2 + De1 parenleftbigu1 ? triangleinvtau1 -tau1 ? triangleinvu1 -(triangleinvu1)T ? tau1parenrightbig =?2 + De2 partialdiffpartialdifft ?2 + De2 parenleftbigu1 ? triangleinv ?1 - ?1 ? triangleinvu1 -(triangleinvu1)T ? ?1parenrightbig , (1.77)with boundary conditionspsi2xvextendsinglevextendsingle(x,0) = 0, psi2yvextendsinglevextendsingle(x,0) = -sin2(x -t)Wprimeprime(0)W(0) ,psi2xvextendsinglevextendsingle(x,hk) = 0, psi2yvextendsinglevextendsingle(x,hk) = U2. (1.78)The zero force condition states thatangbracketlefttau2xyangbracketright = 0. (1.79)Inspecting the second order equations of motion reveals that the right-handside of equation (1.77) is composed of leading order terms exclusively. To solve,let us put the right-hand side of the momentum equations in terms of the com-plexified variables and look only at the xy-componentparenleftbigg1 + De1 partialdiffpartialdifftparenrightbiggtau2xy -parenleftbigg1 + De2 partialdiffpartialdifftparenrightbigg?2xy =Rfracturbraceleftbigg(De2 -De1)2(1 -iDe1)bracketleftbig?uasteriskmath1 ? triangleinv??gamma1 - ??gamma1 ? triangleinv?uasteriskmath1 -(triangleinv?uasteriskmath1)T ? ??gamma1bracketrightbigbracerightbiggxy++Rfracturbraceleftbigg(De2 -De1)2(1 -iDe1) e2i(x-t) bracketleftbig?u1 ? triangleinv??gamma1 - ??gamma1 ? triangleinv?u1 -(triangleinv?u1)T ? ??gamma1bracketrightbigbracerightbiggxy(1.80)This allows us to utilize equation (1.79) as a simplification. To apply thiscondition, let us take the x-average of equation (1.80). Then, by writing therate of strain component in terms of the stream function, and observing the factthat angbracketleftpsi2xxangbracketright = 12pipsi2xvextendsinglevextendsingle2pi0 = - 12pivvextendsinglevextendsingle2pi0 = 0 due to the assumed periodicity of thesolution, we obtainangbracketleftpsi2angbracketrightyy = Rfracturbraceleftbigg(De1 -De2)2(1 -iDe1)bracketleftbig?uasteriskmath1 ? triangleinv??gamma1 - ??gamma1 ? triangleinv?uasteriskmath1 -(triangleinv?uasteriskmath1)T ? ??gamma1bracketrightbigxybracerightbigg. (1.81)22Chapter 1. IntroductionTo solve the problem and calculate the swimming speed we must integratethe above expression for angbracketleftpsi2angbracketrightyy once to obtain an expression for angbracketleftpsi2angbracketrighty and applythe boundary conditions. To be able to apply the boundary conditions, we mustput them into the same averaged form as the equation. For the conditions onpsix, we obtain trivial averages due to the periodicity of the problem, angbracketleftpsi2xangbracketright =12pipsi2vextendsinglevextendsingle2pi0 = 0, and thus no information. However, for the other two boundaryconditions we obtainangbracketleftpsi2yangbracketrightvextendsinglevextendsingley=0 = -Wprimeprime(0)2W(0), (1.82)angbracketleftpsi2yangbracketrightvextendsinglevextendsingley=hk = U2. (1.83)By integrating we find that, writing W(0) = W0 and Wprimeprime(0) = Wprimeprime0 for simplicity,angbracketleftpsi2angbracketrighty = De1(De2 -De1)2(1 + De21)W2bracketleftbigWWprimeprime + Wprime2bracketrightbig - Wprimeprime02W0(1 + De1De2)(1 + De21) . (1.84)Finally, the swimming speed is obtained,U2 = angbracketleftpsi2angbracketrightyvextendsinglevextendsingley=hk = - Wprimeprime02W0(1 + De1De2)(1 + De21) . (1.85)Note that W(hk) = 0 and Wprime(hk) = 0.1.4.5 Limiting CasesDistance Between Wall and Swimmer Approaches InfinityWe would like to assure ourselves that this problem agrees with the infiniteproblem considered Lauga, in [25]. To do this, we take the limit as the distancebetween the wall and the swimmer approaches infinity. It is known thatpsi1 = W(y)W(0) sin(x -t),and so the limit as h -arrowrightinfinity is easily calculated. We find thatlimh-arrowrightinfinityW(y)W(0) = (1 + y)e-y.Thus,limh-arrowrightinfinitypsi1 = (1 + y)e-y sin(x -t),23Chapter 1. Introductionwhich is exactly the infinite solution as stated in [25].It is obvious that if the limit of the streamfunction agrees with the infinitecase that the limit of the swimming speed will as well, as psi2y(hk) = U2. Whencalculated we find thatlimh-arrowrightinfinityU2 = 12 (1 + De1De2)(1 + De21),which again agrees with the Lauga?s result.Similarly, we can calculate the limit for the averaged rate of work done bythe swimmer. We find thatlimh-arrowrightinfinitybraceleftBigg(1 + De1De2)(1 + De21)W20integraldisplay hk0bracketleftbig4Wprime2 + (W + Wprimeprime)2bracketrightbig dybracerightBigg= 2(1 + De1De2)(1 + De21),which again agrees with the previous work done.14So, it is clear that the above calculations for the problem of an organismswimming in an Oldroyd-B viscoelastic fluid near a wall are consistent withthose of the problem of an organism swimming in an infinite fluid of the sametype.Swimmer Approaches the WallNext, we would like to consider the limit as the swimmer and the wall approacheach other. There are two ways of thinking of this limit. The first is keeping thewall fixed at y = h and increasing the amplitude of the swimmer until it reachesthe wall. This line of thinking contradicts the assumption that ak << 1, as itimplies that O(a) arrowright O(h) = O(1/k). This would make our present solutioninvalid. The second way of thinking about this limit is to think of the wallapproaching the swimmer, whose amplitude and wavelength remains the fixed.This implies that O(h) arrowrightO(a), still maintaining the assumption that ak << 1.In fact, in this manner of thinking we get that hk << 1 which is the sameassumption that is applied when making the lubrication approximation. Takingthis limit in this case, however, is not the same as applying the lubricationapproximation fully. It is just a limiting case of the problem under discussion,and will suggest to us to fully apply the lubrication approximation to understand14Please see Appendix A for a full derivation of this problem and the calculation of thework done.24Chapter 1. Introductionthis limit better. We have then, to leading order of this limiting case,U = -epsilon2 (1 + De1De2)2(1 + De21)F(hk), (1.86)whereF(hk) = Wprimeprime0W0 =(hk)2 + sinh2(hk)(hk)2 -sinh2(hk). (1.87)In our limit, hk arrowrightak = epsilon and so we have thatF(epsilon) = epsilon2 + sinh2 epsilonepsilon2 -sinh2 epsilon. (1.88)Expanding sinhepsilon = epsilon + 16epsilon3 + ... and simplifying gives usF(epsilon) = 2epsilon2 + +13epsilon4 + ...-13epsilon4 -... . (1.89)Thus, in the limit when the wall approaches the swimmer, we find that theswimming speed isU = 3(1 + De1De2)(1 + De21), (1.90)which is significantly larger order than the original solution.15 Thus, this prob-lem suggests a more in depth analysis of the problem with the lubrication ap-proximation specifically applied from the start of the formulation of the problem.For a more comprehensive, detailed discussion of this problem please seeAppendix A.1.5 Thesis OutlineChapter Two presents the problem of a infinite swimming sheet with the lubri-cation approximation applied. The sheet shape here is assumed. This problemis first solved for a Newtonian fluid in a preliminary sense, and is then solved fora viscoelastic fluid represented by the Oldroyd-8 constitutive model. ChapterThree presents the problem of a infinite swimming elastic sheet with the lu-brication approximation applied. In this problem, there is a prescribed forcing15Note that in this limit we do not take into account what is happening to the shear rates,which will be getting significantly larger. This is just a motivational example to show thatwhen the swimmer approaches the wall, the swimming speed is not longer proportional to theamplitude of the swimmer squared and something more interesting is happening.25Chapter 1. Introductionalong the length of the sheet. This problem is solved for a Newtonian fluid forforcing of any strength, and then for a viscoelastic fluid, also represented bythe Oldroyd-8 constitutive model, with forcing of moderate strength. The finalchapter presents a summary of the findings and concluding remarks, includingfurther work that may be done.26Chapter 2Swimming with LubricationTheoryAs motivated by the preliminary problem discussed above, we now considerthe problem of a two dimensional swimmer very near a wall. Particularly, weapply the lubrication approximation and assume that the mean distance of theswimmer from the wall, h, compared to the wavelength of the swimmer is verysmall.16Figure 2.1: Lubrication geometry2.1 Mathematical FormulationLet us consider a swimmer in the presence of a wall with a fluid in between.The swimmer shape is assumed to be known. We put it to bey = asin(kx -omegat). (2.1)We consider this wave form and not a more general one because it has beenshown by Katz in [22] that for the Newtonian lubrication problem, horizontalcontributions to the wave shape only affect the higher orders of the solution.For the viscoelastic problem of no wall, Lauga in [25] considered horizontal16For a mathematical discussion of viscoelastic lubrication theory, we refer the reader to [4]27Chapter 2. Swimming with Lubrication Theoryoscillations of the same order as the amplitude in the sheet and found that theydid not affect the leading order solution. Thus, it is very reasonable to assumethat horizontal contributions to the wave form will not affect our leading ordersolution, and so they will be neglected.Here, we are applying the lubrication approximation. Mathematically thismeans thathk lessmuch1. (2.2)The equations for this problem, using the general assumptions stated inChapter One, aretriangleinv?u = 0, (2.3)triangleinvp = triangleinv? tau, (2.4)f(tau, ?) = 0, (2.5)which are the conservation of mass, momentum and the general constitutiveequations respectively.17 This problem is also subject to the constraint thatangbracketlefttauxyangbracketright = 0, (2.6)from the zero force condition. As well, we require that the periodicity in pressurebe considered explicitly for this problem, thus we will includeangbracketleftpxangbracketright = 0, (2.7)in our formulation.The boundary conditions are, for the horizontal velocity componentu(x,h) = U, (2.8)u(x,Y ) = 0, (2.9)and for the vertical velocity componentv(x,h) = 0, (2.10)v(x,Y ) = Yt, (2.11)17Note that due to the low Reynold?s number assumption, the inertial terms are neglectedin the conservation of momentum equation.28Chapter 2. Swimming with Lubrication Theoryas we can see in (2.1). Recall that Y is the shape of the swimmer.Let us non-dimensionalize these equations and apply the lubrication approx-imation. We will use the following non-dimensional variables [42]:t = 1omega ?, ? = omega?,x = 1k ?, y = h?,u = c?,(p,tau) = ?ch (1epsilon ?, ?), ? = Ucwhere c = omegak is the speed of the travelling wave.Upon substituting in the above variables (dropping hats), we find that theconservation of mass equation remains unchangedux + vy = 0. (2.12)Substitution of the non-dimensional variables into the conservation of momen-tum yields (again, hats dropped),px = epsilontauxx,x + tauyx,y,py = epsilon2tauxy,x + epsilontauyy,y.Thus, at leading order we have,px = tauyx,y, (2.13)py = 0. (2.14)This tells us that pressure is independent of y and that the xy-component of thestress tensor is the dominant component. The shape of the swimmer becomesY = ah sin(x -t). (2.15)We can then note that, Yt = -Yx, when in dimensionless variables.Observe that if we integrate the equation px = tauyx,yintegraldisplay yYpxdy =integraldisplay yYtauyx,ydy,we get thattauxy = tau0 + (y -Y )px, (2.16)29Chapter 2. Swimming with Lubrication Theorywhere,tau0 = tauxy|y=Y ,is the stress on the swimmer. If we observe then thatangbracketlefttauxyangbracketright = angbracketlefttau0 + (y -Y )pxangbracketright = angbracketlefttau0 -Y pxangbracketright + yangbracketleftpxangbracketrightand apply the periodic pressure condition, angbracketleftpxangbracketright = 0 which remains unchangedin dimensionless variables, the zero force condition simplifies toangbracketlefttau0 -pxY angbracketright = 0. (2.17)Note that since we havepx = tau1 -tau01 -Y ,where tau1 = tauxy|y=1, it is the case thattau0 -Y px = tau1 -px.Thus, we have thatangbracketlefttau0 -Y pxangbracketright = angbracketlefttau1angbracketright = 0, (2.18)since we are assuming that angbracketleftpxangbracketright = 0.So, our final equations are thenux + vy = 0, (2.19)tauxy = tau0 + (y -Y )px, (2.20)fparenleftBig?, ?parenrightBig= 0, (2.21)angbracketleftpxangbracketright = 0, (2.22)angbracketlefttau1angbracketright = 0, (2.23)and non-dimensional boundary conditionsu(x,1) = U, (2.24)u(x,Y ) = 0, (2.25)30Chapter 2. Swimming with Lubrication Theoryand for the vertical velocity componentv(x,1) = 0, (2.26)v(x,Y ) = -Yx. (2.27)2.2 Newtonian ProblemDue to this problem?s complexity, let us first solve it with a Newtonian fluid.Not only will this give us some insight into the more general problem of anon-Newtonian fluid, but the results will be used as an initial guess for thenumerical method used to solve the viscoelastic lubrication problem, detailed inthe following section. For the viscous problem we have the following:ux + vy = 0, (2.28)px = uyy, (2.29)py = 0, (2.30)angbracketleftuy|y=Y -pxY angbracketright = 0, (2.31)angbracketleftpxangbracketright = 0, (2.32)and boundary conditions given by equations (2.24), (2.25), (2.26), and (2.27).To solve, let us first integrate px = uyy twice in y exploiting the fact that pxis independent of y and apply the boundary conditions on u to obtainu = 12px(y -Y )(y -1) + U (y -Y )(1 -Y ). (2.33)Then, look at the continuity equation. Integrate in y first to obtainintegraldisplay 1Yuxdy + Yx = 0.Pulling the x derivative out of the integral, noting that u|y=Y = 0, and inte-grating in x gives integraldisplay1Yudy + Y = q,where q is the mass flux. Upon calculating this integral and solving for px we31Chapter 2. Swimming with Lubrication Theoryfindpx = 12 (1 -q)(1 -Y )3 + 6 (U -2)(1 -Y )2 . (2.34)Thus, it is left only to determine q and U. To do this, we will simply applythe constraints angbracketleftpxangbracketright = 0, angbracketleftuy|y=Y - pxY angbracketright = 0 and solve simultaneously. Thisinvolves solving integrals of the formIasteriskmathj =integraldisplay 2pi0dx(1 -Y )j j = 1,2,3 (2.35)This gives us thatU = q = 6parenleftbigIasteriskmath1 Iasteriskmath -Iasteriskmath2parenrightbig(4Iasteriskmath1 Iasteriskmath3 -3Iasteriskmath22 ). (2.36)Solving the integrals specifically givesIasteriskmath1 = 1(1 -parenleftbigahparenrightbig2)1/2, (2.37)Iasteriskmath2 = 1(1 -parenleftbigahparenrightbig2)3/2, (2.38)Iasteriskmath3 = (parenleftbigahparenrightbig2 + 2)(1 -parenleftbigahparenrightbig2)5/2, (2.39)which yieldsU = q = 3parenleftbig ahparenrightbig2parenleftBig2parenleftbig ahparenrightbig2 + 1parenrightBig. (2.40)Dimensionally,Uc =3parenleftbig ahparenrightbig2parenleftBig2parenleftbig ahparenrightbig2 + 1parenrightBig. (2.41)This is exactly the result found by Katz, in [22]. By writing this in the formUc =32 + parenleftbighaparenrightbig2, (2.42)we can easily see that the swimming speed is bounded by the wave speed.Substituting equation (2.41) into the expression for px, u and tau1 solves theproblem entirely. Of particular interest are the formulas for the pressure gradientand the stress on the wall. The formula for is given in in equation (2.34), and32Chapter 2. Swimming with Lubrication Theorythe formula for the stress is given bytau1 = (4U -6)(1 -Y ) - 6(q -1)(1 -Y )2 . (2.43)Now that we have a basic procedure for solving the problem, let us apply itto the viscoelastic problem.2.3 Viscoelastic ProblemWith the Newtonian problem solved fully above, we are now in a position toconsider the lubrication problem with a viscoelastic fluid in the gap between theswimmer and the wall. When substituting a viscoelastic fluid into the problemfor a Newtonian fluid, we have to pick a constitutive model for the viscoelasticfluid. As detailed briefly in Chapter One, this is not necessarily a straightforward task. We use the Oldroyd-8 model because it is the most general weakly-nonlinear constitutive model that includes all possible quadratically nonlinearterms involving the deformation rates allowed by the symmetries of the problem[32]. It is given by,tau + lambda1tautriangleinv + 12?0(trtau)? - 12?1 (tau ? ? + ? ? tau) + 12nu1(tau : ?)delta =?parenleftbigg? + lambda2 ?triangleinv -?2( ? ? ?) + 12(? : ?)deltaparenrightbigg, (2.44)where delta is the two-dimensional identity tensor. Recall the upper convectedderivative of a tensor A,Atriangleinv = partialdiffApartialdifft + u ? partialdiffA -parenleftbig(triangleinvu)T ? A + A ? triangleinvuparenrightbig . (2.45)The substitution of the non-dimensional variables into the Oldroyd-8 con-stitutive model is non-trivial. The leading order equations component wise are:tauxx -parenleftBig2De1 + ?1ch - nu1chparenrightBiguytauxy +parenleftBig2De2 + ?2ch - nu2chparenrightBigu2y = 0,(2.46)tauxy -De1uytauyy + 12parenleftBig?0ch -?1chparenrightBig(tauxx + tauyy)uy -uy = 0,(2.47)tauyy -parenleftBig?1ch -nu1chparenrightBiguytauxy +parenleftBig?2ch -nu2chparenrightBigu2y = 0.(2.48)Note that the time derivatives and advective terms are negligible here, just as33Chapter 2. Swimming with Lubrication Theorythey are in the conservation of momentum equations.The component of the stress tensor that we want to solve for is tauxy, sincethis is the component that appears in the equations of motion. Solving theequations (2.46) through (2.48), we find thattauxy = 1 + alphau2y1 + betau2y uy, (2.49)wherealpha = parenleftbignu2ch - ?2ch parenrightbig parenleftbigDe1 - 12 parenleftbig?0ch - ?1ch parenrightbigparenrightbig + 12 parenleftbig?0ch - ?1ch parenrightbig parenleftbig2De2 + ?2ch - nu2ch parenrightbig,andbeta = parenleftbignu1ch - ?1ch parenrightbig parenleftbigDe1 - 12 parenleftbig?0ch - ?1ch parenrightbigparenrightbig + 12 parenleftbig?0ch - ?1ch parenrightbig parenleftbig2De1 + ?1ch - nu1ch parenrightbig.We are using this model only as an example. Much more general modelscan also be dealt with, so long as we know tauxy(uy), and that the details of theconstitutive model do not interfere with the key simplification of the lubricationtheory, namely that the law reduces to that for a steady-state shear flow for thethin-gap geometry with dominant shear stress.2.3.1 Effective ViscosityObserve that the deviatoric stress equation,tauxy = 1 + alphau2y1 + betau2y uy, (2.50)can be written astauxy = ?(uy;alpha,beta)uy, (2.51)where we can think of?(uy;alpha,beta) = 1 + alphau2y1 + betau2yas being an effective viscosity. In the Newtonian case we have that ?(uy;alpha,beta) =1.We can see in figure (2.2), detailed in table (2.1), that there are five fluid typeregimes, including the Newtonian case, depending on the parameters alpha and beta.Two of these regimes are physically unrealistic because the viscosity diverges andreverses sign for increasing shear rate. We will focus our subsequent analysis on34Chapter 2. Swimming with Lubrication Theorythe more realistic regimes of Newtonian, shear thinning and thickening fluids.180 2 4 6 8 1000.511.522.533.544.55uyEffective viscosityEffective Viscosity versus Shear RateFigure 2.2: Effective Viscosities. Red: Newtonian, Solid Blue: Shear Thin-ning, Dashed Blue: Unrealistic Shear Thinning, Solid Green: Shear Thickening,Dashed Green: Unrealistic Shear ThickeningFluid TypesNewtonian alpha = betaShear Thinning 0 < alpha < betaShear Thickening 0 < beta < alphaUnrealistic Thinning alpha < 0 < betaUnrealistic Thickening beta < 0 < alphaTable 2.1: Fluid types and the corresponding parameter values.The specific parameter values chosen to create figure (2.2) are listed below,in table (2.2). These are also the parameter values that are used in the numericalmethod detailed in the following section. Obviously, any parameter values thatsatisfy the restrictions in table (2.1) could be used. These are just chosen as a18Notice that if we reduce the Oldroyd-8 model to the Oldroyd-B model, we find thatalpha = beta = 0, returning to us the Newtonian problem.35Chapter 2. Swimming with Lubrication Theoryspecific example.Fluid Types Parametersalpha betaNewtonian 1 1Shear Thinning 1 2Shear Thickening 2 1Unrealistic Thinning -1 1Unrealistic Thickening 1 -1Table 2.2: The specific parameter values chosen for the various fluid types.2.3.2 Numerical FormulationIn this section, we will simplify the viscoelastic lubrication equations in a fashionthat allows us to solve them numerically with relative ease. Once the equationsare simplified, we will discretize them and outline the numerical method applied.We will use the method employed in the Newtonian problem as a guideline.Since we cannot solve for u explicitly in the viscoelastic problem, let us firstconsider the continuity equation ux + vy = 0. As we did in the Newtonian case,let us integrate twice to obtainintegraldisplay 1Yudy + Y = q. (2.52)We will exploit the fact thatintegraldisplay 1Yuydy = U, (2.53)and so can write integraldisplay 1Yudy = U -integraldisplay 1Yyuydy,which can be found using integration by parts. We have thenU + Y -integraldisplay 1Yyuydy = q.However, we do not know u explicitly, so let us write it in terms of the deviatoric36Chapter 2. Swimming with Lubrication Theorystress. DefineT(uy) = tau0 + (y -Y )px = 1 + alphau2y1 + betau2y . (2.54)Solving for y givesy = T -tau0 + Y pxpxand so we have that dyduy =Tprimepx . (2.55)For simplification, we will define? = uy, (2.56)and?1 = ?|y=1, ?Y = ?|y=Y . (2.57)Substituting this givesintegraldisplay 1Yyuydy = 1p2xbracketleftBiggintegraldisplay ?gamma1?YTTprime ?d? -(tau0 -Y px)integraldisplay ?gamma1?YTprime ?d?bracketrightBigg. (2.58)Let us callI1(?Y , ?1) =integraldisplay ?gamma1?YTTprime ?d?, (2.59)andI0(?Y , ?1) =integraldisplay ?gamma1?YTprime ?d?. (2.60)Note thatI0(?Y , ?1) = Upx (2.61)Our equations are then as follows:I1(?Y , ?1) -p2x (U + Y -q -UY ) -Utau0px = 0, (2.62)I0(?Y , ?1) -Upx = 0, (2.63)angbracketleftpxangbracketright = 0, (2.64)angbracketlefttau0 -Y pxangbracketright = 0, (2.65)37Chapter 2. Swimming with Lubrication Theorywhere we observe from T(?) = tau0 + (y -Y )px thattau0 = T(?Y ), (2.66)px = T(?gamma1) -T(?gammaY )1 -Y , (2.67)and the unknowns are ?1, ?Y , U, and q.To solve this set of equations, we will discretize and apply Newton?s Method.The full discretized system isI1(?jY , ?j1) -parenleftBiggT(?j1) -T(?jY )1 -Y jparenrightBigg2 parenleftbigU + Y j -q -UY jparenrightbig --UT(?jY )parenleftBiggT(?j1) -T(?jY )1 -Y jparenrightBigg= 0, (2.68)I0(?Y , ?1) -UparenleftBiggT(?j1) -T(?jY )1 -Y jparenrightBigg= 0, (2.69)for j = 1,...,N, andNsummationdisplayj=1bracketleftBiggT(?j1) -T(?jY )1 -Y jbracketrightBiggdeltaxi = 0, (2.70)Nsummationdisplayj=1bracketleftBiggT(?jY ) -Y jparenleftBiggT(?j1) -T(?jY )1 -Y jparenrightBiggbracketrightBiggdeltaxi = 0, (2.71)wheredeltaxi = (xij+1 -xij). (2.72)Thus there are 2N + 2 unknowns in the systemU, q, ?1..N1 , ?1..NY ,from which we can easily calculate the pressure gradient and the stress on theswimmer.Figure (2.3) shows the resulting non-dimensional swimming speeds for afixed propagating wave speed.38Chapter 2. Swimming with Lubrication Theory0 0.2 0.4 0.6 0.8 100.20.40.60.811.21.4Non-Dimensional Swimming Speed a/hUFigure 2.3: Non-Dimensional Swimming Speeds. Red: Newtonian. Blue:Shearing Thinning. Green: Shear Thickening.Solving this problem for a fixed propagating wave speed is not necessarilyilluminating. It is not clear which fluid type is optimal for the swimmer, exceptfor a small range of ah. What we would like to do then, is formulate a differentmethod of comparison that allows us to determine the optimal fluid type for theswimmer for the entire range of values of ah. This method is detailed below.2.3.3 Fixed Dimensional WorkTo obtain a clear result of which fluid is optimal for the swimmer, let us comparethe different realistic fluid regimes when the rate of work per wavelength is fixed.This allows the speed of the propagating wave to vary. We can calculate the rateof work per wavelength of the swimmer, W, by calculating the volume integralofangbracketlefttau : ?angbracketright. (2.73)Non-dimensionally and at leading order, this gives usW =integraldisplay 1Yintegraldisplay 2pi02uytauxydxdy. (2.74)39Chapter 2. Swimming with Lubrication TheoryThis can be rewritten in a way that can utilize the numerics detailed above.That is, we can putW = 1pxintegraldisplay 2pi0integraldisplay ?gamma1?Y2?T(?)Tprime(?)d?dxi. (2.75)What we would like to do now is fix the dimensional work and look at the di-mensional swimming speeds for different fluid types and compare. Particularly,we will be looking at Newtonian, shear thickening and shear thinning fluids.We know that the dimensional swimming speed, let us call it ?, is given by?U0,1,2 = c0,1,2U0,1,2, (2.76)where U is the non-dimensional swimming speed, c is the speed of the propa-gating waves and the superscripts 0, 1, 2, represent the fluid types Newtonian,shear thinning and shear thickening respectively. For a fixed dimensional work,the speed of the propagating waves has to vary in order to compensate. Byputting things into dimensional form we know thatWfixed = 2?c0,1,2hkparenleftBigg2pxintegraldisplay 2pi0integraldisplay ?gamma1?Y?T(?)Tprime(?)d?dxiparenrightBigg0,1,2. (2.77)We can easily solve for the speed of the propagating waves and therefore thedimensional swimming speed.From figure (2.4) we can see that the shear thinning fluid yields the fastestdimensional swimming speed for all values of ah. From figure (2.4) we can alsosee that, contrary to what figure (2.3) might suggest, it is not optimal to swimnear the wall.40Chapter 2. Swimming with Lubrication Theory0 0.2 0.4 0.6 0.8 1012x 10-4 Dimensional Swimming Speeda/hU*cFigure 2.4: Dimensional Swimming Speeds. Red: Newtonian. Blue: ShearingThinning. Green: Shear Thickening.2.3.4 Limiting CasesIn the above analysis, there are two cases that may be of analytical interestwhen thinking of Y = ah sinxi, that when a/h -arrowright1, or the swimmer approachesthe wall and that when a/h -arrowright 0. To calculate these limits we assume that? -arrowrightinfinity and that ? -arrowright0 respectively. Let us consider the former limit first.With this assumption of the leading order component of the strain rateapproaching infinity, we have thattauxy = T(?) similar alphabeta ?. (2.78)This simplifies the integrals I0 and I1 significantlyI0 = beta2alpha parenleftbigtau21 -tau20 parenrightbig , (2.79)I1 = beta3alpha parenleftbigtau31 -tau30 parenrightbig . (2.80)We can then solve the system of equations (3.23) through (3.26) explicitly for41Chapter 2. Swimming with Lubrication TheoryU. We find thatU = 6parenleftbigIasteriskmath1 Iasteriskmath -Iasteriskmath2parenrightbig4Iasteriskmath1 Iasteriskmath3 -3Iasteriskmath22 , (2.81)where Iasteriskmathj are defined above in equation (2.35). We can see that this swimmingspeed is independent of the parameters of the problem, and in fact agrees withthe Newtonian solution. Since we are assuming that the swimmer is approachingthe wall here, let us putY = (1 -epsilon2)sin(xi). (2.82)If we assume that the dominant contribution to the integrals Iasteriskmathj are at the pointalong the swimmer that is nearest the wall, [8], that is at xi = pi/2, then we haveU -arrowright1 as a/h -arrowright1. (2.83)This result was not necessarily expected, due to figure (2.3). However, in figure(2.3), we can see that the general trend is, in fact, towards one. The behaviournear one in the figure can be attributed to a numerical singularity due to thepresence of 1/(1 -Y ).Solving for the remainder of the variables gives thatq = (2U -3)Iasteriskmath13Iasteriskmath2 + 1, (2.84)px = alphabeta 6(U -2)(1 -Y )2 - alphabeta 12(q -1)(1 -Y )3 , (2.85)tau1 = alphabeta 4U -6)(1 -Y ) - alphabeta 6(q -1)(1 -Y )2 . (2.86)Notice that U and q are independent of the parameters of the problem alpha and beta,but px and tau1 are not. Assuming the shape of the swimmer to be that as statedin equation (2.82) gives us to leading orderq = 1, (2.87)px = alphabeta -6(1 -sin(x -t))2 , (2.88)tau1 = alphabeta -2(1 -sin(x -t)) + alphabeta 12(1 -sin(x -t))2 . (2.89)Similarly, we have when a/h -arrowright0, assuming this implies ? -arrowright0, thattauxy = T(?) similar ?. (2.90)42Chapter 2. Swimming with Lubrication TheoryThus the problem simplifies to the Newtonian low amplitude problem. Here,the integrals I0 and I1 simplify toI0 = 12 parenleftbigtau21 -tau20 parenrightbig , (2.91)I1 = 13 parenleftbigtau31 -tau30 parenrightbig . (2.92)The solution to this Newtonian problem can be seen above in equations (2.41),(2.34) and (2.43).However, in this situation it is the case that the amplitude of the swimmeris very small. Let us putY = epsilonsin xi. (2.93)We find that to leading orderU = 3epsilon2. (2.94)This agrees quantitatively with other low amplitude results [37, 22, 25]. This re-sult also agrees nicely with figure (2.3), in which the non-dimensional swimmingspeed appears to be approaching zero quadratically.Solving for the remainder of the variables to leading order gives us thatq = U = 3epsilon2, (2.95)px = 12Y, (2.96)tau1 = 6Y. (2.97)Thus, we can see that these two different limiting cases simplify into problemsthat are solvable analytically. In the latter limit, we found that the problemreduces to the Newtonian problem. In the former limit, we found that theproblem reduced to one similar to the Newtonian, but with a simple dependenceon the viscoelastic parameters, alpha and beta.2.4 Bingham ProblemAs stated above, much more general models can be dealt with in this formulationof the lubrication problem, so long as we know tauxy(uy), and that the detailsof the constitutive model do not interfere with the key simplification of thelubrication theory, namely that the law reduces to that for a steady-state shearflow for the thin-gap geometry with dominant shear stress. Let us consider the43Chapter 2. Swimming with Lubrication Theoryspecific example of a Bingham fluid.In the previous discussion of a Oldroyd fluid, the thin-layer approxima-tion washed away many details of the more complicated viscoelastic rheologicalmodel, leaving only the effect of shear thinning or thickening. Considering aBingham fluid leads to a more interesting lubrication problem from the rheo-logical perspective. As well, the addition of yield stress is biologically relevantto moving organisms [8, 26].A Bingham fluid has the following constitutive lawtauxy = tauysgn(uy) + ?uy, |tauxy| > tauy. (2.98)If |tauxy| <=< tauy then it is the case that uy = 0 and the fluid does not yield.As stated above, the standard set of non-dimensional equations governingthis lubrication problem are:ux + vy = 0, (2.99)tauxy = tau1 -(1 -y)px, (2.100)angbracketleftpxangbracketright = 0, (2.101)angbracketlefttau0 -pxY angbracketright = angbracketlefttau1angbracketright = 0, (2.102)and the shape of the swimmer isY = ah sin(x -t). (2.103)Boundary conditions are unchanged, given by equations (2.24), (2.25), (2.26),and (2.27). If we non-dimensionalize the Bingham constitutive model accordingto the lubrication approximation we obtain,tauxy = Bsgn(uy) + uy |tauxy| > B, (2.104)where B = h?ctauy is the Bingham number. If |tauxy| <=< B then it is the case thatuy = 0 and the fluid does not yield. Note that this implies thattauxy,y = uyy, (2.105)which tells us, recalling that px = tauxy,y,px = uyy. (2.106)44Chapter 2. Swimming with Lubrication TheoryLet us begin solving this by equating our two expressions for tauxy. This yields,tau1 -(1 -y)px = Bsgn(uy) + uy. (2.107)Rearranging to solve for uy gives,uy = tau1 -Bsgn(uy) -(1 -y)px. (2.108)If |tauxy| <=< B then it is the case that uy = 0 and the fluid does not yield. Asis evident by equation (2.108), the flow profile is going to be dependent on thesign of uy. The different values of sgn(uy) represent different regions of the flowprofile. Let us classify these regions as illustrated in table (2.3). We will definey = Y+ as the upper yield surface, and y = Y- as the lower yield surface. It isthe case that the fluid yields if y > Y+ or y < Y-. If Y- < y < Y+, it is the casethat uy = 0 and the fluid moves at a constant speed. Let us call this speed up.It is always true that,Y+ > Y-. (2.109)Outer Regions Inner RegionLocation y < Y- or y > Y+ Y- < y < Y+Stress |tauxy| > B |tauxy(Y?)| = BStrain Rate |uy| > 0 uy = 0Flow Type parabolic plug-likeTable 2.3: Bingham Fluid Flow Profile Regions. In the outer regions, we seeparabolic flow profiles typically characterized by Newtonian flows. In the innerregion, we see what we refer to as pseudo-plug flow. In this region the velocityprofile is constant.Basically, what we expect is a flow profile that looks similar the Newtonianflow profile, see figure (2.5), but with a flattened region around where uy = 0.45Chapter 2. Swimming with Lubrication Theory-0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9-0.500.51uyNewtonian Flow Profile for a = 0.5Figure 2.5: Newtonian Flow Profiles for ah = 0.5. The speed of the swimmer,U, is depicted by the vertical red line.We assume continuity of the flow profile. This implies thatu|(y=Y-) = up = u|(y=Y+) (2.110)where, up = Const, is the speed of the plug-like flow in the inner region. Wecall this region the pseudo-plug. Definesigma = sgn(uy|y>Y+). (2.111)Following immediately from this is,-sigma = sgn(uy|y<Y-). (2.112)As we can see from table (2.3), this implies that,tauxy(Y?) = ?sigmaB (2.113)46Chapter 2. Swimming with Lubrication TheoryApplying this to equation (2.108) and rearranging gives us,tau1 -sigmaB = (1 -Y+)px, (2.114)tau1 + sigmaB = (1 -Y-)px. (2.115)Combining these equations to remove sigma gives,tau1 = px2 (2 -Y+ -Y-), (2.116)and combining them to remove tau1 gives,px = 2sigmaB(Y+ -Y-). (2.117)The fact that, Y+ > Y-, tells us thatsgn(px) = sigma. (2.118)Integrating the continuity equation across the gap, applying the boundaryconditions on v, and then integrating in x givesq = Y +integraldisplay 1Yudy, (2.119)since we know that the conditions on u are independent of x.Let us now concern ourselves with the location of Y+ and Y-. There are fivepossible cases. These are detailed in table (2.4).19Case Location Pressure Pseudo-Plug SpeedI Y- < Y+ < Y < 1 px > 0 up < 0II Y- < Y < Y+ < 1 px > 0 up = 0III Y < Y- < Y+ < 1 No restrictions No restrictionsIV Y < Y- < 1 < Y+ px < 0 up = UV Y- < Y+ < Y < 1 px < 0 up > 0Table 2.4: Location of Yield Surfaces.We consider each one of these cases separately, and obtain different equations19For cases I and V, it is shown belwo that sgn(up) = -sgn(px). Due to the fact that theflow profile in the yielded regions is necessarily parabolic and that we are assuming U > 0,we find that up < 0 and up > 0 for cases I and V respectively.47Chapter 2. Swimming with Lubrication Theoryfor each. These equations are summarized below.20 For case I,Y+ = 12(1 + Y ) - Upx(1 -Y ), (2.120)Y- = Y+ - 2B|px|, (2.121)q = Y + U2 (1 -Y ) - px12(1 -Y )3. (2.122)For case II,Y+ = 1 -radicalBigg2Upx , (2.123)Y- = Y+ - 2B|px|, (2.124)q = Y + px6 (1 -Y+)3. (2.125)For case III,Y+ =parenleftBig1 -Y 2 - 2Upx -4 B|px|parenleftBigY + B|px|parenrightBigparenrightBig2(1 -Y -2 B|px| ) , (2.126)Y- = Y+ - 2B|px|, (2.127)q = Y + px6 parenleftbig(1 -Y+)3 + (Y- -Y )3 -3(1 -Y )(Y- -Y )2parenrightbig .(2.128)For case IV,Y- = Y +radicalBigg-2Upx , (2.129)Y+ = Y- + 2B|px|(2.130)q = Y + px6 (Y- -Y )3 + U(1 -Y ), (2.131)20For the full analysis of each case, the reader is referred to Appendix B.48Chapter 2. Swimming with Lubrication Theoryand for case V,Y- = 12(1 + Y ) - Upx(1 -Y ), (2.132)Y+ = Y- + 2B|px|, (2.133)q = Y + U2 (1 -Y ) - px12(1 -Y )3. (2.134)In each case, we have three equations and five unknowns, px, Y+, Y-, U,and q. But, we still have the constraints,angbracketlefttau1angbracketright = 0, (2.135)angbracketleftpxangbracketright = 0. (2.136)We can now solve the system using Newton?s Iteration. The trick is, at each stepin the iteration, checking which case is relevant, and then using the appropriateequations. This is done by checking the locations of Y+ and Y-, as well as thesign of px (for cases I, II, IV and V). Please refer to the appendices for the fullcode used to solve this problem.Figure (2.6) below shows the pressure gradient and the yield surfaces foundfor ah = 1/2 and B = 1.49Chapter 2. Swimming with Lubrication Theory-1.5 -1 -0.5 0 0.5 1 1.5-4-20246810121416xpx(a) Pressure Gradients-1.5 -1 -0.5 0 0.5 1 1.5-0.4-0.200.20.40.60.81xY+ and Y-(b) Pseudo-PlugFigure 2.6: Pressure Gradients and Yield Surfaces for B = 1 and ah = 1/2.The Newtonian pressure gradient is blue, the Bingham red. In the lower figure,Y+ is blue, Y- is red, and the shaded green represents the pseudo-plug. Thesurface of the swimmer is represented by the black curve and the pink verticalline shows where the pressure gradient changes sign.Figure (2.7) shows the velocity profiles characteristic to each case discussedabove.50Chapter 2. Swimming with Lubrication Theory0 0.1 0.2 0.3 0.4 0.50.40.50.60.70.80.91uy(a) Case I-0.1 0 0.1 0.2 0.3 0.4 0.50.50.60.70.80.91uy(b) Case II0 0.1 0.2 0.3 0.4 0.5 0.600.10.20.30.40.50.60.70.80.91uy(c) Case III0 0.1 0.2 0.3 0.4 0.5 0.60.20.30.40.50.60.70.80.91uy(d) Case IV0 0.1 0.2 0.3 0.4 0.5 0.6 0.70.20.30.40.50.60.70.80.91uy(e) Case VFigure 2.7: Characteristic Bingham Velocity Profiles for B = 1 and ah = 1/2.Blue: Yielded Region. Red: Pseudo-plug Region.As the Bingham number, B = h?ctauy, is increased, the size of the pseudo-plugregion increases until, eventually, it fills up most of the flow domain with viscousboundary layers next to the walls and around the location where the pressuregradient changes sign. As the Bingham number approaches zero, the size of thepseudo-plug region approaches zero and the solution approaches the Newtoniansolution, where Y+ = Y-, see figure (2.8).51Chapter 2. Swimming with Lubrication Theory-1.5 -1 -0.5 0 0.5 1 1.5-10-50510152025303540xpx(a) Pressure Gradients for B = 5-1.5 -1 -0.5 0 0.5 1 1.5-0.4-0.200.20.40.60.81xY+ and Y-(b) Pseudo-Plug for B = 5-1.5 -1 -0.5 0 0.5 1 1.5-2024681012xpx(c) Pressure Gradients for B = 0.1-1.5 -1 -0.5 0 0.5 1 1.5-0.4-0.200.20.40.60.81xY+ and Y-(d) Pseudo-Plug for B = 0.1Figure 2.8: Pressure Gradients and Yield Surfaces for large and small B, ah =1/2. In the left figures, the Newtonian pressure gradient is blue, the Binghamred. In the right figures, Y+ is blue, Y- is red, and the shaded green representsthe pseudo-plugs. The surface of the swimmer is represented by the black curveand the pink vertical line shows where the pressure gradient changes sign.Figure (2.9) shows the how the swimming speed and mass flux vary withthe Bingham number. The mass flux is a monotonically decreasing function ofthe Bingham number, which is expected due to the fact that the yielded regiondecreases in size as the yield stress increases. Whereas, the swimming speed isa non-monotonic function of the Bingham number with a clear maximum value,located at B similarequal 1. This suggests that is it beneficial for an organism to swim ina fluid that has a moderate yield stress.2121It has been observed that snails and slugs move over a fluid that exhibits a yield stress[8].52Chapter 2. Swimming with Lubrication Theory0 1 2 3 4 50.4850.490.4950.50.5050.510.515BU(a) U versus B for ah = 1/20 1 2 3 4 50.460.4650.470.4750.480.4850.490.4950.50.505Bq(b) q versus B for ah = 1/2Figure 2.9: Swimming speed and mass flux for varying Bingham number. Thecorresponding Newtonian values are shown in red.2.5 DiscussionIn this chapter we presented the problem of a lubricated swimmer by a wall.We formulated the problem initially with a Newtonian fluid, and returned theresults of previous work done by others [22]. Next, we considered the problemwith a viscoelastic fluid in the gap.We found that the speed of low amplitude swimmers close to walls is signif-icantly larger than those away from walls. This is also observed experimentally[20]. We also found that a shear thinning fluid yields the optimal speed forthe swimmer given a fixed rate of work when compared to a swimmer in aNewtonian fluid and a shear thickening fluid. It is the case that viscoelasticitycan decrease or increase the swimming speed compared to the Newtonian case,depending on the properties of the specific viscoelastic fluid in question.Motivated by the results shown in figure (2.3), we took the analytical limitas the swimmer approaches the wall. We found that the non-dimensional swim-ming speed approaches the value of one, regardless of the fluid it is in.22 Inthe opposite limiting case, when the ratio of the amplitude of the swimmer onthe gap width goes to zero, we have the Newtonian problem and its relevant22This contradicts the limiting case of the introductory problem, where the limit was foundto depend on the parameters of the problem. This is not of concern due to the fact that, inthe introductory problem, no thought was given to what the strain rates were doing in thelimit. The introductory problem was just a preliminary example used to motivate the use oflubrication theory. In fact, if we reduce the Oldroyd-8 constitutive model to the Oldroyd-Bmodel above, the Newtonian problem is returned, hence no contradiction when the limit istaken properly.53Chapter 2. Swimming with Lubrication Theoryresults returned. That is, we found that the non-dimensional swimming speedapproaches zero quadratically in a.Due to the formulation of this problem, theoretically it should be easy toconsider different constitutive models and wave shapes, thus allowing the rangeof applications of this problem to broaden.For example, perhaps we would like to consider a different swimmer shapein Chapter Two. Let us consider the toy example of a swimmer shape givenby Y = tanh(12sin xi), as shown in figure (2.10a). The analysis of the problemdoes not change. The resulting swimming speeds of such a swimmer are shownin figure (2.10b) and were found using the same code as was used for the resultsfound above. It is evident that different swimmer shapes can yield differentresults, depending on the application one has in mind. The mathematical setup we chose here lends itself nicely to generality, and exploring these potentialdifferent applications it certainly something that should be looked into further.0 1 2 3 4 5 6-1-0.8-0.6-0.4-0.200.20.40.60.81xiYSquare Wave Swimmer(a) Wave Shape0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-0.500.511.522.53 x 10-4 Dimensional Swimming Speeda/hU*c(b) Swimming SpeedsFigure 2.10: Example of alternate wave shape and resulting swimming speeds.Red: Newtonian. Blue: Shearing Thinning. Green: Shear Thickening.In the spirit of broadening the applications of this problem, we then consid-ered the problem with a Bingham fluid. We found that in order to solve theproblem, we had to consider five different possible cases that depended on thelocation of the yield surfaces. For a fixed swimmer amplitude, ah = 1/2, wesolved for the pressure gradient and the location of the yield surfaces. At thisfixed amplitude, we found that it is potentially advantageous for an organismto swim in a Bingham fluid with a moderate yield stress, as it results in a fasterswimming speed than a Newtonian fluid.54Chapter 3Elastic Swimming withLubrication TheoryIn the previous chapter we looked at the swimming problem from the perspectiveof the swimmer having a known motion or shape. We then analyzed the resultingswimming speed with the swimmer subjected to a fixed dimensional work orrate of work per unit wavelength. There is, however, a different approach onecould use to look at this problem. Rather than a prescribed shape, we canconsider a prescribed force distribution along the length of the swimmer. In thissituation we have to think more carefully about the swimmer itself. Previously,we only were really concerned with the swimming speed. Now, we must thinkabout how this applied force affects the swimmer. What kind of propertiesdoes this swimmer have to govern how it will respond to this forcing? Thisis a question that has a large number of potential answers. The specifics offlagellum movement are very complex, and so if we were to model the swimmeras a realistic flagella the problem would be very involved [23, 22]. However, herewe are more interested in the fluid dynamics and kinematics of the problem,and so the answer to the question posed above is that we will assume that theswimmer responds to the applied force as an elastic beam. This is a simpler,more practical choice than trying to incorporate all the intricacies of a realisticflagellum. In fact, this assumption makes our problem more general with awider range of potential applications [2, 41, 8, 39].Potentially, we could have considered a different mechanism of wave prop-agation along the swimmer. Rather than having forcing distributed along theentire length of the swimmer, we could consider the situation of a passive elasticswimmer with one fixed end that moves in a vertical fashion thus propagatingwaves down the length of the swimmer. This postulation is interesting from amathematical perspective, and perhaps has some interesting alternate applica-tions, but it has been shown that observed wave forms of spermatozoa cannot55Chapter 3. Elastic Swimming with Lubrication Theorybe accounted for with this mechanism [29].Formally then, this is the problem of an elastic swimmer subject to a knownforcing along its length in a lubrication setting. Its shape and swimming speedare unknown.Our goal in this chapter is to discover how an elastic swimmer reacts whensubjected to a known distribution of applied torques whilst in a viscoelasticfluid. However, the problem of how such a swimmer acts in a Newtonian fluidproves to be interesting itself.3.1 Mathematical FormulationLet us consider the following: a thin, isotropic elastic beam in two-dimensionsimmersed in a fluid very close to a rigid wall. This beam is then subjected toknown internal forces distributed along its length, f(x,t) = asin k(x -ct) thatcause it to bend into an unknown shape, Y . This shape then produces a thrust,propelling it through the fluid at an unknown, constant speed U. For the sake ofsimplicity, we will assume that the applied force is sinusoidal. We use a similarset up for our problem as was used in [2], where the forces applied are referredto as actively distributed torques.Figure 3.1: Geometry of Lubrication Problem. This time, the swimmer shapein unknown as well as the swimming speedWhat we want to discover here is what the unknown swimmer shape is andthe swimming speed it produces.To begin with, let us consider the boundary conditions. In the frame movingwith the swimmer, our boundary conditions remain unchanged for this problemfrom the last. That is, non-dimensionally they are for the horizontal velocity56Chapter 3. Elastic Swimming with Lubrication Theorycomponent,u(x,h) = U, (3.1)u(x,Y ) = 0, (3.2)and for the vertical velocity component v we havev(x,h) = 0, (3.3)v(x,Y ) = -Yx, (3.4)recalling that Yt = -Yx. It is interesting to note here that if we remained inthe original frame with the wall fixed, the boundary condition on the swimmerwould actually be a kinematic or free boundary condition due to the fact thatthe swimmer is now considered elastic. But, in the frame moving with theswimmer there is no horizontal velocity associated with the beam, therefore thekinematic condition reduces to the one stated above.Now, let us consider the fluid. It can be modelled as was done previously,applying the general assumptions detailed in Chapter One, bytriangleinvp = triangleinv? tau, (3.5)triangleinv? u = 0, (3.6)and general constitutive modelf(tau, ?) = 0. (3.7)Here we have an elastic beam, thus we must include another equation for thevertical force balance. This is the beam equation including an internally appliedtorque [2, 28]. In full dimensional terms this equation is, neglecting gravity,0 = -partialdiffxx (f + DpartialdiffxxY ) + p, (3.8)where D is the bending stiffness of the beam, and p, Y , f the pressure, beamshape and applied forcing respectively. We enforce periodicity in the swimmerup to and including its third derivative. Mathematically it is required for thehigher derivatives due to the fact that, once the beam equation is incorporated,we end up with a fourth order ordinary differential equation, thus requiringthe four conditions. Physically, only this condition and the continuity of the57Chapter 3. Elastic Swimming with Lubrication Theoryswimmer are required.23 In truth, we are interested in the pressure gradientonly, thus we take the derivative of the above equation and usepx = partialdiffxxx (f + DpartialdiffxxY ). (3.9)We also have the constraintsangbracketleftpxangbracketright = 0, (3.10)angbracketlefttau1angbracketright = 0, (3.11)angbracketleftY angbracketright = 0. (3.12)When looking at equation (3.9), we can see that in order for angbracketleftpxangbracketright = 0 to holdwe require that angbracketleftYxxxxxangbracketright = 0. This is the same thing as saying Y is periodic inits fourth derivative. This requirement is due to the choice of f(x,t). Thus wecan enforce the constraint on the pressure gradient by enforcing periodicity onthe swimmer shape up to and including its fourth derivative, and from here onout not include the constraint specifically. Observe that the condition angbracketleftpxangbracketright = 0is also utilized in the formulation of the zero force on swimmer condition whenit is put into the form angbracketlefttau1angbracketright = 0, as noted in Chapter One.Non-dimensionalizing as we did in Chapter Two, using the following non-dimensional variables,t = 1omega ?, ? = omega?,x = 1k ?, y = h?,u = c?, ? = epsilonh2k4?c D(p,tau) = ?ch (1epsilon ?, ?), ? = epsilonk2h?c aapplying the lubrication approximation and simplifying algebraically yields,dropping hats,I1(?Y , ?1) -p2x (U + Y -q -UY ) -Utau0px = 0, (3.13)I0(?Y , ?1) -Upx = 0, (3.14)angbracketleftY angbracketright = 0, (3.15)angbracketlefttau1angbracketright = 0. (3.16)23This condition?s specific use was not required previously because the shape of the swimmerwas assumed.58Chapter 3. Elastic Swimming with Lubrication TheoryRecall thatI0 =integraldisplay ?gamma1?YTprime ?d?, (3.17)I1 = TTprime ?d?, (3.18)T = tauxy. (3.19)Similarly then, we also utilize the fact thatpx = T(?gamma1) -T(?gammaY )1 -Y , (3.20)to remove px from the above equations. As well, recall thattau1 = T(?1). (3.21)Putting the beam equation into non-dimensional terms yieldspx = partialdiffxxx (f + DpartialdiffxxY ). (3.22)So, our unknowns are now ?1, ?Y , U, and q and the new unknown, theswimmer shape, Y . Summarizing the full set of non-dimensional equationsrelevant to solving the problem gives usI1(?Y , ?1) -p2x (U + Y -q -UY ) -Utau0px = 0, (3.23)I0(?Y , ?1) -Upx = 0, (3.24)fparenleftBig?, ?parenrightBig= 0, (3.25)px -partialdiffxxx (f + DpartialdiffxxY ) = 0, (3.26)angbracketleftY angbracketright = 0, (3.27)angbracketlefttau1angbracketright = 0, (3.28)along with periodic conditions on the derivatives of Y up to and including thefourth derivative which can be written asangbracketleftYxangbracketright = angbracketleftYxxangbracketright = angbracketleftYxxxangbracketright = angbracketleftYxxxxangbracketright = angbracketleftYxxxxxangbracketright = 0 (3.29)These equations are very similar to what we had before with only the additionof the beam equation and the periodicity conditions on the swimmer and itsfirst five derivatives. Again, the conditions on the swimmer were true before59Chapter 3. Elastic Swimming with Lubrication Theorysince we were imposing the shape of the swimmer, but now we must rigorouslyenforce it.As we have seen above, the layout of this problem and the last one has thefluid properties included simply in the form of the equation tauxy = T(?). This isan elegant and general way of looking at the problem and allows us to look atdifferent constitutive equations with relative ease. However, the complexity ofthe constitutive equation can still heavily affect the ease of solving the systemboth analytically and numerically.As done above, let us first consider the Newtonian problem in order to gainsome insight into the problem before considering more complex fluid models.3.2 Newtonian ProblemSince, with the exception of the inclusion of the beam equation and the con-straints on the swimmer, the equations are the same for this problem as theyare for the Newtonian problem for a swimmer of prescribed shape, we can skipa few steps. We have already solved equations (3.23), (3.24), and (3.28) for aNewtonian fluid in the previous chapter. However, since Y is unknown, we canonly utilize the equations for the pressure gradient and stress on the wall, asthe calculation of the swimming speed and mass flux relied on calculating theintegrals Iasteriskmathj . Not all is lost though, we can use still use the following equationsfor px and tau1px = 6(U -2)(1 -Y )2 - 12(q -1)(1 -Y )3 , (3.30)tau1 = (4U -6)(1 -Y ) - 6(q -1)(1 -Y )2 . (3.31)These are most useful since the pressure gradient and unknown swimmer shapeappear in the beam equation and the stress on the swimmer appears in a con-straint. First, let us substitute the algebraic equation for px, equation (3.30),into the beam equation, equation (3.26). This gives uspartialdiffxxx (asin(x -t) + DpartialdiffxxY ) = 6(U -2)(1 -Y )2 - 12(q -1)(1 -Y )3 . (3.32)If we introduce a new variable,zeta = (1 -Y ), (3.33)60Chapter 3. Elastic Swimming with Lubrication Theoryand utilize the dependence on xi = (x-t) in the problem, we obtain the equationzeta(5) = - 1Dparenleftbiggacos xi + 6(U -2)zeta2 - 12(q -1)zeta3parenrightbigg, (3.34)where derivatives are with respect to xi. We can then solve for zeta and U and qby applying the constraints. To do this, we must first put the constraints intoterms of zeta. The first constraint, < Y >= 0 is simple. It becomesangbracketleftzetaangbracketright = 1. (3.35)The second, angbracketlefttau1angbracketright = 0 is less straight forward. However, due to our abovealgebraic manipulation, we already know what tau1 is in terms of Y , U, and q.Thus it is a simple matter of putting this in terms of zetatau1 = (4U -6)zeta - 6(q -1)zeta2 , (3.36)and so the constraint becomes(4U -6)angbracketleftzeta-1angbracketright -6(q -1)angbracketleftzeta-2angbracketright = 0. (3.37)To solve this analytically requires solving a cubic, so rather than doing that, letus convert the problem to something which we can apply numerics.3.2.1 Numerical FormulationWe can think of equation (3.34) as an ordinary differential equation with eigen-values U and q and constraints given by equations (3.35) and (3.37). To solvethis numerically, let us define a new variable r such that rprime = xi, thus giving usr =integraldisplay xi0zetadxi. (3.38)Similarly define the variable s so that sprime = parenleftbig(4U -6)xi-1 -6(q -1)xi-2parenrightbig, whichgives uss =integraldisplay xi0parenleftbig(4U -6)xi-1 -6(q -1)xi-2parenrightbig dxi. (3.39)We have then thatr(0) = 0, r(2pi) = 2pi, (3.40)61Chapter 3. Elastic Swimming with Lubrication Theoryands(0) = 0, s(2pi) = 0, (3.41)as conditions on these new variables.Thus, our system of equations to solve iszeta(5) = - 1Dparenleftbiggacos xi + 6(U -2)zeta2 - 12(q -1)zeta3parenrightbigg, (3.42)r =integraldisplay xi0zetadxi, (3.43)s =integraldisplay xi0parenleftbig(4U -6)xi-1 -6(q -1)xi-2parenrightbig dxi, (3.44)where we think of U and q as eigenvalues, with conditionszeta(0) = zeta(2pi), zetaprime(0) = zetaprime(2pi), zetaprimeprime(0) = zetaprimeprime(2pi),zeta(3)(0) = zeta(3)(2pi), zeta(4)(0) = zeta(4)(2pi), (3.45)andr(0) = 0, r(2pi) = 2pi, (3.46)s(0) = 0, s(2pi) = 0, (3.47)for unknowns zeta, zetaprime, zetaprimeprime, zeta(3), zeta(4), zeta(5), r and s with eigenvalues U and q. Todo this we use the MATLAB solver bvp4c. To see the code please refer toappendices.Figure (3.2) shows the resulting shape of the swimmer as the amplitude ofthe forcing increases. It shows us that as the amplitude gets larger the swimmershape develops a plateau region followed by a sharp spike region. By lookingat the phase portraits of this large amplitude solution shown in figure (3.3) wecan see that there exists a inner region that reminiscent of a fixed point. Thisexplains the section of the swimmer shape that is nearly constant that existsnear xi similar pi/2.62Chapter 3. Elastic Swimming with Lubrication Theory0 1 2 3 4 5 6-2.5-2-1.5-1-0.500.5Swimmer Shape for D=0.2 and a=0..25xi=x-tY5 10 15 20 250.20.40.60.8a(a) Newtonian Solution for D = 0.20 1 2 3 4 5 6-2.5-2-1.5-1-0.500.5Swimmer Shape for D=1.5 and a=100xi=x-tY20 40 60 80 1000.20.40.60.8a(b) Newtonian Solution for D = 1.5Figure 3.2: Newtonian Solution of Elastic Lubrication Problem. In the bottomleft hand corner of each figure we see the swimming speed (red) and the massflux (green) versus the forcing amplitude. The magenta horizontal lines are thecalculated plateau values for each case.-3-2-101-4-2024-10010Y `Phase Space for D = 0.2YY ``(a) Phase Portrait for D=0.2-3-2-101-4-2024-10010Y `Phase Space for D = 1.5YY``(b) Phase Portrait for D=1.5Figure 3.3: Phase portraits for Newtonian Solution of Elastic Lubrication Prob-lem with large amplitude forcing.We would like to get an approximation of this plateau value. However, itis not immediately obvious how we may scale things in such a way that givesus a rigorous asymptotic solution to this inner region. Let us assume that weare inside this inner region for zeta, near xi = pi/2, and expand in the standardasymptotic fashionzeta = zetap + zeta1 + ... (3.48)where zetap is the plateau constant and |zeta1| << 1. Substituting this into equation63Chapter 3. Elastic Swimming with Lubrication Theory(3.34), and rearranging slightly we obtainDad5dxi5 zeta1 = -cosparenleftBigpi2 + sigmaparenrightBig- 6(U -2)a(zetap + zeta1)+ 12(q -1)a(zetap + zeta1), (3.49)where we are assuming that sigma is some small number. It is easy to see then thatat leading order we obtainzetap = 2(q -1)U -2 . (3.50)To solve this fully analytically is certainly a challenge, but if we take this formulafor the plateau value, recalling that zeta = 1 -Y , and plot it using our numericalvalues for U and q we obtain a very remarkable fit, as can be seen as the magentaline in figure (3.2). Note that the plateau value is independent of the bendingstiffness. These results suggest that there are some truly interesting dynamicsat work in this large amplitude regime which most certainly should be analyzedat greater depth, including a full analysis of the outer region.We can also see from figure (3.2b) that the plateau approaches y = 1. Thisobservation is substantiated by the formulaYp = 1 - 2(q -1)U -2 , (3.51)if we observe that q arrowright1 as the amplitude gets large. As well, we can note fromfigure (3.2) that there is an obvious optimal forcing amplitude for the swimmingspeed and after that, as the swimmer gets too close to the wall and the stressesget large, the speed decreases.3.2.2 Low Amplitude SolutionLet us consider the situation when the amplitude of the applied forcing or torqueon the swimmer is very small. That is, puta lessmuch1, (3.52)where the forcing is given by f = asin(x - t). We can then assume that thisimplies that the amplitude of the swimmer shape is very small|Y | lessmuch1. (3.53)64Chapter 3. Elastic Swimming with Lubrication TheoryAt this point we can then notice that this is a similar situation to the one dis-cussed in Chapter Two. Therefore we can use the low amplitude approximationsfor px and tau1 as we found in Chapter Two. We found initially thatpx = 6(U -2)(1 -Y )2 - 12(q -1)(1 -Y )3 , (3.54)tau1 = (4U -6)(1 -Y ) - 6(q -1)(1 -Y )2 . (3.55)Except here, Y is unknown. However, this does not stop us from being ableto get a low amplitude approximation for these equations. We obtain the sameleading order approximations as we did beforepx = 12Y, (3.56)tau1 = 6Y. (3.57)Notice that the constraint on tau1 holds automatically if we enforce the constrainton Y . Applying the low amplitude approximation to equation (3.9)24 and sub-stituting into it equation (3.56) gives, after some simple rearranging,DYxxxxx -12Y = fxxx, (3.58)with constraintangbracketleftY angbracketright = 0, (3.59)and periodicity in Y up to and including its fourth derivative. Solving this wefind that the homogeneous solutions must be ignored due to the fact that theydo not satisfy the periodicity conditions. Thus it is only the particular solutionthat is relevant. We find the particular solution to beY = aD2 + 122 (-12cos xi + D sinxi) . (3.60)From this low amplitude solution we can see that for a moderate value of Dall terms will balance, but for large or small values of D we will have differentterms balancing. For small amplitude forcing and a small bending stiffness we24Note that the low amplitude approximation does not change the beam equation since allterms in it become small and therefore all balance65Chapter 3. Elastic Swimming with Lubrication Theorywould expect the shape of the swimmer to look likeY = -12aD2 + 122 cos xi (3.61)and for small amplitude forcing and large bending stiffness we expect the shapeof the swimmer to look likeY = aDD2 + 122 sin xi (3.62)3.3 Viscoelastic ResultsLet us now tackle the elastic lubrication problem with a viscoelastic fluid withinthe gap. Here again, we choose to use the Oldroyd-8 constitutive model. Non-dimensionalizing the constitutive model, simplifying, and solving for tauxy, as wedid in the previous chapter, yieldstauxy = T(?) = (1 + alpha?gamma2)(1 + beta ?2) ?gamma. (3.63)To solve this numerically, we employ a similar method as was used in ChapterTwo. That is, we discretize the system of equations and apply Newton?s Method.The system of discretized equations is as follows:I0(?jY , ?j1) -Upjx = 0, (3.64)I1(?jY , ?j1) -pj2x parenleftbigU + Y j -q -UY jparenrightbig -UT(?jY )pjx = 0, (3.65)pjx -fjxxx -DY jxxxxx = 0, (3.66)?Nj=1T(?j1) = 0, (3.67)?Nj=1Y j = 0. (3.68)There are 3N + 2 unknowns, ?1..NY , ?1..N1 , Y 1..N, U, and q. And we havethat, pjx = T(?gammaj1)-T(?gammajY )1-Y j . Note that fxxx is given, as f is given. The Yxxxxxterm is found using a difference formula. Please refer to appendices for the fullnumerical method used to solve this problem.Figure (3.4) shows the swimmer shape varying for changing ah values and theresulting swimming speeds. We see that for low bending stiffness, the differencesbetween the three fluid types is noticeable. Particularly the swimming speed,where you can see that a shear thinning fluid yields the fastest speed and a shear66Chapter 3. Elastic Swimming with Lubrication Theorythickening fluid yields the slowest. But, as we increase the bending stiffness, asseen in figures (3.5) and (3.6), both the swimmer?s shape and speed becomesless differentiated between the different fluid types. Notice the asymmetry herefor certain values of the bending stiffness and the similarities between the shapeof the swimmer for Newtonian, shear thinning, shear thickening.2 4 6-0.1-0.0500.05alpha=beta=1 and D=0.2xY0.2 0.4 0.6 0.8 1 1.20.010.020.030.04aUU vs a for D=0.22 4 6-0.1-0.0500.050.1alpha=1, beta=2, D=0.2xY2 4 6-0.1-0.0500.05alpha=2, beta=1, D=0.2xYFigure 3.4: Solution for viscoelastic elastic lubrication problem for various fluidtypes with D=0.2. Red: Newtonian. Blue: Shearing Thinning. Green: ShearThickening.67Chapter 3. Elastic Swimming with Lubrication Theory2 4 6-0.0500.05alpha=beta=1 and D=7xY0.2 0.4 0.6 0.8 1 1.20.0050.010.0150.020.025aUU vs a for D=72 4 6-0.0500.05alpha=1, beta=2, D=7xY2 4 6-0.0500.05alpha=2, beta=1, D=7xYFigure 3.5: Solution for viscoelastic elastic lubrication problem for variousfluid types with D=7. Red: Newtonian. Blue: Shearing Thinning. Green:Shear Thickening.68Chapter 3. Elastic Swimming with Lubrication Theory2 4 6-0.04-0.0200.020.04alpha=beta=1 and D=27xY0.2 0.4 0.6 0.8 1 1.212345x 10-3aUU vs a for D=272 4 6-0.04-0.0200.020.04alpha=1, beta=2, D=27xY2 4 6-0.04-0.0200.020.04alpha=2, beta=1, D=27xYFigure 3.6: Solution for viscoelastic elastic lubrication problem for various fluidtypes with D=27. Red: Newtonian. Blue: Shearing Thinning. Green: ShearThickening.3.3.1 Low Amplitude SolutionLet us again consider the situation when the amplitude of the applied forcingor torque on the swimmer is very small. Above we assumed that the smallforcing amplitude yields a small swimmer amplitude. Here, we will assume thesame, but also assume that this very small forcing and resulting small swimmeramplitude implies a very small shear rate| ?| lessmuch1. (3.69)This small shear rate affects the form of the constitutive equation for the stress.In particular, we have that the constitutive equation becomesT(?) = ?, (3.70)reducing the problem to the Newtonian one. Therefore, the solution foundabove for the low amplitude Newtonian problem is also the solution for the low69Chapter 3. Elastic Swimming with Lubrication Theoryamplitude viscoelastic problem. Recall that we foundY = aD2 + 122 (-12cos xi + D sinxi) . (3.71)Therefore, we can see that for small values of the bending stiffness the swim-mer will look like a negative cosine curve, and for large values of the bendingstiffness it will look like a sine curve. As we can see from figure (3.7) theseanalytical results agree very nicely with the numerical viscoelastic results forsmall amplitude forcing.1 2 3 4 5 6-505x 10-4 alpha=beta=1 and D=1xY1 2 3 4 5 6-505x 10-4 alpha=1, beta=2, D=1xY1 2 3 4 5 6-505x 10-4 alpha=2, beta=1, D=1xYFigure 3.7: Comparison between analytical and numerical results for low am-plitude swimming problem. Numerical results are solid lines, analytical resultsare shown by *. Red: Newtonian. Blue: Shearing Thinning. Green: ShearThickening.3.4 DiscussionIn this chapter we presented the problem of a lubricated elastic swimmer by awall. We considered first the Newtonian problem and discovered that it yieldsvery interesting results. For large amplitude forcing we found that a plateauregion develops along the swimmer. This plateau approaches the wall as the70Chapter 3. Elastic Swimming with Lubrication Theoryamplitude of the forcing becomes large. We referred to this region as the innerregion as it corresponds to a region in the phase portrait that is reminiscent ofa fixed point. Upon expanding around this inner region, we found a value ofthis plateau analytically that agreed very well with numerical results. Follow-ing this plateau region is a large spike type region. Unfortunately, solving forthis region analytically is beyond the scope of this thesis. These results are notentirely irrelevant to the physical situation under consideration as large ampli-tude, asymmetrical beat patterns are observed in hyperactivated spermatozoa[35], [17], [36], [34]. These results are also extremely interesting from a math-ematical perspective and it is the hope of this author that at some point, inthe not so distant future, they will be analyzed in greater depth. As well, wesaw that there is a definite forcing amplitude that yields an optimal swimmingspeed. The shape of the optimal swimmer at this amplitude can be seen infigure (3.8).0 1 2 3 4 5 6-1.4-1.2-1-0.8-0.6-0.4-0.200.20.40.6Swimmer Shape for D=0.2 and a=0..6xi=x-tYFigure 3.8: Optimal swimmer shape for D=0.2.As we can see in figure (3.8), this optimal shape is neither a true sinusoid,nor does have the flattened, plateau region characteristic of the large amplitudelimit.We next calculated the low amplitude analytical solution of the Newtonian71Chapter 3. Elastic Swimming with Lubrication Theoryproblem. We found that both the pressure gradient and the stress on the wallare proportional to the swimmer shape, and that the swimmer shape isY = aD2 + 122 (-12cos xi + D sinxi) , (3.72)a solution that is valid for any value of the bending stiffness D. This resultwas returned to us when we calculated the low amplitude solution for a elasticswimmer in a viscoelastic fluid. It was found to agree with the numerical resultsfor the viscoelastic problem extremely well.We then took the next logical step and considered a viscoelastic fluid. Wefound that for moderate bending stiffness, a shear thinning fluid yielded thefastest swimming speed and that as the bending stiffness was increased, theshear thinning and shear thickening solutions collapsed towards the Newtoniansolution. When considering this fact and the Newtonian solution, it seems thatthe most interesting parameter regime for problem as a whole is moderate orsmall bending stiffness coupled with a large amplitude forcing.It would be most desirable for us to find a numerical or analytical solution forthe viscoelastic problem with a large amplitude forcing. This is not straightfor-ward, as the code that currently solves the viscoelastic problem cannot handlelarge amplitude forcing. Finding an analytical solution is also not straightfor-ward. To see why this is first note that from the Newtonian problem we havethatpx = uyy, (3.73)giving us thatuy = tau1 -(1 -y)px. (3.74)We can see from figure (3.9) that the strain rates in the Newtonian problem aregetting larger, but it is not obvious that we can make the blanket assumptionthat they get large enough to simplify the constitutive relation.25 Even if theconstitutive equation can be simplified, we would assume that the solution stillhas dependence on the viscoelastic parameters alpha and beta, as did the solution forthe large amplitude swimmer in chapter two.25Please recall that in the Newtonian problem, strain rates and stresses are equal non-dimensionally, tau = ?72Chapter 3. Elastic Swimming with Lubrication Theory0 1 2 3 4 5 6-40-30-20-1001020304050xi=x-ttau1: Black, tau0: CyanStresses\Strain RatesFigure 3.9: Strain rates on the swimmer (cyan) and on the wall (black) fora=100 and D=1.5. Note that the largest strain rates occur in the region wherethe plateau ends and the spike begins.73Chapter 4ConclusionsThis thesis has presented an in depth discussion of swimmers near walls. Severaldifferent problems were analyzed and within these problems, different fluid typeswere considered. The mathematical setup of these problems can be applied toa broad spectrum of applications as it easily adapts to many similar situations.For example, similar problems that have different boundary conditions, walllocations, shapes, or constitutive models all require only simple changes to thebasic formulation.26 Viscoelastic lubrication problems arise in many places, suchas physiological and industrial applications, as well as other types of swimmingproblems [30, 1, 40, 41].4.1 Summary of ResultsTo begin with, we motivated the lubrication problem by looking at the prelim-inary problem of a swimmer in the presence of a wall. We then considered theproblem of a lubricated swimmer with a prescribed shape. It was concludedthat shear thinning fluids are optimal for swimmers in this situation when com-pared to shear thickening and Newtonian fluids. This was found by consideringthe problem with a fixed rate of work and looking at the resulting dimensionalswimming speeds. The two limiting cases of the distance between the wall andthe swimmer approaching zero and infinity were considered analytically. It wasfound that the Newtonian swimming speed was returned in both situations. Inthe former limit, the Newtonian problem was returned entirely. In the latterlimit, the solution, aside from the swimming speed and mass flux, was found todepend on the viscoelastic parameters alpha and beta in a simple fashion.We then considered the lubrication problem with a swimmer of prescribedshape in a Bingham fluid. Here we fixed the amplitude of the swimmer and26There are some restrictions to the scope of these differences, for example with changes inthe constitutive model, we have to know tauxy(uy), and that the law reduces to a steady-stateshear flow for the thin-gap geometry with dominant shear stress.74Chapter 4. Conclusionssolved for the location of the yield surfaces in the problem. The yield surfaces,y = Y+ and y = Y-, were defined to be the upper and lower bounds of theregion of the fluid where |tauxy| <=< tauy. We found that there exist regions alongthe swimmer where the plug-flow reaches the swimmer (case II), reaches thewall (case IV), and regions where the fluid is entirely yielded (cases I and V).When varying the Bingham number, it was found that there exists a maximumswimming speed at B similarequal 1, which is greater than the Newtonian swimmingspeed.Combining the results from Chapter Two, we may speculate that it is themost advantageous for a swimmer to be moving through a fluid that exhibits ayields stress, and once yielded, behaves in the manner of a shear thinning fluid.An example of such a fluid is a Herschel-Bulkley fluid with n < 1.The problem of an elastic lubricated swimmer was then considered. Theswimmer was modelled with the beam equation with an addition of an internalapplied torque as was done in [2]. When looking at the problem with a New-tonian fluid, it was found that there exists an amplitude of forcing that yieldsan optimal swimming speed. As the amplitude of the forcing was increased,it was found that the swimmer shape began to develop a plateau region thatapproached the wall, followed by a large downward spike region. The value ofthe plateau region was successfully calculated analytically. We then introduceda viscoelastic fluid. It was again found that a shear thinning fluid is optimal fora swimmer as it yields the fastest swimming speeds. As the bending stiffnessof the swimmer was increased, it was observed that the problem collapsed intothe Newtonian problem.In general, we found that a shear thinning fluid or a fluid that exhibits ayield stress is theoretically optimal for swimmers near walls. Supporting thistheory is the fact that polymeric fluids, the type of fluid found in most biologicalsettings, are not often shear thickening and often exhibit a yield stress [5].In the case of the elastic swimmer, we found that it is not optimal for aswimmer to continue increasing its energy output, as it does not continue toincrease the swimmer?s speed. This is because as the forcing along the swimmerincreases, so does the stress within the fluid and so there exists a point in whichthe stresses in the fluid are so great that they begin to hamper movement. Thissuggests that in nature, swimmers do not exert themselves to their maximum,but rather moderately in such a way that maximizes speed.75Chapter 4. Conclusions4.2 Further WorkThe results of this thesis suggest further work. Due to the formulation of theproblem, and in particular, the simplification of the constitutive model whennon-dimensionalized, it is theoretically straightforward to consider numerousdifferent constitutive equations and fluid types. This suggests developing amore robust numerical method of solution.From a biological perspective, there are many extensions of the analysis thatcan be done. For example, it is the case that in a biological setting, swimmersare of finite length. Thus, it may be of interest to consider the problems with afinite swimmer and relax the periodicity conditions.There are numerous extensions that can be done with the elastic swimmerproblem. Due to the interesting results of the Newtonian elastic swimmer prob-lem, this author suggests a full analysis of the dynamics of this problem. It isunfortunate that such an analysis is beyond the scope of this thesis. As well,due to the large amplitude of the spike region encountered in this problem, onemight consider the swimmer in a channel, therefore restricting the amplitude ofthe spike. There is potential for things to get quite interesting in this case.One might also be interested in looking into the large amplitude forcingviscoelastic problem, as it is not entirely clear that this will reduce into theNewtonian problem as it did in the similar case considered in Chapter Two.As mentioned throughout this thesis, we are by no means limited to thespecific choices we have made. This is true for the choice of constitutive models,swimmer shapes (Chapter Two) and applied torques (Chapter Three). Themathematical set up we chose here lends itself nicely to generality, and exploringthese potential alternate applications, via different choices of the constitutivemodel, swimmer shape and applied torques, is certainly something that shouldbe looked into further.76Bibliography[1] F. Talay Akyildiz and H. Bellout. Viscoelastic lubrication with phan-thein-tanner fluid. Journal of Tribology, 126:288?291, 2004.[2] M. Argentina, J. Skotheim, and L. Mahadevan. Settling and swimmingof flexible fluid-lubricated foils. Physical Review Letters, 99(22):224503.1?224503.4, 2007.[3] N. Balmforth and R. Craster J. Bush. Roll waves on flowing cornstarchsuspensions. Physics Letters A, 338:479?484, 2005.[4] G. Bayada, L. Chupin, and S. Martin. Viscoelastic fluids in thin domains.Quarterly of Applied Mathematics, 65:625?651, 2007.[5] R. B. Bird, R. C. Armstrong, and O. Hassager. Dynamics of PolymericLiquids. Fluid Mechanics. John Wiley and Sons, second edition, 1977.[6] J. R. Blake. Infinite models for ciliary propulsion. Journal of Fluid Me-chanics, 49:209?222, 1971.[7] C. J. Brokaw. Non-sinusoidal bending waves of sperm flagella. Journal ofExperimental Biology, 43:155?169, 1965.[8] B. Chan, N. J. Balmforth, and A. E. Hosoi. Building a better snail: Lu-brication and adhesive locomotion. Physics of Fluids, 17(11):113101.1?113101.10, 2005.[9] T. K. Chaudhury. On swimming in a visco-elastic liquid. Journal of FluidMechanics, 95(1):189?197, 1978.[10] R. H. Dillon, L. J. Fauci, and C. Omoto. Mathematical modeling of ax-oneme mechanics and fluid dynamics in ciliary and sperm motility. Dynam-ics of Continuous, Discrete and Impulsive Systems, 10(5):745?757, 2003.[11] E. Z. Drobnis and D. F. Katz. On Mammalian Gamete Biology, pages269?300. CRC Press, NY, 1990.77Bibliography[12] L. J. Fauci and A. McDonald. Sperm motility in the presence of boundaries.Bulletin of Mathematical Biology, 57(5):679?699, 1995.[13] H. C. Fu, T. R. Powers, and C. W. Wolgemuth. Theory of swimmingfilaments in viscoelastic media. Physical Review Letters, 99(25):258101.1?258101.4, 2007.[14] J. Gray and G. J. Hancock. The propulsion of sea-urchin spermatozoa.Proceedings of the Royal Society of London. Series A, Mathematical andPhysical Sciences, 209(1099):802?814, 1955.[15] G. J. Hancock. The self-propulsion of microscopic organisms through liq-uids. Proceedings of the Royal Society of London. Series A, Mathematicaland Physical Sciences, 217(1128):96?121, 1953.[16] J. J. L. Higdon. A hydrodynamic analysis of flagellar propulsion. Journalof Fluid Mechanics, 90(4):685?711, 1979.[17] H. C. Ho and S. S. Suarez. Hyperactivation of mammalian spermatozoa:function and regulation. Reproduction, 122:519?526, 2001.[18] A. E. Hosoi and L. Mahadevan. Peeling, healing and bursting in a lubricatedelastic sheet. Physical Review Letters, 93(13):137802, 2004.[19] R. E. Johnson and C. J. Brokaw. Flagellar hydrodynamics. a comparisonbetween resistive-force theory and slender-body theory. Biophysical Jour-nal, 23(1), 1979.[20] D . F. Katz. Movement of bull spermatozoa in cervical mucus. Biology ofReproduction, 25:931?937, 1981.[21] D . F. Katz. Characteristics of sperm motility. New York Academy ofSciences. Annals, 637:409?423, 1991.[22] D. F. Katz. On the propulsion of micro-organisms near solid boundaries.Journal of Fluid Mechanics, 64(1):33?41, 1974.[23] D. F. Katz and R. Yanagimachi. Movement characteristics of hamsterspermatozoa within the oviduct. Biology of Reproduction, 22:759?764, 1980.[24] E. Lauga. Floppy swimming: Viscous locomotion of actuated elastica.Physical review. B, Solid state, 75(4):041916, 2007.78Bibliography[25] E. Lauga. Propulsion in a viscoelastic fluid. Physics of Fluids,19(8):083104.1?083104.13, 2007.[26] E. Lauga and A. E. Hosoi. Tuning gastropod locomotion: Modelling theinfluence of mucus rheology on the cost of crawling. Physics of Fluids,18(11):113102.1?113101.9, 2006.[27] M. J. Lighthill. Flagellar hydrodynamics. SIAM Review, 18:161?230, 1976.[28] A. E. H. Love. A Treatise on the Mathematical Theory of Elasticity. GeneralPublishing Company, fourth edition, 1944.[29] K. E. Machin. Wave propagation along flagella. Journal of ExperimentalBiology, 35:796?806, 1958.[30] A. Meziane, B. Bou-Sad, and John Tichy. Modelling human hip jointlubrication subject to walking cycle. Lubrication Science, 20(3):205?222,2005.[31] J. G. Oldroyd. On the formulation of rheological equations of state. Proceed-ings of the Royal Society of London. Series A, Mathematical and PhysicalSciences, 200(1063):523?541, 1950.[32] J. G. Oldroyd. Non-newtonian effects in steady motion of some idealizedelastico-viscous liquid. Proceedings of the Royal Society of London. SeriesA, Mathematical and Physical Sciences, 245(1241):278?297, 1958.[33] A. J. Reynolds. The swimming of minute organisms. Journal of FluidMechanics, 23:241?260, 1965.[34] S. S. Suarez and X. Dai. Hyperactivation enhances mouse sperm capacityfor penetrating viscoelastic media. Biology of Reproduction, 46:686?691,1992.[35] S. S. Suarez and H. C. Ho. Hyperactivated motility in sperm. Reproductionin Domestic Animals, 38:119?124, 2003.[36] S. S. Suarez, D. F. Katz, D. H. Owen, J. B. Andrew, and R. L. Powell.Evidence for the function of hyperactivated motility in sperm. Biology ofReproduction, 44:375?381, 1991.79Bibliography[37] G. I. Taylor. Analysis of the swimming of microscopic organisms. Proceed-ings of the Royal Society of London. Series A, Mathematical and PhysicalSciences, 209(1099):447?461, 1951.[38] G. I. Taylor. The action of waving cylindrical tails in propelling micro-scopic organisms. Proceedings of the Royal Society of London. Series A,Mathematical and Physical Sciences, 211(1105):225?239, 1952.[39] J. Wilkening and A. E. Hosoi. Shape optimization of a sheet swimmingover a thin liquid layer. Journal of Fluid Mechanics, 601:25?61, 2008.[40] H. Winet. Ciliary propulsion of objects in tubes: Wall drag on swimmingtetrahymena (ciliata) in the presence of mucin and other long chain poly-mers. Journal of Experimental Biology, 64:283?302, 1976.[41] C. Wolgemuth, E. Hoiczyk, D. Kaiser, and G. Oster. How myxobacteriaglide. Current Biology, 12:369?377, 2002.[42] Y. L. Zhang, O. K. Matar, and R. V. Craster. Surfactant spreading on athin weakly viscoelastic film. Journal of Non-Newtonian Fluid Mechanics,105:53?78, 2002.80Appendix AFull Analysis of ParentProblemIn this section we discuss a fundamental problem, that of a swimmer in thepresence of a wall and a viscoelastic fluid without the lubrication approximation.Here we will assume that the amplitude of the waves is small in comparisonto the mean distance to the wall. This is a nice warm up problem as it isa combination of Katz?s biharmonic problem and Lauga?s viscoelastic Taylorproblem [22], [25].A.1 Mathematical FormulationWe assume that the sheet has prescribed motion, that of waves propagating tothe right. In the frame moving with the swimmer, these waves look likey = asin(kx -omegat). (A.1)The equations of motion for this problem aretriangleinv?u = 0, (A.2)triangleinvp = triangleinv? tau, (A.3)tau + lambda1tautriangleinv = ?parenleftBig? + lambda2 ?triangleinvparenrightBig, (A.4)which are the conservation of mass, momentum and the Oldroyd-B constitutiveequation respectively. We also have the zero force condition27angbracketlefttauxyangbracketright. = 0 (A.5)27This equation still holds even though there is no lubrication approximation here due tothe low amplitude compared to mean height approximation81Appendix A. Full Analysis of Parent ProblemThe next step is then to non-dimensionalize the equations as followst = 1omega ?, (x,y) = 1k(?, ?),u = c?, ? = omega?,(p,tau) = ?omega(?, ?),where c = omegak is the speed of the travelling wave. The conservation of momentumand conservation of mass equations do not change:triangleinvp = triangleinv? tau, (A.6)triangleinv? u = 0, (A.7)and the constitutive equation becomes:tau + De1tautriangleinv = ? + De2 ? triangleinv, (A.8)where De1 and De2 are Deborah numbers.It remains to non-dimensionalize the boundary conditions. Using the lowamplitude approximation, we putepsilon = ak lessmuch1, (A.9)and find that (dropping hats):y = epsilonsin(x -t). (A.10)From which we can easily take the time derivative to finddy/dt = -epsiloncos(x -t) (A.11)The boundary conditions become, on the wall y = hk 28u = U, (A.12)v = 0, (A.13)28note that we just write U rather than (1/c)U for simplicity82Appendix A. Full Analysis of Parent Problemand on the swimmer y = epsilonsin(x -t),u = 0, (A.14)v = -epsiloncos(x -t). (A.15)A.2 Stream Function and AsymptoticsSince the problem at hand is an incompressible, two dimensional problem, it ismost useful to consider a stream function psi such thatu = partialdiffpsipartialdiffy , (A.16)v = -partialdiffpsipartialdiffx. (A.17)Automatically we have conservation of mass holding sincetriangleinv? u = ux + vy = partialdiffpartialdiffxparenleftbiggpartialdiffpsipartialdiffyparenrightbigg+ partialdiffpartialdiffyparenleftbigg-partialdiffpsipartialdiffxparenrightbigg= 0. (A.18)It is convenient to observe that the boundary conditions can be written conciselyas:triangleinvpsivextendsinglevextendsingle(x,epsilonsin(x-t)) = epsiloncos(x -t)i (A.19)triangleinvpsivextendsinglevextendsingle(x,hk) = Uj (A.20)where i = (1,0) and j = (0,1) are the unit vectors in the x and y directionsrespectively. This form of the boundary conditions is only practical for substi-tuting in the below perturbation expansions and Taylor expanding.To solve this system, consider low amplitude solutions of the formU = epsilonU1 + epsilon2U2 + ... (A.21)psi = epsilonpsi1 + epsilon2psi2 + ... (A.22)tau = epsilontau1 + epsilon2tau2 + ... (A.23)That is, regular perturbation expansions in epsilon for the swimming speed, the streamfunction and the stress.Since it is assumed that the speed of the swimmer is constant, it must be83Appendix A. Full Analysis of Parent Problemthe case that each Uk in the expansion of U is also constant. That isUk = Const universalk = 1,2,3... (A.24)Notice that since? =parenleftBigg2ux uy + vxuy + vx 2vyparenrightBigg=parenleftBigg2psiyx psiyy -psixxpsiyy -psixx -2psixyparenrightBigg, (A.25)it is convenient to write? = epsilon ?1 + epsilon2 ?2 + ... (A.26)Similarly, we writeu = epsilonu1 + epsilon2u2 + ... (A.27)In the spirit of asymptotics, the next step is to substitute the perturbationexpansions into the equations of motion and the boundary conditions and sortout the first and second order equations. Before this is done, it is most usefulto observe thattriangleinv?triangleinvp = 0, (A.28)and so pressure can be removed from the conservation of momentum equation.Summarizing the equations of motion then givestriangleinv?triangleinv?tau = 0 (A.29)tau + De1tautriangleinv = ? + De2 ? triangleinv. (A.30)Substituting in the perturbation expansions yields at O(epsilon):triangleinv?triangleinv?tau1 = 0, (A.31)tau1 + De1 partialdiffpartialdiffttau1 = ?1 + De2 partialdiffpartialdifft ?1, (A.32)and at O(epsilon2):triangleinv?triangleinv?tau2 = 0, (A.33)tau2 + De1 partialdiffpartialdiffttau2 + De1 parenleftbigu1 ? triangleinvtau1 -tau1 ? triangleinvu1 -(triangleinvu1)T ? tau1parenrightbig =?2 + De2 partialdiffpartialdifft ?2 + De2 parenleftbigu1 ? triangleinv ?1 - ?1 ? triangleinvu1 -(triangleinvu1)T ? ?1parenrightbig . (A.34)Now it is time to look at the boundary conditions. Substituting in the regular84Appendix A. Full Analysis of Parent Problemperturbation expansion for psi and expanding about y = 0 yields at O(epsilon):triangleinvpsi1vextendsinglevextendsingle(x,0) = cos(x -t)i, (A.35)triangleinvpsi1vextendsinglevextendsingle(x,hk) = U1j, (A.36)and at O(epsilon2):triangleinvpsi2vextendsinglevextendsingle(x,0) + triangleinvparenleftbiggpartialdiffpsi1partialdiffyparenrightbiggvextendsinglevextendsinglevextendsingle(x,0)sin(x -t) = 0, (A.37)triangleinvpsi2vextendsinglevextendsingle(x,hk) = U2j. (A.38)Substituting the expansions into the constraint angbracketlefttauxyangbracketright = 0 yields thatangbracketlefttaukxyangbracketright = 0 k = 1,2,3... (A.39)A.3 Leading Order EquationsIf one stares at the leading order equations for a moment,triangleinv?triangleinv?tau1 = 0, (A.40)tau1 + De1 partialdiffpartialdiffttau1 = ? 1 + De2 partialdiffpartialdifft ? 1, (A.41)it becomes apparent that the easiest thing to do is to take the divergence andthen the curl of the second equation thus eliminating leading order stress. Bydoing this we obtainparenleftbigg1 + De2 partialdiffpartialdifftparenrightbiggtriangleinv?triangleinv? ? 1 = 0. (A.42)When substituting in? 1 = ? 1(psi1xx,psi1xy,psi1yy),we actually get parenleftbigg1 + De2 partialdiffpartialdifftparenrightbiggtriangleinv4 psi1 = 0 (A.43)wheretriangleinv4psi = psixxxx + 2psixxyy + psiyyyy.85Appendix A. Full Analysis of Parent ProblemWriting down the boundary conditions in a more practical fashion givespartialdiffpsi1partialdiffxvextendsinglevextendsingle(x,0) = 0,partialdiffpsi1partialdiffyvextendsinglevextendsingle(x,0) = 0,partialdiffpsi1partialdiffxvextendsinglevextendsingle(x,hk) = cos(x -t),partialdiffpsi1partialdiffyvextendsinglevextendsingle(x,hk) = U1.(A.44)One way to solve this is to simply use separation of variables, cleverly guess-ing an initial form of the solution based on the boundary conditionspsi1 = f(y)sin(x -t), (A.45)where f(y) is to unknown. Substituting this into the leading order equation forpsi yieldsd4fdy4 -2d2fdy2 + f = 0. (A.46)It is easily found thatf(y) = (a + by)e-y + (c + dy)ey, (A.47)where constants a, b, c, and d are found by applying the boundary conditionsand will depend on h. When applying the boundary condition psi1yvextendsinglevextendsingle(x,hk) = U1 itbecomes apparent that in order for the assumption that U = constant to hold,it must be the case thatU1 = 0 (A.48)or else U1 would depend on x and t, which is unacceptable.This is a result that we could have easily predicted based on the work donepreviously by others, who all found the leading order swimming speed to bezero [37], [22], [25].Alternately, one can solve this using hyperbolic sines and cosines.The procedure is the same as above, only this time a different form for f(y)will be used. Putpsi1(y) = [(A + BH)cosh(H) + (C + DH)sinh(H)]sin(x -t).This givespsi1x = [(A + BH)cosh(H)++(C + D(H))sinh(H)]cos(x -t),psi1y = -[(A + BH + D)sinh(H)++(B + C + DH)cosh(H)]sin(x -t),86Appendix A. Full Analysis of Parent Problemwhere H = (hk -y). Upon applying the boundary conditions, for y = hk:Acos(x -t) = 0,-(B + C)sin(x -t) = U1.Clearly, in order for the assumption, U1 = Const, to hold, it must be the casethat U1 = 0, thus implying that B = -C. Also, it is clear from the firstcondition that A = 0. Thus so far,psi1 = [B[sinh(hk-y)-(hk-y)cosh(hk-y)]+D(hk-y)sinh(hk-y)]sin(x-t).Or, to simplify the x-derivatives, one could consider the complex analogy of theproblem (?complexification?, see below). That is putpsi1 = 12parenleftBig ?psi1ei(x-t) + ?1asteriskmathe-i(x-t)parenrightBig, (A.49)and consider only?psi1ei(x-t).One then gets:?psi1 = [ ?B[sinh(hk -y) -(hk -y)cosh(hk -y)] + ?D(hk -y)sinh(hk -y)].DefinePsi1(y) = (hk -y)sinh(hk -y), (A.50)Psi2(y) = sinh(hk -y) -(hk -y)cosh(hk -y), (A.51)thus giving?psi1 = [ ?DPsi1 + ?BPsi2].Now, the unknown coefficients ? and ? must be determined using the otherboundary conditions. On y = 0:parenleftBig ?psi1ei(x-t)parenrightBigxvextendsinglevextendsinglevextendsingley=0= iei(x-t)[ ?Psi1(0) + ?Psi2(0)] = ei(x-t),parenleftBig?psi1ei(x-t)parenrightBigyvextendsinglevextendsinglevextendsingley=0= ei(x-t)[ ?Psiprime1(0) + ?Psiprime2(0)] = 0.Solving these gives?D = iPsiprime (0)W(0) , ?B = -iPsiprime (0)W(0) ,87Appendix A. Full Analysis of Parent ProblemwhereW(y) = Psiprime1(0)Psi2(y) -Psiprime2(0)Psi1(y). (A.52)Putting this altogether one has?psi1 = -iW(y)W(0) (A.53)or writing down the real component,psi1 = W(y)W(0) sin(x -t) (A.54)Although this second method of solving for the leading order stream functionand swimming speed is longer, it provides a more concise final form. This isuseful since the leading order swimming speed is zero, so the calculation mustbe done to second order. This aesthetically pleasing form of the leading ordersolution will make the second order problem much easier and less algebraicallychallenging to solve.To solve for the leading order stress, we employ the method of ?complexifi-cation? (referred to as Fourier notations in [25]). That is, consider the leadingorder equation?1 + De1partialdifft ?1 = ?1 + De2partialdifft ?1,and think of the stress and rate of strain as being the real portion of a larger,complex problem (note the periodic dependence explicitly given due to the as-sumed periodicity of the problem)tau1 = Rfractur[ ?1ei(x-t)], ?1 = Rfractur[ ?1ei(x-t)]. (A.55)Alternately, one could writetau1 = 12 parenleftbig ?1ei(x-t) + ?1asteriskmathe-i(x-t)parenrightbig , ?1 = 12 parenleftbig ?1ei(x-t) + ?1asteriskmathe-i(x-t)parenrightbig ,(A.56)observing that ?1 and ?1 are functions of y only, and are independent of x andt (note that this is also true for ?1 defined above).This ?complexification? of the problem implies thatpartialdifft = -partialdiffx = -i.88Appendix A. Full Analysis of Parent ProblemTherefore, by exploiting this simplification of time derivatives, we can easilyobtain?1 = 1 -iDe21 -iDe1??gamma1. (A.57)Upon substituting in ?1 = ?1( ?1xx, ?1xy, ?1yy) as defined above, we get?1 = (1 -iDe2)ei(x-t)(1 -iDe1)W(0)parenleftBigg2Wprime -i(Wprimeprime + W)-i(Wprimeprime + W) -2WprimeparenrightBigg. (A.58)Notice that the condition,angbracketlefttau1xyangbracketright = 0,automatically holds due to the fact thattau1xy proportional ?1xy proportional sin(x -t).A.4 Second Order EquationsThe leading order problem has been solved and it has been shown that in order tomaintain the assumption that the speed of the swimmer is constant, at leadingorder it must be equal to zero as predicted. Thus, in order to calculate theswimming speed of the organism, one must move on to the second order problem.This problem is not quite as straightforward as the leading order problem dueto the nonlinear terms of the upper convected derivative that crop up. Let usrestate the second order problem:triangleinv?triangleinv? tau2 = 0,parenleftbigg1 + De1 partialdiffpartialdifftparenrightbiggtau2 -parenleftbigg1 + De2 partialdiffpartialdifftparenrightbigg?2 = De2 parenleftbigu1 ? triangleinv ?1 - ?1 ? triangleinvu1 -(triangleinvu1)T ? ?1parenrightbig --De1 parenleftbigu1 ? triangleinvtau1 -tau1 ? triangleinvu1 -(triangleinvu1)T ? tau1parenrightbig ,with boundary conditions:triangleinvpsi2vextendsinglevextendsingle(x,0) + triangleinvparenleftBigpartialdiffpsi1partialdiffyparenrightBigvextendsinglevextendsinglevextendsingle(x,0)sin(x -t) = 0,triangleinvpsi2vextendsinglevextendsingle(x,hk) = U2j.89Appendix A. Full Analysis of Parent ProblemThe zero force condition states thatangbracketlefttau2xyangbracketright = 0. (A.59)To start things off slowly, consider first the boundary conditions. Writingthings in a more user friendly fashion gives on y = 0:psi2xvextendsinglevextendsingley=0 = -sin(x -t)psi1yx, psi2yvextendsinglevextendsingley=0 = -sin(x -t)psi1yy.But we know psi1 from solving the leading order problem, which gives psi1yx =Wprime(y)W(0) cos(x - t) recalling that W(y) = Psiprime1(0)Psi2(y) - Psiprime2(0)Psi1(y) and Psi1(y) =(hk - y)sinh(hk - y),Psi2(y) = sinh(hk - y) - (hk - y)cosh(hk - y). Thus, weobtainpsi2xvextendsinglevextendsingley=0 = 0, (A.60)psi2yvextendsinglevextendsingley=0 = -sin2(x -t)Wprimeprime(0)W(0) . (A.61)On y = hk:psi2xvextendsinglevextendsingley=hk = 0, (A.62)psi2yvextendsinglevextendsingley=hk = U2. (A.63)Inspecting the second order equations of motion as written above revealsthat the right-hand side of the momentum equation is composed of leading orderterms exclusively. Since it is the case that the leading order stress can be simplyrepresented in terms of the leading order rate of strain when complexified, let usput the right-hand side of the momentum equations in terms of the complexifiedvariablesu1 = 12(?1ei(x-t) + ?asteriskmath1e-i(x-t)), ?1 = 12( ?1ei(x-t) + ?1asteriskmathe-i(x-t)),and of course?1 = 1 -iDe21 -iDe1??gamma1.90Appendix A. Full Analysis of Parent ProblemWe obtainRHS = (De2 -De1)4(1 -iDe1)bracketleftbig?uasteriskmath1 ? triangleinv??gamma1 - ??gamma1 ? triangleinv?uasteriskmath1 -(triangleinv?uasteriskmath1)T ? ??gamma1bracketrightbig + C.C. ++ (De2 -De1)4(1 -iDe1)e2i(x-t) bracketleftbig?1 ? triangleinv?1 - ?1 ? triangleinv?1 -(triangleinv?1)T ? ?1bracketrightbig + C.C.,(A.64)where C.C. refers to the complex conjugate of the term directly preceding it.29To solve this rather long, complicated looking equation, make use of theperiodicity of the problem and look at the x-average of the equation. Notice atthis point one could take the curl and divergence of the problem to remove the tau2terms, but that is rather cumbersome. Rather, only consider the xy-componentand utilize the zero force condition to simplify the problem.Let us first look at the boundary conditions. Note that angbracketleftpsi2xangbracketright = 12pipsi2vextendsinglevextendsingle2pi0 = 0always, so the boundary conditions for psi2x yield no information. However, theother two boundary conditions have non-trivial averages:angbracketleftpsi2yangbracketrightvextendsinglevextendsingley=0 = -Wprimeprime(0)2W(0), (A.66)angbracketleftpsi2yangbracketright > vextendsinglevextendsingley=hk = U2. (A.67)Now let us average the second order constitutive equation and look at the xy-component (recall that partialdifft = -partialdiffx and that angbracketlefte?2i(x-t)angbracketright = 0), subscripts herereferring to components:angbracketlefttau2xyangbracketright -angbracketleft ?2xyangbracketright = Rfracturbraceleftbigg(De2 -De1)2(1 -iDe1)bracketleftbig?uasteriskmath1 ? triangleinv??gamma1 - ??gamma1 ? triangleinv?uasteriskmath1 -(triangleinv?uasteriskmath1)T ? ??gamma1bracketrightbigxybracerightbigg.(A.68)Applying the zero force condition and writing the rate of strain component in29Alternately, one could writeRHS = Rfracturbraceleftbigg(De2 -De1)2(1 -iDe1)bracketleftBig?asteriskmath1 ? triangleinv?1 - ?1 ? triangleinv?asteriskmath1 -(triangleinv?asteriskmath1)T ? ?1bracketrightBigbracerightbigg++ Rfracturbraceleftbigg(De2 -De1)2(1 -iDe1) e2i(x-t)bracketleftBig?1 ? triangleinv?1 - ?1 ? triangleinv?1 -(triangleinv?1)T ? ?1bracketrightBigbracerightbigg(A.65)91Appendix A. Full Analysis of Parent Problemterms of the stream function gives:angbracketleftpsi2yy -psi2xxangbracketright = Rfracturbraceleftbigg(De1 -De2)2(1 -iDe1)bracketleftbig?uasteriskmath1 ? triangleinv??gamma1 - ??gamma1 ? triangleinv?uasteriskmath1 -(triangleinv?uasteriskmath1)T ? ??gamma1bracketrightbigxybracerightbigg.(A.69)But, angbracketleftpsi2xxangbracketright >= 12pipsi2xvextendsinglevextendsingle2pi0 = - 12pivvextendsinglevextendsingle2pi0 = 0 due to the assumed periodicity of thesolution. Since the bounds of the average are not dependent on y, we can putangbracketleftpsi2angbracketrightyy = Rfracturbraceleftbigg(De1 -De2)2(1 -iDe1)bracketleftbig?uasteriskmath1 ? triangleinv??gamma1 - ??gamma1 ? triangleinv?uasteriskmath1 -(triangleinv?uasteriskmath1)T ? ??gamma1bracketrightbigxybracerightbigg. (A.70)Thus, to solve the problem and calculate the swimming speed we must simplyintegrate once to obtain an expression for angbracketleftpsi2angbracketrighty and apply the boundary condi-tions. To to this, the xy-component of the right-hand side must be calculated.By using the fact that both ?1 and ?1 are functions of ?1 = -iW(y)W(0) = -i WW0 , wecan write down the terms on the right-hand side in a relatively concise manner.Namely, we obtainbracketleftbig ?u1asteriskmath ? triangleinv??gamma1bracketrightbigxy =iW20 (WprimeWprimeprime + 2WWprime + WWprimeprimeprime), (A.71)bracketleftbig ??gamma1 ? triangleinv ?u1asteriskmathbracketrightbigxy =iW20 (WprimeWprimeprime + 3WWprime), (A.72)bracketleftbigtriangleinv ?u1asteriskmathT ? ??gamma1bracketrightbigxy =-iW20 (3WprimeWprimeprime + WWprime). (A.73)Substituting these into the governing equation givesangbracketleftpsi2angbracketrightyy = Rfracturbraceleftbigg(De1 -De2)2(1 -iDe1)iW20 [WWprimeprimeprime + 3WprimeWprimeprime]bracerightbigg. (A.74)But, notice that [WWprimeprimeprime + 3WprimeWprimeprime] can be written as [WWprimeprime + Wprime2]prime and so wefind thatangbracketleftpsi2angbracketrightyy = De1(De2 -De1)2(1 + De21)W2bracketleftbigWWprimeprime + Wprime2bracketrightbigprime . (A.75)Integrating and applying the first boundary condition gives- Wprimeprime02W0 =De1(De2 -De1)2(1 + De21)W20bracketleftbigW0Wprimeprime + Wprime2bracketrightbig + Const.Using the fact that Wprime0 = 0 the constant of integration can be easily solved for,Const = - Wprimeprime02W0(1 + De1De2)(1 + De21) .92Appendix A. Full Analysis of Parent ProblemThus,angbracketleftpsi2angbracketrighty = De1(De2 -De1)2(1 + De21)W2bracketleftbigWWprimeprime + Wprime2bracketrightbig - Wprimeprime02W0(1 + De1De2)(1 + De21) , (A.76)and so finally one has the swimming speedU2 = angbracketleftpsi2angbracketrightyvextendsinglevextendsingley=hk = - Wprimeprime02W0(1 + De1De2)(1 + De21) , (A.77)since W(hk) = 0 and Wprime(hk) = 0.A.5 WorkThe leading order rate of work (per unit wavelength) is given by the volumeintegral ofangbracketlefttau1 : ?1angbracketright, (A.78)and is equal to the rate of viscous dissipation. To simplify this calculation, welook again to the ?complexified? variables. That is, puttau1 : ?1 = 14 ( ?1 + ?1asteriskmath) : parenleftbig ?1 + ?1asteriskmathparenrightbig .Let us calculate the integral in x first. It is the case thatangbracketleft ?1 : ?1angbracketright proportional angbracketlefte2i(x-t)angbracketright = 0,angbracketleft ?1asteriskmath : ?1asteriskmathangbracketright proportional angbracketlefte-2i(x-t)angbracketright = 0.Since both tau1 and ?1 are symmetric, we can putangbracketlefttau1 : ?1angbracketright = 12Rfracturbraceleftbigangbracketleft ?1 : ?1asteriskmathangbracketrightbracerightbigKnowing the form of ?1 in terms of ?1 givesangbracketlefttau1 : ?1angbracketright = 12Rfracturbraceleftbiggparenleftbigg1 -iDe21 -iDe1parenrightbiggangbracketleft ?1 : ?1asteriskmathangbracketrightbracerightbiggUsing the dependence of ?1 on psi1, and the fact that ?1 = -i WW0 , gives usangbracketlefttau1 : ?1angbracketright = 1 + De1De21 + De21parenleftbigg4Wprime2 + (W + Wprimeprime)2W20parenrightbigg. (A.79)93Appendix A. Full Analysis of Parent ProblemThus, leading order rate of work per unit wavelength, is given by(1 + De1De2)(1 + De21)W20integraldisplay hk0bracketleftbig4Wprime2 + (W + Wprimeprime)2bracketrightbig dy. (A.80)A.6 Limiting CasesA.6.1 Distance Between Wall and Swimmer ApproachesInfinityTo make sure that this problem agrees with the infinite case considered by Laugain [25], we must take the limit as the distance between the wall and the swimmerapproaches infinity. It is known thatpsi1 = W(y)W(0) sin(x -t),recalling thatW(y) = Psiprime1(0)Psi2(y) -Psiprime2(0)Psi1(y),andPsi1(y) = (hk -y)sinh(hk -y),Psi2(y) = sinh(hk -y) -(hk -y)cosh(hk -y).The limit as h -arrowrightinfinity is easily calculated. It is found thatlimh-arrowrightinfinityW(y)W(0) = (1 + y)e-yThuslimh-arrowrightinfinitypsi1 = (1 + y)e-y sin(x -t)which is exactly the infinite solution as stated in [25].It is obvious that if the limit of the stream function agrees with the infinitecase that the limit of the swimming speed will as well. But to illustrate thisclearly, recall thatU2 = -Wprimeprime(0)2W(0)(1 + De1De2)(1 + De21) .So, the limitlimh-arrowrightinfinityU2 = 12 (1 + De1De2)(1 + De21),which again agrees with the infinite case.94Appendix A. Full Analysis of Parent ProblemSimilarly, we can calculate the limit for the averaged rate of work done bythe swimmer. It is found thatlimh-arrowrightinfinitybraceleftBigg(1 + De1De2)(1 + De21)W20integraldisplay hk0bracketleftbig4Wprime2 + (W + Wprimeprime)2bracketrightbig dybracerightBigg = 2(1 + De1De2)(1 + De21) .We cannot directly compare this limit with the form of the solution found in[25] as it is stated asWorkinfinityWN =(1 + De1De2)(1 + De21) ,and WN is not stated. However, we can easily see that the calculation of workin the infinite case is exactly the same with the exception of the W function. Inthe infinite case one hasWinfinity = (1 + y)e-y,and so it is found that(1 + De1De2)(1 + De21)Winfinity(0)2integraldisplay infinity0bracketleftbig4Wprime2infinity + (Winfinity + Wprimeprimeinfinity)2bracketrightbig dy = 2(1 + De1De2)(1 + De21) ,exactly matching the limit calculated above.So, it is clear that the above calculations for the problem of an organismswimming in an Oldroyd-B viscoelastic fluid near a wall are consistent withthose of the problem of an organism swimming in an infinite fluid of the sametype.A.6.2 Swimmer Approaches the WallNext, we would like to consider the limit as the swimmer and the wall approacheach other. There are two ways of thinking of this limit. The first is keeping thewall fixed at y = h and increasing the amplitude of the swimmer until it reachesthe wall. This line of thinking contradicts the assumption that ak << 1, as itimplies that O(a) arrowright O(h) = O(1/k). This would make our present solutioninvalid. The second way of thinking about this limit is to think of the wallapproaching the swimmer, whose amplitude and wavelength remains the fixed.This implies that O(h) arrowrightO(a), still maintaining the assumption that ak << 1.In fact, in this manner of thinking we get that hk << 1 which is exactly the95Appendix A. Full Analysis of Parent Problemlubrication approximation. We have that, to leading orderU = -epsilon2 (1 + De1De2)2(1 + De21)F(hk), (A.81)whereF(hk) = Wprimeprime0W0 =(hk)2 + sinh2(hk)(hk)2 -sinh2(hk). (A.82)In our limit, hk arrowrightak = epsilon and so we have thatF(epsilon) = epsilon2 + sinh2 epsilonepsilon2 -sinh2 epsilon. (A.83)Expanding sinhepsilon = epsilon + 16epsilon3 + ... and simplifying gives usF(epsilon) = 2epsilon2 + +13epsilon4 + ...-13epsilon4 -... . (A.84)Thus, in the limit when the wall approaches the swimmer, we find that theswimming speed isU = 3(1 + De1De2)(1 + De21). (A.85)96Appendix BDetails of Bingham AnalysisB.1 Case IBeginning with the first case, Y- < Y+ < Y < 1, let us look at equation(2.108). Considering the region where y > Y+, we can simplify this equationusing equation (2.114). Integrating,integraldisplay yYuydy =integraldisplay yY(y -Y+)pxdy, (B.1)givesu = px2 parenleftbig(y -Y+)2 -(Y -Y+)2parenrightbig . (B.2)In this region, y > Y+, so we can see that this flow profile is given by a concave-up parabola in y. Applying the boundary condition at y = 1 and rearrangingto solve for Y+ gives,Y+ = 12(1 + Y ) - Upx(1 -Y ). (B.3)Technically, the plug region is outside the boundaries we are considering, but, forthe purpose of solving this set of equations, let us imagine that the flow profilecontinues through y = Y . Since the flow profile is a concave-up parabola, wecan see that up < 0. As well, applying the condition on y = Y+ gives us thatup = -px2 (Y -Y+)2. (B.4)In turn, this tells ussgn(up) = -sgn(px). (B.5)Thus, we find thatsgn(px) > 0, (B.6)97Appendix B. Details of Bingham Analysisand so sigma = +1. We can now look at equation (2.119). Substituting in u andsimplifying, givesq = Y + U2 (1 -Y ) - px12(1 -Y )3. (B.7)B.2 Case IIFor case II, we have Y- < Y < Y+ < 1. Let us begin with the region whereY- < y < Y+. Since we have that, u|(y=Y ) = 0, and Y- < Y < Y+, we knowautomatically thatup = 0. (B.8)Now, consider the region where Y+ < y < 1. Taking equation (2.108), simplify-ing using equation (2.114), and integrating,integraldisplay yY +uydy =integraldisplay yY +(y -Y+)pxdy, (B.9)gives usu = px2 (y -Y+)2. (B.10)Applying the boundary condition at y = 1 and solving for Y+ gives,Y+ = 1 -radicalBigg2Upx . (B.11)Since we are assuming that U > 0, the above equation impliespx > 0. (B.12)This implies that sigma = +1. Moving on to equation (2.119),q = Y +integraldisplay 1Y +udy. (B.13)Substituting in our expression for u and simplifying, givesq = Y + px6 (1 -Y+)3. (B.14)98Appendix B. Details of Bingham AnalysisB.3 Case IIICase III is the case when the pseudo-plug is contained within the boundariesof the swimmer and wall, Y < Y- < Y+ < 1. In the pseudo-plug region,Y- < y < Y+, we have that u = up, where up is currently unknown. In theupper region, Y+ < y < 1, we can again simplify equation (2.108) using equation(2.114) and integrate,integraldisplay yY +uydy =integraldisplay yY +(y -Y+)pxdy. (B.15)This gives us thatu = up + px2 (y -Y+)2. (B.16)Applying the boundary condition on y = 1 and solving for up yields,up = U - px2 (1 -Y+)2. (B.17)Now, consider the lower yielded region, Y < y < Y-. In this region, we simplifyequation (2.108) using equation (2.115). Integrating,integraldisplay Y -yuydy =integraldisplay Y -y(y -Y-)pxdy, (B.18)givesu = up + px2 (y -Y-)2. (B.19)Applying the boundary condition at y = Y and solving for up gives,up = -px2 (Y -Y-)2. (B.20)Equating these two equations for up and solving for Y+ gives,Y+ =parenleftBig1 -Y 2 - 2Upx -4 B|px|parenleftBigY + B|px|parenrightBigparenrightBig2(1 -Y -2 B|px| ) (B.21)recalling that, Y- = Y+ - 2 B|px| , from equation (2.117). We can now look atequation (2.119),q = Y +integraldisplay 1Y +udy + up(Y+ -Y-) +integraldisplay Y -Yudy. (B.22)99Appendix B. Details of Bingham AnalysisSubstituting in u and simplifying, givesq = Y + px6 parenleftbig(1 -Y+)3 + (Y- -Y )3 -3(1 -Y )(Y- -Y )2parenrightbig . (B.23)B.4 Case IVFor case IV, we have Y < Y- < 1 < Y+. Let us begin with the region whereY- < y < Y+. Since we have that, u|(y=1) = U, and Y- < 1 < Y+, we knowautomatically thatup = U. (B.24)Now, consider the region where Y < y < Y-. Taking equation (2.108), simpli-fying using equation (2.115), and integrating,integraldisplay Y -yuydy =integraldisplay Y -y(y -Y-)pxdy, (B.25)gives usu = U + px2 (y -Y-)2. (B.26)Applying the boundary condition at y = Y and solving for Y- gives,Y- = Y +radicalBigg-2Upx . (B.27)Since we are assuming that U > 0, we know that, in this case,px < 0. (B.28)This implies that sigma = -1. Looking at equation (2.119),q = Y + U(1 -Y-) +integraldisplay Y -Yudy, (B.29)substituting in u and simplifying, givesq = Y + px6 (Y- -Y )3 + U(1 -Y ). (B.30)100Appendix B. Details of Bingham AnalysisB.5 Case VFinally, the last case, Y < 1 < Y- < Y+. Look again at equation (2.108)and consider the region where y < Y-. Simplifying using equation (2.115) andintegrating, integraldisplay yYuydy =integraldisplay yY(y -Y-)pxdy, (B.31)givesu = px2 parenleftbig(y -Y-)2 -(Y -Y-)2parenrightbig . (B.32)In this region, y < Y-, so we can see that this flow profile is given by a concave-down parabola in y. Applying the boundary condition at y = 1 and rearrangingto solve for Y- gives,Y- = 12(1 + Y ) - Upx(1 -Y ). (B.33)As in case I, the plug region is technically outside the boundaries we are consid-ering. But, for the purpose of solving this set of equations, let us again imaginethat the flow profile continues through y = 1. Since the flow profile is a concave-down parabola, we can see that up > 0. As well, may applying the conditionon y = Y- gives us thatup = -px2 (Y -Y-)2. (B.34)Again, this tells us,sgn(up) = -sgn(px), (B.35)and so we find thatsgn(px) < 0, (B.36)which gives us that sigma = -1. We can now look at equation (2.119). Substitutingin u and simplifying, givesq = Y + U2 (1 -Y ) - px12(1 -Y )3. (B.37)Note that for all the above cases, equation (2.117) allows us to put Y- interms Y+, or vice versa,Y- = Y+ - 2B|px|. (B.38)101Appendix CCode For Chapter TwoBelow is the MATLAB code used to obtain the results in Chapter Two.C.1 Viscoelastic Problem Code% NEWTON?S METHODclearN = 50; h = 2*pi/N;xi = 0:h:2*pi; % uniform meshDimWork=1; % Fixed Dimensional Workmu=1; % viscosityhh=0.1; % mean depthfor k=1:199ah = k*0.005; % a/h < 1;disp([?a/h = ?,num2str(ah)])Y = ah*sin(xi); % REGULAR%Y = ah*tanh(12*sin(xi)); % SQUARE WAVEa0=0; b0=0; % Newtonian - could have a=b=anything....a1=1; b1=2; % Thinninga2=2; b2=1; % Thickening[U1, q1, p_x1, g11, gY1, tau_01, G]=SolveSystem(N,xi,ah,Y,a1,b1);[U2, q2, p_x2, g12, gY2, tau_02, G]=SolveSystem(N,xi,ah,Y,a2,b2);U0=G(1);102Appendix C. Code For Chapter Twoq0=G(2);g10=G(N+4:2*N+4);gY0=G(2*N+5:3*N+5);u0(k)=U0; u1(k)=U1; u2(k)=U2;MF0(k)=q0; MF1(k)=q1; MF2(k)=q2;AH(k)=ah;% must calculate dimensionless work....% wi = <<2ui_z\taui_{xz}>>% Newtonian ; T=g ;p_x0 = (g10-gY0)./(1-Y);W0=@(g) 2*g.^2;% Tprime0 = 1;I0 = (2/3)*(g10.^3-gY0.^3)./(p_x0);II0 = I0(1:N);w0(k)=sum(II0.*diff(xi));;% ThinningW1Tprime1=@(g) 2.*g.^2.*(1+a1*g.^2)./(1+b1*g.^2).*(1-b1*g.^2+3*a1*g.^2+a1*b1*g.^4)./(1+b1*g.^2).^2;for j=1:N+1Int1(j)=quadl(W1Tprime1,gY1(j),g11(j));endI1=Int1./(p_x1);II1=I1(1:N);w1(k)=sum(II1.*diff(xi));% ThickeningW2Tprime2=@(g) 2.*g.^2.*(1+a2*g.^2)./(1+b2*g.^2).*(1-b2*g.^2+3*a2*g.^2+a2*b2*g.^4)./(1+b2*g.^2).^2;for j=1:N+1Int2(j)=quadl(W2Tprime2,gY2(j),g12(j));endI2=Int2./(p_x2);II2=I2(1:N);w2(k)=sum(II2.*diff(xi));103Appendix C. Code For Chapter Two% must calculate wave speed ...c0(k)=DimWork*hh^2/(2*mu*w0(k)); % Newtonianc1(k)=DimWork*hh^2/(2*mu*w1(k)); % Thinningc2(k)=DimWork*hh^2/(2*mu*w2(k)); % Thickeninguu0(k)=c0(k)*U0;uu1(k)=c1(k)*U1;uu2(k)=c2(k)*U2;endfigure plot(AH,uu0,?r?,AH,uu1,?b?,AH,uu2,?g?) title([?DimensionalSwimming Speed?],?FontSize?,12) xlabel(?a/h?,?FontSize?,12)ylabel(?U*c?,?FontSize?,12)104Appendix C. Code For Chapter Twofunction [U, q, p_x, gam1, gamY, tau_0, G] =SolveSystem(N,xi,ah,Y,a,b)E = 10^(-4); maxit = 60;% Initial guess ...G = InitialApproxNewt(N,xi,ah,Y);T = @(x) x.*(1+a*x.^2)./(1+b*x.^2); Tprime = @(x)(1-b*x.^2+3*a*x.^2+a*b*x.^4)./(1+b*x.^2).^2;I1prime = @(g)g.^2.*(1+a*g.^2).*(1-b*g.^2+3*a*g.^2+a*b*g.^4)./(1+b*g.^2).^3;% Observing that -dI1/Dg1 = dI1/dgYI0prime = @(g) g.*(1-b*g.^2+3*a*g.^2+a*b*g.^4)./(1+b*g.^2).^2;% Observing that -dI0/Dg1 = dI0/dgYJ = zeros(2*N+4); % most entries are zeroU=G(1); q=G(2); gam1=G(N+4:2*N+4); gamY=G(2*N+5:3*N+5);Guess = [U q gam1 gamY];for i=1:maxit% Function evaluated at current guessPP = (feval(T,gamY)-((feval(T,gam1)-feval(T,gamY))./(1-Y)).*Y);P = PP(1:N);F(1) = sum(P.*diff(xi));QQ = ((feval(T,gam1)-feval(T,gamY))./(1-Y));Q = QQ(1:N);F(2) = sum(Q.*diff(xi));for j=1:N+1I1 = quadl(I1prime,gamY(j),gam1(j));I0 = quadl(I0prime,gamY(j),gam1(j));% F3?sF(j+2) = U*(feval(T,gam1(j))105Appendix C. Code For Chapter Two- feval(T,gamY(j)))./(1-Y(j))*feval(T,gamY(j))- ((feval(T,gam1(j))-feval(T,gamY(j)))./(1-Y(j)))^2*(q-U-Y(j)+U*Y(j))- I1;% F4?sF(j+N+3) =I0 - U*(feval(T,gam1(j))-feval(T,gamY(j)))./(1-Y(j));% JACOBIAN made gam1 = x1, gamY = x2, q = x3, U = x4% dF3(1..N+1)/dxiJ(j+2,1) = (feval(T,gam1(j))- feval(T,gamY(j)))./(1-Y(j))*feval(T,gamY(j))- ((feval(T,gam1(j))-feval(T,gamY(j)))./(1-Y(j)))^2*(Y(j)-1);J(j+2,2) = -((feval(T,gam1(j))-feval(T,gamY(j)))./(1-Y(j)))^2;J(j+2,j+2) = U*feval(T,gamY(j))*feval(Tprime,gam1(j))/(1-Y(j))- feval(I1prime,gam1(j)) - (q-U-Y(j)+U*Y(j))*2*(feval(T,gam1(j))- feval(T,gamY(j)))*feval(Tprime,gam1(j))/(1-Y(j))^2;J(j+2,j+N+3) = U*(feval(T,gam1(j))*feval(Tprime,gamY(j))- 2*feval(T,gamY(j))*feval(Tprime,gamY(j)))/(1-Y(j))+ feval(I1prime,gamY(j)) + (q-U-Y(j)+U*Y(j))*2*(feval(T,gam1(j))-feval(T,gamY(j)))*feval(Tprime,gamY(j))/(1-Y(j))^2;% dF4(1..N+1)/dxiJ(j+N+3,1) = -(feval(T,gam1(j))-feval(T,gamY(j)))./(1-Y(j));J(j+N+3,j+2) = feval(I0prime,gam1(j))-U*feval(Tprime,gam1(j))/(1-Y(j));J(j+N+3,j+N+3) = -feval(I0prime,gamY(j))+U*feval(Tprime,gamY(j))/(1-Y(j));endfor j=1:N% Constraints% dF1/DxiJ(1,j+2) = -Y(j)/(1-Y(j))*feval(Tprime,gam1(j))*(xi(j+1)-xi(j));J(1,j+N+3) = feval(Tprime,gamY(j))+Y(j)/(1-Y(j))*feval(Tprime,gamY(j))*(xi(j+1)-xi(j));% dF2/DxiJ(2,j+2) = 1/(1-Y(j))*feval(Tprime,gam1(j))*(xi(j+1)-xi(j));J(2,j+N+3) = -1/(1-Y(j))*feval(Tprime,gamY(j))*(xi(j+1)-xi(j));endy = -J\F?;106Appendix C. Code For Chapter TwoGuess = Guess + y?;U=Guess(1);q=Guess(2);gam1=Guess(3:N+3);gamY=Guess(N+4:2*N+4);if norm(y,inf) < Edisp([?Method converging in ?,num2str(i),? iterations?])breakendend if i==maxitdisp([?NEWTONS METHOD NOT CONVERGING IN: ?,num2str(maxit),? ITERATIONS - ERROR?])end tau_0 = feval(T,gamY); p_x =(feval(T,gam1)-feval(T,gamY))./(1-Y);107Appendix C. Code For Chapter TwoC.2 Bingham Problem Code%ssdriva = 0.5; B = 0.0;N=400; x=([0:N-1]+1/2)/N*pi-pi/2; px=a*sin(x); X=zeros(N+2,1);X(1:N)=px;for niter=1:30 ssrhs dX=-EJ\E; conno=max(abs(dX)) X=X+dX; ifconno<1e-12 break end endsubplot(211) plot(x,px)UU=X(N+1); qq=X(N+2); [UU qq]a0=a; aa=a0;%for na = 1:1%a = a0+na/20;B=1;for niter=1:30 ssrhs dX=-EJ\E; conno=max(abs(dX)) if niter==1,dX=dX/2;end X=X+dX; if conno<1e-8 break end end%UU=[UU X(2*N+1)];%qq=[qq X(2*N+2)];%aa=[aa a];hold on, plot(x,px,?r?), hold off axis tight xlabel(?x?,?Fontsize?, 22) ylabel(?p_x?, ?Fontsize?, 22)%end% figure(2), plot(aa,UU)% figure(1)108Appendix C. Code For Chapter Twosubplot(212)xx = [x? x(end:-1:1)?]; YY = [Yp? Ym(end:-1:1)?]; fill(xx,YY,?g?)hold on,plot(x,Yp,?b?,x,Ym,?r?,x,Y,?k?,x,x./x,?LineWidth?,2),axis([-pi/2pi/2 min(Y) 1]) xlabel(?x?, ?Fontsize?, 22) ylabel(?Y_+ and Y_-?,?Fontsize?, 22) hold offlegend(?Pseudo-plug?,?Y_+?,?Y_-?,?Y?)109Appendix C. Code For Chapter Two% ssrhsx=([0:N-1]+1/2)?/N*pi-pi/2; Y=a*sin(x);px = X(1:N); U=X(N+1); q=X(N+2);BB=B./abs(px); %BB=B./sqrt(px.^2+1e-6);Yp1 = (1+Y)/2-U./px./(1-Y); Ym1 = Yp1-2*BB; Is1 =(sign(Y-Yp1)+1)/2.*(1+sign(px))/2; Q1 =U/2*(1-Y)-px/12.*(1-Y).^3+Y;Yp2 = 1-sqrt(2*U./abs(px)); Ym2 = Yp2-2*BB; Is2 =(sign(Y-Ym2)+1)/2.*(sign(Yp2-Y)+1)/2.*(1+sign(px))/2; Q2 =px/6.*(1-Yp2).^3+Y;Ym5 = (1+Y)/2-U./px./(1-Y); Yp5 = Ym5+2*BB; Is5 =(sign(Ym5-1)+1)/2.*(1-sign(px))/2; Q5 =U/2*(1-Y)-px/12.*(1-Y).^3+Y;Ym4 = Y+sqrt(2*U./abs(px)); Yp4 = Ym4+2*BB; Is4 =(sign(Yp4-1)+1)/2.*(sign(1-Ym4)+1)/2.*(1-sign(px))/2; Q4 =U*(1-Y)+px/6.*(Ym4-Y).^3+Y;Yp3 = (1-Y.^2-2*U./px-4*(Y+BB).*BB)/2./(1-Y-2*BB); Ym3 = Yp3-2*BB;Is3 = (sign(1-Yp3)+1)/2.*(sign(Ym3-Y)+1)/2; Q3 =px/6.*((1-Yp3).^3+(Ym3-Y).^3-3*(1-Y).*(Ym3-Y).^2)+Y;Yp = Yp1.*Is1+Yp2.*Is2+Yp3.*Is3+Yp4.*Is4+Yp5.*Is5; Ym =Ym1.*Is1+Ym2.*Is2+Ym3.*Is3+Ym4.*Is4+Ym5.*Is5; QQ =Q1.*Is1+Q2.*Is2+Q3.*Is3+Q4.*Is4+Q5.*Is5;tau1 = (2-Yp-Ym).*px/2;E=zeros(size(X)); E = QQ-q; E(N+1)=sum(tau1); E(N+2)=sum(px);EJ=zeros(N+2,N+2);110Appendix C. Code For Chapter Twofor n=1:N+2 XX = X; XX(n) = X(n)+1e-4; px = XX(1:N); U=XX(N+1);q=XX(N+2);BB=B./abs(px); %BB=B./sqrt(px.^2+1e-6);Yp1 = (1+Y)/2-U./px./(1-Y); Ym1 = Yp1-2*BB; Is1 =(sign(Y-Yp1)+1)/2.*(1+sign(px))/2;; Q1 =U/2*(1-Y)-px/12.*(1-Y).^3+Y;Yp2 = 1-sqrt(2*U./abs(px)); Ym2 = Yp2-2*BB; Is2 =(sign(Y-Ym2)+1)/2.*(sign(Yp2-Y)+1)/2.*(1+sign(px))/2; Q2 =px/6.*(1-Yp2).^3+Y;Ym5 = (1+Y)/2-U./px./(1-Y); Yp5 = Ym5+2*BB; Is5 =(sign(Ym5-1)+1)/2.*(1-sign(px))/2; Q5 =U/2*(1-Y)-px/12.*(1-Y).^3+Y;Ym4 = Y+sqrt(2*U./abs(px)); Yp4 = Ym4+2*BB; Is4 =(sign(Yp4-1)+1)/2.*(sign(1-Ym4)+1)/2.*(1-sign(px))/2; Q4 =U*(1-Y)+px/6.*(Ym4-Y).^3+Y;Yp3 = (1-Y.^2-2*U./px-4*(Y+BB).*BB)/2./(1-Y-2*BB); Ym3 = Yp3-2*BB;Is3 = (sign(1-Yp3)+1)/2.*(sign(Ym3-Y)+1)/2; Q3 =px/6.*((1-Yp3).^3+(Ym3-Y).^3-3*(1-Y).*(Ym3-Y).^2)+Y;Yp = Yp1.*Is1+Yp2.*Is2+Yp3.*Is3+Yp4.*Is4+Yp5.*Is5; Ym =Ym1.*Is1+Ym2.*Is2+Ym3.*Is3+Ym4.*Is4+Ym5.*Is5; QQ =Q1.*Is1+Q2.*Is2+Q3.*Is3+Q4.*Is4+Q5.*Is5;tau1 = (2-Yp-Ym).*px/2;EE=zeros(size(X)); EE = QQ-q; EE(N+1)=sum(tau1); EE(N+2)=sum(px);EJ(:,n)=(EE-E)*1e4;end111Appendix DCode For Chapter ThreeBelow is the MATLAB code used to obtain the results in Chapter Three.D.1 Newtonian Problem CodeFirst, we show the code used to solve the Newtonian problem.function sbvp global aa DD sol aa=0.001; DD = 1.5; U=aa^2/48 q=U;lambda(1)=q; lambda(2)=U; x = [0:100]/100*pi*2;options=bvpset(?Vectorized?,?on?,?RelTol?,1e-4,?AbsTol?,1e-7);soli = bvpinit(x,@solinit,lambda); sol =bvp4c(@solode,@solbc,soli,options);figure(1) plot(sol.x,1-sol.y(1,:)) sol.parametersqqq=sol.parameters(1); uuu=sol.parameters(2); aaa=aa; for na=1:25aa=na*4 sol = bvp4cc(@solode,@solbc,sol,options); holdon,plot(sol.x,1-sol.y(1,:)), hold off qqq=[qqq sol.parameters(1)];uuu=[uuu sol.parameters(2)]; aaa=[aaa aa]; end plateau =12*(1-qqq)./(12-6*uuu); figure(1) holdon,plot(sol.x,1-ones(1,length(sol.x))*plateau(na),?m?,?LineWidth?,2),hold off figure (1) title([?Swimmer Shape for D=?,num2str(DD),?and a=?, num2str(aa)],?FontSize?,22)xlabel(?\xi=x-t?,?FontSize?,22) ylabel(?Y?,?FontSize?,22) axistight axes(?position?,[0.2 0.2 0.28 0.28])plot(aaa,uuu,?r?,aaa,qqq,?g?) xlabel(?a?)axis tightfigure(2)plot3(1-sol.y(1,:),-sol.y(2,:),-sol.y(3,:))112Appendix D. Code For Chapter Threeview([-29,-64])xlabel(?Y?,?FontSize?,22)ylabel(?Y ??,?FontSize?,22)zlabel(?Y???,?FontSize?,22)title([?Phase Space for D = ?,num2str(DD)],?FontSize?,22)px = 6*(sol.parameters(2)-2)./(sol.y(1,:)).^2 -12*(sol.parameters(1)-1)./(sol.y(1,:)).^3;t1 = (4*sol.parameters(2)-6)./(sol.y(1,:)) -6*(sol.parameters(1)-1)./(sol.y(1,:)).^2;t0 = t1 - (sol.y(1,:)).*px;figure(3) plot(sol.x,t1,?k?,sol.x,t0,?c?)xlabel(?\xi=x-t?,?FontSize?,12) ylabel(?\tau_1: Black, \tau_0:Cyan?,?FontSize?,12) title(?Stresses\Strain Rates?, ?Fontsize?,12) axis tightfunction yinit = solinit(x) global aayinit = [ 1+aa/12*cos(x)-aa/12*sin(x)-aa/12*cos(x)aa/12*sin(x)aa/12*cos(x)zeros(size(x))zeros(size(x))];function res = solbc(ya,yb,lambda)res = [ ya(1)-yb(1)ya(2)-yb(2)ya(3)-yb(3)ya(4)-yb(4)ya(5)-yb(5)ya(6)ya(7)yb(6)-pi*2113Appendix D. Code For Chapter Threeyb(7) ];function dydx = solode(x,y,lambda)global aa DDq=lambda(1);U=lambda(2);dydx = [ y(2,:)y(3,:)y(4,:)y(5,:)-(aa*cos(x)+(12*(1-q)-y(1,:)*(12-6*U))./y(1,:).^3)/DDy(1,:)6*(1-q)./y(1,:).^2+(4*U-6)./y(1,:)];114Appendix D. Code For Chapter ThreeD.2 Viscoelastic Problem CodeNext, we list the code for the Viscoelastic problem with moderate amplitude.% siniCD =0.2; %7;%27;% Newtonianalf = 1; bet = 1; I1prime = @(g)g.^2.*(1+alf*g.^2).*(1-bet*g.^2+3*alf*g.^2+...alf*bet*g.^4)./(1+bet*g.^2).^3;I0prime = @(g)g.*(1-bet*g.^2+3*alf*g.^2+alf*bet*g.^4)./(1+bet*g.^2).^2;figure(1), subplot(221) sdrivC figure(1), subplot(221)figure(2), hold on, figure(1), hold on, axis tighttitle([?\alpha=\beta=1 and D=?,num2str(D)],?FontSize?,12)xlabel(?x?,?FontSize?,12)ylabel(?Y?,?FontSize?,12)% Shear Thinningalf = 1; bet = 2; I1prime = @(g)g.^2.*(1+alf*g.^2).*(1-bet*g.^2+3*alf*g.^2+...alf*bet*g.^4)./(1+bet*g.^2).^3;I0prime = @(g)g.*(1-bet*g.^2+3*alf*g.^2+alf*bet*g.^4)./(1+bet*g.^2).^2;figure(1), subplot(222) sdrivC figure(1), subplot(222)figure(2), hold on, figure(1), hold on, axis tighttitle([?\alpha=1, \beta=2, D=?,num2str(D)],?FontSize?,12)xlabel(?x?,?FontSize?,12)ylabel(?Y?,?FontSize?,12)% Shear Thickeningalf = 2; bet = 1; I1prime = @(g)g.^2.*(1+alf*g.^2).*(1-bet*g.^2+3*alf*g.^2+...alf*bet*g.^4)./(1+bet*g.^2).^3;I0prime = @(g)g.*(1-bet*g.^2+3*alf*g.^2+alf*bet*g.^4)./(1+bet*g.^2).^2;115Appendix D. Code For Chapter Threefigure(1), subplot(223) sdrivC figure(1), subplot(223)figure(1), hold on, axis tighttitle([?\alpha=2, \beta=1, D=?,num2str(D)],?FontSize?,12)xlabel(?x?,?FontSize?,12)ylabel(?Y?,?FontSize?,12)figure(1), subplot(224)%hold on,plot([0 DD],[0 DD].^2/48,?r:?),hold offxlabel(?a?,?FontSize?,12)ylabel(?U?,?FontSize?,12)title ([?U vs a for D=?,num2str(D)],?FontSize?,12)axis tight116Appendix D. Code For Chapter Three%sdrivCa=0.01 N=55; x=([0:N-1]+1/2)/N*pi*2;sdrivCC for niter=1:10srhsCdX=-EJ\E;conno=max(abs(dX));X=X+dX;if conno<1e-12disp(?Converging 1?)breakendend disp([?conno = ?, num2str(conno)])if (alf == 1) & (bet == 1)plot(x,Y,?r?)endif (alf == 1) & (bet == 2) %shear thinningplot(x,Y,?b?)endif (alf == 2) & (bet == 1) %shear thickeningplot(x,Y,?g?)end%%For comparison between Numerical & Analytical Solutions:%plot(x,Y)% FOR MEDIUM D SMALL A%hold on,%plot(x,D*a*sin(x)/(D^2+12^2)-a*12*cos(x)/(D^2+12^2),?k*?),%hold off% FOR LARGE D SMALL A%hold on,plot(x,D*a*sin(x)/(D^2+12^2),?k*?),hold off% FOR SMALL D SMALL A%hold on,plot(x,-a*12*cos(x)/(D^2+12^2),?k*?),hold offUU=X(3*N+1); qq=X(3*N+2); [UU qq] a0=a; aa=a0; D0=D; DD=D0; for na= 1:12a = a0 +na/10117Appendix D. Code For Chapter Threefor niter=1:15srhsCdX=-EJ\E;conno=max(abs(dX));X=X+dX;if conno<1e-8disp(?converging 2?)breakendenddisp([?conno = ?, num2str(conno)])UU=[UU X(3*N+1)];qq=[qq X(3*N+2)];aa = [aa a];hold onif (alf == 1) & (bet == 1)plot(x,Y,?r?), axis tightendif (alf == 1) & (bet == 2) %shear thinningplot(x,Y,?b?), axis tightendif (alf == 2) & (bet == 1) %shear thickeningplot(x,Y,?g?), axis tightendhold offend figure(1), subplot(224) if (alf == 1) & (bet == 1)figure(2),plot(aa,qq,?r?)figure(1), subplot(224),hold on, plot(aa,UU,?r?), hold offendif (alf == 1) & (bet == 2) %shear thinningfigure(2),plot(aa,qq,?b?)figure(1), subplot(224),hold on, plot(aa,UU,?b?), hold offendif (alf == 2) & (bet == 1) %shear thickeningfigure(2),plot(aa,qq,?g?)118Appendix D. Code For Chapter Threefigure(1), subplot(224), hold on, plot(aa,UU,?g?) , hold offendfigure(1)119Appendix D. Code For Chapter Three% srhsCgamY = X(1:N); gam1 = X(N+1:2*N); Y = X(2*N+1:3*N); U=X(3*N+1);q=X(3*N+2);x=([0:N-1]+1/2)?/N*pi*2;D3f = -a*cos(x);TY = gamY.*(1+alf*gamY.^2)./(1+bet*gamY.^2); T1 =gam1.*(1+alf*gam1.^2)./(1+bet*gam1.^2); TpY =(1+3*alf*gamY.^2)./(1+bet*gamY.^2)-...2*bet*gamY.^2.*(1+alf*gamY.^2)./(1+bet*gamY.^2).^2;Tp1 = (1+3*alf*gam1.^2)./(1+bet*gam1.^2)-...2*bet*gam1.^2.*(1+alf*gam1.^2)./(1+bet*gam1.^2).^2;I0 = gam1.^2/2-gamY.^2/2; I1 = gam1.^3/3-gamY.^3/3;px = (T1-TY)./(1-Y); pg1 = Tp1./(1-Y); pgY = -TpY./(1-Y);pY =(T1-TY)./(1-Y).^2;DiffYE=zeros(size(X)); EJ=zeros(3*N+2,3*N+2);% F4E(3*N+1)=sum(T1);% F5E(3*N+2)=sum(Y);for n=1:NI0(n)=quadl(I0prime,gamY(n),gam1(n));%F1E(n)=I0(n)-U*px(n); I1(n)=quadl(I1prime,gamY(n),gam1(n));%F2E(N+n)=I1(n)-px(n)^2.*(U+Y(n)-q-U*Y(n))-U*TY(n)*px(n);%F3E(2*N+n) = px(n)-D3f(n)-D*DY5(n);120Appendix D. Code For Chapter Three% DF1EJ(n,n)=-gamY(n)*TpY(n)-U*pgY(n); EJ(n,n+N)=gam1(n)*Tp1(n)-U*pg1(n);EJ(n,2*N+n)=-U*pY(n);EJ(n,3*N+1)=-px(n);% DF2EJ(N+n,n) =-TY(n)*TpY(n)*gamY(n)-2*px(n)*pgY(n)*(U+Y(n)-q-U*Y(n))-U*TpY(n)*px(n)-U*TY(n)*pgY(n);EJ(N+n,N+n) =T1(n)*Tp1(n)*gam1(n)-2*px(n)*pg1(n)*(U+Y(n)-q-U*Y(n))-U*TY(n)*pg1(n);EJ(N+n,2*N+n) =-px(n)^2*(1-U)-2*px(n)*pY(n)*(U+Y(n)-q-U*Y(n))- U*TY(n)*pY(n);EJ(N+n,3*N+1) = -px(n)^2.*(1-Y(n))-px(n)*TY(n);EJ(N+n,3*N+2) =px(n)^2;% DF3EJ(2*N+n,n) = pgY(n); EJ(2*N+n,N+n) = pg1(n);EJ(2*N+n,2*N+n) = pY(n); % Must add DY5 derivatives, see below% DF4EJ(3*N+1,N+n) = Tp1(n);% DF5EJ(3*N+2,2*N+n) = 1;end% Must do dF3/dY separately due to discretization% A is differentiation matrix for DY5EJ(2*N+1:3*N,2*N+1:3*N) =EJ(2*N+1:3*N,2*N+1:3*N)-D*A;121Appendix D. Code For Chapter ThreeThis code calculated the fifth derivative of a function given by a specific setof points numerically. It assumes periodicity of the function. It is assumed thatthe mesh used does not include its endpoints.%DiffYh = x(2)-x(1); %uniform meshA =5*diag(ones(N-1,1),1)-4*diag(ones(N-2,1),2)+diag(ones(N-3,1),3);A = A-A?;A(1,N-2) = -1; A(1,N-1) = 4; A(1,N) = -5; A(2,N-1) = -1; A(2,N) =4; A(3,N) = -1;A(N,1) = 5; A(N,2) = -4; A(N,3) = 1; A(N-1,1) = -4; A(N-1,2) = 1;A(N-2,1) =1;A = (1/(2*h^5))*A; DY5 = A*Y;122Appendix D. Code For Chapter ThreeThis is the code that formulates the initial guess for the Newton?s Iterationused to solve the viscoelastic problem.% srdrivCCx=([0:N-1]+1/2)/N*pi*2; px=-a*cos(x); gamY=-px/2; gam1=px/2;X=zeros(2*N+2,1); X(1:N)=gamY; X(N+1:2*N)=gam1;for niter=1:10 srhsCC dX=-EJ\E; conno=max(abs(dX)); X=X+dX; ifconno<1e-12 break end endfor niter=1:10 srhsCC dX=-EJ\E; conno=max(abs(dX)); X=X+dX; ifconno<1e-8 break end end% Updating X into a 4*N+2 length vector.X(1:N)=gamY; X(N+1:2*N)=gam1; X(2*N+1:3*N)=Y; X(3*N+1)=U;X(3*N+2)=q;123Appendix D. Code For Chapter Three% srhsCCgamY = X(1:N); gam1 = X(N+1:2*N); U=X(2*N+1); q=X(2*N+2);x=([0:N-1]+1/2)?/N*pi*2; px=-a*cos(x);TY = gamY.*(1+alf*gamY.^2)./(1+bet*gamY.^2); T1 =gam1.*(1+alf*gam1.^2)./(1+bet*gam1.^2); TpY =(1+3*alf*gamY.^2)./(1+bet*gamY.^2)-...2*bet*gamY.^2.*(1+alf*gamY.^2)./(1+bet*gamY.^2).^2;Tp1 = (1+3*alf*gam1.^2)./(1+bet*gam1.^2)-...2*bet*gam1.^2.*(1+alf*gam1.^2)./(1+bet*gam1.^2).^2;I0 = gam1.^2/2-gamY.^2/2; I1 = gam1.^3/3-gamY.^3/3;Y = 1-(T1-TY)./px; Yp1 = -Tp1./px; YpY = TpY./px;E=zeros(size(X)); EJ=zeros(2*N+2,2*N+2);E(2*N+1)=sum(T1); E(2*N+2)=sum(Y);for n=1:NI0(n)=quadl(I0prime,gamY(n),gam1(n));E(n)=I0(n)-U*px(n); I1(n)=quadl(I1prime,gamY(n),gam1(n));E(N+n)=I1(n)-px(n)^2.*(U+Y(n)-q-U*Y(n))-U*TY(n)*px(n);EJ(n,n)=-gamY(n)*TpY(n); EJ(n,n+N)= gam1(n)*Tp1(n);EJ(n,2*N+1)=-px(n); EJ(N+n,n) =-TY(n)*TpY(n)*gamY(n)-px(n)^2*(1-U)*YpY(n)-U*TpY(n)*px(n);EJ(N+n,N+n) = T1(n)*Tp1(n)*gam1(n)-px(n)^2*(1-U)*Yp1(n);EJ(N+n,2*N+1) = -TY(n)*px(n)-px(n)^2.*(1-Y(n)); EJ(N+n,2*N+2) =px(n)^2; EJ(2*N+1,N+n) = Tp1(n); EJ(2*N+2,n) = YpY(n);EJ(2*N+2,N+n) = Yp1(n); end124

Cite

Citation Scheme:

    

Usage Statistics

Country Views Downloads
United States 3 0
China 2 7
Canada 1 0
Philippines 1 0
City Views Downloads
Ashburn 2 0
Beijing 2 0
Unknown 1 0
Mountain View 1 0
Vancouver 1 0

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}
Download Stats

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>
                        
                    
IIIF logo 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-0066549/manifest

Comment

Related Items