Microswimming with inertia
Jayant Pande∗1,2 , Kristina Pickl∗2,3 , Oleg Trosman1,2 , Ulrich Rüde1 and Ana-Sunčana
Smith2,3,4†
2 Cluster of Excellence: EAM, Friedrich-Alexander University Erlangen-Nürnberg, Nägelsbachstraße 49b,

91054 Erlangen, Germany
3 PULS Group, Department of Physics, Friedrich-Alexander University Erlangen-Nürnberg, Nägelsbachstraße

49b, 91054 Erlangen, Germany

arXiv:1603.04633v2 [cond-mat.soft] 6 Nov 2016

1 Chair for System Simulation, Friedrich-Alexander University Erlangen-Nürnberg, Cauerstraße 11, 91058

Erlangen, Germany
4 Division of Physical Chemistry, Rud̄er Bošković Institute, Bijenička cesta 54, Zagreb, Croatia
∗ These authors contributed equally to the work.
† Author for correspondence. E–mail: smith@physik.fau.de.

Abstract
Microswimmers, especially in theoretical treatments, are generally taken to be completely inertia-free, since inertial effects on their motion are typically small and assuming
their absence simplifies the problem considerably. Yet in nature there is no discrete break
between swimmers for which inertia is negligibly small and for which it is detectable.
Here we study a microswimming model for which the effect of inertia is calculated explicitly in the regime of transition between the Stokesian and the non-Stokesian flow
limits, which we term the intermediate regime. The model in the inertialess limit is the
bead-spring swimmer. We first show that in the intermediate regime a mechanical microswimmer exhibits damped inertial coasting like an underdamped harmonic oscillator.
We then calculate analytically the swimmer’s velocity by including a mass-acceleration
term in the equations of motion which are otherwise based on the Stokes flow. We show
that this hybrid treatment combining aspects of underdamped and overdamped dynamics provides an accurate description of the motion in the intermediate regime, as verified
here by comparison to simulations using the lattice Boltzmann method, and is a significant improvement over the results from the inertialess theory when either the mass of the
swimmer or the forces driving its motion is/are large enough.
Keywords: microswimming, bead-spring swimmer, inertial effects, lattice-Boltzmann
method, analytical calculation, perturbation theory

1

Introduction

The motion of microswimmers, whether natural ones like E. coli [1] and C. reinhardtii [2]
or artificial ones like various proposed micro-machines [3–7], is dominated by the viscous
forces exerted by the surrounding fluid, with the inertia of the swimmers playing practically
no role. For these organisms the Reynolds number (Re) is typically very small, usually 10−2
or smaller. Because of this, a common approach in the literature to model their motion is to
assume that the Re vanishes [8–12]. In this limit the non-linear terms in the Navier-Stokes
equation for the fluid flow can be eliminated, resulting in the Stokes equation, which makes
the analysis of the motion much easier [9]. The swimmers are then said to be in the “Stokes
regime”.
1

Larger organisms such as small insects swim in what may be termed the “intermediate
regime”, at which both the viscous forces and the inertial forces are important in directing the
motion [13]. Their Reynolds numbers typically range between 1 and a few hundred. Finally,
organisms on the scale of meters, such as large fish, swim at Re & 103 , where the inertial
forces dominate the motion, and the coasting of the swimmer during the non-active part of
each cycle is an integral part of the swimming strategy [14]. This is the domain of the “nonStokes regime”.
The above demarcation of the three regimes is only loosely defined, with the specific
details of a swimmer controlling its Reynolds number as well as the importance of inertial
effects in determining its motion. In practice Re < 1 is often taken to be a sufficient condition
for assuming the Stokes regime, but the validity of such an assumption is not established.
Moreover, inherently non-Stokesian methods such as lattice Boltzmann (LB) simulations are
often used to simulate microswimmers with Re < 1, and any inertial effects if present in the
simulation are typically ignored in the analysis.
In this work we consider the motion of microswimmers beyond the bounds of the Stokes
regime, i.e. when effects of inertia are no longer negligible. As discussed above, the assumption of Stokes flow simplifies the swimming problem a lot, as the inconvenient non-linear
and time-dependent terms in the Navier-Stokes equation disappear. This assumption results
unavoidably in some non-physical effects, such as a direct proportionality of the swimming
speed to the applied force and an instantaneous response of the fluid, up till infinity, to any
disturbance within it. These effects can usually be ignored, firstly in the interest of the linearity of the resulting Stokes equation which is quite evidently crucial for analytic treatments of
the swimming problem, and secondly because these effects turn out not to cause significant
errors in the flow description.
Inclusion of the inertial contributions to the flow is a very daunting task in general. To
bring it within the realm of tractability, we restrict our attention to the lowest order effects of
inertia on our swimmer. In other words, we study the swimming regime where inertial or nonStokesian effects first emerge and impart quantitative, and possibly qualitative, differences
to the motion when compared to that in purely Stokes flow. Our approach will be to take
inertia into account modifying the methods of a Stokes flow calculation. We will do this
in two ways. Firstly, we will provide a heuristic scheme for identifying the non-Stokesian
regime of motion by considering the relaxation of a swimmer from an initially stably-moving
configuration. In the Stokesian approximation, this relaxation is instantaneous, and the degree
of non-instantaneity in it will speak to the non-Stokesian aspect of the motion. We will show
that the larger the mass of the swimmer, the slower does it have to be for the relaxation to
display an inertial component. Secondly, we will perform a full calculation of the swimming
velocity of our bead spring swimmer model with non-negligible inertia, by including a mass
acceleration term in the equations of motion governing the swimming. This shall confirm the
finding that for larger swimmer masses the inertial or the non-Stokes regime begins at smaller
swimming speeds (i.e. smaller driving forces). In both cases the theoretical predictions will
be supported by lattice Boltzmann simulations.

