# Swimming of an assembly of rigid spheres at low Reynolds number

## Abstract

A matrix formulation is derived for the calculation of the swimming speed and the power required for swimming of an assembly of rigid spheres immersed in a viscous fluid of infinite extent. The spheres may have arbitrary radii and may interact with elastic forces. The analysis is based on the Stokes mobility matrix of the set of spheres, defined in low Reynolds number hydrodynamics. For small amplitude swimming optimization of the swimming speed at given power leads to an eigenvalue problem. The method allows straightforward calculation of the swimming performance of structures modeled as assemblies of interacting rigid spheres.

###### pacs:

47.15.G-, 47.63.mf, 47.63.Gd, 45.40.Ln## I Introduction

There is considerable recent interest in the construction of microswimmers capable of locomotion in a viscous fluid 1 ()-3 (). In the case of microorganisms the periodic changes of shape resulting in overall swimming motion are generated internally 4 (). In some robotic microswimmers the distortions of the body are generated from the outside 5 ()-7 (). The theory of swimming at low Reynolds number is nontrivial, and it is desirable to have simple models for which questions of principle can be elucidated and explicit calculations performed.

One such model is the three-sphere model first introduced by Najafi and Golestanian 8 (). In the model three spheres are aligned and motions along the line are considered. Hydrodynamic interactions were approximated by the Oseen tensor. Detailed analytic results were derived by Golestanian and Ajdari 9 (). The model and generalizations thereof were studied in computer simulation by Earl et al. 10 (). Vladimirov 10A () studied the generalization to collinear spheres by means of a two-timing method. The effect of elastic direct interactions was studied in a model with the motion generated by actuating forces 11 ().

In the following we study swimming of more general sphere assemblies. It is shown that the swimming velocity and the rate of dissipation are conveniently calculated in a matrix formalism. As a starting point we use the exact mobility matrix, based on the equations of Stokes flow, for an arbitrary number of spheres immersed in infinite viscous fluid. Since the mobility matrix depends only on relative coordinates, the kinematics of motion can be prescribed in -dimensional relative coordinate space. It is shown that a cyclic path in relative space leads to a net translational shift of the assembly in ordinary space. Matrix expressions are derived for the mean swimming velocity, defined as the shift divided by the period of time, and for the mean rate of dissipation.

For small amplitude motion the mean swimming velocity and the mean rate of dissipation are given by bilinear expressions in terms of the displacements in relative space. These can be written as expectation values of two hermitian matrices for a complex -dimensional displacement vector characterizing the swimming stroke. The ratio of the two expressions provides a natural measure of the efficiency of the stroke. Optimization of the efficiency leads to an eigenvalue problem.

Numerical calculations have shown that the equations of Stokesian dynamics for a set of spheres lead to numerical instability of the solution, unless the spheres are subject to attractive direct interactions keeping them together 11 (). For simplicity one can assume harmonic elastic interactions between the spheres, bound harmonically to sites on a given equilibrium structure. The matrix formulation is extended to include the case where the forces on the spheres are a sum of harmonic interactions and internal or external actuating forces. The method allows straightforward analysis of the swimming performance of structures modeled as assemblies of freely rotating rigid spheres.

As an illustration we consider the three-sphere model of Najafi and Golestanian 8 (). The method works straightforwardly. In particular it can be used to optimize the swimming stroke for the model. The calculations can be performed for improved approximations to the hydrodynamic interactions, beyond the Oseen monopole approximation.

## Ii Displacement and swimming velocity

We consider a set of rigid spheres of radii immersed in a viscous incompressible fluid of shear viscosity . The fluid is of infinite extent in all directions. At low Reynolds number and on a slow time scale the flow velocity and the pressure satisfy the Stokes equations 12 ()

(1) |

The flow velocity is assumed to satisfy the no-slip boundary condition on the surface of the spheres. The fluid is set in motion by time-dependent motions of the spheres. At each time the velocity field tends to zero at infinity, and the pressure tends to the constant ambient pressure . We shall study periodic relative motions which lead to swimming motion of the collection of spheres.

