Effect of fluid inertia on swimming of a sphere in a viscous incompressible fluid

# Effect of fluid inertia on swimming of a sphere in a viscous incompressible fluid

## Abstract

Swimming of a sphere in a viscous incompressible fluid is studied on the basis of the Navier-Stokes equations for wave-type distortions of the spherical shape. At sizable values of the dimensionless scale number the mean swimming velocity is the result of a delicate balance between the net time-averaged flow generated directly by the surface distortions and the flow generated by the mean Reynolds force density. Depending on the stroke, this can lead to a surprising dependence of the mean swimming velocity on the kinematic viscosity of the fluid. The net flow pattern is calculated as a function of kinematic viscosity for axisymmetric strokes of the swimming sphere.

###### pacs:
47.15.G-, 47.63.mf, 47.63.Gd, 87.17.Jj

## I Introduction

The theory of swimming of nearly spherical microorganisms was developed first by Lighthill 1 () and Blake 2 (). The Reynolds number for a microorganism is small and the theory can be based on the time-independent Stokes equations 3 (),4 (). Work in this area was reviewed by Lauga and Powers 5 ().

In the swimming of larger bodies inertia of the fluid can no longer be neglected. The Reynolds stress causes a reactive flow which modifies the Stokes flow. In earlier work we discussed the effect of fluid inertia on swimming on the basis of the complete set of Navier-Stokes equations 6 (). As the simplest geometry we studied the swimming of a sphere with no-slip boundary condition due to time-harmonic surface distortions 7 (). The swimming velocity was evaluated to second order in the amplitude of distortions as an average over a period of time. In terms of the dimensionless ratio , where is the amplitude of stroke and is the ratio of the undistorted sphere, the mean swimming velocity and the mean rate of dissipation were calculated to order .

The effect of fluid inertia may be characterized by the dimensionless ratio , where is the frequency of the stroke, is the mass density of the fluid, and is the shear viscosity. For fixed and one may consider the full range of kinematic viscosity . The Stokes limit corresponds to . In this limit the Reynolds stress vanishes. For larger values of the scale number the Reynolds stress becomes important.

In earlier work we calculated the mean swimming velocity and the mean rate of dissipation of a swimming sphere in the full range of scale number 8 (). For finite values of the theory employs mode functions which become singular in the limit . This causes mathematical difficulties which obscure the relation to the Stokes limit. In our earlier work 8 () we used a modified set of mode functions for small to facilitate the study of low frequency behavior. We show here that the theory can be improved by a matrix representation based on the surface modes of the Stokes limit theory.

In the following we study the mean swimming velocity for a choice of stroke which is the same for all values of the scale number . It turns out that the swimming velocity has a surprising dependence on . For certain strokes the swimming velocity changes sign as increases.

The calculations show that for larger values of there is a delicate balance between the direct effect of surface distortion and the indirect effect of the Reynolds stress. The two contributions to the mean swimming velocity, generated directly and indirectly, nearly cancel for large .

We view the effect of swimming as force-free convection of a body in self-generated fluid flow. It is of interest to study the net time-averaged flow as a superposition of the direct and indirect contributions. Both contributions depend on the scale number . We show in examples that for fixed stroke the net flow pattern varies significantly as a function of scale number. The calculations are exact to order .

The effect of Reynolds stress on the translational velocity of a -active particle, a so-called squirmer, was studied by Wang and Ardekani 9 (), by Khair and Chisholm 10 (), and by Chisholm et al. 11 (). Spelman and Lauga 12 () studied the translational velocity of a squirmer in the inertia-dominated limit to order by the method of matched asymptotic expansion. Wang and Ardekani 13 () studied the effect of fluid inertia on swimming of small organisms via an approximate equation of motion.

## Ii Swimming of a sphere

