A numerical investigation of an Immersed Boundary (IB) model of an effectively inextensible, finite swimmer in a Stokesian Oldroyd-B flow is presented. The swimmer model is a two-dimensional sheet of finite extent and its gait is generated by an elastic force which penalizes deviations from a target shape. A non-stiff IB method is employed to remove the impeding time step limitation induced by strong tangential forces on the swimmer. It is found that for a swimmer with a prescribed gait its mean propulsion speed decreases with increasing Deborah number De toward an apparent asymptotic minimal value. However, as the swimmer is allowed to deviate more from the target shape, the monotonic locomotion behavior with De is broken. For a sufficiently flexible swimmer, viscoelasticity can enhance locomotion but the swimmer in the viscoelastic fluid always remains slower than when it is propelling in a Newtonian fluid. Remarkably, the addition of viscoelastic stress diffusion dramatically alters the swimmer propulsion and can lead to a speed-up over the swimmer in the Newtonian fluid.

For spatially periodic swimmers with a prescribed sinusoidal beating pattern in a non-Newtonian Stokesian (zero-Reynolds number) flow, analytic asymptotic results indicate that viscoelasticity always hinders locomotion; there is a monotonic decrease of the mean propulsion of the swimmer as the Deborah number De (the ratio of the viscoelastic relaxation time and the characteristic time scale of the swimming motion) increases.9,10,13,14 On the other hand, experiments7,16,20 suggest that for a swimmer of a finite extent, the propulsion seems to be sensitive to the details of the motion (undulatory versus helical, flexible versus fixed gait, small amplitude versus large amplitude, etc.). For example, Shen and Arratia20 found that for C. elegans nematodes in shear thinning polymeric fluids, viscoelasticity hinders locomotion while the experiments of Liu, Powers, and Breuer16 of a helical swimmer in a Boger fluid (a viscoelastic fluid with a shear-rate independent viscosity) showed that the swimmer’s propulsion could be increased by viscoelastic response. More recently, the experiments of Espinosa-Garcia, Lauga, and Zenit7 of flexible swimmers with a magnetically actuated head suggest that viscoelasticity always enhances locomotion. There is also a numerical study of an immersed boundary (IB) model of a finite swimmer in a Stokesian Oldroyd B flow by Teran, Fauci, and Shelley25 which indicates that viscoelasticity response can lead to a speed-up of the swimmer with respect to the swimmer in the Newtonian flow. In a subsequent, related numerical investigation, Thomases and Guy26 argue that the combination of gait asymmetry and the ability of the swimmer to deviate from a prescribed gait to respond to viscoelastic stresses are the main mechanisms that contribute to such a speed-up. Refs. 6 and 23 provide recent theoretical and experimental overviews of this locomotion problem.

Here, inspired by the study of Teran, Fauci, and Shelley25 and that of Thomases and Guy,26 we conduct a detailed numerical investigation of the same underlying immersed boundary model of a finite swimmer, exploring a wide range of the model parameters, and examining the important effects of stress diffusion. Another model of immersed boundary type has been proposed by Shirgaonkar, MacIver, and Patankar.21 

We find that viscoelastic stress diffusion, either numerically or explicitly added into the model, is a decisive factor for a speed-up of the swimmer in the viscoelastic fluid over that in the Newtonian one. In the absence of viscoelastic stress diffusion, our study shows that for a swimmer with a prescribed gait, its mean propulsion speed always decreases with increasing De toward an apparent asymptotic minimal value, in agreement with the asymptotic results for spatially periodic small amplitude swimmers.9,10,13,14 However, as the swimmer is allowed to deviate more from the target shape, this monotonic locomotion behavior with De is broken. For a sufficiently flexible swimmer, with large deviations from the target shape, viscoelasticity can enhance locomotion but the swimmer in the viscoelastic fluid always remains slower than that in the Newtonian fluid. Remarkably, the addition of viscoelastic stress diffusion enhances drastically the propulsion of the finite swimmer and can indeed lead to a speed-up. We argue that this peculiar effect is connected to the highly concentrated, singular forces at the swimmer’s end points for this particular model. At the same time, our study highlights the importance of diffusive transport for zero Reynolds number swimming.

The swimmer is modeled as an immersed finite sheet whose inextensibility and gait generation are controlled by appropriate interfacial forces obtained from a discrete elastic and tensile energy. The coupling of the swimmer and the flow is done within the Immersed Boundary (IB) framework.8,18 We employ an efficient, semi-implicit IB method2,4 that overcomes the severe time step restriction induced by the enormous interfacial forces needed to generate the swimmer’s gait and to accurately enforce the inextensibility constraint. We present a systematic numerical study to examine the effect of the various parameters on the model’s solution. This study identifies that high resolution is critical for this problem and that a low accuracy treatment of the viscoelastic stress advection term can lead to a striking change in the simulated locomotion.

The rest of this article is organized as follows. The model is fully described in Section II and the numerical methodology is reviewed in Section III. A numerical validation against analytic, asymptotic results is given in Section IV and a convergence study is given in the  Appendix. Section V is devoted to the numerical results. Finally, some concluding remarks are given in Section VI.

We model the swimmer as an inextensible sheet which undergoes imposed lateral waves in a two-dimensional, incompressible, Stokesian flow. The lateral waving motion gives rise to stresses that cause the forward propulsion of the swimmer. This model goes back to a seminal paper of Taylor24 on the swimming of microscopic organisms. Our study focuses on the particular case in which the sheet is of finite extent (i.e., spatially non-periodic) and the viscoelastic fluid is modeled by an Oldroyd B constitutive equation with a single relaxation time.

An Oldroyd B fluid is a dilute solution of elastic polymer chains modeled as Hookean dumbbells in a Newtonian solvent. The equations governing the flow can be written as1,12

p=ηs2u+τp+F,
(1)
u=0,
(2)

