A theoretical expression is derived for the mechanical contribution to the mean swimming speed of a vibrating two-sphere, consisting of two spheres connected by an elastic spring, and immersed in a viscous incompressible fluid. The spring provides a harmonic potential for oscillations about a mean distance between centers. The system is made to oscillate at a chosen frequency by activating forces which sum to zero. The mechanical contribution to the resulting mean swimming velocity is calculated from the mechanical equations of motion and the corresponding impedance matrix of linear response. The frequency-dependent pair friction coefficients are found from approximate expressions derived earlier. The mechanical contribution is calculated to second order in the amplitude of stroke as a function of the scaling number, a dimensionless combination of size, frequency, and kinematic viscosity. Retarded friction and added mass determine the functional behavior.

The swimming velocity of a deformable body in a fluid can be decomposed as a sum of two contributions. The first contribution follows from the equations of classical mechanics' for the individual body parts. These equations involve friction and added mass effects. The second contribution follows from the Navier–Stokes equations for fluid velocity and pressure. Both contributions are a result of rectification due to nonlinear interactions. We call the first contribution mechanical and the second one hydrodynamic. In the present article, we study the first mechanical contribution for a vibrating two-sphere, consisting of two spheres, connected by an elastic spring, and immersed in a viscous incompressible fluid. A picture of a two-sphere is shown in Fig. 1.

FIG. 1.

Picture of a typical two-sphere.

FIG. 1.

Picture of a typical two-sphere.

Close modal

The work of Klotsa et al.1 established the vibrating two-sphere as an important model in the theory of swimming. The model is surprisingly complex but sufficiently simple to allow some mathematical analysis. Its partial solution can improve understanding and elucidate principles. This will be helpful in the analysis of more sophisticated swimmers.

In the case of the two-sphere, the body parts are the two spheres. The pair mobilities can be evaluated with various degrees of accuracy.2 Hubert et al.3 calculated the mean swimming speed of a two-sphere with Oseen type interactions for large distance between centers. For small amplitude of stroke, we generalized this result to shorter distance4 and included added mass effects in dipole approximation. We used a kinematic approach in which the relative distance between centers is prescribed to oscillate harmonically. The mean speed was obtained as a function of the oscillation frequency.

In the following, we study the two-sphere chain with harmonic elastic interaction, set into motion by periodic activating forces, which sum to zero. First, we consider instantaneous hydrodynamic interactions without specifying the distance dependence and derive a novel expression for the mean speed of small amplitude swimming in terms of the impedance matrix of linear response theory. For Oseen type hydrodynamic interactions, with or without dipolar added mass effects, the expression for the mean speed reduces to that derived earlier by expansion in harmonics.4 

Next, we generalize the theory to frequency-dependent friction coefficients, as first developed by Stokes for a single sphere on the basis of the linearized Navier–Stokes equations.5 It was realized by Basset6 that this corresponds to retarded friction, as described by a memory term in the equation of motion for the sphere.6,7 We use also retarded hydrodynamic pair interactions8 in the impedance matrix.

The second contribution to the swimming velocity results from flow generated by the swimmer via the movement of the spheres and via Reynolds stress. Since we need only the time-averaged swimming velocity, it suffices to solve the steady-state Stokes equations.9 The generated net flow is called the effect of steady streaming.10 We postpone the study of its contribution to the mean swimming velocity to a separate article.

Here, we consider specifically two different systems, as described by four models of increasing sophistication. In the first system, the two spheres are neutrally buoyant. This is indicated by N. In the second system, the larger sphere is again neutrally buoyant, but the smaller sphere is loaded such that it has twice the mass of the larger sphere. This system is indicated by L. In the first model, the hydrodynamic interactions between spheres are treated in Oseen-approximation. This is indicated by O and is the model studied by Hubert et al.3 In the second model, we take approximate account of added mass effects by including the single sphere added mass in the equations of motion. The model is called the Oseen*-model and is indicated by S. The third model takes account of pair added mass effects in dipole approximation.11 This model is indicated by D. In the fourth model, the friction coefficient of a single sphere is replaced by the frequency-dependent coefficient introduced by Stokes,5 and the pair hydrodynamic interactions are also made frequency-dependent, corresponding to retarded hydrodynamic interactions.8 In addition, the single sphere friction coefficient is corrected for a pair effect such that the added masses are identical with those of the dipole model. This is called the retarded friction model and is indicated by R. Thus, NR refers to the neutrally buoyant system described by the retarded friction model.

We consider two spheres of radii a and b and mass densities ρa and ρb 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 are located on the z axis of a fixed Cartesian system of coordinates. The dynamics of the system is governed by a harmonic interaction potential V int = 1 2 k ( z 2 z 1 d ) 2, depending on the instantaneous configuration of centers z 1 ( t ) , z 2 ( t ), and by actuating forces E 1 ( t ) and E 2 ( t ) in the z direction and varying harmonically in time with frequency ω. We assume in this section that friction with the fluid can be described by a 2 × 2 friction matrix with elements which depend on distance z ( t ) = z 2 ( t ) z 1 ( t ) but are independent of frequency ω. Added mass effects are described by a 2 × 2 mass matrix m ( z ), calculated from potential theory.

In the course of time, the motion of the system becomes periodic with period T = 2 π / ω and with a steady mean swimming velocity U ¯ s w, and the average being over a period. In the following, we derive a simple expression for U ¯ s w for the case of small amplitude harmonic variation of the activating forces E 1 ( t ) = φ 1 cos ( ω t ) and E 2 ( t ) = E 1 ( t ). The mean velocity will be calculated to the second order in the amplitude φ 1.