We consider a flexible sphere of radius immersed in a viscous incompressible fluid of shear viscosity and mass density . The fluid is set in motion by time-dependent distortions of the sphere. We shall study axisymmetric periodic distortions which lead to a translational swimming motion of the sphere in the direction in a Cartesian system of coordinates. The surface displacement is written as

 \boldmathξ(θ,t)=Re[\boldmathξω(θ)e−iωt], (1)

with polar angle and complex amplitude . The corresponding first order flow velocity and pressure are given by

 \boldmathv(1)(\boldmathr,t)=Re[\boldmathvω(\boldmathr)e−iωt],p(1)(\boldmath% r,t)=Re[pω(\boldmathr)e−iωt], (2)

with amplitude functions which satisfy the linearized Navier-Stokes equations

 η[∇2\boldmathvω−α2\boldmathvω]−∇pω=0,∇⋅\boldmathvω=0, (3)

with the variable

 α=(−iωρ/η)1/2=(1−i)(ωρ/2η)1/2. (4)

The solution of Eq. (2.3) can be expressed as a linear superposition of modes 14 ()

 \boldmathvl(\boldmathr,α) = 2πeαa[(l+1)kl−1(αr)\boldmath% Al(^\boldmathr)+lkl+1(αr)\boldmathBl(^\boldmathr)], \boldmathul(\boldmathr) = −(ar)l+2\boldmathBl(^\boldmathr),pl(\boldmathr,α)=ηα2a(ar)l+1Pl(cosθ), (5)

with modified spherical Bessel functions 15 () and vector spherical harmonics defined by 16 ()

 \boldmathAl = ^\boldmathAl0=lPl(cosθ)\boldmather−P1l(cosθ)\boldmatheθ, \boldmathBl = ^\boldmathBl0=−(l+1)Pl(cosθ)% \boldmather−P1l(cosθ)\boldmatheθ, (6)

with Legendre polynomials and associated Legendre functions in the notation of Edmonds 17 (). The corresponding surface displacement function is expanded as

 \boldmathξω(^\boldmathr)=−ia∞∑l=1[κl\boldmathvl(\boldmaths,α)+μl% \boldmathul(\boldmaths)], (7)

with and complex coefficients . The first order fluid velocity at the surface is given by .

We define the complex multipole moment vector as the one-dimensional array

 \boldmathψ=(κ1,μ1,κ2,μ2,....). (8)

The absence of uniform displacement implies the constraint . We indicate arrays satisfying this constraint by a hat, as . The mean swimming velocity and the mean rate of dissipation are bilinear in the vector and can be expressed as

 ¯¯¯¯¯¯U2=12ωa(^\boldmathψ|^B|^\boldmathψ),¯¯¯¯¯¯D2=8πηω2a3(^\boldmathψ|^A|^% \boldmathψ). (9)

With truncation at maximum -value the truncated matrices and are -dimensional. The truncated matrices correspond to swimmers obeying the constraint that all multipole coefficients for vanish.

The matrices and are calculated from integrals with integrands which are bilinear in the mode functions defined in Eq. (2.5). The matrix is diagonal in , and the matrix is tridiagonal in . The matrices are frequency-dependent via the variable . We write

 αa=(1−i)s,s=a√ωρ2η (10)

with dimensionless scale number . In Appendix B of Ref. 8 we provided explicit expressions for the matrix elements of and up to .

## Iii Stokes representation

The matrices and are singular at which causes difficulties in numerical calculations and in the discussion of the relation to swimming in the Stokes limit. In our earlier work 8 () we therefore used for small a different set of mode solutions of Eq. (2.3). However, alternatively we can choose a more convenient matrix representation by expanding the surface displacement in terms of a different set of vector functions defined on the surface of the sphere . It is of particular interest to use the set of functions found as limiting values on the sphere surface of the modes defined in the Stokes limit 18 (). The mode functions in the Stokes limit are the same as in Eq. (2.5), but the functions are changed to

 \boldmathv0l(\boldmathr) = (ar)l[(l+1)Pl(cosθ)\boldmather+l−2lP1l(cosθ)\boldmatheθ] (11) = (ar)l[2l+2l(2l+1)% \boldmathAl−2l−12l+1\boldmathBl].