in a domain Ω, where p is the pressure, ηs is the viscosity of the solvent, u is the velocity field, F is the force on the sheet resulting from the gait generation mechanism (described below) and the inextensibility constraint, and τp the polymer stress given by Kramers formula1 

τp=nkBTI+nHQQ.
(3)

In (3), n is the number density of polymer chains, kB is the Boltzmann constant, T is the temperature, and H is the spring constant. The tensor QQ is the second moment of the probability density function for finding a dumbbell with end-to-end displacement Q. The polymeric stress τp satisfies the constitutive equation1 

τp+λHτp=ηpu+uT,
(4)

where λH is the relaxation time constant of the Hookean dumbbells, ηp = nkBH the polymeric viscosity, and

τp=τpt+uτpuτp+τpuT
(5)

is the upper convective derivative of τp (with the last two terms understood as matrix products). Substituting (3) into (4) and defining the dimensionless stress

Sp=(kBT/H)1QQ
(6)

(kBT/H is proportional to the polymer chain’s radius of gyration), we get a simplified constitutive equation12 

λHSp=ISp.
(7)

We select the characteristic flow time scale to be the inverse of the beating pattern frequency ω−1 and the characteristic length scale is its wavelength. The Deborah number is therefore De = λHω and the flow equations in dimensionless form can be written as

p=2u+ηpηs1DeSp+F,
(8)
u=0,
(9)
Sp=1DeISp,
(10)

where p and F have been redefined appropriately. Throughout this study, we take ηpηs1De=1/2 and keep this value fixed, as done in Ref. 25. As the Deborah number is varied, the viscosity ratio is adjusted appropriately to maintain this constant.

For future reference when describing the numerical methodology, we write (10) as

Spt=G(u,Sp),
(11)

where

G(u,Sp)=uSp+uSp+SpuT+1De(ISp).
(12)

The swimmer (sheet) Γ is represented in parametric form as X(s, t), where s is a Lagrangian variable and t stands for time. We use the IB framework18 to describe the interaction of the flow and the swimmer. The interfacial force F is written as a delta distribution of a force density f(X(s, t)) by

F(x,t)=Γf(X(s,t))δxX(s,t)ds
(13)

and the kinematic condition expressing the continuity of velocity at the swimmer is also expressed in terms of a delta distribution,

dX(s,t)dt=Ωu(x,t)δxX(s,t)dx.
(14)

The swimmer’s gait is generated by penalizing deviations from a prescribed, time periodic target κ̄(s,t). We follow the approach proposed by Fauci and Peskin8 to generate the waving motion and to enforce the inextensibility condition through a density force given by

f(X(s,t))=δE[X(s,t)]δX,
(15)

where E[X] is the elastic energy

E[X]=S12ΓX(s,t)s12ds+S22Γκ(s,t)κ̄(s,t)2ds.
(16)

Here κ is the sheet’s mean curvature. We take S1 large enough (as specified in Sections IV and V) to accurately enforce the swimmer inextensibility and we vary S2 over a wide range to allow for small to large deviations from the target curvature.

Some authors have proposed linking the model parameter S2 to the flexural or bending rigidity of some type of micro-swimmers.15,17 This is an interesting idea but there are some important caveats in this association. First, S2 is needed to produce the gait. If reduced too much, no net propulsion will be observed. In contrast, in a real swimmer the flexural rigidity measures the swimmer’s elastic response to forces but there is a separate mechanism to generate the gait. Second, dimensionally S2 does not match the units of the flexural rigidity. Third, direct measurements of the flagellar bending rigidity for human sperm are not available,11 there is a wide range (several orders of magnitude) variation of this elastic physiological coefficient among mammalians, and it is not constant; it varies greatly (about one order of magnitude) with arc-length.11 Nevertheless, if we take the cross-sectional area to match dimensionally S2 with the flexural rigidity we obtain that S2 = 10 corresponds to a flexural rigidity of 10−20 N m2, which is within the range of that elastic parameter for mammalian sperm.

Henceforth we will refer to S2 as the bending stiffness parameter of the model.

The computational domain is the square Ω = [0, C] × [0, C], which is discretized with a uniform Cartesian grid Gh with N × N computational cells of size h × h, where h = C/N. Periodic boundary conditions are assumed for the flow variables and Sp. The swimmer Γ is of extent L and is discretized with a uniform grid Ghs in s, consisting of the points sj = jhs, j = 0, …, Nb with hs = L/Nb.

Let Dh and Lh be the standard centered second-order finite difference operators to approximate the gradient and Laplacian, respectively. Given the excess polymeric stress Spn at the current time tn, we discretize the Stokes equations and the swimmer’s evolution equation as

Dhpn+1Lhun+1=ηpηs1DeDhSpn+Snfhs(Xn+1),
(17)
Dhun+1=0,
(18)
Xn+1XnΔt=Snun+1,
(19)

where Δt is the time step size. Sn and Sn are the IB Method spreading and interpolation operators,

(Sng)(x)=sjGhsg(sj)δh(xXn(sj))hs,
(20)
(Snw)(s)=xijGhw(xij)δh(xijXn(s))h2,
(21)

where δh(x)=φxφy, with x = (x, y) and

φ(r)=14h1+cosπr2h,r2h,0,r>2h
(22)

is an approximation to the 2D delta distribution.18 

To compute the force density fhs(X), we employ the following discretization of the energy:

Ehs[X]=S12j=0NbD̄hsXj12hs+S22j=1Nb1κjκ̄j2hs,
(23)

where D̄hs stands for the centered finite difference for j = 1…Nb − 1 and for the forward and backward finite differences for j = 0 and j = Nb, respectively. denotes the trapezoidal rule sum with a factor of 1/2 for the end points. The mean curvature κj is computed from the swimmer configuration as

κj=Dhs2YjDhsXjDhs2XjDhsYj(DhsXj)2+(DhsXj)23/2,j=1,,Nb1,
(24)

where Dhs and Dhs2 are the second order, centered, finite difference approximations to the first and second derivative, respectively.

