@prefix vivo: . @prefix edm: . @prefix ns0: . @prefix dcterms: . @prefix dc: . @prefix skos: . vivo:departmentOrSchool "Science, Faculty of"@en, "Mathematics, Department of"@en ; edm:dataProvider "DSpace"@en ; ns0:degreeCampus "UBCV"@en ; dcterms:creator "Pachmann, Sydney"@en ; dcterms:issued "2008-08-26T16:50:33Z"@en, "2008"@en ; vivo:relatedDegree "Master of Science - MSc"@en ; ns0:degreeGrantor "University of British Columbia"@en ; dcterms:description """The purpose of this thesis is to study the problem of a low Reynolds number swimmer 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 length in one direction, creating a thrust and therefore propelling it in the opposite direction. We model the swimmer as an infinite, inextensible waving sheet. We consider two main cases of this swimming sheet problem. In the first case, the type of wave being propagated down the length of the swimmer is specified. We compare the swimming speeds of viscoelastic shear thinning, shear thickening and Newtonian fluids for a fixed propagating wave speed. We then compare the swimming speeds of these same fluids for a fixed rate of work per wavelength. In the latter situation, we find that a shear thinning fluid always yields the fastest swimming speed regardless of the amplitude of the propagating waves. We conclude that a shear thinning fluid is optimal for the swimmer. Analytical results are obtained for various limiting cases. Next, we consider the problem with a Bingham fluid. Yield surfaces and flow profiles are obtained. 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 a Newtonian fluid. Large amplitude forcing yields a swimmer shape that has a plateau region following by a large spike region. It is found that there exists an optimal forcing that will yield a maximum swimming speed. Next, we solve the problem for moderate forcing amplitudes for viscoelastic shear thickening and shear thinning fluids. For a given forcing, it is found that a shear thinning fluid yields the fastest swimming speed when compared to a shear thickening fluid and a Newtonian fluid. The difference in swimming speeds decreases as the bending stiffness of the swimmer increases."""@en ; edm:aggregatedCHO "https://circle.library.ubc.ca/rest/handle/2429/1503?expand=metadata"@en ; dcterms:extent "815483 bytes"@en ; dc:format "application/pdf"@en ; skos:note "Swimming in Slime by Sydney Pachmann B.Sc. First Class Honours, The University of Calgary, 2006 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF Master of Science in The Faculty of Graduate Studies (Mathematics) The University Of British Columbia (Vancouver) August, 2008 c© Sydney Pachmann 2008 Abstract The purpose of this thesis is to study the problem of a low Reynolds number swimmer 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 length in one direction, creating a thrust and therefore propelling it in the opposite direction. We model the swimmer as an infinite, inextensible waving sheet. We consider two main cases of this swimming sheet problem. In the first case, the type of wave being propagated down the length of the swimmer is specified. We compare the swimming speeds of viscoelastic shear thinning, shear thickening and Newtonian fluids for a fixed propagating wave speed. We then compare the swimming speeds of these same fluids for a fixed rate of work per wavelength. In the latter situation, we find that a shear thinning fluid always yields the fastest swimming speed regardless of the amplitude of the propagating waves. We conclude that a shear thinning fluid is optimal for the swimmer. Analytical results are obtained for various limiting cases. Next, we consider the problem with a Bingham fluid. Yield surfaces and flow profiles are obtained. 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 a Newtonian fluid. Large amplitude forcing yields a swimmer shape that has a plateau region following by a large spike region. It is found that there exists an optimal forcing that will yield a maximum swimming speed. Next, we solve the problem for moderate forcing amplitudes for viscoelastic shear thickening and shear thinning fluids. For a given forcing, it is found that a shear thinning fluid yields the fastest swimming speed when compared to a shear thickening fluid and a Newtonian fluid. The difference in swimming speeds decreases as the bending stiffness of the swimmer increases. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Overview of Non-Newtonian Fluids . . . . . . . . . . . . . . . . 5 1.2.1 Constitutive Models . . . . . . . . . . . . . . . . . . . . . 8 1.3 Our Problem: Main Assumptions . . . . . . . . . . . . . . . . . 11 1.3.1 Zero Force on Swimmer Condition . . . . . . . . . . . . . 15 1.4 Preliminary Problem . . . . . . . . . . . . . . . . . . . . . . . . 17 1.4.1 Mathematical Formulation . . . . . . . . . . . . . . . . . 17 1.4.2 Stream Function and Asymptotics . . . . . . . . . . . . . 19 1.4.3 Leading Order Equations . . . . . . . . . . . . . . . . . . 20 1.4.4 Second Order Equations . . . . . . . . . . . . . . . . . . 21 1.4.5 Limiting Cases . . . . . . . . . . . . . . . . . . . . . . . . 23 1.5 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2 Swimming with Lubrication Theory . . . . . . . . . . . . . . . . 27 2.1 Mathematical Formulation . . . . . . . . . . . . . . . . . . . . . 27 2.2 Newtonian Problem . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.3 Viscoelastic Problem . . . . . . . . . . . . . . . . . . . . . . . . 33 2.3.1 Effective Viscosity . . . . . . . . . . . . . . . . . . . . . . 34 iii Table of Contents 2.3.2 Numerical Formulation . . . . . . . . . . . . . . . . . . . 36 2.3.3 Fixed Dimensional Work . . . . . . . . . . . . . . . . . . 39 2.3.4 Limiting Cases . . . . . . . . . . . . . . . . . . . . . . . . 41 2.4 Bingham Problem . . . . . . . . . . . . . . . . . . . . . . . . . . 43 2.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3 Elastic Swimming with Lubrication Theory . . . . . . . . . . . 55 3.1 Mathematical Formulation . . . . . . . . . . . . . . . . . . . . . 56 3.2 Newtonian Problem . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.2.1 Numerical Formulation . . . . . . . . . . . . . . . . . . . 61 3.2.2 Low Amplitude Solution . . . . . . . . . . . . . . . . . . 64 3.3 Viscoelastic Results . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.3.1 Low Amplitude Solution . . . . . . . . . . . . . . . . . . 69 3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.1 Summary of Results . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.2 Further Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 Appendices A Full Analysis of Parent Problem . . . . . . . . . . . . . . . . . . 81 A.1 Mathematical Formulation . . . . . . . . . . . . . . . . . . . . . 81 A.2 Stream Function and Asymptotics . . . . . . . . . . . . . . . . . 83 A.3 Leading Order Equations . . . . . . . . . . . . . . . . . . . . . . 85 A.4 Second Order Equations . . . . . . . . . . . . . . . . . . . . . . 89 A.5 Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 A.6 Limiting Cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 A.6.1 Distance Between Wall and Swimmer Approaches Infinity 94 A.6.2 Swimmer Approaches the Wall . . . . . . . . . . . . . . . 95 B Details of Bingham Analysis . . . . . . . . . . . . . . . . . . . . . 97 B.1 Case I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 B.2 Case II . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 B.3 Case III . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 B.4 Case IV . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 iv Table of Contents B.5 Case V . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 C Code For Chapter Two . . . . . . . . . . . . . . . . . . . . . . . . 102 C.1 Viscoelastic Problem Code . . . . . . . . . . . . . . . . . . . . . 102 C.2 Bingham Problem Code . . . . . . . . . . . . . . . . . . . . . . . 108 D Code For Chapter Three . . . . . . . . . . . . . . . . . . . . . . . 112 D.1 Newtonian Problem Code . . . . . . . . . . . . . . . . . . . . . . 112 D.2 Viscoelastic Problem Code . . . . . . . . . . . . . . . . . . . . . 115 v List of Tables 2.1 Fluid types and the corresponding parameter values. . . . . . . . 35 2.2 The specific parameter values chosen for the various fluid types. . 36 2.3 Bingham Fluid Flow Profile Regions. In the outer regions, we see parabolic flow profiles typically characterized by Newtonian flows. In the inner region, we see what we refer to as pseudo-plug flow. In this region the velocity profile is constant. . . . . . . . . 45 2.4 Location of Yield Surfaces. . . . . . . . . . . . . . . . . . . . . . 47 vi List of Figures 1.1 Geometry of Problem . . . . . . . . . . . . . . . . . . . . . . . . 12 1.2 Geometry of Problem in Moving Frame of Reference . . . . . . . 14 2.1 Lubrication geometry . . . . . . . . . . . . . . . . . . . . . . . . 27 2.2 Effective Viscosities. Red: Newtonian, Solid Blue: Shear Thin- ning, Dashed Blue: Unrealistic Shear Thinning, Solid Green: Shear Thickening, Dashed Green: Unrealistic Shear Thickening . 35 2.3 Non-Dimensional Swimming Speeds. Red: Newtonian. Blue: Shearing Thinning. Green: Shear Thickening. . . . . . . . . . . . 39 2.4 Dimensional Swimming Speeds. Red: Newtonian. Blue: Shear- ing Thinning. Green: Shear Thickening. . . . . . . . . . . . . . . 41 2.5 Newtonian Flow Profiles for ah = 0.5. The speed of the swimmer, U , is depicted by the vertical red line. . . . . . . . . . . . . . . . 46 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. The surface of the swimmer is repre- sented by the black curve and the pink vertical line shows where the pressure gradient changes sign. . . . . . . . . . . . . . . . . . 50 2.7 Characteristic Bingham Velocity Profiles for B = 1 and ah = 1/2. Blue: Yielded Region. Red: Pseudo-plug Region. . . . . . . . . . 51 2.8 Pressure Gradients and Yield Surfaces for large and small B, a h = 1/2. In the left figures, the Newtonian pressure gradient is blue, the Bingham red. In the right figures, Y+ is blue, Y− is red, and the shaded green represents the pseudo-plugs. The surface of the swimmer is represented by the black curve and the pink vertical line shows where the pressure gradient changes sign. . . . 52 2.9 Swimming speed and mass flux for varying Bingham number. The corresponding Newtonian values are shown in red. . . . . . . 53 vii List of Figures 2.10 Example of alternate wave shape and resulting swimming speeds. Red: Newtonian. Blue: Shearing Thinning. Green: Shear Thick- ening. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.1 Geometry of Lubrication Problem. This time, the swimmer shape in unknown as well as the swimming speed . . . . . . . . . . . . 56 3.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 values for each case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.3 Phase portraits for Newtonian Solution of Elastic Lubrication Problem with large amplitude forcing. . . . . . . . . . . . . . . . 63 3.4 Solution for viscoelastic elastic lubrication problem for various fluid types with D=0.2. Red: Newtonian. Blue: Shearing Thin- ning. Green: Shear Thickening. . . . . . . . . . . . . . . . . . . . 67 3.5 Solution for viscoelastic elastic lubrication problem for various fluid types with D=7. Red: Newtonian. Blue: Shearing Thin- ning. Green: Shear Thickening. . . . . . . . . . . . . . . . . . . . 68 3.6 Solution for viscoelastic elastic lubrication problem for various fluid types with D=27. Red: Newtonian. Blue: Shearing Thin- ning. Green: Shear Thickening. . . . . . . . . . . . . . . . . . . . 69 3.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: Shearing Thinning. Green: Shear Thickening. . . . . . . . . . . . . . . . . 70 3.8 Optimal swimmer shape for D=0.2. . . . . . . . . . . . . . . . . . 71 3.9 Strain rates on the swimmer (cyan) and on the wall (black) for a=100 and D=1.5. Note that the largest strain rates occur in the region where the plateau ends and the spike begins. . . . . . . . 73 viii Acknowledgements I would like to thank my supervisors, Neil Balmforth and Daniel Coombs, for introducing me to the field of fluid dynamics and allowing me to work on this project. I am grateful for all the help they have provided. I would like to thank the other mathematical faculty at UBC from whom I took courses from. I feel I have learned a lot from you all. I would also like to thank another UBC graduate student, Samara Pillay, for her support and help throughout my time here at UBC. I would also like to make a special thanks to all the students in the IAM for just being great and making my experience here something truly fantastic. Finally, I express my sincere appreciation to all the helpful staff of the Mathematics Department and the IAM. This research was supported by an NSERC PGS M. ix Dedication To my parents, for always being there for me, my sister and brother, who try very hard find what I do interesting, and to Tyler, for putting up with me throughout. x Chapter 1 Introduction Years ago, the question of how microscopic organisms propel themselves through fluid was asked. The standard mechanism of propulsion for larger bodies was understood to rely on the shedding of vortices in an inviscid fluid. It was something entirely new to consider an organism propelling itself through fluid with a very small Reynolds number [37]. The mechanism of vortex shedding no longer 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 as well as the fact that the mechanisms of spermatozoa locomotion are significantly more amenable to mathematical modelling than other microorganisms. Many flagellated microorganisms are propelled by a rather complicated spiralling or corkscrew motion of their flagella. However, the flagellar motion of spermatozoa tends 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 as Taylor, Reynolds, Gray and Hancock, with continual contributions from leading researchers to this day [37, 33, 14, 15, 22, 25]. The purpose of this thesis is to study the problem of a low Reynolds number swimmer 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 length in one direction, creating a thrust and therefore propelling it in the opposite direction. We model the swimmer as an infinite, inextensible waving sheet, a standard approach in this field [37, 33, 22, 25]. Of course, realistic microorganisms have a body or head and flagellum of finite length and width, but this simplification to a swimming sheet provides useful insight into the general problem of swimming. This formulation of the problem also proves to have a wider range of applications [8, 26, 39, 2, 40, 18]. We consider two main cases of the swimming sheet problem. In the first problem, we suppose that the type of wave being propagated down the length 1 Chapter 1. Introduction of the swimmer is specified. In the second problem, the forcing along the length of 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 considering the lubrication approximation, introducing a wall and a non-Newtonian fluid and 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 preparation for the following discussions and mathematical analysis. We will then present the general problem under consideration, emphasizing the assumptions that are made and used in all following chapters, and the physical rationalizations for these assumptions. Finally, we will detail the problem of a two dimensional swimmer in the presence of a wall without the lubrication approximation to motivate the use of the thin gap theory. 1.1 Literature Review In the literature, there are several trains of thought about how to approach the low 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 swimming sheet, 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 that of a low amplitude wave. The solution is expanded as a Taylor series in powers of the non-dimensional wave amplitude [37, 33, 9, 22, 25]. A common result is that the leading order swimming speed is proportional to the amplitude of the swimmer 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 Newtonian fluid. The sheet was fixed to be a sinusoid, which he gave dimensionally by y0 = b sin(kx− σt), (1.1) where b is the amplitude, k the wave number, and σ is the frequency. No dynamical considerations were made for the sheet. He assumed that the waves were low amplitude and found that at first order, the swimming speed was zero. Continuing the calculation to second order, he found a non-zero swimming 2 Chapter 1. Introduction speed, which is henceforth referred to as the leading order swimming speed. This swimming speed, dimensionally, was calculated to be V c = 1 2 b2k2, (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. As is evident, this leading order swimming speed is proportional to the amplitude of the swimmer squared. In this paper, it was observed that real microscopic organisms are in contact with fluid on both sides. When adding a fluid to both sides of the swimmer, only the rate of dissipation of energy changed, from W to 2W . The speed of the swimmer did not. When referring to Taylor’s problem below, we are referring to the problem of an infinite, inextensible waving sheet with 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 commonly referred 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 resistive force theory is slender-body theory, which was introduced when incorporating the head of a microorganism into the analysis. It provides more accuracy in this instance than resistive force theory [16, 19]. In this swimming cylinder approach, the assumption of low amplitude waves is not necessarily made. However, when the assumption is made, the results agree qualitatively with those of the swim- ming sheet approach, that the swimming speed is proportional to the amplitude of 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 a point on the sheet surface and (xm, 0) be the mean position of (x0, y0), x0 = xm + a cos k[xm + (c− V )t] + d sin k[xm + (c− V )t], y0 = b sin k[xm + (c− V )t, ] (1.3) where c is the speed of the propagating wave and V the swimming speed. This wave 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 to the amplitude of the swimmer squared. Therefore, introducing this more general 3 Chapter 1. Introduction wave form did not change the result qualitatively from those that assumed no walls. In the limit as one wall of the channel goes to infinity, Katz’s results return the results of Blake, in [6], who considered a swimmer in the presence of one wall with a general wave form. One thing to note about this general wave form is that if the horizontal wave motion is large enough, the sheet can actually 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 close to the wall. It was found that the leading order swimming speed, V0, is larger than in the case with the wall being further away. The dimensional result was V0 c = 3 2 + ( h b )2 , (1.4) where h is the gap width. This implies that swimming speed is bounded by the speed of the propagating waves. For this lubrication problem, the velocities due to horizontal motion of the sheet are of a higher order and do not contribute to the leading order solution. Therefore, the sheet cannot reverse directions in this problem, as it could in the problem without the lubrication approximation applied. 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 a Newtonian fluid with walls that are both rigid and elastic. Here, the results were very comparable to real observations of sperm movement. The problem of a swimmer in a viscoelastic fluid has also been discussed in various contexts. Chaudhury, in [9], modified Taylor’s problem to include a second-order viscoelastic fluid. He found that viscoelasticity augments propul- sion in low Reynolds numbers, but as Reynolds number increases, there is a range in which the viscoelasticity hampers it. More recently, a paper written by Fu et al. [13] considers the swimming infinite cylinder problem using an upper convected Maxwell model as the constitutive equation for the viscoelastic fluid. They found that viscoelasticity tends to decrease the swimming speed and that swimming speed can reverse direction. This negative effect of viscoelasticity is in direct contrast with the work done by Chaudhury in the two dimensional problem. Lauga, in [25], considered Taylor’s problem in numerous different vis- coelastic fluids. For all cases, he found that the swimming speed was dependent on the viscoelasticity of the fluid and that it was always less that the correspond- 4 Chapter 1. Introduction ing Newtonian swimming speed. He also found that the result was unchanged if more general wave shapes, with both normal and tangential motion, were considered. In all the above viscoelastic problems, the assumption that the amplitude of the swimmer is small compared the wavelength of the swimmer was employed and so swimming speed is technically referring to leading order swimming 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 primary importance. Some related work has been done in this area. In [10], the authors considered a single elastic swimmer in a channel with a Newtonian fluid. They solved the problem using Peskin’s immersed boundary method. Their results were comparable to observed flagella and ciliary motion. There also exists some work done with elastic swimmers in slightly different situations. Lauga in [24] considered a lone, finite elastic swimmer in a Newtonian fluid. This problem of a low amplitude swimmer was solved using resistive force theory and analytical formulae were obtained. Work on finite elastic sheets in the presence of walls has been done as well, namely in the context of settling elastic sheets, or elastic sheets with one end clamped and oscillating vertically [2], [18]. Aside from the formulation of the force balances in these papers, particularly in [2], the results are relatively unrelated to our current discussions as they are for sheets and swimmers of finite length. In this thesis, we use analytic methods, rather than numerics so that we may gain more understanding of the problem, elucidate the dependence on the physical parameters, and determine the importance of the non-Newtonian medium. 1.2 Overview of Non-Newtonian Fluids When modelling fluids, we arm ourselves with equations that describe how the fluid flows. The Navier-Stokes equations governing Newtonian flow in vector form, ρ ( ∂u ∂t + u · ▽u ) = −▽ p+ µ▽2 u+ F. (1.5) They represent the conservation of momentum.1 Here ρ is the density of the fluid, u = (u, v) is the velocity field, p is the pressure within the fluid, µ is the 1Note that we are considering the two dimensional equations here, not the full three di- mensional ones. 5 Chapter 1. Introduction viscosity of the fluid and F represents any body forces acting. These equations must be supplemented with continuity or conservation of mass, ▽ · u = 0. (1.6) This is equivalent to assuming the fluid is incompressible or that density is constant. The question is, how do we describe the flow fluids that have behaviour that deviates from that of Newtonian fluid? To begin with, we must consider the equations of motion in a more general form. Still assuming an incompressible fluid, the conservation of mass equation remains unchanged. However, let us write the conservation of momentum equa- tions in a more general manner, referred to as Cauchy’s equations of motion: ρ ( ∂u ∂t + u · ▽u ) = −▽ p+▽ · τ + F. (1.7) Here τ is the deviatoric stress tensor. In the above Navier-Stokes equations for a Newtonian fluid, the relationship between the deviatoric stress and rate of strain is assumed to be τ = µγ̇, (1.8) thus, giving the µ▽2 u term in equation (1.5). This constitutive relationship between stress and rate of strain is unable to describe the properties of non- Newtonian fluids. Let us put for now the constitutive relation between stress and rate of strain in a general form: f(τ , γ̇) = 0. (1.9) To get an idea of how a non-Newtonian fluid can differ from a Newtonian fluid, let us look at some standard examples of non-Newtonian fluids and related behaviours exhibited by these fluids. To begin with, let us consider a cornstarch solution. This is an interesting example of what is called a shear thickening fluid.2 Anyone that has ever made up their own home solution of cornstarch and water knows that when you play with it, in particular if you hit it or smack it with your palm, it immediately gets hard like a solid. But, when you move your hand through it gently, it behaves again like a liquid. Shear thickening means that the viscosity of the 2It can be shown that cornstarch solutions have behaviour that is more complicated that that which would be expected of a simple shear thickening fluid [3] 6 Chapter 1. Introduction fluid increases with the shear rate, or deformation rate. Therefore, when you smack 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 in fairly concentrated suspensions of small particles. Shear thickening is not the only interesting behaviour that non-Newtonian fluids can exhibit. The opposite effect has also been observed, that of shear thinning. In this instance, the viscosity decreases as the strain rate increases. Examples of fluids that exhibit shear thinning are molten polyethylene and polypropylene, as well as solutions of caryboxymethylcellulose in water and many others. There also exists what are called viscoplastic fluids. These fluids exhibit a yield 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 exhibited by 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 a cylinder. Consider a viscoelastic fluid that starts at rest and is then subject to a pressure gradient, causing it to flow. This pressure gradient is then removed. What is observed is that the fluid actually recoils or moves backwards towards its 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 strong recoil, while with others it is less prominent. It is easy to see why this is often referred to as fading memory. A fully elastic material would return fully to its original state. Two common examples of phenomena observed in viscoelastic fluids are jet swell and rod climbing. After exiting an orifice in a downward direction, a Newtonian fluid has a stream of a set diameter that is the same as that of the jet. In comparison, viscoelastic fluids have a diameter that increases twofold from the orifice diameter. This is also called the Barus or Merrington effect. Rod climbing, or the Weissenburg effect, refers to the situation where a thin solid cylinder is rotating in a viscoelastic fluid. It is observed that the fluid climbs up the rod, moving in the opposite direction to a Newtonian fluid in the same situation. A Newtonian fluid in this situation, besides moving in a different direction, does not climb the rod. In fact, a small indentation at the 7 Chapter 1. Introduction base of the rod is observed. Due to the composition of cervical mucous, namely that it is composed of various 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 a polymeric fluid are very dependent on the molecular weight distribution of the macromolecules that comprise it. That is, the rheological properties of a poly- meric fluid depend not only on the number of molecules composing it, but the sizes and locations of those molecules as well. Polymeric fluids can exhibit many of the behaviours discussed above, including shear thinning and thickening, jet swell and rod climbing, and, obviously, viscoelasticity.5 1.2.1 Constitutive Models There are many different ways to model non-Newtonian fluids ranging from relatively simple, to immensely complex, depending on the properties of the specific fluid one is interested in. On the simpler end of the scale there are generalized Newtonian fluid models and, for viscoelastic fluids, general linear viscoelastic fluid models. On the complex end of things there are corotational and codeformational models. To begin with, let us look at generalized Newtonian fluids. This is the simple case when you take the constitutive model for a Newtonian fluid, as seen in equation (1.8), and relax the assumption that the viscosity, µ, is constant. The constant viscosity is replaced with a viscosity that varies with the shear rate, γ̇, or even sometimes with |τxy|. A standard of this type of constitutive law 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 viscoplastic 3Cervical mucous is a gel that consists of a network of glycoprotein macromolecules and a plasma 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 dilute solutions, the macromolecules are far apart from each other and therefore have negligible influence on each other. 5For a full, in depth discussion on polymeric fluids and the associated phenomena, see [5] 8 Chapter 1. Introduction fluid is a Bingham fluid, which is modelled by µ = ∞, τ ≤ τy, µ = µ0 + τy γ̇ , otherwise. (1.11) Another example of a viscoplastic fluid is a Herschel-Bulkley fluid, which also incorporates shear thinning, given by µ = ∞, τ ≤ τy, µ = (τy + µ|γ̇| n)sgn(γ̇), otherwise. (1.12) Generalized linear viscoelastic fluid models are simple models of viscoelastic fluids that display some of the time dependent properties (i.e. recoil or viscoelas- ticity). These, however, are only applicable to situations that have low strains and shear rates. The most standard example of a generalized linear viscoelastic model is Maxwell’s model τ + λ ∂ ∂t τ = µγ̇, (1.13) where λ is a time constant and µ has units of viscosity. Another generalized linear viscoelastic model is the Jeffreys model, τ + λ1 ∂ ∂t τ = µ ( γ̇ + ∂ ∂t γ̇ ) . (1.14) As you can see, this is just Maxwell’s model with an extra term added in on the right-hand side, a time derivative of the strain rate and another time constant, λ2. The Jeffreys model is important because it has been the starting point for other, more complicated rheological equations of state. Both of these types of models are limited in scope. Ideally, we want a model that describes all flow phenomena, at least qualitatively. Two types of models that do this are corotational models and codeformational models. These are basically just versions of the linear viscoelastic models revised so that they are valid in situations when the fluid gets largely deformed from its initial position. Corotational models are formulated in a reference frame that corotates with the fluid. That is, the coordinates, call them θj , have a time rate of change given by ∂ ∂t θj = − 1 2 ω · θj , (1.15) where ω = ▽× u is the fluid vorticity. Observers in this frame will not detect 9 Chapter 1. Introduction rigid rotation. Codeformational models are formulated in codeforming coordi- nates, call them φj , with time rate of change ∂ ∂t φj = −φj · ▽u. (1.16) The models being used in this thesis are codeformational. For a more in depth discussion 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 both codeformational and corotational and are based on the Jeffreys model. We will restrict our discussion here to the codeformational type, since that is the type that will be used in this thesis. To begin with, let us look at the Oldroyd-B model. This model is in fact the Jeffreys model generalized in codeformational coordinates.6 It is given by, τ + λ1τ ▽ = µ ( γ̇ + λ2γ̇ ▽ ) . (1.17) Here, λ1 is the polymer relaxation time, and λ2, another relaxation time, rep- resents the rate of decay of the residual rate of strain when the fluid is stress free, or the retardation time scale. The upper convected derivative of a tensor A, is defined to be A▽ = ∂A ∂t + u · ▽A− ( (▽u)T ·A+A · ▽u ) . (1.18) It is possible to eliminate one relaxation time and write it in terms of the other by thinking of the fluid as being composed of a Newtonian solvent and a polymeric solute. This allows one to write the viscosities and deviatoric stresses as sums. We put µ = µs + µp and τ = τ s + τ p. In this situation, it is generally assumed that µ = λ1λ2µs. This assumption allows us to remove one relaxation number and also allows us to get rid of the convected strain rate term. Applying these assumptions to the above Oldroyd-B constitutive equation simplifies it to the upper convected Maxwell model for the polymeric component, τ p + λ1τ ▽ p = µpγ̇. (1.19) 6Technically, there are two possible generalizations of the Jeffreys model in codeformational coordinates, the Oldroyd-A model and the Oldroyd-B model. Oldroyd-A is very similar to Oldroyd-B, with the difference being in the convected derivative. Oldroyd-A is given by τ + λ1 ( Dτ Dt + τ · (▽u)T +▽u · τ ) = µ ( γ̇ + λ2 ( Dγ̇ Dt + γ̇ · (▽u)T +▽u · γ̇ )) . 10 Chapter 1. Introduction Another Oldroyd model is the Oldroyd-8 model, aptly named as there are eight parameters in it. This model is also a generalization of the Jeffreys model, in which all possible quadratically nonlinear terms involving the deformation rates allowed by the symmetries of the problem are included. This model is given by, τ + λ1τ ▽ + 1 2 µ0(trτ )γ̇ − 1 2 µ1 (τ · γ̇ + γ̇ · τ ) + 1 2 ν1(τ : γ̇)δ = µ ( γ̇ + λ2γ̇ ▽ − µ2(γ̇ · γ̇) + 1 2 (γ̇ : γ̇)δ ) . (1.20) The eight constants and extra terms in this model allow for a much larger variety in rheological response compared to the Jeffreys model and the Oldroyd- B model [5]. This model is also convenient as it contains several other models with in it that can be obtained by merely setting some of the constants equal to zero, or multiples of each other. For example, it contains the Oldroyd-B model by setting all constants to zero except λ1 and λ2. It also contains the upper convected Maxwell model, the second-order fluid model, and the Oldroyd four and six constant models, to name a few. Codeformational constitutive models can get significantly more complicated than those discussed above. However, a full discussion these rheological equa- tions of state is not necessary here, as the only codeformational models we will be using are Oldroyd-B and Oldroyd-8 constitutive models. 1.3 Our Problem: Main Assumptions Now that we are somewhat comfortable with viscoelastic fluids, we are in a position to set up our main problem mathematically. Consider a two-dimensional, infinite periodic swimmer in the presence of a wall. We would like to consider the case when the fluid between the wall and the swimmer is viscoelastic. This is because the physical properties of cervical mucous are viscoelastic. We are assuming that the swimmer propagates waves down its length therefore propelling itself in the opposite direction. The shape of the swimmer is referred to as Y and is centered about y = 0. The wall is fixed 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 〈Yx〉 = 0 and we assume that 11 Chapter 1. Introduction 〈Y 〉 = 0, so that h is the mean gap width. Here 〈. . .〉 = 1 2pi ∫ 2pi/k 0 (...)dx. (1.21) Figure 1.1: Geometry of Problem The equations governing such a flow are ▽ ·u = 0, (1.22) ρ ( ∂u ∂t + u · ▽u ) = −▽ p+▽ · τ + F, (1.23) f(τ , γ̇) = 0. (1.24) Recall that the constitutive relation between stress and rate of strain in a general form. It is assumed that the swimmer propagates waves down its length in the positive x-direction, therefore propelling itself in the negative x-direction. The swimming speed is referred to as U , thus its velocity is −U . We assume that U > 0. Hence, we have boundary conditions as follows. For the horizontal velocity component, u(x, h) = 0, (1.25) u(x, Y ) = −U, (1.26) and for the vertical velocity component, assuming simple vertical motion rather 12 Chapter 1. Introduction than 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 make several simplifying assumptions to help ourselves along. Due to the physical situation under consideration, a low Reynolds number regime is considered 0 < Re≪ 1, (1.29) allowing the inertial terms or the left hand side of Cauchy’s equation of motion to be neglected. It will be assumed that there are no external forces acting in the problem F = 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 small compared 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 formulated mathematically by stating that 0 < ( maxx∈[0,2pi]|Y | ) k = ε≪ 1, (1.32) where k is the wavenumber associated with Y . It is also assumed that the problem in its entirety is periodic in (kx − (ω − Uk)t). Of specific use will be the periodicity of the pressure, 〈px〉 = 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 is small 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, more asymmetrical waves [36]. 13 Chapter 1. Introduction That is, put x = x̃− Ut, (1.33) dx dt = dx̃ dt − U. (1.34) This implies that the solutions will now be periodic in (kx̃ − ωt). For the remainder of this thesis, we will assume that it is understood we are in the moving frame, and so will drop the tildes on the above changed variables. Figure 1.2: Geometry of Problem in Moving Frame of Reference This simplifies the boundary conditions. In the new frame of reference, for the horizontal velocity component, we have u(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 the new frame of reference have been specified. 14 Chapter 1. Introduction 1.3.1 Zero Force on Swimmer Condition In addition to the above kinematical assumptions, we will also introduce a dy- namical one. We assume that the average force per wavelength on the swimmer is zero. This is a reasonable assumption due to the fact that it has already been assumed 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 minus the total force on the fluid by the swimmer (per unit length). Thus, we can formulate this condition by considering the total force on the fluid. It known that the total stress tensor, σ, is given by σ = −pI+ τ = ( −p+ τxx τxy τyx −p+ τyy ) , and it is found that the unit outward normal to the swimmer is 1√ 1 + Y 2x ( Yx −1 ) . We can write the total force on the fluid per unit length as ∫ 2pi 0 ( −p+ τxx τxy τyx −p+ τyy )( Yx −1 ) dx. Do to the force balance, we know that this is equal to minus the force on the swimmer, which is assumed to be zero. Thus, by writing things in component form, one can see that (recall the stress tensor is symmetric) 〈(−p+ τxx)Yx − τxy〉 = 0, 〈τxyYx − (−p+ τyy)〉 = 0, where the top equation is the tangential force, the bottom the normal force, recalling that 〈. . .〉 denotes the x-average. Now, consider the conservation of momentum equations in component form, slightly rearranged (p− τxx)x = τxy,y, (p+ τyy)y = τxy,x. 15 Chapter 1. Introduction Concerning ourselves with only the first equation 9 and by applying the average in x while maintaining the assumed periodicity of the problem, we discover 〈τxy〉y = 0, and by integrating 〈τxy〉 = Const. Now, consider the same momentum equation. Rather than applying an average in x, integrate in y. This gives ∫ hk Y (p− τxx)xdy = τxy(hk)− τxy(Y ). Pulling the x-derivative out of the integral, ∂x ∫ hk Y (p− τxx)dy + (p− τxx) ∣∣ y=Y Yx = τxy(hk)− τxy(Y ), and again average in x to obtain, 〈τxy(hk)〉 = 〈τxy(Y )〉+ 〈(p− τxx) ∣∣ y=Y Yx〉. 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 force is equal to zero leaves 〈τxy(hk)〉 = 0, and since it was calculated above that 〈τxy〉 = Const, we finally have 〈τxy〉 = 0. (1.39) This is a constraint on the system that will be used in both Chapters Two and Three. 9At this time we will ignore the vertical force balance due to the fact that either the swimmer is assumed to be fixed vertically by some external force or a more complex vertical force balance is required utilizing elasticity theory, as will be seen in Chapter Three. 16 Chapter 1. Introduction 1.4 Preliminary Problem In this section, we discuss the fundamental problem of a swimmer in the presence of a wall and a viscoelastic fluid. We do not apply the lubrication approximation here. We will assume that the amplitude of the waves is small in comparison to the mean distance to the wall and that the fluid between the swimmer and the wall is an Oldroyd-B viscoelastic fluid. This is a nice warm up problem as it is a combination of Katz’s biharmonic problem and Lauga’s viscoelastic version of the Taylor problem [22, 25]. The method of solution applied here follows closely with that used in [25]. See figure (1.2) for geometry. 1.4.1 Mathematical Formulation We assume that the sheet has prescribed motion, that of waves propagating to the right, see figure (1.2). In the frame moving with the swimmer, these waves look like y = a sin(kx− ωt). (1.40) The equations of motion for this problem are ▽ ·u = 0, (1.41) ▽p = ▽ · τ , (1.42) τ + λ1τ ▽ = µ ( γ̇ + λ2γ̇ ▽ ) , (1.43) which are the conservation of mass, momentum and the Oldroyd-B constitutive equation respectively. 10 Recall that the upper convected derivative of a tensor A, is defined to be A▽ = ∂A ∂t + u · ▽A− ( (▽u)T ·A+A · ▽u ) . (1.44) We also have the zero force condition11 〈τxy〉 = 0. (1.45) 10Recall that λ1 is the polymer relaxation time, and λ2 represents the rate of decay of the residual 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 to the low amplitude compared to mean height approximation. 17 Chapter 1. Introduction Non-dimensionalizing the equations as follows t = 1ω t̂, (x, y) = 1 k (x̂, ŷ), u = cû, γ̇ = ωˆ̇γ, (p, τ ) = µω(p̂, τ̂ ), Û = Uc where, c = ωk , is the speed of the travelling wave, and dropping the hat notation, yields ▽ p = ▽ · τ , (1.46) ▽ · u = 0, (1.47) which are unchanged, and τ +De1τ ▽ = γ̇ +De2γ̇ ▽, (1.48) the constitutive equation, where De1 and De2 are Deborah numbers. The Deborah number is the ratio of the characteristic time for the fluid over the characteristic time for the flow system. Here, we have two characteristic times for the fluid, and hence two Deborah numbers. It remains to non-dimensionalize the boundary conditions. Using the low amplitude approximation, we put ε = ak ≪ 1, (1.49) and find that (dropping hats): y = ε sin(x− t). (1.50) The boundary conditions become, for the horizontal velocity component u(x, hk) = U, (1.51) u(x, Y ) = 0, (1.52) and for the vertical velocity component v(x, hk) = 0, (1.53) v(x, Y ) = −ε cos(x− t). (1.54) 18 Chapter 1. Introduction 1.4.2 Stream Function and Asymptotics Since the problem at hand is an incompressible, two-dimensional problem, it is most useful to consider a stream function ψ such that 12 u = ∂ψ∂y , v = − ∂ψ ∂x . (1.55) To solve this system, consider low amplitude solutions of the form U = εU1 + ε 2U2 + ... (1.56) ψ = εψ1 + ε 2ψ2 + ... (1.57) τ = ετ1 + ε 2τ2 + ... (1.58) That is, regular perturbation expansions in ε for the swimming speed, the stream function and the stress. Since it is assumed that the speed of the swimmer is constant, it must be the case that each Uk in the expansion of U is also constant. Notice that, since the velocity field components can be written in terms of the stream function and γ̇ can be written in terms of the velocity field components, we can put γ̇ = εγ̇1 + ε 2γ̇2 + ... (1.59) u = εu1 + ε 2u2 + ... (1.60) Removing the pressure term from the conservation of momentum equations and then summarizing the equations of motion gives ▽×▽ ·τ = 0, (1.61) τ +De1τ ▽ = γ̇ +De2γ̇ ▽. (1.62) The next step is to then substitute in the perturbation expansions and equate orders of ε. 12Conservation of mass holds automatically since ▽ · u = ux + vy = ∂ ∂x ( ∂ψ ∂y ) + ∂ ∂y ( − ∂ψ ∂x ) = 0. 19 Chapter 1. Introduction 1.4.3 Leading Order Equations At O(ε) we have: ▽×▽ ·τ1 = 0, (1.63) τ1 +De1 ∂ ∂t τ1 = γ̇1 +De2 ∂ ∂t γ̇1. (1.64) We find the boundary conditions at this order by substituting in the regular perturbation expansion for ψ and expanding about y = 0. This gives ψ1x ∣∣ (x,0) = 0, ψ1y ∣∣ (x,0) = 0, ψ1x ∣∣ (x,hk) = cos(x− t), ψ1y ∣∣ (x,hk) = U1. (1.65) The zero force constraint 〈τxy〉 = 0 yields that 〈τ1xy〉 = 0. (1.66) If one stares at these leading order equations for a moment, it becomes apparent that the easiest thing to do is to take the divergence and then the curl of equation (1.64),thus eliminating leading order stress. By doing this we obtain ( 1 +De2 ∂ ∂t ) ▽4 ψ1 = 0, (1.67) where ▽4ψ = ψxxxx + 2ψxxyy + ψyyyy. To solve this equation, we simply use separation of variables with hyperbolic sines and cosines. Assume that the form of the solution is ψ1(y) = [(A+BH) cosh(H) + (C +DH) sinh(H)] sin(x− t), where H = hk − y. We find then that U1 = 0, (1.68) and ψ1 = W (y) W (0) sin(x− t), (1.69) 20 Chapter 1. Introduction where Ψ1(y) = (hk − y) sinh(hk − y), (1.70) Ψ2(y) = sinh(hk − y)− (hk − y) cosh(hk − y), (1.71) and W (y) = Ψ′1(0)Ψ2(y)−Ψ ′ 2(0)Ψ1(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 calculate the 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 portion of a larger, complex problem that satisfies the same equations.13 We put τ1 = ℜ[τ̃1e i(x−t)], γ̇1 = ℜ[ ˜̇γ1e i(x−t)]. (1.73) By employing this method of solution, we can easily obtain τ̃1 = 1− iDe2 1− iDe1 ˜̇γ1, (1.74) and upon substituting in ˜̇γ1, we get τ1 = ℜ { (1− iDe2)e i(x−t) (1− iDe1)W (0) ( 2W ′ −i(W ′′ +W ) −i(W ′′ +W ) −2W ′ )} . (1.75) Notice that the condition, 〈τ1xy〉 = 0, automatically holds due to the fact that τ1xy ∝ γ̇1xy ∝ sin(x− t). 1.4.4 Second Order Equations The leading order problem has been solved and it has been shown that in order to maintain the assumption that the speed of the swimmer is constant, at leading order it must be equal to zero. Thus, in order to calculate the swimming speed of the organism, we must move on to the second order problem. This problem is 13Note the periodic dependence explicitly given due to the assumed periodicity of the prob- lem 21 Chapter 1. Introduction not quite as straightforward as the leading order problem due to the nonlinear terms of the upper convected derivative. At O(ε2) we have ▽×▽ ·τ2 = 0, (1.76) τ2 +De1 ∂ ∂t τ2 +De1 ( u1 · ▽τ1 − τ1 · ▽u1 − (▽u1) T · τ1 ) = γ̇2 +De2 ∂ ∂t γ̇2 +De2 ( u1 · ▽γ̇1 − γ̇1 · ▽u1 − (▽u1) T · γ̇1 ) , (1.77) with boundary conditions ψ2x ∣∣ (x,0) = 0, ψ2y ∣∣ (x,0) = − sin2(x− t)W ′′(0) W (0) , ψ2x ∣∣ (x,hk) = 0, ψ2y ∣∣ (x,hk) = U2. (1.78) The zero force condition states that 〈τ2xy〉 = 0. (1.79) Inspecting the second order equations of motion reveals that the right-hand side 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-component ( 1 +De1 ∂ ∂t ) τ2xy − ( 1 +De2 ∂ ∂t ) γ̇2xy = ℜ { (De2 −De1) 2(1− iDe1) [ ũ∗1 · ▽ ˜̇γ1 − ˜̇γ1 · ▽ũ ∗ 1 − (▽ũ ∗ 1) T · ˜̇γ1 ]} xy + +ℜ { (De2 −De1) 2(1− iDe1) e2i(x−t) [ ũ1 · ▽ ˜̇γ1 − ˜̇γ1 · ▽ũ1 − (▽ũ1) T · ˜̇γ1 ]} xy (1.80) This allows us to utilize equation (1.79) as a simplification. To apply this condition, let us take the x-average of equation (1.80). Then, by writing the rate of strain component in terms of the stream function, and observing the fact that 〈ψ2xx〉 = 1 2piψ2x ∣∣2pi 0 = − 12piv ∣∣2pi 0 = 0 due to the assumed periodicity of the solution, we obtain 〈ψ2〉yy = ℜ { (De1 −De2) 2(1− iDe1) [ ũ∗1 · ▽ ˜̇γ1 − ˜̇γ1 · ▽ũ ∗ 1 − (▽ũ ∗ 1) T · ˜̇γ1 ] xy } . (1.81) 22 Chapter 1. Introduction To solve the problem and calculate the swimming speed we must integrate the above expression for 〈ψ2〉yy once to obtain an expression for 〈ψ2〉y and apply the boundary conditions. To be able to apply the boundary conditions, we must put them into the same averaged form as the equation. For the conditions on ψx, we obtain trivial averages due to the periodicity of the problem, 〈ψ2x〉 = 1 2piψ2 ∣∣2pi 0 = 0, and thus no information. However, for the other two boundary conditions we obtain 〈ψ2y〉 ∣∣ y=0 = − W ′′(0) 2W (0) , (1.82) 〈ψ2y〉 ∣∣ y=hk = U2. (1.83) By integrating we find that, writingW (0) =W0 andW ′′(0) =W ′′0 for simplicity, 〈ψ2〉y = De1(De2 −De1) 2(1 +De21)W 2 0 [ WW ′′ +W ′2 ] − W ′′0 2W0 (1 +De1De2) (1 +De21) . (1.84) Finally, the swimming speed is obtained, U2 = 〈ψ2〉y ∣∣ y=hk = − W ′′0 2W0 (1 +De1De2) (1 +De21) . (1.85) Note that W (hk) = 0 and W ′(hk) = 0. 1.4.5 Limiting Cases Distance Between Wall and Swimmer Approaches Infinity We would like to assure ourselves that this problem agrees with the infinite problem considered Lauga, in [25]. To do this, we take the limit as the distance between the wall and the swimmer approaches infinity. It is known that ψ1 = W (y) W (0) sin(x− t), and so the limit as h −→∞ is easily calculated. We find that lim h−→∞ W (y) W (0) = (1 + y)e−y. Thus, lim h−→∞ ψ1 = (1 + y)e −y sin(x− t), 23 Chapter 1. Introduction which is exactly the infinite solution as stated in [25]. It is obvious that if the limit of the streamfunction agrees with the infinite case that the limit of the swimming speed will as well, as ψ2y(hk) = U2. When calculated we find that lim h−→∞ U2 = 1 2 (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 by the swimmer. We find that lim h−→∞ { (1 +De1De2) (1 +De21)W 2 0 ∫ hk 0 [ 4W ′2 + (W +W ′′)2 ] dy } = 2 (1 +De1De2) (1 +De21) , which again agrees with the previous work done.14 So, it is clear that the above calculations for the problem of an organism swimming in an Oldroyd-B viscoelastic fluid near a wall are consistent with those of the problem of an organism swimming in an infinite fluid of the same type. Swimmer Approaches the Wall Next, we would like to consider the limit as the swimmer and the wall approach each other. There are two ways of thinking of this limit. The first is keeping the wall fixed at y = h and increasing the amplitude of the swimmer until it reaches the wall. This line of thinking contradicts the assumption that ak << 1, as it implies that O(a) → O(h) = O(1/k). This would make our present solution invalid. The second way of thinking about this limit is to think of the wall approaching the swimmer, whose amplitude and wavelength remains the fixed. This implies that O(h)→ O(a), still maintaining the assumption that ak << 1. In fact, in this manner of thinking we get that hk << 1 which is the same assumption that is applied when making the lubrication approximation. Taking this limit in this case, however, is not the same as applying the lubrication approximation fully. It is just a limiting case of the problem under discussion, and will suggest to us to fully apply the lubrication approximation to understand 14Please see Appendix A for a full derivation of this problem and the calculation of the work done. 24 Chapter 1. Introduction this limit better. We have then, to leading order of this limiting case, U = −ε2 (1 +De1De2) 2(1 +De21) F (hk), (1.86) where F (hk) = W ′′0 W0 = (hk)2 + sinh2(hk) (hk)2 − sinh2(hk) . (1.87) In our limit, hk → ak = ε and so we have that F (ε) = ε2 + sinh2 ε ε2 − sinh2 ε . (1.88) Expanding sinh ε = ε+ 16ε 3 + . . . and simplifying gives us F (ε) = 2ε2 ++ 13ε 4 + . . . − 13ε 4 − . . . . (1.89) Thus, in the limit when the wall approaches the swimmer, we find that the swimming speed is U = 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 see Appendix A. 1.5 Thesis Outline Chapter Two presents the problem of a infinite swimming sheet with the lubri- cation approximation applied. The sheet shape here is assumed. This problem is first solved for a Newtonian fluid in a preliminary sense, and is then solved for a viscoelastic fluid represented by the Oldroyd-8 constitutive model. Chapter Three presents the problem of a infinite swimming elastic sheet with the lu- brication approximation applied. In this problem, there is a prescribed forcing 15Note 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 that when the swimmer approaches the wall, the swimming speed is not longer proportional to the amplitude of the swimmer squared and something more interesting is happening. 25 Chapter 1. Introduction along the length of the sheet. This problem is solved for a Newtonian fluid for forcing of any strength, and then for a viscoelastic fluid, also represented by the Oldroyd-8 constitutive model, with forcing of moderate strength. The final chapter presents a summary of the findings and concluding remarks, including further work that may be done. 26 Chapter 2 Swimming with Lubrication Theory As motivated by the preliminary problem discussed above, we now consider the problem of a two dimensional swimmer very near a wall. Particularly, we apply the lubrication approximation and assume that the mean distance of the swimmer from the wall, h, compared to the wavelength of the swimmer is very small.16 Figure 2.1: Lubrication geometry 2.1 Mathematical Formulation Let 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 be y = a sin(kx− ωt). (2.1) We consider this wave form and not a more general one because it has been shown by Katz in [22] that for the Newtonian lubrication problem, horizontal contributions to the wave shape only affect the higher orders of the solution. For the viscoelastic problem of no wall, Lauga in [25] considered horizontal 16For a mathematical discussion of viscoelastic lubrication theory, we refer the reader to [4] 27 Chapter 2. Swimming with Lubrication Theory oscillations of the same order as the amplitude in the sheet and found that they did not affect the leading order solution. Thus, it is very reasonable to assume that horizontal contributions to the wave form will not affect our leading order solution, and so they will be neglected. Here, we are applying the lubrication approximation. Mathematically this means that hk ≪ 1. (2.2) The equations for this problem, using the general assumptions stated in Chapter One, are ▽ ·u = 0, (2.3) ▽p = ▽ · τ , (2.4) f(τ , γ̇) = 0, (2.5) which are the conservation of mass, momentum and the general constitutive equations respectively.17 This problem is also subject to the constraint that 〈τxy〉 = 0, (2.6) from the zero force condition. As well, we require that the periodicity in pressure be considered explicitly for this problem, thus we will include 〈px〉 = 0, (2.7) in our formulation. The boundary conditions are, for the horizontal velocity component u(x, h) = U, (2.8) u(x, Y ) = 0, (2.9) and for the vertical velocity component v(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 neglected in the conservation of momentum equation. 28 Chapter 2. Swimming with Lubrication Theory as 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 = 1ω t̂, γ̇ = ω ˆ̇γ, x = 1k x̂, y = hŷ, u = cû, (p, τ ) = µch ( 1 ε p̂, τ̂ ), Û = U c where c = ωk is the speed of the travelling wave. Upon substituting in the above variables (dropping hats), we find that the conservation of mass equation remains unchanged ux + vy = 0. (2.12) Substitution of the non-dimensional variables into the conservation of momen- tum yields (again, hats dropped), px = ετxx,x + τyx,y, py = ε 2τxy,x + ετyy,y. Thus, at leading order we have, px = τyx,y, (2.13) py = 0. (2.14) This tells us that pressure is independent of y and that the xy-component of the stress tensor is the dominant component. The shape of the swimmer becomes Y = a h sin(x− t). (2.15) We can then note that, Yt = −Yx, when in dimensionless variables. Observe that if we integrate the equation px = τyx,y ∫ y Y pxdy = ∫ y Y τyx,ydy, we get that τxy = τ0 + (y − Y )px, (2.16) 29 Chapter 2. Swimming with Lubrication Theory where, τ0 = τxy|y=Y , is the stress on the swimmer. If we observe then that 〈τxy〉 = 〈τ0 + (y − Y )px〉 = 〈τ0 − Y px〉+ y〈px〉 and apply the periodic pressure condition, 〈px〉 = 0 which remains unchanged in dimensionless variables, the zero force condition simplifies to 〈τ0 − pxY 〉 = 0. (2.17) Note that since we have px = τ1 − τ0 1− Y , where τ1 = τxy|y=1, it is the case that τ0 − Y px = τ1 − px. Thus, we have that 〈τ0 − Y px〉 = 〈τ1〉 = 0, (2.18) since we are assuming that 〈px〉 = 0. So, our final equations are then ux + vy = 0, (2.19) τxy = τ0 + (y − Y )px, (2.20) f ( τ̂ , ˆ̇γ ) = 0, (2.21) 〈px〉 = 0, (2.22) 〈τ1〉 = 0, (2.23) and non-dimensional boundary conditions u(x, 1) = U, (2.24) u(x, Y ) = 0, (2.25) 30 Chapter 2. Swimming with Lubrication Theory and for the vertical velocity component v(x, 1) = 0, (2.26) v(x, Y ) = −Yx. (2.27) 2.2 Newtonian Problem Due 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 a non-Newtonian fluid, but the results will be used as an initial guess for the numerical method used to solve the viscoelastic lubrication problem, detailed in the following section. For the viscous problem we have the following: ux + vy = 0, (2.28) px = uyy, (2.29) py = 0, (2.30) 〈uy|y=Y − pxY 〉 = 0, (2.31) 〈px〉 = 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 px is independent of y and apply the boundary conditions on u to obtain u = 1 2 px(y − Y )(y − 1) + U (y − Y ) (1− Y ) . (2.33) Then, look at the continuity equation. Integrate in y first to obtain ∫ 1 Y uxdy + Yx = 0. Pulling the x derivative out of the integral, noting that u|y=Y = 0, and inte- grating in x gives ∫ 1 Y udy + Y = q, where q is the mass flux. Upon calculating this integral and solving for px we 31 Chapter 2. Swimming with Lubrication Theory find px = 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 apply the constraints 〈px〉 = 0, 〈uy|y=Y − pxY 〉 = 0 and solve simultaneously. This involves solving integrals of the form I∗j = ∫ 2pi 0 dx (1− Y )j j = 1, 2, 3 (2.35) This gives us that U = q = 6 ( I∗1 I ∗ 3 − I ∗2 2 ) (4I∗1 I ∗ 3 − 3I ∗2 2 ) . (2.36) Solving the integrals specifically gives I∗1 = 1 (1− ( a h )2 )1/2 , (2.37) I∗2 = 1 (1− ( a h )2 )3/2 , (2.38) I∗3 = ( ( a h )2 + 2) (1− ( a h )2 )5/2 , (2.39) which yields U = q = 3 ( a h )2( 2 ( a h )2 + 1 ) . (2.40) Dimensionally, U c = 3 ( a h )2( 2 ( a h )2 + 1 ) . (2.41) This is exactly the result found by Katz, in [22]. By writing this in the form U c = 3 2 + ( h a )2 , (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 τ1 solves the problem entirely. Of particular interest are the formulas for the pressure gradient and the stress on the wall. The formula for is given in in equation (2.34), and 32 Chapter 2. Swimming with Lubrication Theory the formula for the stress is given by τ1 = (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 it to the viscoelastic problem. 2.3 Viscoelastic Problem With the Newtonian problem solved fully above, we are now in a position to consider the lubrication problem with a viscoelastic fluid in the gap between the swimmer and the wall. When substituting a viscoelastic fluid into the problem for a Newtonian fluid, we have to pick a constitutive model for the viscoelastic fluid. As detailed briefly in Chapter One, this is not necessarily a straight forward task. We use the Oldroyd-8 model because it is the most general weakly- nonlinear constitutive model that includes all possible quadratically nonlinear terms involving the deformation rates allowed by the symmetries of the problem [32]. It is given by, τ + λ1τ ▽ + 1 2 µ0(trτ )γ̇ − 1 2 µ1 (τ · γ̇ + γ̇ · τ ) + 1 2 ν1(τ : γ̇)δ = µ ( γ̇ + λ2γ̇ ▽ − µ2(γ̇ · γ̇) + 1 2 (γ̇ : γ̇)δ ) , (2.44) where δ is the two-dimensional identity tensor. Recall the upper convected derivative of a tensor A, A▽ = ∂A ∂t + u · ∂A− ( (▽u)T ·A+A · ▽u ) . (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: τxx − ( 2De1 + µ1c h − ν1c h ) uyτxy + ( 2De2 + µ2c h − ν2c h ) u2y = 0,(2.46) τxy −De1uyτyy + 1 2 (µ0c h − µ1c h ) (τxx + τyy)uy − uy = 0,(2.47) τyy − (µ1c h − ν1c h ) uyτxy + (µ2c h − ν2c h ) u2y = 0.(2.48) Note that the time derivatives and advective terms are negligible here, just as 33 Chapter 2. Swimming with Lubrication Theory they are in the conservation of momentum equations. The component of the stress tensor that we want to solve for is τxy, since this is the component that appears in the equations of motion. Solving the equations (2.46) through (2.48), we find that τxy = 1 + αu2y 1 + βu2y uy, (2.49) where α = ( ν2c h − µ2c h ) ( De1 − 1 2 ( µ0c h − µ1c h )) + 12 ( µ0c h − µ1c h ) ( 2De2 + µ2c h − ν2c h ) , and β = ( ν1c h − µ1c h ) ( De1 − 1 2 ( µ0c h − µ1c h )) + 12 ( µ0c h − µ1c h ) ( 2De1 + µ1c h − ν1c h ) . We are using this model only as an example. Much more general models can also be dealt with, so long as we know τxy(uy), and that the details of the constitutive model do not interfere with the key simplification of the lubrication theory, namely that the law reduces to that for a steady-state shear flow for the thin-gap geometry with dominant shear stress. 2.3.1 Effective Viscosity Observe that the deviatoric stress equation, τxy = 1 + αu2y 1 + βu2y uy, (2.50) can be written as τxy = µ̃(uy;α, β)uy, (2.51) where we can think of µ̃(uy;α, β) = 1 + αu2y 1 + βu2y as being an effective viscosity. In the Newtonian case we have that µ̃(uy;α, β) = 1. We can see in figure (2.2), detailed in table (2.1), that there are five fluid type regimes, including the Newtonian case, depending on the parameters α and β. Two of these regimes are physically unrealistic because the viscosity diverges and reverses sign for increasing shear rate. We will focus our subsequent analysis on 34 Chapter 2. Swimming with Lubrication Theory the more realistic regimes of Newtonian, shear thinning and thickening fluids.18 0 2 4 6 8 10 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 uy Ef fe ct ive v isc os ity Effective Viscosity versus Shear Rate Figure 2.2: Effective Viscosities. Red: Newtonian, Solid Blue: Shear Thin- ning, Dashed Blue: Unrealistic Shear Thinning, Solid Green: Shear Thickening, Dashed Green: Unrealistic Shear Thickening Fluid Types Newtonian α = β Shear Thinning 0 < α < β Shear Thickening 0 < β < α Unrealistic Thinning α < 0 < β Unrealistic Thickening β < 0 < α Table 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 numerical method detailed in the following section. Obviously, any parameter values that satisfy the restrictions in table (2.1) could be used. These are just chosen as a 18Notice that if we reduce the Oldroyd-8 model to the Oldroyd-B model, we find that α = β = 0, returning to us the Newtonian problem. 35 Chapter 2. Swimming with Lubrication Theory specific example. Fluid Types Parameters α β Newtonian 1 1 Shear Thinning 1 2 Shear Thickening 2 1 Unrealistic Thinning -1 1 Unrealistic Thickening 1 -1 Table 2.2: The specific parameter values chosen for the various fluid types. 2.3.2 Numerical Formulation In this section, we will simplify the viscoelastic lubrication equations in a fashion that allows us to solve them numerically with relative ease. Once the equations are 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 first consider the continuity equation ux+ vy = 0. As we did in the Newtonian case, let us integrate twice to obtain ∫ 1 Y udy + Y = q. (2.52) We will exploit the fact that ∫ 1 Y uydy = U, (2.53) and so can write ∫ 1 Y udy = U − ∫ 1 Y yuydy, which can be found using integration by parts. We have then U + Y − ∫ 1 Y yuydy = q. However, we do not know u explicitly, so let us write it in terms of the deviatoric 36 Chapter 2. Swimming with Lubrication Theory stress. Define T (uy) = τ0 + (y − Y )px = 1 + αu2y 1 + βu2y . (2.54) Solving for y gives y = T − τ0 + Y px px and so we have that dy duy = T ′ px . (2.55) For simplification, we will define γ̇ = uy, (2.56) and γ̇1 = γ̇|y=1, γ̇Y = γ̇|y=Y . (2.57) Substituting this gives ∫ 1 Y yuydy = 1 p2x [∫ γ̇1 γ̇Y TT ′γ̇dγ̇ − (τ0 − Y px) ∫ γ̇1 γ̇Y T ′γ̇dγ̇ ] . (2.58) Let us call I1(γ̇Y , γ̇1) = ∫ γ̇1 γ̇Y TT ′γ̇dγ̇, (2.59) and I0(γ̇Y , γ̇1) = ∫ γ̇1 γ̇Y T ′γ̇dγ̇. (2.60) Note that I0(γ̇Y , γ̇1) = Upx (2.61) Our equations are then as follows: I1(γ̇Y , γ̇1)− p 2 x (U + Y − q − UY )− Uτ0px = 0, (2.62) I0(γ̇Y , γ̇1)− Upx = 0, (2.63) 〈px〉 = 0, (2.64) 〈τ0 − Y px〉 = 0, (2.65) 37 Chapter 2. Swimming with Lubrication Theory where we observe from T (γ̇) = τ0 + (y − Y )px that τ0 = T (γ̇Y ), (2.66) px = T (γ̇1)− T (γ̇Y ) 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 is I1(γ̇ j Y , γ̇ j 1)− ( T (γ̇j1)− T (γ̇ j Y ) 1− Y j )2 ( U + Y j − q − UY j ) − −UT (γ̇jY ) ( T (γ̇j1)− T (γ̇ j Y ) 1− Y j ) = 0, (2.68) I0(γ̇Y , γ̇1)− U ( T (γ̇j1)− T (γ̇ j Y ) 1− Y j ) = 0, (2.69) for j = 1, ..., N , and N∑ j=1 [ T (γ̇j1)− T (γ̇ j Y ) 1− Y j ] ∆ξ = 0, (2.70) N∑ j=1 [ T (γ̇jY )− Y j ( T (γ̇j1)− T (γ̇ j Y ) 1− Y j )] ∆ξ = 0, (2.71) where ∆ξ = (ξj+1 − ξj). (2.72) Thus there are 2N + 2 unknowns in the system U, q, γ̇1..N1 , γ̇ 1..N Y , from which we can easily calculate the pressure gradient and the stress on the swimmer. Figure (2.3) shows the resulting non-dimensional swimming speeds for a fixed propagating wave speed. 38 Chapter 2. Swimming with Lubrication Theory 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Non−Dimensional Swimming Speed a/h U Figure 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 necessarily illuminating. It is not clear which fluid type is optimal for the swimmer, except for a small range of ah . What we would like to do then, is formulate a different method of comparison that allows us to determine the optimal fluid type for the swimmer for the entire range of values of ah . This method is detailed below. 2.3.3 Fixed Dimensional Work To obtain a clear result of which fluid is optimal for the swimmer, let us compare the 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 rate of work per wavelength of the swimmer, W , by calculating the volume integral of 〈τ : γ̇〉. (2.73) Non-dimensionally and at leading order, this gives us W = ∫ 1 Y ∫ 2pi 0 2uyτxydxdy. (2.74) 39 Chapter 2. Swimming with Lubrication Theory This can be rewritten in a way that can utilize the numerics detailed above. That is, we can put W = 1 px ∫ 2pi 0 ∫ γ̇1 γ̇Y 2γ̇T (γ̇)T ′(γ̇)dγ̇dξ. (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 Ũ0,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. By putting things into dimensional form we know that W fixed = 2 µc0,1,2 hk ( 2 px ∫ 2pi 0 ∫ γ̇1 γ̇Y γ̇T (γ̇)T ′(γ̇)dγ̇dξ )0,1,2 . (2.77) We can easily solve for the speed of the propagating waves and therefore the dimensional swimming speed. From figure (2.4) we can see that the shear thinning fluid yields the fastest dimensional swimming speed for all values of ah . From figure (2.4) we can also see that, contrary to what figure (2.3) might suggest, it is not optimal to swim near the wall. 40 Chapter 2. Swimming with Lubrication Theory 0 0.2 0.4 0.6 0.8 1 0 1 2 x 10−4 Dimensional Swimming Speed a/h U* c Figure 2.4: Dimensional Swimming Speeds. Red: Newtonian. Blue: Shearing Thinning. Green: Shear Thickening. 2.3.4 Limiting Cases In the above analysis, there are two cases that may be of analytical interest when thinking of Y = ah sin ξ, that when a/h −→ 1, or the swimmer approaches the wall and that when a/h −→ 0. To calculate these limits we assume that γ̇ −→∞ and that γ̇ −→ 0 respectively. Let us consider the former limit first. With this assumption of the leading order component of the strain rate approaching infinity, we have that τxy = T (γ̇) ∼ α β γ̇. (2.78) This simplifies the integrals I0 and I1 significantly I0 = β 2α ( τ21 − τ 2 0 ) , (2.79) I1 = β 3α ( τ31 − τ 3 0 ) . (2.80) We can then solve the system of equations (3.23) through (3.26) explicitly for 41 Chapter 2. Swimming with Lubrication Theory U . We find that U = 6 ( I∗1 I ∗ 3 − I ∗2 2 ) 4I∗1 I ∗ 3 − 3I ∗2 2 , (2.81) where I∗j are defined above in equation (2.35). We can see that this swimming speed is independent of the parameters of the problem, and in fact agrees with the Newtonian solution. Since we are assuming that the swimmer is approaching the wall here, let us put Y = (1− ε2) sin(ξ). (2.82) If we assume that the dominant contribution to the integrals I∗j are at the point along the swimmer that is nearest the wall, [8], that is at ξ = pi/2, then we have U −→ 1 as a/h −→ 1. (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 behaviour near one in the figure can be attributed to a numerical singularity due to the presence of 1/(1− Y ). Solving for the remainder of the variables gives that q = (2U − 3)I∗1 3I∗2 + 1, (2.84) px = α β 6(U − 2) (1− Y )2 − α β 12(q − 1) (1− Y )3 , (2.85) τ1 = α β 4U − 6) (1− Y ) − α β 6(q − 1) (1− Y )2 . (2.86) Notice that U and q are independent of the parameters of the problem α and β, but px and τ1 are not. Assuming the shape of the swimmer to be that as stated in equation (2.82) gives us to leading order q = 1, (2.87) px = α β −6 (1− sin(x− t))2 , (2.88) τ1 = α β −2 (1− sin(x− t)) + α β 12 (1− sin(x− t))2 . (2.89) Similarly, we have when a/h −→ 0, assuming this implies γ̇ −→ 0, that τxy = T (γ̇) ∼ γ̇. (2.90) 42 Chapter 2. Swimming with Lubrication Theory Thus the problem simplifies to the Newtonian low amplitude problem. Here, the integrals I0 and I1 simplify to I0 = 1 2 ( τ21 − τ 2 0 ) , (2.91) I1 = 1 3 ( τ31 − τ 3 0 ) . (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 swimmer is very small. Let us put Y = ε sin ξ. (2.93) We find that to leading order U = 3ε2. (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 swimming speed appears to be approaching zero quadratically. Solving for the remainder of the variables to leading order gives us that q = U = 3ε2, (2.95) px = 12Y, (2.96) τ1 = 6Y. (2.97) Thus, we can see that these two different limiting cases simplify into problems that are solvable analytically. In the latter limit, we found that the problem reduces to the Newtonian problem. In the former limit, we found that the problem reduced to one similar to the Newtonian, but with a simple dependence on the viscoelastic parameters, α and β. 2.4 Bingham Problem As stated above, much more general models can be dealt with in this formulation of the lubrication problem, so long as we know τxy(uy), and that the details of the constitutive model do not interfere with the key simplification of the lubrication theory, namely that the law reduces to that for a steady-state shear flow for the thin-gap geometry with dominant shear stress. Let us consider the 43 Chapter 2. Swimming with Lubrication Theory specific 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 rheological model, leaving only the effect of shear thinning or thickening. Considering a Bingham fluid leads to a more interesting lubrication problem from the rheo- logical perspective. As well, the addition of yield stress is biologically relevant to moving organisms [8, 26]. A Bingham fluid has the following constitutive law τxy = τysgn(uy) + µuy, |τxy| > τy. (2.98) If |τxy| ≤ τy then it is the case that uy = 0 and the fluid does not yield. As stated above, the standard set of non-dimensional equations governing this lubrication problem are: ux + vy = 0, (2.99) τxy = τ1 − (1− y)px, (2.100) 〈px〉 = 0, (2.101) 〈τ0 − pxY 〉 = 〈τ1〉 = 0, (2.102) and the shape of the swimmer is Y = a h 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 according to the lubrication approximation we obtain, τxy = Bsgn(uy) + uy |τxy| > B, (2.104) where B = hµcτy is the Bingham number. If |τxy| ≤ B then it is the case that uy = 0 and the fluid does not yield. Note that this implies that τxy,y = uyy, (2.105) which tells us, recalling that px = τxy,y, px = uyy. (2.106) 44 Chapter 2. Swimming with Lubrication Theory Let us begin solving this by equating our two expressions for τxy. This yields, τ1 − (1− y)px = Bsgn(uy) + uy. (2.107) Rearranging to solve for uy gives, uy = τ1 −Bsgn(uy)− (1− y)px. (2.108) If |τxy| ≤ B then it is the case that uy = 0 and the fluid does not yield. As is evident by equation (2.108), the flow profile is going to be dependent on the sign of uy. The different values of sgn(uy) represent different regions of the flow profile. Let us classify these regions as illustrated in table (2.3). We will define y = Y+ as the upper yield surface, and y = Y− as the lower yield surface. It is the case that the fluid yields if y > Y+ or y < Y−. If Y− < y < Y+, it is the case that 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 Region Location y < Y− or y > Y+ Y− < y < Y+ Stress |τxy| > B |τxy(Y±)| = B Strain Rate |uy| > 0 uy = 0 Flow Type parabolic plug-like Table 2.3: Bingham Fluid Flow Profile Regions. In the outer regions, we see parabolic flow profiles typically characterized by Newtonian flows. In the inner region, we see what we refer to as pseudo-plug flow. In this region the velocity profile is constant. Basically, what we expect is a flow profile that looks similar the Newtonian flow profile, see figure (2.5), but with a flattened region around where uy = 0. 45 Chapter 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.5 0 0.5 1 u y Newtonian Flow Profile for a = 0.5 Figure 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 that u|(y=Y−) = up = u|(y=Y+) (2.110) where, up = Const, is the speed of the plug-like flow in the inner region. We call this region the pseudo-plug. Define σ = sgn(uy|y>Y+). (2.111) Following immediately from this is, −σ = sgn(uy|y Y−, tells us that sgn(px) = σ. (2.118) Integrating the continuity equation across the gap, applying the boundary conditions on v, and then integrating in x gives q = Y + ∫ 1 Y udy, (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 five possible cases. These are detailed in table (2.4).19 Case Location Pressure Pseudo-Plug Speed I Y− < Y+ < Y < 1 px > 0 up < 0 II Y− < Y < Y+ < 1 px > 0 up = 0 III Y < Y− < Y+ < 1 No restrictions No restrictions IV Y < Y− < 1 < Y+ px < 0 up = U V Y− < Y+ < Y < 1 px < 0 up > 0 Table 2.4: Location of Yield Surfaces. We consider each one of these cases separately, and obtain different equations 19For cases I and V, it is shown belwo that sgn(up) = −sgn(px). Due to the fact that the flow 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. 47 Chapter 2. Swimming with Lubrication Theory for each. These equations are summarized below.20 For case I, Y+ = 1 2 (1 + Y )− U px(1− Y ) , (2.120) Y− = Y+ − 2B |px| , (2.121) q = Y + U 2 (1− Y )− px 12 (1− Y )3. (2.122) For case II, Y+ = 1− √ 2U px , (2.123) Y− = Y+ − 2B |px| , (2.124) q = Y + px 6 (1− Y+) 3. (2.125) For case III, Y+ = ( 1− Y 2 − 2Upx − 4 B |px| ( Y + B|px| )) 2(1− Y − 2 B|px| ) , (2.126) Y− = Y+ − 2B |px| , (2.127) q = Y + px 6 ( (1− Y+) 3 + (Y− − Y ) 3 − 3(1− Y )(Y− − Y ) 2 ) .(2.128) For case IV, Y− = Y + √ −2U px , (2.129) Y+ = Y− + 2B |px| (2.130) q = Y + px 6 (Y− − Y ) 3 + U(1− Y ), (2.131) 20For the full analysis of each case, the reader is referred to Appendix B. 48 Chapter 2. Swimming with Lubrication Theory and for case V, Y− = 1 2 (1 + Y )− U px(1− Y ) , (2.132) Y+ = Y− + 2B |px| , (2.133) q = Y + U 2 (1− Y )− px 12 (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, 〈τ1〉 = 0, (2.135) 〈px〉 = 0. (2.136) We can now solve the system using Newton’s Iteration. The trick is, at each step in the iteration, checking which case is relevant, and then using the appropriate equations. This is done by checking the locations of Y+ and Y−, as well as the sign of px (for cases I, II, IV and V). Please refer to the appendices for the full code used to solve this problem. Figure (2.6) below shows the pressure gradient and the yield surfaces found for ah = 1/2 and B = 1. 49 Chapter 2. Swimming with Lubrication Theory −1.5 −1 −0.5 0 0.5 1 1.5 −4 −2 0 2 4 6 8 10 12 14 16 x p x (a) Pressure Gradients −1.5 −1 −0.5 0 0.5 1 1.5 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 x Y + an d Y − (b) Pseudo-Plug Figure 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. The surface of the swimmer is represented by the black curve and the pink vertical line shows where the pressure gradient changes sign. Figure (2.7) shows the velocity profiles characteristic to each case discussed above. 50 Chapter 2. Swimming with Lubrication Theory 0 0.1 0.2 0.3 0.4 0.5 0.4 0.5 0.6 0.7 0.8 0.9 1 u y (a) Case I −0.1 0 0.1 0.2 0.3 0.4 0.5 0.5 0.6 0.7 0.8 0.9 1 u y (b) Case II 0 0.1 0.2 0.3 0.4 0.5 0.6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 u y (c) Case III 0 0.1 0.2 0.3 0.4 0.5 0.6 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 u y (d) Case IV 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 u y (e) Case V Figure 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µcτy, is increased, the size of the pseudo-plug region increases until, eventually, it fills up most of the flow domain with viscous boundary layers next to the walls and around the location where the pressure gradient changes sign. As the Bingham number approaches zero, the size of the pseudo-plug region approaches zero and the solution approaches the Newtonian solution, where Y+ = Y−, see figure (2.8). 51 Chapter 2. Swimming with Lubrication Theory −1.5 −1 −0.5 0 0.5 1 1.5 −10 −5 0 5 10 15 20 25 30 35 40 x p x (a) Pressure Gradients for B = 5 −1.5 −1 −0.5 0 0.5 1 1.5 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 x Y + an d Y − (b) Pseudo-Plug for B = 5 −1.5 −1 −0.5 0 0.5 1 1.5 −2 0 2 4 6 8 10 12 x p x (c) Pressure Gradients for B = 0.1 −1.5 −1 −0.5 0 0.5 1 1.5 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 x Y + an d Y − (d) Pseudo-Plug for B = 0.1 Figure 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 Bingham red. In the right figures, Y+ is blue, Y− is red, and the shaded green represents the pseudo-plugs. The surface of the swimmer is represented by the black curve and the pink vertical line shows where the pressure gradient changes sign. Figure (2.9) shows the how the swimming speed and mass flux vary with the Bingham number. The mass flux is a monotonically decreasing function of the Bingham number, which is expected due to the fact that the yielded region decreases in size as the yield stress increases. Whereas, the swimming speed is a non-monotonic function of the Bingham number with a clear maximum value, located at B ≃ 1. This suggests that is it beneficial for an organism to swim in a fluid that has a moderate yield stress.21 21It has been observed that snails and slugs move over a fluid that exhibits a yield stress [8]. 52 Chapter 2. Swimming with Lubrication Theory 0 1 2 3 4 5 0.485 0.49 0.495 0.5 0.505 0.51 0.515 B U (a) U versus B for a h = 1/2 0 1 2 3 4 5 0.46 0.465 0.47 0.475 0.48 0.485 0.49 0.495 0.5 0.505 B q (b) q versus B for a h = 1/2 Figure 2.9: Swimming speed and mass flux for varying Bingham number. The corresponding Newtonian values are shown in red. 2.5 Discussion In this chapter we presented the problem of a lubricated swimmer by a wall. We formulated the problem initially with a Newtonian fluid, and returned the results of previous work done by others [22]. Next, we considered the problem with 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 for the swimmer given a fixed rate of work when compared to a swimmer in a Newtonian fluid and a shear thickening fluid. It is the case that viscoelasticity can 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 limit as 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 In the opposite limiting case, when the ratio of the amplitude of the swimmer on the gap width goes to zero, we have the Newtonian problem and its relevant 22This contradicts the limiting case of the introductory problem, where the limit was found to depend on the parameters of the problem. This is not of concern due to the fact that, in the introductory problem, no thought was given to what the strain rates were doing in the limit. The introductory problem was just a preliminary example used to motivate the use of lubrication theory. In fact, if we reduce the Oldroyd-8 constitutive model to the Oldroyd-B model above, the Newtonian problem is returned, hence no contradiction when the limit is taken properly. 53 Chapter 2. Swimming with Lubrication Theory results returned. That is, we found that the non-dimensional swimming speed approaches zero quadratically in a. Due to the formulation of this problem, theoretically it should be easy to consider different constitutive models and wave shapes, thus allowing the range of applications of this problem to broaden. For example, perhaps we would like to consider a different swimmer shape in Chapter Two. Let us consider the toy example of a swimmer shape given by Y = tanh(12 sin ξ), as shown in figure (2.10a). The analysis of the problem does not change. The resulting swimming speeds of such a swimmer are shown in figure (2.10b) and were found using the same code as was used for the results found above. It is evident that different swimmer shapes can yield different results, depending on the application one has in mind. The mathematical set up we chose here lends itself nicely to generality, and exploring these potential different 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.2 0 0.2 0.4 0.6 0.8 1 ξ Y Square Wave Swimmer (a) Wave Shape 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.5 0 0.5 1 1.5 2 2.5 3 x 10−4 Dimensional Swimming Speed a/h U* c (b) Swimming Speeds Figure 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 the problem, we had to consider five different possible cases that depended on the location of the yield surfaces. For a fixed swimmer amplitude, ah = 1/2, we solved for the pressure gradient and the location of the yield surfaces. At this fixed amplitude, we found that it is potentially advantageous for an organism to swim in a Bingham fluid with a moderate yield stress, as it results in a faster swimming speed than a Newtonian fluid. 54 Chapter 3 Elastic Swimming with Lubrication Theory In the previous chapter we looked at the swimming problem from the perspective of the swimmer having a known motion or shape. We then analyzed the resulting swimming speed with the swimmer subjected to a fixed dimensional work or rate of work per unit wavelength. There is, however, a different approach one could use to look at this problem. Rather than a prescribed shape, we can consider a prescribed force distribution along the length of the swimmer. In this situation we have to think more carefully about the swimmer itself. Previously, we only were really concerned with the swimming speed. Now, we must think about how this applied force affects the swimmer. What kind of properties does this swimmer have to govern how it will respond to this forcing? This is a question that has a large number of potential answers. The specifics of flagellum movement are very complex, and so if we were to model the swimmer as a realistic flagella the problem would be very involved [23, 22]. However, here we 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 the swimmer 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 realistic flagellum. In fact, this assumption makes our problem more general with a wider 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 the entire length of the swimmer, we could consider the situation of a passive elastic swimmer with one fixed end that moves in a vertical fashion thus propagating waves down the length of the swimmer. This postulation is interesting from a mathematical perspective, and perhaps has some interesting alternate applica- tions, but it has been shown that observed wave forms of spermatozoa cannot 55 Chapter 3. Elastic Swimming with Lubrication Theory be accounted for with this mechanism [29]. Formally then, this is the problem of an elastic swimmer subject to a known forcing along its length in a lubrication setting. Its shape and swimming speed are unknown. Our goal in this chapter is to discover how an elastic swimmer reacts when subjected to a known distribution of applied torques whilst in a viscoelastic fluid. However, the problem of how such a swimmer acts in a Newtonian fluid proves to be interesting itself. 3.1 Mathematical Formulation Let us consider the following: a thin, isotropic elastic beam in two-dimensions immersed in a fluid very close to a rigid wall. This beam is then subjected to known internal forces distributed along its length, f(x, t) = a sin k(x− ct) that cause 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 of simplicity, we will assume that the applied force is sinusoidal. We use a similar set up for our problem as was used in [2], where the forces applied are referred to as actively distributed torques. Figure 3.1: Geometry of Lubrication Problem. This time, the swimmer shape in unknown as well as the swimming speed What we want to discover here is what the unknown swimmer shape is and the swimming speed it produces. To begin with, let us consider the boundary conditions. In the frame moving with the swimmer, our boundary conditions remain unchanged for this problem from the last. That is, non-dimensionally they are for the horizontal velocity 56 Chapter 3. Elastic Swimming with Lubrication Theory component, u(x, h) = U, (3.1) u(x, Y ) = 0, (3.2) and for the vertical velocity component v we have v(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 in the original frame with the wall fixed, the boundary condition on the swimmer would actually be a kinematic or free boundary condition due to the fact that the swimmer is now considered elastic. But, in the frame moving with the swimmer there is no horizontal velocity associated with the beam, therefore the kinematic 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, by ▽ p = ▽ · τ , (3.5) ▽ · u = 0, (3.6) and general constitutive model f(τ , γ̇) = 0. (3.7) Here we have an elastic beam, thus we must include another equation for the vertical force balance. This is the beam equation including an internally applied torque [2, 28]. In full dimensional terms this equation is, neglecting gravity, 0 = −∂xx (f +D∂xxY ) + p, (3.8) where D is the bending stiffness of the beam, and p, Y , f the pressure, beam shape and applied forcing respectively. We enforce periodicity in the swimmer up to and including its third derivative. Mathematically it is required for the higher derivatives due to the fact that, once the beam equation is incorporated, we end up with a fourth order ordinary differential equation, thus requiring the four conditions. Physically, only this condition and the continuity of the 57 Chapter 3. Elastic Swimming with Lubrication Theory swimmer are required.23 In truth, we are interested in the pressure gradient only, thus we take the derivative of the above equation and use px = ∂xxx (f +D∂xxY ) . (3.9) We also have the constraints 〈px〉 = 0, (3.10) 〈τ1〉 = 0, (3.11) 〈Y 〉 = 0. (3.12) When looking at equation (3.9), we can see that in order for 〈px〉 = 0 to hold we require that 〈Yxxxxx〉 = 0. This is the same thing as saying Y is periodic in its fourth derivative. This requirement is due to the choice of f(x, t). Thus we can enforce the constraint on the pressure gradient by enforcing periodicity on the swimmer shape up to and including its fourth derivative, and from here on out not include the constraint specifically. Observe that the condition 〈px〉 = 0 is also utilized in the formulation of the zero force on swimmer condition when it is put into the form 〈τ1〉 = 0, as noted in Chapter One. Non-dimensionalizing as we did in Chapter Two, using the following non- dimensional variables, t = 1ω t̂, γ̇ = ω ˆ̇γ, x = 1k x̂, y = hŷ, u = cû, D̂ = εh 2k4 µc D (p, τ ) = µch ( 1 ε p̂, τ̂ ), â = εk2h µc a applying the lubrication approximation and simplifying algebraically yields, dropping hats, I1(γ̇Y , γ̇1)− p 2 x (U + Y − q − UY )− Uτ0px = 0, (3.13) I0(γ̇Y , γ̇1)− Upx = 0, (3.14) 〈Y 〉 = 0, (3.15) 〈τ1〉 = 0. (3.16) 23This condition’s specific use was not required previously because the shape of the swimmer was assumed. 58 Chapter 3. Elastic Swimming with Lubrication Theory Recall that I0 = ∫ γ̇1 γ̇Y T ′γ̇dγ̇, (3.17) I1 = TT ′γ̇dγ̇, (3.18) T = τxy. (3.19) Similarly then, we also utilize the fact that px = T (γ̇1)− T (γ̇Y ) 1− Y , (3.20) to remove px from the above equations. As well, recall that τ1 = T (γ̇1). (3.21) Putting the beam equation into non-dimensional terms yields px = ∂xxx (f +D∂xxY ) . (3.22) So, our unknowns are now γ̇1, γ̇Y , U , and q and the new unknown, the swimmer shape, Y . Summarizing the full set of non-dimensional equations relevant to solving the problem gives us I1(γ̇Y , γ̇1)− p 2 x (U + Y − q − UY )− Uτ0px = 0, (3.23) I0(γ̇Y , γ̇1)− Upx = 0, (3.24) f ( τ̂ , ˆ̇γ ) = 0, (3.25) px − ∂xxx (f +D∂xxY ) = 0, (3.26) 〈Y 〉 = 0, (3.27) 〈τ1〉 = 0, (3.28) along with periodic conditions on the derivatives of Y up to and including the fourth derivative which can be written as 〈Yx〉 = 〈Yxx〉 = 〈Yxxx〉 = 〈Yxxxx〉 = 〈Yxxxxx〉 = 0 (3.29) These equations are very similar to what we had before with only the addition of the beam equation and the periodicity conditions on the swimmer and its first five derivatives. Again, the conditions on the swimmer were true before 59 Chapter 3. Elastic Swimming with Lubrication Theory since we were imposing the shape of the swimmer, but now we must rigorously enforce it. As we have seen above, the layout of this problem and the last one has the fluid properties included simply in the form of the equation τxy = T (γ̇). This is an elegant and general way of looking at the problem and allows us to look at different constitutive equations with relative ease. However, the complexity of the constitutive equation can still heavily affect the ease of solving the system both analytically and numerically. As done above, let us first consider the Newtonian problem in order to gain some insight into the problem before considering more complex fluid models. 3.2 Newtonian Problem Since, 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 they are for the Newtonian problem for a swimmer of prescribed shape, we can skip a few steps. We have already solved equations (3.23), (3.24), and (3.28) for a Newtonian fluid in the previous chapter. However, since Y is unknown, we can only utilize the equations for the pressure gradient and stress on the wall, as the calculation of the swimming speed and mass flux relied on calculating the integrals I∗j . Not all is lost though, we can use still use the following equations for px and τ1 px = 6(U − 2) (1− Y )2 − 12(q − 1) (1− Y )3 , (3.30) τ1 = (4U − 6) (1− Y ) − 6(q − 1) (1− Y )2 . (3.31) These are most useful since the pressure gradient and unknown swimmer shape appear 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 us ∂xxx (a sin(x− t) +D∂xxY ) = 6(U − 2) (1− Y )2 − 12(q − 1) (1− Y )3 . (3.32) If we introduce a new variable, ζ = (1− Y ), (3.33) 60 Chapter 3. Elastic Swimming with Lubrication Theory and utilize the dependence on ξ = (x−t) in the problem, we obtain the equation ζ(5) = − 1 D ( a cos ξ + 6(U − 2) ζ2 − 12(q − 1) ζ3 ) , (3.34) where derivatives are with respect to ξ. We can then solve for ζ and U and q by applying the constraints. To do this, we must first put the constraints into terms of ζ. The first constraint, < Y >= 0 is simple. It becomes 〈ζ〉 = 1. (3.35) The second, 〈τ1〉 = 0 is less straight forward. However, due to our above algebraic manipulation, we already know what τ1 is in terms of Y , U , and q. Thus it is a simple matter of putting this in terms of ζ τ1 = (4U − 6) ζ − 6(q − 1) ζ2 , (3.36) and so the constraint becomes (4U − 6)〈ζ−1〉 − 6(q − 1)〈ζ−2〉 = 0. (3.37) To solve this analytically requires solving a cubic, so rather than doing that, let us convert the problem to something which we can apply numerics. 3.2.1 Numerical Formulation We 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 solve this numerically, let us define a new variable r such that r′ = ξ, thus giving us r = ∫ ξ 0 ζdξ. (3.38) Similarly define the variable s so that s′ = ( (4U − 6)ξ−1 − 6(q − 1)ξ−2 ) , which gives us s = ∫ ξ 0 ( (4U − 6)ξ−1 − 6(q − 1)ξ−2 ) dξ. (3.39) We have then that r(0) = 0, r(2pi) = 2pi, (3.40) 61 Chapter 3. Elastic Swimming with Lubrication Theory and s(0) = 0, s(2pi) = 0, (3.41) as conditions on these new variables. Thus, our system of equations to solve is ζ(5) = − 1 D ( a cos ξ + 6(U − 2) ζ2 − 12(q − 1) ζ3 ) , (3.42) r = ∫ ξ 0 ζdξ, (3.43) s = ∫ ξ 0 ( (4U − 6)ξ−1 − 6(q − 1)ξ−2 ) dξ, (3.44) where we think of U and q as eigenvalues, with conditions ζ(0) = ζ(2pi), ζ ′(0) = ζ ′(2pi), ζ ′′(0) = ζ ′′(2pi), ζ(3)(0) = ζ(3)(2pi), ζ(4)(0) = ζ(4)(2pi), (3.45) and r(0) = 0, r(2pi) = 2pi, (3.46) s(0) = 0, s(2pi) = 0, (3.47) for unknowns ζ, ζ ′, ζ ′′, ζ(3), ζ(4), ζ(5), r and s with eigenvalues U and q. To do this we use the MATLAB solver bvp4c. To see the code please refer to appendices. Figure (3.2) shows the resulting shape of the swimmer as the amplitude of the forcing increases. It shows us that as the amplitude gets larger the swimmer shape develops a plateau region followed by a sharp spike region. By looking at the phase portraits of this large amplitude solution shown in figure (3.3) we can see that there exists a inner region that reminiscent of a fixed point. This explains the section of the swimmer shape that is nearly constant that exists near ξ ∼ pi/2. 62 Chapter 3. Elastic Swimming with Lubrication Theory 0 1 2 3 4 5 6 −2.5 −2 −1.5 −1 −0.5 0 0.5 Swimmer Shape for D=0.2 and a=0..25 ξ=x−t Y 5 10 15 20 25 0.2 0.4 0.6 0.8 a (a) Newtonian Solution for D = 0.2 0 1 2 3 4 5 6 −2.5 −2 −1.5 −1 −0.5 0 0.5 Swimmer Shape for D=1.5 and a=100 ξ=x−t Y 20 40 60 80 100 0.2 0.4 0.6 0.8 a (b) Newtonian Solution for D = 1.5 Figure 3.2: Newtonian Solution of Elastic Lubrication Problem. In the bottom 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 values for each case. −3 −2 −1 0 1 −4 −2 0 2 4 −10 0 10 Y ‘ Phase Space for D = 0.2 Y Y ‘‘ (a) Phase Portrait for D=0.2 −3 −2 −1 0 1 −4 −2 0 2 4 −10 0 10 Y ‘ Phase Space for D = 1.5 Y Y‘ ‘ (b) Phase Portrait for D=1.5 Figure 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, it is not immediately obvious how we may scale things in such a way that gives us a rigorous asymptotic solution to this inner region. Let us assume that we are inside this inner region for ζ, near ξ = pi/2, and expand in the standard asymptotic fashion ζ = ζp + ζ1 + . . . (3.48) where ζp is the plateau constant and |ζ1| << 1. Substituting this into equation 63 Chapter 3. Elastic Swimming with Lubrication Theory (3.34), and rearranging slightly we obtain D a d5 dξ5 ζ1 = − cos (pi 2 + σ ) − 6(U − 2) a(ζp + ζ1) + 12(q − 1) a(ζp + ζ1) , (3.49) where we are assuming that σ is some small number. It is easy to see then that at leading order we obtain ζp = 2(q − 1) U − 2 . (3.50) To solve this fully analytically is certainly a challenge, but if we take this formula for the plateau value, recalling that ζ = 1− Y , and plot it using our numerical values for U and q we obtain a very remarkable fit, as can be seen as the magenta line in figure (3.2). Note that the plateau value is independent of the bending stiffness. These results suggest that there are some truly interesting dynamics at work in this large amplitude regime which most certainly should be analyzed at greater depth, including a full analysis of the outer region. We can also see from figure (3.2b) that the plateau approaches y = 1. This observation is substantiated by the formula Yp = 1− 2(q − 1) U − 2 , (3.51) if we observe that q → 1 as the amplitude gets large. As well, we can note from figure (3.2) that there is an obvious optimal forcing amplitude for the swimming speed and after that, as the swimmer gets too close to the wall and the stresses get large, the speed decreases. 3.2.2 Low Amplitude Solution Let us consider the situation when the amplitude of the applied forcing or torque on the swimmer is very small. That is, put a≪ 1, (3.52) where the forcing is given by f = a sin(x − t). We can then assume that this implies that the amplitude of the swimmer shape is very small |Y | ≪ 1. (3.53) 64 Chapter 3. Elastic Swimming with Lubrication Theory At 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 approximations for px and τ1 as we found in Chapter Two. We found initially that px = 6(U − 2) (1− Y )2 − 12(q − 1) (1− Y )3 , (3.54) τ1 = (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 able to get a low amplitude approximation for these equations. We obtain the same leading order approximations as we did before px = 12Y, (3.56) τ1 = 6Y. (3.57) Notice that the constraint on τ1 holds automatically if we enforce the constraint on 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 constraint 〈Y 〉 = 0, (3.59) and periodicity in Y up to and including its fourth derivative. Solving this we find that the homogeneous solutions must be ignored due to the fact that they do not satisfy the periodicity conditions. Thus it is only the particular solution that is relevant. We find the particular solution to be Y = a D2 + 122 (−12 cos ξ +D sin ξ) . (3.60) From this low amplitude solution we can see that for a moderate value of D all terms will balance, but for large or small values of D we will have different terms balancing. For small amplitude forcing and a small bending stiffness we 24Note that the low amplitude approximation does not change the beam equation since all terms in it become small and therefore all balance 65 Chapter 3. Elastic Swimming with Lubrication Theory would expect the shape of the swimmer to look like Y = −12a D2 + 122 cos ξ (3.61) and for small amplitude forcing and large bending stiffness we expect the shape of the swimmer to look like Y = aD D2 + 122 sin ξ (3.62) 3.3 Viscoelastic Results Let us now tackle the elastic lubrication problem with a viscoelastic fluid within the gap. Here again, we choose to use the Oldroyd-8 constitutive model. Non- dimensionalizing the constitutive model, simplifying, and solving for τxy, as we did in the previous chapter, yields τxy = T (γ̇) = (1 + αγ̇2) (1 + βγ̇2) γ̇. (3.63) To solve this numerically, we employ a similar method as was used in Chapter Two. That is, we discretize the system of equations and apply Newton’s Method. The system of discretized equations is as follows: I0(γ̇ j Y , γ̇ j 1)− Up j x = 0, (3.64) I1(γ̇ j Y , γ̇ j 1)− p j2 x ( U + Y j − q − UY j ) − UT (γ̇jY )p j x = 0, (3.65) pjx − f j xxx −DY j xxxxx = 0, (3.66) ΣNj=1T (γ̇ j 1) = 0, (3.67) ΣNj=1Y j = 0. (3.68) There are 3N + 2 unknowns, γ̇1..NY , γ̇ 1..N 1 , Y 1..N , U , and q. And we have that, pjx = T (γ̇j 1 )−T (γ̇j Y ) 1−Y j . Note that fxxx is given, as f is given. The Yxxxxx term is found using a difference formula. Please refer to appendices for the full numerical method used to solve this problem. Figure (3.4) shows the swimmer shape varying for changing ah values and the resulting swimming speeds. We see that for low bending stiffness, the differences between 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 shear 66 Chapter 3. Elastic Swimming with Lubrication Theory thickening fluid yields the slowest. But, as we increase the bending stiffness, as seen in figures (3.5) and (3.6), both the swimmer’s shape and speed becomes less differentiated between the different fluid types. Notice the asymmetry here for certain values of the bending stiffness and the similarities between the shape of the swimmer for Newtonian, shear thinning, shear thickening. 2 4 6 −0.1 −0.05 0 0.05 α=β=1 and D=0.2 x Y 0.2 0.4 0.6 0.8 1 1.2 0.01 0.02 0.03 0.04 a U U vs a for D=0.2 2 4 6 −0.1 −0.05 0 0.05 0.1 α=1, β=2, D=0.2 x Y 2 4 6 −0.1 −0.05 0 0.05 α=2, β=1, D=0.2 x Y Figure 3.4: Solution for viscoelastic elastic lubrication problem for various fluid types with D=0.2. Red: Newtonian. Blue: Shearing Thinning. Green: Shear Thickening. 67 Chapter 3. Elastic Swimming with Lubrication Theory 2 4 6 −0.05 0 0.05 α=β=1 and D=7 x Y 0.2 0.4 0.6 0.8 1 1.2 0.005 0.01 0.015 0.02 0.025 a U U vs a for D=7 2 4 6 −0.05 0 0.05 α=1, β=2, D=7 x Y 2 4 6 −0.05 0 0.05 α=2, β=1, D=7 x Y Figure 3.5: Solution for viscoelastic elastic lubrication problem for various fluid types with D=7. Red: Newtonian. Blue: Shearing Thinning. Green: Shear Thickening. 68 Chapter 3. Elastic Swimming with Lubrication Theory 2 4 6 −0.04 −0.02 0 0.02 0.04 α=β=1 and D=27 x Y 0.2 0.4 0.6 0.8 1 1.2 1 2 3 4 5 x 10−3 a U U vs a for D=27 2 4 6 −0.04 −0.02 0 0.02 0.04 α=1, β=2, D=27 x Y 2 4 6 −0.04 −0.02 0 0.02 0.04 α=2, β=1, D=27 x Y Figure 3.6: Solution for viscoelastic elastic lubrication problem for various fluid types with D=27. Red: Newtonian. Blue: Shearing Thinning. Green: Shear Thickening. 3.3.1 Low Amplitude Solution Let us again consider the situation when the amplitude of the applied forcing or torque on the swimmer is very small. Above we assumed that the small forcing amplitude yields a small swimmer amplitude. Here, we will assume the same, but also assume that this very small forcing and resulting small swimmer amplitude implies a very small shear rate |γ̇| ≪ 1. (3.69) This small shear rate affects the form of the constitutive equation for the stress. In particular, we have that the constitutive equation becomes T (γ̇) = γ̇, (3.70) reducing the problem to the Newtonian one. Therefore, the solution found above for the low amplitude Newtonian problem is also the solution for the low 69 Chapter 3. Elastic Swimming with Lubrication Theory amplitude viscoelastic problem. Recall that we found Y = a D2 + 122 (−12 cos ξ +D sin ξ) . (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 bending stiffness it will look like a sine curve. As we can see from figure (3.7) these analytical results agree very nicely with the numerical viscoelastic results for small amplitude forcing. 1 2 3 4 5 6 −5 0 5 x 10−4 α=β=1 and D=1 x Y 1 2 3 4 5 6 −5 0 5 x 10−4 α=1, β=2, D=1 x Y 1 2 3 4 5 6 −5 0 5 x 10−4 α=2, β=1, D=1 x Y Figure 3.7: Comparison between analytical and numerical results for low am- plitude swimming problem. Numerical results are solid lines, analytical results are shown by *. Red: Newtonian. Blue: Shearing Thinning. Green: Shear Thickening. 3.4 Discussion In this chapter we presented the problem of a lubricated elastic swimmer by a wall. We considered first the Newtonian problem and discovered that it yields very interesting results. For large amplitude forcing we found that a plateau region develops along the swimmer. This plateau approaches the wall as the 70 Chapter 3. Elastic Swimming with Lubrication Theory amplitude of the forcing becomes large. We referred to this region as the inner region as it corresponds to a region in the phase portrait that is reminiscent of a fixed point. Upon expanding around this inner region, we found a value of this plateau analytically that agreed very well with numerical results. Follow- ing this plateau region is a large spike type region. Unfortunately, solving for this region analytically is beyond the scope of this thesis. These results are not entirely 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, in the not so distant future, they will be analyzed in greater depth. As well, we saw that there is a definite forcing amplitude that yields an optimal swimming speed. The shape of the optimal swimmer at this amplitude can be seen in figure (3.8). 0 1 2 3 4 5 6 −1.4 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 Swimmer Shape for D=0.2 and a=0..6 ξ=x−t Y Figure 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 amplitude limit. We next calculated the low amplitude analytical solution of the Newtonian 71 Chapter 3. Elastic Swimming with Lubrication Theory problem. We found that both the pressure gradient and the stress on the wall are proportional to the swimmer shape, and that the swimmer shape is Y = a D2 + 122 (−12 cos ξ +D sin ξ) , (3.72) a solution that is valid for any value of the bending stiffness D. This result was returned to us when we calculated the low amplitude solution for a elastic swimmer in a viscoelastic fluid. It was found to agree with the numerical results for the viscoelastic problem extremely well. We then took the next logical step and considered a viscoelastic fluid. We found that for moderate bending stiffness, a shear thinning fluid yielded the fastest swimming speed and that as the bending stiffness was increased, the shear thinning and shear thickening solutions collapsed towards the Newtonian solution. When considering this fact and the Newtonian solution, it seems that the most interesting parameter regime for problem as a whole is moderate or small bending stiffness coupled with a large amplitude forcing. It would be most desirable for us to find a numerical or analytical solution for the viscoelastic problem with a large amplitude forcing. This is not straightfor- ward, as the code that currently solves the viscoelastic problem cannot handle large amplitude forcing. Finding an analytical solution is also not straightfor- ward. To see why this is first note that from the Newtonian problem we have that px = uyy, (3.73) giving us that uy = τ1 − (1− y)px. (3.74) We can see from figure (3.9) that the strain rates in the Newtonian problem are getting larger, but it is not obvious that we can make the blanket assumption that they get large enough to simplify the constitutive relation.25 Even if the constitutive equation can be simplified, we would assume that the solution still has dependence on the viscoelastic parameters α and β, as did the solution for the large amplitude swimmer in chapter two. 25Please recall that in the Newtonian problem, strain rates and stresses are equal non- dimensionally, τ = γ̇ 72 Chapter 3. Elastic Swimming with Lubrication Theory 0 1 2 3 4 5 6 −40 −30 −20 −10 0 10 20 30 40 50 ξ=x−t τ 1 : Bl ac k, τ 0 : Cy an Stresses\\Strain Rates Figure 3.9: Strain rates on the swimmer (cyan) and on the wall (black) for a=100 and D=1.5. Note that the largest strain rates occur in the region where the plateau ends and the spike begins. 73 Chapter 4 Conclusions This thesis has presented an in depth discussion of swimmers near walls. Several different problems were analyzed and within these problems, different fluid types were considered. The mathematical setup of these problems can be applied to a broad spectrum of applications as it easily adapts to many similar situations. For example, similar problems that have different boundary conditions, wall locations, shapes, or constitutive models all require only simple changes to the basic formulation.26 Viscoelastic lubrication problems arise in many places, such as physiological and industrial applications, as well as other types of swimming problems [30, 1, 40, 41]. 4.1 Summary of Results To 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 the problem of a lubricated swimmer with a prescribed shape. It was concluded that shear thinning fluids are optimal for swimmers in this situation when com- pared to shear thickening and Newtonian fluids. This was found by considering the problem with a fixed rate of work and looking at the resulting dimensional swimming speeds. The two limiting cases of the distance between the wall and the swimmer approaching zero and infinity were considered analytically. It was found that the Newtonian swimming speed was returned in both situations. In the former limit, the Newtonian problem was returned entirely. In the latter limit, the solution, aside from the swimming speed and mass flux, was found to depend on the viscoelastic parameters α and β in a simple fashion. We then considered the lubrication problem with a swimmer of prescribed shape in a Bingham fluid. Here we fixed the amplitude of the swimmer and 26There are some restrictions to the scope of these differences, for example with changes in the constitutive model, we have to know τxy(uy), and that the law reduces to a steady-state shear flow for the thin-gap geometry with dominant shear stress. 74 Chapter 4. Conclusions solved 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 the region of the fluid where |τxy| ≤ τy. We found that there exist regions along the swimmer where the plug-flow reaches the swimmer (case II), reaches the wall (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 maximum swimming speed at B ≃ 1, which is greater than the Newtonian swimming speed. Combining the results from Chapter Two, we may speculate that it is the most advantageous for a swimmer to be moving through a fluid that exhibits a yields 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. The swimmer was modelled with the beam equation with an addition of an internal applied 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 yields an optimal swimming speed. As the amplitude of the forcing was increased, it was found that the swimmer shape began to develop a plateau region that approached the wall, followed by a large downward spike region. The value of the plateau region was successfully calculated analytically. We then introduced a viscoelastic fluid. It was again found that a shear thinning fluid is optimal for a swimmer as it yields the fastest swimming speeds. As the bending stiffness of the swimmer was increased, it was observed that the problem collapsed into the Newtonian problem. In general, we found that a shear thinning fluid or a fluid that exhibits a yield stress is theoretically optimal for swimmers near walls. Supporting this theory is the fact that polymeric fluids, the type of fluid found in most biological settings, 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 a swimmer to continue increasing its energy output, as it does not continue to increase the swimmer’s speed. This is because as the forcing along the swimmer increases, so does the stress within the fluid and so there exists a point in which the stresses in the fluid are so great that they begin to hamper movement. This suggests that in nature, swimmers do not exert themselves to their maximum, but rather moderately in such a way that maximizes speed. 75 Chapter 4. Conclusions 4.2 Further Work The results of this thesis suggest further work. Due to the formulation of the problem, and in particular, the simplification of the constitutive model when non-dimensionalized, it is theoretically straightforward to consider numerous different constitutive equations and fluid types. This suggests developing a more robust numerical method of solution. From a biological perspective, there are many extensions of the analysis that can be done. For example, it is the case that in a biological setting, swimmers are of finite length. Thus, it may be of interest to consider the problems with a finite swimmer and relax the periodicity conditions. There are numerous extensions that can be done with the elastic swimmer problem. 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 is unfortunate 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, one might consider the swimmer in a channel, therefore restricting the amplitude of the spike. There is potential for things to get quite interesting in this case. One might also be interested in looking into the large amplitude forcing viscoelastic problem, as it is not entirely clear that this will reduce into the Newtonian problem as it did in the similar case considered in Chapter Two. As mentioned throughout this thesis, we are by no means limited to the specific choices we have made. This is true for the choice of constitutive models, swimmer shapes (Chapter Two) and applied torques (Chapter Three). The mathematical set up we chose here lends itself nicely to generality, and exploring these potential alternate applications, via different choices of the constitutive model, swimmer shape and applied torques, is certainly something that should be looked into further. 76 Bibliography [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 swimming of 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 cornstarch suspensions. 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 Polymeric Liquids. 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 of Experimental 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 Fluid Mechanics, 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, pages 269–300. CRC Press, NY, 1990. 77 Bibliography [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 swimming filaments 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 and Physical 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, Mathematical and Physical Sciences, 217(1128):96–121, 1953. [16] J. J. L. Higdon. A hydrodynamic analysis of flagellar propulsion. Journal of 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 lubricated elastic sheet. Physical Review Letters, 93(13):137802, 2004. [19] R. E. Johnson and C. J. Brokaw. Flagellar hydrodynamics. a comparison between 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 of Reproduction, 25:931–937, 1981. [21] D . F. Katz. Characteristics of sperm motility. New York Academy of Sciences. 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 hamster spermatozoa 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. 78 Bibliography [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 the influence 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. General Publishing Company, fourth edition, 1944. [29] K. E. Machin. Wave propagation along flagella. Journal of Experimental Biology, 35:796–806, 1958. [30] A. Meziane, B. Bou-Sad, and John Tichy. Modelling human hip joint lubrication 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 Physical Sciences, 200(1063):523–541, 1950. [32] J. G. Oldroyd. Non-newtonian effects in steady motion of some idealized elastico-viscous liquid. Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 245(1241):278–297, 1958. [33] A. J. Reynolds. The swimming of minute organisms. Journal of Fluid Mechanics, 23:241–260, 1965. [34] S. S. Suarez and X. Dai. Hyperactivation enhances mouse sperm capacity for penetrating viscoelastic media. Biology of Reproduction, 46:686–691, 1992. [35] S. S. Suarez and H. C. Ho. Hyperactivated motility in sperm. Reproduction in 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 of Reproduction, 44:375–381, 1991. 79 Bibliography [37] G. I. Taylor. Analysis of the swimming of microscopic organisms. Proceed- ings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 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 swimming over a thin liquid layer. Journal of Fluid Mechanics, 601:25–61, 2008. [40] H. Winet. Ciliary propulsion of objects in tubes: Wall drag on swimming tetrahymena (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 myxobacteria glide. Current Biology, 12:369–377, 2002. [42] Y. L. Zhang, O. K. Matar, and R. V. Craster. Surfactant spreading on a thin weakly viscoelastic film. Journal of Non-Newtonian Fluid Mechanics, 105:53–78, 2002. 80 Appendix A Full Analysis of Parent Problem In this section we discuss a fundamental problem, that of a swimmer in the presence of a wall and a viscoelastic fluid without the lubrication approximation. Here we will assume that the amplitude of the waves is small in comparison to the mean distance to the wall. This is a nice warm up problem as it is a combination of Katz’s biharmonic problem and Lauga’s viscoelastic Taylor problem [22], [25]. A.1 Mathematical Formulation We assume that the sheet has prescribed motion, that of waves propagating to the right. In the frame moving with the swimmer, these waves look like y = a sin(kx− ωt). (A.1) The equations of motion for this problem are ▽ ·u = 0, (A.2) ▽p = ▽ · τ , (A.3) τ + λ1τ ▽ = µ ( γ̇ + λ2γ̇ ▽ ) , (A.4) which are the conservation of mass, momentum and the Oldroyd-B constitutive equation respectively. We also have the zero force condition27 〈τxy〉. = 0 (A.5) 27This equation still holds even though there is no lubrication approximation here due to the low amplitude compared to mean height approximation 81 Appendix A. Full Analysis of Parent Problem The next step is then to non-dimensionalize the equations as follows t = 1ω t̂, (x, y) = 1 k (x̂, ŷ), u = cû, γ̇ = ωˆ̇γ, (p, τ ) = µω(p̂, τ̂ ), where c = ωk is the speed of the travelling wave. The conservation of momentum and conservation of mass equations do not change: ▽ p = ▽ · τ , (A.6) ▽ · u = 0, (A.7) and the constitutive equation becomes: τ +De1τ ▽ = γ̇ +De2γ̇ ▽, (A.8) where De1 and De2 are Deborah numbers. It remains to non-dimensionalize the boundary conditions. Using the low amplitude approximation, we put ε = ak ≪ 1, (A.9) and find that (dropping hats): y = ε sin(x− t). (A.10) From which we can easily take the time derivative to find dy/dt = −ε cos(x− t) (A.11) The boundary conditions become, on the wall y = hk 28 u = U, (A.12) v = 0, (A.13) 28note that we just write U rather than (1/c)U for simplicity 82 Appendix A. Full Analysis of Parent Problem and on the swimmer y = ε sin(x− t), u = 0, (A.14) v = −ε cos(x− t). (A.15) A.2 Stream Function and Asymptotics Since the problem at hand is an incompressible, two dimensional problem, it is most useful to consider a stream function ψ such that u = ∂ψ ∂y , (A.16) v = − ∂ψ ∂x . (A.17) Automatically we have conservation of mass holding since ▽ · u = ux + vy = ∂ ∂x ( ∂ψ ∂y ) + ∂ ∂y ( − ∂ψ ∂x ) = 0. (A.18) It is convenient to observe that the boundary conditions can be written concisely as: ▽ ψ ∣∣ (x,ε sin(x−t)) = ε cos(x− t)i (A.19) ▽ψ ∣∣ (x,hk) = U j (A.20) where i = (1, 0) and j = (0, 1) are the unit vectors in the x and y directions respectively. 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 form U = εU1 + ε 2U2 + ... (A.21) ψ = εψ1 + ε 2ψ2 + ... (A.22) τ = ετ1 + ε 2τ2 + ... (A.23) That is, regular perturbation expansions in ε for the swimming speed, the stream function and the stress. Since it is assumed that the speed of the swimmer is constant, it must be 83 Appendix A. Full Analysis of Parent Problem the case that each Uk in the expansion of U is also constant. That is Uk = Const ∀k = 1, 2, 3... (A.24) Notice that since γ̇ = ( 2ux uy + vx uy + vx 2vy ) = ( 2ψyx ψyy − ψxx ψyy − ψxx −2ψxy ) , (A.25) it is convenient to write γ̇ = εγ̇1 + ε 2γ̇2 + ... (A.26) Similarly, we write u = εu1 + ε 2u2 + ... (A.27) In the spirit of asymptotics, the next step is to substitute the perturbation expansions into the equations of motion and the boundary conditions and sort out the first and second order equations. Before this is done, it is most useful to observe that ▽×▽p = 0, (A.28) and so pressure can be removed from the conservation of momentum equation. Summarizing the equations of motion then gives ▽×▽ ·τ = 0 (A.29) τ +De1τ ▽ = γ̇ +De2γ̇ ▽. (A.30) Substituting in the perturbation expansions yields at O(ε): ▽×▽ ·τ1 = 0, (A.31) τ1 +De1 ∂ ∂t τ1 = γ̇1 +De2 ∂ ∂t γ̇1, (A.32) and at O(ε2): ▽×▽ ·τ2 = 0, (A.33) τ2 +De1 ∂ ∂t τ2 +De1 ( u1 · ▽τ1 − τ1 · ▽u1 − (▽u1) T · τ1 ) = γ̇2 +De2 ∂ ∂t γ̇2 +De2 ( u1 · ▽γ̇1 − γ̇1 · ▽u1 − (▽u1) T · γ̇1 ) . (A.34) Now it is time to look at the boundary conditions. Substituting in the regular 84 Appendix A. Full Analysis of Parent Problem perturbation expansion for ψ and expanding about y = 0 yields at O(ε): ▽ ψ1 ∣∣ (x,0) = cos(x− t)i, (A.35) ▽ψ1 ∣∣ (x,hk) = U1j, (A.36) and at O(ε2): ▽ ψ2 ∣∣ (x,0) +▽ ( ∂ψ1 ∂y ) ∣∣∣ (x,0) sin(x− t) = 0, (A.37) ▽ψ2 ∣∣ (x,hk) = U2j. (A.38) Substituting the expansions into the constraint 〈τxy〉 = 0 yields that 〈τkxy〉 = 0 k = 1, 2, 3... (A.39) A.3 Leading Order Equations If one stares at the leading order equations for a moment, ▽×▽ ·τ 1 = 0, (A.40) τ 1 +De1 ∂ ∂t τ 1 = γ̇1 +De2 ∂ ∂t γ̇1, (A.41) it becomes apparent that the easiest thing to do is to take the divergence and then the curl of the second equation thus eliminating leading order stress. By doing this we obtain ( 1 +De2 ∂ ∂t ) ▽×▽ ·γ̇1 = 0. (A.42) When substituting in γ̇1 = γ̇1(ψ1xx, ψ1xy, ψ1yy), we actually get ( 1 +De2 ∂ ∂t ) ▽4 ψ1 = 0 (A.43) where ▽4ψ = ψxxxx + 2ψxxyy + ψyyyy. 85 Appendix A. Full Analysis of Parent Problem Writing down the boundary conditions in a more practical fashion gives ∂ψ1 ∂x ∣∣ (x,0) = 0, ∂ψ1∂y ∣∣ (x,0) = 0, ∂ψ1 ∂x ∣∣ (x,hk) = cos(x− t), ∂ψ1∂y ∣∣ (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 conditions ψ1 = f(y)sin(x− t), (A.45) where f(y) is to unknown. Substituting this into the leading order equation for ψ yields d4f dy4 − 2 d2f dy2 + f = 0. (A.46) It is easily found that f(y) = (a+ by)e−y + (c+ dy)ey, (A.47) where constants a, b, c, and d are found by applying the boundary conditions and will depend on h. When applying the boundary condition ψ1y ∣∣ (x,hk) = U1 it becomes apparent that in order for the assumption that U = constant to hold, it must be the case that U1 = 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 done previously by others, who all found the leading order swimming speed to be zero [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. Put ψ1(y) = [(A+BH) cosh(H) + (C +DH) sinh(H)] sin(x− t). This gives ψ1x = [(A+BH) cosh(H)+ +(C +D(H)) sinh(H)] cos(x− t), ψ1y = −[(A+BH +D) sinh(H)+ +(B + C +DH) cosh(H)] sin(x− t), 86 Appendix A. Full Analysis of Parent Problem where H = (hk − y). Upon applying the boundary conditions, for y = hk: A cos(x− t) = 0, −(B + C) sin(x− t) = U1. Clearly, in order for the assumption, U1 = Const, to hold, it must be the case that U1 = 0, thus implying that B = −C. Also, it is clear from the first condition that A = 0. Thus so far, ψ1 = [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 the problem (‘complexification’, see below). That is put ψ1 = 1 2 ( ψ̃1e i(x−t) + ψ̃1 ∗ e−i(x−t) ) , (A.49) and consider only ψ̃1e i(x−t). One then gets: ψ̃1 = [B̃[sinh(hk − y)− (hk − y) cosh(hk − y)] + D̃(hk − y) sinh(hk − y)]. Define Ψ1(y) = (hk − y) sinh(hk − y), (A.50) Ψ2(y) = sinh(hk − y)− (hk − y) cosh(hk − y), (A.51) thus giving ψ̃1 = [D̃Ψ1 + B̃Ψ2]. Now, the unknown coefficients B̃ and D̃ must be determined using the other boundary conditions. On y = 0: ( ψ̃1e i(x−t) ) x ∣∣∣ y=0 = iei(x−t)[D̃Ψ1(0) + B̃Ψ2(0)] = e i(x−t),( ψ̃1e i(x−t) ) y ∣∣∣ y=0 = ei(x−t)[D̃Ψ′1(0) + B̃Ψ ′ 2(0)] = 0. Solving these gives D̃ = i Ψ′2(0) W (0) , B̃ = −i Ψ′1(0) W (0) , 87 Appendix A. Full Analysis of Parent Problem where W (y) = Ψ′1(0)Ψ2(y)−Ψ ′ 2(0)Ψ1(y). (A.52) Putting this altogether one has ψ̃1 = −i W (y) W (0) (A.53) or writing down the real component, ψ1 = W (y) W (0) sin(x− t) (A.54) Although this second method of solving for the leading order stream function and swimming speed is longer, it provides a more concise final form. This is useful since the leading order swimming speed is zero, so the calculation must be done to second order. This aesthetically pleasing form of the leading order solution will make the second order problem much easier and less algebraically challenging 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 leading order equation τ̃1 +De1∂tτ̃1 = ˜̇γ1 +De2∂t ˜̇γ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) τ1 = ℜ[τ̃1e i(x−t)], γ̇1 = ℜ[ ˜̇γ1e i(x−t)]. (A.55) Alternately, one could write τ1 = 1 2 ( τ̃1e i(x−t) + τ̃1 ∗e−i(x−t) ) , γ̇1 = 1 2 ( ˜̇γ1e i(x−t) + ˜̇γ1 ∗e−i(x−t) ) , (A.56) observing that τ̃1 and ˜̇γ1 are functions of y only, and are independent of x and t (note that this is also true for ψ̃1 defined above). This ‘complexification’ of the problem implies that ∂t = −∂x = −i. 88 Appendix A. Full Analysis of Parent Problem Therefore, by exploiting this simplification of time derivatives, we can easily obtain τ̃1 = 1− iDe2 1− iDe1 ˜̇γ1. (A.57) Upon substituting in ˜̇γ1 = ˜̇γ1(ψ̃1xx, ψ̃1xy, ψ̃1yy) as defined above, we get τ̃1 = (1− iDe2)e i(x−t) (1− iDe1)W (0) ( 2W ′ −i(W ′′ +W ) −i(W ′′ +W ) −2W ′ ) . (A.58) Notice that the condition, 〈τ1xy〉 = 0, automatically holds due to the fact that τ1xy ∝ γ̇1xy ∝ sin(x− t). A.4 Second Order Equations The leading order problem has been solved and it has been shown that in order to maintain the assumption that the speed of the swimmer is constant, at leading order it must be equal to zero as predicted. Thus, in order to calculate the swimming 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 due to the nonlinear terms of the upper convected derivative that crop up. Let us restate the second order problem: ▽×▽ · τ2 = 0, ( 1 +De1 ∂ ∂t ) τ2 − ( 1 +De2 ∂ ∂t ) γ̇2 = De2 ( u1 · ▽γ̇1 − γ̇1 · ▽u1 − (▽u1) T · γ̇1 ) − −De1 ( u1 · ▽τ1 − τ1 · ▽u1 − (▽u1) T · τ1 ) , with boundary conditions: ▽ψ2 ∣∣ (x,0) +▽ ( ∂ψ1 ∂y ) ∣∣∣ (x,0) sin(x− t) = 0, ▽ψ2 ∣∣ (x,hk) = U2j. 89 Appendix A. Full Analysis of Parent Problem The zero force condition states that 〈τ2xy〉 = 0. (A.59) To start things off slowly, consider first the boundary conditions. Writing things in a more user friendly fashion gives on y = 0: ψ2x ∣∣ y=0 = − sin(x− t)ψ1yx, ψ2y ∣∣ y=0 = − sin(x− t)ψ1yy. But we know ψ1 from solving the leading order problem, which gives ψ1yx = W ′(y) W (0) cos(x − t) recalling that W (y) = Ψ ′ 1(0)Ψ2(y) − Ψ ′ 2(0)Ψ1(y) and Ψ1(y) = (hk − y) sinh(hk − y),Ψ2(y) = sinh(hk − y) − (hk − y) cosh(hk − y). Thus, we obtain ψ2x ∣∣ y=0 = 0, (A.60) ψ2y ∣∣ y=0 = − sin2(x− t) W ′′(0) W (0) . (A.61) On y = hk: ψ2x ∣∣ y=hk = 0, (A.62) ψ2y ∣∣ y=hk = U2. (A.63) Inspecting the second order equations of motion as written above reveals that the right-hand side of the momentum equation is composed of leading order terms exclusively. Since it is the case that the leading order stress can be simply represented in terms of the leading order rate of strain when complexified, let us put the right-hand side of the momentum equations in terms of the complexified variables u1 = 1 2 (ũ1e i(x−t) + ũ∗1e −i(x−t)), γ̇1 = 1 2 ( ˜̇γ1e i(x−t) + ˜̇γ1 ∗e−i(x−t)), and of course τ̃1 = 1− iDe2 1− iDe1 ˜̇γ1. 90 Appendix A. Full Analysis of Parent Problem We obtain RHS = (De2 −De1) 4(1− iDe1) [ ũ∗1 · ▽ ˜̇γ1 − ˜̇γ1 · ▽ũ ∗ 1 − (▽ũ ∗ 1) T · ˜̇γ1 ] + C.C.+ + (De2 −De1) 4(1− iDe1) e2i(x−t) [ ũ1 · ▽ ˜̇γ1 − ˜̇γ1 · ▽ũ1 − (▽ũ1) T · ˜̇γ1 ] + C.C., (A.64) where C.C. refers to the complex conjugate of the term directly preceding it.29 To solve this rather long, complicated looking equation, make use of the periodicity of the problem and look at the x-average of the equation. Notice at this point one could take the curl and divergence of the problem to remove the τ2 terms, but that is rather cumbersome. Rather, only consider the xy-component and utilize the zero force condition to simplify the problem. Let us first look at the boundary conditions. Note that 〈ψ2x〉 = 1 2piψ2 ∣∣2pi 0 = 0 always, so the boundary conditions for ψ2x yield no information. However, the other two boundary conditions have non-trivial averages: 〈ψ2y〉 ∣∣ y=0 = − W ′′(0) 2W (0) , (A.66) 〈ψ2y〉 > ∣∣ y=hk = U2. (A.67) Now let us average the second order constitutive equation and look at the xy- component (recall that ∂t = −∂x and that 〈e ±2i(x−t)〉 = 0), subscripts here referring to components: 〈τ2xy〉 − 〈γ̇2xy〉 = ℜ { (De2 −De1) 2(1− iDe1) [ ũ∗1 · ▽ ˜̇γ1 − ˜̇γ1 · ▽ũ ∗ 1 − (▽ũ ∗ 1) T · ˜̇γ1 ] xy } . (A.68) Applying the zero force condition and writing the rate of strain component in 29Alternately, one could write RHS = ℜ { (De2 −De1) 2(1− iDe1) [ ũ ∗ 1 · ▽ ˜̇γ1 − ˜̇γ1 · ▽ũ ∗ 1 − (▽ũ ∗ 1) T · ˜̇γ1 ]} + + ℜ { (De2 −De1) 2(1− iDe1) e2i(x−t) [ ũ1 · ▽ ˜̇γ1 − ˜̇γ1 · ▽ũ1 − (▽ũ1) T · ˜̇γ1 ]} (A.65) 91 Appendix A. Full Analysis of Parent Problem terms of the stream function gives: 〈ψ2yy − ψ2xx〉 = ℜ { (De1 −De2) 2(1− iDe1) [ ũ∗1 · ▽ ˜̇γ1 − ˜̇γ1 · ▽ũ ∗ 1 − (▽ũ ∗ 1) T · ˜̇γ1 ] xy } . (A.69) But, 〈ψ2xx〉 >= 1 2piψ2x ∣∣2pi 0 = − 12piv ∣∣2pi 0 = 0 due to the assumed periodicity of the solution. Since the bounds of the average are not dependent on y, we can put 〈ψ2〉yy = ℜ { (De1 −De2) 2(1− iDe1) [ ũ∗1 · ▽ ˜̇γ1 − ˜̇γ1 · ▽ũ ∗ 1 − (▽ũ ∗ 1) T · ˜̇γ1 ] xy } . (A.70) Thus, to solve the problem and calculate the swimming speed we must simply integrate once to obtain an expression for 〈ψ2〉y 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 = −i W (y) W (0) = −i W W0 , we can write down the terms on the right-hand side in a relatively concise manner. Namely, we obtain [ ũ1 ∗ · ▽ ˜̇γ1 ] xy = i W 20 (W ′W ′′ + 2WW ′ +WW ′′′) , (A.71) [ ˜̇γ1 · ▽ũ1 ∗ ] xy = i W 20 (W ′W ′′ + 3WW ′) , (A.72) [ ▽ũ1 ∗T · ˜̇γ1 ] xy = −i W 20 (3W ′W ′′ +WW ′) . (A.73) Substituting these into the governing equation gives 〈ψ2〉yy = ℜ { (De1 −De2) 2(1− iDe1) i W 20 [WW ′′′ + 3W ′W ′′] } . (A.74) But, notice that [WW ′′′+3W ′W ′′] can be written as [WW ′′+W ′2]′ and so we find that 〈ψ2〉yy = De1(De2 −De1) 2(1 +De21)W 2 0 [ WW ′′ +W ′2 ]′ . (A.75) Integrating and applying the first boundary condition gives − W ′′0 2W0 = De1(De2 −De1) 2(1 +De21)W 2 0 [ W0W ′′ 0 +W ′2 0 ] +Const. Using the fact that W ′0 = 0 the constant of integration can be easily solved for, Const = − W ′′0 2W0 (1 +De1De2) (1 +De21) . 92 Appendix A. Full Analysis of Parent Problem Thus, 〈ψ2〉y = De1(De2 −De1) 2(1 +De21)W 2 0 [ WW ′′ +W ′2 ] − W ′′0 2W0 (1 +De1De2) (1 +De21) , (A.76) and so finally one has the swimming speed U2 = 〈ψ2〉y ∣∣ y=hk = − W ′′0 2W0 (1 +De1De2) (1 +De21) , (A.77) since W (hk) = 0 and W ′(hk) = 0. A.5 Work The leading order rate of work (per unit wavelength) is given by the volume integral of 〈τ1 : γ̇1〉, (A.78) and is equal to the rate of viscous dissipation. To simplify this calculation, we look again to the ‘complexified’ variables. That is, put τ1 : γ̇1 = 1 4 (τ̃1 + τ̃1 ∗) : ( ˜̇γ1 + ˜̇γ1 ∗ ) . Let us calculate the integral in x first. It is the case that 〈τ̃1 : ˜̇γ1〉 ∝ 〈e 2i(x−t)〉 = 0, 〈τ̃1 ∗ : ˜̇γ1 ∗〉 ∝ 〈e−2i(x−t)〉 = 0. Since both τ1 and γ̇1 are symmetric, we can put 〈τ1 : γ̇1〉 = 1 2 ℜ { 〈τ̃1 : ˜̇γ1 ∗〉 } Knowing the form of τ̃1 in terms of ˜̇γ1 gives 〈τ1 : γ̇1〉 = 1 2 ℜ {( 1− iDe2 1− iDe1 ) 〈 ˜̇γ1 : ˜̇γ1 ∗〉 } Using the dependence of γ̇1 on ψ1, and the fact that ψ̃1 = −i W W0 , gives us 〈τ1 : γ̇1〉 = 1 +De1De2 1 +De21 ( 4W ′2 + (W +W ′′)2 W 20 ) . (A.79) 93 Appendix A. Full Analysis of Parent Problem Thus, leading order rate of work per unit wavelength, is given by (1 +De1De2) (1 +De21)W 2 0 ∫ hk 0 [ 4W ′2 + (W +W ′′)2 ] dy. (A.80) A.6 Limiting Cases A.6.1 Distance Between Wall and Swimmer Approaches Infinity To make sure that this problem agrees with the infinite case considered by Lauga in [25], we must take the limit as the distance between the wall and the swimmer approaches infinity. It is known that ψ1 = W (y) W (0) sin(x− t), recalling that W (y) = Ψ′1(0)Ψ2(y)−Ψ ′ 2(0)Ψ1(y), and Ψ1(y) = (hk − y) sinh(hk − y), Ψ2(y) = sinh(hk − y)− (hk − y) cosh(hk − y). The limit as h −→∞ is easily calculated. It is found that lim h−→∞ W (y) W (0) = (1 + y)e−y Thus lim h−→∞ ψ1 = (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 infinite case that the limit of the swimming speed will as well. But to illustrate this clearly, recall that U2 = − W ′′(0) 2W (0) (1 +De1De2) (1 +De21) . So, the limit lim h−→∞ U2 = 1 2 (1 +De1De2) (1 +De21) , which again agrees with the infinite case. 94 Appendix A. Full Analysis of Parent Problem Similarly, we can calculate the limit for the averaged rate of work done by the swimmer. It is found that lim h−→∞ { (1 +De1De2) (1 +De21)W 2 0 ∫ hk 0 [ 4W ′2 + (W +W ′′)2 ] dy } = 2 (1 +De1De2) (1 +De21) . We cannot directly compare this limit with the form of the solution found in [25] as it is stated as Work∞ WN = (1 +De1De2) (1 +De21) , and WN is not stated. However, we can easily see that the calculation of work in the infinite case is exactly the same with the exception of the W function. In the infinite case one has W∞ = (1 + y)e −y, and so it is found that (1 +De1De2) (1 +De21)W∞(0) 2 ∫ ∞ 0 [ 4W ′2∞ + (W∞ +W ′′ ∞) 2 ] 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 organism swimming in an Oldroyd-B viscoelastic fluid near a wall are consistent with those of the problem of an organism swimming in an infinite fluid of the same type. A.6.2 Swimmer Approaches the Wall Next, we would like to consider the limit as the swimmer and the wall approach each other. There are two ways of thinking of this limit. The first is keeping the wall fixed at y = h and increasing the amplitude of the swimmer until it reaches the wall. This line of thinking contradicts the assumption that ak << 1, as it implies that O(a) → O(h) = O(1/k). This would make our present solution invalid. The second way of thinking about this limit is to think of the wall approaching the swimmer, whose amplitude and wavelength remains the fixed. This implies that O(h)→ O(a), still maintaining the assumption that ak << 1. In fact, in this manner of thinking we get that hk << 1 which is exactly the 95 Appendix A. Full Analysis of Parent Problem lubrication approximation. We have that, to leading order U = −ε2 (1 +De1De2) 2(1 +De21) F (hk), (A.81) where F (hk) = W ′′0 W0 = (hk)2 + sinh2(hk) (hk)2 − sinh2(hk) . (A.82) In our limit, hk → ak = ε and so we have that F (ε) = ε2 + sinh2 ε ε2 − sinh2 ε . (A.83) Expanding sinh ε = ε+ 16ε 3 + . . . and simplifying gives us F (ε) = 2ε2 ++ 13ε 4 + . . . − 13ε 4 − . . . . (A.84) Thus, in the limit when the wall approaches the swimmer, we find that the swimming speed is U = 3 (1 +De1De2) (1 +De21) . (A.85) 96 Appendix B Details of Bingham Analysis B.1 Case I Beginning 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 equation using equation (2.114). Integrating, ∫ y Y uydy = ∫ y Y (y − Y+)pxdy, (B.1) gives u = px 2 ( (y − Y+) 2 − (Y − Y+) 2 ) . (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 rearranging to solve for Y+ gives, Y+ = 1 2 (1 + Y )− U px(1− Y ) . (B.3) Technically, the plug region is outside the boundaries we are considering, but, for the purpose of solving this set of equations, let us imagine that the flow profile continues through y = Y . Since the flow profile is a concave-up parabola, we can see that up < 0. As well, applying the condition on y = Y+ gives us that up = − px 2 (Y − Y+) 2. (B.4) In turn, this tells us sgn(up) = −sgn(px). (B.5) Thus, we find that sgn(px) > 0, (B.6) 97 Appendix B. Details of Bingham Analysis and so σ = +1. We can now look at equation (2.119). Substituting in u and simplifying, gives q = Y + U 2 (1− Y )− px 12 (1− Y )3. (B.7) B.2 Case II For case II, we have Y− < Y < Y+ < 1. Let us begin with the region where Y− < y < Y+. Since we have that, u|(y=Y ) = 0, and Y− < Y < Y+, we know automatically that up = 0. (B.8) Now, consider the region where Y+ < y < 1. Taking equation (2.108), simplify- ing using equation (2.114), and integrating, ∫ y Y+ uydy = ∫ y Y+ (y − Y+)pxdy, (B.9) gives us u = px 2 (y − Y+) 2. (B.10) Applying the boundary condition at y = 1 and solving for Y+ gives, Y+ = 1− √ 2U px . (B.11) Since we are assuming that U > 0, the above equation implies px > 0. (B.12) This implies that σ = +1. Moving on to equation (2.119), q = Y + ∫ 1 Y+ udy. (B.13) Substituting in our expression for u and simplifying, gives q = Y + px 6 (1− Y+) 3. (B.14) 98 Appendix B. Details of Bingham Analysis B.3 Case III Case III is the case when the pseudo-plug is contained within the boundaries of 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 the upper region, Y+ < y < 1, we can again simplify equation (2.108) using equation (2.114) and integrate, ∫ y Y+ uydy = ∫ y Y+ (y − Y+)pxdy. (B.15) This gives us that u = up + px 2 (y − Y+) 2. (B.16) Applying the boundary condition on y = 1 and solving for up yields, up = U − px 2 (1− Y+) 2. (B.17) Now, consider the lower yielded region, Y < y < Y−. In this region, we simplify equation (2.108) using equation (2.115). Integrating, ∫ Y− y uydy = ∫ Y− y (y − Y−)pxdy, (B.18) gives u = up + px 2 (y − Y−) 2. (B.19) Applying the boundary condition at y = Y and solving for up gives, up = − px 2 (Y − Y−) 2. (B.20) Equating these two equations for up and solving for Y+ gives, Y+ = ( 1− Y 2 − 2Upx − 4 B |px| ( Y + B|px| )) 2(1− Y − 2 B|px| ) (B.21) recalling that, Y− = Y+ − 2 B |px| , from equation (2.117). We can now look at equation (2.119), q = Y + ∫ 1 Y+ udy + up(Y+ − Y−) + ∫ Y− Y udy. (B.22) 99 Appendix B. Details of Bingham Analysis Substituting in u and simplifying, gives q = Y + px 6 ( (1− Y+) 3 + (Y− − Y ) 3 − 3(1− Y )(Y− − Y ) 2 ) . (B.23) B.4 Case IV For case IV, we have Y < Y− < 1 < Y+. Let us begin with the region where Y− < y < Y+. Since we have that, u|(y=1) = U , and Y− < 1 < Y+, we know automatically that up = U. (B.24) Now, consider the region where Y < y < Y−. Taking equation (2.108), simpli- fying using equation (2.115), and integrating, ∫ Y− y uydy = ∫ Y− y (y − Y−)pxdy, (B.25) gives us u = U + px 2 (y − Y−) 2. (B.26) Applying the boundary condition at y = Y and solving for Y− gives, Y− = Y + √ −2U px . (B.27) Since we are assuming that U > 0, we know that, in this case, px < 0. (B.28) This implies that σ = −1. Looking at equation (2.119), q = Y + U(1− Y−) + ∫ Y− Y udy, (B.29) substituting in u and simplifying, gives q = Y + px 6 (Y− − Y ) 3 + U(1− Y ). (B.30) 100 Appendix B. Details of Bingham Analysis B.5 Case V Finally, 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) and integrating, ∫ y Y uydy = ∫ y Y (y − Y−)pxdy, (B.31) gives u = px 2 ( (y − Y−) 2 − (Y − Y−) 2 ) . (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 rearranging to solve for Y− gives, Y− = 1 2 (1 + Y )− U px(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 imagine that 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 condition on y = Y− gives us that up = − px 2 (Y − Y−) 2. (B.34) Again, this tells us, sgn(up) = −sgn(px), (B.35) and so we find that sgn(px) < 0, (B.36) which gives us that σ = −1. We can now look at equation (2.119). Substituting in u and simplifying, gives q = Y + U 2 (1− Y )− px 12 (1− Y )3. (B.37) Note that for all the above cases, equation (2.117) allows us to put Y− in terms Y+, or vice versa, Y− = Y+ − 2B |px| . (B.38) 101 Appendix C Code For Chapter Two Below is the MATLAB code used to obtain the results in Chapter Two. C.1 Viscoelastic Problem Code % NEWTON’S METHOD clear N = 50; h = 2*pi/N; xi = 0:h:2*pi; % uniform mesh DimWork=1; % Fixed Dimensional Work mu=1; % viscosity hh=0.1; % mean depth for k=1:199 ah = k*0.005; % a/h < 1; disp([’a/h = ’,num2str(ah)]) Y = ah*sin(xi); % REGULAR %Y = ah*tanh(12*sin(xi)); % SQUARE WAVE a0=0; b0=0; % Newtonian - could have a=b=anything.... a1=1; b1=2; % Thinning a2=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); 102 Appendix C. Code For Chapter Two q0=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));; % Thinning W1Tprime1=@(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+1 Int1(j)=quadl(W1Tprime1,gY1(j),g11(j)); end I1=Int1./(p_x1); II1=I1(1:N); w1(k)=sum(II1.*diff(xi)); % Thickening W2Tprime2=@(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+1 Int2(j)=quadl(W2Tprime2,gY2(j),g12(j)); end I2=Int2./(p_x2); II2=I2(1:N); w2(k)=sum(II2.*diff(xi)); 103 Appendix C. Code For Chapter Two % must calculate wave speed ... c0(k)=DimWork*hh^2/(2*mu*w0(k)); % Newtonian c1(k)=DimWork*hh^2/(2*mu*w1(k)); % Thinning c2(k)=DimWork*hh^2/(2*mu*w2(k)); % Thickening uu0(k)=c0(k)*U0; uu1(k)=c1(k)*U1; uu2(k)=c2(k)*U2; end figure plot(AH,uu0,’r’,AH,uu1,’b’,AH,uu2,’g’) title([’Dimensional Swimming Speed’],’FontSize’,12) xlabel(’a/h’,’FontSize’,12) ylabel(’U*c’,’FontSize’,12) 104 Appendix C. Code For Chapter Two function [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/dgY I0prime = @(g) g.*(1-b*g.^2+3*a*g.^2+a*b*g.^4)./(1+b*g.^2).^2; % Observing that -dI0/Dg1 = dI0/dgY J = zeros(2*N+4); % most entries are zero U=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 guess PP = (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+1 I1 = quadl(I1prime,gamY(j),gam1(j)); I0 = quadl(I0prime,gamY(j),gam1(j)); % F3’s F(j+2) = U*(feval(T,gam1(j)) 105 Appendix 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’s F(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)/dxi J(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)/dxi J(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)); end for j=1:N % Constraints % dF1/Dxi J(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/Dxi J(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)); end y = -J\\F’; 106 Appendix C. Code For Chapter Two Guess = Guess + y’; U=Guess(1); q=Guess(2); gam1=Guess(3:N+3); gamY=Guess(N+4:2*N+4); if norm(y,inf) < E disp([’Method converging in ’,num2str(i),’ iterations’]) break end end if i==maxit disp([’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); 107 Appendix C. Code For Chapter Two C.2 Bingham Problem Code %ssdriv a = 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; if conno<1e-12 break end end subplot(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) 108 Appendix C. Code For Chapter Two subplot(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/2 pi/2 min(Y) 1]) xlabel(’x’, ’Fontsize’, 22) ylabel(’Y_+ and Y_-’, ’Fontsize’, 22) hold off legend(’Pseudo-plug’,’Y_+’,’Y_-’,’Y’) 109 Appendix C. Code For Chapter Two % ssrhs x=([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); 110 Appendix C. Code For Chapter Two for 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; end 111 Appendix D Code For Chapter Three Below is the MATLAB code used to obtain the results in Chapter Three. D.1 Newtonian Problem Code First, 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.parameters qqq=sol.parameters(1); uuu=sol.parameters(2); aaa=aa; for na=1:25 aa=na*4 sol = bvp4cc(@solode,@solbc,sol,options); hold on,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) hold on,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) axis tight axes(’position’,[0.2 0.2 0.28 0.28]) plot(aaa,uuu,’r’,aaa,qqq,’g’) xlabel(’a’) axis tight figure(2) plot3(1-sol.y(1,:),-sol.y(2,:),-sol.y(3,:)) 112 Appendix D. Code For Chapter Three view([-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 tight function yinit = solinit(x) global aa yinit = [ 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*2 113 Appendix D. Code For Chapter Three yb(7) ]; function dydx = solode(x,y,lambda) global aa DD q=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)/DD y(1,:) 6*(1-q)./y(1,:).^2+(4*U-6)./y(1,:)]; 114 Appendix D. Code For Chapter Three D.2 Viscoelastic Problem Code Next, we list the code for the Viscoelastic problem with moderate amplitude. % siniC D =0.2; %7;%27; % Newtonian alf = 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 tight title([’\\alpha=\\beta=1 and D=’,num2str(D)],’FontSize’,12) xlabel(’x’,’FontSize’,12) ylabel(’Y’,’FontSize’,12) % Shear Thinning alf = 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 tight title([’\\alpha=1, \\beta=2, D=’,num2str(D)],’FontSize’,12) xlabel(’x’,’FontSize’,12) ylabel(’Y’,’FontSize’,12) % Shear Thickening alf = 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; 115 Appendix D. Code For Chapter Three figure(1), subplot(223) sdrivC figure(1), subplot(223) figure(1), hold on, axis tight title([’\\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 off xlabel(’a’,’FontSize’,12) ylabel(’U’,’FontSize’,12) title ([’U vs a for D=’,num2str(D)],’FontSize’,12) axis tight 116 Appendix D. Code For Chapter Three %sdrivC a=0.01 N=55; x=([0:N-1]+1/2)/N*pi*2; sdrivCC for niter=1:10 srhsC dX=-EJ\\E; conno=max(abs(dX)); X=X+dX; if conno<1e-12 disp(’Converging 1’) break end end disp([’conno = ’, num2str(conno)]) if (alf == 1) & (bet == 1) plot(x,Y,’r’) end if (alf == 1) & (bet == 2) %shear thinning plot(x,Y,’b’) end if (alf == 2) & (bet == 1) %shear thickening plot(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 off UU=X(3*N+1); qq=X(3*N+2); [UU qq] a0=a; aa=a0; D0=D; DD=D0; for na = 1:12 a = a0 +na/10 117 Appendix D. Code For Chapter Three for niter=1:15 srhsC dX=-EJ\\E; conno=max(abs(dX)); X=X+dX; if conno<1e-8 disp(’converging 2’) break end end disp([’conno = ’, num2str(conno)]) UU=[UU X(3*N+1)]; qq=[qq X(3*N+2)]; aa = [aa a]; hold on if (alf == 1) & (bet == 1) plot(x,Y,’r’), axis tight end if (alf == 1) & (bet == 2) %shear thinning plot(x,Y,’b’), axis tight end if (alf == 2) & (bet == 1) %shear thickening plot(x,Y,’g’), axis tight end hold off end 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 off end if (alf == 1) & (bet == 2) %shear thinning figure(2),plot(aa,qq,’b’) figure(1), subplot(224),hold on, plot(aa,UU,’b’), hold off end if (alf == 2) & (bet == 1) %shear thickening figure(2),plot(aa,qq,’g’) 118 Appendix D. Code For Chapter Three figure(1), subplot(224), hold on, plot(aa,UU,’g’) , hold off end figure(1) 119 Appendix D. Code For Chapter Three % srhsC gamY = 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; DiffY E=zeros(size(X)); EJ=zeros(3*N+2,3*N+2); % F4 E(3*N+1)=sum(T1); % F5 E(3*N+2)=sum(Y); for n=1:N I0(n)=quadl(I0prime,gamY(n),gam1(n)); %F1 E(n)=I0(n)-U*px(n); I1(n)=quadl(I1prime,gamY(n),gam1(n)); %F2 E(N+n)=I1(n)-px(n)^2.*(U+Y(n)-q-U*Y(n))-U*TY(n)*px(n); %F3 E(2*N+n) = px(n)-D3f(n)-D*DY5(n); 120 Appendix D. Code For Chapter Three % DF1 EJ(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); % DF2 EJ(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; % DF3 EJ(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 % DF4 EJ(3*N+1,N+n) = Tp1(n); % DF5 EJ(3*N+2,2*N+n) = 1; end % Must do dF3/dY separately due to discretization % A is differentiation matrix for DY5 EJ(2*N+1:3*N,2*N+1:3*N) =EJ(2*N+1:3*N,2*N+1:3*N)-D*A; 121 Appendix D. Code For Chapter Three This code calculated the fifth derivative of a function given by a specific set of points numerically. It assumes periodicity of the function. It is assumed that the mesh used does not include its endpoints. %DiffY h = x(2)-x(1); %uniform mesh A = 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; 122 Appendix D. Code For Chapter Three This is the code that formulates the initial guess for the Newton’s Iteration used to solve the viscoelastic problem. % srdrivCC x=([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; if conno<1e-12 break end end for niter=1:10 srhsCC dX=-EJ\\E; conno=max(abs(dX)); X=X+dX; if conno<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; 123 Appendix D. Code For Chapter Three % srhsCC gamY = 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:N I0(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); end 124"@en ; edm:hasType "Thesis/Dissertation"@en ; vivo:dateIssued "2008-11"@en ; edm:isShownAt "10.14288/1.0066549"@en ; dcterms:language "eng"@en ; ns0:degreeDiscipline "Mathematics"@en ; edm:provider "Vancouver : University of British Columbia Library"@en ; dcterms:publisher "University of British Columbia"@en ; dcterms:rights "Attribution-NonCommercial-NoDerivatives 4.0 International"@en ; ns0:rightsURI "http://creativecommons.org/licenses/by-nc-nd/4.0/"@en ; ns0:scholarLevel "Graduate"@en ; dcterms:title "Swimming in slime"@en ; dcterms:type "Text"@en ; ns0:identifierURI "http://hdl.handle.net/2429/1503"@en .