2

Swimmer model

Our swimmer, based on the Najafi-Golestanian design [4], consists of three rigid spheres
of equal masses m and radii r connected by two linear springs of equal stiffness constants
2

Figure 1: Three-sphere swimmer model, based on [4].

k and mean rest lengths l. In the simulations, additionally, we include angular springs of
high stiffness constants to ensure linear and rotational alignment of the swimmer [15]. Timeirreversibility in the stroke, as required by the Scallop Theorem for net propulsion at zero
Reynolds number [3], is achieved via the sinusoidal driving forces
Fd1 (t ) = − ( Fd2 (t ) + Fd3 (t )) ẑ,
Fd2 (t ) = −a sin (ωt ) ẑ, and
Fd3 (t ) = b sin (ωt + α) ẑ, for t ≥ (π − α) /ω,

(2.1)
(2.2)
(2.3)

that are applied to the centers of mass of the spheres along the main axis of the swimmer, the
x-axis. The parameters a and b denote the amplitudes of the sinusoidal forces, ω is the force
frequency, T is the cycle period (so that T = 2π/ω), α is the relative phase shift in the forces
and t is the instantaneous time. The force Fd3 on the right sphere is kept zero for the first
(π − α) /ω time steps, so that it starts continuously from a value of zero as that aids in the
stability of simulations; this has no bearing on the model itself. Due to Eq. (2.1) the criterion
of force neutrality is always satisfied. As a result of the forces, the trailing and leading arms
perform sinusoidal motion with oscillation amplitudes d1 and d2 , respectively (Fig. 1).

3

Simulation framework

The simulation system consists of two parts, the parallel software framework waLBerla [16–
18] for simulating the fluid, and the rigid body engine pe [19] for the bodies within it.
waLBerla is based on the lattice Boltzmann method (LBM) [20, 21] and uses a D3Q19
model [22] for three-dimensional space discretisation employing 19 particle distribution functions (PDFs) per lattice site. For the collision operator, Ginzburg’s two-relaxation time model
[23] is applied, which splits the PDFs, the used incompressible equilibrium distribution function, and the relaxation parameter into symmetric and antisymmetric parts. For the simulations, the best accuracy at the solid walls is achieved with the symmetric relaxation parameter
set to 1/τ and the antisymmetric one set to 8(2 − 1/τ)/(8 − 1/τ) [23], where τ denotes the
relaxation time. The software utilizes the message passing interface (MPI) and is optimized
for scalable and efficient execution on the fastest supercomputers available [24].
The massively parallel rigid body engine pe [19] handles the dynamics of the bodies using
Newton’s equations of motion, and includes mechanisms for resolving frictional rigid body
collisions and modelling external forces like gravity or the sinusoidal forces driving the swimmer. Potential contacts are handled in this work using the parallel Discrete Element Method
(DEM) [25]. In addition to these soft-contact models it also supports, e.g., a contact resolution
3

based on hard-contact models [26]. The algorithm is augmented with an MPI communication
strategy that can handle general pairwise spring-damper systems [15]. To avoid problems
resulting from non-local communication among processes resulting from extended springs,
communication is restricted to those pairs of processes on which objects interact. Global
communication is avoided since it might lead to a deterioration of parallel performance and
scalability on modern supercomputers.
The interactions between the swimmer and the fluid, as well as between different bodies
within the swimmer, are modelled by a four-way coupling scheme [15, 27, 28]. The rigid
bodies of the swimmer overlap with the cells of the LBM grid and are marked as obstacles
within the fluid. Those that interface to the fluid are specified to have a moving boundary
condition [29–31]. The fluid couples to a rigid body via the momentum exchange method [29,
30, 32] that takes into account the instantaneously acting hydrodynamic forces from the fluid
on the rigid body. A detailed description of the algorithms and methods used in the simulation
systems can be found in Pickl et al. [15, 27].

4

Stokes regime: velocity and stroke

To first test the suitability of the LBM simulation system to reproduce swimming within the
Stokes regime, we compare the simulation results with purely Stokesian theory. We have
described this theory in previous work [33], where we have determined the velocity of the
swimmer by considering the effect on each sphere of the different forces it faces, these being
the spring forces, the driving forces and the hydrodynamic forces. For small driving forces, the
oscillations of the spheres are small too, and one obtains a coupled system of ordinary differential equations for the sphere positions as functions of time, due to the linear velocity-force
relationship in Stokes flow. This system is solved in a perturbative manner with the oscillation
of the swimmer arms being the variable of perturbation, and this leads to a swimming velocity
given by (in the notation employed in this paper) [33]
h 



