1 Introduction
In the natural environment, a variety of organisms negotiate their way through complex fluids ranging on a microscopic physiological scale from bacteria and sperm cells immersed in mucus (Fulford, Katz & Powell Reference Fulford, Katz and Powell1998; Bansil et al. Reference Bansil, Celli, Hardcastle and Turner2013), to the more familiar scale at which mudskippers swim through marshes and tidal flats and snails and slugs slide over mucus trails (Denny Reference Denny1980; McInroe et al. Reference McInroe, Astley, Gong, Kawano, Schiebel, Rieser, Choset, Blob and Goldman2016). These organisms must cope with the complications introduced by the ambient fluid rheology and, as a consequence, can develop special strategies for locomotion. Likewise, a number of artificial machines are designed to move through a non-Newtonian environment, including microbots targeted at diagnosis or drug delivery within the human body (Valdastri, Simi & Webster Reference Valdastri, Simi and Webster2012), and screw vehicles built to traverse challenging muddy terrain (Neumeyer & Jones Reference Neumeyer and Jones1965).
The fluid mechanics of locomotion through viscous fluids was pioneered by Taylor and Lighthill over half a century ago, with Taylor’s (Reference Taylor1951) flexible sheet providing one of the simplest and most idealized models of swimming through viscous fluids. Within the confines of this model, Taylor exploited a regular perturbation expansion in the limit of low-amplitude wavy motions on the flexible sheet to construct analytically the swimming speed. In a number of more recent works, the model has been generalized to explore swimming through viscoelastic, two-phase and generalized Newtonian fluids (Chaudhury Reference Chaudhury1979; Lauga Reference Lauga2007; Lauga & Powers Reference Lauga and Powers2009; Elfring, Pak & Lauga Reference Elfring, Pak and Lauga2010; Fu, Shenoy & Powers Reference Fu, Shenoy and Powers2010; Espinosa-Garcia, Lauga & Zenit Reference Espinosa-Garcia, Lauga and Zenit2013; Vélez-Cordero & Lauga Reference Vélez-Cordero and Lauga2013; Riley & Lauga Reference Riley and Lauga2014; Li & Ardekani Reference Li and Ardekani2015; Elfring & Goyal Reference Elfring and Goyal2016). The goal of the present work is to explore the consequence of a yield stress in the rheology of the ambient fluid.
The complications in generalizing Taylor’s model to a yield-stress fluid are significant because a finite stress is required to yield the fluid and permit propulsion. However, as the stress decays away from the swimmer, the fluid must remain rigid sufficiently far from the surface, leading to flow localization around the swimmer. The resultant yield surface between rigid and yielded fluid indicates that one must deal with a type of nonlinear free-boundary-value problem, in which the simplifications of Taylor’s low-amplitude perturbation method do not immediately apply. A similar situation is encountered for locomotion through granular media (Hosoi & Goldman Reference Hosoi and Goldman2015).
To tackle this complication, we offer a combined numerical and asymptotic approach to the problem. As in Taylor’s original problem, we consider an infinite two-dimensional sheet that deforms under prescribed wave-like motions. We consider two specific waveforms: a pure transverse wave (cf. Taylor Reference Taylor1951), discussed in § 3, and a pure longitudinal wave (cf. Tuck Reference Tuck1968), discussed in § 4. In each case, we examine the flow pattern, swimming speed and energetics for a range of different wave amplitudes and yield stresses. As in the Newtonian limit, transverse waves generate retrograde locomotion (propulsion in the opposite direction to the wave), while longitudinal waves lead to prograde motion (swimming in the same direction as the wave).
In the limit of low yield stress, the solutions approach the corresponding Newtonian results, with the flow field deviating only in the vicinity of a relatively distant yield surface. In the opposite limit of large yield stress, the leading-order formulation for transverse waves instead reduces to a problem in ideal plasticity and the stress solution can be found using the classical method of sliplines (the characteristic curves of the stress field). The construction for low wave amplitudes turns out to be closely related to Prandtl’s solution for the indentation of a punch into a plastic half-space, which is a classical problem in plasticity theory (Prager & Hodge Reference Prager and Hodge1951). For longitudinal waves, flow in the limit of large yield stress is instead dominated by a viscoplastic boundary layer (Balmforth et al. Reference Balmforth, Craster, Hewitt, Hormozi and Maleki2017) against the swimmer. We find the curious result that the speed in the plastic limit is exactly twice that for a swimmer in a Newtonian fluid, for both types of waves. Our numerical calculations provide solutions across a range of wave amplitudes and bridge the gap between these Newtonian and plastic limits.
2 Formulation
2.1 Governing equations
Consider a two-dimensional, incompressible Bingham fluid in an unbounded domain, containing a flexible sheet (a ‘swimmer’). Wave-like motions are introduced on the sheet, forcing the viscoplastic fluid to flow and generating locomotion. We consider a periodic section of the sheet, idealizing the swimmer as infinitely long with no significant effect from the head or tail. We neglect gravity and inertial forces, and focus only on fluid motion above the swimmer, assuming a reflectional symmetry about the
$x$
-axis. Using the Cartesian coordinate system
$(x,y)$
sketched in figure 1, the governing equations of incompressibility and force balance can be written in the dimensionless form,
$$\begin{eqnarray}u_{x}+v_{y}=0,\end{eqnarray}$$
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}_{xx}}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}_{xy}}{\unicode[STIX]{x2202}y},\quad \frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}y}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}_{xy}}{\unicode[STIX]{x2202}x}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}_{xx}}{\unicode[STIX]{x2202}y},\end{eqnarray}$$
in terms of the velocity
$(u,v)$
, pressure
$p$
and deviatoric-stress components
$(\unicode[STIX]{x1D70F}_{xx},\unicode[STIX]{x1D70F}_{xy})$
. Note that, except in the case of the stress components, we use subscripts of
$x$
and
$y$
(and, below, of
$t$
) to denote partial derivatives. The rheology is given by the Bingham constitutive relation
$$\begin{eqnarray}\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D70F}_{xx}\\ \unicode[STIX]{x1D70F}_{xy}\end{array}\right)=\left(1+\frac{Bi}{\dot{\unicode[STIX]{x1D6FE}}}\right)\left(\begin{array}{@{}c@{}}2u_{x}\\ v_{x}+u_{y}\end{array}\right)\quad \text{for }\unicode[STIX]{x1D70F}\equiv \sqrt{\unicode[STIX]{x1D70F}_{xx}^{2}+\unicode[STIX]{x1D70F}_{xy}^{2}}>Bi,\end{eqnarray}$$
and
$u_{x}=u_{y}+v_{x}=0$
otherwise, where
$\dot{\unicode[STIX]{x1D6FE}}=\sqrt{(v_{x}+u_{y})^{2}+4u_{x}^{2}}$
. To arrive at this scaled system, we have used the characteristic wavenumber of the swimming pattern
$k$
and wavespeed
$c$
to remove the dimensions of length and velocity; the stresses and pressure are scaled by
$\unicode[STIX]{x1D707}kc$
, resulting in the Bingham number or characteristic ratio of the yield and viscous shear stresses,
$$\begin{eqnarray}Bi=\frac{\unicode[STIX]{x1D70F}_{Y}}{\unicode[STIX]{x1D707}kc},\end{eqnarray}$$
where
$\unicode[STIX]{x1D707}$
is the (plastic) viscosity and
$\unicode[STIX]{x1D70F}_{Y}$
the yield stress.
Sketch of the geometry in the frame of the swimmer, for (a) retrograde transverse waves and (b) prograde longitudinal waves.
We consider either transverse or longitudinal wave motions of the sheet. In the frame of the swimmer, travelling at speed
$U$
, we locate the surface of the sheet by
$$\begin{eqnarray}\left(\begin{array}{@{}c@{}}x\\ y\end{array}\right)=\left[\begin{array}{@{}c@{}}x_{0}+X(x_{0},t)\\ Y(x_{0},t)\end{array}\right],\end{eqnarray}$$
in terms of two waveform functions
$X(x_{0},t)$
and
$Y(x_{0},t)$
, where
$x_{0}$
is a Lagrangian coordinate of a material point on the surface. We assume that the swimmer is freely extensible (cf. Blake Reference Blake1971; Katz Reference Katz1974), such that the corresponding surface velocity is
$$\begin{eqnarray}\left(\begin{array}{@{}c@{}}u\\ v\end{array}\right)=\left(\begin{array}{@{}c@{}}X_{t}\\ Y_{t}\end{array}\right)\,.\end{eqnarray}$$
Note that Taylor originally modelled the opposite limit of a perfectly inextensible swimmer, for which the boundary conditions reduce to the same limit when
$a\ll 1$
but differ when the amplitude is large (see § 3.1).
For transverse waves, we set
$X=0$
and
$x=x_{0}$
, so that the surface
$[x,Y(x,t)]$
moves purely vertically. For longitudinal waves, we instead set
$Y=0$
so that material surface points shift only horizontally (see figure 1). We assume sinusoidal waveforms that travel to the left and set
$$\begin{eqnarray}[X,Y]=[0,a\sin \unicode[STIX]{x03C0}(x+t)]\quad \text{or}\quad [X,Y]=[a\sin \unicode[STIX]{x03C0}(x_{0}+t),0],\end{eqnarray}$$
respectively, where
$a$
is the wave amplitude measured in units of wavelength. Note that, in either case, the dimensionless wavelength is
$2$
, leading to the periodic domain
$-1\leqslant x\leqslant 1$
.
Sufficiently far above the swimmer, the fluid is unyielded and, in this frame, translates to the left at the swimming speed:
$u(x,y\rightarrow \infty )=-U$
(measured in units of wavespeed). The swimming speed is not known a priori but must be determined as part of the solution. To calculate the speed, we impose the constraints that there is no pressure drop across each periodic section of the channel,
$$\begin{eqnarray}0=\left[p(x,y)\right]_{x=-1}^{1}\equiv \left[\unicode[STIX]{x1D70F}_{xx}(x,y)\right]_{x=-1}^{1}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}\int _{-1}^{1}\unicode[STIX]{x1D70F}_{xy}(x,y)\,\text{d}x,\end{eqnarray}$$
along any horizontal level
$y=y_{I}>Y$
, and there is no net horizontal force on the swimmer,
$$\begin{eqnarray}0=\int _{-1}^{1}\left[\unicode[STIX]{x1D70F}_{xy}+Y_{x}(p-\unicode[STIX]{x1D70F}_{xx})\right]_{y=Y}\text{d}x\equiv \int _{-1}^{1}\unicode[STIX]{x1D70F}_{xy}(x,y_{I})\,\text{d}x.\end{eqnarray}$$
The last relation follows from integrating (2.2a
) over the area enclosed by
$y=Y(x,t)$
and
$y=y_{I}$
, turning the double integrals into surface integrals using Gauss’s theorem, and cancelling the contributions from the vertical sides of the domain in view of the periodic boundary conditions and the vanishing pressure drop. Provided
$\unicode[STIX]{x1D70F}_{xx}$
is periodic, equation (2.9) implies (2.8), and therefore does not need to be imposed independently.
The solution must also satisfy the power integral,
$$\begin{eqnarray}\displaystyle \langle \dot{\unicode[STIX]{x1D6FE}}(Bi+\dot{\unicode[STIX]{x1D6FE}})\rangle & \equiv & \displaystyle {\mathcal{P}}\equiv \int _{-1}^{1}\left[\boldsymbol{u}\boldsymbol{\cdot }\left(\unicode[STIX]{x1D749}-p\boldsymbol{I}\right)\boldsymbol{\cdot }\boldsymbol{n}\right]_{y=Y}\sqrt{1+Y_{x}^{2}}\,\text{d}x\nonumber\\ \displaystyle & \equiv & \displaystyle \int _{-1}^{1}\left[X_{t}\unicode[STIX]{x1D70F}_{xy}+X_{t}Y_{x}(p-\unicode[STIX]{x1D70F}_{xx})-Y_{t}Y_{x}\unicode[STIX]{x1D70F}_{xy}-Y_{t}(p+\unicode[STIX]{x1D70F}_{xx})\right]_{y=Y}\,\text{d}x,\end{eqnarray}$$
where the angular brackets denote an area integral over the entire domain and
$\boldsymbol{n}$
is the (upward) normal vector to the sheet. The left-hand side of (2.10) is the dissipation rate of the yield and viscous stress, while
${\mathcal{P}}$
is the total power input from the swimming stroke. In the frame of the wave (where there is no normal velocity for either transverse or longitudinal waveforms), the latter can alternatively be written as
$$\begin{eqnarray}{\mathcal{P}}=\int _{0}^{S}w\unicode[STIX]{x1D70F}_{sn}\;\text{d}s,\end{eqnarray}$$
where
$\unicode[STIX]{x1D70F}_{sn}=[(1-Y_{x}^{2})\unicode[STIX]{x1D70F}_{xy}-2Y_{x}\unicode[STIX]{x1D70F}_{xx}]_{y=Y}/(1+Y_{x}^{2})$
and
$w$
are the local shear stress and tangential speed at the swimmer surface, respectively, and
$s$
is the arc length of that boundary, as measured from
$x=-1$
, with total perimeter
$S$
. For transverse waves,
$w=\sqrt{1+Y_{x}^{2}}$
and
$\text{d}s=\sqrt{1+Y_{x}^{2}}\;\text{d}x$
. For longitudinal waves,
$w=X_{t}-1$
and
$\text{d}s=\text{d}x$
.
Given the power input
${\mathcal{P}}$
, we can formulate a measure of the efficiency
${\mathcal{E}}$
of the swimming stroke,
$$\begin{eqnarray}{\mathcal{E}}=\frac{U^{2}}{{\mathcal{P}}}.\end{eqnarray}$$
(cf. Blake Reference Blake1971), which provides a ratio of the characteristic thrusting power to the total power expended. Note that some authors have questioned the use of
$U^{2}$
to estimate the thrusting power, preferring to define the speed per unit power
$U/{\mathcal{P}}$
as a measure of ‘swimming economy’ (cf. Krieger, Spagnolie & Powers Reference Krieger, Spagnolie and Powers2014).
Finally, we note that sufficiently far above the swimmer, the stress in the fluid must fall below the yield stress to form a rigid plug. Unlike in a Newtonian fluid, therefore, rigid boundaries in the fluid have no effect on the flow if they are located above this yield distance. The model presented here is thus equally applicable to locomotion in the vicinity of a sufficiently distant rigid boundary. We exploit this feature in our numerical simulations in order to simulate the flow in a finite domain.
2.2 Perfect plasticity: slipline theory
For
$Bi\gg 1$
, one anticipates that regions of perfectly plastic flow may arise over which viscous stresses
$\dot{\unicode[STIX]{x1D6FE}}_{ij}$
are negligible relative to plastic stresses
$Bi\dot{\unicode[STIX]{x1D6FE}}_{ij}/\dot{\unicode[STIX]{x1D6FE}}$
and the stress invariant is held just above the yield stress,
$\unicode[STIX]{x1D70F}\sim Bi$
. In this situation, equation (2.2) reduces to the hyperbolic equations of perfect plasticity, which can be solved using the method of characteristics, or ‘sliplines’, in the more commonly used terminology of this field (e.g. Prager & Hodge Reference Prager and Hodge1951).
Within the region of perfectly plastic flow, the characteristics of the stress field are given by two orthogonal families of sliplines, which we denote by
$\unicode[STIX]{x1D6FC}$
and
$\unicode[STIX]{x1D6FD}$
. To incorporate the plastic constraint
$\unicode[STIX]{x1D70F}=Bi$
, we represent the stress field by
$$\begin{eqnarray}(\unicode[STIX]{x1D70F}_{xx},\unicode[STIX]{x1D70F}_{xy})=Bi(-\!\sin 2\unicode[STIX]{x1D719},\cos 2\unicode[STIX]{x1D719}),\end{eqnarray}$$
in which the angle
$\unicode[STIX]{x1D719}$
can be related directly to the properties of the sliplines. More specifically, the characteristics of (2.2) comprise
$\unicode[STIX]{x1D6FC}-$
lines, which make an angle
$\unicode[STIX]{x1D719}$
with the
$x$
-axis and along which the quantity
$p+2Bi\unicode[STIX]{x1D719}$
is invariant, and
$\unicode[STIX]{x1D6FD}-$
lines, which are orthogonal to the
$\unicode[STIX]{x1D6FC}-$
lines and have the invariant
$p-2Bi\unicode[STIX]{x1D719}$
.
The plastic flow lines (i.e. the characteristics of the velocity field) coincide with the sliplines. They can be calculated from the constraint of incompressibility (2.1), which reduces to the so-called Geiringer equations,
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}v_{\unicode[STIX]{x1D6FC}}}{\unicode[STIX]{x2202}s_{\unicode[STIX]{x1D6FC}}}=v_{\unicode[STIX]{x1D6FD}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}s_{\unicode[STIX]{x1D6FC}}},\quad \text{and}\quad \frac{\unicode[STIX]{x2202}v_{\unicode[STIX]{x1D6FD}}}{\unicode[STIX]{x2202}s_{\unicode[STIX]{x1D6FD}}}=-v_{\unicode[STIX]{x1D6FC}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}s_{\unicode[STIX]{x1D6FD}}},\end{eqnarray}$$
along the
$\unicode[STIX]{x1D6FC}$
and
$\unicode[STIX]{x1D6FD}$
lines, respectively. Here,
$(v_{\unicode[STIX]{x1D6FC}},v_{\unicode[STIX]{x1D6FD}})$
denote the velocity components along the sliplines and
$(s_{\unicode[STIX]{x1D6FC}},s_{\unicode[STIX]{x1D6FD}})$
denote their arc-length coordinates.
2.3 Numerical method
To satisfy mass conservation (2.1), we introduce a streamfunction
$\unicode[STIX]{x1D713}$
with convention
$(u,v)=(\unicode[STIX]{x1D713}_{y},-\unicode[STIX]{x1D713}_{x})$
. After eliminating the pressure in (2.2)–(2.3), the equations reduce to
$$\begin{eqnarray}\frac{\unicode[STIX]{x1D6FB}^{4}\unicode[STIX]{x1D713}}{Bi}+4\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}y}\left(\frac{\unicode[STIX]{x1D713}_{xy}}{\dot{\unicode[STIX]{x1D6FE}}}\right)+\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}x^{2}}-\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}y^{2}}\right)\left(\frac{\unicode[STIX]{x1D713}_{xx}-\unicode[STIX]{x1D713}_{yy}}{\dot{\unicode[STIX]{x1D6FE}}}\right)=0;\quad \dot{\unicode[STIX]{x1D6FE}}^{2}=4\unicode[STIX]{x1D713}_{xy}^{2}+(\unicode[STIX]{x1D713}_{xx}-\unicode[STIX]{x1D713}_{yy})^{2},\end{eqnarray}$$
provided
$\unicode[STIX]{x1D70F}>Bi$
, and
$\dot{\unicode[STIX]{x1D6FE}}=0$
otherwise. We work in the frame of the swimmer and freeze time at
$t=0$
because it does not explicitly enter the problem. As noted above, we fix the upper boundary of the domain at
$y=H$
, for some constant
$H$
that is always selected to lie above the yield surface and thus plays no role.
The boundary conditions on the sheet are
$$\begin{eqnarray}\unicode[STIX]{x1D713}=-a\sin \unicode[STIX]{x03C0}x,\quad \unicode[STIX]{x1D713}_{y}=0,\quad \text{on }y=Y(x)=a\sin \unicode[STIX]{x03C0}x,\end{eqnarray}$$
for transverse waves and
$$\begin{eqnarray}\unicode[STIX]{x1D713}=0,\quad \unicode[STIX]{x1D713}_{y}=a\unicode[STIX]{x03C0}\cos \unicode[STIX]{x03C0}x_{0},\quad \text{on }x=x_{0}+a\sin \unicode[STIX]{x03C0}x_{0}\text{ and }y=0,\end{eqnarray}$$
for longitudinal waves. In either case, the flow is periodic in the
$x$
direction over the domain
$-1\leqslant x\leqslant 1$
and we impose the conditions on the upper surface,
$$\begin{eqnarray}\unicode[STIX]{x1D713}=K,\quad \unicode[STIX]{x1D713}_{yy}=0,\quad \text{on }y=H,\end{eqnarray}$$
where the unknown constant
$K$
in (2.18) determines the total horizontal flux through the domain. We determine
$K$
by calculating the pressure drop across the domain (2.9) and iterating the value of
$K$
using an interval-bisection algorithm until both the pressure drop and the corresponding uncertainty in
$K$
fall below a prescribed tolerance (for most computations these are
$10^{-4}Bi$
and
$10^{-5}$
, respectively). For each value of
$K$
, we use an augmented Lagrangian scheme (Fortin & Glowinski Reference Fortin and Glowinski2006) to solve (2.15) numerically. The basic details of this scheme are outlined in the Appendix, while further information can be found in Liu et al. (Reference Liu, Balmforth, Hormozi and Hewitt2016) and Balmforth et al. (Reference Balmforth, Craster, Hewitt, Hormozi and Maleki2017). We confirmed explicitly that the net horizontal force along the swimmer (2.8) vanishes for the final converged value of
$K$
. Finally, we read off the swimming speed from
$U\equiv -\unicode[STIX]{x1D713}_{y}(x,H)$
, given that the fluid is fully plugged at the top boundary.
3 Transverse wave solutions
3.1 Numerical observations
Figures 2–3 show sample numerical solutions with a transverse waveform (as in figure 1
a) across a range of values of
$Bi$
and
$a$
. For small
$Bi$
, the solutions take the form of a relatively deep yielded zone over which the flow is much like that in the Newtonian problem (left-hand panels in figures 2–3). For small amplitude
$a$
(figure 2), the streamlines (in the frame of the wave, wherein the swimmer has a fixed shape) follow the topography of the swimmer, implying no net entrainment of fluid. At higher amplitudes (figure 3), the streamlines develop recirculation zones in the topographic troughs, which carry certain fluid regions along with the swimmer. In all these near-Newtonian solutions, the yield stress only becomes important sufficiently far above the swimmer, causing the fluid to plug up at a relatively distant yield surface.
As the yield stress is increased, the yielded region becomes more localized to the swimmer (central panels in figures 2–3). With the highest yield stresses (right-hand panels in figures 2–3), the flow field develops distinctive structure, containing a melange of sharp boundary layers, wider ‘plastic’ regions with relatively low (but finite) strain rate and disconnected rigidly rotating plugs.
Solutions for low and moderate amplitude
$a$
, showing density maps of
$\log _{10}\dot{\unicode[STIX]{x1D6FE}}$
with streamlines in the wave frame superposed (blue). (a)–(c)
$a=0.1$
and (d)–(f)
$a=0.4$
, with
$Bi=0.25$
on (a,d),
$Bi=4$
in (b,e) and
$Bi=2048$
on (c,f). The swimmer is shaded grey, and black regions show the unyielded plugs.
Solutions for moderate and high amplitude
$a$
, showing density maps of
$\log _{10}\dot{\unicode[STIX]{x1D6FE}}$
with streamlines in the wave frame superposed (blue). (a)–(c)
$a=0.5$
and (d)–(f)
$a=1$
, with
$Bi=0.25$
on (a,d),
$Bi=4$
in (b,e) and
$Bi=2048$
on (c,f). The swimmer is shaded grey, and black regions show the unyielded plugs. The black dashed line in (f) shows the centreline of the boundary layer
$y=Y+\unicode[STIX]{x1D6E5}$
computed from (3.18).
For a fixed, high yield stress, the solutions pass through an interesting change in the flow pattern as the wave amplitude is increased. For low
$a$
, the swimmer is sheathed by a relatively wide (order-one) region of predominantly plastic deformation (figure 2
c,f). As the wave amplitude grows, rigidly rotating plugs within this region expand until the topography of the swimmer exceeds the natural thickness of the nearly plastic region. At that stage, there is a switch to a different type of flow pattern with a more boundary-layer-like form (figure 3
c,f). At first glance, one might assume that the boundary layers are thin viscoplastic regions, in which viscous stresses enter the force balance (cf. Balmforth et al.
Reference Balmforth, Craster, Hewitt, Hormozi and Maleki2017). However, a dissection of the numerical solutions demonstrates that the viscous part of the stress does not play a significant role in the force balance anywhere for
$Bi\gg 1$
. This observation is confirmed by measurements of the boundary-layer width in figure 4(f), discussed below. In other words, the boundary-layer-like regions against the swimmer at high amplitude (figure 3
c,f) are actually slender regions of plastic deformation, indicating that solutions undergo a transition from a plastic indentation flow at low
$a$
to a plastic boundary-layer flow at high
$a$
(see §§ 3.3–3.4).
(a,c,f) Data as a function of
$Bi$
for different fixed amplitudes:
$a=0.025$
,
$a=0.1$
,
$a=0.2$
,
$a=0.3$
,
$a=0.5$
,
$a=1$
and
$a=2$
. (b,d) Data as a function of
$a$
for different fixed Bingham numbers
$Bi$
:
$Bi=1/16$
,
$Bi=1/4$
,
$Bi=1$
,
$Bi=16$
,
$Bi=128$
and
$Bi=2048$
. For each plot, the series with the lowest (highest) value of the parameter (
$a$
or
$Bi$
) is marked with stars (open circles) for clarity. (a,b) Swimming speed
$U$
; (c) scaled power consumption
${\mathcal{P}}/a$
and (d) unscaled power consumption
${\mathcal{P}}$
; (e) efficiency
${\mathcal{E}}=U^{2}/{\mathcal{P}}$
as a density map in
$(Bi,a)$
space, together with three contours of constant economy
$U/{\mathcal{P}}$
; and (f) distance
$\unicode[STIX]{x1D6FF}$
from the top of the topography at
$x=0.5$
to the yield surface. Short-dashed (blue) lines and long-dashed (red) lines show low-amplitude predictions (
$a\ll 1$
) in the limits
$Bi\ll 1$
and
$Bi\gg 1$
, respectively. The green dot-dashed line in (b) shows the results of Sauzade, Elfring & Lauga (Reference Sauzade, Elfring and Lauga2011) for an inextensible swimmer in a Newtonian fluid, while that in (d) shows the high-
$Bi$
, high-
$a$
, asymptotic prediction (3.15) for
$Bi=2048$
.
Data extracted from the numerical calculations are shown in figure 4. For given wave amplitude, the swimming speed
$U$
increases with
$Bi$
away from its Newtonian limit, before converging to a limit
$U\rightarrow U_{\infty }(a)$
for
$B\rightarrow \infty$
(figure 4
a). At low
$a$
,
$U$
progresses monotonically up to
$U_{\infty }$
, whereas for higher amplitudes,
$U$
attains a maximum at a finite
$Bi$
and then decreases towards
$U_{\infty }$
. Plotting the swimming speed against
$a$
for given
$Bi$
(figure 4
b), shows that the speed converges to Taylor’s limit for Newtonian fluid
$U\sim \unicode[STIX]{x03C0}^{2}a^{2}/2$
in the limit
$(a,Bi)\ll 1$
. In the limit
$a\ll 1$
and
$Bi\gg 1$
, we instead observe that
$U_{\infty }\sim \unicode[STIX]{x03C0}^{2}a^{2}$
. In other words, the swimming speed for low-amplitude waves driving motion through an ideal plastic is exactly double that for motion through a Newtonian fluid. These observations are rationalized in the following subsections.
The power consumption
${\mathcal{P}}$
in (2.10) increases monotonically with both
$a$
and
$Bi$
(figure 4
c,d), reflecting how more fluid is forced to flow and a greater yield stress must be overcome, respectively. In the plastic limit,
${\mathcal{P}}$
increases linearly with the combination
$aBi$
. Given the power consumption and swimming speed, we can quantify the effectiveness of the swimming strategy using Blake’s definition of efficiency
${\mathcal{E}}=U^{2}/{\mathcal{P}}$
(2.12), as shown in figure 4(e). The efficiency initially increases with both
$a$
and
$Bi$
due to the rise of the swimming speed as the swimmer gains increased traction on the ambient fluid. The increase in speed, however, is limited at higher
$a$
and
$Bi$
, and eventually the gradual rise in the power expenditure reduces the efficiency. An optimal yield stress and amplitude for swimming thereby arises for
$(a,Bi)\approx (0.5,4)$
(figure 4
e), suggesting a potential advantage for biological organisms to tune their locomotion strategy to the rheology of their environment, as has previously been proposed for swimming in viscoelastic fluid. Note, however, that the swimming economy
$U/{\mathcal{P}}$
defined by Krieger et al. (Reference Krieger, Spagnolie and Powers2014) (contours in figure 4
e) is maximized in the Newtonian limit with
$a\rightarrow 0$
.
Figure 4(f) shows the thickness
$\unicode[STIX]{x1D6FF}$
of the yielded layer above the crest of the waves (
$x=1/2$
), and illustrates the transition in the flow structure at high yield stress as the wave amplitude is increased. For low
$a$
, this distance approaches
$1/\sqrt{2}$
, as predicted in § 3.3 below. For
$a\approx 0.4$
and higher, however, the thickness instead approaches another, much smaller limit, indicative of the overlying boundary layer. Crucially,
$\unicode[STIX]{x1D6FF}$
is independent of
$Bi$
, which is only possible if viscous resistance is negligible and the boundary layers are slender regions of almost plastic flow.
Finally, figure 4(b) also includes a set of data taken from Sauzade et al. (Reference Sauzade, Elfring and Lauga2011), who computed swimming speeds for Taylor’s sheet in a Newtonian fluid up to comparable values of
$a$
as in our calculations. While their data agree with our low-
$Bi$
results when the wave amplitude is small, the predictions diverge for large
$a$
. This disagreement results because Sauzade et al. (Reference Sauzade, Elfring and Lauga2011) demand that the swimmer is inextensible, as in Taylor’s original model, whereas the swimmer is freely extensible in our model. In reality, both models are crude idealizations of any real swimmer when the swimming stroke is large; we explore solutions for
$a\gg 1$
chiefly to understand the nature of the limiting behaviour rather than to provide a direct model of high-amplitude locomotion.
3.2 Low-amplitude waves in nearly Newtonian fluid (
$a\ll 1$
,
$Bi\ll 1$
)
In the low-amplitude Newtonian limit, the streamfunction satisfies the biharmonic equation and the swimmer becomes almost flat, so that the boundary condition can be imposed perturbatively about
$y=0$
. This device allowed Taylor (Reference Taylor1951) to find the swimming speed for the Newtonian problem. We repeat this analysis here, but for the weakly viscoplastic problem.
3.2.1 Near-field Newtonian region
For
$a\ll 1$
and
$Bi=O(a^{2})$
, a regular perturbation expansion of (2.15) yields Taylor’s biharmonic problem at leading order. With this solution in hand, one can then continue on to next order and simultaneously include the yield stress and terms arising from the shift in the position of the swimmer. The result is
$$\begin{eqnarray}\unicode[STIX]{x1D713}=-\left(a+\frac{Bi}{\unicode[STIX]{x03C0}^{2}}\right)(1+\unicode[STIX]{x03C0}y)\text{e}^{-\unicode[STIX]{x03C0}y}\sin \unicode[STIX]{x03C0}x+\frac{Bi}{\unicode[STIX]{x03C0}^{2}}\sin \unicode[STIX]{x03C0}x-\frac{1}{2}\unicode[STIX]{x03C0}^{2}a^{2}y(1-\text{e}^{-2\unicode[STIX]{x03C0}y}\cos 2\unicode[STIX]{x03C0}x)+O(a^{3}).\end{eqnarray}$$
If
$Bi=0$
, the
$O(a^{2})$
contribution to this solutions gives the far-field Newtonian swimming speed of
$U=\unicode[STIX]{x03C0}^{2}a^{2}/2$
(see figure 4
b). With this expression for
$\unicode[STIX]{x1D713}$
, one can also calculate the leading-order strain rate and dissipation rate (2.10),
$$\begin{eqnarray}\dot{\unicode[STIX]{x1D6FE}}\sim 2\unicode[STIX]{x03C0}^{3}ay\text{e}^{-\unicode[STIX]{x03C0}y}\quad \text{and}\quad {\mathcal{P}}\sim 2\unicode[STIX]{x03C0}^{3}a^{2}.\end{eqnarray}$$
Both the plastic dissipation
$Bi\langle \dot{\unicode[STIX]{x1D6FE}}\rangle$
and the
$O(Bi)$
correction to the viscous dissipation furnish a correction to
${\mathcal{P}}$
that is of order
$aBi$
. Thus,
${\mathcal{P}}\sim 2\unicode[STIX]{x03C0}^{3}a^{2}+O(aBi)$
, as seen in figure 4(d).
(a) Numerical solution with a flat base (i.e. boundary conditions imposed on
$y=0$
rather than
$y=Y(x)$
, implying zero swimming speed), for
$Bi=1/8$
and
$a=0.05$
, showing a density plot of
$\log _{10}\dot{\unicode[STIX]{x1D6FE}}$
together with streamlines (blue, solid) and the leading-order prediction from (3.1) (white dashed). Streamlines are shown in the swimmer frame. (b) A magnification of the streamlines for the same solution near the yield surface, where the flow deviates from (3.1) and is given instead by the parameter-free equation (3.4). (c) Yield surfaces, after subtracting their mean height
$D$
(see figure 6), for a set of solutions of the full problem (i.e. non-flat bottom boundary) for different values of
$Bi$
and
$a$
(
$Bi=2^{-n}$
,
$n=1,2,3,4$
;
$a=0.025$
and
$a=0.05$
).
3.2.2 Far-field viscoplastic region
In view of the non-decaying yield-stress term
$\unicode[STIX]{x03C0}^{-2}Bi\sin \unicode[STIX]{x03C0}x$
, the regular solution in (3.1) becomes disordered for
$$\begin{eqnarray}y\sim O(y_{B})\quad \text{with }y_{B}=\frac{1}{\unicode[STIX]{x03C0}}\log \left(\frac{\unicode[STIX]{x03C0}^{3}ay_{B}}{Bi}\right),\end{eqnarray}$$
and the prediction in (3.1) breaks down (see figure 5
a). Above this level a new expansion is needed in which the yield-stress features in the leading-order balance of terms in (2.15). On defining
$\unicode[STIX]{x1D6F9}(x,\check{y})=(\unicode[STIX]{x1D713}+\unicode[STIX]{x03C0}^{2}a^{2}y/2)/Bi$
, with
$\check{y}=y-y_{B}$
(which eliminates the swimming-speed term
$-\unicode[STIX]{x03C0}^{2}a^{2}y/2$
that dominates the streamfunction above
$y_{B}$
but gives no contribution to the strain rate), we find that
$\unicode[STIX]{x1D6F9}$
satisfies the far-field problem,
$$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{4}\unicode[STIX]{x1D6F9}+4\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}\check{y}}\left(\frac{\unicode[STIX]{x1D6F9}_{x\check{y}}}{\unicode[STIX]{x1D6E4}}\right)+\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}x^{2}}-\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}\check{y}^{2}}\right)\left(\frac{\unicode[STIX]{x1D6F9}_{xx}-\unicode[STIX]{x1D6F9}_{\check{y}\check{y}}}{\unicode[STIX]{x1D6E4}}\right)=0,\quad \unicode[STIX]{x1D6E4}^{2}=4\unicode[STIX]{x1D6F9}_{x\check{y}}^{2}+(\unicode[STIX]{x1D6F9}_{xx}-\unicode[STIX]{x1D6F9}_{\check{y}\check{y}})^{2},\end{eqnarray}$$
subject to the matching condition
$$\begin{eqnarray}\unicode[STIX]{x1D6F9}\rightarrow \frac{1}{\unicode[STIX]{x03C0}^{2}}\text{e}^{-\unicode[STIX]{x03C0}\check{y}}\sin \unicode[STIX]{x03C0}x\quad \text{for }\check{y}\rightarrow O(-y_{B}),\end{eqnarray}$$
corresponding to
$\unicode[STIX]{x1D713}\sim a\unicode[STIX]{x03C0}y\text{e}^{-\unicode[STIX]{x03C0}y}\sin \unicode[STIX]{x03C0}x$
for
$y=O(1)$
.
Equations (3.4)–(3.5) comprise a parameter-free nonlinear problem that corrects (3.1) to allow the yield stress to assert itself and plug up the fluid. The yield surface is thus located at a height
$y=O(y_{B})$
(see figure 6). Moreover, aside from a scaling, the flow in the far field adopts a universal form, as illustrated in figure 5(b). In particular, the yield surfaces collapse once shifted by
$y_{B}$
(figure 5
c).
The average distance
$D$
above
$y=0$
of the yield surface, for
$a=0.1$
(black stars),
$a=0.05$
(red circles) and
$a=0.025$
(blue squares). The solid (black) line shows
$D=y_{B}$
, as defined implicitly in (3.2) and the dashed (blue) line shows the prediction of the perfectly plastic theory
$D=0.643$
, given by the average height of the yield surface in figure 7.
3.3 Low-amplitude waves in the perfectly plastic limit (
$a\ll 1$
,
$Bi\gg 1$
)
The low-amplitude limit can also be reconsidered in the opposite limit,
$Bi\gg 1$
, in which case the leading-order problem corresponds to the indentation of a perfectly plastic half-space as considered in the early plasticity literature by Prandtl (see Prager & Hodge Reference Prager and Hodge1951). As in that literature, we use slipline theory (see § 2.2) to construct solutions in this limit.
3.3.1 Leading-order solutions
Guided by the numerical solutions, we observe that the leading-order stress field can be decomposed spatially into three parts in the limit
$Bi\gg 1$
and
$a\ll 1$
. These parts are illustrated in figure 7 and comprise two fans (
$F_{\pm }$
) centred at
$[x,y]=[\pm 1/2,0]$
, buffered from one another by triangles of constant stress (
$T_{\pm }$
), all beneath an overlying plug (
$P$
).
Low-amplitude plastic slipline solution, showing the plug
$P$
, triangular checkerboards
$T_{\pm }$
and fans
$F_{\pm }$
, together with a selection of sliplines (red and blue for the
$\unicode[STIX]{x1D6FC}$
and
$\unicode[STIX]{x1D6FD}-$
lines, respectively).
Within the triangle
$T_{+}$
above
$-1/2<x<1/2$
, the sliplines are inclined at
$45^{\circ }$
, such that their local arc-length coordinates
$(s_{\unicode[STIX]{x1D6FC}},s_{\unicode[STIX]{x1D6FD}})=(\unicode[STIX]{x1D709},\unicode[STIX]{x1D702})\equiv (x-y,x+y)/\sqrt{2}$
. The stress field is then given by
$(\unicode[STIX]{x1D719},p,\unicode[STIX]{x1D70F}_{xx},\unicode[STIX]{x1D70F}_{xy})=(-\unicode[STIX]{x03C0}/4,\mathit{Bi}\unicode[STIX]{x03C0}/2,\mathit{Bi},0)$
, with a suitable choice for the background pressure. For the triangle
$T_{-}$
above
$-3/2<x<-1/2$
, the sliplines are instead
$(s_{\unicode[STIX]{x1D6FC}},s_{\unicode[STIX]{x1D6FD}})=(\unicode[STIX]{x1D702},-\unicode[STIX]{x1D709})$
with
$(\unicode[STIX]{x1D719},p,\unicode[STIX]{x1D70F}_{xx},\unicode[STIX]{x1D70F}_{xy})=(\unicode[STIX]{x03C0}/4,-\mathit{Bi}\unicode[STIX]{x03C0}/2,-\mathit{Bi},0)$
.
The fan
$F_{-}$
from
$(-1/2,0)$
contains
$\unicode[STIX]{x1D6FC}-$
sliplines with the form of circular arcs of radius
$r=\sqrt{(x-1/2)^{2}+y^{2}}=\text{const.}$
,
$0<r<\sqrt{2}$
, and
$\unicode[STIX]{x1D6FD}-$
lines that are the radial spokes,
$\unicode[STIX]{x1D703}=\tan ^{-1}[y/(x-1/2)]=\text{const.}$
,
$\unicode[STIX]{x03C0}/4<\unicode[STIX]{x1D703}<3\unicode[STIX]{x03C0}/4$
. Thence,
$(\text{d}s_{\unicode[STIX]{x1D6FC}},\text{d}s_{\unicode[STIX]{x1D6FD}})\equiv (-r\,\text{d}\unicode[STIX]{x1D703},\text{d}r)$
,
$\unicode[STIX]{x1D719}=\unicode[STIX]{x1D703}-\unicode[STIX]{x03C0}/2$
and
$p=-2Bi\unicode[STIX]{x1D719}=Bi(\unicode[STIX]{x03C0}-2\unicode[STIX]{x1D703})$
. Likewise, the fan
$F_{+}$
from
$(1/2,0)$
has circular
$\unicode[STIX]{x1D6FD}-$
lines, with
$(\text{d}s_{\unicode[STIX]{x1D6FC}},\text{d}s_{\unicode[STIX]{x1D6FD}})\equiv -(\text{d}r,r\,\text{d}\unicode[STIX]{x1D703})$
,
$\unicode[STIX]{x1D719}=\unicode[STIX]{x1D703}-\unicode[STIX]{x03C0}$
and
$p=Bi(2\unicode[STIX]{x1D703}-\unicode[STIX]{x03C0})$
.
All the preceding predictions for the stress components match satisfyingly with the low-amplitude numerical solutions when
$Bi\gg 1$
. However, the normal-stress component along
$y=0$
in this construction is
$\unicode[STIX]{x1D70F}_{yy}=-\unicode[STIX]{x1D70F}_{xx}=\mp Bi$
in
$T_{\pm }$
, which is inconsistent with the boundary data: since
$u_{x}=0$
along
$y=0$
, we expect
$\unicode[STIX]{x1D70F}_{yy}=0$
there. Evidently, the constraints of periodicity and no pressure drop or net horizontal force on the swimmer dictate the punch-like stress pattern, but the velocity boundary condition then demands the insertion of a stress discontinuity across
$y=0$
. Along this discontinuity, both
$\unicode[STIX]{x1D70F}_{yy}$
and
$p$
jump by
$2Bi$
across
$y=0$
, thereby enforcing
$\unicode[STIX]{x1D70F}_{yy}=0$
at
$y=0$
.
The slipline solution (figure 7) therefore predicts that the yield surface forms the arc of a circle of radius
$1/\sqrt{2}$
, in agreement with the numerical solutions (figure 2
c and figure 6). Moreover, the leading-order power input (2.10) is found to be
$$\begin{eqnarray}{\mathcal{P}}\sim \int _{-1}^{1}\unicode[STIX]{x03C0}a\cos \unicode[STIX]{x03C0}x\;[p+\unicode[STIX]{x1D70F}_{xx}]_{y=0}\;\text{d}x\rightarrow (2+\unicode[STIX]{x03C0})\unicode[STIX]{x03C0}aBi\int _{-1/2}^{1/2}\cos \unicode[STIX]{x03C0}x\;\text{d}x=2(2+\unicode[STIX]{x03C0})aBi,\end{eqnarray}$$
again in agreement with the numerical data in figure 4(c).
Having constructed the stress field, the velocity over the triangles
$T_{\pm }$
follows directly from the constraint
$\unicode[STIX]{x1D70F}_{xy}=0$
, or
$\unicode[STIX]{x1D713}_{xx}=\unicode[STIX]{x1D713}_{yy}$
. In view of the leading-order boundary condition
$\unicode[STIX]{x1D713}=-a\sin \unicode[STIX]{x03C0}x$
on
$y=0$
, we thus find
$$\begin{eqnarray}\unicode[STIX]{x1D713}=-{\textstyle \frac{1}{2}}a[\sin \unicode[STIX]{x03C0}(x+y)+\sin \unicode[STIX]{x03C0}(x-y)]=-a\sin (\unicode[STIX]{x03C0}x)\cos (\unicode[STIX]{x03C0}y)\quad \text{in }T_{\pm }.\end{eqnarray}$$
Over the fans
$F_{\pm }$
, the velocity field can be found using Geiringer’s equations (2.14). Throughout the fan
$F_{-}$
,
$v_{\unicode[STIX]{x1D6FD}}=0$
and
$\unicode[STIX]{x2202}v_{\unicode[STIX]{x1D6FC}}/\unicode[STIX]{x2202}\unicode[STIX]{x1D703}=0$
, which implies a differential rotation about
$(x,y)=(-1/2,0)$
. Similarly, throughout
$F_{+}$
,
$v_{\unicode[STIX]{x1D6FC}}=\unicode[STIX]{x2202}v_{\unicode[STIX]{x1D6FD}}/\unicode[STIX]{x2202}\unicode[STIX]{x1D703}=0$
, and fluid rotates around
$(x,y)=(1/2,0)$
. Thus, by continuity of angular velocity at the borders of the fans,
$$\begin{eqnarray}\unicode[STIX]{x1D713}=\pm {\textstyle \frac{1}{2}}a[1+\cos (\sqrt{2}\,\unicode[STIX]{x03C0}r_{\pm })]\quad \text{in }F_{\mp },\end{eqnarray}$$
where
$r_{\pm }^{2}=(x\pm 1/2)^{2}+y^{2}$
.
3.3.2 The swimming speed
The leading-order streamfunction calculated above does not predict a swimming speed; locomotion, as in Taylor’s classical problem, is captured only by the contributions at
$O(a^{2})$
. Curiously, the numerical solutions suggest that there is no change in the stress field at
$O(a^{2})$
, which is apparently achieved by suitably shifting the stress discontinuity off the swimmer surface (as highlighted by a curve of zero shear rate in figure 2
b,c). The slipline field then remains unchanged and we need only consider the corrections to the streamfunction of this order. Writing
$\unicode[STIX]{x1D713}=a\unicode[STIX]{x1D713}_{1}+a^{2}\unicode[STIX]{x1D713}_{2}+O(a^{3})$
, with
$\unicode[STIX]{x1D713}_{1}$
the leading-order solution described above, then over the triangles
$T_{\pm }$
we still have
$\unicode[STIX]{x1D713}_{2xx}=\unicode[STIX]{x1D713}_{2yy}$
. However, the boundary conditions (2.16) at this order demand that
$$\begin{eqnarray}\unicode[STIX]{x1D713}_{2}(x,0)+\unicode[STIX]{x1D713}_{1y}(x,0)\sin \unicode[STIX]{x03C0}x=0,\quad \unicode[STIX]{x1D713}_{2y}(x,0)+\unicode[STIX]{x1D713}_{1yy}(x,0)\sin \unicode[STIX]{x03C0}x=0,\end{eqnarray}$$
or
$$\begin{eqnarray}\unicode[STIX]{x1D713}_{2}(x,0)=0,\quad \unicode[STIX]{x1D713}_{2y}(x,0)=-\unicode[STIX]{x03C0}^{2}\sin ^{2}\unicode[STIX]{x03C0}x.\end{eqnarray}$$
Hence,
$$\begin{eqnarray}\unicode[STIX]{x1D713}_{2}={\textstyle \frac{1}{8}}\unicode[STIX]{x03C0}[\sin 2\unicode[STIX]{x03C0}(x+y)-\sin 2\unicode[STIX]{x03C0}(x-y)]-{\textstyle \frac{1}{2}}\unicode[STIX]{x03C0}^{2}y=1/4\unicode[STIX]{x03C0}\sin (2\unicode[STIX]{x03C0}y)\cos (2\unicode[STIX]{x03C0}x)-{\textstyle \frac{1}{2}}\unicode[STIX]{x03C0}^{2}y,\end{eqnarray}$$
in
$T_{\pm }$
. At the tops of the triangles, where
$(x,y)=(n,1/2)$
for
$n=0,\pm 1,\ldots$
, we therefore have the
$O(a^{2})$
corrections to the velocity field,
$$\begin{eqnarray}(u_{2},v_{2})=(\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{2}/\unicode[STIX]{x2202}y,-\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{2}/\unicode[STIX]{x2202}x)=(-\unicode[STIX]{x03C0}^{2},0),\end{eqnarray}$$
which implies a plug, or swimming, speed of
$U=-a^{2}u_{2}=a^{2}\unicode[STIX]{x03C0}^{2}$
. This is exactly twice Taylor’s prediction for Newtonian fluid, as suggested by the numerical data (figure 4
b).
To complete the second-order velocity field, we fill out the solution inside the fans
$F_{\pm }$
: along the fan borders, where
$x=1/2\mp y$
, the correction to the velocity field is
$$\begin{eqnarray}\left(\begin{array}{@{}c@{}}u_{2}\\ v_{2}\end{array}\right)=\left(\begin{array}{@{}c@{}}-\unicode[STIX]{x03C0}^{2}\\ 0\end{array}\right)+1/4\unicode[STIX]{x03C0}^{2}(1-\cos 4\unicode[STIX]{x03C0}y)\left(\begin{array}{@{}c@{}}1\\ \pm 1\end{array}\right),\end{eqnarray}$$
such that, aside from a rigid translation in the
$x$
direction at the swimming speed, the velocity into and out of the fans is again purely angular, with no radial component. Since there is no change to the stress field at this order, Geiringer’s equations (2.14) for
$(u_{2},v_{2})$
remain identical to those at leading order. Hence, the flow throughout the fans combines the rigid translation in
$x$
with differential rotation, and by continuity of angular velocity,
$$\begin{eqnarray}\unicode[STIX]{x1D713}_{2}=-\unicode[STIX]{x03C0}^{2}y+\frac{\unicode[STIX]{x03C0}^{2}}{2\sqrt{2}}r_{\pm }-\frac{\unicode[STIX]{x03C0}}{8}\sin (2\sqrt{2}\,\unicode[STIX]{x03C0}r_{\pm })\quad \text{in }F_{\mp },\end{eqnarray}$$
where, again,
$r_{\pm }^{2}=(x\pm 1/2)^{2}+y^{2}$
. The corrected streamfunction
$\unicode[STIX]{x1D713}=a\unicode[STIX]{x1D713}_{1}+a^{2}\unicode[STIX]{x1D713}_{2}$
is compared with numerical solutions in figure 8.
Streamlines (blue solid) overlying a density map of
$\log _{10}\dot{\unicode[STIX]{x1D6FE}}$
for a full numerical solution with
$Bi=2048$
and
$a=0.025$
. The streamlines are overlain by the prediction
$a\unicode[STIX]{x1D713}_{1}+a^{2}\unicode[STIX]{x1D713}_{2}$
of perfect plasticity theory (white dashed; as calculated from § 3.3.2). (a) Flow in the frame of the swimmer; (b) flow in the frame of the wave.
3.4 High-amplitude waves in the plastic limit (
$a\gg 1$
,
$Bi\gg 1$
)
For high amplitude and yield stress, the flow pattern switches to become dominated by narrow regions of high dissipation adjacent to the swimmer (see figure 3
c,f). In this situation, one can estimate the power input by arguing that the local shear stress
$\unicode[STIX]{x1D70F}_{sn}$
dominates the stress state; i.e.
$\unicode[STIX]{x1D70F}_{sn}\sim Bi$
. Thence, from (2.11),
$$\begin{eqnarray}{\mathcal{P}}\sim Bi\int _{0}^{S}w\;\text{d}s=Bi\int _{-1}^{1}(1+Y_{x}^{2})\;\text{d}x=(2+\unicode[STIX]{x03C0}^{2}a^{2})Bi,\end{eqnarray}$$
where
$w=\sqrt{1+Y_{x}^{2}}$
is the tangential speed at the swimmer surface in the frame of the wave. For
$Bi\gg 1$
, this estimate maximizes the total possible power input (since
$|\unicode[STIX]{x1D70F}_{sn}|\leqslant Bi$
for perfectly plastic flow), and closely reproduces the numerical data at large yield stress, as shown in figure 4(d).
We can generate a prediction of the asymptotic swimming speed
$U=U_{\infty }(a)$
for
$Bi\gg 1$
by comparing the power input (3.15) with the total dissipation rate. In an
$(s,n)$
local coordinate system aligned with the swimmer surface in the frame of the wave, the velocity jump across the boundary layer is
$$\begin{eqnarray}\left[w+\frac{U_{\infty }-1}{\sqrt{1+Y_{x}^{2}}},\frac{(U_{\infty }-1)Y_{x}}{\sqrt{1+Y_{x}^{2}}}\right].\end{eqnarray}$$
Hence, the dissipation rate must be approximately
$$\begin{eqnarray}\langle \dot{\unicode[STIX]{x1D6FE}}(Bi+\dot{\unicode[STIX]{x1D6FE}})\rangle \sim Bi\int _{0}^{S}\left[\left(w+\frac{U_{\infty }-1}{\sqrt{1+Y_{x}^{2}}}\right)^{2}+\frac{4(U_{\infty }-1)^{2}Y_{x}^{2}}{1+Y_{x}^{2}}\right]^{1/2}\text{d}s+\langle \dot{\unicode[STIX]{x1D6FE}}^{2}\rangle .\end{eqnarray}$$
Given that the dissipation must balance (3.15), (3.17) appears to indicate that
$U\rightarrow 1$
in the plastic limit
$Bi\gg 1$
.
Curiously, the data in figure 4(a,b) demonstrate that the swimmer is able to reach speeds in excess of this limit. A simple explanation for this observation is that the power input
${\mathcal{P}}$
in (3.15) arises from an integral over the swimmer surface. On the other hand, the leading-order dissipation rate for
$w\gg (U_{\infty }-1)$
in (3.17) is given more accurately by an integral of the strain rate along the centreline of the boundary layer. As is clear from figure 3(f), the shape of the boundary layer ensures that this centreline is shorter than the perimeter of the swimmer. In other words, the thickening of the boundary layer in the trough of the wave ensures that the power input acts over a longer curve than the effective dissipation rate. To provide a crude estimate of the effect of this foreshortening, we locate the centreline of the boundary layer by
$$\begin{eqnarray}y=Y(x)+\unicode[STIX]{x1D6E5}(x)=\frac{\displaystyle \int _{Y}^{\infty }y\dot{\unicode[STIX]{x1D6FE}}\;\text{d}y}{\displaystyle \int _{Y}^{\infty }\dot{\unicode[STIX]{x1D6FE}}\;\text{d}y},\end{eqnarray}$$
where the function
$\unicode[STIX]{x1D6E5}\ll a$
can be diagnosed from the numerical computations, as illustrated in figure 3(f). Thence, when
$U_{\infty }-1\ll a$
, equating (3.15) and (3.17) gives
$$\begin{eqnarray}U_{\infty }\sim 1-\frac{1}{2}\int _{-1}^{1}\unicode[STIX]{x1D6E5}_{x}Y_{x}\;\text{d}x=1-\frac{\unicode[STIX]{x03C0}^{2}a}{2}\int _{-1}^{1}\unicode[STIX]{x1D6E5}\sin \unicode[STIX]{x03C0}x\;\text{d}x.\end{eqnarray}$$
For the computed curve in figure 3(f), we find
$U_{\infty }\approx 1.51$
, in comparison to the numerical result of
$1.23$
for that value of
$Bi$
and
$a$
. The estimate in (3.19) also predicts the interesting result that
$U_{\infty }$
could continue to increase above unity as
$a$
is increased provided
$\unicode[STIX]{x1D6E5}$
varies more weakly with amplitude than
$\unicode[STIX]{x1D6E5}\sim 1/a$
, which is a trend that is certainly suggested by the data in figure 4(a,b).
A collage of numerical solutions for the longitudinal wave problem, showing density maps of
$\log _{10}\dot{\unicode[STIX]{x1D6FE}}$
overlain by streamlines (blue) in a stationary frame. (a–c)
$a=0.1$
, (d–f)
$a=0.3$
, and (a,d)
$Bi=1/8$
, (b,e)
$Bi=2$
and (c, f)
$Bi=512$
. The inset in (c) shows the thin viscoplastic boundary layer above the swimmer on a magnified scale.
4 Longitudinal waves
For the longitudinal waveform described by (2.5)–(2.6), material points in the surface of the sheet shift horizontally along the swimmer and the vertical surface velocity vanishes (see figure 1
b). With this model, the amplitude of the longitudinal wave cannot exceed
$a=\unicode[STIX]{x03C0}^{-1}$
in order that the sheet does not fold back over itself. We again begin our discussion of the problem by describing observations from the numerical computations, and then follow up by constructing asymptotic solutions to explain these observations.
4.1 Numerical observations
Numerical solutions for different wave amplitudes and yield stresses are shown in figure 9, while a summary of the compiled data from the simulations is given in figure 10. The flow with a longitudinal waveform adopts a qualitatively similar form irrespective of the wave amplitude and comprises a set of rotating cells (figure 9). As the yield stress is raised, the flow develops an increasingly clear structure, comprising rigidly rotating plugs at the centre of each cell, buffered from one another, and from the overlying plug, by regions of weaker shear. A region of higher shear builds up against the swimmer, below the rotating plugs. In the limit
$Bi\gg 1$
, the boundary layer against the swimmer sharpens significantly, and the rotating plugs above become bordered by weakly sheared plastic zones. Note how the centres of the rotating plugs are located slightly above the swimmer, and the plastic zones are composed of circular segments and rectangular constant-stress regions (figure 9
c,f), but the upper part of the overall pattern appears much like the Prandtl punch solution of § 3.3.
Data for longitudinal wave motion, for (a,c,f) different fixed amplitudes
$a=0.025$
,
$a=0.05$
,
$a=0.1$
,
$a=0.2$
and
$a=0.3$
, and for (b,d) different fixed Bingham numbers
$Bi=1/8$
,
$Bi=2$
,
$Bi=32$
,
$Bi=512$
and
$Bi=2048$
. For each plot, the series with the lowest (highest) value of the parameter (
$a$
or
$Bi$
) is marked with stars (open circles) for clarity. (a,b) Swimming speed
$-U$
; (c) scaled power consumption
${\mathcal{P}}/a$
and (d) unscaled power consumption
${\mathcal{P}}$
; (e) a density plot of the efficiency
${\mathcal{E}}=U^{2}/{\mathcal{P}}$
in
$(Bi,a)$
space, together with three contours of constant economy
$U/{\mathcal{P}}$
; and (f) the thickness
$\unicode[STIX]{x1D6FF}$
of the yielded boundary layer above the sheet at
$x=0$
. Short-dashed (blue) lines lines show the theoretical predictions of § 4.2 for
$a\ll 1$
and
$Bi\ll 1$
. The long-dashed (red) lines show the
$Bi\gg 1$
predictions in (4.14) and (4.15) (with
$Bi=2048$
in (d)); the triangles indicate the limits for
$a\ll 1$
.
Unlike for transverse waves, but as for the Newtonian problem (Tuck Reference Tuck1968), the swimming speed (figure 10
a,b) is negative for longitudinal waves. For the physical range of wave amplitudes (
$a<\unicode[STIX]{x03C0}^{-1}$
), the speed increases monotonically from the Newtonian limit (
$Bi\ll 1$
) up to another plastic limit (
$Bi\gg 1$
). Interestingly, with low wave amplitudes
$a\ll 1$
, the swimming speed approaches the same limits as for transverse waves:
$|U|\sim \unicode[STIX]{x03C0}^{2}a^{2}/2$
for
$Bi\ll 1$
and
$|U|\sim \unicode[STIX]{x03C0}^{2}a^{2}$
for
$Bi\gg 1$
(figure 10
b). The power input
${\mathcal{P}}$
again increases with both
$a$
and
$Bi$
(figure 10
c,d), scaling roughly with
$4aBi$
for higher yield stress (which is lower than the corresponding limit
${\mathcal{P}}\sim (2+\unicode[STIX]{x03C0}^{2}a^{2})Bi$
for transverse waves; see figure 4
d). The efficiency
${\mathcal{E}}$
now increases monotonically as
$a$
is increased or as
$Bi$
is decreased (figure 10
e), with higher-amplitude longitudinal strokes in low-yield-stress fluid being more efficient than the optimal transverse waves (by a factor of about two; compare with figure 4
e). The economy
$U/{\mathcal{P}}$
(contours in figure 10
e) shows qualitatively the same behaviour. The thickness of the boundary layer against the swimmer at
$x=0$
, denoted
$\unicode[STIX]{x1D6FF}$
in figure 10(f), decreases with
$Bi$
, and for sufficiently large yield stress scales with
$Bi^{-1/2}$
. These features are rationalized below.
4.2 Low-amplitude solutions in the nearly Newtonian limit (
$a\ll 1$
,
$Bi\ll 1$
)
In the limit
$a\ll 1$
, the boundary conditions (2.17) imply
$$\begin{eqnarray}\unicode[STIX]{x1D713}_{x}(x,0)=0\quad \text{and}\quad \unicode[STIX]{x1D713}_{y}(x,0)\sim a\unicode[STIX]{x03C0}\sin \unicode[STIX]{x03C0}x+{\textstyle \frac{1}{2}}\unicode[STIX]{x03C0}^{2}a^{2}(1+\cos 2\unicode[STIX]{x03C0}x)+O(a^{3}).\end{eqnarray}$$
Hence, following Tuck (Reference Tuck1968) and Blake (Reference Blake1971) for the Newtonian problem (
$Bi=0$
), we find
$$\begin{eqnarray}\unicode[STIX]{x1D713}\sim a\unicode[STIX]{x03C0}y\text{e}^{-\unicode[STIX]{x03C0}y}\cos \unicode[STIX]{x03C0}x+{\textstyle \frac{1}{2}}\unicode[STIX]{x03C0}^{2}a^{2}y(1+\text{e}^{-2\unicode[STIX]{x03C0}y}\cos 2\unicode[STIX]{x03C0}x)+O(a^{3}),\end{eqnarray}$$
which gives the leading-order power consumption and strain rate,
$$\begin{eqnarray}{\mathcal{P}}\sim 2\unicode[STIX]{x03C0}^{3}a^{2},\quad \dot{\unicode[STIX]{x1D6FE}}\sim \dot{\unicode[STIX]{x1D6FE}}_{0}=2a\unicode[STIX]{x03C0}^{2}\text{e}^{-\unicode[STIX]{x03C0}y}|1-\unicode[STIX]{x03C0}y|.\end{eqnarray}$$
This strain rate vanishes at
$y=\unicode[STIX]{x03C0}^{-1}$
, which suggests that plugs may appear in the vicinity of that level when
$Bi>0$
. While such plugs could conceivably interrupt the decay of the Newtonian streamfunction and stress towards a far-field form above this level, we will proceed under the assumption that they do not, and verify this assumption below. The swimming speed predicted by (4.1) is therefore again
$U\sim -\unicode[STIX]{x03C0}^{2}a^{2}/2$
for small yield stress (as confirmed in figure 10
b), and a regular expansion of (2.15) for
$Bi\ll 1$
and
$a\gg Bi\gg a^{2}$
leads to a correction at
$O(Bi)$
to the leading-order solution above:
$$\begin{eqnarray}\unicode[STIX]{x1D713}=a\unicode[STIX]{x03C0}y\text{e}^{-\unicode[STIX]{x03C0}y}\cos \unicode[STIX]{x03C0}x-\unicode[STIX]{x03C0}^{-2}Bi\;G(y)\cos \unicode[STIX]{x03C0}x+O(a^{2}),\end{eqnarray}$$
with
$G(y)$
being the solution to
$$\begin{eqnarray}\left(\frac{\text{d}}{\text{d}y}-\unicode[STIX]{x03C0}^{2}\right)^{2}G=\unicode[STIX]{x03C0}^{4}\text{sgn}(\unicode[STIX]{x03C0}y-1),\quad G(0)=G^{\prime }(0)=0,\quad G\rightarrow 0\text{ as }y\rightarrow \infty .\end{eqnarray}$$
The solution to (4.5) has limiting behaviour
$G(y)\sim 1+O(\text{e}^{-\unicode[STIX]{x03C0}y})$
for
$y\gg 1$
, which, just as for the transverse wave solutions in § 3.2, fails to decay with
$y$
, implying the emergence of another far-field viscoplastic region for
$y\sim y_{B}$
(see figure 11
a). Indeed, if the plugs around the level
$y=\unicode[STIX]{x03C0}^{-1}$
have no significant effect, one expects the same universal solution
$\unicode[STIX]{x1D6F9}(x,y-y_{B})=(\unicode[STIX]{x1D713}+\unicode[STIX]{x03C0}^{2}a^{2}y/2)/Bi$
as in § 3.2 to describe the plugging up of the Newtonian solution, except for a shift of
$1/2$
in the horizontal position (the limiting leading-order streamfunction is proportional to
$\cos \unicode[STIX]{x03C0}x$
rather than the
$\sin \unicode[STIX]{x03C0}x$
in (3.5)). This prediction is confirmed in figure 11(b) which compares shifted yield surfaces and the local streamlines from numerical solutions with
$(a,Bi)\ll 1$
for both types of waves.
(a) The mean distance
$D$
from
$y=0$
to the upper yield surface, across a range of values of
$a$
and
$Bi$
for longitudinal waveforms (black stars), and compared with the same data for transverse waveforms from figure 6 (red squares). The solid line again shows
$D=y_{B}$
from (3.2). (b) The upper yield surface for longitudinal waveforms (black solid) translated by its mean distance
$D$
, compared with the same data for transverse waveforms (red dashed), shifted horizontally by
$1/2$
, for
$Bi=1/8$
and
$a=0.025$
. A few (equally spaced) streamlines in the laboratory frame are also shown for the longitudinal (black) and transverse (red) waves.
Returning now to the zero-strain-rate line at
$y=\unicode[STIX]{x03C0}^{-1}$
, equation (4.4) implies
$$\begin{eqnarray}\dot{\unicode[STIX]{x1D6FE}}^{2}\approx \dot{\unicode[STIX]{x1D6FE}}_{0}\left[\dot{\unicode[STIX]{x1D6FE}}_{0}+Bi(\unicode[STIX]{x1D6E4}+\unicode[STIX]{x1D6F6}\cos 2\unicode[STIX]{x03C0}x)\right]+O(Bi^{2}),\end{eqnarray}$$
where
$$\begin{eqnarray}\unicode[STIX]{x1D6E4}=[2\unicode[STIX]{x03C0}^{-1}G^{\prime }-G-\unicode[STIX]{x03C0}^{-2}G^{\prime \prime }]_{\unicode[STIX]{x03C0}y=1}\quad \text{and}\quad \unicode[STIX]{x1D6F6}=-[2\unicode[STIX]{x03C0}^{-1}G^{\prime }+G+\unicode[STIX]{x03C0}^{-2}G^{\prime \prime }]_{\unicode[STIX]{x03C0}y=1}.\end{eqnarray}$$
This suggests that the curve along which
$\dot{\unicode[STIX]{x1D6FE}}=0$
shifts to the positions,
$$\begin{eqnarray}y=\unicode[STIX]{x03C0}^{-1}\pm \frac{Bi}{a}\text{Max}(0,-\unicode[STIX]{x1D6E4}-\unicode[STIX]{x1D6F6}\cos 2\unicode[STIX]{x03C0}x)\end{eqnarray}$$
(given
$\dot{\unicode[STIX]{x1D6FE}}_{0}>0$
, the curve can only be shifted away from
$y=\unicode[STIX]{x03C0}^{-1}$
where
$\unicode[STIX]{x1D6E4}+\unicode[STIX]{x1D6F6}\cos 2\unicode[STIX]{x03C0}x<0$
). As illustrated in figure 12, equation (4.8) predicts that over two localized windows around the points
$x=\pm 1$
and 0, the leading-order zero-shear-rate line
$y=\unicode[STIX]{x03C0}^{-1}$
splits into two curves. The region enclosed by these two curve is then potentially below the yield stress and may plug up. Of course, over a narrow region surrounding
$y=\unicode[STIX]{x03C0}^{-1}$
, equation (4.4) cannot remain valid as the asymptotic sequence for the strain-rate solution becomes disordered and another, boundary-layer-like solution is needed. Nevertheless, as also shown in figure 12, localized plugs do appear in our numerical solutions in the vicinity of the points
$(\pm 1,\unicode[STIX]{x03C0}^{-1})$
and
$(0,\unicode[STIX]{x03C0}^{-1})$
, with a vertical extent of order
$a^{-1}Bi$
. The shape of the bordering yield surfaces is, however, qualitatively different to the prediction in (4.8).
The shifted positions where
$\dot{\unicode[STIX]{x1D6FE}}=0$
to
$O(Bi^{2})$
. A plug is predicted in the region between the two curves (shaded). (a) Predictions overlain by the plugs that appear in linearized numerical calculations (in which
$x=x_{0}$
on the boundary, such that the swimming speed is zero) for
$Bi=2^{-n}$
,
$n=3,4,5,6$
(blue, red, green and black dashed respectively). (b) Predictions overlain by the plugs that appear in full numerical calculations for
$Bi=1/8$
, for
$a=0.0125$
(blue),
$a=0.025$
(red) and
$a=0.05$
(green).
Note that the perturbed streamfunction connects directly across the leading-order zero-strain-rate line outside the windows where plugs form for both (4.4) and in the numerical solutions. Thus, by a process of continuation in
$x$
and
$y$
, one can circumnavigate the plugs that do form, implying that the leading-order flow is not interrupted by a rigid layer spanning the domain. This observation justifies the use of (4.1) and (4.4) both below and above the level
$y=\unicode[STIX]{x03C0}^{-1}$
, and thus the prediction of the far-field swimming speed
$U\sim -\unicode[STIX]{x03C0}^{2}a^{2}/2$
.
4.3 Viscoplastic boundary-layer solutions for
$Bi\gg 1$
Irrespective of the amplitude
$a$
, in the limit
$Bi\gg 1$
the numerical solutions develop narrow boundary layers against the swimmer of thickness
$\unicode[STIX]{x1D716}\equiv Bi^{-1/2}$
(figures 9
c,f and 10
f). Over these layers, the velocity field is largely horizontal and
$$\begin{eqnarray}[u,v]\rightarrow [{\mathcal{U}}(x,\unicode[STIX]{x1D701})-U,\unicode[STIX]{x1D716}V(x,\unicode[STIX]{x1D701})],\quad \text{with }\unicode[STIX]{x1D701}=\unicode[STIX]{x1D716}^{-1}y,\end{eqnarray}$$
where
${\mathcal{U}}$
denotes the horizontal boundary-layer speed in the stationary laboratory frame, which jumps from
$$\begin{eqnarray}{\mathcal{U}}(x,0)\equiv U+X_{t}=U+a\unicode[STIX]{x03C0}\cos \unicode[STIX]{x03C0}(x_{0}+t),\end{eqnarray}$$
to
$O(\unicode[STIX]{x1D716})$
at the edge of the boundary layer, which we denote by
$y=\unicode[STIX]{x1D716}\unicode[STIX]{x1D701}=\unicode[STIX]{x1D716}Z(x)$
. In addition, the shear stress dominates the stress state and the combination of force balance and the constitutive equation demand
$$\begin{eqnarray}\unicode[STIX]{x1D70F}_{xy}\sim Bi\;\text{sgn}({\mathcal{U}}_{\unicode[STIX]{x1D701}})+\unicode[STIX]{x1D716}^{-1}{\mathcal{U}}_{\unicode[STIX]{x1D701}},\quad P_{x}\sim \unicode[STIX]{x1D716}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}\unicode[STIX]{x1D70F}_{xy}\sim {\mathcal{U}}_{\unicode[STIX]{x1D701}\unicode[STIX]{x1D701}},\quad P_{\unicode[STIX]{x1D701}}\sim 0,\end{eqnarray}$$
where
$P=Bi^{-1}p$
. Above the boundary layer, the strain rate must fall sufficiently to allow a match to the overlying rigidly rotating plugs. Hence,
${\mathcal{U}}_{\unicode[STIX]{x1D701}}=O(\unicode[STIX]{x1D716})$
at
$\unicode[STIX]{x1D701}=Z$
. This match also demands that
$p=O(Bi)$
in view of the Riemann invariants of the sliplines within the plastic regions, which guides the preceding choice of scaling for
$P$
(and rationalizes the observed boundary-layer width
$\unicode[STIX]{x1D716}=Bi^{-1/2}$
). The boundary-layer solution is therefore
$$\begin{eqnarray}{\mathcal{U}}=\left(1-\frac{\unicode[STIX]{x1D701}}{Z}\right)^{2}[U+a\unicode[STIX]{x03C0}\cos \unicode[STIX]{x03C0}(x_{0}+t)],\quad \text{with }[U+a\unicode[STIX]{x03C0}\cos \unicode[STIX]{x03C0}(x_{0}+t)]=\frac{1}{2}P_{x}Z^{2},\end{eqnarray}$$
which compares well with numerical solutions, as illustrated in figure 13(a).
Boundary-layer horizontal velocity profiles in a stationary frame for
$a=0.3$
and
$Bi=512$
. Profiles, scaled by the velocity of the wave
${\mathcal{U}}(x,0)$
, are taken at
$x=0$
,
$0.1$
,
$0.2$
,
$0.3$
(black) and
$x=0.7$
,
$0.8$
,
$0.9$
,
$1$
(blue) and compared with the theoretical prediction in (4.12) (red circles). (a) A large-scale plot of the boundary layer against the swimmer; (b) the same data on a larger scale, showing the gently rotating plugs above the boundary layer (where the scalings do not collapse all the curves).
From (4.12) we may deduce the leading-order shear stress on the swimmer,
$$\begin{eqnarray}\unicode[STIX]{x1D70F}_{xy}(x,0)=-Bi\;\text{sgn}[U+a\unicode[STIX]{x03C0}\cos \unicode[STIX]{x03C0}(x_{0}+t)],\end{eqnarray}$$
which enables us to calculate the swimming speed, given that there is no net force exerted by this stress:
$$\begin{eqnarray}0=\int _{-1}^{1}\text{sgn}[U+a\unicode[STIX]{x03C0}\cos \unicode[STIX]{x03C0}(x_{0}+t)]\;\text{d}x=\int _{-1}^{1}(1+\unicode[STIX]{x03C0}a\cos \unicode[STIX]{x03C0}z)\;\text{sgn}[U+a\unicode[STIX]{x03C0}\cos \unicode[STIX]{x03C0}z]\;\text{d}z.\end{eqnarray}$$
In other words, the swimming speed is selected so that there are equal areas of positive and negative shear stress, or equivalently translation speed, along the swimmer (with the symmetry of the problem placing the sign switches at
$x=\pm 1/2$
). A short calculation with (4.14) shows that
$U\sim -\unicode[STIX]{x03C0}^{2}a^{2}$
for
$a\ll 1$
, while the speed for larger wave amplitudes can be calculated numerically from (4.14) and is plotted, together with numerical data, in figure 10(b). The leading-order power input also now follows as
$$\begin{eqnarray}{\mathcal{P}}=Bi\int _{-1}^{1}|U+a\unicode[STIX]{x03C0}\cos \unicode[STIX]{x03C0}z|\;(1+a\unicode[STIX]{x03C0}\cos \unicode[STIX]{x03C0}z)\;\text{d}z,\end{eqnarray}$$
equivalent to the leading-order plastic dissipation over the boundary layer, in accordance with the constraint (2.10). The power input has the limit
${\mathcal{P}}\sim 4aBi$
for
$a\ll 1$
, and is again shown for higher
$a$
and compared with numerical data in figure 10(d).
Finally, mass conservation across the boundary layer demands that
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\int _{0}^{Z}{\mathcal{U}}\,\text{d}\unicode[STIX]{x1D701}\equiv \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left[1/3Z(U+a\unicode[STIX]{x03C0}\cos \unicode[STIX]{x03C0}(x_{0}+t))\right]=-V[x,Z(x)],\end{eqnarray}$$
with
$V(x,Z)$
specified by the match to the two overlying rigidly rotating plugs and intervening plastic regions. Unfortunately, without solving for the detailed flow pattern above the boundary layer, we cannot prescribe
$V(x,Z)$
, and so the integral of (4.16) simply relates the boundary-layer thickness
$Z$
to an unknown flux function. Although a diagnosis of both
$Z$
and
$V(x,Z)$
from the numerical computations is consistent with (4.16), we continue no further with the analysis.
5 Concluding remarks
In this paper we have considered the locomotion of a flexible sheet in a yield-stress fluid driven by either transverse retrograde waves (i.e. waves travelling opposite to the direction of motion) or longitudinal prograde waves (so the motion is in the wave direction). Both swimming problems are classical and popular models for locomotion through Newtonian fluid, and our focus here was to examine the effect of a yield stress in each case. For the task, we provided combined numerical and asymptotic analyses, treating the ambient as a Bingham fluid.
For both types of wave motion, the swimming speed increases with the dimensionless yield stress (except for the highest-amplitude transverse waves), reflecting how the swimmer gains increased traction on the ambient fluid due to its rheology. While the gain in speed is at most by a factor of order unity over that expected for Newtonian fluid, the power expenditure required for locomotion increases steadily with yield stress. We also calculated a measure of the efficiency of each swimming strategy, following Blake (Reference Blake1971), allowing a quantitative comparison for different wave amplitude and yield stress. For both strokes, the steady power increase for high yield-stress fluid always leads to inefficient swimming, with transverse waves being slightly less efficient owing to a higher power expenditure. Interestingly, however, the gain in speed from the fluid rheology for transverse wave motion initially overcomes the increased power expenditure, leading to a optimal finite yield stress for locomotion. Thus, there may be potential for an organism to tune its swimming strategy to its rheological environment, as proposed previously for locomotion in viscoelastic fluid. Such an optimum is not present in longitudinal wave motion, where high-amplitude Newtonian swimming is the most efficient. Another measure of the effectiveness of the locomotion strategy, the swimming speed per unit power (the economy of Krieger et al. Reference Krieger, Spagnolie and Powers2014), shows qualitatively the same behaviour, although swimming is optimized for transverse waves at low amplitude in Newtonian fluid.
In the limit of high yield stress (or slow wave speed;
$Bi\gg 1$
) and low wave amplitude, transverse wave motions give rise to an ideal plastic flow that we have constructed using the sliplines (characteristics of the stress field) of plasticity theory. The problem resembles Prandtl’s classical analysis of the indentation of a punch into a plastic half-space, although, curiously, this solution is not what one would expect in view of the boundary conditions against the swimmer. In fact, the boundary conditions are reconciled by inserting a stress discontinuity along the boundary. Although this feature is presumably dictated by the global constraint of force balance, we have no satisfying explanation other than observations based on the numerical computations. For longitudinal prograde waves, the high-yield-stress limit can instead be analysed using viscoplastic boundary-layer theory (Balmforth et al.
Reference Balmforth, Craster, Hewitt, Hormozi and Maleki2017). Interestingly, for both waveforms and despite the solutions having rather different constructions, the swimmer is exactly twice as fast in this limit as it would be in a Newtonian fluid.
Another interesting feature of our analysis concerns the manner in which the introduction of a small yield stress modifies the Newtonian flow pattern. In the flow problems we have considered, the stresses exerted on the fluid decay away from the moving boundary, which raises the question of how, and where, the yield stress must be included to correct the leading-order Newtonian solution and allow the fluid to plug up. The mathematical fashion in which this occurs furnishes an exercise in matched asymptotics, some details of which we have exposed here.
The differences between the two swimming strokes naturally raise the question of whether a more general wave motion might be more advantageous, such as a mix of longitudinal and transverse flexures or non-sinusoidal waveforms. Furthermore, while we have fixed the swimming stroke in this study, in an actual physical situation the swimmer will instead impose a force, which has been shown to have a significant impact for locomotion above a thin viscoplastic film (Pegler & Balmforth Reference Pegler and Balmforth2013). Of course, the simplistic geometry of the two-dimensional flexible sheet is itself problematic for an immediate application of these results to real biological organisms, and an obvious next step is to consider the effect of a yield stress on the propulsion of helical waves travelling down a cylindrical flexible filament (e.g. Lauga & Powers Reference Lauga and Powers2009). We leave such considerations for future work.
Appendix. Numerical details
We solve (2.15) numerically using an augmented Lagrangian (AL) scheme (Fortin & Glowinski Reference Fortin and Glowinski2006), which is widely used to compute viscoplastic flow (e.g. Vinay, Wachs & Agassant Reference Vinay, Wachs and Agassant2005). The scheme solves the nonlinear partial differential equation (2.15) iteratively by the introduction of ‘dummy’ strain rate and plastic stress tensors, with each step in the iteration process requiring only the solution of a linear fourth-order elliptic equation and a nonlinear algebraic expression.
For transverse waves, the lower boundary at
$y=Y(x)$
is not flat. We deal with this complication by simulating a domain which does extends down to a flat surface, including the region occupied by the swimmer. We impose the desired velocity field over this region by setting
$\check{\unicode[STIX]{x1D713}}=\unicode[STIX]{x1D713}+a\sin \unicode[STIX]{x03C0}x$
and then forcing
$\check{\unicode[STIX]{x1D713}}$
to satisfy the constitutive law. By artificially increasing the yield stress below the surface of the swimmer,
$y\leqslant Y(x)$
, we can then engineer
$\check{\unicode[STIX]{x1D713}}=\check{\unicode[STIX]{x1D713}}_{y}=0$
over that region, thereby enforcing the correct boundary conditions at the swimmer surface. Note that some care is required with the procedure, since it adds an unwanted contribution to the magnitude of the stress at each stage of the iteration which must be carefully accounted for.
With the domain expanded to a rectangular region, the equations are solved using a Fourier transform in the
$x$
-direction and second-order finite difference discretization in the
$y$
-direction (leading to a standard band-diagonal matrix inversion). In most simulations we used grid spacings
$(\unicode[STIX]{x0394}x,\unicode[STIX]{x0394}y)$
in the
$(x,y)$
directions of approximately
$(\unicode[STIX]{x0394}x,\unicode[STIX]{x0394}y)=(2^{-10},2^{-10})$
, although we confirmed that solutions were well resolved by comparison with a few calculations at a higher resolution. In general, simulations with higher
$a$
and higher
$Bi$
required a higher resolution. In particular, the spectral method in the
$x$
-direction required a large number of grid points to reduce inaccuracies associated with Gibbs phenomena at plug boundaries. We deemed solutions to have converged in the iterative scheme once the error in the average magnitude of the dummy strain-rate tensor had fallen below
$10^{-7}$
, although in some cases we confirmed that solutions with a much smaller tolerance gave indistinguishable solutions. We used an iteration parameter in the AL scheme of
$r=\max (1,3Bi/4)$
.
We iterate the augmented Lagrangian (AL) scheme to convergence at each step of another iteration to find the constant
$K$
associated with a vanishing pressure drop (see § 2.3). At each step of this second iteration, we use the converged solution from the previous step as an initial guess for the new AL iteration. This strategy significantly reduces the number of steps required for convergence of the AL scheme, and thus for the whole nested iteration: initial AL iterations take several thousand steps (and increasingly many for larger
$a$
and
$Bi$
), whereas once
$K$
is close to convergence, AL iterations take several hundred steps or fewer.