We summarize the positions of centers in the 2-dimensional vector R = ( z 1 , z 2 ) and the sphere momenta in p = ( p 1 , p 2 ). The momenta are related to the velocities U = ( U 1 , U 2 ) by p = m · U with a symmetric 2 × 2 mass matrix m, which depends on distance z. The dynamics of the system is assumed to be governed by the approximate equations of motion,12 
d R d t = U , d p d t = K R ζ · U V int R + E ,
(2.1)
where the kinetic energy K is given by
K = 1 2 p · w · p ,
(2.2)
with inverse mass matrix w = m 1. The friction matrix ζ and the interaction potential Vint are invariant under translations.

We abbreviated E = ( E 1 , E 2 ). We assume E 2 ( t ) = E 1 ( t ), thus ensuring that the total applied force vanishes at any time. Another choice would imply that the center of mass of the swimmer is swung back and forth by an applied force during its motion. The activating forces cause internal motions of the body without acceleration of its center of mass.

The total sphere momentum P = p 1 + p 2 can be expressed as P = u · p with vector u = ( 1 , 1 ). The momentum varies in time due to friction with the fluid as
d P d t = u · ζ · U .
(2.3)
In the asymptotic periodic swimming situation, we, therefore, have
u · ζ · U ¯ = 0 ,
(2.4)
where the overline bar indicates the average over a period. This shows that in periodic swimming, the mean drag vanishes.
We perform a formal expansion in powers of the amplitude of the activating forces. To second order in the amplitude Eq. (2.4) reads
u · ζ ( 1 ) · U ( 1 ) ¯ + u · ζ 0 · U ( 2 ) ¯ = 0 ,
(2.5)
where ζ 0 is the friction matrix of the rest configuration, which is time-independent. The sphere velocities U = ( U 1 , U 2 ) can be expressed as
U = U u + d ̇ ,
(2.6)
where U = ( U 1 + U 2 ) / 2 is the center velocity and d = ( d 1 , d 2 ) are the deviations from the co-moving rest positions. In the periodic swimming situation, we have
U ¯ = U ¯ u .
(2.7)
From Eq. (2.5), we find for the mechanical mean second order swimming velocity
U ( 2 ) ¯ M = 1 Z 0 u · ζ ( 1 ) · U ( 1 ) ¯ ,
(2.8)
where
Z 0 = u · ζ 0 · u
(2.9)
is the friction coefficient of the rest configuration.
The average on the right-hand side of Eq. (2.8) can be evaluated from linear response theory. The first order friction matrix can be expressed as
ζ ( 1 ) = d ( 1 ) · ζ | 0 ,
(2.10)
where = ( / z 1 , / z 2 ). With activating forces E j ( t ) = Re [ E j ω exp ( i ω t ) ] , ( j = 1 , 2 ), the displacements in linear response take the following form:
d j ( j ) ( t ) = d j c cos ( ω t ) + d j s sin ( ω t ) , j = ( 1 , 2 ) ,
(2.11)
where d j c = Re [ d j ω ] , d j s = Im [ d j ω ]. The velocities are
U j ( 1 ) ( t ) = U j c cos ( ω t ) + U j s sin ( ω t ) , j = ( 1 , 2 ) .
(2.12)
The bilinear averages over a period are
d j ( 1 ) U k ( 1 ) ¯ = 1 2 ( d j c U k c + d j s U k s ) .
(2.13)

We evaluate these quantities from the linear response.

From the equations of motion Eq. (2.1), we find for the one-dimensional situation with coordinates q 1 = z 1 , q 2 = z 2 d and velocities U 1 , U 2 after linearization
d q 1 d t = U 1 , d q 2 d t = U 2 , m 11 d U 1 d t + m 12 d U 2 d t + ζ 11 U 1 + ζ 12 U 2 + k ( q 1 q 2 ) = E 1 ( t ) , m 12 d U 1 d t + m 22 d U 2 d t + ζ 12 U 1 + ζ 22 U 2 k ( q 1 q 2 ) = E 2 ( t ) ,
(3.1)
with constant coefficients taken at the equilibrium distance z = d. For oscillating activating forces at frequency ω, it is convenient to use complex notation. The equation for the amplitudes ψ ω = ( q 1 ω , q 2 ω , U 1 ω , U 2 ω ) then takes the form
M ( ω ) · ψ ω = φ ω ,
(3.2)
with matrix
M ( ω ) = ( i ω 0 1 0 0 i ω 0 1 k k i ω m 11 + ζ 11 i ω m 12 + ζ 12 k k i ω m 12 + ζ 12 i ω m 22 + ζ 22 ) ,
(3.3)
variable vector
ψ ω = ( q 1 ω , q 2 ω , U 1 ω , U 2 ω ) ,
(3.4)
and force vector
φ ω = ( 0 , 0 , E 1 ω , E 1 ω ) .
(3.5)
The solution of Eq. (3.2) is given by
ψ ω = G ( ω ) · φ ω , G ( ω ) = M ( ω ) 1 .
(3.6)
It is convenient to introduce the 2 × 2 impedance matrix
z ( ω ) = ζ i ω m .
(3.7)
From the solution vector ψ ω, we can construct the corresponding displacements
d j ω = q j ω + U ω i ω , ( j = 1 , 2 ) ,
(3.8)
where U ω = ( U 1 ω + U 2 ω ) / 2. It follows that
2 i ω d 1 ω = U 2 ω U 1 ω , 2 i ω d 2 ω = U 1 ω U 2 ω .
(3.9)
Substituting Eq. (2.10) into Eq. (2.8), we see that the mechanical mean swimming velocity is linear in the four correlations d j ( 1 ) U k ( 1 ) ¯. From the solution Eq. (3.6), we find the equalities
d 1 ( 1 ) U 1 ( 1 ) ¯ = d 1 ( 1 ) U 2 ( 1 ) ¯ = d 2 ( 1 ) U 1 ( 1 ) ¯ = d 2 ( 1 ) U 2 ( 1 ) ¯ .
(3.10)
Using this, we obtain
U ( 2 ) ¯ M = 2 Z 0 d Z 0 d 1 ( 1 ) U 1 ( 1 ) ¯ ,
(3.11)
with derivative with respect to d,
Z 0 d = Z 0 / d ,
(3.12)
taken at z = d. In complex notation, we use that for two real quantities A(t) and B(t) oscillating at frequency ω,
A ( t ) = Re [ A ω e i ω t ] , B ( t ) = Re [ B ω e i ω t ] ,
(3.13)
the time-average of the product is given by
A ( t ) B ( t ) ¯ = 1 2 Re [ A ω * B ω ] .
(3.14)