i
7r 2 a2 − b2 κr − ab κ2 + 12r2 sin α
ẑ .
(4.1)
vforce =
24l2 π2 ν2 ρ2 ω (κ4 + 40κ2 r2 + 144r4 )
Here vforce denotes the swimming velocity calculated using the forces acting upon the swimmer, ν is the kinematic viscosity of the fluid, ρ is its density, and κ is a constant defined for
convenience as κ = k/(πνρω). This calculation assumes that there is no slip between the
spheres and the fluid, and that the distances between the spheres are much larger than the other
length scales in the problem.
In response to the driving forces, the swimming stroke induced is sinusoidal, and can be
expressed as
L1 (t ) = l + d1 cos(ωt + δ1 ),
L2 (t ) = l + d2 cos(ωt + δ2 ),

(4.2)

where Li is the instantaneous length of the ith arm, and di and δi are the amplitude and the
phase of its oscillation. Once this form of the swimming stroke is adopted, and the earliermentioned conditions of no fluid slip, large sphere separations and small arm-length oscillations are assumed, then the swimmer’s velocity can be written as [34]
vstroke = Gd1 d2 ω sin (δ1 − δ2 ) ẑ.
4

(4.3)

Here vstroke is the swimming velocity (with the subscript referring to the stroke-dependence
of its calculation), and G is a geometrical constant which for equal sphere radii r and equal
mean arm-lengths l is given by G = 7r/(24l2 ).
The two approaches of finding the swimming velocity–i.e. assuming known driving forces
and known strokes–may be reconciled by determining the various stroke parameters for the
assumed force protocol, and then using the now-known strokes to find the velocity as in
Eq. (4.3). The stroke parameters are found to be
v
t
h
i
a2 κ2 + 4r2 4a2 + b2 − 4ab cos α + 4abκr sin α
,
d1 =
π2 ν2 ρ2 ω2 (κ4 + 40κ2 r2 + 144r4 )
(4.4)
v
t
h
i
2
2
2
2
2
b κ + 4r a + 4b − 4ab cos α − 4abκr sin α
d2 =
.
π2 ν2 ρ2 ω2 (κ4 + 40κ2 r2 + 144r4 )
The stroke phase difference sin(δ1 − δ2 ) is presented in the Supplementary Information (SI)
due to the heft of the expression involved.
Recasting now the stroke-based velocity form of Eq. (4.3) into the force-based one of
Eq. (4.1), we find
7(4l − 3r )(4l − 7r )
vstroke ≈
vforce .
(4.5)
4l(28l − 45r )
Clearly, to the lowest order in r/l, the two velocities match each other perfectly.

5

Relaxation of inertial swimmers

We wish to observe and identify non-inertial effects in our swimmer’s motion by considering
its relaxation in the fluid when all forces on it cease. For this purpose, we make the beads
in the swimmer much heavier than the surrounding fluid, and let the swimming motion be
much faster, with the expectation that both these contributions aid in violating the assumptions
of Stokes flow. We therefore run simulations of swimmers with spherical beads using the
waLBerla-pe framework, where the radius of each sphere is 6∆x and the rest length of each
arm is 32∆x. The cycle period is kept at 8000 time steps to aid simulation accuracy. The
phase shift α is set to 0.5π.
The simulations are run in four sets, with the mass of each sphere in the four sets being
20000, 30000, 40000 and 50000 lattice units, respectively. Note that these masses imply that
the sphere densities are about 22 to 55 times larger than the density of the surrounding fluid,
respectively, in the four sets, meaning that the beads are not neutrally buoyant. This is not a
problem since gravity plays no role in our simulations.
Within each simulation set we increase the driving force amplitudes successively and
check first the swimming veolcity of the swimmer in the steady state, and then its relaxation to zero velocity when the driving forces are suddenly switched off. The range of driving
force amplitudes checked is slightly different in the different sets, because (as we shall show)
the critical force amplitude values where non-inertial effects first become visible depend on
the swimmer’s mass. The force ranges chosen in each simulation set are those at which the
Stokes regime, the non-Stokes regime, and an intermediate regime each become visible.
We find that for all the simulation sets, the Stokesian theory of section 4 gives the correct
swimmer velocities in the steady state only for low driving force amplitudes, and diverges
5

[10-4 l. u.]

3

2

1

(A = B) [10-1 l. u.]
Figure 2: Velocity of a swimmer with inertia for different force amplitudes, from simulations and theory.

from the simulations as the forces increase. In each case, the velocities found from the simulations are higher than those predicted by the theory. As an example, Fig. 2 shows the velocity
curves from the Stokesian theory and the simulations for a sphere mass of 40000 lattice units.
This matches the expectation of a continuous transition from the Stokes regime to the nonStokes one, which we have loosely marked as blue and yellow regions respectively in Fig. 2.
In between lies an intermediate regime, marked in pink in Fig. 2, where the transition occurs.
We now locate this transition more precisely.

5.1

Underdamped relaxation of the swimmer