We assume that the motion is caused by time-dependent periodic forces which satisfy the condition that their sum vanishes at any time. The forces are transmitted by the spheres to the fluid. The spheres can rotate freely, so that they exert no torques on the fluid. Hence the rotational velocities can be ignored. The translational velocities are linearly related to the forces,

(2) |

with translational mobility tensors . The tensors have many-body character and depend in principle on the positions of all particles 13 ()-15 (). By translational invariance only relative distance vectors occur in the functional dependence. We abbreviate Eq. (2.2) as

(3) |

with a symmetric mobility matrix . Conversely

(4) |

with friction matrix . The friction matrix is the inverse of the mobility matrix, , and is also symmetric. Henceforth we omit the superscripts for brevity.

The positions of the centers change as a function of time. The equations of motion of Stokesian dynamics read

(5) |

The explicit time-dependence on the right originates in the time-dependence of the forces . We assume that the forces are periodic in time with period , so that . As mentioned, we impose the condition that at no time there is a net force acting on the set of spheres, so that

(6) |

We shall show that the solution of the set of nonlinear differential equations (2.5) takes the form

(7) |

where the positions describe the rigid body motion of a configuration with mean swimming velocity , and the displacements are periodic in time, . Hence is the net shift of a configuration in period .

The condition Eq. (2.6) allows us to eliminate one of the forces, for example . Correspondingly the velocity can also be eliminated. Since the mobility matrix depends only on relative coordinates, we can then reduce the spatial dimension of the algebraic problem. Eliminating we define reduced velocities

(8) |

with reduced mobility tensors

(9) |

Correspondingly we denote and , and define the corresponding reduced mobility matrix and friction matrix by

(10) |

The matrices and are not symmetric. The velocity of the -th sphere is given by

(11) |

with dimensionless transfer tensors given by

(12) |

It is convenient to define relative coordinates as

(13) |

From the corresponding differentials we define displacements for as the solution of the equations

(14) |

Via the tensors the displacements depend on the relative coordinates , and therefore are uniquely defined at each accessible point of the -dimensional relative configuration space as linear combinations of the differentials . In abbreviated notation Eq. (2.14) and its inverse read

(15) |

with and with -dimensional matrices and which depend on the relative coordinates. We define in addition

(16) |

In the kinematic formulation of swimming the relative coordinates and their time derivatives are prescribed as periodic functions of time with frequency . One cycle corresponds to a closed loop in space. It follows from Eqs. (2.14) and (2.16) that after a period the integral of the displacements over a cycle is independent of the label , so that we may define the net displacement as

(17) |

The mean swimming velocity in Eq. (2.7) is given by

(18) |

with period . The displacement clearly is a geometrical property associated with radii and the Stokes equations Eq. (2.1). It may be evaluated for any given closed path in -space and is independent of the nature of the prescribed time-dependence on the path.

For given time-dependence on a chosen closed path one may define the time-dependent displacements

(19) |

and the corresponding mean value

(20) |

with positive weights . In particular one may choose corresponding to the center of mass, or corresponding to the size distribution. The time-derivative may be defined as the instantaneous swimming velocity . These quantities clearly depend on the choice of the weights , but the net displacement and the mean swimming velocity do not.

We have defined the displacement for the mobility matrix corresponding to the set of radii and the no-slip boundary condition. For any closed path in -space the displacement may also be calculated for an approximation to the mobility matrix, for example, as given by the Oseen-interaction. The same procedure may be followed for any real -dimensional matrix depending on for which the inverse of the corresponding matrix exists.

## Iii Dissipation and efficiency

As the relative configuration of spheres varies as runs through a cycle, the corresponding flow pattern in three-dimensional space also varies. The integral of the local rate of dissipation in the fluid at any time yields the total rate of dissipation . The efficiency of swimming for the given cycle may be defined as the ratio of the net displacement and the mean rate of dissipation , where the overline indicates the average over a period.