For the Oseen* and Oseen-Dipole models, Eq. (3.11) agrees with the expression for the mean swimming velocity derived by a different method.4 The mean swimming velocity must be supplemented with the additional contribution from the motion of the fluid. First, we show how the above derivation must be modified when the friction coefficients depend on frequency and find a general expression for the moment in Eq. (3.11).

It was first shown by Stokes5 on the basis of the linearized Navier–Stokes equations that the flow pattern about an oscillating sphere differs from that about a steadily moving sphere. Correspondingly, the friction coefficient of a single sphere depends on frequency. For a sphere of radius a with no-slip boundary condition, it takes the value ζ 0 = 6 π η a at zero frequency, but at higher frequency, the complex coefficient ζ ( ω ) differs. It was shown by Basset6 that the frequency-dependent friction corresponds to a retardation effect in the hydrodynamic force exerted by the fluid. If the fluid is viscoelastic, then that is an additional reason for retardation. In the following, we discuss how the retardation effect can be taken into account in the theory of swimming.

We modify the equations of motion Eq. (2.1) to include retarded friction and replace the frictional term ζ · U by Basset-type forces B ( t ) with
B i ( t ) = t β i j ( z 2 ( t ) z 1 ( t ) , t t ) U j ( t ) d t , ( i , j ) = ( 1 , 2 ) ,
(4.1)
with memory kernel β i j ( z , τ ). The momentum balance Eq. (2.3) is modified correspondingly to
d P d t = u · 0 β ( z ( t τ ) , τ ) · U ( t τ ) d τ .
(4.2)
A δ ( t τ )-contribution to the kernel takes account of the added mass effect.7 
We consider the asymptotic periodic swimming motion. Integrating Eq. (4.2) over a period T = 2 π / ω, we obtain
0 T d t u · 0 β ( z ( t τ ) , τ ) · U ( t τ ) d τ = 0.
(4.3)
To lowest order, the swimmer is at rest with positions z 0 and vanishing velocities U ( 0 ) = 0. To first order, the motion is oscillatory, and Eq. (4.3) is satisfied because
0 T d t U ( 1 ) ( t τ ) = 0.
(4.4)
To second order Eq. (4.3) yields an equation for the mechanical mean swimming velocity, in generalization of Eq. (2.8). The equation reads
U ( 2 ) ¯ M = 1 Z 0 0 T d t u · 0 β ( 1 ) ( z 0 , τ ) · U ( 1 ) ( t τ ) d τ .
(4.5)
This shows that the mechanical mean second order swimming velocity can be calculated from the linear response. The equations are the same as before, except that the friction matrix has become frequency-dependent.
In the linear theory, the equations take the form (3.2) with frequency-dependent complex friction coefficients ζ i j ( ω ) and masses mij. The added masses are defined by
m i j = m i δ i j + m aij ,
(4.6)
where m1 and m2 are the bare masses. These definitions imply that the friction coefficients do not have an imaginary part proportional to ω at high frequency. Since only the impedance matrix z ( ω ) is relevant, we could alternatively use the bare masses only and include the added masses in the friction coefficients.
The property
z i j ( ω ) = z i j ( ω ) * , ( i , j ) = ( 1 , 2 )
(4.7)
leads to the same mathematical structure as in Sec. III with a 4 × 4 matrix M ( ω ) with the property (3.7). As a consequence, for real activating forces E 1 ( t ) = E 1 c cos ( ω t )   and   E 2 ( t ) = E 1 ( t ) , the solution ψ ( t ) takes the form
ψ ( t ) = ψ ω c cos ( ω t ) + ψ ω s sin ( ω t ) ,
(4.8)
with real vectors ψ ω c = Re [ ψ ω ] and ψ ω s = Im [ ψ ω ] which follow from the solution Eq. (3.6) with the new matrix M ( ω ) and with φ ω determined from E 1 , E 2 = E 1. We put E 1 c = φ k a, where φ is a dimensionless amplitude.
The complex velocity amplitudes are found from the linear equations as
U 1 ω = i ω φ ka Δ ( k , ω ) ( z 12 + z 22 ) , U 2 ω = i ω φ ka Δ ( k , ω ) ( z 12 z 11 ) ,
(4.9)
with denominator
Δ ( k , ω ) = k Z ( ω ) i ω | z | ,
(4.10)
where | z | is the determinant of the complex impedance matrix z. The corresponding displacement amplitudes are
d 1 ω = d 2 ω = φ kaZ ( ω ) 2 Δ ( k , ω ) .
(4.11)
The expression for the moment in Eq. (3.11) becomes
d 1 ( 1 ) U 1 ( 1 ) ¯ = 1 4 φ 2 k 2 a 2 ω N ( ω ) D ( k , ω ) ,
(4.12)
with numerator
N ( ω ) = z 11 ( z 12 + z 22 ) + z 12 ( z 11 z 22 ) + z 22 ( z 12 + z 11 ) ,
(4.13)
and denominator
D ( k , ω ) = | Δ ( k , ω ) | 2 .
(4.14)
The numerator N ( ω ) can be expressed as
N ( ω ) = 1 2 u · z * ( ω ) · σ y · z ( ω ) · u ,
(4.15)
where σ y is the Pauli spin-matrix,12,
σ y = ( 0 i i 0 ) .
(4.16)
The sign of N ( ω ) determines the direction of swimming. The absolute value | N ( ω ) | equals twice the area of the triangle in the z z plane made up by the points P 0 = ( z 12 , z 12 ) , P 1 = ( z 11 , z 11 ) , and   P 2 = ( z 22 , z 22 ), as given by the so-called shoelace formula.13,14 We call the triangle as the impedance triangle. With 3-vectors P 0 = ( P 0 , 0 ) , P 1 = ( P 1 , 0 ) ,   and   P 2 = ( P 2 , 0 ), and difference vectors S 1 = P 1 P 0 and S 2 = P 2 P 0, the value N ( ω ) is the third component of the cross product N ( ω ) = S 1 × S 2. The vector N ( ω ) can be regarded as a generalized torque corresponding to the impedance matrix z ( ω ).
We wish to choose the amplitude factor φ in such a way that the mean square displacement does not vary with the scale number s. It follows from Eq. (3.9) that d 2 ( 1 ) = d 1 ( 1 ), so that it suffices to consider ( d 1 ( 1 ) ) 2 ¯. For E 1 ( t ) = φ k a cos ( ω t ), the mean square displacement is calculated as
( d 1 ( 1 ) ) 2 ¯ = 1 8 φ 2 k 2 a 2 | Z ( ω ) | 2 D ( k , ω ) .
(4.17)
If we choose
φ 2 = ε 2 D ( k , ω ) k 2 | Z ( ω ) | 2 ,
(4.18)
then ( d 1 ( 1 ) ) 2 ¯ takes the constant value ε 2 a 2 / 8, and the moment in Eq. (4.12) becomes
d 1 ( 1 ) U 1 ( 1 ) ¯ = 1 4 ε 2 a 2 ω N ( ω ) | Z ( ω ) | 2 .
(4.19)
In combination with Eq. (3.11), we, therefore, find for the mechanical contribution to the mean swimming velocity
U ( 2 ) ¯ M = ε 2 a 2 ω Z 0 d 2 Z 0 N ( ω ) | Z ( ω ) | 2 .
(4.20)
We note that this is independent of the elasticity coefficient k. The dependence on k has been fully absorbed in the expression (4.18) for the amplitude factor φ.
In the same way, we calculate the averages ( U 1 ( 1 ) ) 2 ¯ , ( U 2 ( 1 ) ) 2 ¯, and U 1 ( 1 ) U 2 ( 1 ) ¯. We find
( U j ( 1 ) ) 2 ¯ = 1 2 ε 2 a 2 ω 2 N j ( ω ) | Z ( ω ) | 2 , U 1 ( 1 ) U 2 ( 1 ) ¯ = 1 2 ε 2 a 2 ω 2 N 3 ( ω ) | Z ( ω ) | 2 ,
(4.21)
with numerators
N 1 ( ω ) = ( z 12 + z 22 ) 2 + ( z 12 + z 22 ) 2 , N 2 ( ω ) = ( z 12 + z 11 ) 2 + ( z 12 + z 11 ) 2 , N 3 ( ω ) = ( z 12 + z 11 ) ( z 12 + z 22 ) + ( z 12 + z 11 ) ( z 12 + z 22 ) .
(4.22)
Again the averages are independent of the elasticity coefficient.
In this section, we study the frequency-dependence of the mean swimming velocity in more detail. We use the dimensionless number R = a 2 ω ρ / η = 2 s 2 as a measure of frequency. We have R = ω τ v, where τ v = a 2 ρ / η is the viscous relaxation time. In the notation of Derr et al.16, R = M 2. We define the dimensionless effectiveness V(R) as
U s w ¯ = ε 2 a ω V ( R ) .
(5.1)
In words, ε 2 V ( R ) = U s w ¯ / ( a ω ) is the ratio of the distance traveled during a period T = 2 π / ω to the length 2 π a, when swimming with a stroke of amplitude ε at a scale number R. Clearly, one wishes to optimize this value. We recall that the efficiency, as defined by Shapere and Wilczek,16 involves both the mean swimming velocity and the power dissipated during a period. The effectiveness of the stroke is defined regardless of the power.
The mechanical contribution to the effectiveness V M ( R ) in the limit ε 0 is found by combination of Eqs. (4.20) and (5.1) as
V M ( R ) = a Z 0 d 2 Z 0 N ( ω ) | Z ( ω ) | 2 .
(5.2)
This takes a simple form for models with frequency-independent friction matrix. In that case, we find from Eq. (4.15)
N ( ω ) = N ω ,
(5.3)
with coefficient
N = i 2 u · ( m · σ y · ζ ζ · σ y · m ) · u ,
(5.4)
given by
N = m 11 ( ζ 12 + ζ 22 ) m 12 ( ζ 11 ζ 22 ) m 22 ( ζ 12 + ζ 11 ) = 1 2 ( Z 0 Δ m M 0 Δ ζ ) ,
(5.5)
with differences
Δ m = m 11 m 22 , Δ ζ = ζ 11 ζ 22 .
(5.6)