We denote the corresponding set of superposition coefficients as and the corresponding Stokes representation of the matrices as and . The Stokes limit is denoted as and .

In our earlier work 8 () we gave the relation between the two sets of mode coefficients and , found from a comparison of the two expansions of the surface displacement on the sphere. More generally we find by the same method the linear relation

 ^\boldmathψ=T⋅^\boldmathψI, (12)

with a transformation matrix . The matrix is block-diagonal, as given by a factor , with a 2-dimensional at order given by the relations

 κl = πl(2l+1)ezkl−1(z)κIl, μl = 12l+1[2l−1+2kl+1(z)kl−1(z)]κIl+μIl,z=(1−i)s. (13)

At we have simply . The relation between the two sets of matrices is

 ^AI=T†⋅^A⋅T,^BI=T†⋅^B⋅T, (14)

where is the Hermitian conjugate of . The mean swimming velocity and the mean rate of dissipation can be expressed alternatively as

 ¯¯¯¯¯¯U2=12ωa(^\boldmathψI|^BI|^\boldmathψI),¯¯¯¯¯¯D2=8πηω2a3(^\boldmathψI|^AI|^\boldmathψI). (15)

The explicit expression for the matrix up to order reads

 ^AI13=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝300000AI22185000185600000AI4450700050710⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠, (16)

with elements

 AI22 = 3109+18s+18s2+2s31+2s+2s2, AI44 = 221531+1062s+1062s2+708s3+244s4+8s59+18s+18s2+12s3+4s4. (17)

The matrix is simpler than . It differs only in the - and -elements from the matrix , which can be read off from Eq. (7.17) in Ref. 18.

The matrix is a sum of two terms , corresponding to contributions from surface displacements directly and from a bulk term originating in the Reynolds stress tensor. The nonvanishing elements are given by

 BIS12 = BI∗S21=3+(3−3i)s+4is25i+(5+5i)s, BIS13 = BI∗S31=−3i, BIS24 = BI∗S42=6359+18s+(18−14i)s2−(6+10i)s3−(8−8i)s43i+6is+(2+6i)s2+(2+2i)s3, BIS25 = BI∗S52=635−15i+(15−15i)s−4s21+(1+i)s, BIS34 = BI∗S43=63545+(45−45i)s−26is2+(4+4i)s33i+(3+3i)s+2s2, BIS35 = BI∗S53=−6i. (18)

The matrix is simpler than .

The nonvanishing elements are given by

 BIB12 = BI∗B21=s25−i−(1+i)s+s2−(1−i)s3−2is4F−i+(1+i)s, BIB24 = BI∗B42=s2630[3i+6is+(2+6i)s2+(2+2i)s3]× [ −216i−432is−(234+432i)s2+(198+378i)s3 − 12i(7−56i+96F2−45F++21F−)s4+(6+6i)(−20+33i+90F+−42iF−)s5, + (−51+124i−384F2+585F++255F−)s6+(1+i)(55+6i−225iF++87F−)s7 + (16+3i−117iF++3iF−)s8−(1−i)(−3+15iF++F−)s9−6F+s10], BIB25 = BI∗B52=s2420(1+(1+i)s)[18+(18+18i)s+6is2+(6−6i)s3−9s4 + 11(1+i)s5+i(1−24F+)s6+(1−i)s7−2F+s8], BIB34 = BI∗B43=s21260[3i+3+3i)s+2s2]× (19) [ −450i−(450+450i)s−222s2+(78−78i)s3+81is4 − (83+83i)s5−(1−168F−)s6+(1−i)s7+2iF−s8],

with the abbreviations

 F+=F(s+is),F−=F(s−is),F2=F(2s), (20)

where the function with complex variable is defined by

 F(z)=ezE1(z)=∫∞0e−uz+udu. (21)

