Effect of inertia on laminar swimming and flying of an assembly of rigid spheres in an incompressible viscous fluid

# Effect of inertia on laminar swimming and flying of an assembly of rigid spheres in an incompressible viscous fluid

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

A mechanical model of swimming and flying in an incompressible viscous fluid in the absence of gravity is studied on the basis of assumed equations of motion. The system is modeled as an assembly of rigid spheres subject to elastic direct interactions and to periodic actuating forces which sum to zero. Hydrodynamic interactions are taken into account in the virtual mass matrix and in the friction matrix of the assembly. An equation of motion is derived for the velocity of the geometric center of the assembly. The mean power is calculated as the mean rate of dissipation. The full range of viscosity is covered, so that the theory can be applied to the flying of birds, as well as to the swimming of fish or bacteria. As an example a system of three equal spheres moving along a common axis is studied.

###### pacs:
47.15.G-, 47.63.mf, 47.63.Gd, 45.50.Jf

## I Introduction

The swimming of fish and the flying of birds continue to pose challenging theoretical problems. The physics of bird flight was first studied in detail by Otto Lilienthal in the nineteenth century 1 (). Since then, significant progress has been made in many years of dedicated research 2 ()-5 ().

The goal of theory is to calculate the time-averaged speed and power for given periodic shape variations of the body, at least for a simple model system. It is assumed that the motion of the fluid is well described by the Navier-Stokes equations for an incompressible viscous fluid. On average over a period the force exerted by the body on the fluid vanishes, so that thrust and drag cancel. In early work by Lighthill 6 () and Wu 7 () the thrust and power were calculated approximately as functions of the speed on the basis of potential flow theory for a planar strip. This work and subsequent developments have been reviewed by Childress 3 (), by Wu 8 (),9 (), and by Sparenberg 10 (). However, an independent calculation of the mean speed for given periodic shape variations is still lacking. Measurement of the power consumption has led to a surprisingly small friction coefficient, much smaller than that of an inert body, as was first observed by Gray 11 ().

It was first shown by Taylor 12 () that in the slow swimming of a microorganism the calculation of thrust can be circumvented. In this limiting case one can use the time-independent Stokes equations. The mean swimming velocity and mean rate of dissipation then follow from a purely kinematic calculation 13 (),14 (). For small amplitude swimming both quantities are quadratic in the amplitude of the stroke to lowest order. For a simple system, where the body is modeled as an assembly of rigid spheres held together by direct interaction forces and subject to periodic actuating forces which sum to zero, we have shown that in the high viscosity limit the swimming velocity and power can be calculated for any amplitude of stroke from kinematics alone 15 (),16 ().

In the following we investigate questions of thrust, velocity, and power for swimming or flying in a fluid of any viscosity, including the limit of low viscosity, for the same mechanical model as before. We assume for simplicity that the spheres experience Stokes friction. In addition we incorporate hydrodynamic interactions via virtual mass effects, as found from potential flow theory. We use Hamilton’s equations of motion with added damping terms. In the limit of high viscosity, where resistive forces dominate, the earlier results are recovered. The model provides valuable insight also in the limit of low viscosity, where reactive forces dominate. In that regime the motion is dominated by virtual mass effects. Bernoulli forces and modified linear friction should be taken into account in a more realistic model. Nonetheless, the principle of the calculation, which exploits elimination of the fluid degrees of freedom, remains valid.

The flow is assumed to be laminar at all times. It is now realized that the boundary layer of swimming fish is laminar even at high Reynolds number 9 (). Virtual mass effects were discussed earlier by Lighthill 17 (). The numerical modeling of animal swimming and flight was reviewed by Deng et al. 18 ().

As an example a system of three equal spheres moving along a common axis is studied. For this simple system the mean swimming speed and mean power to second order in the amplitude of stroke can be evaluated analytically. The solution to a corresponding eigenvalue problem provides the optimal stroke to this order, as we found elsewhere in the resistive regime 15 ().

In our model the mean thrust, i.e. the frictional force exerted on the fluid averaged over a period in periodic swimming, vanishes identically. We find that the velocity of the geometric center of the assembly is driven by a different force, which we call the impetus. It has both a reactive and a resistive component. The impetus determines the center velocity with retardation. The mean impetus does not vanish.

It is known for small amplitude swimming in the resistive regime that the mean power is directly proportional to the mean velocity. We find for our example that the relation between mean power and mean velocity is nearly linear also for large amplitude swimming. Presumably the near linearity holds also for other systems in the whole regime of viscosity. If true, this would resolve the so-called Gray paradox 9 (), which is based on the mistaken notion that the power is quadratic in the velocity, as in Stokes friction.