This shows an interesting symmetry between mass and friction.

We find that V M ( R ) takes the form4 
V M 2 ( R ) = V x M 2 R x M R R x M 2 + R 2 ,
(5.7)
with parameters RxM and VxM. The function V M 2 ( R ) has maximum value VxM at RxM, and the curve has half-width Δ R = 2 3 R x M. We call a curve in the VR-plane with behavior described by Eq. (5.7) a VR2-curve.
From the behavior at small and large R
V M 2 ( R ) C R as    R 0 , V M 2 ( R ) L / R as    R ,
(5.8)
we find
R x M = L / C , V x M = 1 2 L C .
(5.9)
From Eq. (4.20), we find the values
C = a Z 0 d 2 Z 0 η | N | a 2 ρ Z 0 2 , L = a Z 0 d 2 Z 0 a 2 ρ | N | η M 2 ,
(5.10)
where M = u · m · u = m 11 + 2 m 12 + m 22 is the total mass. Therefore,
R x M = a 2 ρ Z 0 η M , V x M = a Z 0 d 4 Z 0 | N | M Z 0 .
(5.11)
The aforementioned expressions can be cast in the form4 
V M 2 ( R ) = Im χ 0 1 i ω τ M ,
(5.12)
with
τ M = τ v R x M = M Z 0 , χ 0 = 2 V x M R x M .
(5.13)
This shows that for such models, the dependence on frequency is quite simple.
In earlier work,4 we derived explicit expressions for the parameters RxM and VxM for three specific models. The Oseen-model is based on the Oseen-approximation for the hydrodynamic pair mobility, and added mass effects are neglected. The Oseen*-model uses the same approximation for the pair mobility and includes the single sphere added mass effect. The Oseen-Dipole model describes in addition the pair added mass in dipole approximation. We note that for the Oseen*-model M = m 1 * + m 2 *, where m j * = m j + 1 2 m f j with m f j = ( 4 π / 3 ) ρ a j 3. For the Oseen-Dipole model, M = m 11 + 2 m 12 + m 22 with mij given below. The factors Z 0 d and Z0 are the same for the three models. The Oseen mobility matrix corresponds to the friction matrix,
ζ = 12 π η z 4 z 2 9 a b ( 2 a z 3 a b 3 a b 2 b z ) ,
(5.14)
where z = z 2 z 1. In the Oseen-Dipole model, one uses a dipole approximation to the mass matrix given by Felderhof11,
m = 2 π 3 z 6 3 a 3 b 3 ( a 3 [ z 6 ( ρ + 2 ρ a ) + 2 a 3 b 3 ( ρ ρ a ) ] 3 a 3 b 3 z 3 ρ 3 a 3 b 3 z 3 ρ b 3 [ z 6 ( ρ + 2 ρ b ) + 2 a 3 b 3 ( ρ ρ b ) ] ) .
(5.15)
In an expansion in inverse powers of z to terms of order z 6, this agrees with the expression given by Lamb.21 The diagonal elements include the single sphere added mass effect.
The mechanical contribution to the effectiveness, given by Eq. (5.2), can be evaluated straightforwardly for the three models by use of the above expressions for ζ and m. We find agreement with those found before,4 except for a copying error in the coefficient rb in Eq. (6.16) of Ref. 4. This should read instead
r b = 2 ξ 2 ( ξ δ 2 ) ( 2 δ 3 ξ ) ( δ 4 + δ 2 ξ + ξ 2 ) ,
(5.16)
where ξ = b / a and δ = d / a.
In particular, we find for the Oseen-model
V M O ( R ) = 27 4 s 2 ( 2 δ 3 ) ( 2 δ 3 ξ ) δ 3 ξ + δ ξ ( 2 δ 3 ) σ a ( 2 δ 3 ξ ) ξ 2 σ b 81 δ 2 ( δ 3 ξ + δ ξ ) 2 + ( 4 δ 2 9 ξ ) 2 ( σ a + ξ 3 σ b ) 2 s 4 ,
(5.17)
where σ a = ρ a / ρ and σ b = ρ b / ρ. We compare with the expression for the effectiveness in the Oseen-model, as calculated by Hubert et al.3 for large d/a. Their Eq. (18) reads in our notation
V M ¯ S = 3 ξ 2 2 ( 1 + ξ ) 3 δ 2 ω ( τ a τ b ) 1 + ω 2 τ M O 2 ,
(5.18)
where τ a = m 1 / ( 6 π η a ) , τ b = m 2 / ( 6 π η b ) ,   and   τ M O = ( m 1 + m 2 ) / ( 6 π η ( a + b ) ). This agrees to order δ 2 with Eq. (5.17). For the Oseen*-model, the same equations hold with mi replaced by m i * = m i + ( 2 π / 3 ) ρ a i 3 or with σi replaced by σ i * = σ i + 1 2.