Alternatively the instantaneous rate of dissipation may be calculated from the forces and velocities as

(21) |

For the motion we calculate velocities from Eq. (2.15) as

(22) |

The corresponding forces follow from Eq. (2.10). The velocity follows from Eq. (2.11) and is given by . The expression Eq. (3.1) may be cast in the form

(23) |

where the tilde denotes the transpose, is the identity matrix, is a -dimensional matrix constructed from a -dimensional matrix with the tensors on the diagonal and zeros elsewhere, and is a -dimensional matrix given by

(24) |

in Gibbs dyadic notation, where the symbol denotes a -dimensional vector with 1 on the positions, 0 on the positions, and cyclic. The antisymmetric part of the matrix drops out in Eq. (3.3) and we may use instead the symmetrized part

(25) |

where is the radius of the largest sphere. The matrix elements of are dimensionless. The rate of dissipation becomes

(26) |

Since the rate of dissipation is positive definite, the matrix can be used to define a Riemannian metric in the accessible part of -space with distance given by

(27) |

We consider in particular harmonically varying differentials of the form

(28) |

with infinitesimal factor and constant vectors and , which do not depend on , and which characterize the stroke. Then at each accessible point the mean swimming velocity and the mean dissipation are at least of second order in ,

(29) |

with values which depend on and . From Eqs. (2.18) and (3.8) we find

(30) |

The dimensionless efficiency for a particular stroke is defined as 16 (),17 ()

(31) |

This can be maximized with respect to the stroke at fixed values of the parameters. It was pointed out by Shapere and Wilczek 16 () that in the bilinear theory of swimming the above definition is preferable to Lighthill’s efficiency 19 (), which is essentially the ratio of the square of velocity and the dissipation. In the bilinear theory the ratio of speed and power in Eq. (3.11) is independent of the amplitude of the stroke. At any amplitude, the stroke which maximizes the ratio of speed and power for chosen values of the parameters also maximizes Lighthill’s efficiency, and minimizes the swimming drag coefficient of Avron et al. 20 ().

## Iv Optimizing efficiency

The quest for optimum efficiency of stroke at a point leads to an eigenvalue problem for the vectors and . We abbreviate

(32) |

with a matrix for . We have chosen sphere 1 for the definition, but this choice is arbitrary. Only the antisymmetric part of the matrix contributes to the swimming velocity in Eq. (3.10), so that we define the antisymmetric matrices by

(33) |

The swimming velocity in Eq. (3.10) can then be expressed as

(34) |

The associated eigenvalue problem is expressed conveniently in complex notation. We introduce the complex dimensionless vector

(35) |

and the dimensionless hermitian matrices

(36) |

We ask for the stroke with maximum swimming velocity in a class of strokes with equal rate of dissipation for fixed values of the geometric parameters, fixed frequency , and fixed viscosity . This leads to the eigenvalue problem

(37) |

The eigenvalues are real. The real and imaginary parts of the eigenvectors yield the corresponding and . The eigenvalues occur in pairs with eigenvectors which are complex conjugate, corresponding to swimming in opposite directions. It is convenient to use the notation

(38) |

with the scalar product

(39) |

The maximum efficiency for motion in direction is given by the maximum eigenvalue as

(40) |

The set depends on the choice of Cartesian coordinate system. Further optimization may be possible by a rotation of axes. In particular cases a natural choice of axes will suggest itself.

We recall that the matrices and depend on the point in -space under consideration. At a chosen point the maximum eigenvalue for optimized choice of axes yields the maximum swimming velocity for given power, and the corresponding eigenvector yields the nature of the corresponding stroke. The amplitude factor in Eq. (3.9) implies that the amplitude must be small. It is natural to consider larger amplitudes by use of Eq. (3.8) and calculate the mean swimming velocity and rate of dissipation by use of Eqs. (2.18) and (3.6) for the corresponding cyclic path centered about .