Our bead-spring swimmer can be viewed as a system of connected harmonic oscillators, which
are coupled through both the middle sphere and the surrounding fluid, as well as damped by
the fluid. Since the drag force is dominant in Stokes flow, we postulate that the Stokes and
the non-Stokes regimes are characterised by overdamping and underdamping, respectively,
in the motion of this oscillator system. For an underdamped driven harmonic oscillator, the
trajectory is given by
x(t ) = x0 e−γt cos (ωr t − αr ) ,
(5.1)
where x0 is the maximum amplitude, γ is the damping constant, and ωr and αr are constants of
oscillation. For underdamped motion of a coupled system such as ours, and of any mechanical
microswimmer in general, it is difficult to specify the different parameters in Eq. (5.1), yet the
damping coefficient γ may still be identified. To do this, consider the swimmer’s relaxation
from an initially steady state when the driving forces are switched off. In the steady state, the
swimmer faces the Stokes drag force FSt = −6πηreff u where reff is its effective hydrodynamic
radius. When the driving forces vanish, then the body stops instantly if inertial effects are
discounted, but in the presence of inertia the body exhibits coasting, and its velocity decreases

6

A = B = 0.280
1.0

[10-3 l.u.]

1.0

[10-4 l.u.]

[10-5 l.u.]

A = B = 0.056

0.5

A = B = 0.560

-1.0

Figure 3: Relaxation of a swimmer for different driving force amplitudes.

continuously as
FSt = mdu/dt = −6πηreff u.
6πηreff
⇒ u = Ce−γt , with γ =
and C a constant.
m

(5.2)

Due to the use of the Stokes drag force FSt in obtaining γ, Eq. (5.1) with said γ describes
the relaxation only in the intermediate regime between the Stokes and the non-Stokes ones.
Therefore, by fitting the relaxation curve of the swimmer with Eq. (5.1), the intermediate
regime can be identified. This can then also be used to determine the swimmer’s effective
hydrodynamic radius reff , by using Eq. (5.2) to find reff once γ has been identified from the
fit to Eq. (5.1). This procedure can be used for other mechanical microswimmers such as in
[3, 6] which can be viewed as oscillators.

5.2

Identification of intermediate regime

We now check the relaxation curves obtained from the simulations of our swimmers when the
driving forces are turned off after the steady state has been reached. We find that Eq. (5.1) at
times underestimates and at times overestimates the damping seen from the simulations. As
illustration, Fig. 3 shows the relaxation obtained from simulation (blue dotted curve) and the
fit to this curve using Eq. (5.1) (orange solid curve) for three different force amplitude values
for the simulation set with sphere mass = 40000 lattice units. In each case, the fit parameters
are chosen such that they minimise the error to the simulation curve in the initial part of the
relaxation, i.e. from the initial point on the left to the first local minimum. It may be seen
that the middle simulation, with a force amplitude of 0.28 lattice units, shows near-perfect
agreement with the theoretically predicted relaxation curve based on Eq. (5.1). In contrast, in
the plot on the left, with a force amplitude of 0.056 lattice units, the swimmer’s damping is
underestimated by the fitting curve, while in the plot on the right with a force amplitude of
0.56 lattice units, the fitting curve from Eq. (5.1) overestimates the damping.
The different extents of agreement between the simulations and the theoretical fitting
curves for different force amplitudes are understandable, since the theoretical fit is only expected to work well in the regime which shows characteristics of both Stokesian and nonStokesian motion (due respectively to Eq. (5.2) which depends on the Stokes drag law and
Eq. (5.1) which assumes underdamped motion which is a non-Stokesian effect). This means
that the relaxation from the simulation should match the theory the most perfectly in the intermediate regime. Such a reasoning shows the three simulations plotted in Fig.3 to lie, from
left to right, in the Stokes, intermediate, and non-Stokes regimes, respectively.
7

[10 l. u.]

(A = B) [10-1 l. u.]
Figure 4: Effective radius of a swimmer for different force amplitudes, as determined by relaxation after
force cut-off. Only the simulations in the intermedite regime are expected to reproduce the theoretically
predicted value of 12.3 lattice units (marked by a dashed red line).

There is another way to check such a determination of the three regimes, by checking
the effective radius reff of our swimmer, calculated in Appendix A. In all the simulations that
we are here considering, we have λi = 6∆x and li = 32∆x, which gives a theoretically
expected reff value of reff = 12.3∆x. The corresponding values for each simulation can
be found by combining Eqs. (5.1) and (5.2). Focussing first on the case of sphere mass =
40000, we plot the reff value obtained from each simulation in Fig. 4, with the black circles.
The theoretical value of 12.3∆x is marked with the horizontal dashed red line. It is clear
that there is an excellent agreement between the theoretical and the simulation values, for a
force amplitude of 0.28 lattice units, i.e. precisely the force amplitude which we identified as
marking the intermediate regime on the basis of the accuracy of the underdamped oscillator
description of the swimmer’s relaxation. In general, the reff value found from simulations
increases monotonically as the driving forces increase, because the parameter γ found from
Eqs. (5.1) and (5.2) initially underestimates and then overestimates the actual damping in the
simulations. This is a good check of our reasoning to identify the intermediate regime.
We now consider similar curves obtained from all four sets of simulations. Fig. 5 shows
that in each case, the reff values given by the simulations increase monotonically. (The only
exceptions to this are the simulations with the very highest force amplitude values probed–see
for instance the orange curve for sphere mass = 30000 lattice units–but that is a relic of the
spheres almost colliding with each other at these high driving force values.) Moreover, as
the sphere mass increases, the curves shift to the left, and so do the respective force ranges
marking the intermediate regime. In each case, the simulations where Eq. (5.1) best fits the
swimmer relaxation are the same ones where the reff curves are closest to the theoretical
value, marked by the horizontal red dashed line in Fig. 5. What this means is that as the
mass of the swimmer increases, the non-inertial effects, which first make their presence felt in
the intermediate regime, become visible at smaller driving force amplitudes, i.e. for smaller

