Swimming of a linear chain with a cargo in an incompressible viscous fluid with inertia

Swimming of a linear chain with a cargo 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

An approximation to the added mass matrix of an assembly of spheres is constructed on the basis of potential flow theory for situations where one sphere is much larger than the others. In the approximation the flow potential near a small sphere is assumed to be dipolar, but near the large sphere it involves all higher order multipoles. The analysis is based on an exact result for the potential of a magnetic dipole in the presence of a superconducting sphere. Subsequently, the approximate added mass hydrodynamic interactions are used in a calculation of the swimming velocity and rate of dissipation of linear chain structures consisting of a number of small spheres and a single large one, with account also of frictional hydrodynamic interactions. The results derived for periodic swimming on the basis of a kinematic approach are compared with bilinear theory, valid for small amplitude of stroke, and with the numerical solution of the approximate equations of motion. The calculations cover the whole range of scale number between the friction-dominated Stokes limit and the inertia-dominated regime.

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

I Introduction

In earlier work we have studied the swimming of linear chain structures, consisting of a number of small spheres and a single large one, in the Stokes limit where inertia of the spheres and the fluid is neglected 1 (). A model of this kind with two little spheres pushing a big one was first studied by Golestanian 2 (). It is of interest to investigate the effects of inertia on the swimming of bodies immersed in an incompressible viscous fluid. Many microorganisms are so small that the Stokes limit provides a valid approximation, but for larger bodies inertia becomes important. The corresponding reactive effects are caused by the dependence of the added mass of the whole structure on its instantaneous configuration.

In models where the structure is approximated as an assembly of spheres the added mass hydrodynamic interactions are manifested as a dependence of the mass matrix on the instantaneous positions of centers. In principle the dependence can be found from the theory of potential flow. Elsewhere we have used a dipole approximation to calculate the mass matrix of an assembly of spheres 3 ()4 (). The approximation is useful when all spheres are far apart, but fails for configurations where some small spheres are close to a large one. For such configurations an approximate method of calculation of the Stokes friction matrix was proposed by Ekiel-Jeżewska and Felderhof 5 (). In the following we derive a similar approximation to the mass matrix on the basis of an exact expression for the field of a magnetic point dipole in the presence of a superconducting sphere, as derived by Palaniappan 6 (). The corresponding scalar potential can be used in hydrodynamics if account is taken of convection of the spheres by the fluid. In Sec. II of this article we present the derivation of the approximate mass matrix in some detail.

Subsequently we apply the approximation in a study of the longitudinal swimming of linear chain structures with a single large sphere and a number of small ones. As an example, we calculate in Sec. III the mean swimming velocity and power for swimming of a linear chain of three small spheres and one large one on the basis of kinematics, taking into account both friction and added mass effects 7 (). We assume that the stroke is proportional to the one that is optimal at small amplitude in the Stokes limit, as calculated earlier 1 ().

For a linear chain of three identical spheres the optimal stroke, as derived from bilinear theory valid to second order in the amplitude of swimming, is nearly independent of the scale number which characterizes the effect of fluid inertia 4 (). In Sec. IV we find that the same is true for chains of three small spheres and an additional large one. We show also that the bilinear theory provides a good approximation to the results derived from the kinematic approach of Sec. III.

In Sec. V we calculate the corresponding actuating forces for a chain with harmonic elastic interactions and compare with the swimming of the chain with application of a cargo constraint, implying that the actuating force on the large sphere vanishes. We use bilinear theory to optimize the actuating forces subject to the cargo constraint. It turns out that the latter does not lead to a significant decrease of swimming efficiency compared to the optimum value without the constraint. Finally, we study transient effects in swimming by solving the approximate equations of motion for the spheres numerically, for a chain subject to actuating forces with the cargo constraint, starting from the rest situation.

Ii Added mass in small particle approximation

As a preliminary we study in this section added mass hydrodynamic interactions for situations where one sphere is much larger than the others. We consider an -body system consisting of a sphere of radius , labeled , and much smaller particles of radii , all immersed in a viscous incompressible fluid of shear viscosity and mass density . The fluid is of infinite extent in all directions. The whole system is considered to be at rest for . At time the system is caused to move by instantaneous forces of short duration

 \boldmathEj(t)=\boldmathSjδ(t)j=0,1,...,N, (1)