As we shall see, the contributions from the matrix are important for large . In the limit the matrix vanishes, and the matrix tends to the matrix which can be read off from Eq. (7.11) in Ref. 18.

## Iv Simple swimmers

We study the effect of fluid inertia on swimming performance by calculating the mean swimming velocity as a function of scale number for fixed surface displacement . This implies a fixed set of coefficients . In order to compare different swimmers we define the dimensionless reduced swimming velocity

 Ured(s)=(^ψI|^BI(s)|^ψI)(^ψI|^A0|^ψI). (22)

The denominator provides a measure of the intensity of surface agitation. For a chosen fixed set of coefficients it is independent of scale number .

In the following we consider the simplest swimmers involving only modes of orders . For such swimmers the reduced swimming velocity can be evaluated from the explicit expressions for the matrix elements of given in Sec. III. The matrix can be read off from Eq. (7.17) in Ref. 18, or from Eq. (3.7) with .

The simplest swimmer is the potential one, superposing a dipolar and a quadrupolar flow field, corresponding to moments , , and all other moments vanishing. For this swimmer the reduced swimming velocity is , independent of . We have optimized the ratio of the two moments. With three optimized moments the swimming velocity increases to .

Next we study the so-called -swimmer, as defined in terms of the modes introduced by Lighthill 1 () and Blake 2 (). In terms of the mode coefficients defined above

 μI1=B1,κI2=−μI2,μI2=13iB2. (23)

In our scheme the first order swimming velocity vanishes, so that Lighthill’s coefficient equals . The constraint causes the first order component of the flow to be tangential to the sphere. We consider in particular the case , the same ratio as for the -active particle studied by Ishikawa et al. 19 (). In Fig. 1 we plot the reduced swimming velocity as a function of . In the Stokes limit . For large scale number .

In Fig. 1 we compare with the reduced swimming velocity of the 12-swimmer with optimal velocity in the Stokes limit, as given by the solution of the generalized eigenvalue problem 18 () with matrices and . This corresponds to moments . In the Stokes limit . For large scale number .

There are no diagonal elements in Eqs. (3.8) and (3.9). This implies that the stroke must contain at least two different modes of the chosen type. Another simple swimmer has only the coefficients and different from zero. Optimizing again in the Stokes limit we obtain the values . In Fig. 2 we plot the reduced swimming velocity for this swimmer. In the Stokes limit . For large scale number . Remarkably, the reduced swimming velocity changes sign as a function of . The swimmer is not very efficient.

In Fig. 3 we plot the reduced swimming velocity for a 123-swimmer with so-called combined stroke 6 (). The optimized coefficients are

 μI1=1,κI2=53√230413i,μI2=0,κI3=−2759,μI3=0, (24)

as calculated from the Stokes generalized eigenvalue problem with constraints . Again, in this case the reduced swimming velocity changes sign as a function of . In the Stokes limit . For large scale number .

The sign change occurs also for a 23-swimmer with moments

 μI1=0,κI2=1,μI2=0,κI3=92√7295i,μI3=0, (25)

as shown in Fig. 4. Again we optimized the moments in the Stokes limit. The corresponding reduced velocity is . For large the velocity tends to .

Finally we consider the 123-swimmer with optimized moments in the Stokes limit. The moments are

 μI1=1,κI2=−1.553i,μI2=1.824i,κI3=1.373,μI3=−1.440, (26)

corresponding to Stokes reduced velocity . For large the velocity tends to . In Fig. 5 we show the reduced velocity as a function of .

## V Analysis of mean swimming velocity

In this section we analyze the results shown above in some more detail. The potential swimmer needs no further comment. The reduced velocity of the -swimmer can be expressed as

 Ured(s)=4β90+5β212+24s+24s2+4s3−s4+2s6−is6(1+s−is)F++is6(1+s+is)F−]1+2s+2s2, (27)

where . In particular for and for

 Ured(0)=48β90+5β2,Ured(∞)=72β90+5β2. (28)

