Swimming at small Reynolds number of a planar assembly of spheres in an incompressible viscous fluid with inertia

# Swimming at small Reynolds number of a planar assembly of spheres in an incompressible viscous fluid with inertia

B. U. Felderhof Institut für Theorie der Statistischen Physik
RWTH Aachen University
Templergraben 55
52056 Aachen
Germany
July 5, 2019
###### Abstract

Translational and rotational swimming at small Reynolds number of a planar assembly of identical spheres immersed in an incompressible viscous fluid is studied on the basis of a set of equations of motion for the individual spheres. The motion of the spheres is caused by actuating forces and forces derived from a direct interaction potential, as well as hydrodynamic forces exerted by the fluid as frictional and added mass hydrodynamic interactions. The translational and rotational swimming velocities of the assembly are deduced from momentum and angular momentum balance equations. The mean power required during a period is calculated from an instantaneous power equation. Expressions are derived for the mean swimming velocities and the mean power, valid to second order in the amplitude of displacements from the relative equilibrium positions. Hence these quantities can be evaluated for prescribed periodic displacements. Explicit calculations are performed for three spheres interacting such that they form an equilateral triangle in the rest configuration.

###### pacs:
47.15.G-, 47.63.mf, 47.63.Gd, 47.63.M-

## I Introduction

The mechanism of swimming of a body immersed in an infinite viscous fluid is not yet fully understood. Much progress has been made in the study of the Stokes limit 1 (),2 (), pertaining to the swimming of microorganisms, with neglect of inertia of the body and the fluid. The effect of fluid inertia can be characterized by a dimensionless viscosity 3 (), or equivalently by a scale number 4 (), given by the ratio of a typical body dimension to the range of viscous dissipation during a period. The Stokes limit corresponds to infinite dimensionless viscosity and vanishing scale number. In the following we study a model swimmer in the full range of scale number. The Reynolds number is assumed to be small, so that vortex shedding is neglected and the fluid flow is laminar. In the swimming of fish and humans and in the flying of birds the Reynolds number is large and vortex shedding is relevant 5 (),6 (),6A ().

One can gain much insight by modeling the body as an assembly of spheres interacting directly, for example via harmonic springs, and with the effect of fluid motion embodied in frictional and added mass hydrodynamic interactions. The added mass of human bodies has been measured experimentally 7 (). In swimming the body is set in motion by actuating forces which vary periodically in time.

So far most model calculations in the Stokes limit have been made for linear chain structures with longitudinal motions along the chain. In particular, the linear three-sphere chain has been studied in detail 8 (),9 (). In earlier work we studied the effect of body and fluid inertia on the swimming of a linear three-sphere chain 10 (),10A (). In the following we consider the more complicated situation of planar assemblies with motions restricted to a plane. The particular case of a triangular assembly in the Stokes limit was studied by Vladimirov 11 ().

In Sec. II the dynamics of the resistive-reactive model is formulated on the basis of Hamiltonian equations of motion with added frictional forces. The frictional hydrodynamic interactions can in principle be calculated in the Stokes limit 12 (). In practice one uses an Oseen 13 () or Rotne-Prager approximation 14 (). The added mass hydrodynamic interactions can in principle be found from potential theory 15 (),16 (). In practice one uses a dipole approximation 10 (). The spheres are also subject to direct forces, as specified in an interaction potential, and to actuating forces, driving the motion of the assembly. The total actuating force and torque are assumed to vanish.

The translational and rotational swimming velocity are found from momentum and angular momentum balance equations for the assembly. In Sec. III these equations are derived from the equations of motion for the individual spheres. In addition we derive an instantaneous power equation 17 ().

For simplicity we restrict attention to small amplitude motion, which implies small Reynolds number. The latter is a combination of the amplitude and the scale number. In Secs. IV and V we derive expressions for the mean translational and rotational swimming velocities, valid to second order in the amplitude of displacements from the equilibrium positions. In Secs. VI and VII we apply the theory to an assembly of three spheres, with rest configuration in the form of an equilateral triangle. In Sec. VIII we calculate the power, averaged over a period, to second order in the amplitude of displacements and optimize the mean translational swimming velocity for given power with respect to the stroke. For short mean distances between the spheres the optimal efficiency of swimming of the triangular configuration is substantially larger than that of a comparable three-sphere linear chain.

## Ii Dynamics of spheres centered in a plane

