Efficient swimming of an assembly of rigid spheres at low Reynolds number
aa r X i v : . [ phy s i c s . f l u - dyn ] M a y Efficient swimming of an assembly of rigid spheres at lowReynolds number
B. U. Felderhof ∗ Institut f¨ur Theorie der Statistischen PhysikRWTH Aachen UniversityTemplergraben 5552056 AachenGermany (Dated: February 18, 2018)
Abstract
The swimming of an assembly of rigid spheres immersed in a viscous fluid of infinite extentis studied in low Reynolds number hydrodynamics. The instantaneous swimming velocity andrate of dissipation are expressed in terms of the time-dependent displacements of sphere centersabout their collective motion. For small amplitude swimming with periodically oscillating dis-placements, optimization of the mean swimming speed at given mean power leads to an eigenvalueproblem involving a velocity matrix and a power matrix. The corresponding optimal stroke per-mits generalization to large amplitude motion in a model of spheres with harmonic interactions andcorresponding actuating forces. The method allows straightforward calculation of the swimmingperformance of structures modeled as assemblies of interacting rigid spheres. A model of threecollinear spheres with motion along the common axis is studied as an example.
PACS numbers: 47.15.G-, 47.63.mf, 47.63.Gd, 45.50.Jf ∗ Electronic address: [email protected] . INTRODUCTION In earlier work [1] we presented a method to analyze the performance of a microswim-mer modeled as an assembly of N rigid spheres immersed in a viscous incompressible fluidof infinite extent, with a no-slip boundary condition on the surface of each sphere. Themotion of the whole system is determined by the Stokes equations of low Reynolds numberhydrodynamics. The swimming motion of such a system was discussed earlier by Alougeset al. [2],[3]. The particular case of collinear spheres was studied by Vladimirov [4] using atwo-timing method.For small displacements of the spheres from fixed positions in the collective rest framethe time-averaged swimming velocity and rate of dissipation can be evaluated in terms of a(3 N − × (3 N −
3) velocity matrix and a (3 N − × (3 N −
3) power matrix, which can beconstructed from the mobility matrix for each relative rest configuration [1]. Optimizationof the velocity at fixed power leads to a generalized eigenvalue problem involving the twomatrices. Optimal efficiency corresponds to the maximum eigenvalue.In a model with harmonically interacting spheres the optimal stroke of small amplitudemotion can be used to calculate a set of corresponding actuating forces. Large amplitudemotion can be studied by solving the equations of Stokesian dynamics for the same actuatingforces multiplied by a factor. The mean swimming velocity and the mean power of the largeamplitude motion can then be determined numerically from the limit cycle of the solution.In the following we present an alternative method based on a purely kinematic point ofview. Expressions are derived for the instantaneous swimming velocity and power in termsof the sphere displacements from the center and their instantaneous time derivative. Thisallows calculation of the mean swimming velocity and mean power for given periodic strokeof any amplitude. The present method also provides an alternative derivation of the velocitymatrix and power matrix of small amplitude motion.For large amplitude swimming the present method is more straightforward than the earlierone [1], since it does not require numerical solution of the equations of Stokesian dynamics.A large amplitude stroke may be determined by amplifying the optimal stroke found fromthe eigenvalue problem of the small amplitude theory for a given equilibrium structure. Theinstantaneous swimming velocity and power are then determined from explicit expressionsin terms of the given displacements. Subsequently the mean swimming velocity and meanpower can be found by integration over a period.Both methods are tested on a model of three collinear spheres with motion along thecommon axis, as formulated by Najafi and Golestanian [5] and studied in detail by Golesta-nian and Ajdari [6]. The two methods of calculation lead to similar numerical results for awide range of amplitude.
2. DISPLACEMENT AND SWIMMING VELOCITY
We consider a set of N rigid spheres of radii a , ..., a N immersed in a viscous incompressiblefluid of shear viscosity η . The fluid is of infinite extent in all directions. At low Reynoldsnumber and on a slow time scale the flow velocity v and the pressure p satisfy the Stokesequations [7] η ∇ v − ∇ p = 0 , ∇ · v = 0 . (2.1)2he flow velocity v is assumed to satisfy the no-slip boundary condition on the surface ofthe spheres. The fluid is set in motion by time-dependent motions of the spheres. At eachtime t the velocity field v ( r , t ) tends to zero at infinity, and the pressure p ( r , t ) tends tothe constant ambient pressure p . We shall study periodic relative motions which lead toswimming motion of the collection of spheres.We assume that the motion is caused by time-dependent periodic forces F ( t ) , ..., F N ( t )which satisfy the condition that their sum vanishes at any time. The forces are transmittedby the spheres to the fluid. The spheres can rotate freely, so that they exert no torques onthe fluid. Hence the rotational velocities Ω ( t ) , ..., Ω N ( t ) can be ignored. The translationalvelocities U , ..., U N are linearly related to the forces, U j = N X k =1 µ ttjk · F k , j = 1 , ..., N, (2.2)with translational mobility tensors µ ttjk . The tensors have many-body character and dependin principle on the positions of all particles [8]-[10]. By translational invariance only relativedistance vectors { R i − R j } occur in the functional dependence. We abbreviate eq. (2.2) as U = µ · F , (2.3)with a symmetric 3 N × N mobility matrix µ . Conversely F = ζ · U , (2.4)with friction matrix ζ . The friction matrix is the inverse of the mobility matrix, ζ = µ − ,and is also symmetric.The positions of the centers change as a function of time. The equations of motion ofStokesian dynamics read d R j dt = U j ( R , ..., R N , t ) , j = 1 , ..., N. (2.5)The explicit time-dependence on the right originates in the time-dependence of the forces F ( t ). In the swimming motion the forces are periodic in time with period T , so that F ( t + T ) = F ( t ). As mentioned, we impose the condition that at no time is there a net force acting onthe set of spheres, so that N X j =1 F j ( t ) = 0 . (2.6)We look for a solution of eq. (2.5) corresponding to swimming motion, of the form R j ( t ) = S j + Z t U ( t ′ ) dt ′ + δ j ( t ) , j = 1 , ..., N, (2.7)where the first two terms describe the collective motion of the configuration S =( S , ..., S N ) with swimming velocity U ( t ) caused by the displacements { δ j ( t ) } . We re-quire that the latter are periodic with period T , and exclude uniform displacements, so thatthe 3 N -dimensional vector d ( t ) = { δ ( t ) , ..., δ N ( t ) } satisfies d ( t ) · u α = 0 , ( α = x, y, z ) , (2.8)3here the symbol u x denotes a 3 N -dimensional vector with 1 on the x positions, 0 on the y, z positions, and cyclic. Periodicity implies U ( t + T ) = U ( t ) , d ( t + T ) = d ( t ) . (2.9)The mean swimming velocity is defined as U sw = 1 T Z T U ( t ) dt. (2.10)We require that d ( t ) is purely oscillating, so that Z T d ( t ) dt = . (2.11)We show in the following that the instantaneous swimming velocity U ( t ) can be calcu-lated from the displacement vector d ( t ) and its time derivative ˙ d ( t ). Later we compare thepresent kinematic description to a dynamical model, in which the forces are decomposedinto actuating forces and elastic restoring forces.
3. SWIMMING VELOCITY AND DISSIPATION
By substitution of eq. (2.7) into eqs. (2.4) and (2.5) one finds F = ζ · ( U β u β + ˙ d ) , (3.1)where summation over repeated greek indices is implied. The condition (2.6) can be ex-pressed as u α · F = 0, so that Z αβ U β = − u α · ζ · ˙ d (3.2)with friction tensor Z αβ = u α · ζ · u β . (3.3)Hence we obtain the swimming velocity U α = − M αβ u β · ζ · ˙ d , (3.4)where M αβ is the inverse of the friction tensor. The 3 N × N friction matrix ζ depends onlyon the instantaneous relative positions. Therefore the friction tensor Z and the mobilitytensor M depend on the displacement vector d , but not on the central coordinates R α = u α · R /N .By series expansion of the mobility tensor M and the friction matrix ζ in powers of thedisplacement vector d we obtain a corresponding expansion of the swimming velocity U = U (1) + U (2) + U (3) + ..., (3.5)with first order term U (1) α = − M αβ u β · ζ · ˙ d , (3.6)with mobility tensor M αβ and friction matrix ζ calculated for the configuration S . Byperiodicity of d ( t ) the time average of the first order swimming velocity vanishes, U (1) = .4e introduce the friction vectors f α = u α · ζ = ζ · u α , (3.7)where we have used the symmetry of the friction matrix ζ . The vectors are related to thefriction tensor by u α · f β = u β · f α = Z αβ . (3.8)From the Taylor series expansion of eq. (3.4) we find that the second order instantaneousswimming velocity can be expressed as U (2) α = − d · V α (cid:12)(cid:12) · ˙ d , (3.9)with matrix V α given by V α = ∇ (cid:2) M αβ f β (cid:3) , (3.10)where ∇ is the gradient vector in 3 N -dimensional configuration space. The notation (cid:12)(cid:12) ineq. (3.9) indicates that the matrix-function is to be evaluated at R = S .The expression on the right of eq. (3.10) may be written as a sum of two terms, V α = ( ∇ M αβ ) f β + M αβ D β , (3.11)with derivative friction matrix D β = ∇ f β . (3.12)We introduce the gradient vectors g βγ = D β · u γ = ∇ Z βγ , (3.13)and use the identity Z αγ M γβ = δ αβ (3.14)to show that ∇ M αβ = − M αγ g γδ M δβ . (3.15)Then eq. (3.11) may be expressed alternatively as V α = M αβ ˘ D β , (3.16)with reduced derivative friction matrix˘ D β = D β − g βγ M γδ f δ . (3.17)This matrix has the property ˘ D β · u α = 0 . (3.18)From the fact that ζ depends only on relative coordinates it follows that u α · ∇ ζ = , andhence u α · D β = , u α · g βγ = 0 . (3.19)As a consequence u α · V β = , V α · u β = . (3.20)5he time-dependent rate of dissipation can be expressed in the same matrix formalism.The rate of dissipation is given by D = F · U = F · ˙ d , (3.21)since F · u α = 0 on account of the condition eq. (2.6). Substituting eq. (3.1) we find D = ˙ d · ζ · ˙ d + U α ˙ d · f α . (3.22)It follows from eq. (3.4) that the rate of dissipation is at least of second order in d and ˙ d .To second order, by use of eq. (3.6), D (2) = ˙ d · P · ˙ d (3.23)with matrix P = ζ − M αβ f α f β . (3.24)The matrix is symmetric and has the properties u α · P = , P · u α = . (3.25)The properties eq. (3.20) and (3.25) allow us to reduce the dimension of the matrix descrip-tion by three by the introduction of center and relative coordinates.
4. VELOCITY MATRIX VECTOR AND POWER MATRIX
The center of the assembly is given by R = 1 N N X j =1 R j = 1 N e α u α · R (4.1)with Cartesian unit vectors e α . We define relative coordinates { r j } as r = R − R , r = R − R , ..., r N − = R N − R N − , j = 1 , ..., N − , (4.2)and the corresponding (3 N − r = ( r , ..., r N − ). The 3 N -vector ( R , r ) is related tothe vector R by a transformation matrix T according to( R , r ) = T · R (4.3)with explicit form given by eqs. (4.1) and (4.2).The matrices V α and P are transformed to V αT = T · V α · T − , P T = T · P · T − . (4.4)The first three rows of T consist of u α /N and the first three columns of T − consist of u α .It follows from the properties eq. (3.20) and (3.25) that the first three rows and columns ofthe transformed matrices V αT and P T vanish identically. Hence in this representation we candrop the center coordinates and truncate the matrices by erasing the first three rows and6olumns. We denote the truncated (3 N − × (3 N − V αT and ˆ P T and definedisplacements ξ in relative space by ( , ξ ) = T · d . (4.5)With this notation the second order swimming velocity and rate of dissipation are given by U (2) α = ξ · C T · ˆ V αT · ˙ ξ , D (2) = ˙ ξ · C T · ˆ P T · ˙ ξ , (4.6)with the matrix C T = [ g T − · T − ] ˆ . (4.7)This (3 N − × (3 N −
3) dimensional matrix consists of numerical coefficients and is obtainedfrom the corresponding 3 N × N matrix by truncation, as indicated by the final hat symbol.We consider in particular harmonically varying displacements of the form d ( t ) = d s sin ωt + d c cos ωt, (4.8)with a corresponding expression for ξ ( t ). The time-averaged second order swimming velocityand rate of dissipation are then given by U (2) α = 12 ω (cid:2) ξ s · C T · ˆ V αT (cid:12)(cid:12) · ξ c − ξ c · C T · ˆ V αT (cid:12)(cid:12) · ξ s (cid:3) , D (2) = 12 ω (cid:2) ξ s · C T · ˆ P T · ξ s + ξ c · C T · ˆ P T · ξ c (cid:3) . (4.9)We introduce the complex dimensionless vector ξ c = 1 b ( ξ c + i ξ s ) , (4.10)where b is a typical length scale. With the definitions B α = 12 ib (cid:0) C T · ˆ V αT (cid:12)(cid:12) − ^ C T · ˆ V αT (cid:12)(cid:12) (cid:1) , A = 1 bη C T · ˆ P T , (4.11)and the scalar product ( ξ c | η c ) = N − X j =1 ξ c ∗ j · η cj (4.12)the mean swimming velocity and mean rate of dissipation can then be expressed as U (2) α = 12 ωb ( ξ c | B α | ξ c ) , D (2) = 12 ηω b ( ξ c | A | ξ c ) . (4.13)We have normalized such that the matrix elements of B α and A are dimensionless. We call B α the velocity matrix and A the power matrix.We ask for the stroke with maximum swimming velocity in a class of strokes with equalrate of dissipation for fixed values of the geometric parameters, fixed frequency ω , and fixedviscosity η . This leads to the generalized eigenvalue problem B α ξ c = λ α A ξ c . (4.14)7he eigenvalues { λ α } are real. The maximum efficiency for motion in direction α is givenby the maximum eigenvalue as E αT max = λ αmax . (4.15)The set { E xT max , E yT max , E zT max } depends on the choice of Cartesian coordinate system. Fur-ther optimization may be possible by a rotation of axes. In particular cases a natural choiceof axes will suggest itself.In the formulation of the mobility matrix in Eq. (2.2) the nature of the forces { F j } need not be specified. In an earlier calculation [11] we considered microswimmers withinternal harmonic interactions, driven by actuating forces. In matrix form the forces maybe expressed as F = E + H · ( R − S ) , (4.16)with a real symmetric matrix H with the property H · u α = 0. The actuating forces { E j ( t ) } are assumed to satisfy N X j =1 E j ( t ) = 0 . (4.17)They can be generated internally or externally.
5. THREE-SPHERE SWIMMER
The simplest application of the theory is to a three-sphere swimmer with three spheresaligned on the x axis, as studied by Golestanian and Ajdari [6]. The spheres move along the x axis, and the y and z coordinates can be ignored. There are only two relative coordinates r = x − x and r = x − x , and the relevant parts of the matrices B x and A are two-dimensional. The elements of the 3 × µ ttjk = 16 πη (cid:20) a j δ jk + 32 | x j − x k | (1 − δ jk ) (cid:21) . (5.1)In the bilinear theory we consider a point r in r -space with coordinates ( d , d ), corre-sponding to the configuration S of the rest system. As an example we consider the case ofequal-sized spheres with a = a = a = a and equal distances between centers d = d = d .For this case the explicit expressions for the matrices B x and A are identical to those derivedearlier by a different method [1]. Explicit expressions for the eigenvectors ξ ± and eigenvalues λ ± of the two-dimensional eigenvalue problem B x · ξ = λ A ξ , as functions of the ratio d/a ,were derived in ref. 1.In the bilinear theory, corresponding to small ε , the orbit ( r ( t ) , r ( t )) = ( x ( t ) − x ( t ) , x ( t ) − x ( t )) in relative space is given by r ( t ) = r + ξ ( t ) with r = ( d, d ) and ξ ( t ) = εa Re ξ + exp( − iωt ) , (5.2)with amplitude factor ε and eigenvector ξ + = (1 , ξ + ) corresponding to the largest eigenvalue.In fig. 1 of ref. 1 we have shown the elliptical orbit in relative space for d = 5 a and ε = 0 . d ( t ) = T − · (cid:18) ξ ( t ) (cid:19) , T =
13 13 13 − − . (5.3)8n fig. 1 we show the reduced mean swimming velocity U sw / ( ε ωa ) as a function of ε for d = 5 a , as calculated from eq. (3.4). In fig. 2 we show the reduced mean rate of dissipation D / ( ε ηω a ), as calculated from eq. (3.22). In fig. 3 we show the efficiency E T = ηωa U sw / D as a function of ε . The efficiency increases monotonically with the amplitude factor.It is of interest to compare the above results with values obtained by the numericalsolution of the Stokesian equations of motion eq. (2.5) with hydrodynamic interactionsgiven by eq. (5.1) and prescribed oscillating actuating forces. We use harmonic interactionsgiven by the 3 × H = k − − − (5.4)with elastic constant k . This corresponds to nearest neighbor interactions of equal strength k between the three spheres. The stiffness of the swimmer is characterized by the dimensionlessnumber σ defined by σ = kπηaω . (5.5)In general, the first order forces F (1)0 ( t ) corresponding to the displacement vector d ( t ) andthe corresponding first order swimming velocity U (1)0 ( t ), calculated from eq. (3.6), followfrom eq. (3.1) as F (1)0 = ζ · ( U (1)0 β u β + ˙ d ) . (5.6)In the present case only the x components are relevant. The corresponding actuating forces E ( t ) are found from eq. (4.16) as E ( t ) = F (1)0 ( t ) − H · d ( t ) . (5.7)These have the property u α · E ( t ) = 0, so that the sum of actuating forces vanishes. Wechoose initial conditions for the x coordinates x (0) = 0 , x (0) = d + εa, x (0) = 2 d + εa + εa Re ξ + . (5.8)In fig. 4 we show the numerical solution of the equations of Stokesian dynamics eq. (2.5)with forces given by F ( t ) = E ( t ) + H · ( R ( t ) − S ) (5.9)for d = 5 a , stiffness σ = 1, and amplitude factor ε = 2 for the first ten periods. We comparethe orbit with the ellipse given by eq. (5.2). The mean swimming velocity and mean power,calculated as time-averages over the last period for values of the amplitude factor in therange 0 < ε <
2, are shown in figs. 1 and 2. The corresponding efficiency is shown infig. 3. The dashed curves in figures 1 − ε calculated by the kinematic methodis always larger than that calculated by the dynamic method from the limit cycle withactuating forces. However, we must compare the mean swimming velocity for two differentstrokes of the same mean power. In fig. 5 we plot the power as a function of ε in the range1 . < ε < D = 52 ηω a of themean power occurs at ε k = 1 .
949 in the kinematic method, and at ε d = 1 .
970 in the dynamic9ethod. For these values the mean swimming velocity is found to be U sw = 0 . ωa forthe elliptical orbit of the kinematic method, and U sw = 0 . ωa for the limit cycle of thedynamical method. Thus in the present case the elliptical orbit is the most efficient of thetwo. This does not exclude that for the same power an orbit with yet higher speed can befound.At ε = 1 .
38 and for d = 5 a we have U sw ≈ . ωa from eq. (3.4) and D ≈ . ηω a from eq. (3.22) for the orbit given by eq. (5.3). This can be compared with the numericalcalculation of Alouges et al. [2],[3] on the basis of a Stokes solver. The authors used radius a = 0 .
05 mm, and period T = 1 s. For viscosity of water η = 0 .
01 poise our calculationyields ∆ = U sw T ≈ D T ≈ . × − J . The latter value is somewhatless than the one given in table 1 of ref. 3, and the displacement agrees well with the value0 .
01 mm of Alouges et al..Finally we consider the efficiency calculated from eqs. (3.4) and (3.22) for displace-ment in relative space of the form eq. (5.2), but with the eigenvector ξ + replaced by ξ = (1 , A exp( iδ )) with absolute value A and phase δ . The values of A and δ can be relatedto the Stokes parameters of the elliptical orbit [12]. In fig. 6 we show the efficiency foramplitude factor ε = 2 and ratio d/a = 5 as a function of A and δ . The maximum is notvery pronounced.
6. DISCUSSION
The swimming performance of an assembly of spheres as a function of the amplitude ofa chosen stroke can be studied in a purely kinematic formulation. From eq. (3.4) we findthe instantaneous swimming velocity, and from eq. (3.22) we find the instantaneous rate ofdissipation or power. The mean swimming velocity and the mean power follow by averagingover a period. The ratio of these two quantities yields the efficiency of the stroke.Alternatively one may use a dynamic approach [1],[11] in which the swimmer is modeled asa set of spheres bound harmonically to equilibrium positions and with harmonic interactions.The spheres are subject to actuating forces which sum to zero. The corresponding swimmingmotion may be found as the limit cycle of the solution of the equations of Stokesian dynamics.The mean swimming velocity and the mean power may be found numerically from the limitcycle.We have shown in sect. 5 that for a collinear three-sphere swimmer the two methodslead to similar results over a wide range of amplitude, provided that for small amplitude theactuating forces correspond to the chosen kinematic stroke. We have chosen the latter to bethe optimal one at small amplitude, as determined from the velocity matrix and the powermatrix of the bilinear theory.The kinematic method is the more straightforward one, since it does not require numericalsolution of the equations of Stokesian dynamics. The dynamic approach has the advantagethat it provides a physical model of the swimmer. It will be of interest to explore thedifference in efficiency for given stroke or given actuating forces as a function of amplitudefactor for more sophisticated model swimmers, with actuating forces chosen to agree withthe optimal stroke at small amplitude. 10
1] B. U. Felderhof, Eur. Phys. J. E ,110 (2014).[2] F. Alouges, A. DeSimone, and A. Lefebvre, J. Nonlinear Sci. , 277 (2008).[3] F. Alouges, A. DeSimone, and A. Lefebvre, Eur. Phys. J. E , 279 (2009).[4] V. A. Vladimirov, J. Fluid Mech. , R1-1 (2013).[5] A. Najafi and R. Golestanian, Phys. Rev. E , 062901 (2004).[6] R. Golestanian and A. Ajdari, Phys. Rev. E , 036308 (2008).[7] J. Happel and H. Brenner, Low Reynolds number hydrodynamics (Noordhoff, Leyden, 1973).[8] B. Cichocki, B. U. Felderhof, K. Hinsen, E. Wajnryb, and J. Blawzdziewicz, J. Chem. Phys. , 3780 (1994).[9] B. Cichocki, M. L. Ekiel-Je˙zewska, and E. Wajnryb, J. Chem. Phys. , 3265 (1999).[10] M. L. Ekiel-Je˙zewska and E. Wajnryb, in
Theoretical Methods for Micro Scale Viscous Flows ,edited by F. Feuillebois and A. Sellier (Transworld Research Network, Kerala, 2009).[11] B. U. Felderhof, Phys. Fluids , 063101 (2006).[12] C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (Wiley, New York, 1983). igure captionsFig. 1 Plot of the reduced mean swimming velocity U sw / ( ε ωa ) for d = 5 a as a function ofthe amplitude ε as calculated by the kinematic method (solid curve), and by the dynamicmethod with stiffness parameter σ = 1 (dashed curve). Fig. 2
Plot of the reduced mean swimming power D / ( ε ηω a ) for d = 5 a as a function ofthe amplitude ε as calculated by the kinematic method (solid curve), and by the dynamicmethod with stiffness parameter σ = 1 (dashed curve). Fig. 3
Plot of the efficiency E T = ηωa U sw / D for d = 5 a as a function of the amplitude ε ascalculated by the kinematic method (solid curve), and by the dynamic method with stiffnessparameter σ = 1 (dashed curve). Fig. 4
Plot of the orbit in the r r plane calculated from the equations of Stokesian dynamicsfor d = 5 a, ε = 2 , σ = 1 for ten periods. The initial values correspond to Eq. (5.8) andthe forces follow from eq. (5.9). We also plot the elliptical orbit for d = 5 a, ε = 2 (dashedcurve). Fig. 5
Plot of the mean swimming power D / ( ηω a ) for d = 5 a as a function of the amplitude ε in the range 1 . < ε < σ = 1 (dashed curve). Fig. 6