## V Harmonic interactions and actuating forces

In the formulation of the mobility matrix in Eq. (2.2) the nature of the forces need not be specified. In an earlier calculation 11 () we have considered microswimmers with internal harmonic interactions, driven by actuating forces. The forces are of the form

(41) |

with displacements

(42) |

where are equilibrium positions at time . The constant harmonic interaction tensors are symmetric and satisfy . The actuating forces are assumed to satisfy

(43) |

They can be internal or generated from the outside. As seen from Eq. (2.13) the differences depend linearly on the relative coordinates , so that for an equilibrium configuration with relative distances the first forces may be expressed as

(44) |

with a real matrix . It is convenient to rewrite the second order mean swimming velocity as a bilinear expression in terms of the displacements and forces . To simplify notation we agree that the subscripts comprise both particle label and Cartesian index. The subscripts are reserved for particle labels. It follows from Eqs. (2.15) and (2.17) that the mean second order swimming velocity can be expressed as

(45) |

where summation over repeated indices is implied. The first order equations of motion in relative coordinate space are given by

(46) |

Substituting into Eq. (5.5) we obtain for the swimming velocity

(47) |

where the matrix is independent of time and has matrix elements

(48) |

The matrix is independent of the particle label . Similarly, from Eq. (3.6) the second order mean rate of dissipation may be expressed as

(49) |

with matrix , given by Eq. (3.5), calculated at .

Substituting Eq. (5.4) into Eq. (5.6) and assuming harmonically varying actuating forces one obtains in complex notation

(50) |

with the solution

(51) |

This relation may be used to express the swimming velocity in Eq. (5.7) and the mean rate of dissipation in Eq. (5.9) in terms of the actuating forces . The resulting dynamic expression for the efficiency is maximized by the method of Sec. IV and calculation of the corresponding actuating force amplitudes and , defined in analogy to Eq. (3.8), from Eq. (5.11).

## Vi Three-sphere swimmer

The simplest application of the theory is to a three-sphere swimmer with three spheres aligned on the axis, as studied by Golestanian and Ajdari 9 (). The arms mentioned by these authors do not appear in the theory, and the derived equations apply to three spheres in the absence of direct interactions. The authors define the instantaneous swimming velocity as the mean . As seen above, this definition is somewhat arbitrary. One could also consider the center of mass or the center of sizes. As we have shown, for periodic motion the mean velocity of any sphere, averaged over a period, provides the mean swimming velocity.

In the example the and coordinates can be ignored. There are only two relative coordinates and , and the various matrices are two-dimensional. The elements of the mobility matrix are approximated by use of the Oseen interaction as 12 ()

(52) |

In the bilinear theory we consider a point in -space with coordinates . It is not necessary to assume the ratios to be small, as done by Golestanian and Ajdari. The more complete expressions will provide a better approximation to the situation with exact hydrodynamic interactions. Similarly one could use the more complicated Rotne-Prager expressions 18 () instead of Eq. (6.1).

As an example we consider the case of equal-sized spheres with and equal distances between centers . The matrix in Eq. (4.5) takes the form

(53) |

with element

(54) |

We note that the corresponding swimming velocity given by Eq. (4.7) with when expanded to lowest order in the ratio agrees with Eq. (12) of ref. 9. The matrix takes the form

(55) |

These expressions show a singularity at , which is less than the minimum distance . The matrix , defined in Eq. (2.15), is singular at . From Eq. (4.6) one finds the eigenvalues

(56) |

as well as the corresponding eigenvectors with

(57) |

normalized to . The maximum efficiency, corresponding to by Eq. (4.9), tends to zero monotonically as the ratio tends to infinity.

