The mechanism of propulsion of a model microswimmer in a viscoelastic fluid next to a solid boundary
Thomas R. Ives Alexander Morozov
physics.flu-dyn cond-mat.soft





SUPA, School of Physics and Astronomy, The University of Edinburgh,
    James Clerk Maxwell Building, Peter Guthrie Tait Road, Edinburgh, EH9 3FD,
    United Kingdom


alexander.morozov@ph.ed.ac.uk
SUPA, School of Physics and Astronomy, The University of Edinburgh,
    James Clerk Maxwell Building, Peter Guthrie Tait Road, Edinburgh, EH9 3FD,
    United Kingdom


In this paper we study swimming of a model organism, the so-called Taylor's
swimming sheet, in a viscoelastic fluid close to a solid boundary. This
situation comprises natural habitats of many swimming microorganisms, and while
previous investigations have considered the effects of both swimming next to a
boundary and swimming in a viscoelastic fluid, seldom have both effects been
considered simultaneously. We re-visit the small wave amplitude result obtained
by Elfring and Lauga (Gwynn J. Elfring and Eric Lauga, in Saverio E. Spagnolie,
editor, Complex Fluids in Biological Systems, Springer New York, New York, NY,
2015), and give a mechanistic explanation to the decoupling of the effects of
viscoelasticity, which tend to slow the sheet, and the presence of the boundary,
which tends to speed up the sheet.  We also develop a numerical spectral method
capable of finding the swimming speed of a waving sheet with an arbitrary
amplitude and waveform.  We use it to show that the decoupling mentioned above
does not hold at finite wave amplitudes and that for some parameters the
presence of a boundary can cause the viscoelastic effects to increase the
swimming speed of microorganisms.


The mechanism of propulsion of a model microswimmer in a viscoelastic
    fluid next to a solid boundary
    Alexander Morozov
    December 30, 2023
========================================================================================================





§ INTRODUCTION


Many microorganisms are able to propel themselves through fluid environments by
deforming their bodies. The small size of these organisms, ranging from a few
micrometers in the case of most bacteria, to tens or hundreds of micrometers in
the case of eukaryotes, and their relatively small propulsion speeds dictate
that their swimming typically occurs in the low-Reynolds-number regime, and that
the fluid flow around them obeys the Stokes equation <cit.>. As was
pointed out by Purcell <cit.>, this poses a severe restriction on
how microorganisms move since they have to break the intrinsic
time-reversibility of the Stokes equation – a result commonly known as
Purcell's scallop theorem. In order to propel, microorganisms have to either
deform their bodies or move parts of their bodies in a non-time-reversible
fashion, and a vast number of studies considered how various modes of propulsion
work and what the resulting properties of microorganisms' motion are (see
<cit.> and references therein).

Arguably the most influential model of swimming at low Reynolds numbers is
Taylor's swimming sheet model <cit.> that guided the later studies
of microorganism propulsion.  It comprises an infinite inextensible two
dimensional sheet, periodic in space, that can change its shape by propagating a
wave with a speed c along its waveform. From the point of view of an external
observer, its shape traces out the curve y_s in the xy-plane, given by

    y_s(x, t) = f(x-(c-U)t),

where the wave is travelling in the positive x-direction, and the organism
swims at a speed U along the same axis. The waveform f(x) is a periodic
function with the period 2π/k, where k is the associated wavenumber. For
sheets with sinusoidal waveforms, f(x) = bsin(kx), and small amplitudes, bk
≪ 1, Taylor demonstrated that the sheet moves with the speed U_BN =
c b^2k^2/2 in the direction opposite to the direction of wave propagation
<cit.>.

This result was extended by Katz <cit.> who considered a Taylor sheet
swimming next to a solid boundary and showed that to lowest order in bk its
swimming speed is given by

    U_N = cb^2k^2/2(
            sinh^2(hk)+h^2k^2/sinh^2(hk)-h^2k^2),

where h is the distance between the middle line of the organism and the
boundary. Notably, this speed is larger than U_BN, the swimming speed in the
bulk, for any finite value of h, although this conclusion relies on the
assumption that the organism keeps the same kinematics in the bulk and next to a
wall.  The works by Taylor <cit.> and Katz <cit.> were
instrumental in guiding later studies of low-Reynolds propulsion of various
model swimmers, both in the bulk
<cit.>
and close to surfaces
<cit.>.

Another important extension of the Taylor's result was derived by Lauga
<cit.> who studied a waving sheet swimming in the bulk of a
viscoelastic fluid. Lauga showed <cit.> that for a range of
constitutive models, the small-amplitude swimming speed is given by

    U_B = cb^2k^2/2(1 + β^2/1+^2).

Here, β = η_s/(η_s + η_p), with η_s and η_p being the
viscosity of the solvent and the polymer components, respectively, and =
λ c k is the Deborah number of the problem, where λ is the
longest relaxation time of the fluid. For the fixed kinematics of the organism,
this result suggests that viscoelasticity reduces the propulsion speed of a
small-amplitude sheet compared to its Newtonian value, reaching for large
Deborah numbers the limit β U_BN. These conclusions
were extended to other swimmers
<cit.>
or fluids with different rheological properties
<cit.>,
and were used as a motivation for experimental studies
<cit.>.