8

Figure 5: Effective radius of a swimmer for different force amplitudes and different mass values, as determined by relaxation after force cut-off.

swimming strokes. This is understandable, since in the Reynolds number it is the product of
the velocity and the density which comes together, and increasing either sufficiently breaks
the bounds of the Stokes flow assumption. Note that to find the true Reynolds number of a
swimmer which is not neutrally buoyant, one should replace the density of the fluid in the
definition of the Reynolds number by that of the swimmer itself.

6

Theory for inertial swimming

Having phenomenologically studied in the previous section the motion of our swimmer when
inertial effects start becoming visible, we now present a basic theoretical treatment of the
situation. The calculations are based on our theory for a swimmer with rigid beads swimming
in Stokes flow [33], and non-Stokesian effects are included in the model by adding a mass
acceleration term to the governing equations of motion. As we shall show, the theory provides
results which are in good agreement with the simulations of the swimmer with non-negligible
mass, and also confirms the three regimes and their characteristics discussed in the previous
section.

6.1

Calculation of velocity for a swimmer with inertia

We wish to calculate to the lowest non-negligible order the inertia-induced effect on the motion of our bead-spring swimmer. To begin with, we adopt the same model as in section 2,
with the sinusoidal driving forces specified by Eqs. (2.1)-(2.3). The equation of motion of the

9

ith bead can in this case be written as
Fdi (t ) + Fsi (t )

=

3
X

Qi j ( t ) v j ( t ) ,

(6.1)

j=1

where Fdi (t ) and Fsi (t ) are respectively the driving and the spring force on the ith bead, and the
introduced variables Qi j (t ) are functions of the differences R j (t ) − Rk (t ) ( j, k = 1, 2, 3) of
the bead positions. The left hand side in Eq. (6.1) is the sum of the forces on the ith bead not
counting the force applied by the fluid, and in the Stokesian description this sum is exactly
balanced by the hydrodynamic force of the fluid which opposes the motion of each bead. This
is manifestly not the case in the non-Stokes regime, and we postulate that in the latter regime
the effect of the applied forces (the driving and the spring forces) on each bead is to accelerate
the bead, in addition to neutralising the opposing force provided by the fluid which is given
by the right hand side of Eq. (6.1). In other words,
Fdi (t ) + Fsi (t )

=

3
X

Qi j (t )vmass
(t) + mi v̇mass
(t ),
j
i

(6.2)

j=1

where mi is the mass of the ith bead and the velocity of the ith bead is now written as vmass
(t ),
i
in order to highlight its mass-dependence.
Since the sum of the driving and the spring forces over the entire swimmer is still zero,
the above equation yields
3 X
3
X

Qi j (t )vmass
(t ) +
j

i=1 j=1

3
X

mi v̇mass
(t) = 0.
i

(6.3)

i=1

Eq. (6.3) represents a set of coupled homogenous ordinary differential equations for the velocities vmass
(t) of the three beads, which can be solved if a suitable form for the deformations
i
of the lengths of the two arms of the swimmer is given, a specification which allows one to
treat the coefficients Qi j (t ) as known functions of time. It may be noted that this approach is
stroke-based, unlike our force-based approach for the Stokes regime case in section 4. The
form for the armlength deformations that we assume is sinusoidal, given by
R2 (t ) − R1 (t ) = L1 (t )ẑ = (l1 + d1 cos(ωt + δ1 )) ẑ, and
R3 (t ) − R2 (t ) = L2 (t )ẑ = (l2 + d2 cos(ωt + δ2 )) ẑ.

(6.4)

Here ω is the frequency of the swimming cycles (which equals the frequency of the driving
forces). The above form for the arm length deformations is the same as that adopted by the
arms of the swimmer in the non-inertial (Stokesian) case, where it emerges as a response to the
sinusoidal driving forces, and also identical to the form assumed in [34] in the stroke-based
formulation of the swimmer model in the non-inertial case (Eq. (4.2)). In the inertial case the
armlengths may be expected to have a dependence on the masses mi , but comparison with the
armlength trajectories obtained in simulations suggests that the functions in Eq. (6.4) describe
the armlengths well even in the inertial swimming case if the three beads are identical.
The above equations can be solved fully for unequal masses mi , but here we present the
result only for the case mi = m for the sake of brevity. Using the set of equations (6.4) to
specify the coefficients Qi j (t ), Eqs. (6.3) can be decoupled to
v̇mass
(t) + β(t)vmass
(t) + γi (t) = 0,
i
i
10

