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.

## I. INTRODUCTION

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, experiments^{7,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 Arratia^{20} found that for *C. elegans nematodes* in shear thinning polymeric fluids, viscoelasticity hinders locomotion while the experiments of Liu, Powers, and Breuer^{16} 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 Zenit^{7} 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 Shelley^{25} 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 Guy^{26} 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 Shelley^{25} 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 method^{2,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.

## II. THE MODEL

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 Taylor^{24} 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 as^{1,12}

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 formula

^{1}

In (3), *n* is the number density of polymer chains, *k _{B}* 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 equation

^{1}

where *λ _{H}* is the relaxation time constant of the Hookean dumbbells,

*η*=

_{p}*nk*the polymeric viscosity, and

_{B}Tλ_{H} 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

($kBT/H$ is proportional to the polymer chain’s radius of gyration), we get a simplified constitutive equation^{12}

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

where *p* and **F** have been redefined appropriately. Throughout this study, we take $\eta p\eta 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

where

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 framework^{18} 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

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

The swimmer’s gait is generated by penalizing deviations from a prescribed, time periodic target $\kappa \u0304(s,t)$. We follow the approach proposed by Fauci and Peskin^{8} to generate the waving motion and to enforce the inextensibility condition through a density force given by

where *E*[**X**] is the elastic energy

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

Some authors have proposed linking the model parameter *S*_{2} 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, *S*_{2} 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 *S*_{2} 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 *S*_{2} with the flexural rigidity we obtain that *S*_{2} = 10 corresponds to a flexural rigidity of 10^{−20} N m^{2}, which is within the range of that elastic parameter for mammalian sperm.

Henceforth we will refer to *S*_{2} as the bending stiffness parameter of the model.

## III. THE NUMERICAL METHOD

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 **S**_{p}. The swimmer Γ is of extent *L* and is discretized with a uniform grid $Ghs$ in *s*, consisting of the points *s _{j}* =

*jh*,

_{s}*j*= 0, …,

*N*with

_{b}*h*=

_{s}*L*/

*N*.

_{b}Let **D**_{h} and *L _{h}* 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

*t*, we discretize the Stokes equations and the swimmer’s evolution equation as

_{n} where Δ*t* is the time step size. $Sn$ and $Sn\u2217$ are the IB Method spreading and interpolation operators,

where $\delta h(x)=\phi x\phi y$, with **x** = (*x*, *y*) and

is an approximation to the 2D delta distribution.^{18}

To compute the force density **f**_{hs}(**X**), we employ the following discretization of the energy:

where $D\u0304hs$ stands for the centered finite difference for *j* = 1…*N _{b}* − 1 and for the forward and backward finite differences for

*j*= 0 and

*j*=

*N*, respectively. $\u2211\u2033$ denotes the trapezoidal rule sum with a factor of 1/2 for the end points. The mean curvature

_{b}*κ*is computed from the swimmer configuration as

_{j} where *D*_{hs} and $Dhs2$ are the second order, centered, finite difference approximations to the first and second derivative, respectively.

The elastic force density is given by

This discretization satisfies $\u2211kfhs(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, **u**^{n+1} and **X**^{n+1}. To solve this system we use the methodology developed in Refs. 2 and 4, which we briefly outline next.

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 **X**^{n+1} only,

where $Mn$ the *flow-structure* operator is given by $Mn=\Delta tSn\u2217LhSn$ and

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 **X**^{n+1} is computed, **u**^{n+1} is obtained from (26) and the excess stress is updated. This is done using **u**^{n+1} in an explicit, second order, total variation diminishing (TVD) Runge-Kutta method^{22}

Here **G**_{h} is the discrete version of (12) obtained by using standard second order spatial finite differences for all the terms except for the advection **u** ⋅ ∇**S**_{p}, 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.

## IV. NUMERICAL VALIDATION

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*) = *a*sin(*kx* + *ωt*), the ratio of an Oldroyd B swimmer (mean) speed *U _{V}* to the Newtonian fluid swimmer speed (Taylor

^{24})

satisfies

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}*/

*η*= 1/2 and

_{s}*De*= 1.

Domain Ω | [0, 1] × [0, 1] |

Mesh size h = 1/N | N = 64, 128, 256, 512 |

Lagrangian mesh size h _{s} | 0.5h |

Time step Δt | h _{s} |

Inextensibility parameter S_{1} | 10^{6} |

Bending stiffness parameter S_{2} | 10^{4} |

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 h _{s} | 0.5h |

Time step Δt | h _{s} |

Inextensibility parameter S_{1} | 10^{6} |

Bending stiffness parameter S_{2} | 10^{4} |

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*(*ak*^{2}), 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).

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

## V. NUMERICAL RESULTS

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

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

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 (*x*_{0}(*s*), *y*_{0}(*s*)) we solve the system,

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

by the composite trapezoidal rule quadrature for *j* = 0, …, *N _{b}*. Here

*N*= [

_{b}*L*/

*h*] and

_{s}*s*=

_{j}*jh*.

_{s}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.

Lagrangian mesh size h _{s} | 0.5h |

Time step Δt | 0.5h _{s} |

Inextensibility parameter S_{1} | As specified in each simulation |

Bending stiffness parameter S_{2} | As specified in each simulation |

Swimmer’s length L | 0.6 |

Wave number k | $2\pi L$ |

Frequency ω | 2π |

Final Time | 20 |

Amplitude A | 0.1 |

Viscosity ratio polymer/solvent η/_{p}η _{s} | $12$ |

Lagrangian mesh size h _{s} | 0.5h |

Time step Δt | 0.5h _{s} |