While the previous studies provide understanding of how individual effects
influence microswimming (with a notable exception of
<cit.>), the actual ecological situation of many
microorganisms often comprises a combination of these effects. Examples range
from sperm moving in a mucus along the cervix wall <cit.> to
bacterial pathogens invading biofilms of different bacterial species
<cit.>. The simplest model to study such systems should include both
viscoelasticity of the suspending fluid and the presence of a solid boundary,
i.e. be a combination of the effects discussed above. The first step in this
direction was taken by Elfring and Lauga <cit.> who calculated the
swimming speed of a small-amplitude Taylor sheet swimming next to a boundary in
an Oldroyd-B fluid. Surprisingly, the swimming speed they obtain,

    U = cb^2k^2/2(1 + β^2/1+^2)
        (sinh^2(hk)+h^2k^2/sinh^2(hk)-h^2k^2).

is simply a combination of the swimming speeds U_BN, U_N and U_B,
i.e. the effects of viscoelasticity and
the boundary factorise. While <ref>, <ref>, and
<ref> are often cited, there is currently no simple
understanding of the corresponding effects.

The purpose of this work is to provide a mechanistic explanation of the
interplay between viscoelasticity of the fluid and the presence of a solid wall.
Our paper is organised as follows. In <ref> we consider
a small-amplitude Taylor sheet model swimming next to a boundary in an Oldroyd-B
fluid. We re-derive the result obtained by Elfring and Lauga <cit.>,
<ref>, and obtain explicit expressions for the velocity and
stress fields around the swimmer which we will use to develop a
    small-amplitude physical mechanism in <ref>.  In
<ref> we develop a numerical method based on a spectral
representation of hydrodynamic fields to calculate the swimming speed of a
Taylor sheet of any amplitude and shape, and apply it to the situation discussed
above. Finally, in <ref>, we use the velocity and stress
fields calculated in the previous Sections to explain the origin of
<ref>, and <ref>, and discuss the emerging
mechanism of propulsion.




§ SMALL-AMPLITUDE SWIMMING: ANALYTICAL SOLUTION





In this Section we consider a Taylor's waving sheet swimming in a viscoelastic
fluid next to a boundary in the small-amplitude limit, see
<ref>.  In the low-Reynolds-number limit, the flow of the fluid
around the organism is governed by the Stokes equation,
∇·Σ=0, where Σ is the total stress in
the fluid given by

    Σ = -p1 + 2η_s𝐃 + τ.

Here, p is the pressure, 1 is the identity matrix, 𝐃 =
(∇ u + ∇ u^T)/2 is the symmetric strain rate tensor, u is the velocity of the fluid, η_s is the solvent viscosity, and
(…)^T denotes the transpose. The polymeric contribution to the total
stress, τ, arises due to the polymers being stretched and oriented by
local velocity gradients. Here we use one of the simplest viscoelastic
constitutive equations, the Oldroyd-B model <cit.>,
that develops large normal stresses in shear flows that are responsible for many
non-trivial effects exhibited by viscoelastic fluids <cit.>;
but the model does not have any shear thinning effects, i.e. its material
properties are independent of local velocity gradients. Combined together, the
governing equations are given by

    
    
    -∇ p + 2η_s ∇·𝐃 +
            ∇·τ = 0, 
    ∇·u  = 0, 
    τ + λτ = 2 η_p 𝐃,


where η_p is the polymeric contribution to the fluid's viscosity, λ
is the longest relaxation time of the solution and we assumed the fluid to be
incompressible. The upper-convected Maxwell derivative is given by

    τ = ∂_t τ + u·∇τ - ∇u^T·τ - τ·∇u.

The boundary conditions are given by the no-slip boundary conditions at the sheet
and the wall:

    
    
    .u|_y=y_s = u_s, 
    .u|_y=h = u_w,


where u_s and u_w are the velocity of the material points of the
sheet and the wall, respectively.

To address the situation sketched in <ref>, we solve a slightly
more general problem of the sheet in a channel with walls both above and below
it, placed at distances h_+ and h_- from the centreline, respectively.  We
then return to the original single-wall problem by taking h_+→ h and
h_-→∞.

We start by simultaneously introducing dimensionless variables and removing any
explicit time dependence with the help of the following transformation to
starred quantities

    x^* = k(x - ct), y^* = k y, h_±^* = h_± k, 
    
        U^* = U/c, u^* = 1 + u/c, v^* = v/c, 
    
        p^* = p/η c k, τ^* = τ/η_p c k,
         Σ^* = Σ/η c k.

In these coordinates the velocity of the walls is u_w^*= -_x and the
shape of the sheet is approximately fixed in time, such that y^*_s(x^*) =
ϵsin(x^*) + 𝒪(ϵ^3) where ϵ = bk is the
dimensionless wave amplitude.  Our goal then is to find the steady velocity
field surrounding the sheet and, from this, to calculate the sheet's swimming
speed.  From now on we will drop the ^*s.

We expand the fields p, u and τ into Taylor series about
ϵ = 0.  For example, the pressure field is given by

    p = ∑_n=0^∞ p^(n)ϵ^n.

where p^(n) is the “nth-order” contribution to the pressure field.

The velocity field is the solution to <ref>, subject to the
no-slip boundary conditions of <ref>. This solution is
periodic in x, reflecting the symmetry of the underlying problem.  As we are
only interested in the lowest order of the small-ϵ expansion of the
swimming speed, we only need the boundary conditions to the lowest order in
ϵ.  From Taylor's original paper <cit.>, we can show that,
to the lowest order, the velocity of the material points of the sheet,
u_s, are given in our coordinates, where we are co-moving with the
    wave, by


    u_s = -1 + 𝒪(ϵ^2), 
    
        v_s = -ϵcos(x)+ 𝒪(ϵ^2).