## Ii Equations of motion

We consider a set of rigid spheres of radii and masses , centered at positions , and immersed in an incompressible viscous fluid of shear viscosity and mass density . The fluid is of infinite extent in all directions. The flow velocity and pressure of the fluid are assumed to satisfy the Navier-Stokes equations

 ρ[∂\boldmathv∂t+\boldmathv% ⋅∇\boldmathv]=η∇2\boldmathv−∇p,∇⋅\boldmathv=0. (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 .

As the spheres move in the fluid they experience a frictional force. In addition there may be applied forces and direct interaction forces which depend on the relative positions of sphere centers. We shall assume that the sum of applied forces vanishes, so that

 N∑j=1\boldmathEj(t)=0. (2)

The sum of direct interaction forces vanishes owing to Newton’s third law. We assume that the frictional forces are linear in the sphere velocities, as given by low Reynolds number hydrodynamics on the slow timescale 19 (). The spheres are freely rotating so that there are no frictional torques.

The forces exerted by pressure gradients resist instantaneous acceleration and give rise to virtual mass effects 20 (). For a single sphere immersed in infinite fluid the virtual mass would be , where is the mass of fluid displaced by the sphere. The virtual mass effect for a collection of spheres is embodied in a -dimensional mass matrix . This can be derived from potential flow theory by considering the irrotational flow pattern generated instantaneously by a set of sudden impulses from a state of rest. The total corresponding kinetic energy is

 K=12U⋅m⋅U, (3)

where is the set of sphere velocities. The kinetic energy is a sum of two contributions, , with

 Kp=12N∑j=1mpjU2j,Kf=12ρ∫Vf(∇ϕ)2d\boldmathr, (4)

where the integration is over the part of space occupied by fluid and is the velocity potential corresponding to the set of instantaneous velocities . The potential is linear in the sphere velocities , so that the kinetic energy is a quadratic form in . The contribution to the total kinetic energy defines the virtual mass. This depends parametrically on the positions , leading to corresponding hydrodynamic interactions.

Besides the mass matrix it will be useful to define also the mobility matrix , the friction matrix , and the inverse mass matrix ,

 \boldmathμ=\boldmathζ−1,w=m−1. (5)

The four matrices are symmetric and depend on the relative positions of the sphere centers. The sphere momenta, including the virtual mass contribution, are given by

 p=m⋅U,U=w⋅p. (6)

Correspondingly we postulate the equations of motion

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

where is given by and is the potential of direct interaction forces. The partial derivative is taken at constant momenta . It is clear that the equations of motion (2.7) have only limited validity. The frictional forces are assumed to be linear in the sphere velocities, and Basset memory forces are neglected. Nonetheless it is of interest to explore the consequences of the equations as they stand.

We note that it follows from Eq. (2.7) that total dressed sphere momentum is not conserved but changes due to friction with the fluid. If an impulse is imparted at time to the spheres in a state of rest due to applied forces , then part of the momentum is transferred instantaneously to the fluid, reducing the sphere velocities to , as can be seen by integration over an infinitesimal time interval about and use of Eq. (2.6). The remaining momentum of the spheres is transferred to the fluid in the course of time by friction. The total momentum of spheres and fluid is conserved at all times.

## Iii Impetus, center velocity, and rate of dissipation

We are looking for a solution of the equations of motion where the center of the assembly moves on average with uniform translational velocity and the individual spheres perform periodic motions about the moving center. The velocity of the geometric center is defined by

 Uα=1NU⋅uα,α=(x,y,z) (8)

where the symbol denotes a -dimensional vector with on the positions, on the positions, and cyclic. The -dimensional displacement vector is defined by

 \boldmathRj(t)=\boldmathRj0+∫t0% \boldmathU(t′)dt′+\boldmathδj(t),j=1,...,N, (9)

with the property

 d(t)⋅uα=0,α=(x,y,z), (10)

and with a set of equilibrium positions for which the direct interaction forces vanish. Correspondingly the velocity vector is decomposed as

 U=Uβuβ+˙d. (11)

Summation over repeated greek indices is implied. By substitution into Eqs. (2.6) and (2.7) we find

 ddt[m⋅(Uβuβ+˙d)]+∂K∂R+\boldmathζ⋅(Uβuβ+˙d)+∂Vint∂R=E(t). (12)

Multiplying from the left with the vector we obtain

 (13)

with mass tensor and friction tensor defined by

 Mαβ=uα⋅m⋅uβ,Zαβ=uα⋅\boldmathζ⋅uβ. (14)

We have used the fact that depends only on relative coordinates, so that . Also we have used Newton’s third law and Eq. (2.2). We note that in Eq. (3.6) the center velocity occurs only in the two places explicitly shown.

We rewrite Eq. (3.6) as

 (15)

with time-dependent impetus given by

 Iα(t)=−ddt(uα⋅m⋅˙d)−uα⋅\boldmathζ⋅˙d. (16)

The mean impetus, averaged over a period ,

 ¯¯¯¯¯¯Iα=1τ∫τ0Iα(t)dt=−uα⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathζ⋅˙d, (17)

does not vanish, even though the total mean force exerted by the fluid on the assembly does vanish. The drag exerted on the assembly by the fluid is

 Dα=−uα⋅\boldmathζ⋅U=−ZαβUβ−uα⋅\boldmathζ⋅˙d. (18)

The thrust is equal and opposite to the drag, . There are frictional forces on the spheres, but on time average the total drag vanishes, , as follows from Eqs. (2.2), (2.7) and Newton’s third law for the interaction forces. We use the fact that the integral of over a period vanishes by periodicity, as well as the property , mentioned below Eq. (3.7). From it follows that the mean impetus can also be expressed as

 ¯¯¯¯¯¯Iα=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ZαβUβ. (19)

We may regard Eq. (3.8) as a balance between central and internal motion. The equation must be solved for the center velocity for given impetus , which can be calculated from the displacement vector . In the resistive limit the total drag vanishes at any time 21 () and the reactive forces vanish, so that then .

The displacement vector may be calculated from displacements in relative space by using a transformation to center and relative coordinates. The geometric center of the assembly is given by

 \boldmathR=1NN∑j=1\boldmathRj=1N\boldmatheαuα⋅R (20)

with Cartesian unit vectors . We define relative coordinates with as

 \boldmathr1 = \boldmathR2−\boldmathR1,% \boldmathr2=\boldmathR3−\boldmathR2,..., \boldmathrN−1 = \boldmathRN−\boldmathRN−1, (21)

and denote the corresponding -vector . The -vector is related to the vector by a transformation matrix according to

 (\boldmathR,r)=T⋅R (22)

with explicit form given by Eqs. (3.13) and (3.14). The displacement vector is derived from the displacement in relative space as

 d=T−1⋅(\boldmath0,\boldmathξ). (23)

Therefore the impetus is determined by displacements in relative space.

The time-dependent rate of dissipation can be expressed in the same matrix formalism. The rate of dissipation is given by

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

Once the center velocity has been calculated for known displacements we can also calculate the time-dependent rate of dissipation by use of Eq. (3.4). However, we derive an alternative expression which will be useful in the following. We can solve for the product from the equation of motion

 ddt(m⋅U)+\boldmathζ⋅U=F, (25)

where is the vector of mechanical forces acting on the spheres,

 F=−∂K∂R−∂Vint∂R+E, (26)

which has the property . Using this property we find for the rate of dissipation

 D=U⋅F−U⋅ddt(m⋅U)=˙d⋅F−U⋅ddt(m⋅U). (27)

Using here Eq. (3.18) again we can rewrite this as

 D=˙d⋅\boldmathζ⋅˙d+Uα˙d⋅fα−Uαuα⋅ddt(m⋅U), (28)

with friction vector . This generalizes an expression derived earlier in the resistive limit 16 ().

## Iv Small amplitude swimming

For vanishing displacements Eq. (3.8) has the solution . By formal series expansion in powers of the displacement vector we obtain a corresponding expansion of the center velocity

 \boldmathU=\boldmathU(1)+\boldmathU(2)+% \boldmathU(3)+.... (29)

The first order velocity satisfies the equation

 M0αβdU(1)βdt+Z0αβU(1)β=−uα⋅m0⋅¨d−uα⋅\boldmathζ0⋅˙d, (30)

where the superscript indicates that the quantity is calculated for the configuration . In particular for oscillating displacement, in complex notation

 d(t)=Re[dωe−iωt], (31)

we find correspondingly

 U(1)αω=[−iω\boldmathM0+% \boldmathZ0]−1αβuβ⋅(ω2m0+iω\boldmathζ0)⋅dω. (32)

In this situation the first order velocity oscillates in time and vanishes on time average.

Multiplying Eq. (3.8) by the mobility tensor and expanding to second order we obtain a more complicated equation for the second order velocity . It suffices to derive an expression for the time-average over a period ,

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathU(2)=1τ∫τ0% \boldmathU(2)(t)dt. (33)

Using periodicity and the first order equation Eq. (4.2) we obtain the expression

 ¯¯¯¯¯¯¯¯¯¯U(2)α=−M0αβuβ⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathζ(1)⋅˙d+¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯M(1)αβZ0βγU(1)γ. (34)

Alternatively the expression can be derived directly from Eq. (3.12) by use of . Here we can use

 \boldmathζ(1)=d⋅\boldmath∇% \boldmathζ∣∣0,\boldmathM(1)=d⋅\boldmath∇\boldmathM∣∣0, (35)

where indicates the gradient operator in configuration space, and the notation indicates that the value of the gradient is taken at . The time average in the first term in Eq. (4.6) can be expressed as

 uβ⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathζ(1)⋅˙d=12Re[−iωd∗ω⋅Dβ∣∣0⋅dω], (36)

with derivative friction matrix

 Dβ=\boldmath∇fβ,fβ=\boldmathζ⋅uβ, (37)

as introduced earlier 16 ().

In the second term in Eq. (4.6) we use the identity

 ZαγMγβ=δαβ (38)

to show that

 \boldmath∇Mαβ=−MαγgγδMδβ (39)

 gβγ=\boldmath∇Zβγ=Dβ⋅uγ. (40)

The first order velocity is eliminated by use of Eq. (4.4). Then the second order mean swimming velocity can be expressed as

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

with frequency-dependent matrix given by

 Vα(ω)=Mαβ˘Dβ(ω), (42)

with reduced derivative friction matrix

 ˘Dβ(ω)=Dβ−gβγYγδ(ω)fδ(ω), (43)

 \boldmathY(ω)=[−iω\boldmathM+%\boldmath$Z$]−1, (44)

and impedance vector

 fδ(ω)=(−iωm+\boldmathζ)⋅uδ. (45)

The matrix has the property

 ˘Dβ(ω)⋅uα=0. (46)

Since depends only on relative coordinates , and hence

 uα⋅Dβ=0,uα⋅gβγ=0. (47)

As a consequence

 uα⋅Vβ(ω)=0,Vα(ω)⋅uβ=0. (48)

These properties generalize those derived earlier at zero frequency 16 ().

To second order in the displacements the last term in Eq. (3.21) can be written as

 Uαuα⋅ddt(m⋅U)≈12ddt(U(1)αuα⋅m0⋅U(1)βuβ)+U(1)αuα⋅m0⋅¨d. (49)

In the average over a period the first term on the right does not contribute. Hence for periodic motion the mean second order rate of dissipation can be expressed as

 ¯¯¯¯¯¯¯¯¯D(2)=1τ∫τ0[˙d⋅\boldmathζ0⋅˙d+U(1)α˙d⋅f0α−U(1)αuα⋅m0⋅¨d]dt. (50)

Substituting Eq. (4.4) we find

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

with the complex matrix

 P(ω)=\boldmathζ0−Y0αβ(ω)f0α(ω)f0β(ω). (52)

The matrix is symmetric and has the properties

 uα⋅Re[P(ω)]=0,Re[P(ω)]⋅uα=0. (53)

The properties Eq. (4.20) and (4.25) allow us to reduce the dimension of the matrix description by three by introducing center and relative coordinates.

## V Velocity matrices and power matrix

The transformation given by Eqs. (3.13) and (3.14) can be used to reduce the calculation of the second order swimming velocity and rate of dissipation to one in relative space. The matrices and are transformed to

 VαT(ω)=T⋅Vα(ω)⋅T−1,PT(ω)=T⋅P(ω)⋅T−1. (54)

The first three rows of consist of and the first three columns of consist of . It follows from the properties Eq. (4.20) and (4.25) that the first three rows and columns of the transformed matrices and vanish identically. Hence in this representation we can drop the center coordinates and truncate the matrices by erasing the first three rows and columns. We denote the truncated -matrices as and and define displacements in relative space by

 (\boldmath0,\boldmathξω)=T⋅dω. (55)

With this notation the mean second order swimming velocity and rate of dissipation are given by

 ¯¯¯¯¯¯¯¯¯¯U(2)α = 12Reiω\boldmathξ∗ω⋅CT⋅^VαT(ω)⋅%\boldmath$ξ$ω, ¯¯¯¯¯¯¯¯¯D(2) = 12ω2Re\boldmathξ∗ω⋅CT⋅^PT(ω)⋅% \boldmathξω, (56)

with the matrix

 CT=[˜T−1⋅T−1]% \boldmath^. (57)

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

We rewrite the expressions in Eq. (5.3) in a more convenient form with vectors and matrices consisting of dimensionless numbers. We introduce the complex dimensionless vector

 \boldmathξc=1b\boldmathξω, (58)

where is a typical length scale, and define

 Bα = b[(CT⋅i^VαT(ω)∣∣0)s′+i(CT⋅i^VαT(ω)∣∣0)a′′], A = 1bη[(CT⋅^PT(ω)∣∣0)s′+i(CT⋅^PT(ω)∣∣0)a′′], (59)

where the superscript indicates the symmetric part, the superscript the antisymmetric part, the single prime the real part, and the double prime the imaginary part. With the scalar product

 (\boldmathξc|\boldmathηc)=N−1∑j=1% \boldmathξc∗j⋅\boldmathηcj (60)

the mean swimming velocity and mean rate of dissipation can then be expressed as

 ¯¯¯¯¯¯¯¯¯¯U(2)α=12ωb(\boldmathξc|Bα|\boldmathξc),¯¯¯¯¯¯¯¯¯D(2)=12ηω2b3(\boldmathξc|A|% \boldmathξc). (61)

The matrices and are hermitian. We call the velocity matrix and the power matrix.

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 given values of viscosity and mass density . Maximizing the quadratic form with Lagrange multiplier we obtain the generalized eigenvalue problem

 Bα\boldmathξc=λαA% \boldmathξc. (62)

The eigenvalues are real. The maximum efficiency for motion in direction is given by the maximum eigenvalue as

 EαTmax=λαmax. (63)

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.

## Vi Power and dissipation

We view the swimmer or flyer as a dynamical system in periodic motion, driven by actuating forces satisfying and . In an expansion in powers of the actuating forces the first order displacements are given by Eq. (4.3). To second order the mean swimming velocity and mean rate of dissipation are given by Eq. (5.8). In this section we relate the mean power supplied by the actuating forces to the mean rate of dissipation.

The instantaneous power supplied by the actuating forces is given by

 P(t)=E(t)⋅U(t). (64)

From Eq. (2.7)

 ddt(K+Vint)=−U⋅\boldmathζ⋅U+E⋅U. (65)

Hence we find for periodic motion

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯U⋅\boldmathζ⋅U=¯¯¯¯¯¯¯¯¯¯¯E⋅U. (66)

This shows that on time average over a period the power is fully dissipated by friction.

Since the mean thrust vanishes the usual definition of energy wastage 2 () makes no sense here. Instead we define the energy wastage as the difference

 E=¯¯¯¯P−¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathU⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯% \boldmathI. (67)

In the theory of fish swimming and bird flight the energy wastage has been associated with energy being lost to the vortical wake 9 (). We define the corresponding Froude efficiency as

 ηF=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathU⋅¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\boldmathI¯¯¯¯P. (68)

This concept may be useful for the comparison of different strokes.

## Vii Hydrodynamic interactions

In order to perform explicit calculations we must specify the form of the hydrodynamic interactions appearing in the friction matrix and the mass matrix. In practice one uses approximate expressions which are presumed to be reasonably accurate in the range of distances considered.

The friction matrix can be calculated from an approximation to the mobility matrix based on Oseen’s pair interaction 19 (). In this approximation the pair mobility tensor for the pair is given by

 \boldmathμjk=16πηaj\boldmath1δjk+18πη[\boldmath1|Rj−Rk|+(\boldmathRj−\boldmathRk)(\boldmathRj−\boldmathRk))|Rj−Rk|3](1−δjk). (69)