where is the impulse imparted to the sphere or particle. In order to satisfy the kinematic boundary condition that sphere and particles are impenetrable, the fluid moves with a flow velocity satisfying Laplace’s equation,

 \boldmathv=−∇ϕ,∇2ϕ=0. (2)

The hydrodynamic interactions are embodied in the inverse mass matrix which relates the translational velocities to the impulses imparted to the bodies via 4 ()

 \boldmathUj=N∑k=0wjk⋅\boldmathSk,(j=0,...,N), (3)

where each tensor depends on the positions of all centers. Correspondingly the mass matrix also depends on all center positions. We shall derive approximate simplified expressions for the inverse mass tensors, based on an approximation to the irrotational flow pattern and the corresponding convective effects. In the approximation the element is independent of position, elements and with depend only on via a pair hydrodynamic interaction, and elements , with depend only on and via a three-body hydrodynamic interaction. In order to derive the expression for the pair hydrodynamic interactions and it suffices to consider the sphere and a single small particle. In order to derive the three-body hydrodynamic interaction we must consider the sphere and two small particles.

We use Cartesian coordinates with origin at the center of the large sphere. We recall first that in the absence of particles, the sphere made to move with translational velocity generates a potential 8 ()

 ϕ(\boldmathr)=12b3\boldmathrr3⋅\boldmathU0,r>b, (4)

corresponding to the dipole moment

 \boldmathq0=12b3\boldmathU0. (5)

The corresponding Poisson flow pattern is given by

 \boldmathvP(\boldmathr)=F0(\boldmathr% )⋅\boldmathq0,F0(\boldmathr)=−I+3^\boldmathr^\boldmathrr3, (6)

where is the dipolar tensor with unit tensor . The dipole moment is related to the impulse by

 \boldmathq0=β0\boldmathS0,β0=b32m∗0, (7)

where is the effective mass

 m∗0=m0+12mf0,mf0=4π3ρb3. (8)

Here is the mass of displaced fluid and is the added mass.

Next we consider a sphere of radius and a single small particle, both immersed in the fluid. The pair hydrodynamic interaction between particle and sphere is embodied in the six-dimensional inverse mass matrix relating the sphere velocity and the particle velocity to the impulses and imparted to the sphere and particle according to

 \boldmathU0 = w00⋅\boldmathS0+w01⋅\boldmathS1, \boldmathU1 = w10⋅\boldmathS0+w11⋅\boldmathS1. (9)

We consider first the situation where an impulse is imparted to the sphere, but . In our approximation the particle is subjected to the dipolar flow pattern generated by , but its reaction is neglected, except for its convective motion. Hence the parts and of the inverse mass matrix are given simply by

 w00=1m∗0I,w10=γ1β0F0(\boldmathr1), (10)

with convective coefficient

 γ1=3mf12m∗1=4πρβ1, (11)

and relative vector . The convective coefficient with and follows from linear response theory 9 (). The value is consistent with the derivations of Landau and Lifshitz 10 () and of Batchelor 11 ().

In order to derive an approximate expression for the part of the inverse mass matrix in Eq. (2.9) we consider a flow situation where an impulse is applied to the particle, but the sphere is convected such that it exerts no force on the fluid. In the absence of the sphere and in infinite fluid the impulse would generate a dipolar flow with dipole moment . We regard this as an unperturbed flow acting on the freely moving sphere. The resulting sphere velocity is 9 ()

 \boldmathU01=γ0¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯F0(\boldmathr% −\boldmathr1)V⋅\boldmathq1, (12)

where the overline indicates an average over the volume of the sphere. Here

 γ0=3mf02m∗0=4πρβ0, (13)

as in Eq. (2.11). The average in Eq. (2.12) corresponds to the field at due to a uniformly polarized sphere. This is dipolar in , so that we have

 \boldmathU01=γ0F0(\boldmathr1)⋅\boldmathq1, (14)

corresponding to

 w01=γ0β1F0(\boldmathr1). (15)