The elastic force density is given by

fhs(X)k=EhsXk,EhsYk,k=0,1,,Nb.
(25)

This discretization satisfies kfhs(X)k=0, which is a compatibility condition for Stokes equations.

Equations (17)-(19) are an implicit system for the future velocity and sheet configuration, un+1 and Xn+1. To solve this system we use the methodology developed in Refs. 2 and 4, which we briefly outline next.

The solution to the discrete Stokes equations (17) and (18) can be written as

un+1=LhSnfhs(Xn+1)+ηpηs1DeDhSpn,
(26)

where Lh is the fluid solver operator whose application produces the velocity field given the force density and the excess polymeric stress. We compute Lh using the Fast Fourier Transform (FFT). Substituting this into (19), we obtain a system of equations for Xn+1 only,

Xn+1=Mnfhs(Xn+1)+bn,
(27)

where Mn the flow-structure operator is given by Mn=ΔtSnLhSn and

bn=Xn+ΔtDeηpηsSnLhDhSpn.
(28)

The non-linear system (27) is solved efficiently with Newton’s method combined with an expedited evaluation of the matrix representation of Mn as detailed in Ref. 4. Once Xn+1 is computed, un+1 is obtained from (26) and the excess stress is updated. This is done using un+1 in an explicit, second order, total variation diminishing (TVD) Runge-Kutta method22 

Sp=Spn+ΔtGh(un+1,Spn),
(29)
Sp=Sp+ΔtGh(un+1,Sp),
(30)
Spn+1=Sp+Sp2.
(31)

Here Gh is the discrete version of (12) obtained by using standard second order spatial finite differences for all the terms except for the advection u ⋅ ∇Sp, which is computed with third order ENO (Essentially Non-Oscillatory) scheme,22 unless otherwise stated. As noted in Ref. 3, the ENO, which does upwinding plus a high order correction with a sliding stencil, provides a localized high order regularization which is effective for capturing thin and potentially singular stress structures in the Oldroyd B model.

While the numerical methodology is not new,2–4 it has not been applied to the type of interfacial forces arising in the model under consideration. In addition to performing extensive resolution and convergence studies, we carry out as a validation test a comparison with the asymptotic, small amplitude result of Lauga.13 He showed that for a swimmer of the form y(x, t) = asin(kx + ωt), the ratio of an Oldroyd B swimmer (mean) speed UV to the Newtonian fluid swimmer speed (Taylor24)

U=12ωk(ak)2+O(ak)4
(32)

satisfies

UVU=1+ηsηs+ηpDe21+De2+O(ak)2.
(33)

We perform simulations for various resolutions up to the time equal to five wave periods. All the parameters used for this test are provided in Table I with ηp/ηs = 1/2 and De = 1.

TABLE I.

Parameters for the spatially periodic swimmer.

Domain Ω [0, 1] × [0, 1] 
Mesh size h = 1/N N = 64, 128, 256, 512 
Lagrangian mesh size hs 0.5h 
Time step Δt hs 
Inextensibility parameter S1 106 
Bending stiffness parameter S2 104 
Wave number k 2π 
Frequency ω 2π 
Amplitude a 0.001, 0.003, 0.006, 0.01, 0.02, 0.03, 0.04, 0.05 
Domain Ω [0, 1] × [0, 1] 
Mesh size h = 1/N N = 64, 128, 256, 512 
Lagrangian mesh size hs 0.5h 
Time step Δt hs 
Inextensibility parameter S1 106 
Bending stiffness parameter S2 104 
Wave number k 2π 
Frequency ω 2π 
Amplitude a 0.001, 0.003, 0.006, 0.01, 0.02, 0.03, 0.04, 0.05 

Figure 1 shows that the numerical results follow closely the parabolic growth of Lauga’s asymptotic projection and there is excellent agreement for amplitudes of up to 2% of the wavelength. For larger amplitudes, there is an evident gap between the computed speed and the corresponding asymptotic value. This gap is of O(ak2), consistent with the asymptotic approximation (33), even though higher resolution would be needed to obtain a more accurate speed for the largest amplitude values (0.04 and 0.05).

FIG. 1.

Numerical validation. Comparison of the computed swim speed for different numerical resolutions with Lauga’s asymptotic projection (continuous curve), De = 1 and ηp/ηs = 1/2.

FIG. 1.

Numerical validation. Comparison of the computed swim speed for different numerical resolutions with Lauga’s asymptotic projection (continuous curve), De = 1 and ηp/ηs = 1/2.

Close modal

We have performed similar simulations for a Newtonian fluid and they agree well with Taylor’s formula (32).

We now consider a finite swimmer in an Stokesian Oldroyd-B fluid and two asymmetric swimming motions. One in which the gait amplitude increases from head to tail and the target curvature is given by

κ̄T(s,t)=(1s)Ak2sin(ks+ωt),s[0,L],
(34)

this is the same wave form as that in Ref. 25, and the other in which the gait amplitude increases from tail to head and the target curvature is

κ̄H(s,t)=(1L+s)Ak2sin(ks+ωt),s[0,L].
(35)

In (34) and (35), A is the maximum amplitude of the swimmer’s gait and L is the swimmer’s length. Note that s = 0 and s = L correspond to the tail and the head, respectively. The wave direction is the same for both beating patterns. All of our numerical studies were performed for both swimmers but we did not find any substantial, qualitative difference in the results. Therefore, we report only on the swimmer with target curvature (34).

To obtain the swimmer’s initial configuration (x0(s), y0(s)) we solve the system,

θ0s(s)=κ̄(s,0),forκ̄=κ̄Tandκ̄=κ̄H,
(36)
x0s(s)=cos(θ0(s)),
(37)
y0s(s)=sin(θ0(s)).
(38)

To do this, we first integrate analytically (36) to get the tangent angle θ0 (we set the constant of integration to zero). We then approximate

x0(sj)=0sjcos(θ0(t))dt,
(39)
y0(sj)=0sjsin(θ0(t))dt,
(40)