The mobility matrix is composed of pair tensors, and the friction matrix is its inverse.

The calculation of the mass matrix is based on potential flow theory. A dipole approximation to the mass matrix can be evaluated on the basis of an expression for the force on a sphere in a uniform flow in potential flow theory, as given by Landau and Lifshitz 22 () and by Batchelor 23 ().

In potential flow theory the flow velocity is expressed as with a scalar potential , which satisfies Laplace’s equation by incompressibility. A sphere of radius , centered at the origin and moving with velocity in a fluid at rest generates a potential

 ϕU(\boldmathr)=12a3\boldmathrr3⋅\boldmathU,r>a, (70)

corresponding to the dipole moment

 \boldmathqU=12a3\boldmathU. (71)

If the sphere is placed in a uniform flow this is modified to 24 ()

 \boldmathq=12a3(\boldmathU−\boldmathv0). (72)

For a collection of spheres in a fluid at rest at infinity the velocities and dipole moments are related in dipole approximation as

 \boldmathUj=2a3j\boldmathqj+∑k≠j\boldmathFjk⋅\boldmathqk,j=1,...,N, (73)

with dipole interaction tensor

 \boldmathFjk=\boldmathF(\boldmathRj−% \boldmathRk),\boldmathF(\boldmathr)=−%\boldmath$1$+3^\boldmathr^\boldmathrr3, (74)