(6.5)

where vmass
(t) denotes the magnitude of vmass
(t), and the coefficients β(t) and γi (t) are given
i
i
by
3

3

1 XX
β(t ) =
Qi j (t ), and
(6.6)
3m i=1 j=1


3
3
P
P


Q2 j (t ) + ( L̇1 (t ) + L̇2 (t ))
Q3 j (t ) + m (2L̈1 (t ) + L̈2 (t )) 
 L̇1 (t )




j=1
j=1

 γ1 (t ) 

3
3

P
P
1 



 .
−L̇1 (t )
Q1 j (t ) + L̇2 (t )
Q3 j (t ) + m (−L̈1 (t ) + L̈2 (t ))
 γ2 (t )  =

3m 
j=1
j=1
γ3 ( t )


3
3
 − ( L̇1 (t ) + L̇2 (t )) P Q1 j (t ) − L̇2 (t ) P Q2 j (t ) − m ( L̈1 (t ) + 2L̈2 (t )) 
j=1

j=1

(6.7)
The equations (6.5)-(6.7) are closed by requiring the motion to be periodic in the steady
state, i.e.
vmass
(t) = vmass
(t + 2π/ω).
(6.8)
i
i
In terms of the velocities of the individual beads, the velocity vmass of the whole swimmer is
mass

v

mass

=v

2π/ω
Z X
3

ω
ẑ =
6π

0

(t)dt,
vmass
i

(6.9)

i=1

where the averaging of the bead velocities has been done over the three beads and over one
swimming cycle.
Eqs. (6.5)-(6.9) are solved for vmass by the integrating factor method, and the final closedform expression comes out to be


2π/ω
R


1
0
0
0
 Zt

−
γ
(
t
)
I
(
t
)
dt
2π/ω
i
Z X
3


I
(
2π/ω
)
ω
1
0
− γ (t0 ) I (t0 )dt0 +
 dt, (6.10)
vmass =
i


1
6π
I
(
t
)


i
=
1
1
−
0
 0

I (2π/ω)
where the function

 t

Z



I (t ) = exp  β( s)ds



(6.11)

0

has been introduced for ease of expression.

6.2

Comparison with simulations

We check the results of our calculation for the same four sets of simulations (for sphere masses
20000, 30000, 40000 and 50000 on the lattice) discussed in section 5. Fig. 6 displays the
velocity of the swimmer for each mass value as a function of the driving forces. The different
velocity curves shown are the one obtained from the simulations (labelled as vsimulation and
marked with orange triangles), the one given by vmass using the theory of section 6.1 (marked
11

Cycle-averaged velocity [10-4 LU]

(a)

Mass = 20000

4.0

3.0

vsimulation

2.0

vmass
vstroke

1.0

0.0

0.1

0.2
0.3
0.4
0.5
Driving force amplitudes (A = B) [LU]

(b)

0.6

Cycle-averaged velocity [10-4 LU]

Mass = 30000
3.0

vsimulation

2.0

vmass
vstroke

1.0

0.0

0.2

0.3
0.4
0.5
Driving force amplitudes (A = B) [LU]

Cycle-averaged velocity [10-4 LU]

(c)

0.6

Mass = 40 000
vsimulation

2.0

vmass
vstroke
1.0

0.0

(d)

0.1

0.5
0.2
0.3
0.4
Driving force amplitudes (A = B) [LU]

Cycle-averaged velocity [10-4 LU]

Mass = 50000
1.5

1.0

vsimulation
vmass
vstroke

0.5

0.0

0.1
0.3
12 0.2
Driving force amplitudes (A = B) [LU]

0.4

Figure 6: Comparison of velocity expressions from simulations (vsimulation ) and from theories which do
(vmass ) and do not (vstroke ) include the effect of inertia, for different mass values.

Mass
Mass
Mass
Mass

Error in velocity

0.15
0.10

=
=
=
=

20000
30000
40000
50000

0.05
0.00
-0.05
-0.10

0.1

0.2

0.3

0.4

0.5

0.6

Driving force amplitudes (A = B) [LU]
Figure 7: Relative error between vsimulation and vmass for different mass values.

with blue circles), and the one given by vstroke in Eq. (4.3) which does not consider the effect
of inertia (marked with green squares).
An immediate observation is that the vmass curves match much better with the vsimulation
ones than the vstroke curves do. The higher the mass the better is the agreement between the
vsimulation curves and the theoretical vmass curves. This shows that at least in the parameter
ranges that we have explored, our velocity calculation for swimming with inertia works well.
A second result is that the vmass curves initially overestimate the simulation velocities and
then, beyond some force amplitude value dependent on the mass value, begin to underestimate
them. This is most clearly seen in parts (a) and (b) of Fig. 6, for masses 20000 and 30000 on
the lattice, but is also the case for the other two mass values. To show this more clearly, we
plot in Fig. 7 the relative errors between the vsimulation and the vmass curves for the different
sphere masses as functions of the driving force amplitudes. Here one sees clearly that the
relative errors, never very large (. 15%), are initially negative and then become positive for
each mass value except for mass 50000, in which case they are small enough (< 5%) to be
within the error bounds of the LBM. Moreover, the errors get progressively smaller as the
mass of the spheres increases.
To sum up, the observations above and in section 5 lead to the following conclusions:
1. The Stokesian approaches for predicting the swimmer’s motion work well for low driving forces (or equivalently, for small swimming strokes).
2. As the driving forces/strokes get large, the swimming begins to diverge from the Stokesian predictions.
3. Combining the Stokesian approach with non-Stokesian elements–such as an underdamped relaxation of the swimmer or a mass acceleration term in the equations of
motion–leads to accurate description of the swimmer’s motion in a narrow, intermediate regime between the Stokes and non-Stokes regimes. To the left of this regime (in
a velocity vs. driving force plot), the combined approach overestimates the swimmer velocity, and to the right of the regime it underestimates the velocity, which is consistent
with the identification of the left and right regions as the Stokesian and non-Stokesian
regimes, respectively.
13