The numerator of the last factor in Eq. (5.17) for fixed σ a , ξ , δ can switch sign as a function of σb. For example, for σ a = 1 , ξ = 1 / 2 , and   δ = 3, the expression switches sign at σ b = 8 / 3. This corresponds to a change of swimming direction.

For the Oseen-model, the maximum value VxMO depends on the mass densities ρ a , ρ b only via the ratio r = ρ b / ρ a = σ b / σ a. The expression is, from Eq. (5.11),
V xMO = 3 ξ 2 ( 2 δ 3 ) ( 2 δ 3 ξ ) ( 3 2 δ + 2 r δ ξ 2 3 r ξ 3 ) 8 δ ( 4 δ 9 ξ 2 ) ( δ 3 ξ + δ ξ ) 2 ( 1 + q ξ 3 ) .
(5.19)
To get the value for the Oseen*-model, replace r by r * = σ b * / σ a *. For the position of the maximum, one finds
R xMO = 18 δ ( δ 3 ξ + δ ξ ) ( 4 δ 2 9 ξ ) ( 1 + r ξ 3 ) 1 σ a ,
(5.20)
and similarly for RxMS. These simple properties do not hold for the Oseen-Dipole model. For typical values of ξ , δ, the value of VxMO ranges from positive to negative values. From Eq. (5.19), one sees that in the Oseen-model, the sign switch occurs at
r 0 = 2 δ 3 ξ 2 ( 2 δ 3 ξ ) .
(5.21)
As a standard system, we choose a two-sphere with parameters b = 1 2 a , d = 3 a in a fluid with viscosity η and mass density ρ. The mass densities of the spheres can be varied and take values ρ a = σ a ρ and ρ b = σ b ρ. We consider first neutrally buoyant spheres with σ a = 1 and σ b = 1. In Fig. 2, we show the effectiveness VM as a function of R for the Oseen-model, the Oseen*-model, and the Oseen-Dipole model. In the notation introduced at the end of Sec. I, the three models are indicated, respectively, as the NO-model, the NS-model, and the ND-model. In Table I, we list the values of the parameters RxM and VxM for the three models.
FIG. 2.