Comparing with Eq. (2.10) we see that , since . There is symmetry, as expected on general grounds.

In order to derive an approximation to the part of the inverse mass matrix we consider the flow generated by the sphere when acted upon by the dipolar flow due to the particle. First we recall an important result derived by Palaniappan 6 () for the field of a magnetic dipole near a superconducting sphere. For a fixed sphere centered at the origin and with a point dipole acting at he derived the scalar potential 6 ()

 Φ(\boldmathr,\boldmathr1)=Φ0(\boldmathr−\boldmathr1)+ΦR(\boldmathr,\boldmathr1), (16)

where is the potential for infinite fluid in the absence of the sphere,

 Φ0(\boldmathr−\boldmathr1)=(\boldmathr−\boldmathr1)⋅\boldmathq1|\boldmathr−\boldmathr1|3, (17)

and is the reflection potential

 ΦR(\boldmathr,\boldmathr1) = b3r31(\boldmathr−¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathr1)⋅(\boldmathq1,⊥−\boldmathq1,∥)¯¯¯d3 (18) − 1br1[r2−(\boldmathr⋅\boldmathr% 1/r1)2](r−\boldmathr⋅(\boldmathr−¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathr1)¯¯¯d)\boldmathr⋅\boldmathq1,⊥,

where the image point is defined by

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathr1=b2r21\boldmathr1, (19)

and is the distance from the field point to the image point,

 ¯¯¯d=|\boldmathr−¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathr1|. (20)

Furthermore, the vector

 \boldmathq1,∥=\boldmathr1\boldmath% r1r21⋅\boldmathq1 (21)

is the part of the dipole moment parallel to and is the part perpendicular to . The potential satisfies Laplace’s equation , as well as the boundary condition

 ∂Φ∂r=0onr=b. (22)

Hence the corresponding flow is tangential to the spherical surface. The flow can be expressed as

 \boldmathv1(\boldmathr)=−∇Φ(\boldmathr,\boldmathr1)=F(\boldmathr,\boldmathr1)⋅\boldmathq1, (23)

with Green tensor . The Green tensor can be decomposed as

 F(\boldmathr,\boldmathr1)=F0(% \boldmathr−\boldmathr1)+FR(\boldmathr,\boldmathr1), (24)

where is the reflection tensor. The tensors have the symmetry

 F0(\boldmathr−\boldmathr1)=F0(\boldmathr1−\boldmathr)T,FR(% \boldmathr,\boldmathr1)=FR(\boldmathr1,\boldmathr)T, (25)

where the superscript indicates the transpose.

In order to obtain the flow generated by in the presence of the freely moving sphere we must add the dipolar flow corresponding to the sphere velocity given by Eq. (2.14). The total flow can be expressed as

 ^\boldmathv1(\boldmathr)=^F(% \boldmathr,\boldmathr1)⋅\boldmathq1, (26)

with Green tensor

 ^F(\boldmathr,\boldmathr1)=F(%\boldmath$r$,\boldmathr1)+V(\boldmathr,% \boldmathr1) (27)

with tensor function given by

 V(\boldmathr,\boldmathr1)=12b3γ0F0(\boldmathr)⋅F0(% \boldmathr1). (28)

The tensor has the symmetry

 V(\boldmathr,\boldmathr1)=V(% \boldmathr1,\boldmathr)T. (29)

The total image dipole of the fixed sphere is

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathq1=−12b3F0(% \boldmathr1)⋅\boldmathq1, (30)

as one sees from the long range behavior of the reflection potential in Eq. (2.18). For a neutrally buoyant sphere , and then this is precisely canceled by the dipole corresponding to the flow caused by sphere velocity . More generally, the tensor tends to

 ^FR(\boldmathr,\boldmathr1)≈(γ0−1)12b3F0(\boldmathr)⋅F0(\boldmathr1),(r>>b) (31)

at large distance from sphere and particle.

We derive an approximate expression for the inverse mass tensor by considering the reflected part of the modified flow at the point . This yields

 w11=1m∗1I+β1γ1^FR(\boldmathr1,\boldmathr1). (32)