where . We abbreviate Eq. (7.5) as

 U=\boldmathA−1⋅q,q=\boldmathA⋅U. (75)

The velocity of sphere after a sudden impulse from a state of rest is given by 22 (),23 ()

 (mpj+12mfj)\boldmathUj=\boldmathSj+32mfj∑k≠j\boldmathFjk⋅% \boldmathqk,j=1,...,N. (76)

Substituting from Eq. (7.7) and solving for the velocities we find by use of Eq. (7.5) and

 U=w⋅S (77)

with matrix

 w=[mp−mf+4πρ\boldmathA]−1, (78)

where the matrices and are diagonal with elements and . The approximate effective mass matrix is

 m=mp−mf+4πρ\boldmathA (79)

If the spheres are neutrally buoyant one has simply

 m=4πρ\boldmathA,w=14πρ\boldmathA−1. (80)

In our application we shall consider this case.

## Viii Three-sphere swimmer or flyer

The simplest application of the theory is to a three-sphere swimmer or flyer with the three spheres aligned on the axis, as studied by Golestanian and Ajdari 25 () in the resistive limit. The spheres move along the axis, and the and coordinates can be ignored 26 (). There are only two relative coordinates and , and the relevant parts of the matrices and are two-dimensional. The relevant parts are denoted as and . In the bilinear theory we consider a point in -space with coordinates , corresponding to the configuration of the rest system. As an example we consider the case of equal-sized spheres with , equal masses , and equal distances between centers .