Inextensibility parameter S_{1} | As specified in each simulation |

Bending stiffness parameter S_{2} | As specified in each simulation |

Swimmer’s length L | 0.6 |

Wave number k | $2\pi 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.

### A. Propulsion as the bending stiffness parameter *S*_{2} is varied

One of the key parameters in the model is the bending stiffness constant *S*_{2}, which penalizes deviations from the target curvature. Effectively, for sufficiently large *S*_{2}, the swimmer will follow closely the given gait generated by the time-periodic target curvature. As *S*_{2} is reduced, larger deviations from the target curvature are allowed in response to the flow. Thus, large and small values of *S*_{2} 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 *S*_{2} 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 *S*_{2} for (a) a Newtonian flow (*De* = 0) and (b) a viscoelastic flow with *De* = 5. The propulsion is drastically affected by the value of *S*_{2}. The swimmer with *S*_{2} = 10^{4} achieves a displacement in excess of 50% larger than that of the swimmer with *S*_{2} = 10. For *S*_{2} = 2, the swimmer propulsion is minimal, with a final displacement after 20 periods of less than 17% its length.

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 *S*_{2} = 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%. *S*_{2} = 10 and *S*_{2} = 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, *S*_{2} = 10^{4}, such deviations are negligible (about 0.4%) and the swimmer effectively follows the given gait.

S_{2}
. | 2 . | 10 . | 10^{2}
. | 10^{3}
. | 10^{4}
. |
---|---|---|---|---|---|

$\kappa \u0304\u2212\kappa \u221e\kappa \u0304\u221e$ | 0.83 | 0.19 | 0.07 | 0.01 | 0.004 |

$count\kappa \u0304\u2212\kappa \kappa \u0304\u22650.1Nb$ | 0.94 | 0.82 | 0.30 | 0.04 | 0.00 |

S_{2}
. | 2 . | 10 . | 10^{2}
. | 10^{3}
. | 10^{4}
. |
---|---|---|---|---|---|

$\kappa \u0304\u2212\kappa \u221e\kappa \u0304\u221e$ | 0.83 | 0.19 | 0.07 | 0.01 | 0.004 |

$count\kappa \u0304\u2212\kappa \kappa \u0304\u22650.1Nb$ | 0.94 | 0.82 | 0.30 | 0.04 | 0.00 |

### B. Swimmer with a prescribed gait

We now look at the dynamics of a stiff finite swimmer with a given gait (*S*_{2} = 10^{4}) 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 *S*_{1} is set to 10^{7} to guarantee an accurate enforcement of this constraint. The size of the total elastic force on the swimmer depends on the ratio *S*_{1}/*S*_{2}, which we find, via numerical experiments, needs to be of *O*(10^{3}) 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 Arratia^{20} for *C. elegans nematode* and of Qin, Gopinath, Yang, Gollub, and Arratia^{19} 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 Arratia^{20} 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.

### C. Flexible swimmer

We look at the locomotion as the bending stiffness parameter *S*_{2} parameter is varied. As discussed above, *S*_{2} 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 *S*_{2}. 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 (*S*_{2} = 2, 10), which is allowed to respond to the flow, and a stiff swimmer (*S*_{2} = 10^{3}, 10^{4}), 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.

### D. The effect of viscoelastic stress diffusion

Based on measurements of the mean square displacement of flow particles, Shen and Arratia^{20} 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 Guy^{26} 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

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** ⋅ ∇**S _{p}** in the viscoelastic stress evolution equation (10) and (b) added explicitly [the

*δ*∇

^{2}

**S**in (41)]. The parameters remain the same as in Table II with a resolution of

_{p}*h*= 1/128 and Δ

*t*= 0.3

*h*. We employ first order upwinding to compute

_{s}**u**⋅ ∇

**S**and to assess the effect of this approximation we use a more accurate, third order ENO scheme.

_{p}We look first at the effect of viscoelastic stress diffusion on a stiff swimmer with a given gait (bending stiffness parameter *S*_{2} = 10^{4}). 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.28*h* = 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.

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 **S**_{p}, 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)].

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 *S*_{2} = 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.

## VI. CONCLUDING REMARKS

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

*S*_{2}), 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 $\varphi 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),

where $FT(Xs,Ys)=12(Xs2+Ys2\u22121)2$. To obtain the corresponding force, the derivative of $EhsT$ with respect to *X _{j}*,

*Y*,

_{j}*j*= 0, 1, …,

*N*needs to be evaluated. For

_{b}*j*= 2…

*N*− 2, there is convergence to the tensile force in the continuum level

_{b} where *D*_{1} stands for the derivative with respect to the first argument. A similar limit holds for $\u2202EhsT\u2202Yj$, for *j* = 2…*N _{b}* − 2. But this is not the case at the end points. For example, for

*j*= 0,

This would yield an unbounded force if *D*_{1}*F ^{T}*(

*X*(

_{s}*s*

_{0}),

*X*(

_{s}*s*

_{0})) does not go to zero fast enough with

*h*. The value of

_{s}*D*

_{1}

*F*(

^{T}*X*(

_{s}*s*

_{0}),

*X*(

_{s}*s*

_{0})) 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.

## Acknowledgments

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.).

### APPENDIX: DOMAIN SIZE AND NUMERICAL RESOLUTION

#### 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 *S*_{2} = 10^{4} and the inextensibility parameter is *S*_{1} = 10^{7}. 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.

#### 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 *S*_{2} to 10^{4} and the inextensibility parameter *S*_{1} to 10^{7}. 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 *S*_{1} and *S*_{2} 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.25*h* and the Lagrangian mesh size is *h _{s}* = 0.5

*h*. 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.