From Eqs. (2.18) and (2.27) we find the explicit expression

 ^FR(\boldmathr1,\boldmathr1) = −2b3(r21−b2)3^\boldmathr1^\boldmathr1−b2+r212r21b3(r21−b2)3(I−^\boldmathr1^% \boldmathr1) (33) + 2b3γ0r61^\boldmathr1^\boldmathr1+b3γ02r61(I−^\boldmathr1^\boldmathr1).

The derivation of the first two terms is best performed in spherical coordinates. The last two terms depend on the mass of the sphere via . The expression satisfies Eq. (2.31), since

 F0(\boldmathr1)⋅F0(\boldmathr1)=4r61^\boldmathr1^\boldmathr1+1r61(I−^\boldmathr1^\boldmathr1). (34)

Note that becomes singular at short distance and shows a rapid decay at large distance. The tensor is obviously symmetric. In our approximation it involves only hydrodynamic interactions between particle 1 and the big sphere. It can be shown that the above results are consistent with the dipole approximation 3 ()4 () and with the calculation of the kinetic energy of irrotational flow for two moving spheres quoted by Lamb 12 ().

We can use the same formalism to derive the approximate expression for the inverse mass tensor , with . It suffices to consider a sphere and two small particles located at and . The linear relation between velocities and impulses becomes

 \boldmathU0 = w00⋅\boldmathS0+w01⋅\boldmathS1+w02⋅\boldmathS2, \boldmathU1 = w10⋅\boldmathS0+w11⋅\boldmathS1+w12⋅\boldmathS2, \boldmathU2 = w20⋅\boldmathS0+w21⋅\boldmathS1+w22⋅\boldmathS2. (35)

Most parts of the inverse mass matrix are given by expressions derived above. Only the parts and require further consideration. In generalization of Eq. (2.32) we have

 w12=γ1β2^F(\boldmathr1,\boldmathr2),w21=γ2β1^F(\boldmathr2,\boldmathr1), (36)

with and . The symmetry relation

 w21=wT12 (37)

is satisfied, as follows from Eqs. (2.25) and (2.29), in agreement with general arguments.

The expression for given by Eq. (2.36) is complicated, but it simplifies for configurations for which and are collinear with the origin, so that . For such configurations we find

 w12 = γ1γ24πρ[F0(%\boldmath$r$1−\boldmathr2)−2b3(r1r2−b2)3^\boldmathr1^\boldmathr1 − b2+r1r22r1r2b3(r1r2−b2)3(I−^\boldmathr1^\boldmathr1)+b3γ02r31r32(I+3^% \boldmathr1^\boldmathr1)].

The first term is the pair hydrodynamic interaction between two small particles. The last three terms depend on both and , and therefore represent a three-body hydrodynamic interaction. Note that the interaction becomes singular as both and tend to the sphere radius . It decays as as both and tend to infinity.

The expressions derived above can be used as approximations in the many-body inverse mass matrix in Eq. (2.3). If we regard as a typical small particle radius, as a typical distance between small particles, and as the minimum separation distance of the center of a small particle from the surface of the sphere, then the ratios , and may be regarded as small parameters. It follows from a consideration of the multipole expansion of the exact mass matrix that our expressions represent the first few terms in a systematic expansion in powers of the small parameters. We call our expressions the small particle approximation to the inverse mass matrix.

Iii Kinematic swimming

As an application of the added mass hydrodynamic interactions derived in the preceding section we consider swimming of a linear chain structure consisting of a big sphere of radius and three small spheres of radius with periodic motions of the four bodies along the axis of a Cartesian system of coordinates. In Fig. 1 we show a sketch of the structure. The small spheres have centers at with and the center of the big sphere is at with . In this section we use a kinematic approach and prescribe the periodic relative motion of the four spheres, where . For given relative motion of the spheres the asymptotic periodic swimming velocity and the periodic rate of dissipation are determined by the added mass and frictional hydrodynamic interactions. We use a mass matrix as derived in the preceding section, and a friction matrix as derived in earlier work with Ekiel-Jeżewska 5 ().