Plot of the mechanical effectiveness VM for a two-sphere with a = 1, b = 1/2, d = 3, and ρ a = ρ b = ρ, with hydrodynamic interactions of Oseen-tye (NO: long dashes), dipole-type added mass (ND: short dashes), and retarded friction (NR: solid curve). The curve for NS is not shown, since it cannot be distinguished from the one for ND.

FIG. 2.

Plot of the mechanical effectiveness VM for a two-sphere with a = 1, b = 1/2, d = 3, and ρ a = ρ b = ρ, with hydrodynamic interactions of Oseen-tye (NO: long dashes), dipole-type added mass (ND: short dashes), and retarded friction (NR: solid curve). The curve for NS is not shown, since it cannot be distinguished from the one for ND.

Close modal
TABLE I.

List of values characterizing the mechanical swimming speed of a two-sphere with radii a = 1 , b = 1 / 2, equilibrium distance between centers d = 3, in a fluid with viscosity η = 1, mass density ρ = 1, for four models: Oseen (O), Oseen*(S), dipole (D), retarded friction (R), and of two systems N = ( ρ , ρ ) or L = ( ρ , 16 ρ ). The listed values are explained in Sec. V. In the case of retarded friction, the relations in Eq. (5.11) do not hold.

( ρ a , ρ b ) Model R 2 x M 10 6 V 2 x M Z 0 / ( η a ) M / ( ρ a 3 )
( ρ , ρ )  Oseen  4.571  2480  21.5423  4.712 
( ρ , ρ )  Osee n *  3.048  2480  21.5423  7.069 
( ρ , ρ )  Dipole  3.072  2538  21.5423  7.012 
( ρ , ρ )  RetFr  5.255  1298  21.7598  7.014 
( ρ , 16 ρ )  Oseen  1.714  –7440  21.5423  12.566 
( ρ , 16 ρ )  Osee n *  1.444  –5874  21.5423  14.923 
( ρ , 16 ρ )  Dipole  1.449  –5879  21.5423  14.866 
( ρ , 16 ρ )  RetFr  4.409  –2827  21.7598  14.868 
( ρ a , ρ b ) Model R 2 x M 10 6 V 2 x M Z 0 / ( η a ) M / ( ρ a 3 )
( ρ , ρ )  Oseen  4.571  2480  21.5423  4.712 
( ρ , ρ )  Osee n *  3.048  2480  21.5423  7.069 
( ρ , ρ )  Dipole  3.072  2538  21.5423  7.012 
( ρ , ρ )  RetFr  5.255  1298  21.7598  7.014 
( ρ , 16 ρ )  Oseen  1.714  –7440  21.5423  12.566 
( ρ , 16 ρ )  Osee n *  1.444  –5874  21.5423  14.923 
( ρ , 16 ρ )  Dipole  1.449  –5879  21.5423  14.866 
( ρ , 16 ρ )  RetFr  4.409  –2827  21.7598  14.868 

For our choice of parameters ξ = 1 / 2 , δ = 3, we have R xMNO = 32 / 7 = 4.571 and V xMNO = 5 / 2016 = 0.002 480. As a function of log 10 R, the curve resembles a Gaussian. For the NS-model, the corresponding numbers read: R xMNS = 64 / 21 = 3.048 and V xMNS = 5 / 2016 = 0.002 480. We compare with the curve V MND ( R ) for the ND-model, characterized by the same friction matrix Eq. (5.14), but with mass matrix given by Eq. (5.15). The function V MND ( R ) has the form Eq. (5.7) with parameter values RxMND = 3.063 and VxMND = 0.00 2538. Thus, the two V R 2 curves are nearly identical.

For the NS-model, the points Pj of the impedance triangle, defined below Eq. (4.16), have coordinates P 0 = ( 5.386 η a , 0 ) , P 1 = ( 21.542 η a , 19.149 ω ρ a 3 ) , and   P 2 = ( 10.771 η a , 2.394 ω ρ a 3 ) at R = R xMO. The value of the coefficient N is N = 21.15 η ρ a 4. For the ND-model, the points are P 0 = ( 5.386 η a , 0.089 ω ρ a 3 ) , P 1 = ( 21.542 η a , 19.308 ω ρ a 3 ) ,   and   P 2 = ( 10.771 η a , 2.413 ω ρ a 3 ), and N = 21.47 η ρ a 4.

As a second example, we consider the L-system with the same geometry, but with mass densities ρ a = ρ and ρ b = 16 ρ, so that m 2 = 2 m 1. For this geometry r 0 = 8 / 3, so that the r-values of the two systems, r = 1 and r = 16, lie on opposite sides of r0. Hence, the two systems swim in opposite directions. We list the corresponding values in Table I.