For this case analytic expressions for the matrices and can be derived, but are too complicated to be presented. In the high viscosity limit the expressions reduce to those derived previously 15 (). The matrices, defined with in Eq. (5.6), depend only on the ratio and the dimensionless viscosity 27 ()

 η∗=ηωa2ρ. (81)

The two eigenvalues , as well as the corresponding eigenvectors with , also depend only on these two variables. The dependence on the viscosity is surprisingly slight over the whole range of values. The absolute value is close to unity over the whole range. In Fig. 1 we show the variation of and with for . The argument of increases slightly from at to at .

In the bilinear theory the optimal orbit in relative space is given by with and

 \boldmathξ0(t)=εaRe[\boldmathξ% +exp(−iωt)], (82)

with amplitude factor . The corresponding displacement vector in configuration space is given by

 (83)

In Fig. 2 we show a snapshot of the spheres and their velocities in the instantaneous rest frame at for and . In Fig. 1 of Ref. 15 we showed the elliptical orbit in relative space for and in the limit of high viscosity. In Fig. 3 we compare this with the corresponding orbit in the limit of zero viscosity, corresponding to small . The two plots are indistinguishable on the scale of the figure.

From the periodic displacement we can calculate the instantaneous swimming velocity as a series of harmonics from Eq. (3.8). The zeroth harmonic yields the mean swimming velocity . From and we can calculate the time-dependent rate of dissipation , and hence the mean , by use of Eq. (3.17).