The asymptotic periodic swimming velocity can be calculated from the periodic total mass and total friction coefficient by use of an expression derived earlier 7 (). From the given relative motion and the calculated swimming velocity the periodic sphere velocities in the laboratory frame can be evaluated. The velocities of the individual spheres are a sum of the swimming velocity and a displacement velocity. The latter is a component of the vector

 ˙d=T−1⋅(0,drdt), (39)

where is the matrix relating center and relative coordinates to the Cartesian coordinates . In the present case this is given explicitly by

 T=⎛⎜ ⎜ ⎜ ⎜⎝14141414−11000−11000−11⎞⎟ ⎟ ⎟ ⎟⎠. (40)

The periodic rate of dissipation then follows from the expression , where , with , is the four-vector comprising the -components of the sphere velocities. It follows from Eq. (3.1) that .

The reduced speed is defined by , where the overhead bar indicates the time average over a period , and the reduced rate of dissipation is defined by . For small amplitude swimming we denote the reduced quantities as and , where is an amplitude factor and the subscript 2 indicates that the calculation is performed to second order in the amplitude. We study these quantities as functions of the dimensionless scale number defined by 13 ()

 s2=a2ωρ2η. (41)

For the swimming is dominated by added mass effects, whereas corresponds to the Stokes limit, dominated by friction. The efficiency of the stroke is defined by

 ET=^U^D. (42)

The time-dependent periodic swimming velocity and rate of dissipation are calculated from an expansion in Fourier series 7 ().

In earlier work we have studied the swimming of a collinear chain in the Stokes limit 1 () and determined the optimal stroke for small amplitude swimming of a chain of harmonically linked spheres with equilibrium relative distances . The optimal stroke is independent of the strength of direct interactions, since it is determined solely by kinematics. The optimal stroke follows from an eigenvalue problem with a matrix determining and a matrix determining , where is the complex vector of relative displacements. The maximum eigenvalue of the eigenvalue problem corresponds to the optimal ratio . For the case we found with corresponding complex eigenvector

 ^\boldmathξc0=(0.720,−0.121+0.598i,−0.311−0.112i). (43)

The vector is normalized to unity. The corresponding reduced speed and power are

 ^U2=0.000792,^D2=8.523,(s=0). (44)

The eigenvector determines the optimal stroke in the Stokes limit. The relative motion during the optimal stroke is given by

 r(t)=r0+εaRe[^\boldmath% ξc0e−iωt]. (45)

with amplitude factor . In Fig. 2 we show the motion of the four spheres for in the center system, where the geometrical center of the four spheres is at rest. For comparison, the maximum eigenvalue of the 3-chain without the big sphere in the Stokes limit is , as given by Eq. (6.5) in Ref. 14. The corresponding values of reduced speed and power are and .

For the system with inertia we use the same stroke, given by Eq. (3.5). We have found for the chain with three identical spheres 4 () that the optimal stroke at small amplitude varies little as a function of scale number , and we shall show in Sec. IV that the same is true for the present system.

We assume that the four spheres are neutrally buoyant. In particular we then find numerically for and

 ^U=0.00328,^D=34.381,(s=1,ε=2), (46)

corresponding to efficiency . For and we find

 ^U=0.00317,^D=34.381,(s2=10,ε=2), (47)

corresponding to efficiency . For comparison, in the Stokes limit and to second order in the amplitude

 ε2^U2=0.00317,ε2^D2=34.092,(s=0,ε=2), (48)

corresponding to efficiency .

The Fourier coefficients of the swimming velocity decrease rapidly with increasing order. For the above case with we find for the first five absolute ratios

 {1,∣∣∣Usw,1¯¯¯¯¯¯¯¯¯Usw∣∣∣,∣∣∣Usw,2¯¯¯¯¯¯¯¯¯Usw∣∣∣,∣∣∣Usw,3¯¯¯¯¯¯¯¯¯Usw∣∣∣,∣∣∣Usw,4¯¯¯¯¯¯¯¯¯Usw∣∣∣}={1,73.80,10×10−4,4×10−4,10−5}. (49)

For the case with we find

 {1,∣∣∣Usw,1¯¯¯¯¯¯¯¯¯Usw∣∣∣,∣∣∣Usw,2¯¯¯¯¯¯¯¯¯Usw∣∣∣,∣∣∣Usw,3¯¯¯¯¯¯¯¯¯Usw∣∣∣,∣∣∣Usw,4¯¯¯¯¯¯¯¯¯Usw∣∣∣}={1,77.76,5×10−4,4×10−4,1×10−4}. (50)

