Waving transport and propulsion in a generalized Newtonian fluid
aa r X i v : . [ phy s i c s . f l u - dyn ] M a r Waving transport and propulsion in a generalizedNewtonian fluid
J. Rodrigo V´elez-Cordero a , Eric Lauga b a Instituto de Investigaciones en Materiales, Universidad Nacional Aut´onoma de M´exico, Apdo. Postal70-360 M´exico D.F. 04510. b Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Centre forMathematical Sciences, Wilberforce Road, Cambridge, CB3 0WA, UK.
Abstract
Cilia and flagella are hair-like appendages that protrude from the surface of a variety ofeukaryotic cells and deform in a wavelike fashion to transport fluids and propel cells.Motivated by the ubiquity of non-Newtonian fluids in biology, we address mathemat-ically the role of shear-dependent viscosities on both the waving flagellar locomotionand ciliary transport by metachronal waves. Using a two-dimensional waving sheet asmodel for the kinematics of a flagellum or an array of cilia, and allowing for both nor-mal and tangential deformation of the sheet, we calculate the flow field induced by asmall-amplitude deformation of the sheet in a generalized Newtonian Carreau fluid upto order four in the dimensionless waving amplitude. The net flow induced at far fromthe sheet can be interpreted either as a net pumping flow or, in the frame moving withthe sheet, as a swimming velocity. At leading order (square in the waving amplitude),the net flow induced by the waving sheet and the rate of viscous dissipation is the sameas the Newtonian case, but is di ff erent at the next nontrivial order (four in the wavingamplitude). If the sheet deforms both in the directions perpendicular and parallel tothe wave progression, the shear-dependence of the viscosity leads to a nonzero flow in-duced in the far field while if the sheet is inextensible, the non-Newtonian influence isexactly zero. Shear-thinning and shear-thickening fluids are seen to always induce op-posite e ff ects. When the fluid is shear-thinning, the rate of working of the sheet againstthe fluid is always smaller than in the Newtonian fluid, and the largest gain is obtainedfor antiplectic metachronal waves. Considering a variety of deformation kinematics forthe sheet, we further show that in all cases transport by the sheet is more e ffi ciency in ashear-thinning fluid, and in most cases the transport speed in the fluid is also increased.Comparing the order of magnitude of the shear-thinning contributions with past workon elastic e ff ects as well as the magnitude of the Newtonian contributions, our theoret-ical results, which beyond the Carreau model are valid for a wide class of generalizedNewtonian fluids, suggest that the impact of shear-dependent viscosities on transportcould play a major biological role. Preprint submitted to Elsevier September 21, 2018 . Introduction
The editors of the recently-created journal
Cilia called their inaugural article “Cilia- the prodigal organelle” [1]. This is perhaps an appropriate denotation for a biologicalappendage found on the surface of a variety of eukaryotic cells, with a fascinatingand still not fully understood relationship between internal structure and biologicalfunction [2]. Cilia not only have the capacity to transport surrounding fluids by meansof periodic movements but they also play important sensory roles [1, 3]. Remarkably,their basic morphology and internal structure (the so-called axoneme), which is thesame as in all eukaryotic flagella, has been conserved throughout evolution [4].The active dynamics of cilia plays a crucial role in a major aspect of biology,namely fluid transport [4–6]. The motion of biological fluids induced by the collec-tive beating of an array of cilia is the third major form of mechanical fluid transportdone by the inner organs of superior animals, ranking only behind blood pumping andperistaltic motion. In humans, the flow produced by the deformation of cilia is involvedin the transport of several biological fluids, including the removal of tracheobronchialmucus in the respiratory track [4–6], the transport of ovulatory mucus and the ovumin the oviduct of the female reproductive track [4, 7], and the motion of epididymalfluid in the e ff erent ducts of the male reproductive track [8]. Failure of the transportfunctionality of cilia can lead to serious illness of the respiratory system [4]. In ad-dition, it is believed that cilia found in specialized brain cells (ependymal cilia) areinvolved in the transport of cerebrospinal fluid at small scales [9]. When cilia are an-chored on a free-swimming cell, such as in many protozoa including the oft-studied Paramecium , their collective deformation and active transport of the surrounding fluidleads to locomotion of the cell, typically at very low Reynolds number [5, 10].The study of fluid transport by cilia, and more generally by eukaryotic flagella,requires an understanding of their beating patterns as well as the physical proper-ties of the surrounding fluid. Many biological fluids in the examples above displaystrong non-Newtonian characteristics, including nonzero relaxation times and shear-dependent viscosities [11]. Past work has addressed theoretically [12–18] and experi-mentally [19–21] the e ff ect of fluid viscoelasticity on transport and locomotion, but theimpact of shear-dependance material functions on the performance of cilia and flagellahas yet to be fully quantified [22].In the present paper we mathematically address the role of shear-dependent vis-cosity on the propulsion performance of a two-dimensional undulating surface whoseboundary conditions model the beating strokes of a cilia array or of a waving flagel-lum. The calculation is an extension of the classical results by Taylor [23, 24] to thenon-Newtonian case. In the next section we describe the physics of the cilia beatingpatterns, distinguish individual vs. collective motion, and introduce the mathematicalmodels that have been proposed to quantify cilia dynamics both in Newtonian and non-Newtonian fluids. 2 igure 1: Schematic representation of the stroke of an individual cilium and the envelope model. Top:e ff ective (solid) and recovery (dashed) stroke of a cilium. The ellipse (thin dashed) represents the strokecycle as modeled by the envelope model; in this figure the direction of the dynamics along the ellipse isclockwise, but it could be reversed (see text). Bottom: Representation of the envelope model covering thecilia layer and the propagation of the metachronal wave. II. Background
II.1. Kinematics of an individual cilium
A variety of beating patterns is displayed by cells employing cilia or flagella, de-pending on the surrounding geometry, the purpose of the beating (locomotion, inges-tion, mixing...) as well as genetic factors and biochemical regulation [2, 4, 5, 25, 26]. Adetailed description of the stroke of cilia has been o ff ered in the past [5, 6, 23, 27] andwe summarize the main features here. In general, the beating pattern of an individualcilium displays a two-stroke e ff ective-recovery motion, as illustrated in Fig. 1 (top).During the e ff ective stroke, the cilium extends into the fluid dragging the maximumvolume of fluid forward. In contrast, the recovery (or backward) stroke is executedby bending the cilium towards itself and the nearest boundary, minimizing thereby thedrag on the surrounding fluid in the opposite direction. It is this asymmetry in the beat-ing pattern that induces a net fluid flow in the direction of the e ff ective stroke. In somecells, the cilium not only bends towards the wall during the recovery stroke but alsotilts with respect to the vertical plane described by the e ff ective stroke, hence display-ing a three-dimensional pattern [6, 28]. Cilia in the respiratory track seem, however, toremain two-dimensional [29]. II.2. Collective cilia motion and metachronal waves
If we now turn our attention to the collective motion of cilia array, the cilia areobserved to beat in an organized manner such that an undulating surface will appear toform on top of the layer of cilia, and to deform in a wave-like fashion. These waves,known as metachronal waves, are due to a small phase lag between neighboring cilia,3 arameter ReferenceCilium length ( L ) 6 µ m [31]Beating frequency ( f = ω/ π ) 20 Hz [9]Distance between cilia ( d ) 0.4 µ m [6]Metachronal wavelength ( λ = π/ k ) 30 µ m [31] Table 1: Beating parameters for metachronal waves of tracheobronchial cilia: cilia length ( L ), beating fre-quency ( f ), distance between cilia ( d ), and metachronal wavelength ( λ ). These numbers are found in thereferences indicated in the last column and are akin to waves created in sports stadiums by waving spectators. Many biologicalstudies have been devoted to this collective e ff ect [26, 30]. In Table 1 we reproducesome typical values of the beating parameters measured experimentally for cilia move-ment in the respiratory track (tracheobronchial cilia). The cilia length and the distancebetween cilia at their base are denoted by L and d , while f = ω/ π is the beat frequencyand λ = π/ k is the metachronal wavelength.Di ff erent types of metachronal waves can be classified according to the relation-ship between their dynamics and that of the e ff ective stroke of the constituting cilia [5].When the propagative direction of the metachronal wave is the same as the directionof the e ff ective stroke, the beat coordination is called symplectic. If instead both direc-tions oppose each other then the coordination is termed antiplectic [23]. Other typesof metachronal waves that requires certain three-dimensional coordination have alsobeen identified [6, 27, 28, 32]. Although there is no yet agreed-upon explanation forthe origin of these di ff erent types of coordination, it seems that patterns close to theantiplectic metachrony (including tracheobronchial cilia) are more widely used thanthe symplectic wave [33] possibly because the first one is characterized by the separa-tion of consecutive cilia during the e ff ective stoke, allowing them to propel more fluidvolume and thus increasing the stoke e ffi ciency [34]. II.3. Modeling cilia arrays
Two di ff erent approaches have been proposed to model the periodic motion of cilia,namely the sub-layer and the envelope models [23]. In sub-layer modeling, one keepsthe basic ingredients of the cilia beating cycle by parametrizing each cilium shape alongits stroke period. This approach allows to compute the forces and bending momentsinduced by each cilium to the surrounding fluid and at the same time obtain the meanflow produced above the cilia layer, although the hydrodynamic description can betedious and usually has to be done numerically [10, 34, 35].In contrast, the envelope modeling approach takes advantage of the formation ofmetachronal waves above the cilia layer. The motion of the cilia array is simplified asan undulating surface that covers the cilia layer (Fig. 1, bottom), ignoring the details ofthe sub-layer dynamics. In this case, the metachronal waves are prescribed by settingappropriate boundary conditions on the material points of the envelope. In the casewhere the boundary conditions are nearly inextensible, the dynamics of the flow pro-duced by the cilia array becomes similar to the undulating motion of a two-dimensionalflagellum [5]. The simplification allowed by the envelope model has been useful incomparing, even quantitatively, the swimming velocities obtained theoretically with4hose observed in water for a number of microorganisms [27, 36]. Furthermore, thismodel, which is originally due to Taylor [24], is amenable to perturbation analysis andallows to incorporate certain non-Newtonian e ff ects in a systematic fashion [12, 15].To address the limit of validity of the envelope model, Brennen [27] derived the setof conditions to be satisfied in order for the envelope model to be approximately valid,which are ν/ ( ω L ) ≫ kd and kd <
1, where ν is the kinematic viscosity of the fluid.With the numbers from Table 1, we see that the envelope model is valid for the caseof the cilia array in the respiratory track when the liquid has a viscosity similar to thatof water (and ν will be actually higher if the fluid is a polymeric solution, hence themodel is expected to always be valid). II.4. Cilia in complex fluids
Because of the importance of mucus in human health, numerous theoretical papershave attempted to elucidate the mechanisms of tracheobronchial mucus transport usingboth modeling approaches outlined above [31, 35, 37, 38]. Mucus is a typical exampleof a non-Newtonian fluid, displaying both elastic behavior and stress relaxation as wellas a shear-dependent viscosity [11].In order to incorporate a variation of the viscosity in the transport problem, a two-fluid model has been proposed consisting of a low-viscosity fluid located at the cilialayer (peri-ciliary fluid) with a second high viscosity fluid laying above the first one[6, 35, 38, 39]. This hypothesis of a two-fluid layer, supported by experimental obser-vations of tissue samples, allows to obtain values of the mucus flow similar to thosemeasured in experiments [6, 28].An alternative modeling approach proposed to use a fluid with a viscosity whosevalue changes linearly with the distance from the cell wall [37]. Such approximation,however, is not a formal mechanical model that rigorously takes into account the rheo-logical properties of the fluid. Later work considered formally the elastic properties ofmucus by introducing the convective Maxwell equation in the theoretical formulation[31]. The values of the mean velocity of the fluid in this case were found to be muchsmaller than the experimental ones, indicating that the shear-thinning properties of mu-cus also need to be considered. Elastic e ff ects were the focus on a number of studies onflagellar locomotion, both theoretically [12–18] and experimentally [19–21], and theconsensus so far seems to be that for small-amplitude motion elastic e ff ects do alwayshinder locomotion whereas it can potentially be enhanced for large amplitude motioncorresponding to order-one Weissenberg, or Deborah, numbers.For mucus transport by cilia, it is in fact believed that it is a combination of shear-thinning, which can reduce the fluid flow resistance at the cilia layer, and elastic e ff ects,which maintain a semi-rigid surface at the upper mucus layer, which allow for the op-timal transport of particles in the respiratory track [6, 28, 37, 39, 40]. A constitutiveequation that takes into account a shear-dependent viscosity was formally used withinthe lubrication approximation in the context of gastropod [41, 42] and flagellar loco-motion [43]. These models consider the locomotion due to a tangential or normal (orboth) deformation of a sheet above a thin fluid layer with shear-dependent viscosity. Inthe case of tangential deformation of the sheet, the average locomotion speed is seen togo down for a shear-thinning fluid whereas in the case where normal deformation arealso included the speed increases. This indicates that the impact of shear-dependent5aterial function on transport depends on the particular beating pattern, something ourmodel also shows. Recent studies on the e ff ect of shear-dependent viscosity on wav-ing propulsion showed that locomotion is una ff ected by shear-thinning for nematodes(experiments, [19]), increases in a two-dimensional model flagellum with an increasein flagellar amplitude for shear-thinning fluids (computations, [22]), while decreases inthe case of a two-dimensional swimming-sheet-like swimmer model in shear-thinningelastic fluids (experiments, [21]).Other models focused on the presence of yield stresses and the heterogeneity inthe surrounding environment [42–44]. The presence of a yield stress in the limit ofsmall amplitude oscillations hinders the mean propulsion or transport velocity whereasthe inclusion of obstacles mimicking an heterogeneous polymeric solution enhancespropulsion. II.5. Outline
In this paper we solve for the envelope model in a generalized Newtonian fluidusing a domain perturbation expansion. The fluid is the Carreau model for general-ized Newtonian fluids where the viscosity is an instantaneous nonlinear function of thelocal shear rate. We solve for the flow velocity and rate of work in powers of the am-plitude of the sheet deformation. In particular, we compute analytically the influenceof non-Newtonian stresses on the fluid velocity induced in the far field: that veloc-ity can be interpreted either as a transport velocity in the context of fluid pumping oras a propulsive speed in the context of locomotion. In section III we present the for-mulation of the envelope model and the dynamic equations that govern the fluid flow.Section IV details the asymptotic calculations followed by section V which is devotedto an analysis of our mathematical results and the illustration of the importance of non-Newtonian e ff ects. Finally in section VI we discuss the impact of our results in thecontext of biological locomotion and transport and conclude the paper. In Appendix Athe applicability of our results to a wider class of generalized Newtonian fluid modelsis examined. III. Formulation
III.1. The envelope model
In the envelope model approach, the cilia tips, attached to a stationary base, displaya combination of normal and tangential periodic motions. Depending on the phase shift φ between these two motions and their relative amplitude ( − π ≤ φ ≤ π ), in general thecilia tips will describe an elliptic shape in the xy -plane over one period of oscillation(Fig. 1). We assume that the cilia tips are connected by a continuous waving sheet, andthe variation of the distance between adjacent cilia is therefore idealized as a stretchingor compression of the surface. The positions of material points ( x m , y m ) on the sheetare given by x m = x + Aa cos ( kx − ω t − φ ) , y m = Ab sin ( kx − ω t ) , (1)where 2 Aa and 2 Ab are the maximum displacements of the material points in the x and y directions respectively. The wave amplitude is thus A (dimensions of length)6hile the relative values of the dimensionless parameters a > b > x direction (Fig. 1) and the material points of thesheet can describe a cyclic motion either in a clockwise direction (symplectic), or acounterclockwise direction (antiplectic).We nondimensionalize lengths by k − and times by ω − to get the nondimensional-ized form of eq. 1 as x m = x + ǫ a cos ( x − t − φ ) , y m = ǫ b sin ( x − t ) , (2)where we have defined ǫ = Ak = π A /λ . Assuming that the cilia array propagatesmetachronal waves with a single characteristic amplitude much smaller than its wave-length, we can use a standard domain perturbation method to solve for the pumpingvelocity in the limit ǫ ≪
1. In the remaining text we will use for convenience z = x − t and ξ = x − t − φ . Using the formula cos ( z − φ ) = cos z cos φ + sin z sin φ , we will alsodefine β = a cos φ , γ = a sin φ , and a = β + γ . III.2. Governing equations of the fluid flow and boundary conditions
Neglecting any inertial e ff ects in the fluid, mechanical equilibrium is simply givenby ∇ · σ = where σ is the total stress tensor. Denoting p the pressure and τ the stresstensor component due to fluid deformation (deviatoric) we have ∇ p = ∇ · τ , (3)which is associated with the incompressibility condition, ∇ · u =
0, where u denotesthe velocity field.In our model, the fluid problem is two-dimensional, and we can therefore reduce thenumber of scalars used in the formulation by introducing a streamfunction, ψ ( x , y , t ).The incompressibility condition is satisfied everywhere in the fluid provided we use u = ∂ψ∂ y , v = − ∂ψ∂ x · (4)Using eq. 4 and the time derivatives of the material position of the undulating sheet, u = ∂ x m /∂ t and v = ∂ y m /∂ t , we can identify the velocity components of the fluid onthe sheet assuming the no-slip boundary condition ∇ ψ | ( x m , y m ) = ǫ b cos ( x − t ) e x + ǫ a sin ( x − t − φ ) e y (5)In the far field, y → ∞ , the flow moves at the unknown steady speed, U e x , whosesign and value will depend on the particular movement of the material points on thesheet and on the fluid rheology. The main goal of this paper is to derive analyticallythe non-Newtonian contributions to U . III.3. Constitutive equation: the Carreau model
As discussed above, many biological fluids exhibit non-Newtonian rheology. InFig. 2 we reproduce the data from several published experiments on the shear-dependance7 −6 −4 −2 −2 −1 shear rate (1/s) v i sc o s i t y ( P a ⋅ s ) Figure 2: Shear-thinning viscosity of mucus reproduced from published data: Steady ( (cid:7) ) and oscillatory ( (cid:4) )measurements of human sputum [45]; oscillatory measurements ( N ) of pig small intestine mucus [46]; steadymeasurements ( • ) of human cervicovaginal mucus [47]; micro-rheological measurements ( (cid:7) ) of human lungmucus [48], micro-rheological measurements ( • ) of human cervical mucus [48]. When the data was reportedin terms of oscillatory moduli, the Cox-Merz rule was applied to estimate the steady values [49]. Thecontinuous line denotes a power-law fitting, η = . γ − . . The dashed lines denote a fitting to the Carreaumodel for the lung mucus ( (cid:7) : η o = . , λ t = . , n = .
08) and the cervical human mucus ( • : η o = . , λ t = . , n = . . ˙ γ .
50 s − ,indicates the typical oscillation frequencies of cilia [5]. of mucus viscosity [45–48]. The plot shows the value of the shear viscosity, η , as afunction of the shear rate, ˙ γ . Measurements are either steady-state or, for oscillatorydata, their steady-state values are estimated using the Cox-Merz rule [49]. Clearly, mu-cus is very strong shear-thinning, the viscosity dropping by three orders of magnitudeon the range of shear rates 10 − . ˙ γ . .The dependance of the viscosity with the shear-rate can be modeled assuming aCarreau generalized Newtonian fluid whose viscosity is given by [49] η − η ∞ η o − η ∞ = h + ( λ t | ˙ γ | ) i n − , (6)where η ∞ is the infinite-shear-rate viscosity, η o the zero-shear-rate viscosity, n the di-mensionless power-law index, and λ t a time constant associated with the inverse of theshear-rate at which the viscosity acquires the zero-shear-rate value. In eq. 6, the mag-nitude of the shear-rate tensor, | ˙ γ | , is equal to ( Π / / , where the second invariant of8he shear rate tensor, Π , is given by Π = X i = X j = ˙ γ i j ˙ γ ji = ∂ u ∂ x ! + ∂ u ∂ y ! + ∂ u ∂ y ∂ v ∂ x + ∂ v ∂ x ! + ∂ v ∂ y ! · (7)Based on the data reproduced in Fig. 2 it is clear that η o ≫ η ∞ , and this is the limitwe will consider here. If we nondimensionalize stresses and shear-rates by η o ω and ω ,respectively, we obtain the nondimensional form of the deviatoric stress, τ = h + (Cu | ˙ γ | ) i N ˙ γ , (8)where we have defined the non-Newtonian index N ≡ ( n − / = ωλ t is the Carreau number. When n = N = ff ects of shear-dependence onthe flow. Any other model, provided that it is well-defined in the zero-shear-rate limit,would either give results mathematically equivalent to ours, or would predict an impacton the flow generated by the sheet arising at a higher-order in the sheet amplitude. IV. Asymptotic solution for small wave amplitudes
Having posed the mathematical problem, we proceed in this section to compute itssolution order by order in the dimensionless wave amplitude, ǫ [12, 23, 24]. Our goal isto calculate the leading-order influence of the shear-dependence of the viscosity on theflow field and, in particular, the flow induced at infinity, U e x . We will also computeits influence on the rate of work done by the sheet in order to deform and the resultingtransport e ffi ciency.We therefore write a regular perturbation expansion of the form { u , ψ, τ , p , ˙ γ , | ˙ γ |} = ǫ { u (1) , ψ (1) , τ (1) , p (1) , ˙ γ (1) , | ˙ γ | (1) } + ǫ { u (2) , ψ (2) , τ (2) , p (2) , ˙ γ (2) , | ˙ γ | (2) } + ... (9)Using the definition of the magnitude of the shear-rate tensor, | ˙ γ | = √ Π /
2, the term | ˙ γ | in eq. 8 can be expanded in terms of ǫ as | ˙ γ | = ǫ Π (1) + ǫ | ˙ γ (1) || ˙ γ (2) | + O ( ǫ ) . (10)Writing the viscosity in eq. 8 as η = " + ǫ Cu Π (1) + ǫ Cu | ˙ γ (1) || ˙ γ (2) | + O ( ǫ ) N , (11)9nd Taylor-expanding eq. 11 up to O ( ǫ ), we finally obtain the constitutive relationshipas ǫτ (1) + ǫ τ (2) + ... = " + ǫ N Cu Π (1) ǫ ˙ γ (1) + ǫ ˙ γ (2) + ... i . (12)Note that, because of the ǫ → − ǫ symmetry, the fluid transport at infinity, quantifiedby U e x , has to scale as an even power of ǫ . Since the perturbation in the viscosity,eq. 11, is order ǫ , this means that the flow up to order ǫ is going to be the same asthe Newtonian one, and thus we expect the non-Newtonian e ff ects to come in only atorder ǫ in U (or higher, see Appendix A). As we see below, this means that we haveto compute the flow field everywhere up to order ǫ . IV.1. Solution of the velocity field at O ( ǫ ) and O ( ǫ )Is is straightforward to get from eq. 12 that τ (1 , = ˙ γ (1 , and thus the solutionsat order one and two are the same as the Newtonian ones. The general procedure toderive these solutions consists in taking the divergence and then the curl of eq. 3 toeliminate the pressure, substituting the constitutive relation for the stress, and using thedefinition of the streamfunction, eq. 4, together with the definition of the shear-rate,˙ γ = ∇ u + ∇ u T , to obtain the biharmonic equation satisfied by the streamfunction ∇ ψ ( p ) = , p = , . (13)The general solution of eq. 13 at order O ( ǫ p ) is already known to be [23, 50] ψ ( p ) = U ( p ) y + ∞ X j = h(cid:16) A ( p ) j + B ( p ) j y (cid:17) sin ( jz ) + (cid:16) C ( p ) j + D ( p ) j y (cid:17) cos ( jz ) i e − jy , (14)which already satisfies the boundary conditions in the far field, u = ∂ψ/∂ y | ( x , ∞ ) = U ( p ) and v = − ∂ψ/∂ x | ( x , ∞ ) =
0. The values of the coe ffi cients in eq. 14 are computed usingthe boundary conditions on the sheet, eq. 5, where ∇ ψ is expanded in terms of ǫ usinga Taylor approximation of the gradients around the material points, x m = x and y m = O ( ǫ ) and O ( ǫ ) are then [23] ψ (1) = ( b + β y + by ) e − y sin z − γ ye − y cos z , (15) ψ (2) = (cid:16) b + b β − a (cid:17) y + y e − y (cid:16) γ − β − β b − b (cid:17) cos 2 z − ye − y ( b γ + βγ ) sin 2 z , (16)where the net contribution to the fluid velocity in the far field appears up to O ( ǫ ) as U (2) = (cid:16) b + ab cos φ − a (cid:17) . (17) IV.2. Solution of the velocity field at O ( ǫ )At order three, eq. 12 leads to the non-Newtonian contribution from the constitutiveequation as τ (3) = N Cu Π (1) ˙ γ (1) + ˙ γ (3) . (18)10he values of the tensor invariant, Π (1) , and the components of the shear rate tensor,˙ γ (1) , can be computed using the first order solution, eq. 15, recalling that ˙ γ = ∇ u + ∇ u T ,and using eq. 7 we obtain u (1) = e − y (cid:8)(cid:2) β − ( β + b ) y (cid:3) sin z + γ ( y −
1) cos z (cid:9) , (19) v (1) = − e − y (cid:8)(cid:2) b + ( β + b ) y (cid:3) cos z + γ y sin z (cid:9) , (20)and therefore ˙ γ (3)11 = e − y { [ β − ( β + b ) y ] cos z − γ ( y −
1) sin z } , (21)˙ γ (3)12 = ˙ γ (3)21 = e − y { [( β + b ) y − β ] sin z + γ (1 − y ) cos z } , (22)˙ γ (3)22 = − ˙ γ (3)11 , (23)which allows us to finally compute Π (1) = e − y h b y + ab ( y − y ) cos φ + a ( y − i . (24)After taking the divergence and then the curl of eq. 18 we obtain inhomogeneousbiharmonic equation for the streamfunction at order three given by ∇ ψ (3) = N Cu e − y nh γ b f ( y ) + b βγ g ( y ) + a γ h ( y ) i cos z + h b i ( y ) − b β f ( y ) − b (3 β + γ ) g ( y ) − a β h ( y ) i sin z o , (25)where the y -dependent functions are f ( y ) = y − y + y − , (26) g ( y ) = y − y + y − , (27) h ( y ) = y − y + y − , (28) i ( y ) = − y + y − y . (29)The homogeneous solution of this equation, ψ (3) h , is given by the general solution ineq. 14. The particular solution, ψ (3) p , can be found using the method of variation ofparameters, leading to ψ (3) p = N Cu e − y h { γ b F ( y ) + b βγ G ( y ) + a γ H ( y ) } cos z + { b I ( y ) − b β F ( y ) − b (3 β + γ ) G ( y ) − a β H ( y ) } sin z i , (30)where the y -dependent functions are F ( y ) = y + y + y , (31) G ( y ) = y − y + y + , (32) H ( y ) = y − y + y − , (33) I ( y ) = − y − y − y − · (34)11n order to determine the unknown constants in the homogeneous solution, we needto evaluate the boundary condition on the sheet at O ( ǫ ). Using a Taylor expansion, wehave ∇ ψ (3) | ( x , = − a cos ξ ∂∂ x ∇ ψ (2) | ( x , − b sin z ∂∂ y ∇ ψ (2) | ( x , − a cos ξ ∂ ∂ x ∇ ψ (1) | ( x , − ab cos ξ sin z ∂∂ x ∂∂ y ∇ ψ (1) | ( x , − b sin z ∂ ∂ y ∇ ψ (1) | ( x , . (35)After substituting in eq. 35 the solutions at order one and two we obtain, after sometedious but straightforward algebra, ∂∂ x ψ (3) | ( x , = (cid:16) − b − b β + b β + b γ (cid:17) cos z + (cid:16) b + b β + b β − b γ (cid:17) cos 3 z − (cid:16) b γ + b βγ (cid:17) sin z + (cid:16) b γ + b βγ (cid:17) sin 3 z , (36)and ∂∂ y ψ (3) | ( x , = (cid:16) − b γ − b βγ + a γ (cid:17) cos z + (cid:16) b γ + b βγ + β γ − γ (cid:17) cos 3 z + (cid:16) b + b β + b β − b γ − βγ − β (cid:17) sin z + (cid:16) − b − b β − b β + b γ + βγ − β (cid:17) sin 3 z . (37)Applying these boundary conditions to the general solution, ψ (3) = ψ (3) h + ψ (3) p , we find12he value of the only nonzero coe ffi cients in eq. 14 which are given by A (3)1 = N Cu − ! b + N Cu + ! b β + N Cu + ! b γ − (cid:16) N Cu a + b (cid:17) β, (38) A (3)3 = b + b β + b β − b γ , (39) B (3)1 = N Cu − ! b + N Cu − ! b β + − N Cu ! b β − + N Cu ! b γ − βγ + N Cu a β − β , (40) B (3)3 = − b − b β − b β + b γ − β + βγ , (41) C (3)1 = − N Cu ! b βγ + N Cu a + b ! γ, (42) C (3)3 = − (cid:16) b γ + b βγ (cid:17) , (43) D (3)1 = − N Cu + ! b γ + N Cu − ! b βγ + b + a − N Cu a ! γ, (44) D (3)3 = β γ − γ + b γ + b βγ, (45)and, as expected, the fluid velocity at infinity at this order is U (3) = . (46)In order to illustrate quantitatively the di ff erence between the Newtonian flow in-duced by the traveling wave and the new non-Newtonian terms, we show in Fig. 3the instantaneous streamlines (dotted lines) and contours (colors and numbers) of iso-values of the vorticity at order O ( ǫ ) for two di ff erent values of the phase, φ , and threedi ff erent values of the non-Newtonian term, N Cu , assuming a = b . Changing thephase φ has the e ff ect of tilting the vortical streamlines formed in the flow. Increas-ing the magnitude of N Cu leads to a significant modification of the flow: large strongvortices are created above the waving sheet whose sign changes with the sign of N Cu .The location and strength of the regions of high vorticity near the sheet are also modi-fied. Note that even though the non-Newtonian term is scaling linearly with N , the totalflow, sum of the Newtonian plus the non-Newtonian component, is not anti-symmetricshear-thinning ( N <
0) vs. shear-thickening ( N > IV.3. Net far-field velocity at O ( ǫ )To finish the calculation we have to compute the velocity induced in the far field, U (4) e x , to quantify the role of non-Newtonian stresses on fluid transport, either pumping13 a) N Cu = φ = N Cu = − φ = N Cu = φ = N Cu = φ = π/ N Cu = − φ = π/ N Cu = φ = π/ ǫ obtained for the case a = b for two values of the phase, φ = π/ N Cu = −
10 (shear-thinning fluid, middle column) and +
10 (shear-thickening, right column). The solid lines denote the positionof the traveling wave at t = or propulsion. This can actually be determined without having to compute the entireflow field. Indeed, just like for the previous order, the streamfunction, ψ (4) , will satisfyan inhomogeneous biharmonic equation. The swimming speed, which contributes toa uniform U (4) y term to the streamfunction, does not have to be computed in the farfield but can instead be evaluated using the boundary conditions. The expansion ofthe x -velocity, u , at this order evaluated at ( x ,
0) is obtained by Taylor-expanding the14oundary conditions on the sheet as ∂∂ y ψ (4) | ( x , = − b sin z ∂ ∂ y ψ (3) | ( x , − a cos ξ ∂∂ x ∂∂ y ψ (3) | ( x , − a cos ξ ∂ ∂ x ∂∂ y ψ (2) | ( x , − ab cos ξ sin z ∂∂ x ∂ ∂ y ψ (2) | ( x , − b sin z ∂ ∂ y ψ (2) | ( x , − a cos ξ ∂ ∂ x ∂∂ y ψ (1) | ( x , − a b cos ξ sin z ∂ ∂ x ∂ ∂ y ψ (1) | ( x , − ab cos ξ sin z ∂∂ x ∂ ∂ y ψ (1) | ( x , − b sin z ∂ ∂ y ψ (1) | ( x , . (47)The first term of the right-hand side of eq. 47 is the only one that provides a nonzeronon-Newtonian term after evaluating the derivatives at y =
0. Computing the meanvelocity over one period of oscillation, we obtain a Newtonian plus a non-Newtoniancontribution of the form12 π Z π − b sin z ∂ ∂ y ψ (3) | ( x , dt = N Cu (cid:16) b β − b β − b γ − a b β (cid:17) + (cid:16) b + b β + b β − b γ − a b β (cid:17) . (48)The contribution to the mean velocity of the eight other terms appearing in the right-hand side of eq. 47 are the same as for the Newtonian solution and the final result canbe found in Blake [32]. Adding all the terms and expressing the mean velocity at orderfour in terms of the original parameters a , b , and φ , we are finally able to write the O ( ǫ ) velocity in the far field as the sum of a Newtonian component, U (4) N , and a newnon-Newtonian term, U (4) NN , as U (4) = U (4) N + U (4) NN , (49)with U (4) N = − b + a b − (cid:16) ab + a b (cid:17) cos φ, (50) U (4) NN = N Cu h − a b + (3 ab − a b ) cos φ − a b cos 2 φ i . (51)The Newtonian component, eq. 50, is Blake’s solution [32]. The new term, eq. 51,which quantifies the leading-order e ff ect of a shear-dependent viscosity on the fluidtransport, is the main result of our paper.Two interesting features can be noted at this point by a simple inspection of eq. 51.First we see that if the sheet does not display normal deformation ( b =
0) or no tan-gential deformation ( a =
0) then U (4) NN =
0; in order for the shear-dependent viscos-ity to a ff ect transport, it is this important that both modes of deformation be present, a , , b ,
0. The second interesting point is that the e ff ect scales linearly with N , andthus with the di ff erence between the power-law index, n , and 1. Consequently, whenthey have an influence on transport, shear-thinning and shear-thickening fluids alwaysact in opposite direction and if one type of fluid hinders transport the other fluid willfacilitate it. 15 V.4. Rate of work
We now turn to the calculation of the rate of work, W , done by the sheet to trans-port the fluid (per unit area of the sheet). Since W scales quadratically with the shearrate, the leading order rate of work appears at O ( ǫ ) and is identical to the one for aNewtonian fluid [23, 32]. Using the nondimensionalization W ∼ η o ω / k we have indimensionless variables h W (2) i = ǫ (cid:16) a + b (cid:17) , (52)where we have used h ... i to denote mean over a period of oscillation.In order to compute the rate of work at higher order we recall that, in the absenceof inertia and without the presence of external forces, the rate of work W done by awaving surface S against a fluid is equal to the value of the energy dissipation over theenclosed volume V W = Z S u · σ · n dS = Z V σ : ∇ u dV . (53)Using that the stress tensor and velocity vector have been expanded in powers of thewave amplitude ǫ , the rate of work at third and fourth order are W (3) = Z V (cid:16) σ (1) : ∇ u (2) + σ (2) : ∇ u (1) (cid:17) dV , (54) W (4) = Z V (cid:16) σ (1) : ∇ u (3) + σ (2) : ∇ u (2) + σ (3) : ∇ u (1) (cid:17) dV . (55)At first and second orders in ǫ , the constitutive equations for the stress tensor arethe Newtonian ones, σ (1 , = − p (1 , I + ˙ γ (1 , . The constitutive equation at third orderis in contrast σ (3) = − p (3) I + N Cu Π (1) ˙ γ (1) + ˙ γ (3) . Substituting these into eqs. 54 and55 we see that the pressure at any order multiplies the divergence of the velocity whichis zero by incompressibility. The mean rates of work per unit area of the sheet at order O ( ǫ ) and O ( ǫ ) are thus given by h W (3) i = π Z π Z ∞ (cid:16) ˙ γ (1) : ˙ γ (2) (cid:17) dydz , (56) h W (4) i = π Z π Z ∞ ˙ γ (1) : ˙ γ (3) +
12 ˙ γ (2) : ˙ γ (2) + N Cu Π (1) ˙ γ (1) : ˙ γ (1) ! dydz . (57)Substituting the values of the shear-rate tensor and its second invariant at the respectiveorders in ǫ , we obtain that h W (3) i =
0, which could have been anticipated by symmetry ǫ → − ǫ . The mean rate of work per unit area at O ( ǫ ) is found to be the sum of aNewtonian term plus a non-Newtonian contribution as h W (4) i = h W (4) i N + h W (4) i NN , (58)where h W (4) i N = −
14 ( a + b ) + a b + ( ab + a b ) cos φ − a b cos 2 φ, (59) h W (4) i NN = N Cu (cid:16) a + b + a b − a b cos φ + a b cos 2 φ (cid:17) . (60)16imilarly to the non-Newtonian contribution to the velocity, eq. 51, the e ff ect of theshear-dependent viscosity on W scales linearly with N , and thus shear-thinning andshear-tickening fluids always contribute in opposite sign to the rate of working of thesheet. In contrast however with the fluid velocity, we see that modes with pure normal( a =
0) or tangential ( b =
0) motion do contribute to the rate of work at this order.
IV.5. Transport e ffi ciency We finally consider the e ffi ciency of transport, E , measured as the ratio betweenthe useful work done in the fluid to create the flow at infinity and the total dissipation,written E ∼ U / h W i in dimensionless variables [10]. Given that we have U = ǫ U (2) + ǫ U (4) and h W i = ǫ h W (2) i + ǫ h W (4) i , it is clear what we will obtain E = ǫ E (2) + ǫ E (4) with each term found by Taylor expansion as E (2) = (cid:16) U (2) (cid:17) h W (2) i , E (4) = (cid:16) U (2) (cid:17) h W (2) i U (4) U (2) − h W (4) ih W (2) i ! . (61)The first term, E (2) , is the same as the Newtonian one whereas the second one, E (4) ,includes a non-Newtonian contribution, given by E (4) NN = (cid:16) U (2) (cid:17) h W (2) i U (4) NN U (2) − h W (4) i NN h W (2) i · (62)Non-Newtonian stresses will lead to an increase in the e ffi ciency of transport if theterm on the right-hand side of eq. 62 is positive, which we will investigate below forthe various kinematics considered. V. Do shear-thinning fluids facilitate locomotion and transport?
With our mathematical results we are now ready to address the impact of non-Newtonian stresses on flow transport. We focus on four aspects. We first show that,within our mathematical framework, shear-thinning always decrease the cost of flowtransport. We then demonstrate that the maximum energetic gain is always obtainedfor antiplectic metachronal waves. When the sheet is viewed as a model for flagellarlocomotion in complex fluids, for example spermatozoa, we then show that, at the orderin which we carried out the calculations, there is no influence of the shear-dependentviscosity on the kinematics of swimming. We finally address the more general caseof viewing the sheet as an envelope for the deformation of cilia array, and address theimpact of non-Newtonian stresses on the magnitude of the flow induced in the far fieldfor a variety of sheet kinematics.
V.1. Shear-thinning fluids always decrease the cost of transport
We start by considering the non-Newtonian contribution to the rate of working bythe sheet, eq. 60. We can write h W (4) i NN = N Cu F ( a , b , φ ) where the function F is F ( a , b , φ ) = (cid:16) a + b + a b − a b cos φ + a b cos 2 φ (cid:17) . (63)17e are going to show that this function is always positive, meaning that for shear-thinning fluids ( N < F ≥ F with respect to the phase is given byd F d φ = a b sin φ ( a − b cos φ ) , (64)while its second derivative isd F d φ = a b ( a cos φ − b cos 2 φ ) . (65)The local extrema of F are reached when d F / d φ = φ = , π and, when a < b , cos φ = a / b . When a ≥ b , only the points φ = , π areextrema, and plugging them into eq. 65 we see that the minimum of F is obtained at φ = F min = [7 a + b + a b + a ( a − b )] / a ≥ b ≥
0, and a ≥ b . In the case where a < b it is straightforward to show that the minimum of F occurs at the phase φ satisfying cos φ = a / b . Plugging that value of the phase intoeq. 63, and using that cos 2 φ = a / b − F min = [11 a + b + a b ] / F ≥ h W (4) i NN , always has the same sign as N . For shear-thinningfluids N < U (4) NN ,is equal to zero. As a consequence, the last term in eq. 62, h W (4) i NN / h W (2) i , will al-ways be negative for shear-thinning fluids, indicating a contribution to an increase inthe e ffi ciency. The opposite is true for shear-thickening fluids which always add to theenergetic cost. This result is reminiscent of earlier work on the cost of transport by gas-tropod showing that shear-thinning fluids are advantageous, although for a completelydi ff erent reason (mechanical work versus mucus production) [41].As a final note, it is worth pointing out that, had the flow kinematics been un-changed at all orders compared to the Newtonian ones, then this result of reductionin W would have in fact been obvious. Indeed, if the shear is the same everywherein the fluid, and if the viscosity decreases, then the total dissipation should also de-crease. However, as detailed in § IV.4, the dissipation rate at order four in ǫ dependson the third-order kinematics in the fluid which are di ff erent from what they wouldhave been in the absence of shear-dependence. The fact that our calculations show thatshear-thinning fluids always decrease the cost of transport is therefore not obvious. V.2. Antiplectic metachronal waves o ff er the largest non-Newtonian energy saving Another general result can be gained from a close inspection of the function F in eq. 63. For given values of both amplitudes a and b , it is clear that F it is maxi-mized when cos φ = − φ =
1, which occurs when φ = π . This is thus the18alue of the phase between normal and tangential motion leading to the maximum non-Newtonian decrease in the energy expanded by the sheet to transport the fluid. For thisparticular value of the phase di ff erence, the kinematics of the material point at x = x m = − ǫ a cos t , y m = − ǫ b sin t which means that material points follow elliptical tra-jectories of semi axes ǫ a and ǫ b in the anticlockwise direction. Viewing the sheet as amodel for the motion of cilia tips this is therefore a wave for which the e ff ective strokeis in the − x direction while the recovery stroke occurs in the + x direction. Given thatthe metachronal waves are assumes to propagate in the positive direction, this meansthat antiplectic metachronal waves o ff er the largest non-Newtonian energy saving in ashear-thinning fluid. V.3. The sheet as a model for flagellar locomotion
We have seen in § IV that the presence of non-Newtonian stresses a ff ects the fluidtransported by the waving sheet only if both oscillatory motions, normal ( b ,
0) andtangential ( a , ff erence however between the waving kinematics of a flagellumand the one we have studied in this paper is that a flagellum is inextensible. Conse-quently, the flagellum has O ( ǫ ) motion only in the normal direction, while the tangentialdeformation is O ( ǫ ) as a result from the requirement of inextensibility (material pointson a flagellum describe a typical figure-8 motion). Higher-order kinematics can alsobe described similarly [23, 24, 51].To address the impact of a shear-dependent viscosity on locomotion we have fol-lowed the same procedure as above and computed the swimming velocity systemati-cally up to O ( ǫ ) using the Carreau model and enforcing inextensibility condition. Weobtain formally that U (4) NN =
0: the shear-dependence of the viscosity does not a ff ectthe locomotion speed of the sheet at this order. Inextensible swimmers, such as thoseemploying flagella, and deforming in a constant, small-amplitude waving motion ofthe type modeled here are thus expected to be hardly a ff ected by the shear-dependenceof the fluid. Since the rate of working decreases for shear-thinning fluids, swimmersare thus more e ffi cient than in Newtonian fluids (see eq. 62). Our theoretical resultis consistent with experiments by Shen and Arratia [19] which have shown that theswimming velocity of a nematode ( C. elegans ) waving freely in a shear-thinning fluid(xanthan gum solution) is the same as that observed in a Newtonian fluid with similarviscosity.
V.4. The sheet as a model for fluid transport and pumping
Returning to the sheet as a model for fluid transport by the tips of cilia arrays,we investigate here the e ff ect of the shear-dependent viscosity on the pumping perfor-mance as a function of the values of a , b , and φ . Recall that the magnitude of thenon-Newtonian contribution, U (4) NN , scales with the absolute value of N Cu while itssign is that of N , which is positive if the fluid is shear-thickening, and negative is thefluid is shear-thinning. Among the many combinations of a , b , and φ that we maychoose, we limit the study to three relevant cases: a = b , b ≫ a , and a ≫ b , and in19 igure 4: Representation of the sheet kinematics for the three di ff erent types of motion considered in thiswork: (1) normal plus tangential motion of the same magnitude, a = b ; (2) motion dominated by the normalmode, b ≫ a ; (3) motion dominated by the tangential mode, a ≫ b . In all cases displayed a value of φ = π/ t = π (dashed solid line). In (3) the positions of thematerial points are shown at t = π/ π , and 3 π/ each case consider four illustrative values of φ , namely, φ = φ = ± π/
2, and φ = π . InFig. 4 we display the three di ff erent kinematics of the sheets in the case where the phaseis φ = π/
4. Below we reproduce the results for the fourth-order swimming speed, U (4) ,and discuss its relation to the leading-order Newtonian contribution at order ǫ . Wealso compute the mean rate of working, h W (4) i . The main results are summarized inTable 2 in the discussion section of the paper. V.4.1. Fluid transport with both normal and tangential motion
In the case where the normal and tangential amplitudes are similar, a = b ≡ c , weobtain (i) U (4) = − c N Cu , for φ =
0; the O ( ǫ ) velocity is U (2) = c ; (66)(ii) U (4) = c − N Cu ! , for φ = ± π/ U (2) =
0; (67)(iii) U (4) = c (cid:16) + N Cu (cid:17) , for φ = π ; U (2) = − c . (68)We see that for all the values of the phase φ , the fluid transport is increased, i.e. thesecond order velocity and the non-Newtonian contribution have the same sign, if the20uid is shear-thinning, N <
0. This results remains true for φ = ± π/ O ( ǫ ) is zero as in this case the non-Newtonian contribution has the samesign as the Newtonian term at order four when N < h W (4) i = c (cid:16) + N Cu (cid:17) , for φ =
0; (69)(ii) h W (4) i = c + N Cu ! , for φ = ± π/
2; (70)(iii) h W (4) i = c N Cu , for φ = π ; (71)and the maximum reduction of the rate of work is given when φ = π , as anticipatedabove. In all cases, the transport e ffi ciency in eq. 62 is increased for two reasons: morefluid is being transported, and transport occurs at a lower cost. V.4.2. Motion dominated by the normal mode
In the case where the sheet kinematics is dominated by the normal deformation, b ≫ a , we have U (2) ≈ b /
2. Keeping all terms in a / b up to the leading-order non-Newtonian contribution, we obtain that the velocity at O ( ǫ ) is given by(i) U (4) ≈ − b + ab N Cu − ! , for φ =
0; (72)(ii) U (4) ≈ − b + a b − N Cu ! , for φ = ± π/
2; (73)(iii) U (4) ≈ − b + ab − N Cu ! , for φ = π ; (74)while at leading-order in a / b we have h W (4) i ≈ b N Cu − ! , (75)for all values of the phase. In the cases where φ = ± π/
2, and π , the velocity inducedby the motion of the sheet is increased if the fluid is shear-thinning. In contrast, forin-phase motion, φ =
0, the non-Newtonian term acts in the direction opposite to theleading-order Newtonian term, U (2) . However, the non-Newtonian contribution to thevelocity, eq. 72, scales as ∼ N ab while the contribution to the energetics, eq. 75,scales as N b , and thus in the limit b ≫ a considered here, eq. 62 leads to a positivevalue for E (4) NN : the motion is still more e ffi cient in a shear-thinning fluid. V.4.3. Motion dominated by the tangential mode
The last case we consider is the one where a ≫ b and the kinematics is dominatedby the tangential deformation. In this case, the second order velocity is U (2) ≈ − a / b / a up to the leading-order non-Newtonian contribution in U (4) a ) |N Cu | = 0 . c ) |N Cu | = 10 ( b ) |N Cu | = 1 ab φφ φ Figure 5: In the wave configuration space a / b vs. φ we use colors to indicate the regions in which themagnitude of the net non-Newtonian flow dominates that of the Newtonian component at order O ( ǫ ),i.e. | U (4) NN | > | U (4) N | (from eqs. 50 and 51) for three illustrative values of the non-Newtonian coe ffi cient, N Cu (0.1, 1 and 10). The cut-o ff value of | U (4) NN | / | U (4) N | =
20 is chosen to account for the zeros of U (4) N . Thecolor scheme, shown on the right, gives the iso-values of | U (4) NN | / | U (4) N | . we obtain (i) U (4) ≈ − ba + N Cu ! , for φ =
0; (76)(ii) U (4) ≈ a b − N Cu ! , for φ = ± π/
2; (77)(iii) U (4) ≈ ba + N Cu ! , for φ = π ; (78)while the rate of work is given by h W (4) i ≈ a N Cu − ! · (79)This time, for a shear-thinning fluid, the non-Newtonian contribution acts in the direc-tion of the leading-order Newtonian flow only for φ = π while a flow in the oppositedirection is created for φ = ± π/
2. However, similarly to the discussion in the pre-vious section, since the non-Newtonian e ff ect on the rate of working is ∼ N a (eq. 79)while the impact on the flow speed is ∼ N ba (eq. 76) and ∼ N a b (eq. 77), thetransport e ffi ciency, eq. 62, is systematically decreased in the shear-thinning case. V.4.4. Numerical approach
With the solution for both the net Newtonian flow speed, U (4) N (eq. 50), and thenon-Newtonian one, U (4) NN (eqs. 51), it is straightforward to numerically compare theirvalues. In Fig. 5 we plot the regions in which the ratio in magnitude between bothterms, | U (4) NN / U (4) N | , is above 1 in the wave configuration space a / b vs. φ , for three val-ues of the non-Newtonian term N Cu (0.1, 1, and 10). The figure indicates therefore22he regions where the non-Newtonian contribution dominates the Newtonian one. Thecolor scheme for the iso-values of the velocity ratio is indicated on the right of thefigure, and the cut-o ff of the ratio to account for the zeros of U (4) N is chosen to be 20.When the flow is almost Newtonian ( N Cu = .
1) the non-Newtonian contribution canalmost everywhere be neglected, except near the region where the Newtonian term ex-actly cancels out. As N Cu increases, the non-Newtonian contribution to the net flowbecomes predominant over almost all configuration space; the only region remainingwhite (indicating a ratio of less than 1) is at the bottom of each panel indicating theprevalence of the Newtonian term for waves dominates by normal modes ( a / b ≪ VI. Discussion
In this paper we propose a mathematical modeling approach to address the roleof shear-dependent viscosity on flagellar locomotion and transport of mucus by ciliaarray. We employ the envelope model originated by Taylor [23, 24] and allow for bothnormal and tangential deformation of a two-dimensional waving sheet. We compute theflow field induced by a small-amplitude deformation of the envelope in a generalizedNewtonian Carreau fluid with power index n up to order 4 in the dimensionless wavingamplitude, ǫ . The net flow induced at infinity can be interpreted either as a net pumpingflow or, in the frame moving with the sheet, as a swimming velocity.At leading order, i.e. at order ǫ , the flow induced by the sheet, and the rate ofworking by the sheet to induce this flow, is the same as the Newtonian solution. Thenon-Newtonian contributions in both the induced net flow and the rate of working ap-pear at the next order, i.e. O ( ǫ ). In both cases, the e ff ect is linear in n − n <
1) and shear-thickening fluids ( n >
1) always induceopposite e ff ects. The leading-order non-Newtonian contribution to the flow created atinfinity is non-zero only if the sheet deforms both in the direction normal and tangen-tial to the wave direction. The non-Newtonian contribution to the rate of working isfound to always be negative in the case of shear-thinning fluids, which are thus seen tosystematically decrease the cost of transport independently of the details of the wavekinematics. The maximum gain in dissipated energy is obtained for antiplectic waveswhere the direction of the propagating wave and that of the e ff ective cilia stroke areopposite.The sheet kinematics can be interpreted as a model for the locomotion of a flag-ellated microorganism. In this case, and making the biologically-realistic assumptionthat the flagellum is inextensible, we find that the swimmer is more e ffi cient in a shear-thinning fluid but that the non-Newtonian contribution to the swimming speed of thesheet is exactly equal to zero. This result is consistent with recent experiments [19]showing that the swimming speed of a nematode waving freely in a shear-thinning fluidis same as that observed in a Newtonian fluid. In contrast, computations with swim-mers deforming their flagella with increasing amplitude showed that the locomotionspeed was enhanced by shear thinning. The e ff ect was attributed to a viscosity gradientinduced by the increased beating flagellar amplitude, a result which could potentiallybe tackled theoretically with an approach similar to ours adapted to their kinematics[22]. In contrast, recent experiments with a two-dimensional swimmer model showed23 ormal + tangential Normal mode Tangential mode a = b b ≫ a a ≫ b φ = + velocity: − velocity: − e ffi ciency: + e ffi ciency: + e ffi ciency: + φ = − π/ , π/ + velocity: + velocity: − e ffi ciency: + e ffi ciency: + e ffi ciency: + φ = π velocity: + velocity: + velocity: + e ffi ciency: + e ffi ciency: + e ffi ciency: + Table 2: E ff ect of a shear-thinning viscosity on the mean transport velocity and e ffi ciency for three di ff erentwave kinematics: normal plus tangential deformation, motion dominated by the normal mode, and motiondominated by the tangential mode. Four di ff erent values of the phase between normal and tangential mo-tion are considered. The sign + (resp. − ) indicates a situation where a shear-thinning fluid improves on(resp. hinders) the velocity or e ffi ciency; all signs would be reversed in a shear-thickening fluid. a decrease of the swimming speed in a shear-thinning elastic fluid [21], a result not ex-plained by our theoretical approach but perhaps due to the importance of higher-orderterms for large-amplitude swimming or due to non-negligible linear elastic e ff ects.Viewing the sheet as a model for the transport of non-Newtonian fluids by ciliaarrays, we address three di ff erent types of kinematics to understand the impact of shear-thinning fluids on the flow transport: kinematics dominated by normal deformation,those dominated by tangential motion, and kinematics where normal and tangentialbeating were of the same order of magnitude. In all three cases we consider fourdi ff erent values of the phase between the normal and tangential waving motion, andthe results are summarized in Table 2 where we used “ + ” to denote instances when theshear-thinning fluid improves on the performance (flow velocity or transport e ffi ciency)and “ − ” otherwise. In all cases, a shear-thinning fluid renders the flow transport moree ffi cient, and in eight out of the twelve cases it also increases the transport speed –in particular for the biologically relevant case where the magnitudes of the normaland tangential deformation are of the same order. A further numerical investigation in § V.4.4 also identifies the regions in which the non-Newtonian term is stronger than theNewtonian one.In this paper we have assumed that the fluid has no memory but instead its viscosityis an instantaneous function of the shear rate – a class of models known as generalizedNewtonian fluids. Under that limit, we saw that the shear-thinning property of mucusdoes facilitate its transport by cilia arrays. The constitutive model we used in thiswork then serves as an alternative model to the previously-studied two-fluid modeland allows us to treat the mucus layer as a continuous fluid. We still need, however, toestimate the contribution of the shear-thinning e ff ects. Indeed, they do not play any roleat leading order but only contribute at O ( ǫ ). Elastic stresses, in contrast, do contribute24t order O ( ǫ ) and hinder the fluid transport by the prefactor [12, 19].1 + De η s /η + De , (80)where the Deborah number, De, is defined as De = λ r ω , with λ r being the relaxationtime of the fluid, and η s /η is the ratio of the solvent viscosity with the viscosity ofthe polymeric solution (typically η s /η ≪ λ r ≈
30 s [52],and thus, using the typical frequency in Table 1 for tracheobronchial cilia, we haveDe ≈ × ≫
1. For a solvent with the same viscosity as water, we have η s /η ≈ − at zero shear rate (see Fig. 2) and thus De η s /η ≫
1. Elastic e ff ects lead therefore to adecrease of the transport speed by a factor of η s /η .To compare simply the elastic and shear-thinning contributions, consider a vis-coelastic fluid being deformed by a low-amplitude undulating surface with normaland tangential motion of the same order of magnitude. The magnitude of the shear-dependent contribution to the transport velocity scales as ∼ ǫ | N | Cu while the elas-tic contribution scales as ∼ ǫ η s /η . With | N | ≈ . λ t ≈ ω = π ×
20 Hz (Table 1), the ratio of shear-dependent to elastic e ff ects varies as ∼ ǫ | N | Cu / ( η s /η ) ∼ ǫ . Despite the di ff erent scaling with the motion amplitude ǫ , and although our results are only strictly valid in the mathematical limit ǫ → U (4) NN , tothe Newtonian speed at leading order, U (2) , and we find U (4) NN / U (2) ∼ ǫ | N | Cu /ǫ ∼ ǫ , another potentially very large ratio . It is therefore clear that the rheologicalproperties of the surrounding fluid will play important roles in the fluid flow producedby cilia arrays [6, 28, 40].Beyond the assumption on the fluid rheology, the model proposed in this paper hasbeen carried out under a number of mathematical assumptions. The calculation wastwo-dimensional and as such does not address the flow in the direction perpendicularto the wave propagation. For the sheet as a model for flagellar locomotion, the finite-length of the flagellum is expected to play an important role in non-linear fluids. Andfor cilia, the discrete nature of the cilia in a dense array will also lead to nontrivial flowin the sublayer with contributions to both transport and energetics still unquantified. Itis hoped that these limitations do not prevent the current results to still be relevant toboth locomotion and transport studies, for example mucus transport in the respiratorytrack, but instead suggest potentially interesting directions with future work. Acknowledgements
R. V´elez acknowledges CONACyT-M´exico and UC-MEXUS for their financialsupport during his postdoctoral stay at UCSD. We also thank the support of the NSF We emphasize again that the calculations presented in this paper are only valid in the mathematical limit ǫ →
0. The order-of-magnitude estimates presented above are only formally true when ǫ is small and areused to point out how easy it is for the shear-dependent terms to become predominant. Appendix A. Application to other generalized Newtonian fluids
In this paper we have used a specific empirical model (Carreau fluid) to mathemat-ically describe the relationship between viscosity and shear rate. How applicable areour results to fluids described by other rheological laws, η ( ˙ γ )? In the asymptotic limitassumed this work, we need a fluid model able to describe the behavior of flows in thelimit of low shear rates. We first require the viscosity, η , to be a smooth function ofthe shear rate, ˙ γ , leading a nonzero value of η in the limit ˙ γ →
0, denoted η (0). Thisrequirement rules out the use of the power-law model, valid only in the limit of highshear rates [49]. Due to the ǫ → − ǫ symmetry, the shear-dependent viscosity needsto be an even function of the shear rate. The only even function at order one in ǫ isthe absolute value function, which is not smooth (di ff erentiable) at ˙ γ =
0. Hence, theviscosity necessarily has to be a function of the shear rate square. As described in thetext, we have | ˙ γ | = ǫ Π (1) + O ( ǫ ) . (A.1)In the limit of small wave amplitudes, the viscosity then takes the general form η ( | ˙ γ | ) = η " ǫ Π (1) + O ( ǫ ) = η (0) + ǫ Π (1) d η d ˙ γ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (A.2)using a Taylor expansion near zero. The shear-dependent viscosity will thus give ane ff ect at order O ( ǫ ) at best, and thus at O ( ǫ ) in the flow field. If the derivative ineq. (A.2) is zero then the e ff ect will be at higher order. If instead the derivative isnot defined, then the model is not appropriate to describe flow behavior at small shearrates. In the case of the Carreau model considered above, the derivative d η/ d ˙ γ | takesa finite value, so we indeed obtain a nonzero e ff ect at second order. There are othermodels that could potentially be used as well to fit rheological data. One of such modelis that of a Cross fluid, given in a dimensionless form by η − = + h ( Cu | ˙ γ | ) i − N . (A.3)In that case, it is first straightforward to show that only N ≤ N = − d η/ d ˙ γ | . When N = − N > − N < − ff ects on the fluidflow would appear at a higher order. A similar analysis can be carried out for the Ellisfluid, given by, η − = + h | τ | /τ / i N (A.4)with similar results ( | τ | is the magnitude of the stress tensor and τ (1 / is the shear atwhich η = / ǫ at least). References [1] P. Beales, P.K. Jackson, Cilia-the prodigal organelle, Cilia , 1 (2012).[2] C.B. Lindemann, K.A. Lesich, Flagellar and ciliary beating: the proven and thepossible, J. Cell Sci. , 519 (2010).[3] R.A. Bloodgood, Sensory reception is an attribute of both primary cilia and motilecilia, J. Cell Sci. , 505 (2010).[4] I. Iba˜nez-Tallon, N. Heintz, H. Omran, To beat or not to beat: roles of cilia indevelopment and disease, Hum. Mol. Genet. , R27 (2003).[5] C. Brennen, H. Winet, Fluid Mechanics of Propulsion by Cilia and Flagella, Ann.Rev. Fluid Mech. , 339 (1977).[6] M.A. Sleigh, J.R. Blake, N. Liron, The Propulsion of Mucus by Cilia, Am. Rev.Respir. Dis. , 726 (1988).[7] L.J. Fauci, R. Dillon, Biofluidmechanics of reproduction, Ann. Rev. Fluid Mech. , 371 (2006).[8] T.J. Lardner, W.J. Shack, Cilia transport, Bull. Math. Biophys. , 325 (1972).[9] C. O’Callaghan, K. Sikand, M.A. Chilvers, Analysis of ependymal ciliary beatpattern and beat frquency using high speed imaging: comparison with the photo-multiplier and photodiode methods, Cilia , 8 (2012).[10] E. Lauga, T.R. Powers, The hydrodynamics of swimming microorganisms, Rep.Prog. Phys. , 096601 (2009).[11] S.K. Lai, Y-Y Wang, D. Wirtz, J. Hanes, Micro- and macrorheology of mucus,Adv. Drug Deliv. Rev. , 86 (2009).[12] E. Lauga, Propulsion in a viscoelastic fluid, Phys. Fluids , 083104 (2007).[13] H.C. Fu, C.W. Wolgemuth, T.R. Powers, Swimming speeds of filaments in non-linearly viscoelastic fluids, Phys. Fluids , 033102 (2009).[14] H. C. Fu, T. R. Powers, and H. C. Wolgemuth. Theory of swimming filaments inviscoelastic media. Phys. Rev. Lett. , 99:258101–258105, (2007).[15] E. Lauga. Life at high Deborah number.
Europhys. Lett. , 86:64001, (2009).[16] J. Teran, L. Fauci, and M. Shelley. Viscoelastic fluid response can increase thespeed and e ffi ciency of a free swimmer. Phys. Rev. Lett. , 104:038101, (2010).2717] L. Zhu, M. Do-Quang, E. Lauga, and L. Brandt, Locomotion by tangential defor-mation in a polymeric fluid, Phys. Rev. E, , 011901 (2011)[18] L. Zhu, E. Lauga, and L. Brandt, Self-propulsion in viscoelastic fluids: pushersvs. pullers, Phys. Fluids, , 051902 (2012)[19] X. N. Shen, P.E. Arratia, Undulating swimming in viscoelastic fluids, Phys. Rev.Lett. , 208101 (2011).[20] B. Liu, T. R. Powers, and K. S. Breuer. Force-free swimming of a model helicalflagellum in viscoelastic fluids. Proc. Natl. Acad. Sci. USA , 108:19516–19520,(2011).[21] M. Dasgupta, B. Liu, H. C. Fu, M. Berhanu, K. S. Breuer, T. R. Powers andA. Kudrolli, Speed of a swimming sheet in Newtonian and viscoelastic fluids,Phys. Rev. E, , 013015 (2013)[22] T. D. Montenegro-Johnson, A. A. Smith, D. J. Smith, D. Loghin, and J. R. BlakeModelling the fluid mechanics of cilia and flagella in reproduction and develop-ment. Eur. Phys. J. E , 35: 111, (2012).[23] S. Childress,
Mechanics of Swimming and Flying (Cambridge University Press,Cambridge, England, 1981).[24] G.I. Taylor, Analysis of the swimming of microscopic organisms, Proc. R. Soc.London , 447 (1951).[25] P. Verdugo, Introduction: Mucociliary Function in Mammalian Epithelia, CellMotil. Suppl. , 1 (1982).[26] K. Takahashi, C. Shingyoji, Control of Eukaryotic Flagella and Cilia, Zool. Sci. , 1393 (2002).[27] C. Brennen, An oscillating-boundary-layer theory for ciliary propulsion, J. FluidMech. , 799 (1974).[28] M.J. Sanderson, M.A. Sleigh, Ciliary activity of cultured rabbit tracheal epithe-lium: beat pattern and metachrony, J. Cell Sci.
331 (1981).[29] M.A. Chilvers, C. O’Callaghan, Analysis of ciliary beat pattern and beat fre-quency using digital high speed imaging: comparison with the photomultiplierand photodiode methods, Thorax , 314 (2000).[30] A. Murakami, K. Takahashi, Correlation of electrical and mechanical responsesin nervous control of cilia, Nature , 48 (1975).[31] S.M. Ross, S. Corrsin, Results of an analytical model of mucociliary pumping, J.Appl. Physiol. , 333 (1974).[32] J.R. Blake, Infinite models for ciliary propulsion, J. Fluid Mech. , 209 (1971).2833] E.W. Knight-Jones, Relations between metachornism and the direction of ciliarybeat in metazoa. Q. J. Microsc. Sci. , 503 (1954).[34] J.R. Blake, A model for the micro-structure in ciliated organisms, J. Fluid Mech. , 1 (1972).[35] J.R. Blake, On the movement of mucus in the lungs, J. Biomech. , 179 (1975).[36] J.R. Blake, A spherical envelope approach to ciliary propulsion, J. Fluid Mech. , 199 (1971).[37] C. Barton, S. Raynor, Analytical investigation of cilia-induced mucous flow, Bull.Math. Biophys. , 419 (1967).[38] G.R. Fulford, J.R. Blake, Muco-ciliary transport in the lung, J. Theor. Biol. ,381 (1986).[39] J.G. Widdicombe, Airway Surface Liquid: Concepts and Measurements, in: Air-way Mucus: Basic Mechanisms and Clinical Perspectives , Ed. D.F. Rogers, M.I.Lethem, Birkh¨auser Verlag, Germany (1997).[40] J. Sad´e, N. Eliezer, A. Silberberg, A.C. Nevo, The role of mucus in transport bycilia, Am. Rev. Respir. Dis. , 48 (1970).[41] E. Lauga, A.E. Hosoi, Tuning gastropod locomotion: Modeling the influence ofmucus rheology on the cost of crawling, Phys. Fluids , 113102 (2006).[42] B. Chan, N.J. Balmforth, A.E. Hosoi, Building a better snail: Lubrication andadhesive locomotion, Phys. Fluids , 113101 (2005).[43] N.J. Balmforth, D. Coombs, S. Pachmann, Microelastohydrodynamics of Swim-ming Organisms Near Solid Boundaries in Complex Fluids, Q. J. MechanicsAppl. Math. , 267 (2010).[44] A.M. Leshansky, Enhanced low-Reynolds-number propulsion in heterogeneousviscous environments, Phys. Rev. E , 051911 (2009).[45] M. Dawson, D. Wirtz, J. Hanes, Enhanced Viscoelasticity of Human Cystic Fi-brotic Sputum Correlates with Increasing Microheterogeneity in Particle Trans-port, J. Biol. Chem. , 50393 (2003).[46] L.A. Sellers, A. Allen, E.R. Morris, S.B. Ross-Murphy, The rheology of pig smallintestinal and colonic mucus: weakening of gel structure by non-mucin compo-nents, Biochim. Biophys. Acta , 174 (1991).[47] S. K. Lai, D. E. O’Hanlon, S. Harrold, S. T. Man, Y. Y. Wang, R. Cone, J. Hanes,Rapid transport of large polymeric nanoparticles in fresh undiluted human mucus,Proc. Natl. Acad. Sci. , 1482 (2007).[48] S.H. Hwang, M. Litt, W.C. Forsman, Rheological properties of mucus, Rheol.Acta , 438 (1969). 2949] R. B. Bird, R.C. Armstrong, O. Hassager, Dynamics of Polymeric Liquids, Vol. 1 (John Wiley & Sons, USA, 1987).[50] L. G. Leal,
Advances Transport Phenomena (Cambridge University Press, Cam-bridge, England, 2007).[51] M. Sauzade, G. J. Elfring, and E. Lauga Taylor’s swimming sheet: Analysis andimprovement of the perturbation series
Physica D. , 240:1567, (2011).[52] A. Gilboa, A. Silberberg, In-situ rheological characterization of epithelial mucus,Biorheology13