It is of interest to compare the above analytical results for small amplitude motion with values obtained by numerical solution of the equations of motion Eq. (2.5) with hydrodynamic interactions given by Eq. (6.1) and prescribed oscillating forces. As we have noted earlier 11 (), the numerical solution is unstable unless one introduces harmonic elastic forces as in Eq. (5.4), keeping the particles together. With suitable harmonic forces the motion quickly tends to a limit cycle. The latter corresponds to the periodic motion studied by Golestanian and Ajdari 9 (). In our numerical work we use harmonic interactions given by the matrix

(58) |

with elastic constant . This corresponds to nearest neighbor interactions of equal strength between the three spheres. The stiffness of the swimmer is characterized by the dimensionless number defined by

(59) |

We consider actuating forces oscillating at frequency with given by Eq. (5.11) with , and . This corresponds to maximum efficiency in the bilinear theory. We choose initial conditions

(60) |

In the bilinear theory, corresponding to small , the orbit in relative space is given by with and

(61) |

independent of the elastic constant .

In Fig. 1 we show the orbit for and ratio for one period. In the figure the numerical solution cannot be distinguished from the elliptical orbit given by Eq. (6.10). By use of the Stokes parameters 21 () calculated from one finds that the tilt angle in the plane is , independent of the ratio , and obtains an expression for the ellipticity. In the present case the ellipticity equals 0.66.

In Fig. 2 we show the orbit for , ratio , and stiffness for . The mean swimming velocity and the mean rate of dissipation are calculated as time averages over the last period , corresponding to the limit cycle. It turns out that in the range both quantities vary approximately in proportion to . In Fig. 3 we show the reduced mean swimming velocity as a function of for . In Fig. 4 we show the reduced mean rate of dissipation and in Fig. 5 we show the efficiency as functions of . Interestingly, the efficiency increases monotonically with the amplitude factor.

At we have and . This can be compared with the numerical calculation of Alouges et al. 22 (),23 () on the basis of a Stokes solver. The authors used radius mm and period s. For viscosity of water poise our calculation yields 0.0072 mm and . The latter value is about the same as the one given in Table 1 of ref. 27, and the displacement agrees well with the value mm of Alouges et al.. A more precise comparison would require a calculation with the same hydrodynamic interactions in both procedures. In our calculation the efficiency for given amplitude could be maximized numerically by variation of the vector in the neighborhood of and variation of the stiffness . Also one can vary the equilibrium situation used as a starting point of the bilinear theory. The present method allows fast and straightforward design of an efficient large amplitude mechanical swimmer corresponding to an approximate form of the mobility matrix.

## Vii Discussion

The matrix-formulation presented above provides insight into the mathematical problem of the swimming of assemblies of spheres at low Reynolds number. It allows straightforward calculation of the swimming performance of assemblies of interest. A practical procedure for a particular swimmer would be to model the rest shape by a set of spheres, identifying the positions of the centers with the equilibrium sites of a structure with harmonic elastic interactions. By use of an approximation to the mobility matrix one can then evaluate the hermitian matrix , which appears in the expression for the mean swimming velocity, and the real and symmetric matrix , which appears in the expression for the mean rate of dissipation. The swimming velocity and required power can be calculated for a chosen stroke of small amplitude with harmonic time variation as expectation values of the two hermitian matrices. Large amplitudes can also be handled in principle. These require a parametric integration along a chosen closed path in the space of relative positions.

For the simple example of three collinear spheres with Oseen-type interactions, discussed in Sec. VI, the calculations can be performed in analytic form. For more complicated structures and more accurate hydrodynamic interactions the algebra rapidly becomes cumbersome, but it is straightforward to derive numerical results. Elsewhere we have applied the method to the analysis of three- and four-sphere swimmers pushing a cargo sphere 24 ().

As we have shown, the optimization of swimming at small amplitude leads to an eigenvalue problem. This allows determination of the optimum stroke yielding the largest swimming speed at given power. For an assembly of spheres with elastic interactions the required actuating forces which lead to optimal speed can be evaluated from the eigenvector with maximum eigenvalue. It is then of interest to study the swimming speed and power for the same set of force ratios as functions of an amplitude factor.