by the composite trapezoidal rule quadrature for j = 0, …, Nb. Here Nb = [L/hs] and sj = jhs.

We fix the length of the swimmer to L = 0.6, following Ref. 25, and the amplitude is chosen to be A = 0.1. This is the largest value for which we could afford an accurate numerical resolution of the coupled viscoelastic flow-swimmer system. The model and numerical parameters are presented in Table II.

TABLE II.

Parameters for a finite swimmer in Stokesian Oldroyd-B flow.

Lagrangian mesh size hs 0.5h 
Time step Δt 0.5hs 
Inextensibility parameter S1 As specified in each simulation 
Bending stiffness parameter S2 As specified in each simulation 
Swimmer’s length L 0.6 
Wave number k 2πL 
Frequency ω 2π 
Final Time 20 
Amplitude A 0.1 
Viscosity ratio polymer/solvent ηp/ηs 12 
Lagrangian mesh size hs 0.5h 
Time step Δt 0.5hs 
Inextensibility parameter S1 As specified in each simulation 
Bending stiffness parameter S2 As specified in each simulation 
Swimmer’s length L 0.6 
Wave number k 2πL 
Frequency ω 2π 
Final Time 20 
Amplitude A 0.1 
Viscosity ratio polymer/solvent ηp/ηs 12 

Given that this is a finite swimmer and we employ periodic boundary conditions as in Refs. 25 and 26, we have conducted a study to examine the effects of the domain size. This is presented in Subsection 1 of the Appendix. Based on this study we choose [0, 4] × [0, 4] as an adequate computational domain, with negligible domain size effect (0.5% variation in the mean propulsion speed with respect to that using [0, 8] × [0, 8]), to perform all the simulations. The mean propulsion speed is calculated by averaging the horizontal (x) component of the center of mass of the swimmer over one time period.

We have also performed a convergence and resolution study which is summarized in Subsection 2 of the Appendix. This study shows that significantly high resolution is needed to accurately capture the enormous singularly supported forces on the swimmer as well as small scale viscoelastic structures.

One of the key parameters in the model is the bending stiffness constant S2, which penalizes deviations from the target curvature. Effectively, for sufficiently large S2, the swimmer will follow closely the given gait generated by the time-periodic target curvature. As S2 is reduced, larger deviations from the target curvature are allowed in response to the flow. Thus, large and small values of S2 yield two very distinctive types of swimming dynamics, corresponding to a “stiff” and a “soft” swimmer, respectively.26 

We examine the propulsion of the swimmer as the bending stiffness parameter S2 is varied over wide range of values, spanning four orders of magnitude, while keeping all other parameters fixed (as in Table II with a spatial resolution of h = 1/256).

Figure 2 shows the swimmer’s mean x location as a function of time for various values of S2 for (a) a Newtonian flow (De = 0) and (b) a viscoelastic flow with De = 5. The propulsion is drastically affected by the value of S2. The swimmer with S2 = 104 achieves a displacement in excess of 50% larger than that of the swimmer with S2 = 10. For S2 = 2, the swimmer propulsion is minimal, with a final displacement after 20 periods of less than 17% its length.

FIG. 2.

Swimmer’s mean x position for values of the bending stiffness parameter S2 = 2, 101, 102, 103, and 104: (a) De = 0 and (b) De = 5.

FIG. 2.

Swimmer’s mean x position for values of the bending stiffness parameter S2 = 2, 101, 102, 103, and 104: (a) De = 0 and (b) De = 5.

Close modal

The final shape and position of the swimmer after 20 periods are shown in Fig. 3 for De = 1. There are notorious differences in the amplitude, curvature, and overall displacement of the swimmer for the various values of the bending stiffness parameter. The deviations from the target curvature are quantified in Table III. For S2 = 2, these deviations are up to 83% of the maximum of the target curvature and the majority of the points on the swimmer (94%) experience a departure from the target curvature greater or equal than 10%. S2 = 10 and S2 = 100 yield a swimmer whose curvature deviates as much as 19% and 7% from the target curvature, respectively, and there is a sizable curvature variation for a substantial number of points on the swimmer. For the largest value, S2 = 104, such deviations are negligible (about 0.4%) and the swimmer effectively follows the given gait.

FIG. 3.

Shape and location of the swimmer after 20 periods for De = 1 and different values of the bending stiffness parameter S2.

FIG. 3.

Shape and location of the swimmer after 20 periods for De = 1 and different values of the bending stiffness parameter S2.

Close modal
TABLE III.

Relative deviation from the target curvature and percentage of marker points on the swimmer where such deviation is more than 10% (last row), for various values of the bending stiffness parameter S2 and De = 1.

S2210102103104
κ̄κκ̄ 0.83 0.19 0.07 0.01 0.004 
countκ̄κκ̄0.1Nb 0.94 0.82 0.30 0.04 0.00 
S2210102103104
κ̄κκ̄ 0.83 0.19 0.07 0.01 0.004 
countκ̄κκ̄0.1Nb 0.94 0.82 0.30 0.04 0.00 

We now look at the dynamics of a stiff finite swimmer with a given gait (S2 = 104) for different Deborah numbers. Based on our resolution study, we fix h = 1/256 and set all the other parameters as in Table II. The inextensibility parameter S1 is set to 107 to guarantee an accurate enforcement of this constraint. The size of the total elastic force on the swimmer depends on the ratio S1/S2, which we find, via numerical experiments, needs to be of O(103) to adequately penalize extensibility variations and maintain the swimmer’s structure. The use of such large stiffness parameters with practical time steps is made possible by the semi-implicit IB method.4 