In our coordinate system, the swimming speed of the sheet can be found by
averaging the velocity field along the length of the sheet.  Thus, up to the
second order in ϵ, the swimming speed of the sheet is given by

    U = ⟨.u^(1)|_y=y_s⟩ϵ +
        ⟨.u^(2)|_y=0⟩ϵ^2 +
        𝒪(ϵ^3).

To the lowest order in ϵ, averages over the line
y=y_s=ϵsin(x) are equal to averages over the line y=0, thus we take
the ⟨⟩ above as simple x-averages.

We find the first and second order velocity fields by substituting the Taylor
expansion of each of the fields into <ref> and consider each
power of ϵ separately.  To the zeroth order, this procedure yields the
following set of equations

    -∇ p^(0) + β∇^2u^(0) +
         (1 - β) ∇·τ^(0) = 0, 
    ∇·u^(0) = 0, 
    τ^(0) = 𝐃^(0), 
    .u^(0)|_y=y_s =
        .u^(0)|_y=h_+ =
        .u^(0)|_y=-h_- = -_x,

which has the trivial solution p^(0) = 0, τ^(0) = 0, and
u^(0) = -_x.  Note that the zeroth order velocity field does not
contribute to the swimming speed of the sheet as the latter is given by the
difference between the average velocities of the fluid at the sheet and at the
wall, which vanishes at the zeroth order.

The first order velocity field is in fact the same for an Oldroyd-B fluid as for
a Newtonian one <cit.>.  To demonstrate this, we consider the first
order equations:

    
    
    -∇ p^(1) + ∇·(2β𝐃^(1) +
            (1 - β) τ^(1)) = 0,
            
    ∇·u^(1) = 0, 
    
            (1 -  ∂_x)τ^(1) =
            𝐃^(1), 
    .u^(1)|_y=y_s = -cos(x)_y, .u^(1)|_y=h_+ =
            .u^(1)|_y=-h_- = 0.


Here we have used the previous solution, u^(0) = -1, in
<ref>, and we have re-arranged <ref> using
∇^2u^(1) = 2∇·𝐃^(1). Let ℒ be the
linear operator defined by

    ℒ(a) = (1- ∂_x)∇·ℰ·a,

where ℰ = _x_y - _y_x.  Applying ℒ to
<ref>, we obtain

    ∇·ℰ·∇·(
            (1- ∂_x)(2β𝐃^(1) +
            (1 - β) τ^(1))
        ) =2∇·ℰ·∇·(
            1 - β ∂_x
        ) 𝐃^(1)= 0,

where we have used the commutativity of differential operators and
<ref> to remove τ^(1).  This equation is satisfied
either by a 𝐃^(1) for which
∇·ℰ·∇·𝐃^(1) = 0, or by a
𝐃^(1) for which  (1 - β ∂_x)𝐃^(1) = 0.
The first of these conditions is satisfied by the Newtonian solution, while the
second has no non-trivial solutions which are periodic in x.  Moreover, since
the boundary conditions are the same as in the Newtonian case, we conclude that
the first order velocity field in an Oldroyd-B fluid is that same as its
Newtonian counterpart.

In his original analysis of a Taylor sheet swimming next to a wall, Katz
<cit.> showed that the first order velocity field, u^(1)_±
in the Newtonian fluid above (+) and below (-) the sheet is given by

    
    
    u^(1)_± = (1 + A_± - B_± y)sin(x)sinh(y) +
            A_± y sin(x)cosh(y), 
    
            v^(1)_± = -(A_± y + B_±)cos(x)sinh(y) -
            (1 - B_± y)cos(x)cosh(y),


where

    A_± = sinh^2(h_±)/sinh^2(h_±) - h_±^2, 
            B_± = ±sinh(h_±)cosh(h_±) ± h_±/sinh^2(h_±) - h_±^2.

The contribution of this field to the swimming speed is given by the first term of
<ref>, which reads

    ⟨.u^(1)|_y=y_s⟩ϵ =
        ϵ^2 (A + 1/2) + 𝒪(ϵ^3).

Here we have dropped the explicit ± notation, but we note that this
contribution to the swimming speed is different for each region of the fluid.
Below, we ensure that the total swimming speed is the same regardless of whether
we use the fluid velocity above or below the sheet; but first we calculate the
second term in <ref>.

The second order set of governing equations is given by

    
    
    -∇ p^(2) + β∇^2u^(2) +
            (1 - β) ∇·τ^(2) = 0,
            
    ∇·u^(2) = 0, 
    
            (1 -  ∂_x)τ^(2) =
            𝐃^(2) - [
            u^(1)·∇τ^(1) -
            (∇u^(1))^T·τ^(1) -
            τ^(1)·∇u^(1)].