In both cases the swimming velocity is essentially just the sum of the mean and the first harmonic, and the amplitude of the first harmonic is larger than the mean.

Iv Bilinear theory

Next we study small amplitude swimming in a fluid with inertia. We consider collinear spheres with dynamics dominated by frictional and added mass hydrodynamic interactions. The bilinear theory yields results for mean swimming velocity and power valid to second order in the amplitude of the stroke. These quantities are evaluated from matrices and operating in the space of relative displacements as

 ¯¯¯¯¯¯¯¯¯¯U(2)=12ωa(\boldmathξc|B|\boldmathξc),¯¯¯¯¯¯¯¯¯D(2)=12ηω2a3(\boldmathξc|A|\boldmathξc), (51)

with a typical sphere radius . The -dimensional vector and the matrices and have dimensionless elements. The expressions given earlier 4 () for the matrices need correction, as detailed below. The numerical calculations in the earlier work 4 () were performed with the correct expressions.

It was shown earlier 4 () that for harmonic displacements the second order mean swimming velocity is given by

 ¯¯¯¯¯¯¯¯¯¯U(2)=12Re[iωd∗ω⋅V(ω)∣∣0⋅dω], (52)

with frequency-dependent matrix given by

 V(ω)=μ˘D(ω), (53)

with mobility , where , and reduced derivative friction matrix

 ˘D(ω)=D−Y(ω)gf(ω), (54)

 Y(ω)=[−iωM+Z]−1, (55)

and vectors

 g=D⋅u,f(ω)=(−iωm+\boldmathζ)⋅u. (56)

Here , and the derivative friction matrix is given by

 D=\boldmath∇f,f=\boldmath% ζ⋅u, (57)

with gradient operator . In Eq. (4.2) the matrix is calculated in the time-independent equilibrium configuration. In the expression the product can be replaced by its hermitian part

 [iωV(ω)]h=iωVa,Va=12(V−V†). (58)

Similarly the second order mean rate of dissipation is given by

 ¯¯¯¯¯¯¯¯¯D(2)=12ω2Re[d∗ω⋅P(ω)∣∣0⋅dω], (59)

with the matrix

 P(ω)=−iωm+\boldmathζ−Y(ω)f(ω)f(ω). (60)

In Eq. (4.9) the matrix can be replaced by its hermitian part .

The matrices and have the properties

 u⋅V(ω) = 0,V(ω)⋅u=0, u⋅P(ω) = 0,P(ω)⋅u=0, (61)

and this allows to reduce the calculation to one in relative space. We define the transformed matrices

 VaT=T⋅Va⋅T−1,PhT=T⋅Ph⋅T−1, (62)

where is the matrix relating center and relative coordinates in generalization of Eq. (3.2). The first row and column of the matrices and vanish identically on account of the properties Eq. (4.11). The truncated matrices obtained by erasing the first row and column are denoted as and . Finally, the dimensional matrices and in Eq. (4.1) are given by

 B=iaCT⋅^VaT,A=1ηaCT⋅^PhT, (63)

with the matrix

 CT=[(T−1)T⋅T−1]\boldmath^. (64)

This dimensional matrix consists of numerical coefficients and is obtained from the corresponding matrix by truncation, as indicated by the final hat symbol.

In the present context the bilinear theory can be used for two purposes. On the one hand we can study the optimal stroke, as determined from the eigenvalue problem

 B|\boldmathξc)=λA|\boldmathξc), (65)

as the solution with maximum eigenvalue , and show that for neutrally buoyant spheres this depends only weakly on the scale number , defined in Eq. (3.3). On the other hand, for a calculated or assumed periodic relative motion we can obtain a quick estimate of mean swimming velocity and power as functions of the scale number. For example, we can use the stroke given by Eq. (3.5), or the long-time nearly periodic motion found from the numerical solution of the equations of motion with cargo constraint in Sec. V, approximating this by its first harmonic. Furthermore, it is of interest to show that the bilinear theory provides a good approximation over a wide range of amplitudes.