In Fig. 4 we show the reduced mean swimming velocity as a function of for and . In Fig. 5 we show the reduced mean power vs. the reduced mean swimming velocity in the range . In Fig. 6 we show the efficiency as a function of . The efficiency increases monotonically with the amplitude factor. In Fig. 7 we show the time-dependence of the impetus and the center velocity as functions of time during a period for . We also plot separately the resistive contribution to the impetus. At this value of this is much smaller than the reactive part. It is noteworthy that the variations in time of the center velocity are much smaller than those of the impetus. The center velocity follows the impetus with an after-effect. In Fig. 8 we show the absolute value of the Fourier coefficients of the velocity , normalized to , for . This shows that only a small number of harmonics contribute appreciably.

The equality of mean power and mean rate of dissipation given by Eq. (6.3) can be checked numerically. The Froude efficiency , defined in Eq. (6.5), for these three values of is , respectively.

It is of interest to compare the above results with values obtained from the numerical solution of the equations of motion Eq. (2.7) with hydrodynamic interactions given by Eqs. (7.1) and (7.12) and with prescribed oscillating actuating forces. In order to stabilize the system we consider microswimmers with internal harmonic interactions 28 (). In matrix form the forces may be expressed as

 F=H⋅(R−R0)+E, (84)

with a real symmetric matrix with the property . The actuating forces can be chosen to correspond to the eigenvector with maximum eigenvalue in the problem Eq. (5.9). The first term in Eq. (8.4) represents a harmonic approximation to the direct interactions. We use harmonic interactions given by the -matrix

 H=k⎛⎜⎝−1101−2101−1⎞⎟⎠ (85)

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

 σ=kπηaω. (86)