4. As the mass of the swimmer increases, the intermediate regime becomes visible at
smaller driving forces (or swimming strokes).
It may be observed that there is a small difference between the values of the driving force
amplitudes in Figs.5 and 7 at which, respectively, the reff curves equal the theoretical value
and the errors between the vsimulation and vmass curves become zero. This is to be expected
since the methods of sections 5 and 6 are quite different, with the one based on the force-free
relaxation of the swimmer, and the other on the force-induced swimming. Nevertheless, the
fact that both the methods lead to the above listed conclusions points to the validity of the
general approach in trying to identify inertial features in the motion.

7

Appendix A

For negligibly small oscillation amplitudes, i.e. to the zeroth order in di /l j (from Eq. (4.2)),
the reff for the swimmer over a cycle is the same as that of the equilibrium configuration
moving rigidly with a constant velocity. To obtain such a configuration, we let the three beads
have effective friction coefficients λ1 , λ2 and λ3 , and let constant forces λ1 F, λ2 F and λ3 F
respectively act on the three beads, with F = F ẑ. We will first show that in the steady state
these forces indeed result in a rigid equilibrium configuration, i.e. the springs are unextended.
Assume that in the steady state both the springs are extended, with the respective extensions being ∆l1 and ∆l2 . Then, the forces on the left, middle and right spheres, apart from the
hydrodynamic forces, are respectively (λ1 F + k∆l1 ), (λ2 F + k∆l2 − k∆l1 ), and (λ3 F − k∆l2 ),
all in the +z-direction. Since we are considering the steady state, all the three spheres are asFid+s
. Therefore, we have
sumed to be moving with the same speed, which equals
6πηλi

(λ1 F + k∆l1 )
(λ2 F + k∆l2 − k∆l1 )
(λ3 F − k∆l2 )
=
=
.
6πηλ1
6πηλ2
6πηλ3
k∆l1
k∆l2 k∆l1
k∆l2
=F+
−
=F−
.
λ1
λ2
λ2
λ3
∆l1
∆l2 − ∆l1
∆l2
⇒
=
=−
.
λ1
λ2
λ3
⇒F +

(7.1)

(7.2)
(7.3)

If ∆l1 is positive (negative), then ∆l2 must be negative (positive), and then the middle term
is negative (positive), resulting in a contradiction. Therefore, the only solution to the above
equation is
∆l1 = ∆l2 = 0
(7.4)
which is what we wanted to show.
Then the velocity of the swimmer in the steady state is
vforce

1
=
3T

ZT X
3
0

1
=
3T

vi dt

i=1



ZT X
3 
3

X
 1
λi F +
T (Ri − R j ) · λ j F dt.

6πηλi
0

i=1

j,i

14

(7.5)

As the springs are always at their rest lengths, all the terms in the above expression are timeindependent, and therefore can be taken out of the time integral. So


3 3

1 X X
F
+ 
vforce =
λ j T (Ri − R j ) · F
6πη 3
i=1 j,i

=

1
F
+ (λ1 T12 + λ1 T13 + λ2 T21 + λ2 T23 + λ3 T31 + λ3 T32 ) · F,
6πη 3

(7.6)

where Ti j is shorthand for T (Ri − R j ).
Now,
1
λ1 (T12 + T13 ) · F =
8πη
1
λ2 (T21 + T23 ) · F =
8πη
1
λ3 (T31 + T32 ) · F =
8πη

!
1
1
+
(2λ1 F).
l1 l1 + l2
!
1
1
+
(2λ2 F).
l2 l1
!
1
1
+
(2λ3 F).
l1 + l2 l2