We have left out the boundary conditions, which we will deal with later.  Note
that in the second term of <ref>, the x-average commutes with
the y-substitution, as the substitute is independent of x.  Thus, we only
have to solve the x-average of <ref> to find the second order
swimming speed.  Considering the x-averages of the x-component of
<ref> and the xy-component of <ref>, we
have

    
    
    β∂_yy⟨ u^(2)⟩ +
            (1-β)∂_y⟨τ_xy^(2)⟩ = 0, 
    ⟨τ_xy^(2)⟩ =
            ∂_y⟨ u^(2)⟩ -
        [
               ⟨ u^(1)∂_xτ_xy^(1)⟩ +
            ⟨ v^(1)∂_yτ_xy^(1)⟩ - a
       ⟨ D_xy^(1)(τ_xx^(1)+τ_yy^(1))⟩ +
            ⟨Ω_xy^(1)(τ_xx^(1)-τ_yy^(1))⟩],


where Ω = (∇u^T-∇u)/2 is the vorticity
tensor, and we have ignored the x-averages of x-derivatives of x-periodic
functions, which must vanish.  Since the first order fields are known,
<ref> is simply an ordinary differential
equation for ⟨ u^(2)⟩, the solution to which is given by

    ⟨ u^(2)_±⟩ = E_± + F_± y +
        (1-β)^2/4(1+^2)[
            G_±cosh(2y) + H_±sinh(2y)
        ],

where

    G_± = (2 + 4 A_± + A_±^2 - B_±^2) -
            4 B_± (1 + A_±) y + 2(A_±^2 + B_±^2) y^2, 
        H_± = 2 A_± (B_± + 2 (1 + A_±) y - 2 B_± y^2).

Here, E_± and F_± are arbitrary constants.











Until now, the solutions in the domains above and below the sheet were
completely independent. Their coupling is now ensured by applying appropriate
boundary conditions, which determine the constants E_± and F_±. Similar
to the solution developed by Katz for swimming of a Taylor sheet next to a wall
in a Newtonian fluid <cit.>, we require that (i) the swimming speeds
we calculate from the upper and lower fluids match, and that (ii) the sheet is,
on average, force-free to second order in the x-direction.  The first order
flow field, which is the same as the first order Newtonian flow field, does not
apply a net-force to the sheet <cit.>, while the second order flow
field contributes a force
.Σ^(2)·n_±|_y=0 where n_± =
±_y + 𝒪(ϵ) is the inward-normal of the sheet in the
upper/lower fluid. To the second order, these boundary conditions are given by

    
    
    .⟨ u^(2)_±⟩|_y=± h_± = 0,
    ⟨. u^(1)_+|_y=y_s⟩ϵ +
            .⟨ u^(2)_+ ⟩|_y=0ϵ^2 =
            ⟨. u^(1)_-|_y=y_s⟩ϵ +
            .⟨ u^(2)_- ⟩|_y=0ϵ^2 +
            𝒪(ϵ^3), 
    .-β∂_y⟨ u^(2)_+⟩|_y=0
            +.(1-β)⟨τ^(2)_xy,+⟩|_y=0 =
            .-β∂_y⟨ u^(2)_-⟩|_y=0
            +.(1-β)⟨τ^(2)_xy,-⟩|_y=0,


which result in

    
    E_± = (1 + β^2)/(1 + ^2)h_±(A_∓ - A_±)/(h_+ + h_-)
             + (1 - β)^2/(1 + ^2)(A_±^2 - B_±^2), 
    
            F_± = (1 + β^2)/(1 + ^2)(A_∓ - A_±)/(h_+ + h_-)
             + (1 - β)^2/h_±(1 + ^2)(
            J_±cosh^2(h_±) + K_±sinh^2(h_±)
            ),


where

    J_± = (1 - 2h_± B_±) + h_±^2 (A_±^2 - B_±^2), 
        K_± = (1 + 2A_±)(1 - 2h_± B_±) +
            (1 - h_±^2)(A_±^2 - B_±^2).


Substituting the first and second order velocity fields into
<ref>, we finally arrive at

    U    = ϵ^2(1 + β^2)/(h_+ + h_-)(1 + ^2)[
        h_-(A_+ + 1/2) + h_+(A_- + 1/2)
        ] 
    U   = ϵ^2(1 + β^2)/2(h_+ + h_-)(1 + ^2)[
        h_-(sinh^2(h_+) + h^2_+/sinh^2(h_+) - h^2_+) +
        h_+(sinh^2(h_-) + h^2_-/sinh^2(h_-) - h^2_-)
        ].

In the limit of h_+→ hk and h_-→∞, we recover
<ref> as mentioned above.  The main implication of this result
is the observation that the effects of swimming next to a wall and swimming in
an Oldroyd-B fluid decouple at small wave amplitudes. The mechanism of this
decoupling is discussed in <ref>; but first we develop a
numerical method capable of calculating the swimming speed for any value of the
wave amplitude.










§ LARGE-AMPLITUDE SWIMMING: NUMERICAL METHOD



Here we develop a numerical method to solve <ref> subject to the
boundary conditions of <ref> for an arbitrary wave, with any
amplitude or shape.  As in the previous Section, we in fact solve the more
general problem of the sheet in a channel, with the walls above and below the
sheet a distance h_+ and h_- from the centreline, respectively.  We perform
the following transformation to the starred variables

    x^* = k(x - (c-U)t), y^* = k y, h_±^* = h_± k, 
    
        U^* = U/c, u^* = 1 - U^* + u/c, v^* = v/c, 
    
        p^* = p/η c k, τ^* = τ/η_p c k,
         Σ^* = Σ/η c k,

that render the problem dimensionless. This transformation is different from the
small-amplitude  transformation because in this frame of reference the shape of
the sheet is exactly independent of time, as opposed to being independent of
time only in the limit of small wave amplitudes.  Again, we drop the ^*s
in what follows.