In the literature a Stokes active particle with has been called a puller, and one with a pusher 19 (). For a Stokes -swimmer, as defined here, the names are inappropriate, since one deals with the same swimmer swimming in opposite directions. The reduced velocity of the swimmer is maximal for at any . In the Stokes limit . The distinction between Stokes active particle and Stokes swimmer was discussed by one of us 20 ().

The same type of factorization as shown in Eq. (5.1) occurs for the swimmer of type considered in Eq. (4.4). More generally we put

 μI1=0,κI2=1,μI2=0,κI3=iγ,μI3=0. (29)

The reduced swimming velocity takes the form

 Ured(s)=216γ567+1180γ2G(s), (30)

where the function takes the limiting values

 G(0)=1,G(∞)=−559. (31)

The prefactor takes its maximum at the value given in Eq. (4.4).

In general the reduced swimming velocity can be written as a sum of two terms

 Ured(s)=US(s)+UB(s), (32)

corresponding to the decomposition mentioned above Eq. (3.8). The bulk term vanishes at , as is evident from Eq. (3.9). Including terms up to the normalization factor can be expressed as

 (^\boldmathψI|^A0|^%\boldmath$ψ$I) = 3|μI1|2+2710|κI2|2+365R(κI2μI∗2)+6|μI2|2 (33) + 11821|κI3|2+1007R(κI∗3μI3)+10|μI3|2.

The numerator in Eq. (4.1) at can be expressed similarly as

 (^\boldmathψI|^BI(0)|^\boldmath% ψI)=635I[7μI∗1κI2+35μI∗1μI2+6κI2κI∗3+30κI∗2μI3+30μI∗2κI3+70μI∗2μI3]. (34)

The behavior at large can be expressed as

 (^\boldmathψI|^BIS(s)|^\boldmathψI) = Ps+QS+O(1/s), (^\boldmathψI|^BIB(s)|^\boldmathψI) = −Ps+QB+O(1/s), (35)

showing that the terms linear in precisely cancel. The constant term is given by

 QS=235I[35μI∗1κI2+105μI∗1μI2+54κI∗2κI3+102κI∗2μI3+114μI∗2κI3+210μI∗2μI3], (36)

and the term is given by

 QB=−835I[14μI∗1κI2+41κI∗2κI3+33κI∗2μI3+14μI∗2κI3]+1635R(κI∗2κI3), (37)

The sum is

 QS+QB = 235I[−21μI∗1κI2+105μI∗1μI2−110κI∗2κI∗3−30κI∗2μI3+58μI∗2κI3+210μI∗2μI3] (38) + 1635R(κI∗2κI3).

It is clear that this expression can have either sign, depending on the coefficients, and that the sign is independent of that in the expression Eq. (5.8). One can check that the velocity

 Ured(∞)=QS+QB(^\boldmathψI|^A0|^\boldmathψI), (39)

in our examples takes the values given in Sec. IV. It is easy to extend the expressions in Eqs. (5.7) and (5.12) to include higher order mode coefficients.

The two expressions on the left in Eq. (5.9) originate in the mean convective flow generated directly by the surface distortions of the sphere, and the reactive flow generated by the Reynolds force density, respectively. Apparently there is a delicate balance between these two quantities. It follows from Eq. (3.9) that the bulk matrix has vanishing -elements. Although the matrix originates in the Reynolds stress of the bulk flow, there is a contribution only from the boundary layer of thickness . In the limit of large the boundary layer becomes very thin, but its effect on the sphere velocity depends on the details of its structure, as evident from Eq. (5.12).

Spelman and Lauga 12 () studied the boundary layer problem for an active particle in the inertia-dominated limit by the method of matched asymptotic expansions. For the case of radial surface displacements they find a quite complicated expression for the swimming speed involving Gaunt coefficients, which is not easy to compare with our result Eq. (5.13) for . They also find the possibility of sign changes depending on the choice of mode combinations.