We consider a system of identical spheres of radius and mass density immersed in a viscous incompressible fluid of shear viscosity and mass density . The fluid is of infinite extent in all directions. We assume that at all times the centers of the spheres are located in the plane of a Cartesian system of coordinates, so that the centers constitute a horizontal planar configuration. The dynamics of the system is governed by an interaction potential , depending on the instantaneous configuration of centers, by actuating forces , directed in the plane and summing to zero total force and zero total torque at any time , and by hydrodynamic forces exerted by the fluid. We assume that the spheres are freely rotating, so that the torque on any sphere vanishes. Further we assume that the hydrodynamic interactions can be approximated by Stokes friction, calculated from the Stokes equations, and by added mass effects, calculated from potential theory. We summarize the positions of centers in the -dimensional vector , and the sphere momenta in the -dimensional vector . The momenta are related to the velocities by

 p=m⋅U, (1)

where is the mass matrix, which depends on the relative positions of the spheres, so that it is invariant under translations of the whole assembly. Its tensor character determines the transformation under rotations of the assembly. The dynamics of the system is assumed to be governed by the approximate equations of motion 10 ()

 dRdt=U,dpdt=−∂K∂R−\boldmathζ⋅U−∂Vint∂R+E, (2)

where the kinetic energy is given by

 K=12p⋅w⋅p, (3)

with inverse mass matrix . The derivative with respect to positions in Eq. (2.2) is to be taken at constant momenta . The friction matrix and the interaction potential are invariant under translations of the assembly. The interaction potential is invariant under rotations and the friction matrix transforms according to its tensor character. We abbreviated . In spite of the fact that the total actuating force and torque vanish, the assembly can experience a net translation and rotation which is identified as translational and rotational swimming.

In the absence of actuating forces the system comes to rest due to friction with the fluid. The rest situation corresponds to a solution of Eq. (2.2) with constant configuration , which is a minimum of the potential energy . In the rest configuration the center and relative positions are

 \boldmathC0=1NN∑j=1\boldmathR0j,\boldmathc0j=\boldmathR0j−\boldmathC0,j=1,...,N. (4)

In shorthand notation . From the definitions it follows that

 uα⋅c0=0,(α=x,y), (5)

where and . With respect to fixed axes with unit vectors

 \boldmathc0j=bjx\boldmathex+bjy\boldmathey,j=1,...,N, (6)

where the constants are determined by the potential energy. We call the basic configuration. By isotropy a configuration rotated about the center by angle with vectors

 \boldmathcj=bjx\boldmatheρ+bjy\boldmath% eφ,j=1,...,N, (7)

with rotated axes

 \boldmatheρ = cosφ\boldmathex+sinφ\boldmath% ey, \boldmatheφ = −sinφ\boldmathex+cosφ% \boldmathey (8)

is also an equilibrium configuration.

The positions of the centers at time may be decomposed as

 \boldmathRj(t)=\boldmathC0+∫t0\boldmath% U(t′)dt′+\boldmathcj(t)+\boldmathdj(t),j=1,...,N, (9)

where is the center velocity at time , and is the additional displacement from the equilibrium structure. The instantaneous positions have center , given by the first two terms in Eq. (2.9). In shorthand notation

 R=C+c+d. (10)

From the definitions it follows that

 uα⋅d(t)=0,(α=x,y). (11)

The orientation of the equilibrium structure may be defined conveniently from the instantaneous positions . We choose to define it such that the squared distance is minimal for the orientation .

The time-derivative of Eq. (2.7) yields

 d\boldmathcjdt=Ω\boldmathez×%\boldmath$c$j,j=1,...,N, (12)

with angular velocity

 Ω=dφdt. (13)

In matrix notation this can be expressed as

 dcdt=−ΩX⋅c, (14)

 X=⎛⎜ ⎜ ⎜⎝0100−1000000100−10⎞⎟ ⎟ ⎟⎠. (15)

The generalization to arbitrary is obvious. With this notation the time-derivative of Eq. (2.10) reads

 U=Uβuβ−ΩX⋅c+˙d. (16)

Substituting this into Eq. (2.2) and requiring that the total actuating force and torque vanish we obtain equations of motion for , involving also the time-derivatives and of the displacements .

## Iii Momentum, angular momentum, and power

The requirement that the sum of actuating forces vanishes reads in abbreviated notation

 uα⋅E=0,(α=x,y). (17)