To exploit the two-dimensional nature of the problem, we introduce a
stream-function ψ(x, y), which is defined via its relationship to the flow
field u: u = ∂_y ψ and v = -∂_x ψ.  This
substitution satisfies <ref> for any choice of ψ. To
re-formulate <ref> in terms of the stream-function, we take the
curl and divergence of  <ref> to obtain the complete set of
governing equations given by
 
    
    
    β∇^4 ψ - (1-β) [
                ∂_xy(τ_yy - τ_xx) +
                ^2 τ_xy] = 0, 
    ∇^2 p - (1-β) [
                ∂_xxτ_xx +
                2∂_xyτ_xy +
                ∂_yyτ_yy] = 0, 
    τ_xx - 2∂_xyψ + [
                (∂_yψ∂_x- ∂_xψ∂_y) τ_xx -
                2τ_xx∂_xyψ -
                2τ_xy∂_yyψ] = 0, 
    τ_xy + ^2 ψ  + [
                (∂_yψ∂_x- ∂_xψ∂_y) τ_xy +
                τ_xx∂_xxψ -
                τ_yy∂_yyψ] = 0,  
    τ_yy + 2∂_xyψ + [
                (∂_yψ∂_x- ∂_xψ∂_y) τ_yy +
                2τ_xy∂_xxψ +
                2τ_yy∂_xyψ]= 0,


where ^2 = ∂_xx - ∂_yy, β = η_s/(η_s +
η_p) is the viscosity ratio, and = λ c k is the Deborah number of
the fluid. There are five differential equations, three of which are
non-linear, with five fields to solve for (ψ, p, τ_xx,
τ_xy, τ_yy).

To solve <ref>  numerically, we developed a spectral
method adapted for our geometry, where the two-dimensional stream-function,
pressure and polymeric stress fields are represented by Fourier-Chebyshev
series <cit.>. Since convergence properties of the
Fourier-Chebyshev basis are optimal in rectangular domains, we perform two
independent coordinate transformations, one for the fluid above the sheet and
the other for the fluid below, that project the corresponding domains onto
rectangular strips, periodic in one direction. These transformations from the
original coordinates (x, y) to the new ones (η_±,ξ_±) are given by

    
    
    η_± = x, 
    ξ_± = 1 - 2± h_± - y/± h_± - f(x),


where “+” and “-” denote the regions above and below the sheet,
respectively.  In each domain, ξ_± = 1 corresponds to the domain's wall,
while ξ_± = -1 corresponds to the sheet, i.e. the lower domain has been
flipped.  The two domains can be treated equivalently and from now on we will
drop the ± unless it is necessary. The solutions in these domains only
couple to each other through the boundary conditions at the sheet.

In each deformed domain (η, ξ) ∈ [0,2π) × [-1, 1], the
hydrodynamic variables are represented by truncated Fourier-Chebyshev series.
For example, the pressure field, p, is given by

    p(η, ξ) = ∑_n=0^N-1∑_m=0^M-1 p^(nm)
        F_n(η) T_m(ξ),

where T_m(ξ) = cos(marccos(ξ)) is the mth Chebyshev polynomial, and

    F_n(η) = sin(n + 1/2η)     n odd
    cos(n/2η)     n even,
 is the nth Fourier mode.  We choose the resolution (N, M) such that the
error on truncating the series in <ref> is small; for each
set of physical parameters this is assessed by increasing the resolution (N,
M) until the value of the swimming speed of the sheet does not change by more
than 0.5% between the two highest resolutions. Typically, this precision is
achieved by N=33 and M=80.

The spatial derivaties ∂_η and ∂_ξ of the
Fourier-Chebyshev representations are calculated by multiplying vectors
containing spectral coefficients with the NM × NM spectral derivative
matrices <cit.>.  The spatial
derivatives in the original (x, y)-space are then trivially constructed with
the help of <ref>, giving for each domain

    
    
    (∂/∂ x)_± =
            ∂/∂η_± +
            (ξ_± - 1) f'(η_±)/± h_± - f(η_±)∂/∂ξ_±,
            
    (∂/∂ y)_± =
            2/± h_±-f(η_±)∂/∂ξ_±.


To calculate products of the fields represented in the spectral space, we use
the Fast Fourier transform <cit.> with the
following collocation points

    
    
    η_n_c = 2π n_c /N_c, 
    ξ_m_c = cos(π m_c/M_c-1),


to evaluate the fields in the real space, calculate their product, and transform
the result back to the spectral space. Here, n_c ∈[0,N_c), m_c
∈[0,M_c), and the collocation resolution (N_c, M_c) is selected
to satisfy N_c > 1.5 N and M_c > 1.5 M in order to avoid aliasing issues
<cit.>.

Representing five governing equations in the truncated Fourier-Chebyshev basis
for each fluid domain yields a set of 10NM non-linear algebraic equations that
need to be complemented by the boundary conditions.  By using Fourier modes, we
have implicitly imposed periodic boundary conditions in the η-direction,
which correctly reflects the symmetry of the underlying problem.  We still need,
however, six boundary conditions (four for ψ and two for p) along the
lines ξ=±1.  These boundary conditions are expanded in the Fourier basis
(as they are functions of η), generating 12N discretised boundary
conditions to substitute into the original set of 10NM discretised
governing equations.