In this section, we study the effect of frequency-dependent friction. We do not aim at a complete calculation based on exact frequency-dependent friction coefficients but limit ourselves to the dominant effects for not too small distance between centers. At very large distance, the interaction tends to zero, and the friction matrix is diagonal with single sphere terms given by the frequency-dependent friction coefficient first calculated by Stokes.5 

The mutual mobility will be approximated by a result,8 based on a single Green function between the spheres and on Faxén theorems. We have used the approximation earlier in a study of velocity relaxation of a pair of spheres.17 In the limit of zero frequency, the approximation reduces to that of Rotne and Prager,18 as generalized to spheres of different sizes.19 The Rotne–Prager approximation is popular in soft matter modeling.20 The self-mobility is given by the single sphere friction coefficient, corrected by a distance-dependent term, as detailed below, chosen such that the mass matrix is identical with that of the Oseen-Dipole model, given by Eq. (5.15).

The single sphere friction coefficient of spheres A and B is, respectively,
ζ A ( ω ) = 6 π η a A 0 ( α a ) , ζ B ( ω ) = 6 π η b A 0 ( α b ) ,
(6.1)
with the abbreviation7,
A 0 ( λ ) = 1 + λ + 1 9 λ 2 ,
(6.2)
and α = i ω ρ / η , Re [ α ] > 0. The approximate mutual mobility function for two spheres with distance z between centers reads8 
α A B t t ( z , ω ) = B 0 ( α a ) B 0 ( α b ) ( 1 + α z ) e α ( a + b z ) 2 π η α 2 A 0 ( α a ) A 0 ( α b ) z 3 ,
(6.3)
with
B 0 ( λ ) = 1 + λ + 1 3 λ 2 .
(6.4)
The pair mobility matrix will be approximated by
μ ( z , ω ) = ( α A A t t ( z , ω ) α A B t t ( z , ω ) α A B t t ( z , ω ) α B B t t ( z , ω ) ) ,
(6.5)
with the off diagonal terms given by Eq. (6.3). The diagonal terms are given by the single sphere terms with a Van der Waals-type pair correction, corresponding to a dipolar reflection of a dipole field,
    α A A t t ( z , ω ) = 1 6 π η a A 0 ( z , ω ) , A 0 ( z , ω ) = 1 + α a + 1 9 ( 1 6 a 3 b 3 z 6 ) α 2 a 2 .
(6.6)
The approximate pair friction matrix is
ζ ( z , ω ) = μ ( z , ω ) 1 .
(6.7)
The real and imaginary parts of the matrix elements are denoted as ζ i j and ζ i j , respectively.
The pair friction matrix given by Eqs. (6.1)–(6.7) increases with frequency in proportion to i ω. It, therefore, contains added mass terms, as explained in Eq. (4.6). We incorporate all added mass effects in the friction matrix and use the impedance matrix
z ( ω ) = ζ i ω m 0 ,
(6.8)
rather than Eq. (3.7). With ζ ( z , ω ) given by Eqs. (6.1)–(6.7), the added mass matrix is identical with that given by Eq. (5.15).
In the expression (5.2) for the function V M ( R ), we need the pre-factor a Z 0 d / Z 0, where Z 0 = u · ζ ( 0 ) · u is the steady-state friction coefficient and Z 0 d its derivative with respect to d. With the above expressions, we find
Z 0 = 24 π η d 3 ( a + b ) d 3 3 a b d 2 + a b ( a 2 + b 2 ) 4 d 6 9 a b d 4 + 6 a b ( a 2 + b 2 ) d 2 a b ( a 4 + b 4 ) 2 a 3 b 3 .
(6.9)
This is identical with the Rotne-Prager-approximation. For large d,
Z 0 d Z 0 = 3 a b ( a + b ) d 2 9 a b ( a 2 + b 2 ) 2 ( a + b ) 2 d 3 + 3 a b ( a + b ) 3 d 4 ( a 4 + 2 a 3 b 7 a 2 b 2 + 2 a b 3 + b 4 ) + O ( d 5 ) .
(6.10)

To the order shown, this is identical with the exact result.19 

In Fig. 3, we show the mechanical effectiveness V M ( R ) for the RF-model as specified above. For the N-system, the maximum is at R mMNR = 5.573 and has value VmMNR = 0.001 300. We compare with the VR2-curve for these parameter values. It is evident that the retardation of friction has a significant effect. In a plot of V vs R, the half-width of the actual curve is Δ R = 73.28, whereas for the V R 2 curve, Δ R 2 = 2 3 R mMNR = 19.30. The values for the minimum for the L-system with m 2 = 2 m 1 are RmMLR = 4.389 and V mMLR = 0.002 828. The half-width of the spectral curve is Δ R = 64.02, whereas for the V R 2 curve, Δ R 2 = 2 3 R mMLR = 15.20. Hence, in both cases, the retarded friction causes a shift to higher frequency of the extremum and a strong broadening of the spectral curve. In Fig. 4, we show the analogue of Fig. 2 for the L-system, and in Fig. 5, we show the analogue of Fig. 3.

FIG. 3.

Plot of mechanical effectiveness VM for the NR-model (solid curve) compared with curve of the form (5.7) for the same maximum (dashed).

FIG. 3.

Plot of mechanical effectiveness VM for the NR-model (solid curve) compared with curve of the form (5.7) for the same maximum (dashed).

Close modal
FIG. 4.

Plot of the mechanical effectiveness VM for a two-sphere with a = 1, b = 1/2, d = 3, and ρ a = ρ , ρ b = 16 ρ, with hydrodynamic interactions of Oseen-tye (LO: long dashes), dipole-type added mass (LD: short dashes), and retarded friction (LR: solid curve). The curve for LS is not shown, since it cannot be distinguished from the one for LD.

FIG. 4.