Similarly, the requirement that the total torque of actuating forces vanishes can be expressed as

 R⋅X⋅E=0. (18)

We can use these requirements to derive simple balance equations for the following components of total momentum and orbital angular momentum,

 Px=ux⋅p,Py=uy⋅p,Lz=R⋅X⋅p. (19)

These quantities vary due to interaction with the fluid. From Eq. (2.2) we derive

 dPxdt=−ux⋅\boldmathζ⋅U,dPydt=−uy⋅\boldmathζ⋅U,dLzdt=−R⋅X⋅\boldmathζ⋅U. (20)

The kinetic energy and potential energy terms in Eq. (2.2) do not contribute to these expressions on account of translational invariance and isotropy.

The orbital angular momentum can be decomposed as

 Lz = L′z+C⋅X⋅P,P=(Pxux+Pyuy)/N, L′z = (R−C)⋅X⋅p, (21)

where is the angular momentum relative to the center of mass. This has the rate of change

 dL′zdt=−(R−C)⋅X⋅% \boldmathζ⋅U−dCdt⋅X⋅P. (22)

The right-hand side depends only on velocities and relative positions.

By use of in Eq. (3.4) and substitution of Eq. (2.16) we derive coupled equations for the velocity components involving the time-derivatives and . The actuating forces do not occur explicitly in these equations and this allows a kinematic point of view in which the velocities are determined from the equations for prescribed displacements .

We define the Hamiltonian as

 H=K+Vint. (23)

From Eq. (2.2) we find for its time-derivative

 dHdt=−D+E⋅U, (24)

with rate of dissipation

 D=U⋅\boldmathζ⋅U. (25)

This may be called the instantaneous power equation 17 ().

In periodic swimming the time-average of the rate of dissipation over a period equals the power used. We denote the average as

 ¯¯¯¯D=1τ∫τ0D(t)dt, (26)

where is the period. From Eq. (3.8) we see that the mean rate of dissipation equals the power, i.e. the work performed by the actuating forces during a period,

 ¯¯¯¯D=¯¯¯¯¯¯¯¯¯¯¯E⋅U. (27)

In the same way we see from Eqs. (3.4) and (3.6)

 ux⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathζ⋅U=0,uy⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathζ⋅U = 0, ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(R−C)⋅X⋅% \boldmathζ⋅U+¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯dCdt⋅X⋅P = 0. (28)

The first two equations show that in periodic swimming the mean drag vanishes. In the next section we investigate the third equation for small amplitude motion.

## Iv Bilinear theory

In the following we consider an assembly with small deviations from an equilibrium configuration . The averaged Eqs. (3.12) are solved by formal expansion in powers of the displacements . We include terms up to second order. The displacements are assumed to vary harmonically in time at frequency .

To second order the first two equations (3.12) read

 uα⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathζ(1)⋅U(1)+uα⋅\boldmathζ0⋅¯¯¯¯¯¯¯¯¯U(2)=0,(α=x,y), (29)

where is the friction matrix of the equilibrium configuration, which is time-independent. Similarly, to second order the third equation reads

 c0⋅X⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathζ(1)⋅U(1)+¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(R(1)−C(1))⋅X⋅\boldmathζ0⋅U(1)+c0⋅X⋅\boldmathζ0⋅¯¯¯¯¯¯¯¯¯U(2)+¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯dC(1)dt⋅X⋅P(1)=0. (30)

In these equations we can put from Eq. (2.16)

 ¯¯¯¯¯¯¯¯¯U(2)=¯¯¯¯¯¯¯¯¯¯U(2)βuβ−¯¯¯¯¯¯¯¯¯Ω(2)X⋅c0, (31)

where we have used , as follows from

 c(1)=−φ(1)X⋅c0,Ω(1)=dφ(1)/dt. (32)

This shows that the second order mean velocities can be calculated from time-averaged products of first order quantities. We use the abbreviation

 ^c0=a−1c0⋅X=−a−1X⋅c0. (33)

From Eq. (4.1) we find

 Z0αβ¯¯¯¯¯¯¯¯¯¯U(2)β+Z0αca¯¯¯¯¯¯¯¯¯Ω(2)=¯¯¯¯¯¯¯¯¯I(2)Tα(α=x,y), (34)

with friction elements

 Z0αβ=uα⋅\boldmathζ0⋅uβ,Z0αc=uα⋅% \boldmathζ0⋅^c0, (35)