Figures 4 and 5 show a clear monotonic ordering of the mean displacement and locomotion speed with respect to De (after 20 periods) for this fixed gait swimmer. The steady state mean speed of the swimmer decreases with an increase in De, i.e., viscoelasticity hinders locomotion, as observed in the experiments of Shen and Arratia20 for C. elegans nematode and of Qin, Gopinath, Yang, Gollub, and Arratia19 for a more flexible algal flagellum. This occurs for both types of beating patterns, (34) and (35). Note that there is also a convergent trend toward an apparent minimal speed, consistent with the experiments of Shen and Arratia20 and with the asymptotic results for spatially periodic swimmers.9,10,13,14 Before reaching steady state, however, there is a period of time (t < 3, see Fig. 5) during which the monotonic displacement ordering with respect to De is reversed. This transient behavior is an effect of larger relaxation times associated with larger Deborah numbers.

FIG. 4.

Mean x position of a stiff swimmer with a given gait (bending stiffness parameter S2 = 104) versus time for various De. After an initial transient, there is a monotonic decrease in the mean x position with De as plot (b) shows. Viscoelasticity hinders locomotion of a swimmer with a prescribed gait.

FIG. 4.

Mean x position of a stiff swimmer with a given gait (bending stiffness parameter S2 = 104) versus time for various De. After an initial transient, there is a monotonic decrease in the mean x position with De as plot (b) shows. Viscoelasticity hinders locomotion of a swimmer with a prescribed gait.

Close modal
FIG. 5.

Newtonian-normalized mean speed versus time (in periods) for a (stiff) swimmer with a prescribed gait (bending stiffness parameter S2 = 104) for various values of the Deborah number De.

FIG. 5.

Newtonian-normalized mean speed versus time (in periods) for a (stiff) swimmer with a prescribed gait (bending stiffness parameter S2 = 104) for various values of the Deborah number De.

Close modal

We look at the locomotion as the bending stiffness parameter S2 parameter is varied. As discussed above, S2 controls how much the swimmer is allowed to deviate from the target curvature to respond to the flow. All the other parameters are set as in Table II.

Figure 6 shows the average final speed, after 20 periods, normalized by the corresponding Newtonian swimmer speed for four values of S2. The final mean speed is computed by averaging the forward propulsion of the swimmer’s center of mass throughout the final period. We observe a clearly distinct behavior between a soft swimmer (S2 = 2, 10), which is allowed to respond to the flow, and a stiff swimmer (S2 = 103, 104), which follows strictly a prescribed gait. For the stiff swimmer, as observed before (Fig. 5), viscoelasticity monotonically hinders locomotion. However, for a soft swimmer the propulsion speed increases with De. But despite this increasing trend, Fig. 6 shows no indication that the speed of swimmer in the viscoelastic medium would ever surpass that of the swimmer in the Newtonian flow. This is in apparent disagreement with the results in Refs. 25 and 26, where a speed-up with respect to the swimmer in the Newtonian flow is reported for De ≈ 1. As we document below, this seeming discrepancy might be attributed to the details of the model, specially the definition of the elastic force at the end points of the swimmer and whether or not diffusion is introduced in the viscoelastic stress.

FIG. 6.

Newtonian-normalized mean speed after 20 periods versus the Deborah number De for four different values of the bending stiffness parameter S2. Markedly distinct behaviors between a soft swimmer (S2 = 2, 10), which is allowed to respond to the flow, and a stiff swimmer (S2 = 103, 104), which follows a prescribed gait, are observed.

FIG. 6.

Newtonian-normalized mean speed after 20 periods versus the Deborah number De for four different values of the bending stiffness parameter S2. Markedly distinct behaviors between a soft swimmer (S2 = 2, 10), which is allowed to respond to the flow, and a stiff swimmer (S2 = 103, 104), which follows a prescribed gait, are observed.

Close modal

Based on measurements of the mean square displacement of flow particles, Shen and Arratia20 point out that the fluid transport induced by the nematode’s swimming is non-diffusive and that large viscoelastic stresses hinder this transport. Without significant diffusion, these very large stress values occur in highly localized regions (e.g., hyperbolic points), where the polymer molecules can elongate substantially and lead to a noticeable resistance to the locomotion and the flow transport. The numerical studies of Teran et al.25 and of Thomases and Guy26 show noticeable stress diffusion, either explicitly added or numerically induced.

Within the Oldroyd B model, viscoelastic stress diffusion could be linked to the Brownian motion of the center of mass of a polymer molecule (modeled as an elastic, Hookean dumbbell) as pointed out by El-Kareh and Leal.5 Taking such Brownian motion into account the Oldroyd-B constitutive equation becomes

Sp=1DeISp+δ2Sp.
(41)

However, as El-Kareh and Leal show, the dimensionless diffusion coefficient δ (a constant, not to be confused with the delta function approximation in the IB method) is minuscule. For a microscopic swimmer such as mammalian sperm, the simple dimensional analysis in Ref. 5 yields that an upper bound for δ is of order 10−5.

We explore the effects of viscoelastic stress diffusion, (a) introduced numerically via the approximation of the advection term u ⋅ ∇Sp in the viscoelastic stress evolution equation (10) and (b) added explicitly [the δ ∇2Sp in (41)]. The parameters remain the same as in Table II with a resolution of h = 1/128 and Δt = 0.3hs. We employ first order upwinding to compute u ⋅ ∇Sp and to assess the effect of this approximation we use a more accurate, third order ENO scheme.

We look first at the effect of viscoelastic stress diffusion on a stiff swimmer with a given gait (bending stiffness parameter S2 = 104). Figure 7 displays the Newtonian-normalized swimmer’s mean speed, averaged over the final (20th) period, as a function of the Deborah number, for various values of the stress diffusion coefficient δ and the corresponding results using third order ENO scheme and δ = 0. Notably, setting δ = 1.28h = 0.01 results in a drastic enhancement of the mean swimming velocity. For this value of δ, the viscoelastic swimmers are all faster than the corresponding swimmer in the Newtonian flow. A similar behavior of speed-up and an increasing locomotion speed with De are observed in the flexible swimmer experiments of Espinosa-García, Lauga, and Zenit. Viscoelastic swimmer speed-up has also been reported in Refs. 25 and 26, as indicated above.