We use the first order equations of motion,

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

to calculate the actuating forces corresponding to the optimal linear motion given by Eqs. (8.3) and (4.4). These have the property , so that the sum of actuating forces vanishes. We choose initial conditions corresponding to the rest configuration

 x1(0) = 0,x2(0)=d,x3(0)=2d, U1(0) = 0,U2(0)=0,U3(0)=0. (88)

The kinetic energy term in Eq. (2.7) makes the direct numerical solution of the equations of motion time-consuming. Instead we use an iterative procedure, neglecting the kinetic energy term in first approximation. Then the equations can be solved by a fast procedure. From the solution the kinetic forces, defined as , can be calculated as a function of time. In the next step we include the kinetic forces and repeat the procedure. The kinetic forces are small compared to the actuating forces and the solution converges after a few iterations.

In Fig. 9 we show the numerical solution of the equations of motion Eqs. (2.6) and (2.7) with forces given by Eq. (8.4) for , viscosity , stiffness , and amplitude factor for the first fifty periods. In Fig. 10 we show the mean value of the kinetic energy for successive periods. The orbit for the last period hardly differs from the ellipse given by Eq. (8.3), as shown in Fig. 11. The mean swimming velocity and the mean rate of dissipation can be calculated as time-averages over the last period. The efficiency is , equal to the value calculated from the periodic orbit by use of Eq. (3.8) for displacements with , shown in Fig. 6.

## Ix Discussion

In general the performance of a swimmer or flyer can be measured in terms of the dimensionless efficiency , defined as the ratio

 ET=ηωa2¯¯¯¯U¯¯¯¯P, (89)

where is a conveniently chosen length scale, is the mean speed and is the mean power, averaged over a period . We note that the lower bound of the efficiency vanishes, since for a periodic motion in relative phase space which is time-reversible the mean swimming velocity vanishes. A striking result of the present analysis is that for the simple three-sphere model, for which analytic calculations can be performed, the maximum efficiency is nearly independent of the dimensionless viscosity , as shown in Fig. 1, see Eq. (5.10). As a consequence the mean speed is nearly inversely proportional to the shear viscosity for given power. We expect that this is a general feature of the mechanical system defined by equations of motion of the type Eq. (2.7) with hydrodynamic interactions as detailed in Sec. VII. This explains the great advantage that birds in air have over fish in water.

A second result of the analysis is that the mean power is equal to the mean rate of dissipation, as shown in Sec. VI. There is no additional energy loss related to the rate of change of the virtual mass.

The mean speed and mean power can be evaluated for more sophisticated model swimmers or flyers by similar analysis. Elsewhere we studied longer chains with both longitudinal and transverse excitation in the resistive limit 26 (). That analysis can be extended to the full range of viscosity, based on the equations of motion Eqs. (2.6) and (2.7).

For assumed periodic displacements the velocity of the assembly can be derived from the equation of motion Eq. (3.8) by decomposition in harmonics, as demonstrated in the three-sphere model. The mean power is found as a collateral of the calculation. The full range of viscosity is covered, so that the analysis provides an interesting link between the flying of birds, the swimming of fish, and the swimming of bacteria.

## References