## Vi Mean Reynolds force density

It is of interest to consider besides the mean swimming velocity also the net flow pattern of a swimmer. One of us studied the net flow pattern of simple swimmers in the Stokes limit 20 (). Here we extend the analysis to arbitrary values of the scale number . Both the direct contribution from the second order surface velocity and the contribution from the second order Reynolds stress depend on the scale number . For the mean second order flow we need to solve a steady state Stokes problem 6 (). The direct contribution to the net flow can be evaluated fairly straightforwardly as the solution of the Stokes problem with boundary condition given by the second order velocity on the surface of the sphere with radius .

In order to find the contribution from the Reynolds stress we must solve an inhomogeneous Stokes problem with driving term given by the mean Reynolds force density

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathf(2)R(\boldmathr)=−ρT∫T0\boldmathv(1)(\boldmathr,t)⋅∇% \boldmathv(1)(\boldmathr,t)dt, (40)

where the time-average on the right is over a period . The resulting mean flow velocity can be expressed as

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathv(2)R(\boldmathr)=1η∫r′>a\boldmathG(\boldmathr,\boldmathr′)⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathf(2)R(\boldmathr′)d\boldmathr′, (41)

where is the Green function of the Stokes equations with no-slip boundary condition on the surface of a fixed sphere with radius centered at the origin. We shall discuss the Green function in the next section.

In our axisymmetric problem the mean Reynolds force density can be expanded in vector spherical harmonics and ,

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathf(2)R(\boldmathr)=∞∑l=1[fAl(r)\boldmathAl+fBl(r)\boldmathBl]. (42)

As before it is convenient to use complex notation. We therefore express the time-average in Eq. (6.1) as

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathf(2)R(\boldmathr)=−12ρR[\boldmathv∗ω⋅∇\boldmathvω]. (43)

With the expansion corresponding to Eq. (2.7)

 \boldmathvω(\boldmathr)=−ωa∞∑l=1[κl\boldmathvl(\boldmathr,α)+μl% \boldmathul(\boldmathr)], (44)

and the understanding that , the complex function can be expressed as

 \boldmathv∗ω⋅∇\boldmathvω=∞∑l=1(^\boldmathψ|FAl(r,s)|^% \boldmathψ)\boldmathAl+∞∑l=1(^% \boldmathψ|FBl(r,s)|^\boldmathψ)% \boldmathBl, (45)

where the elements of the matrices and can be evaluated from the expansions of the bilinear expressions

 \boldmathv∗j⋅∇\boldmathvk,% \boldmathv∗j⋅∇\boldmathuk,\boldmathu∗j⋅∇\boldmathvk,\boldmathu∗j⋅∇\boldmathuk (46)

in terms of vector spherical harmonics and . These can be evaluated by use of the orthonormality relations for the and which read

 ∫π0\boldmathAk⋅\boldmathAlsinθdθ = 2kδkl,∫π0\boldmathAk⋅\boldmathBlsinθdθ=0, ∫π0\boldmathBk⋅\boldmathBlsinθdθ = (2k+2)δkl. (47)

The expression for the tensor in spherical coordinates is given by Happel and Brenner 3 ().

Finally we can rewrite

 (^\boldmathψ|FAl(r,s)|^% \boldmathψ) = (^\boldmathψI|FIAl(r,s)|^\boldmathψI), (^\boldmathψ|FBl(r,s)|^% \boldmathψ) = (^\boldmathψI|FIBl(r,s)|^\boldmathψI), (48)

by use of the transformation matrix as in Eq. (3.4). A comparison of Eqs. (6.3) and (6.6) yields expressions for the radial functions and . If the coefficient vector is truncated at low order, then only correspondingly low order matrices and need to be calculated.

## Vii Net flow pattern

The mean second order flow velocity is defined in the volume . It tends to at infinity 6 (). The net flow pattern is defined as 20 ()

 \boldmathv′(\boldmathr)=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathv% (<