FIG. 7.

The effect of viscoelastic stress diffusion on a swimmer with a given gait (S2 = 104). Newtonian-normalized velocities after 20 periods for various stress diffusion coefficients δ. All the results were computed using a first order, upwind approximation of u ⋅ ∇Sp except those labeled “ENO” which were obtain using a third order approximation for this advection term and setting δ = 0.

FIG. 7.

The effect of viscoelastic stress diffusion on a swimmer with a given gait (S2 = 104). Newtonian-normalized velocities after 20 periods for various stress diffusion coefficients δ. All the results were computed using a first order, upwind approximation of u ⋅ ∇Sp except those labeled “ENO” which were obtain using a third order approximation for this advection term and setting δ = 0.

Close modal

Our study underlines the decisive role that stress diffusion plays to achieve such locomotion enhancement and the strong dependence on the details of the motion, the model and its approximation. In particular, by comparing with the δ = 0 and the third order ENO results, we observe in Fig. 7 the striking effect of the numerical dissipation introduced by the first order upwind scheme. Moreover, the use of the high order ENO scheme makes it possible to accurately resolve the highly localized and large viscoelastic stress and with that a limiting trend toward a minimal locomotion speed as De → ∞ becomes apparent. It is however, extremely challenging to capture this perceptible limit numerically as larger De simulations become computationally formidable. Highly localized flow structures and viscoelastic instabilities are ubiquitous for large De.3 

Figure 8 presents a plot of the trace of Sp, which is proportional to the mean-square polymer extension (the trace of QQ), and the shape of the swimmer at t = 20 and for De = 1. The strong stress diffusion either explicitly added [Fig. 8(a)] or numerically introduced [Fig. 8(b)] is quite noticeable. Moreover, for δ = 2.5 × 10−5 [Fig. 8(d)], which corresponds to an upper limit of a physically based stress diffusion coefficient for this system as discussed above, the result is qualitatively the same as that with δ = 0 [Fig. 8(c)].

FIG. 8.

The effect of viscoelastic stress diffusion a swimmer with a given gait (S2 = 104). The trace of Sp, a measure of mean square polymer extension, after 20 periods for De = 1. (a) Viscoelastic stress diffusion coefficient δ = 0.01, (b) δ = 0 and first order approximation of u ⋅ ∇Sp, (c) δ = 0 and third order (ENO) approximation of u ⋅ ∇Sp, (d) δ = 2.5 × 10−5 and third order (ENO) approximation of u ⋅ ∇Sp.

FIG. 8.

The effect of viscoelastic stress diffusion a swimmer with a given gait (S2 = 104). The trace of Sp, a measure of mean square polymer extension, after 20 periods for De = 1. (a) Viscoelastic stress diffusion coefficient δ = 0.01, (b) δ = 0 and first order approximation of u ⋅ ∇Sp, (c) δ = 0 and third order (ENO) approximation of u ⋅ ∇Sp, (d) δ = 2.5 × 10−5 and third order (ENO) approximation of u ⋅ ∇Sp.

Close modal

The effect of stress diffusion is even more extraordinary for a soft swimmer. Figure 9 presents the Newtonian-normalized mean speed after 20 periods of a swimmer with S2 = 10 against De for several values of the dimensionless stress diffusion parameter δ. There is a monotonic speed-up over the Newtonian swimmer for δ ≥ 10−3 starting at De ≈ 1. However, as discussed earlier, this critical value of δ is still significantly larger than the diffusion coefficient stemming from the Brownian motion of the polymer’s center of mass for typical viscoelastic media in which the mammalian sperm locomotion takes place.

FIG. 9.

Newtonian-normalized mean speed after 20 periods versus the Deborah number De for a soft swimmer (S2 = 10) and for different values of dimensionless stress diffusion coefficient δ. There is a monotonic speed-up over the Newtonian swimmer for δ ≥ 10−3 starting at De ≈ 1.

FIG. 9.

Newtonian-normalized mean speed after 20 periods versus the Deborah number De for a soft swimmer (S2 = 10) and for different values of dimensionless stress diffusion coefficient δ. There is a monotonic speed-up over the Newtonian swimmer for δ ≥ 10−3 starting at De ≈ 1.

Close modal

We presented a numerical investigation of an IB model of a finite, effectively inextensible swimmer in a two-dimensional Stokesian Oldroyd B flow. The study shows that for this model:

  • When the swimmer follows a given gait (large bending stiffness constant S2), viscoelasticity hinders locomotion monotonically with De and the mean locomotion speed exhibits a trend toward a minimal asymptotic value.

  • When the swimmer is allowed to significantly deviate from the target gait, the monotonically decreasing speed order with De is broken and, for some parameters, can be inverted. However, no speed-up relative to the swimmer in the Newtonian flow was observed in the absence of significant stress diffusion.

  • Stress diffusion consistently enhances locomotion and, if large enough, it can lead to a propulsion speed larger than that in the Newtonian medium for a flexible swimmer with a sufficiently small bending stiffness constant. For some parameters, the mean locomotion speed was observed to increase with De.

These results do not preclude the possibility that speed-up could be achieved under different choices of parameters, force, and stress approximations. We performed simulations with a different approximation of the delta function in the IB method (the function ϕ2C in Ref. 27). While there was increase of about 6%-7% of the propulsion speed, this occurred for both the viscoelastic and the Newtonian flow. Consequently, there was no qualitative change in the behavior of the normalized mean propulsion speed. Nevertheless, other force spreadings could potentially yield qualitatively different results.

The model has two evident limitations. It is two-dimensional and the fluid is Oldroyd B. Swimming of micro-organisms typically occurs in viscoelastic fluids whose rheology is not appropriately captured by the simple Oldroyd B fluid. In particular, the unrestricted extensibility of the polymer chains in the Oldroyd B model could lead to unbounded stresses and non-physical behavior (see e.g., Ref. 3). A less perceptible limitation of the model is the treatment of the finiteness of the swimmer in the computation of the elastic forces via a discrete version of a variational formulation, (15) and (16). In other words, the swimmer’s boundary conditions and the discretization of the elastic energy. To illustrate the issue, let us consider the discrete tensile energy in (23),

