Teaching Stokesian Dynamics to swim, a tutorial
James W. Swan and John F. Brady
This tutorial will provide the tools needed to perform Stokesian Dynamics simulations of the motion of a finite number spherical colloidal particles constrained to form a single self-propelled body.
The code coupled to this tutorial is distributed freely as supplementary material with the article
titled “Modeling hydrodynamic self-propulsion with Stokesian Dynamics. Or teaching Stokesian
Dynamics to swim” which appears in the journal Physics of Fluids. The article itself describes the
fundamentals of the Stokesian Dynamics method and methods employed to achieve self-propulsion,
while this tutorial deals with the details of the accompanying FORTRAN code. Section 1 describes
the installation, compilation and execution of the Stokesian Dynamics code. In section 2, the input
and output files for the code are discussed and a simple simulation of sedimenting rods is performed.
The results may be compared with slender-body theory as a means of verifying the code works correctly. Section 3, discusses some algorithmic details of the code and introduces simulations of the
“pusher/puller” and “spinner” swimming bodies depicted in the POF article.
1
Installation, compilation and execution
This Stokesian Dynamics software package is archived and compressed using the “tar” and “gzip”
formats. Under Windows, the free application, 7-Zip (www.7-zip.org), may be used to extract the
package, while Linux and Mac OS X may extract it natively by navigating to the directory containing the archive via the Terminal application and typing the command: tar xvfz sd-swim.tar.gz
.
Within the extracted directory sd-swim are the FORTRAN code files for the simulation with
suffices .f90 and .f, the simulation input files, pref.in and conf.in, the simulation output
file, pos.out, Makefile used for compiling the code on Linux and Mac OS X and the directory
lubdat which contains data files storing the forces due to hydrodynamic lubrication between pairs
of nearly touching spheres of equal size. A FORTRAN compiler is needed to build this simulation.
It is intended to work with the Intel Fortran Compiler (http://software.intel.com/) though
other compilers may be employed.
Specific compilation directions for Windows machines are eschewed as they will require specification of a visual FORTRAN compiler. Whereas, in Linux and Mac OS X, the Unix underpinnings
make compilation with an arbitrary compiler straightforward. The variable COMP within the file
Makefile specifies the path to the FORTRAN compiler executable and must be made to point to
the appropriate location on each individual machine. When the correct path has been specified, enter the command make in the command line within the “Terminal” application while in the sd-swim
directory. This will compile the Stokesian Dynamics simulation and produce an executable titled
sd.x. The simulation may be run by entering the command ./sd.x into the command line. It is
important that sd.x and the directory lubdat containing important data files are always located
in the same parent directory. Otherwise, the code will fail to execute. For Windows users, the
order in which the FORTRAN files must be compiled is: glob.f90, blas.f, chol.f90, init.f90,
lub.f90, dist.f90, mob.f90, vel.f90, traj.f90, out.f90, sd.f90.
1
2
Input, output and validation
There are two input files for the Stokesian Dynamics simulation – pref.in: in which the simulation
preferences are specified and conf.in in which the initial configuration of particles is specified.
Within the preference file, the number of spherical particles may be specified by changing np. Any
external forces or torques on the constrained, self-propelled body may be specified by changing
fexth i and texth i where “i” specifies the cartesian direction along which the external force or
torque is applied.
Time and distance are dimensionless in these simulations so that a single particle will move a
unit, dimensionless distance in a unit, dimensionless time when subjected to a unit, dimensionless
force. Recall, that a single particle is subject to Stokes drag law so that its dimensional velocity is
F
,
6πηa
U=
(1)
where F is the dimensional force, η is the fluid viscosity and a is the particle’s radius. The
dimensional distance travelled by the particle in a dimensional time ∆t is
∆x =
F ∆t
,
6πηa
(2)
so that ∆x/a = 1 and F ∆t/6πηa2 = 1. Therefore, distances in the simulation are made dimensionless on the particle radius a and times are made dimensionless on the factor F/6πηa2 . If no force
is specified in lieu of specifying some particle velocities, one may take the converse perspective in
which ∆t has dimension T such that all velocities are made dimensionless on a/T . In this way, the
dimensional factor T may be defined relative to the driving velocity.
The final dimensionless running time for the simulation is specified by changing tf while the
time step over which the particle trajectories are integrated is specified by changing dt. Two
final preferences pos time and out mode specify the total number of times to output the particle
positions over the course of the simulation and the mode of output respectively. Selecting “1” for
out mode results in an output data file pos.out formatted as:
x1
x2
...
xn
−−−
y1
y2
z1
z2
θ1
θ2
φ1
φ2
ψ1
ψ2
yn
zn
θn
φn
ψn
while selecting “2” for out mode results in an output data file pos.out formatted as
t1
x1
x2
...
xn
y1
y2
z1
z2
θ1
θ2
φ1
φ2
ψ1
ψ2
yn
zn
θn
φn
ψn
Here, (xi , yi , zi , θi , φi , ψi ), are the Cartesian positions and Euler angles associated with particle i
and tj is the dimensionless time at step j through the simulation. If fewer than np coordinates
2
are specified in conf.in, the code will exit reporting an error reading from file. If more than np
coordinates are specified, the first np are read and used for the simulation. Additionally, all particles
must be a minimum of two radii apart in order to ensure that there are no particle overlaps. In
the event of an overlap, the simulation will exit and report and representative error code. The
parameter drset is specified in init.f90. It is the minimum allowed inter-particle separation
and may be adjusted to admit small overlaps in cases when the finite size time step may prove a
computational limitation.
To measure an approximation of the drag coefficient for a long, slender rod in Stokes flow, open
the file conf.in and enter the following coordinates:
0.0e0
2.01e0
4.02e0
6.03e0
8.04e0
...
0.0e0
0.0e0
0.0e0
0.0e0
0.0e0
0.0e0
0.0e0
0.0e0
0.0e0
0.0e0
where the number of entries is equal to the number of particles comprising the rod. Here, the
particles are aligned with the x-axis and spaced so that there is at least 1/100th of a particle radius
between their surfaces. In pref.in specify the desired number of particles and set fexth x and
fexth y to 1.0e0. This will cause the rod to sediment along and normal to its axis of symmetry.
Set the end time and time step for the simulation, tf and dt, to 1.0e0. Now, the distance travelled
by any particle in the rod over that first and only time step is identical to the instantaneous velocity
of the rod as a whole during that time step. As motion along the rod’s axis and normal to it are
linearly independent, the drag coefficient associated with motion in each direction (the ratio of
rod velocity to total applied force) is straightforward to determine. Execute the code, change the
number of particles and repeat to generate figure 1.
3
Algorithmic details and examples
The main loop through time in the Stokesian Dynamics simulation is located in the the file sd.f90.
The loop is organized into six discrete steps/subroutines:
• check dist in which the distances between all particle pairs are computed,
• form mob in which the far-field approximation to the grand mobility tensor (M∞ ) is formed
and inverted ((M∞ )−1 = R),
• add lub in which the lubrication contributions to the hydrodynamic resistance are added to
the inverted grand mobility tensor,
• new vel in which rigid body connectivity tensor (Σ) is calculated, the velocity of the swimmer
is calculated subject to equation POF-1 in the associated article, and the velocity of the
constituent particles are calculated,
• trajectory in which the velocities of all the particles are integrated over a time-step to
generate a new configuration for the particles,
3
Rigid rod, parallel
Rigid rod, perpendicular
0.7
πη
FN/Ua
0.6
6 :ecrof/yticolev doR
0.5
0.4
0.3
0.2
5
10
15
20
Number of particles: N
Figure 1: A rod composed of N particles with centers separated by 2.01 radii sediments parallel
and perpendicular to its axis of symmetry. The ratio of the resulting sedimentation velocity to the
total applied force on the rod, N F , is plotted.
• and out inter in which the positions of the particles in the new configuration are written to
file.
Of primary concern to the user is the subroutine new vel located in the file vel.f90. It is here
that the kinematics of a particular swimming body may be specified. The relevant location in
the routine begins and ends with the symbol !!!!. Any changes made outside this region are at
the users peril. We will model two swimmers here: the pusher/puller and the spinner, which are
implicit gate and explicit gate swimmers respectively. In both cases, we must begin by specifying
a configuration of two force-free particles separated by a distance r = 2.01a along the x-axis (i.e.
make the pertinent changes to pref.in and conf.in).
In the case of the pusher/puller, we must specify the implicit gate, E. In the code, this is
represented by the array einf which is 5N in length. Every five elements correspond to the implicit
4
gate associated with a particular particle so that the E for particle 1 is stored in einf(1:5),
the E for particle 2 is stored in einf(6:10), etc. As E is a symmetric, traceless tensor, it is
represented uniquely by only five elements which are mapped for particle 1 as einf(1:5) = (Exx −
Ezz , 2Exy , 2Eyz , 2Exz , Eyy − Ezz ). A dipolar implicit gate oriented along the x-axis (the axis of
alignment for the particle pair) takes the form: E = p(I − 3ex ex ), where p is the magnitude of the
dipole, I is the idem tensor and ex is the unit vector parallel to the x-axis. A value of p = 1 will
result in a time step which is normalized by the dipole magnitude. The form of the array einf in
the case where the first particle is subject to this implicit gate while the second particle is static is
einf(1) = -3.0e0
einf(2) = 0.0e0
einf(3) = 0.0e0
einf(4) = 0.0e0
einf(5) = 0.0e0
einf(5) = 0.0e0
einf(6) = 0.0e0
einf(7) = 0.0e0
einf(8) = 0.0e0
einf(9) = 0.0e0
einf(10) = 0.0e0
The velocity of such a swimmer may be calculated in exactly the same manner as in the previous
section: set the final simulation time and time step equal to 1.0, run the simulation and find
the displacement of the swimmer over that time step from the particle position data in pos.out.
Repeat this process for two particle, pusher/puller swimmers with larger inter-particle spacings to
reproduce the plot in figure 2. Additional data about the rate of energy dissipated by the swimmer
is stored in the variable dissipation and is computed at the end of the subroutine new vel.
Output and investigation of this variable is left to the user.
In the case of the spinner, we must specify the explicit gate, US . In the code, this is represented
by the array ugate which is 6N in length. Every six elements correspond to the explicit gate
associated with a particular particle so that the US for particle 1 is stored in ugate(1:6), the
US for particle 2 is stored in ugate(7:12), etc. The explicit gate contains the six components
of the translational and rotational velocities of each particle organized so that ugate(1:6) =
(UxS , UyS , UzS , UθS , UφS , UψS ) where θ, φ, ψ correspond to rotation about the x, y, z axes respectively.
The form of the array ugate in the case where the particles are counter rotating about the y-axis
is
ugate(1) = 0.0e0
ugate(2) = 0.0e0
ugate(3) = 0.0e0
ugate(4) = 0.0e0
ugate(5) = -1.0e0
ugate(6) = 0.0e0
ugate(7) = 0.0e0
ugate(8) = 0.0e0
ugate(9) = 0.0e0
ugate(10) = 0.0e0
ugate(11) = 1.0e0
ugate(12) = 0.0e0
5
2
ap/U :deeps gnimmiws rellup/rehsuP
0.1
9
8
7
6
5
4
3
2
0.01
3
4
5
6 7 8
2
3
0.1
4
5
6 7 8
2
3
4
5
6 7 8
1
Inter-particle separation: r/a - 2
Figure 2: A pair of particles forms a pusher/puller type swimmer with one particle subject to the
dipole characterized by E = p(I − ex ex ). The ratio of the swimmer’s velocity to the magnitude of
the applied dipole is plotted.
6
As with the dipole strength for the explicit gate, the rate of rotation here has been used to make time
dimensionless in the simulation. Calculate the translational velocity of such a swimmer, varying
the inter-particle separation to recreate the plot in figure 3. Note, that the positions of the particles
themselves are stored in the variable x which is 6N in length where every six elements correspond
to the x, y, z coordinates of a particle’s center and the three Euler angles corresponding the angle
of rotation of that particle about its center.
2
0.1
Ω
a
9
/U :deeps gnimmiws rennipS
8
7
6
5
4
3
2
0.01
2
3
4
0.01
5
6 7
2
3
0.1
4
5
6 7
2
3
4
5
6 7
1
Inter-particle separation: r/a - 2
Figure 3: A pair of counter rotating particles forms a spinner-type swimmer. The ratio of the
swimming velocity of the spinner to the rate of rotation of the particles, Ω, is plotted.
4
A final example
Constructing more complicated swimmers is a task left to the user, but some general advice regarding swimmers with explicit gaits may be appropriate. It is often the case that the coordinates
for points on the surface of a swimmer are known: x(t). Consequently, the translational velocity
7
of points on the swimmer’s surface are uS (x(t)) = x(t).
˙
From the perspective of the corresponding
(1)
discretized swimmer, the translational velocity of particle 1 located at x1 (t) is US = uS (x1 (t)).
(1)
The rotational velocity of that same particle is ΩS = 21 ∇ × uS (x1 (t)).
Consider a helical waveform: x = (R sin(ks − ωt), R cos(ks − ωt), s), where R is the helix
radius, k is the helix wavelength, s parameterizes the helix and ω is the rate at which the helix
rotates. The velocity of a point on the helix is uS (s) = (−Rω cos(ks − ωt), Rω sin(ks − ωt), 0)
or equivalently uS (x) = (−ωy, ωx, 0), where we have replaced the parameterization, s, by the
coordinates themselves (x, y, z). Therefore, the translational velocity of particle 1 located at x1 (t)
(1)
which is on the contour of the helix is simply US = (−ωy1 (t), ωx1 (t), 0). The rotational velocity
(1)
of that particle is ΩS = (0, 0, ω). It is the inversion from parameterization to coordinates that
can prove difficult in practice as in the simulation, the positions of the particles are always known.
There is no luxury of an arbitrary parameterization.
Note, that a single helix cannot translate force and torque free in response to rotating. The
corresponding swimming velocity in response to this particular gait will simply be U = −US
and Ω = −ΩS so that there is no net motion. However, it is simple to reconstruct the counterrotating helix example from the corresponding article by applying this same treatment. There, the
counter-rotation results in no net torque on the body, while translation of the swimmer balances the
propulsive thrust generated by that same rotation. The user should attempt this task in preparation
for the study of a novel swimmer.
8