• (1) O. Lilienthal, Der Vogelflug als Grundlage der Fliegekunst (R. Gaertner’s Verlagsbuchhandlung, Berlin, 1889).
• (2) M. J. Lighthill, Mathematical Biofluiddynamics (SIAM , Philadelphia, 1975).
• (3) S. Childress, Mechanics of swimming and flying (Cambridge University Press, Cambridge, 1981).
• (4) C. J. Pennycuick, Modelling the Flying Bird (Academic Press, Burlington MA, 2008).
• (5) H. Tennekes, The Simple Science of Flight (MIT Press, Cambridge MA, 2009).
• (6) M. J. Lighthill, ”Note on the swimming of slender fish”, J. Fluid Mech. 9, 305 (1960).
• (7) T. Y. Wu, ”Swimming of a waving plate”, J. Fluid Mech. 10, 321 (1961).
• (8) T. Y. Wu, ”On Theoretical Modeling of Aquatic and Aerial Animal Locomotion”, Adv. Appl. Mech. 38, 291 (2001).
• (9) T. Y. Wu, ”Fish Swimming and Bird/Insect Flight”, Annu. Rev. Fluid Mech. 43, 25 (2011).
• (10) J. A. Sparenberg, ”Survey of the mathematical theory of fish locomotion”, J. Eng. Math. 44, 395 (2002).
• (11) J. Gray, ”Studies of animal locomotion”, J. Expl. Biol. 13, 192 (1936).
• (12) G. I. Taylor, ”Analysis of the swimming of microscopic organisms”, Proc. Roy. Soc. London A 209, 447 (1951).
• (13) J. Lighthill, ”Flagellar Hydrodynamics: The John von Neumann Lecture, 1975”, SIAM Review 18, 161 (1976).
• (14) E. Lauga and T. R. Powers, ”The hydodynamics of swimming microorganisms”, Rep. Prog. Phys. 72, 09660 (2009).
• (15) B. U. Felderhof, ”Swimming of an assembly of rigid spheres at low Reynolds number”, Eur. Phys. J. E 37, 110 (2014).
• (16) B. U. Felderhof, ”Efficient swimming of an assembly of rigid spheres at low Reynolds number”, arXiv1504.05794[physics.flu-dyn].
• (17) J. Lighthill, ”Large-amplitude elongated body theory of fish locomotion”, Proc. R. Soc. Lond. B 179, 125 (1971).
• (18) H.-B. Deng, Y.-Q. Xu, D.-D. Chen, H. Dai, J. Wian, and F.-B. Tian, ”On numerical modeling of animal swimming and flight”, Comput. Mech. 52, 1221 (2013).
• (19) J. Happel and H. Brenner, Low Reynolds number hydrodynamics (Noordhoff, Leyden, 1973).
• (20) J. Lighthill, An Informal Introduction to Theoretical Fluid Mechanics (Clarendon Press, Oxford, 1986).
• (21) B. U. Felderhof, ”Comment on: ”Circular motion of asymmetric self-propelling particles””, Phys. Rev. Lett. 113, 029801 (2014).
• (22) L. D. Landau and E. M. Lifshitz, Fluid Mchanics (Pergamon, Oxford, 1987).
• (23) G. K. Batchelor, An Introduction to Fluid Dynamics (Cambridge University Press, Cambridge, 1967).
• (24) B. U. Felderhof, ”Virtual mass and drag in two-phase flow”, J. Fluid Mech. 225, 177 (1991).
• (25) R. Golestanian and A. J. Ajdari, ”Analytic results for the three-sphere swimmer at low Reynolds number”, Phys. Rev. E 77. 036308 (2008).
• (26) B. U. Felderhof, ”Optimization of flagellar swimming by a model sperm”, arXiv:1412.3937.
• (27) A. S. Sangani and A. Prosperetti, ”Numerical simulation of the motion of particles at large Reynolds numbers” in Particulate two-phase flow ed. M. C. Roco (Butterworth-Heinemann, Boston, 1993).
• (28) B. U. Felderhof, ”The swimming of animalcules”, Phys. Fluids 18, 063101 (2006).

## Figure captions

### Fig. 1

Plot of the eigenvalue and the absolute value for the corresponding eigenvector of the three-sphere model with as functions of the dimensionless viscosity .

### Fig. 2

Snapshot of the spheres and their velocities in the rest frame at for the motion given by Eq. (8.2) with and .

### Fig. 3

Plot of the elliptical orbit in the plane corresponding to Eq. (8.2) with and (solid curve). We also plot the elliptical orbit for the high viscosity limit (dashed curve). The two curves cannot be distinguished on the scale of the figure.

### Fig. 4

Plot of the reduced mean swimming velocity for and as a function of the amplitude as calculated from Eq. (3.8) for displacements given by Eq. (8.3).

### Fig. 5

Parametric plot of the reduced mean swimming power vs. the reduced mean swimming velocity for , and .

### Fig. 6

As in Fig. 4 for the efficiency .

### Fig. 7

Plot of the impetus (solid curve), the resistive contribution to the impetus (short dashes) and of the center velocity (long dashes) as functions of time during a period for the three-sphere swimmer for , and amplitude factor .

### Fig. 8

Plot of the absolute values of the Fourier coefficients of harmonics of frequency , normalized to , of the center velocity of the three-sphere swimmer for and with amplitude factor .

### Fig. 9

Plot of the positions of the three spheres found by numerical integration of the equations of motion Eq. (2.7) for , for fifty periods of time.

### Fig. 10

Plot of the mean value of the kinetic energy for successive periods corresponding to Fig. 9.

### Fig. 11

Plot of the orbit during the last period of Fig.9 (solid curve) compared with the elliptical orbit given by Eq. (8.2) (dashed curve).