EhsT[X]=S1hsj=0NbFT(DhsXk,DhsYk),
(42)

where FT(Xs,Ys)=12(Xs2+Ys21)2. To obtain the corresponding force, the derivative of EhsT with respect to Xj, Yj, j = 0, 1, …, Nb needs to be evaluated. For j = 2…Nb − 2, there is convergence to the tensile force in the continuum level

limhs01hsEhsTXj=S1ddsD1FT(Xs(sj),Ys(sj)),
(43)

where D1 stands for the derivative with respect to the first argument. A similar limit holds for EhsTYj, for j = 2…Nb − 2. But this is not the case at the end points. For example, for j = 0,

limhs01S1hsEhsTX0=limhs01hsD1FT(Xs(s0),Xs(s0))+34ddsD1FT(Xs(s0),Xs(s0)).
(44)

This would yield an unbounded force if D1FT(Xs(s0), Xs(s0)) does not go to zero fast enough with hs. The value of D1FT(Xs(s0), Xs(s0)) is proportional to the target condition. Thus, the elastic forces are largely concentrated at the end points, particularly when the inextensibility condition is not strictly enforced. Similar issues arise for the bending force, which is intimately linked to the discretization and approximation of the curvature throughout the swimmer. This might explain why stress diffusion has such a dramatic effect in the flow and the locomotion for this model.

The authors would like to thank Mr. Felipe N. Franco for sharing his numerical results, which motivated us to take a closer look at this swimmer model and to Professor Paulo Arratia for an insightful discussion. The authors gratefully acknowledge partial support for this work by the National Science Foundation under Grant Nos. DMS-1016310 (H.D.C. and D.S.) and DMS-1317684 (H.D.C.), CNPq Research Grant Nos. 309433/2011-8 and 303514/2014-0 (A.M.R.), and the Ciência Sem Fronteiras program (H.D.C.).

1. Domain size effects

We examine the effect of the periodic boundary condition on the swimmer’s dynamics by considering different domain sizes at a fixed resolution (h = 2/128). We fix De = 0 for this test. We consider a stiff swimmer with a given gait; the stiffness parameter is S2 = 104 and the inextensibility parameter is S1 = 107. All other model and numerical parameters are listed in Table II.

Figure 10 presents a plot of the mean x position of the swimmer as a function of time for the three domains considered, as well as the final position of the swimmers. There is a difference of about 2% between the [0, 2] × [0, 2] and the [0, 8] × [0, 8] results and the swimmer’s locomotion is visibly affected by that of the periodic images for the smaller domain. The difference between the [0, 4] × [0, 4] and the [0, 8] × [0, 8] results is less than 0.5%. To better afford the required high resolution for the viscoelastic simulations, we take [0, 4] × [0, 4] as our computational domain.

FIG. 10.

Effect of finite domain size on a finite, stiff, inextensible Newtonian swimmer (De = 0). The computational domains are [0, D] × [0, D] with D = 2, 4, 8. (a) Location of mean x position during final 2 periods and (b) final swimmer’s position. There is a variation of 2% ([0, 2] × [0, 2]) and 0.5% ([0, 4] × [0, 4]) in the mean position when compared with that using [0, 8] × [0, 8].

FIG. 10.

Effect of finite domain size on a finite, stiff, inextensible Newtonian swimmer (De = 0). The computational domains are [0, D] × [0, D] with D = 2, 4, 8. (a) Location of mean x position during final 2 periods and (b) final swimmer’s position. There is a variation of 2% ([0, 2] × [0, 2]) and 0.5% ([0, 4] × [0, 4]) in the mean position when compared with that using [0, 8] × [0, 8].

Close modal

2. A convergence and resolution study

We perform a convergence and resolution study with De = 0 and De = 5 (the maximum value of De considered in this work). We fixed the stiffness parameter S2 to 104 and the inextensibility parameter S1 to 107. This choice yields a stiff, effectively inextensible swimmer with a prescribed gait. Of all the cases considered here, this is the most challenging computationally. On one hand, enormous and fast oscillatory forces are generated on the swimmer with this particular values of S1 and S2 and on the other hand, for De = 5, the flow develops large and highly concentrated stresses. Moreover, there is also the potential of viscoelastic instabilities for such large De.

We consider three spatial resolutions, h = 1/64, 1/128, and 1/256, for the [0, 4] × [0, 4] domain. The time step is varied according to Δt = 0.25h and the Lagrangian mesh size is hs = 0.5h. The rest of the parameters are listed in Table II.

Figure 11 displays the mean x position of the swimmer as a function of time for the three numerical resolutions. This figure shows convergence of the swimmer’s mean displacement as the resolution is increased. The lower resolutions significantly underestimate the swimmer’s mean x position, by more than 20% for h = 1/64 (with respect to that with h = 1/256 and De = 0) and by about 7% and 4% for h = 1/128 with De = 0 and De = 5, respectively.

FIG. 11.

Swimmer’s mean x position for spatial resolutions h = 1/64, 1/128, 1/256: (a) De = 0 and (b) De = 5. The mean propulsion speed is significantly underestimated by the lower resolutions.

FIG. 11.

Swimmer’s mean x position for spatial resolutions h = 1/64, 1/128, 1/256: (a) De = 0 and (b) De = 5. The mean propulsion speed is significantly underestimated by the lower resolutions.