and mean second order translational impetus

 ¯¯¯¯¯¯¯¯¯I(2)Tα=−uα⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathζ(1)⋅U(1),(α=x,y). (36)

We write the second term in Eq. (4.2) as

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(R(1)−C(1))⋅X⋅% \boldmathζ0⋅U(1)=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(R(1)−C(1))⋅X⋅\boldmathζ0′⋅U(1)+ζ0¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(R(1)−C(1))⋅X⋅U(1), (37)

where

 \boldmathζ0′=\boldmathζ0−ζ0I,ζ0=6πηa, (38)

with unit matrix . We call the average in the last term in Eq. (4.9) the generalized ellipticity. It is a sum of ellipticities of the single particle displacements, but it also has pair contributions. We write the term in the form

 ζ0¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(R(1)−C(1))⋅X⋅U(1)=ζ0a2ω¯¯¯E=(1−α)Z0cca2¯¯¯¯¯¯¯¯¯Ω(2) (39)

with rotational friction coefficient

 Z0cc=^c0⋅\boldmathζ0⋅^c0. (40)

The coefficient and the second order rotational velocity are to be determined. The remaining terms in Eq. (4.2) yield

 Z0αc¯¯¯¯¯¯¯¯¯¯U(2)α+αZ0cca¯¯¯¯¯¯¯¯¯Ω(2)=¯¯¯¯¯¯¯¯¯I(2)R, (41)

with mean second order rotational impetus

 ¯¯¯¯¯¯¯¯¯I(2)R=−^c0⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯% \boldmathζ(1)⋅U(1)−a−1¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(R(1)−C(1))⋅X⋅\boldmathζ0′⋅U(1)−a−1¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯dC(1)dt⋅X⋅P(1). (42)

In the first term in Eq. (4.13) we have used the symmetry of the friction matrix . From Eqs. (4.6) and (4.13) we can solve for in terms of . Finally the coefficient can be determined from Eq. (4.11).

## V Mean translational and rotational impetus

In this section we study the mean translational and rotational impetus defined in Eqs. (4.8) and (4.14) in more detail. The first order friction matrix can be expressed as

 \boldmathζ(1)=aφ(1)^c0⋅∇\boldmathζ∣∣0+d⋅∇\boldmathζ% ∣∣0, (43)

where is the gradient operator in -dimensional configuration space. From Eq. (2.16) we find for the first order velocity vector and the corresponding position vector

 R(1)=C(1)βuβ+aφ(1)^c0+d,U(1)=U(1)βuβ+aΩ(1)^c0+˙d. (44)

Here we use first order equations derived from Eqs. (3.4) and (3.6). From Eq. (3.4) we find

with mass elements

 M0αβ=uα⋅m0⋅uβ,M0αc=uα⋅m0⋅^c0. (46)

From Eq. (3.6) we find

with mass elements

 M0cβ=^c0⋅m0⋅uβ=M0βc,M0cc=^c0⋅m0⋅^c0. (48)

These equations may be regarded as linear response equations determining the first order velocities in terms of the displacements . The equations can be solved by Fourier analysis. The equations for the complex Fourier coefficients read

 [−iωM0αβ+Z0αβ]U(1)βω+[−iωM0αc+Z0αc]aΩ(1)ω = uα⋅[ω2m0+iω% \boldmathζ0]⋅dω, [−iωM0cβ+Z0cβ]U(1)βω+[−iωM0cc+Z0cc]aΩ(1)ω = ^c0⋅[ω2m0+iω% \boldmathζ0]⋅dω. (49)

We introduce the shorthand notation and corresponding impedance vectors with

 fx,y(ω)=(−iωm0+\boldmathζ0)⋅ux,y,fc(ω)=(−iωm0+\boldmathζ0)⋅^c0. (50)

Then the solution of Eq. (5.7) can be expressed as

 ^U(1)ρω=iωYρσ(ω)fσ(ω)⋅dω, (51)

with admittance matrix , or alternatively as

 ^U(1)ρω=iω\boldmathΨρ(ω)⋅dω,\boldmathΨρ(ω)=Yρσ(ω)fσ(ω). (52)

The elements of the vector are dimensionless.

In the calculation of the time-averaged components of the impetus in Eqs. (4.8) and (4.14) we encounter bilinear expressions. The averages are evaluated conveniently in complex notation. For example

 ¯¯¯¯¯¯¯¯¯¯¯¯¯dU(1)β=12Reiωd∗ω\boldmathΨβ(ω)⋅dω. (53)