The first boundary conditions to consider are the no-slip boundary conditions at
both the sheet and the wall, <ref>, where the velocities of
the material points of the sheet and the walls are given by <cit.>

    
    
    u_s(x) = - Q/√(1 + f'(x)^2),
            
    
            v_s(x) = - Qf'(x)/√(1 + f'(x)^2),
            
    u_w = (U - 1)_x,


and

    Q = ∫_0^2π√(1 + f'(x)^2) x.

The four boundary conditions are, therefore,

    
    
    . ∂_y ψ|_ξ = -1   = u_s,
            
    . -∂_x ψ|_ξ = -1   = v_s,
            
    . ∂_y ψ|_ξ = 1   = U - 1,
            
    . -∂_x ψ|_ξ = 1   = 0.


Note that the x-derivative of the n=0 Fourier mode vanishes, and that the
sheet's swimming speed U, which is unknown, appears in the n=0 mode of
<ref>, thus different sets of boundary conditions are required for
the n=0 and the n0 Fourier modes.  We address this below.  First, we
consider the other two boundary conditions required for the n0 case.

As already mentioned above, we do not directly solve the force balance equation,
<ref>, but instead solve its derivatives (specifically its curl and
divergence, see <ref>). The solutions to both problems may
differ, at most, by curl-free and divergence-free terms. To fix those terms, we
explicitly ensure that the force balance equation is satisfied at the
boundaries. Specifically, at both the sheet and the wall, we require that the
n·∇·Σ = 0, where n is the normal to the
surface. This yields the final two boundary conditions for the n0 Fourier
modes:

    
    [f'(η) ∂_x p -∂_y p +
                β  (f'(η)∂_y - ∂_x) ∇^2 ψ +
                x
    x + f'(η) ∂_x τ_xx +
                (f'(η) ∂_y + ∂_x) τ_xy +
                ∂_y τ_yy]_ξ = -1 = 0,

    
    [ ∂_y p + β∂_x ∇^2 ψ -
                ∂_x τ_xy - ∂_y τ_yy]_ξ = 1 = 0,


where (0, -1)^T is the normal to the wall, and (-f'(η), 1)^T the normal to
the surface of the sheet.

For the n=0 mode we replace
<ref>
with alternative boundary conditions.  First of all, we note that ψ and p
are defined up to a constant as only their derivatives are physical, and we set
those constants to some arbitrary value. The other two boundary conditions
ensure that the average x-force being applied to each of the walls is zero.
And similarly to the small amplitude case, we have to couple the two domains by
requiring that the swimming speed of the sheet as calculated by each domain is
the same and that the sheet is a force-free swimmer.  Thus, we have

    
    
    . p|_ξ=1 = 0, 
    .ψ|_ξ=1 = 0, 
    [ β^2ψ - τ_xy]_ξ=1 = 0,
            
    [ f'(η_+) p_+ - 2β f'(η_+)∂_xyψ_+ -
            β^2ψ_+ - f'(η_+)τ_xx, + + τ_xy, +]_ξ_+=-1
            = 
    [ f'(η_-) p_- - 2β f'(η_-)∂_xyψ_- -
            β^2ψ_- - f'(η_-)τ_xx, - + τ_xy, -]_ξ_-=-1,
            
    .∂_yψ_+|_ξ_+=1 =
        .∂_yψ_-|_ξ_-=1,


where the absence of a ± implies that the boundary condition applies to
both domains.

In the spirit of the Chebyshev-tau method <cit.>, for each Fourier
mode we replace the four highest Chebyshev modes of the discretised
<ref> and the two modes of <ref> with the boundary conditions
presented above. Combining everything together leads to the set of 10NM
non-linear discretised equations, with the structure outlined in
<ref>.  With the solution to this set of equations,  the swimming
speed of the sheet is given by

    U = .∂_y ψ|_ξ=1 + 1.


To actually solve this set of non-linear equations we use the Newton-Raphson
method <cit.> with an analytically calculated Jacobian. In general,
for > 0, starting from an arbitrary initial guess does not lead to
convergence of the Newton-Raphson algorithm, and therefore, we employ a simple
continuation strategy. For each set of parameters, we start from the Newtonian
case, = 0, which is linear and can always be solved, and use its solution
as the initial guess for a slightly higher . This process is continued
until we either reach the target value of  or the algorithm fails to
converge, in which case a smaller step Δ is selected. In practice,
the Δ required for continuation becomes very small at sufficiently large
, leading to unreasonable computation times in which case we only report
the results up to that value of .



We verify that our numerical method correctly reproduces the
small-amplitude prediction <ref>. In <ref> we plot
the swimming speed for a sheet with bk = 0.01 as a function of the Deborah
number  for various distances to the wall and viscosity ratios. As
expected, for this amplitude the numerically computed swimming speeds (symbols)
agree well with the analytical prediction of  <ref> (solid
line), again demonstrating that the effects of swimming next to a wall decouple
from the effects of swimming in a viscoelastic fluid at small amplitudes.



The situation changes significantly for finite values of the wave amplitude. In
<ref> we plot the swimming speed for a sheet with bk=0.5 as a
function of the Deborah number  for various distances to the wall and
viscosity ratios. We observe that the numerical data now deviates significantly
from the small-amplitude prediction <ref>. Despite this
deviation, for most values of h and β the swimming speed follows the
same trend as predicted by <ref>: starting from its Newtonian
value, it decreases monotonically with  and reaches a plateau value at
large Deborah numbers.
However, for sufficiently small h (hk =1.05, 1.1) at β = 0.5 and
    β = 0.9 the swimming speed breaks this trend and exhibits a
    non-monotonic dependence on .  This effect seems to be the stronger for
    larger values of β, which corresponds to more dilute solutions,
    reaching speeds faster than the Newtonian case for β = 0.9.  Also,
    there are indications that at lower values of β the swimming speed
    starts to increase with  at sufficiently large values of the Deborah
    number.
These results are further discussed in the
next Section.





§ DISCUSSION



As we have demonstrated above, at small wave amplitudes the influence of the
polymeric stress on the swimming speed is the same for both swimming in the bulk
and next to a wall. In other words, the effects of the boundary and polymers
decouple and the swimming speed is given by the product of the corresponding
contributions, see <ref>. Let us discuss the mechanism of this
behaviour.



We start by considering the kinematics of a Taylor sheet swimming in the bulk of
a Newtonian fluid. As noted by Taylor <cit.> and Lauga and Powers
<cit.>, at small wave amplitudes every point of the sheet is
oscillating approximately up and down, generating locally a vertical motion of
the surrounding fluid. Along one period of the sheet's waveform, for every point
moving up there is another point moving down with the same speed. Since the
fluid is incompressible, this sets an array of counter-rotating vortices along
the sheet, two vortices per period, see <ref>. We will be
referring to them as the `sheet vortices'. As can be seen from
<ref>, the presence of these vortices implies a velocity component
along the surface of the sheet, which on average drags the sheet along the
horizontal direction.

When this configuration of sheet vortices is placed next to a solid boundary, as
shown in <ref>, it generates a non-zero velocity at the boundary
which obviously does not satisfy the no-slip boundary condition at this
boundary. This velocity is cancelled by the creation of an array of vortices
localised at the boundary with the same periodicity as the sheet vortices.
Along the boundary, the velocity of these `wall vortices' has the same magnitude
but opposing direction of the velocity of the sheet vortices. And thus, the
no-slip boundary condition is satisfied for the total velocity field. In turn,
the wall vortices have a contribution along the surface of the sheet that
effectively increases the speed of the sheet vortices, which in turn causes the
sheet to speed up. Effectively, this implies that a small-amplitude sheet swimming
next to a wall can be viewed as a free-swimming sheet with faster sheet vortices.

This argument is further corroborated by rearranging the first order velocity
field <ref> in the following form, with dimensional
quantities
 
    u^(1) = c/2((A+B)ky-A-1)sin(k(x-ct))exp(-ky) +
    c/2((A-B)ky+A+1)sin(k(x-ct))exp(ky),
        
    
        v^(1) = c/2((A+B)ky+B-1)cos(k(x-ct))exp(-ky) +
    c/2((B-A)ky-B-1)cos(k(x-ct))exp(ky),

where we have dropped the ±, as the fluid domains are equivalent and the
distinction between them is unimportant.  In the upper domain, the terms
proportional to exp(-ky) and exp(ky) correspond to the sets of vortices
which are localised at the sheet and at the wall, respectively.  The centres of
these vortices are located along the lines ky=ω_s and ky=ω_w,
where

    ω_s = A+1/A+B =
        2sinh^2(hk)-h^2k^2/hk+sinh(hk)cosh(hk)+sinh^2(hk), 
    ω_w = A+1/B-A =
        2sinh^2(hk)-h^2k^2/hk+sinh(hk)cosh(hk)-sinh^2(hk).

For hk>1, ω_s≪ hk and ω_w≈ hk, thus justifying sheet and
wall label of the arrays of vortices. Note that the vertical location
of the sheet vortices is moved from its bulk location, ky=0, to the line
ky=ω_s≪ hk.  At small wave amplitudes, the effect of this displacement
is minor, while we suspect it plays an important role in finite-amplitude
swimming, see below.

Now we turn to the effect that the viscoelasticity of the fluid has on the
swimming speed. We have recently studied the mechanism of the polymer-induced
slowing down of a Taylor sheet in the bulk <cit.>, that we briefly
summarise here. As we have seen in <ref>, the
first order velocity field is the same for a Newtonian and Oldroyd-B fluids, and
we start by considering its effect on the polymeric stresses. While this
vortical velocity field is locally a simple shear flow in most parts of the
domain, the regions in-between the vortices correspond to locally extensional
flows around stagnation points, similar to the ones observed in PIV measurements
by Shen and Arratia <cit.>. These extensional flows generate large
normal components of the polymer stresses, τ_xx and τ_yy, that are
advected downstream by the fluid flow and, in turn,  generate significant shear
stresses in the fluid that drive the sheet horizontally, in the direction of the
wave propagation <cit.>. Since the Taylor sheet swims in
the direction opposite to the direction of its wave propagation, this results in
slowing down of the sheet compared to a Newtonian fluid. We note that for
swimmers that travel in the direction of the wave propagation, this effect would
enhance their propulsion speed, see also <cit.>. See also <cit.> for a  discussion of the interaction between the multiple waves propagating in the opposite directions and viscoelasticity.

Next to a boundary, the same mechanism applies at small wave amplitudes, since,
as we have argued above, a Newtonian sheet swimming next to a wall is equivalent
to a free-swimming sheet with faster sheet vortices and, hence, with a larger
swimming speed given by <ref>. This `effective' free-swimming sheet
would experience the same slowing down as discussed above and would swim with
the speed set by <ref>, where the Newtonian swimming speed c b^2
k^2/2 should be replaced with <ref>, arriving finally at
<ref>. This is the fundamental reason behind the factorisation of 
the effects of viscoelasticity and the boundary.

Now we turn to the case of finite-amplitude swimming. As shown in
<ref>, at bk=0.5 the swimming speed deviates significantly from
the small amplitude prediction of <ref>.  Although, in the
majority of cases the trend of the velocity decreasing with the Deborah number
persists, the high- value of the swimming speed U_∞ is larger than
the asymptotic prediction of <ref>, U_∞ > β U_N,
and increases as the boundary is brought closer to the swimmer.  Moreover, in
some cases the swimming speed no longer decreases monotonically and can even
increase to swimming speeds greater than the Newtonian value.

To understand this behaviour we analyse the spatial distribution of the elastic
stresses in the fluid. Additionally, we plot the flow type parameter χ,
defined as χ =
|𝐃|-|Ω|/|𝐃|+|Ω|
<cit.>.  Based on the invariants of the velocity gradient tensor, it
is designed to determine velocity type at every point in space, independent of
its local orientation: χ=1 corresponds to purely extensional,  χ=0 –
to shear, and χ=-1 – to purely rotation flows. Note that the flow type
parameter does not measure the magnitude of the flow, only its topology.



First we consider the case of moderate viscosity ratio and distance to the wall,
β=0.5 and h_+ k = 1.2. In  <ref> we plot the shear stress
τ_xy, the difference between normal stresses τ_xx-τ_yy, and
the flow type parameter for the Newtonian case =0 and two values of the
Deborah number, =0.5 and 3.2, corresponding to the monotonic decrease of
the swimming speed from its Newtonian value. The plots are superimposed with the
isolines of the stream-function ψ (black lines) and the local direction of
the velocity field (black arrows). Note that we plot the total velocity field
but will discuss it as consisting of the sheet and wall vortices, when useful.




In line with the small-amplitude mechanism discussed above, we observe formation
of lines of strong extensional flow that generate large normal stresses
τ_xx-τ_yy, which are advected by the local flow. The main difference
between this case and the small-amplitude one is the fact that the normal
stresses generated in-between the wall vortices are now advected by the vortices
towards the sheet and also contribute to the shear stress τ_xy that
generates an additional average flow that, in turn, drags the sheet in the
direction of the wave. At small Deborah number, =0.5, the negative value
(blue) of the normal stresses τ_xx -τ_yy is rotated into extra
positive (red) τ_xy, which is responsible for the slow down of the sheet
relative to swimming in a Newtonian fluid.  This only happens in the vortices in
the troughs of the sheet, where the vortices are not restricted too much by the
presence of the wall. At larger Deborah number, =3.2, in addition to a
region with extra positive τ_xy, there is a region with extra negative
τ_xy which pushes the sheet in the direction of its swimming, in
competition with the positive region.  The growth of this region of negative
polymeric shear stress is absent from the small amplitude solution (shown
    in <ref>) and is responsible for the increased swimming
speed compared to the small amplitude prediction <ref>.



Now we turn to the case of non-monotonic behaviour of the swimming speed with
the Deborah number.  In <ref> we compare the stress distributions
for β=0.5 and β=0.9 for h k = 1.05 and =2.7. For β=0.5
these values approximately correspond to the local maximum of the swimming
speed, although its value is still smaller than the Newtonian one, while at
these parameters the case with β=0.9 exhibits swimming speeds larger than
U_N, see <ref>. First, we note that now both the trough and crest
vortices are equally close to the wall, somewhat in contrast to
<ref>. This behaviour is also observed in the Newtonian case
=0, not shown. However, in  <ref> there are no significant
differences between the stress distributions for β=0.5 and β=0.9
cases besides the numerical values of the stresses, and we conclude that whether
the swimming speed is larger or smaller than its Newtonian counterpart is
determined by a numerical competition between the wall and sheet stresses that
cannot be deduced from hand-waving arguments.

In conclusion, we have provided a mechanistic explanation for the
small-wave-amplitude swimming speed of a Taylor sheet derived by Elfring and
Lauga <cit.>, and explain why the effects of fluid's viscoelasticity
and the presence of a solid boundary decouple. We also developed a numerical
method with spectral accuracy that allows us to study finite-amplitude sheets of
various waveforms close to and away from solid walls. We observe that at finite
amplitudes the swimming speed is no longer a monotonic function of the Deborah
number, and can even become larger than the corresponding Newtonian value.
Interestingly, this effect seems to be the stronger, the more dilute the
viscoelastic solution is (large values of β), although there are
indications that at lower values of β the swimming speed starts to
increase with  at sufficiently large values of the Deborah number. This
result suggests that even small amounts of polymer, either excreted or naturally
present in the solution, can aid propulsion next to solid boundaries, although
the speed increase reported in this work is minute. Our numerical data is not
sufficient to determine whether this increase would eventually lead to swimming
speeds larger than the Newtonian values for all values of β, at what
values of  this can be achieved, and how significant this speed up might
be. Further study is required to address these questions.



§ ACKNOWLEDGEMENTS


Research outputs generated through the EPSRC grant EP/I004262/1 can be found at
http://dx.doi.org/xxxxxx.

unsrt