Close modal
1.
R. B.
Bird
,
C. F.
Curtiss
,
R. C.
Armstrong
, and
O.
Hassager
,
Dynamics of Polymeric Liquids. Volume 2: Kinetic Theory
(
John Wiley and Sons
,
New York
,
1987
).
2.
H. D.
Ceniceros
and
J. E.
Fisher
, “
A fast, robust, and non-stiff immersed boundary method
,”
J. Comput. Phys.
230
(
12
),
5133
5153
(
2011
).
3.
H. D.
Ceniceros
and
J. E.
Fisher
, “
Peristaltic pumping of a viscoelastic fluid at high occlusion ratios and large Weissenberg numbers
,”
J. Non-Newtonian Fluid Mech.
171
,
31
41
(
2012
).
4.
H. D.
Ceniceros
,
J. E.
Fisher
, and
A. M.
Roma
, “
Efficient solutions to robust, semi-implicit discretizations of the immersed boundary method
,”
J. Comput. Phys.
228
,
7137
7158
(
2009
).
5.
A. W.
El-Kareh
and
L. G.
Leal
, “
Existence of solutions for all Deborah numbers for a non-Newtonian model modified to include diffusion
,”
J. Non-Newtonian Fluid Mech.
33
,
257
287
(
1989
).
6.
G. J.
Elfring
and
E.
Lauga
, “
Theory of locomotion through complex fluids
,” in
Complex Fluids in Biological Systems: Experiment, Theory, and Computation
, edited by
S. E.
Spagnolie
,
Biological and Medical Physics Biomedical Engineering
(
Springer
,
2015
), pp.
283
317
.
7.
J.
Espinosa-Garcia
,
E.
Lauga
, and
R.
Zenit
, “
Fluid elasticity increases the locomotion of flexible swimmers
,”
Phys. Fluids
25
,
031701
(
2013
).
8.
L. J.
Fauci
and
C. S.
Peskin
, “
A computational model of aquatic animal locomotion
,”
J. Comput. Phys.
77
,
85
108
(
1988
).
9.
H. C.
Fu
,
T. R.
Powers
, and
C. W.
Wolgemuth
, “
Theory of swimming filaments in viscoelastic media
,”
Phys. Rev. Lett.
99
,
258101
(
2007
).
10.
H. C.
Fu
,
C. W.
Wolgemuth
, and
T. R.
Powers
, “
Swimming speeds of filaments in nonlinearly viscoelastic fluids
,”
Phys. Fluids
21
,
033102
(
2009
).
11.
E. A.
Gaffney
,
H.
Gadêlha
,
D. J.
Smith
,
J. R.
Blake
, and
J. C.
Kirkman-Brown
, “
Mammalian sperm motility: Observation and theory
,”
Annu. Rev. Fluid Mech.
43
,
501
528
(
2011
).
12.
R. G.
Larson
,
The Structure and Rheology of Complex Fluids
(
Oxford University Press
,
Oxford
,
1999
).
13.
E.
Lauga
, “
Propulsion in a viscoelastic fluid
,”
Phys. Fluids
19
,
083104
(
2007
).
14.
E.
Lauga
, “
Life at large Deborah number
,”
Europhys. Lett.
86
,
64001
(
2009
).
15.
S.
Lim
and
C. S.
Peskin
, “
Simulations of the whirling instability by the immersed boundary method
,”
SIAM J. Sci. Comput.
25
(
6
),
2066
2083
(
2004
).
16.
B.
Liu
,
T. R.
Powers
, and
K. S.
Breuer
, “
Force-free swimming of a model helical flagellum in viscoelastic fluids
,”
Proc. Natl. Acad. Sci. U. S. A.
108
,
19516
19520
(
2011
).
17.
S. D.
Olson
,
S. S.
Suarez
, and
L. J.
Fauci
, “
Coupling biochemistry and hydrodynamics captures hyperactivated sperm motility in a simple flagellar model
,”
J. Theor. Biol.
283
(
1
),
203
216
(
2011
).
18.
C. S.
Peskin
, “
Numerical analysis of blood flow in the heart
,”
J. Comput. Phys.
25
,
220
252
(
1977
).
19.
B.
Qin
,
A.
Gopinath
,
J.
Yang
,
J. P.
Gollub
, and
P. E.
Arratia
, “
Flagellar kinematics and swimming of algal cells in viscoelastic fluids
,”
Sci. Rep.
5
,
9190
(
2015
).
20.
X. N.
Shen
and
P. E.
Arratia
, “
Undulatory swimming in viscoelastic fluids
,”
Phys. Rev. Lett.
106
,
208101
(
2011
).
21.
A. A.
Shirgaonkar
,
M. A.
MacIver
, and
N. A.
Patankar
, “
A new mathematical formulation and fast algorithm for fully resolved simulation of self-propulsion
,”
J. Comput. Phys.
228
,
2366
2390
(
2009
).
22.
C. W.
Shu
and
S.
Osher
, “
Efficient implementation of essentially non-oscillatory shock capturing schemes
,”
J. Comput. Phys.
88
,
439
(
1988
).
23.
J.
Sznitman
and
P. E.
Arratia
, “
Locomotion through complex fluids: An experimental view
,” in
Complex Fluids in Biological Systems: Experiment, Theory, and Computation
, edited by
S. E.
Spagnolie
,
Biological and Medical Physics Biomedical Engineering
(
Springer
,
2015
), pp.
283
317
.
24.
G. I.
Taylor
, “
Analysis of the swimming of microscopic organisms
,”
Proc. R. Soc. London, Ser. A
209
,
447
461
(
1951
).
25.
J.
Teran
,
L.
Fauci
, and
M.
Shelley
, “
Viscoelastic fluid response can increase the speed and efficiency of a free swimmer
,”
Phys. Rev. Lett.
104
(
3
),
38101
(
2010
).
26.
B.
Thomases
and
B.
Guy
, “
Mechanisms of elastic enhancement and hindrance for the finite length undulatory swimmers in viscoelastic fluids
,”
Phys. Rev. Lett.
113
,
098102
(
2014
).
27.
A.-K.
Tornberg
and
B.
Engquist
, “
Numerical approximations of singular source terms in differential equations
,”
J. Comput. Phys.
200
,
462
488
(
2004
).