The leading contribution to the time-average of the translational impetus in Eq. (4.8) takes the form

 −uα⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(d⋅∇\boldmathζ)⋅˙d=12Re[iωd∗ω⋅Dα∣∣0⋅dω], (54)

with derivative friction matrix

 Dα=\boldmath∇fα,fα=\boldmathζ⋅uα=uα⋅\boldmathζ, (55)

as introduced earlier 9 (). We write the complete expression as

 ¯¯¯¯¯¯¯¯¯I(2)Tα=12Re[iωd∗ω⋅˘Dα∣∣0⋅dω], (56)

where the matrix differs from by corrections coming from the remaining terms in Eqs. (5.1) and (5.2). We have

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯φ(1)Ω(1)=0, (57)

so that there are four correction terms. It is convenient to use tensor notation and to define the -array as

 Fjkl=∂ζjl∂xk, (58)

where denotes the components of . We define the corresponding contractions

 Gαβk = uαjFjkluβl,Gαck=uαjFjkl^c0l, Hαl = uαj^c0kFjkl,Hcl=^c0j^c0kFjkl, (59)

where we use Einstein’s summation convention for latin indices. By comparison with Eq. (5.13)

 Gαβk=Dαkluβl,Gαck=Dαkl^c0l. (60)

The mean translational impetus can then be expressed as

 ¯¯¯¯¯¯¯¯¯I(2)Tα=[−Dαkl¯¯¯¯¯¯¯¯¯¯dk˙dl−Hαluβl¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯aφ(1)U(1)β+(Hαk−Gαck)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯aΩ(1)dk−Gαβk¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯dkU(1)β]∣∣∣0, (61)

where it is indicated that finally the coefficients must be evaluated at . The averages can be evaluated in complex notation. Besides Eqs. (5.11) and (5.12) we have

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯aφ(1)U(1)β = 12Re[−iωd∗ω⋅\boldmathΨ∗c(ω)\boldmathΨβ(ω)⋅dω], ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯aΩ(1)dk = 12Re[iωd∗kωΨcl(ω)dlω]. (62)

Substituting from Eq. (5.10) and using complex notation we obtain the matrix in Eq. (5.14) as

 ˘Dα(ω)=Dα+Hα⋅uβ\boldmathΨ∗c(ω)\boldmathΨβ(ω)−Gαβ\boldmathΨβ(ω)+(Hα−Gαc)\boldmathΨc(ω). (63)

Next we consider the time-average of the rotational impetus, defined in Eq. (4.14). In generalization of Eq. (5.13) we need the rotational derivative friction matrix given by

 Dc=\boldmath∇fc,fc=\boldmathζ⋅^c0=^c0⋅\boldmathζ. (64)

In analogy with Eq. (5.14) we have

 ¯¯¯¯¯¯¯¯¯I(2)R=12Re[iωd∗ω⋅˘Dc(ω)∣∣0⋅dω], (65)

with a matrix . From the first term in Eq. (4.14) we obtain correction terms as in Eq. (5.21) with the subscript replaced by . The corresponding -coefficients are

 Gcβk=^c0jFjkluβl=Gβck,Gcck=^c0jFjkl^c0l. (66)

The time-average of the second term in Eq. (4.14) is

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(R(1)−C(1))⋅X⋅% \boldmathζ0′⋅U(1)=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯aφ(1)^c0⋅X⋅\boldmathζ0′⋅U(1)+¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯d⋅X⋅% \boldmathζ0′⋅U(1). (67)

The first term on the right is evaluated to

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯aφ(1)^c0⋅X⋅% \boldmathζ0′⋅U(1)=−a−1c0⋅\boldmathζ0′⋅uβ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯aφ(1)U(1)β+a−1c0⋅\boldmathζ% 0′⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯aΩ(1)d, (68)

with use of and Eq. (5.15). The two averages are given by Eq. (5.20). The second term on the right in Eq. (5.25) is given by

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯d⋅X⋅\boldmathζ0′⋅U(1)=12Re[iωd∗kωXkm(f0′βmΨβl(ω)+f0′cmΨcl(ω)−ζ0′ml)dlω]. (69)