Therefore,
!"
!
!
1
1
1
1
λ1
+
+ λ2
+
l1 l1 + l2
l1 l2
!#
1
1
+
+λ3
l1 + l2 l2
(
"
!
!
F
1
1
1
1
1
=
1 + λ1
+
+ λ2
+
6πη
2
l1 l1 + l2
l1 l2
!#)
1
1
+
.
+λ3
l1 + l2 l2

2F 1
F
+
v=
6πη
3 8πη

(7.7)

By the definition of reff , we have
v=

(λ1 + λ2 + λ3 )F
Ftotal
=
,
6πηreff
6πηreff

(7.8)

from where we get
(
"
!
!
1
1
1
1
1
1
1
=
1 + λ1
+
+ λ2
+
reff
λ1 + λ2 + λ3
2
l1 l1 + l2
l1 l2
!#)
1
1
+ λ3
+
.
l1 + l2 l2

(7.9)

For identical and equidistant beads (for which we have λi = λ and li = l), the above
expression for the swimmer effective radius simplifies to
reff =

6λl
.
2l + 5λ

(7.10)

For identical beads, the expression in Eq. (7.10), correct to the first order in λ/l and the zeroth
order in di /l, can also be obtained from the swimmer effective radius provided in Eq. (40) in
[34].
15

References
[1] H. C. Berg, E. Coli in Motion, Springer (New York), 2004.
[2] J. S. Guasto, K. A. Johnson and J. P. Gollub, Phys. Rev. Lett., 2010, 105, 168102.
[3] E. M. Purcell, Am. J. Phys., 1977, 45, 3.
[4] A. Najafi and R. Golestanian, Phys. Rev. E, 2004, 69, 062901.
[5] R. Dreyfus, J. Baudry, M. L. Roper, M. Fermigier, H. A. Stone and J. Bibette, Nature,
2005, 437, 862.
[6] J. E. Avron, O. Kenneth and D. H. Oaknin, New J. Phys., 2005, 7, 234.
[7] L. Baraban, M. Tasinkevych, M. N. Popescu, S. Sanchez, S. Dietrich and O. G. Schmidt,
Soft Matter, 2012, 8, 48.
[8] J. M. Yeomans, D. O. Pushkin and H. Shum, Eur. Phys. J. Special Topics, 2014, 223,
1771–1785.
[9] J. Happel and H. Brenner, Low Reynolds Number Hydrodynamics, Prentice-Hall Inc.,
1965.
[10] D. L. Koch and G. Subramanian, Annu. Rev. Fluid Mech., 2011, 43, 637–659.
[11] E. Lauga and T. R. Powers, Rep. Prog. Phys., 2009, 72, 096601.
[12] J. Elgeti, R. G. Winkler and G. Gompper, Rep. Prog. Phys., 2015, 78, 056601.
[13] V. Ngo and M. J. McHenry, J. Exp. Biol., 2014, 217, 72740.
[14] S. Childress, Mechanics of swimming and flying, Cambridge University Press, 1981.
[15] K. Pickl, M. Hofmann, T. Preclik, H. Köstler, A.-S. Smith and U. Rüde, Parallel Computing: Accelerating Computational Science and Engineering (CSE), IOS Press, 2014,
vol. 25, pp. 395–404.
[16] C. Feichtinger, S. Donath, H. Köstler, J. Götz and U. Rüde, J. Comput. Sci., 2011, 2,
105.
[17] C. Feichtinger, J. Habich, H. Köstler, G. Hager, U. Rüde and G. Wellein, Parall. Comput., 2011, 37, 536 – 549.
[18] H. Köstler and U. Rüde, IT - Information Technology, 2013, 55, 91–96.
[19] K. Iglberger and U. Rüde, Multibody Syst. Dyn., 2011, 25, 81–95.
[20] C. K. Aidun and J. R. Clausen, Annu. Rev. Fluid Mech., 2010, 42, 439.
[21] S. Succi, The Lattice Boltzmann Equation - For Fluid Dynamics and Beyond, Clarendon
Press, 2001.

16

[22] Y. H. Qian, D. d’Humières and P. Lallemand, Europhys. Lett., 1992, 17, 479.
[23] I. Ginzburg, F. Verhaeghe and D. d’Humières, Commun. Comput. Phys., 2008, 3, 427.
[24] C. Godenschwager, F. Schornbaum, M. Bauer, H. Köstler and U. Rüde, Proceedings of
the International Conference on High Performance Computing, Networking, Storage,
and Analysis, 2013, pp. 35:1–35:12.
[25] P. A. Cundall and O. D. L. Strack, Geotechnique, 1979, 29, 47.
[26] T. Preclik and U. Rüde, Comp. Part. Mech., 2015, 2, 173.
[27] K. Pickl, J. Götz, K. Iglberger, J. Pande, K. Mecke, A.-S. Smith and U. Rüde, J. Comput.
Sci., 2012, 3, 374.
[28] J. Götz, K. Iglberger, M. Stürmer and U. Rüde, 2010 ACM/IEEE International Conference for High Performance Computing, Networking, Storage and Analysis, 2010, p. 1.
[29] A. J. C. Ladd, J. Fluid Mech., 1994, 271, 285.
[30] A. J. C. Ladd, J. Fluid Mech., 1994, 271, 311.
[31] D. Yu, R. Mei, L.-S. Luo and W. Shyy, Prog. Aero. Sci., 2003, 39, 329.
[32] K. Iglberger, N. Thürey and U. Rüde, Comput. Math. Appl., 2008, 55, 1461.
[33] J. Pande and A.-S. Smith, Soft Matter, 2015, 11, 2364.
[34] R. Golestanian and A. Ajdari, Phys. Rev. E, 2008, 77, 036308.

17