## Figure captions

### Fig. 1

Plot of the orbit in the plane for for one period.

### Fig. 2

Plot of the orbit in the plane for for ten periods. The initial values correspond to Eq. (6.9).

### Fig. 3

Plot of the reduced mean swimming velocity for as a function of the amplitude .

### Fig. 4

Plot of the reduced mean swimming power for as a function of the amplitude .

### Fig. 5

Plot of the efficiency for as a function of the amplitude .

### References

- E. Lauga and T. R. Powers, Rep. Prog. Phys. 72, 096601 (2009).
- P. Garstecki and M. Cieplak, J. Phys. C.:Condens. Matter 21, 200301 (2009).
- J. W. Swan, J. F. Brady, R. S. Moore, and ChE 174, Phys. Fluids 23, 071901 (2011).
- D. B. Dusenbery, Living at Micro Scale (Harvard University Press, Cambridge (Mass.), 2009).
- M. Roper, R. Dreyfus, J. Baudry, M. Fermiger, J. Bibette, and H. A. Stone, J. Fluid Mech. 554, 167 (2006).
- M. Roper, R. Dreyfus, J. Baudry, M. Fermiger, J. Bibette, and H. A. Stone, Proc. Roy. Soc. A 464, 877 (2008).
- E. E. Keaveny and M. R. Maxey, J. Fluid Mech. 598, 293 (2008).
- A. Najafi and R. Golestanian, Phys. Rev. E 69, 062901 (2004).
- R. Golestanian and A. Ajdari, Phys. Rev. E 77, 036308 (2008).
- D. J. Earl, C. M. Pooley, J. F. Ryder, I. Bredberg, and J. M. Yeomans, J. Chem. Phys. 126, 064703 (2007).
- V. A. Vladimirov, J. Fluid Mech. 716, R1-1 (2013).
- B. U. Felderhof, Phys. Fluids 18, 063101 (2006).
- J. E. Marsden, Lectures on Mechanics (Cambridge University Press, Cambridge, 1992).
- R. Yang and P. S. Krishnaprasad, Proc. IEEE Conf. on Decision and Control 2, 1632 (1989).
- M. Berry, Proc. Roy. Soc, London A 392, 45 (1984).
- J. Happel and H. Brenner, Low Reynolds number hydrodynamics (Noordhoff, Leyden, 1973).
- B. Cichocki, B. U. Felderhof, K. Hinsen, E. Wajnryb, and J. Blawzdziewicz, J. Chem. Phys. 100, 3780 (1994).
- B. Cichocki, M. L. Ekiel-Jeżewska, and E. Wajnryb, J. Chem. Phys. 111, 3265 (1999).
- M. L. Ekiel-Jeżewska and E. Wajnryb, in Theoretical Methods for Micro Scale Viscous Flows, edited by F. Feuillebois and A. Sellier (Transworld Research Network, Kerala, 2009).
- A. Shapere and F. Wilczek, J. Fluid Mech. 198, 587 (1989).
- B. U. Felderhof and R. B. Jones, Physica A 202, 94 (1994).
- P. J. Zuk, E. Wajnryb, K. A. Mizerski, and P. Szymczak, J. Fluid Mech. 741, R5 (2014).
- M. J. Lighthill, Commun. Pure Appl. Maths. 5, 109 (1952).
- J. E. Avron, O. Gat, and O. Kenneth, Phys. Rev. Lett. 93, 186001 (2004).
- C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (Wiley, New York, 1983).
- F. Alouges, A. DeSimone, and A. Lefebvre, J. Nonlinear Sci. 18, 277 (2008).
- F. Alouges, A. DeSimone, and A. Lefebvre, Eur. Phys. J. E 28, 279 (2009).
- B. U. Felderhof, arxiv/physics.flu-dyn1408.3530.