The time-average of the third term on the right in Eq. (4.14) vanishes on account of

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯dC(1)dt⋅X⋅P(1)=M0xy12NRe[ω2d∗kω(Ψ∗xk(ω)Ψyl(ω)−Ψ∗yk(ω)Ψxl(ω))dlω]=0. (70)

Collecting terms we find that the matrix in Eq. (5.23) is given by

 ˘Dc(ω) = Dc+(Hc−a−2c0⋅% \boldmathζ0′)⋅uβ\boldmathΨ∗c(ω)\boldmathΨβ(ω)+a−1X⋅\boldmathζ0′ (71) − (Gcβ+a−1X⋅f0′β)\boldmathΨβ(ω)+(Hc−a−2c0⋅\boldmathζ0′−Gcc−a−1X⋅f0′c)\boldmathΨc(ω).

The second order mean swimming velocities are calculated from the solution of Eqs. (4.6), (4.11) and (4.13). The expressions are complicated, but in our application they simplify because many of the coefficients vanish by symmetry. We remark that the expressions in Eqs. (5.21) and (5.29) vanish in the absence of frictional hydrodynamic interactions. In that case the -array in Eq. (5.16) is identically zero and the matrix vanishes.

Finally we note that in Eqs. (5.14) and (5.23) the matrices and can be reduced to their hermitian part. Also we can use projectors to take account of the fact that in swimming the displacement vector must be orthogonal to , and .

## Vi Triangular assembly

As an application of the equations derived in the preceding sections we consider swimming of a triangular assembly of spheres of radius and mass density . For simplicity we consider neutrally buoyant spheres. As basic equilibrium configuration we consider the equilateral triangle given by

 c0=(−d2√3,−d2,d√3,0,−d2√3,d2). (72)

The center of this triangle is at the origin and the three sides have length . We have chosen the orientation such that with suitable displacements the swimming is in the direction. The basic vectors and are

 ux = (1,0,1,0,1,0),uy=(0,1,0,1,0,1), ^c0 = da(12,−12√3,0,1√3,−12,−12√3). (73)

We take the potential energy to be given by the expression

 Vint(R)=kd2[(\boldmathr12⋅% \boldmathr12−d2)2+(\boldmathr23⋅\boldmathr% 23−d2)2+(\boldmathr31⋅\boldmathr31−d2)2], (74)

with elastic constant and relative distance vectors

 \boldmathr12=\boldmathR2−\boldmathR1,\boldmathr23=\boldmathR3−\boldmathR2,\boldmathr31=\boldmathR1−\boldmathR3. (75)

The potential energy is positive definite, translation-invariant, isotropic, and it vanishes at configuration .

We assume that the mobility matrix is given by Oseen hydrodynamic interactions 13 (). Hence the friction matrix is evaluated to first order in the ratio . We assume that the inverse mass matrix is evaluated in dipole approximation 10 (). Hence the mass matrix is evaluated to first order in the ratio . The kinetic energy is positive definite, translation-invariant, and isotropic. The friction matrix and the mass matrix depend only on the relative distance vectors given by Eq. (6.4).

We introduce projection operators and defined as

 Pop=13uxux+13uyuy+^c0^c0^c0⋅^c0,Q=I−Pop. (76)

The displacement vector must satisfy

 dω=Q⋅dω, (77)

to exclude rigid body motion. The projected vector has only three independent components. This allows reduction of the matrices and in Eqs. (5.21) and (5.29) to a three-dimensional representation.

In order to construct the reduced matrices we expand the displacement vector in terms of a convenient set of basis vectors. We use the orthonormal set of eigenvectors of the elasticity matrix corresponding to the interaction energy . To second order in deviations from equilibrium we find from Eq. (6.3)

 Vint2=12(R−R0)⋅H⋅(R−R0), (78)

with elasticity matrix given explicitly by

 H=k⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝62√3−6−2√3002√310−2√3−20−8−6−2√3120−62√3−2√3−2042√3−200−62√36−2√30−82√3−2−2√310⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (79)

This has the orthonormal set of eigenvectors

 e1 = (1,0,1,0,1,0)/√3=ux/√3, e2 = (0,1,0,1,0,1)/√3=uy/√3, e3 = (12,−12√3,0,1√3,−12,−12√3)=ad^c0, e4 = (12,−12√3,−12,−12√3,0,1√3), e5 = (−12√3,−12,−12√3,12,1√3,0), e6 = (−12√3,−12,1√3,0,−12√3,12)=1dc0,