We consider first the eigenvalue problem Eq. (4.15). In Fig. 3 we plot the maximum eigenvalue for and neutrally buoyant spheres as a function of . The function is nearly constant, but shows a minimum at . At the minimum . The corresponding eigenvector is

 ^\boldmathξc0=(0.719,−0.127+0.596i,−0.310−0.122i). (66)

The eigenvalue tends to as . The corresponding eigenvector is

 ^\boldmathξc0=(0.722,−0.123+0.595i,−0.310−0.114i). (67)

These values should be compared to Eq. (3.5) corresponding to . As can be seen, the eigenvalue and eigenvector hardly vary over the whole range of . This justifies the use of the Stokes eigenvector Eq. (3.5) in the whole range.

Finally we compare the values of mean swimming velocity and power for the bilinear theory in the inertial regime with those of the exact calculation. The values for the bilinear theory, as calculated for the Stokes eigenvector Eq. (3.5), are

 ε2^U2=0.00320,ε2^D2=34.394,(s→∞,ε=2), (68)

corresponding to efficiency . These values should be compared with Eq. (3.9). The bilinear theory provides a good estimate even at this large amplitude factor.

V Swimming with the cargo constraint

In the calculations of Sec. III the relative motion of the spheres was prescribed. The asymptotic periodic swimming velocity was derived as a solution of the equation of motion

 ddt(MU)+ZU=I (69)

with time-dependent mass and friction coefficient given by

 M=u⋅m⋅u,Z=u⋅% \boldmathζ⋅u,u=(1,...,1), (70)

and impetus given by

 I(t)=−ddt(u⋅m⋅˙d)−u⋅\boldmathζ⋅˙d, (71)

where , given by Eq. (3.1), is the time-derivative of the displacement vector . For collinear spheres with positions and velocities the center velocity is defined as . The equation of motion for the center velocity Eq. (5.1) was derived from the equations of motion for the individual spheres, which read 4 ()

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

with momenta and . The partial derivative is taken at constant momenta . Furthermore is the potential of direct interaction forces. The actuating forces are summarized in . These oscillate at frequency and sum to zero at any time. The cargo constraint implies that the actuating force on the big sphere vanishes. In our case =0, and there are only two independent actuating forces , with .

We assume the direct interaction forces to be harmonic of the form

 −∂Vint∂R=H⋅(R−R0), (73)

with a real and symmetric matrix with the property . Specifically we choose

 H=k⎛⎜ ⎜ ⎜⎝−11001−21001−32002−2⎞⎟ ⎟ ⎟⎠ (74)

with elastic constant . This corresponds to nearest neighbor interactions between the three small spheres and a twice as strong link between the last small sphere and the big sphere. The stiffness of the assembly is characterized by the dimensionless number defined by . The elastic forces must be chosen such that stability in the numerical calculations is ensured. We choose and .

For given actuating forces, oscillating at frequency and with the property , we can solve the first order equations of motion,

 dR(1)dt=U(1),m0⋅dU(1)dt=−\boldmathζ0⋅U(1)+H⋅R(1)+E, (75)

as well as the linearized version of Eq. (5.1), to obtain the oscillating displacement with the property . The cargo constraint limits the class of displacements. The question arises how to select actuating forces, satisfying the cargo constraint, such that the efficiency is maximal. In the bilinear theory we can optimize by solving the eigenvalue problem in a reduced subspace.

Assuming harmonic time-dependence we find from the linearized version of Eq. (5.1) for the amplitude of the first order swimming velocity 4 ()

 U(1)ω=[−iωM0+Z0]−1u⋅(ω2m0+iω\boldmathζ0)⋅dω. (76)

The first order velocities are . Using Eq. (3.1) we find from Eq. (5.7) that the relation is satisfied automatically, and we obtain a linear relation between the amplitude vector in relative space and the force amplitude vector .

In our case the force amplitude vector is assumed to have the form , and substitution of the relation into and leads to a two-dimensional reduction of the speed matrix and the power matrix, acting in the space of complex vectors . The eigenvector with maximum eigenvalue of the two-dimensional eigenvalue problem