Plot of the mechanical effectiveness VM for a two-sphere with a = 1, b = 1/2, d = 3, and ρ a = ρ , ρ b = 16 ρ, with hydrodynamic interactions of Oseen-tye (LO: long dashes), dipole-type added mass (LD: short dashes), and retarded friction (LR: solid curve). The curve for LS is not shown, since it cannot be distinguished from the one for LD.

Close modal
FIG. 5.

Plot of the mechanical effectiveness VM for the LR-model (solid curve) compared with curve of the form (5.7) for the same minimum (dashed).

FIG. 5.

Plot of the mechanical effectiveness VM for the LR-model (solid curve) compared with curve of the form (5.7) for the same minimum (dashed).

Close modal

In the earlier, we analyzed the mechanical contribution to the mean swimming velocity of a two-sphere swimmer. We showed that this contribution can be expressed in terms of the impedance matrix of the two-sphere. The latter embodies the linear response of the two-sphere to applied oscillatory forces acting on either sphere. The response depends strongly on frequency. We studied this dependence for various, more or less realistic choices of the hydrodynamic interactions incorporating the effects of friction and added mass.

Part of the desired information on the mean swimming velocity is missing. We must yet add the hydrodynamic contribution due to fluid flow generated by the oscillatory motion of the two spheres and affected by the Reynolds stress in the fluid. Many aspects of the required calculations were discussed by Derr et al.15 In separate work, we shall compare the hydrodynamic and mechanical contributions. Once a reliable calculation of the mean swimming velocity of the vibrating two-sphere is established, we can apply the method to other systems of interest, for example to a vibrating chain of three spheres.

The authors have no conflicts to disclose.

B. U. Felderhof: Conceptualization (lead); Writing – review & editing (lead).

The data that support the findings of this study are available within the article.

1.
D.
Klotsa
,
K. A.
Baldwin
,
R. J. A.
Hill
,
R. M.
Bowley
, and
M. R.
Swift
, “
Propulsion of a two-sphere swimmer
,”
Phys. Rev. Lett.
115
,
248102
(
2015
).
2.
B.
Cichocki
,
B. U.
Felderhof
, and
R.
Schmitz
, “
Hydrodynamic interactions between two spherical particles
,”
PhysicoChem. Hyd.
10
,
383
(
1988
).
3.
M.
Hubert
,
O.
Trosman
,
Y.
Collard
,
A.
Sukhov
,
J.
Harting
,
N.
Vandewalle
, and
A.-S.
Smith
, “
Scallop theorem and swimming at the mesoscale
,”
Phys. Rev. Lett.
126
,
224501
(
2021
).
4.
B. U.
Felderhof
, “
Optimizing the mean swimming velocity of a model two-sphere swimmer
,”
Phys. Fluids
34
,
081901
(
2022
).
5.
G. G.
Stokes
, “
On the effect of the internal friction of fluids on the motion of pendulums
,”
Trans. Cambridge Philos. Soc.
9
,
8
(
1851
).
6.
A. B.
Basset
, “
On the motion of a sphere in a viscous liquid
,”
Philos. Trans. R. Soc., A
179
,
43
(
1888
).
7.
S.
Kim
and
S. J.
Karrila
,
Microhydrodynamics: Principles and Selected Applications
(
Butterworth-Heinemann
,
Boston
,
1991
).
8.
B. U.
Felderhof
, “
Retarded hydrodynamic interaction between two spheres immersed in a viscous incompressible fluid
,”
Phys. Fluids
31
,
053604
(
2018
).
9.
B. U.
Felderhof
and
R. B.
Jones
, “
Inertial effects in small-amplitude swimming of a finite body
,”
Physica A
202
,
94
(
1994
).
10.
N.
Riley
, “
Steady streaming
,”
Annu. Rev. Fluid Mech.
33
,
43
(
2001
).
11.
B. U.
Felderhof
, “
Effect of inertia on laminar swimming and flying of an assembly of rigid spheres in an incompressible viscous fluid
,”
Phys. Rev. E
92
,
053011
(
2015
).
12.
E.
Merzbacher
,
Quantum Mechanics
(
Wiley
,
New York
,
1970
).
13.
C. E.
Garza-Hume
,
M. C.
Jorge
, and
A.
Olivera
, “
Areas and shapes of planar irregular polygons
,”
Forum Geometricorum
18
,
17
(
2018
).
14.
B.
Polster
,
The Shoelace Book
(
American Mathematical Society, Providence
,
2006
).
15.
N. J.
Derr
,
T.
Dombrowski
,
C. H.
Rycroft
, and
D.
Klotsa
, “
Reciprocal swimming at intermediate Reynolds number
,”
J. Fluid Mech.
952
,
A8
(
2022
).
16.
A.
Shapere
and
F.
Wilczek
, “
Efficiencies of self-propulsion at low Reynolds number
,”
J. Fluid Mech.
198
,
587
(
1989
).
17.
B. U.
Felderhof
, “
Collinear velocity relaxation of two spheres in a viscous incompressible fluid
,”
Phys. Rev. E
101
,
043103
(
2020
).
18.
J.
Rotne
and
S.
Prager
, “
Variational treatment of hydrodynamic interaction in polymers
,”
J. Chem. Phys.
50
,
4831
(
1969
).
19.
B. U.
Felderhof
, “
Hydrodynamic interactions between two spheres
,”
Physica A
89
,
373
(
1977
).
20.
P. J.
Zuk
,
E.
Wajnryb
,
K. A.
Mizerski
, and
P.
Szymczak
, “
Rotne–Prager–Yamakawa approximation for different-sized particles in application to macromolecular bead models
,”
J. Fluid Mech.
741
,
R5
(
2014
).
21.
H.
Lamb
,
Hydrodynamics
(
Dover
,
New York
,
1945
).