Nonlinear dynamics in colloidal
model systems under confined
flow
vorgelegt von
Diplom-Physiker
Sebastian Reddig
geboren in Berlin
Von der Fakultät II - Mathematik und Naturwissenschaften
der Technischen Universität Berlin
zur Erlangung des akademischen Grades
Doktor der Naturwissenschaften (Dr. rer. nat.)
genehmigte Dissertation
Promotionsauschuss:
Vorsitzender: Prof. Dr. Martin Schoen
Erster Gutachter: Prof. Dr. Holger Stark
Zweiter Gutachter: Prof. Dr. Roland Netz
Tag der wissenschaftlichen Aussprache: 20.11.2013
Berlin 2013
D83
Contents
1 Introduction 1
2 Hydrodynamics at low Reynolds number 7
2.1 Hydrodynamic equations . . . . . . . . . . . . . . . . . . . . . . . 9
2.1.1 Continuity equation and incompressibility . . . . . . . . . 9
2.1.2 The Navier-Stokes equation for incompressible fluids . . . 10
2.1.3 The Reynolds number . . . . . . . . . . . . . . . . . . . . 11
2.1.4 The Stokes equations . . . . . . . . . . . . . . . . . . . . . 12
2.2 Solutions of the Stokes equations . . . . . . . . . . . . . . . . . . 13
2.2.1 Green’s function of the Stokes equations . . . . . . . . . . 14
2.2.2 The Oseen tensor . . . . . . . . . . . . . . . . . . . . . . . 14
2.3 Flowpastasphere .......................... 16
2.3.1 The flow past a translating sphere . . . . . . . . . . . . . 16
2.3.2 The flow past a rotating sphere . . . . . . . . . . . . . . . 18
2.3.3 Faxén’stheorem........................ 19
2.4 Hydrodynamic interactions . . . . . . . . . . . . . . . . . . . . . . 21
2.4.1 Hydrodynamic interactions at the Oseen level . . . . . . . 21
2.4.2 Method of reflections . . . . . . . . . . . . . . . . . . . . . 23
2.5 Hydrodynamic interactions between two parallel planar walls in
Poiseuilleflow............................. 25
2.5.1 The two-wall Green tensor . . . . . . . . . . . . . . . . . . 25
2.5.2 Hydrodynamic interactions of point-like particles . . . . . 29
2.5.3 Multipole expansion . . . . . . . . . . . . . . . . . . . . . 30
2.5.4 L-multipole approximation and Lubrication correction . . 34
2.6 Stokesian dynamics simulations . . . . . . . . . . . . . . . . . . . 36
i
2.6.1 Equations of motion . . . . . . . . . . . . . . . . . . . . . 36
2.6.2 Numerical integration scheme . . . . . . . . . . . . . . . . 37
3 Brownian motion 39
3.1 The Langevin equation . . . . . . . . . . . . . . . . . . . . . . . . 41
3.1.1 Fluctuation-dissipation theorem . . . . . . . . . . . . . . . 42
3.1.2 The Wiener process . . . . . . . . . . . . . . . . . . . . . . 43
3.2 The Smoluchowski equation . . . . . . . . . . . . . . . . . . . . . 46
3.3 Brownian dynamics simulation . . . . . . . . . . . . . . . . . . . . 47
3.3.1 Numerical integration scheme . . . . . . . . . . . . . . . . 48
3.3.2 The Cholesky decomposition . . . . . . . . . . . . . . . . . 49
3.3.3 Gaussian random numbers: the Box-Muller transform . . . 50
4 Semiflexible polymer in confined Poiseuille flow 53
4.1 Introduction.............................. 54
4.2 Modeling................................ 57
4.2.1 The worm-like chain model . . . . . . . . . . . . . . . . . . 57
4.2.2 The bead-spring model . . . . . . . . . . . . . . . . . . . . 59
4.2.3 Equations of motion . . . . . . . . . . . . . . . . . . . . . 60
4.2.4 System parameters . . . . . . . . . . . . . . . . . . . . . . 61
4.3 Results................................. 63
4.3.1 Polymer conformation and orientation . . . . . . . . . . . 64
4.3.2 Steady state center-of-mass distribution . . . . . . . . . . . 66
4.4 Kinetic theory for a semiflexible polymer . . . . . . . . . . . . . . 70
4.4.1 Smoluchowski equation for center-of-mass current . . . . . 70
4.4.2 Analysis of the lateral center-of-mass current . . . . . . . . 73
4.5 Summary and conclusion I . . . . . . . . . . . . . . . . . . . . . . 75
5 Nonlinear dynamics of spherical particles in confined Poiseuille
flow 79
5.1 Introduction.............................. 80
5.2 Modeling................................ 82
5.2.1 System geometry . . . . . . . . . . . . . . . . . . . . . . . 82
5.2.2 Equations of motion . . . . . . . . . . . . . . . . . . . . . 83
5.2.3 System parameters . . . . . . . . . . . . . . . . . . . . . . 83
5.3 Results................................. 84
5.3.1 Translational and angular velocity of a sphere . . . . . . . 84
5.3.2 Relative motion of two identical spherical colloids . . . . . 86
5.3.2.1 Unbound and oscillatory bound states . . . . . . 86
5.3.2.2 Dependence on particle size . . . . . . . . . . . . 90
ii
iv
1
Introduction
«Newton hat ganz Recht, wenn
er bemerkt, daß dasjenige, was
wir Gesetz in der Natur nennen,
eigentlich nicht existiert, und daß
es nur Formeln sind, die unserer
Fassungskraft zu Hilfe kommen,
um eine Reihe von Erscheinungen
in der Natur zu erklären.»
Henrich Heine
1
1. INTRODUCTION
In the last decades, technology has evolved considerably. For example, the first
computers easily filled an entire room, whereas nowadays nearly everyone carries
a much more powerful smart device in their pocket. But not only the physical
size of computer chips decreased in the recent years, also the size of microfluidic
devices changed drastically. Since one can create micrometer or sub-micrometer
sized mechanical structures in silicon and other materials it is possible to design
devices that integrate one or several laboratory functions on a single chip of only a
few square millimeters in size. These so called Lab on a chip devices find applica-
tions in many physical, biological or chemical systems, such as DNA sequencing,
polymerase chain reaction, cell sorting, cell culturing, microfabrication, and as
microreactors [1, 2, 3, 4] to name just a few examples. Furthermore, microfluidic
channels provide an ideal tool for basic research. They allow controlled studies
on the influence of confinement or mimic biological systems where confinement
is essential. Examples are the flow of red blood cells through blood vessels or
single actin filaments in the actin network of the cell cortex [5]. Most impor-
tantly, using the pressure driven Poiseuille flow through a microchannel, one can
controllably drive suspended objects out of equilibrium and thereby induce novel
and intriguing dynamic structure formation in complex fluids. One thing that all
these microfluidic devices have in common is that they crucially depend on the
precise manipulation of particles interacting hydrodynamically, such as emulsion
droplets, polystyrene microbeads, biological cells, etc. [6, 7, 8]. Therefore, an es-
sential aspect for developing and understanding such devices is the flow induced
cross-streamline migration of suspended objects in confined geometries.
Already in 1836, Jean Louis Marie Poiseuille revealed in his pioneering work
that blood flow in arterioles and venules features a layer of red blood cell-free
plasma near the vessel wall [9]. More than one hundred years later, Segré and
Silberberg investigated a dilute suspension of neutrally buoyant rigid particles
flowing in circular tube. They found that the particles equilibrates at a distance
of 0.6 times the tube radius from the centerline [10]. So, in addition to depletion
at confining walls reported in earlier studies [9, 11], migration of the particles
away from the channel’s centerline occurs. A fundamental consequence of this
famous Segré-Silberberg effect is that even if the injected particles are uniformly
distributed, they move across streamlines and are well focused between the cen-
terline and the wall of the tube while traveling along the channel. This offers po-
tential applications in filtration and separation of bioparticles, such as red blood
cells [12, 13, 14]. This effect is confirmed by further experiments [15, 16, 17, 18]
and is explained by inertial forces acting on the particles in non-zero Reynolds
number flow [19, 20, 21, 22, 23].
At low Reynolds numbers close to zero, where viscous forces dominate over
inertial forces, this effect does no longer occur. Single spherical particles just
follow the streamlines in a laminar flow as a result of the kinematic reversibility
2
of the Stokes equations [24], which is explicitly demonstrated in [18, 25]. There-
fore, at low Reynolds numbers rigid particles distribute homogeneously across
the channel due to Brownian motion, whereas at finite Reynolds numbers they
migrate to a position between the channel wall and centerline. However, in dense
colloid suspensions the lateral particle profiles also become inhomogeneous [26].
The cross-streamline migration in such dense suspensions gives rise to interesting
dynamical patterns in particulate flow, such as the characteristic phonon modes
of one-dimensional microfluidic crystals in flat microchannels [27, 28, 29] or the
collective dynamics in two-dimensional arrays [30, 31, 32, 33]. To gain a better
understanding for the rheology of such suspensions, their relative motion is es-
sential. Whereas the motion of a pair of rigid particles is well studied in bounded
and unbounded shear flow [34, 35], there are fewer studies in confined Poiseuille
flow.
Cross-streamline migration at low Reynolds numbers is also observed for de-
formable objects [36], such as a polymers [37, 38, 39], emulsion droplets [40, 41]
or red blood cells [42, 43]. Investigations of dilute suspensions of biopolymers in
confined flow, such as DNA or actin filaments [44, 45, 46, 47, 48, 49, 50, 51] reveal
a depletion of polymers near walls in planar shear flow as well as in Poiseuille
flow. In addition, migration away from the channel’s centerline is also observed,
especially for semiflexible polymers [37, 38, 52]. As a result, the concentration of
the polymer becomes bimodally distributed across the channel. Whereas the
migration away from bounding channel walls is fully understood in terms of
hydrodynamic wall-polymer interactions [53, 54], the migration away from the
centerline is explained by a position-dependent diffusivity, which gives rise to a
diffusion current away from the center [37, 52].
The focus of this work is to explore the potential of hydrodynamic interactions
to induce novel motional patterns in confined flow, and therefore to gain a deeper
understanding of relevant physical and biological processes, introduced above.
First, we investigate the cross-streamline migration of semiflexible polymers,
which are confined between two planar walls and are driven out of equilibrium us-
ing the pressure driven Poiseuille flow. To do so, we perform Brownian dynamics
simulations for a bead-spring chain with bending elasticity. Since the dynamics
of the bead-spring chain takes place at very low Reynolds numbers, where viscous
forces dominate over inertial effects, it is predominantly governed by hydrody-
namic interactions between the point-like beads. We take these interactions into
account by the two-wall Green tensor of the Stokes equations [55]. Our simula-
tions reproduce the bimodal distribution only when hydrodynamic interactions
are taken into account. In addition, we derive a Smoluchowski equation for the
center-of-mass distribution and carefully analyze the different contributions to
3
1. INTRODUCTION
the probability current that causes the bimodal distribution. We show that the
flow-induced deterministic drift current mainly determines the migration away
from the centerline, whereas diffusional currents due to a position-dependent dif-
fusivity become less important with increasing polymer stiffness.
Then, we investigate the dynamics of a pair of force- and torque-free spherical
colloids, whose diameter is comparable to the distance of the two walls. For this
situation, we can no longer work with the simplified point particles but have to
take into account their finite size. We determine the grand mobility matrix using
a multipole expansion of the force densities on the spheres and taking into account
the imposed Poseuille flow as well as the no-slip condition at both bounding walls
[55, 56]. We truncate the expansion after the fifth order in the inverse particle
separation and apply a lubrication correction [56] to take nearly touching spheres
into account. We show that the bounded Poiseuille flow gives rise to new classes
of trajectories resulting in cross-streamline migration. Two particles moving on
these new trajectories exhibit either bound or unbound states. In the first case
they oscillate on closed trajectories in the center-of-mass frame. In the second
case, they exhibit cross-swapping trajectories in addition to swapping trajectories
which were already observed in unbounded or bounded linear shear flow [34, 35].
The different classes of trajectories occur depending on the initial positions of the
two particles and their size. We present state diagrams in the lateral positions,
where we categorize the trajectories and color code the oscillation frequencies of
the bound states. Finally, we discuss how the results on the two-particle system
help to understand the stability of particle trains composed of several particles.
This work is organized as follows:
Chapter 2 In this chapter, we introduce the hydrodynamic equa-
tions, which are the basis to describe the interactions of
suspended particles. We first briefly derive the Stokes
equations, which govern the dynamics of a Newtonian
fluid at low Reynolds number. We then introduce
the concepts to describe hydrodynamic interactions of
spherical particles suspended in an unbounded fluid as
well as in a bounded fluid under the influence of a pres-
sure driven Poiseuille flow. At the end of this chapter
we also guide through the Stokesian dynamics simula-
tion method, which is used in the present work.
4
Chapter 3 In this chapter, we introduce the basic concept of Brow-
nian motion. First, we present the Langevin equation, a
stochastic differential equation describing the time evo-
lution of the suspended particles. We then introduce an
equivalent approach to consider the probability distribu-
tion of the particles, which is given by the Smoluchowski
equation. At the end of this chapter we guide through
the method of Brownian dynamics simulation.
Chapter 4 In this chapter, we present and discuss our numerical re-
sults for a bead-spring chain as model for a semiflexible
polymer in confined Poiseuille flow. After introducing
the model system, we analyze the steady-state proba-
bility distribution for the center of mass and study both
conformation and orientational order of the polymer in
the channel. Then, we derive the Smoluchowski equa-
tion for the center-of-mass distribution and investigate
the different contributions to the probability current.
Chapter 5 In this chapter, we investigate the nonlinear dynamics
of a pair of force and torque-free spherical particles in
Poiseuille flow. First we introduce the model and an-
alyze the particle trajectories for different lateral and
axial initial positions of the two particles as well as for
different particle sizes. Finally, we discuss how the re-
sults on the two-particle system help to understand the
stability of particle trains composed of several particles.
Chapter 6 In this chapter, we summarize and conclude our results.
5
1. INTRODUCTION
6
2
Hydrodynamics at low Reynolds number
«... the scientist finds his reward in
what Henri Poincaré calls the joy
of comprehension, and not in the
possibilities of application to which
any discovery of his may lead.»
Albert Einstein
7
2. HYDRODYNAMICS AT LOW REYNOLDS NUMBER
Hydrodynamic interactions play an important role whenever two or more par-
ticles move in a viscous fluid. Due to their long-range nature, the dynamics of
colloidal suspensions are governed by these interactions. Since the suspended
particles are much larger than the solvent molecules, the fluid can be considered
to be continuous on the length scale of the colloids. The mechanical state of the
fluid is described by the local flow velocity u(r, t)at the position rin the fluid
at some time t, the pressure p(r, t)and the mass density ρ(r, t). For our purpose
the mass density is considered to be constant, both spatially and in time. This
incompressibility is a good approximation for a wide range of fluids, such as water.
The remaining two variables which describe the state of the fluid are governed by
the incompressibility condition and the Navier-Stokes equation which is Newton’s
equation of motion for the fluid flow.
So, in principle, if the spatial distribution of the suspended particles is given
along with their motions, one has to solve a well defined boundary value prob-
lem for the Navier-Stokes equation to describe the flow field u(r, t). To do so,
various numerical techniques have been developed in the recent years. One kind
of simulations methods are coarse grained molecular dynamic techniques, such as
the lattice Boltzmann technique [57, 58], dissipative particle dynamics [59, 60],
multiparticle collision dynamics [61, 62] or stochastic rotation dynamics [63]. In
all these concepts the fluid is described by "effective" fluid particles satisfying the
fundamental balance equations of continuum theory. These fluid particles are on
a larger scale than single solvent molecules but still smaller than the suspended
colloids.
In the regime of low Reynolds number, where viscous forces dominate over
inertial forces, the fluid dynamics is governed by the Stokes equations, which are
both linear and independent of time. These equations can be solved using the
method of Green’s functions, which leads to the technique of Stokesian dynamics
simulation [24, 64, 65]. This concept follows a different path. Instead of calculat-
ing the flow field explicitly it solves Newton’s equation of motion for the colloidal
particles only. The hydrodynamic coupling through the fluid flow is taken into
account by the mobility matrix µor its inverse the friction matrix ζ=µ−1, which
describe the relation between the particles velocities and the hydrodynamic forces
and torques.
In this chapter we introduce the essentials of the theory as well as methods
to obtain expressions for mobility matrix in confined Poiseuille flow. At the end
of this chapter we also guide through the numerical simulation method, which is
used in the present work.
8
2.1 Hydrodynamic equations
2.1 Hydrodynamic equations
In this section we introduce the theoretical concepts to describe the fluid dynam-
ics, which are governed by the fundamental equations: the continuity equation
and the Navier-Stokes equation, which can be simplified to the Stokes equations
in the creeping flow regime.
2.1.1 Continuity equation and incompressibility
We consider a arbitrary volume Vwith the boundary surface ∂V , which is fixed
in space. The total mass Mof this volume is given by
M=Z
V
ρ(r, t)dr3(2.1)
where ρ(r, t)is the mass density, which is, in gerneral, a function of space rand
time t. Since matter is neither created nor destroyed the number of fluid molecules
flowing through the boundary ∂V must be equal to the number flowing outwards
through ∂V .
d
dtM=Z
V
∂
∂tρ(r, t)dr3=−I
∂V
j(r, t)·ndr2.(2.2)
Here, j(r, t) = ρ(r, t)u(r, t)is the flux of mass and u(r, t)denotes the flow
field.1. Using Gauss’s law on the right hand side of equation (2.2) we obtain the
well known continuity equation [24]
∂
∂tρ(r, t) + ∇·j(r, t) = 0,(2.3)
where ∇= ( ∂
∂x ,∂
∂y ,∂
∂z )denotes the Nabla operator. We consider the fluid to be
incompressible, so that the mass density is constant, both spatially and in time
ρ(r, t) = ρ. Together with the continuity equation (2.3) we obtain the following
condition for the flow field
∇·u(r, t)=0,(2.4)
which is referred to as incompressibility condition. Note that this equation (2.4)
is not sufficient to determine the flow field u(r, t), it has to be supplemented by
Newton’s equation of motion for the fluid flow, the Navier-Stokes equation.
1In the present work, the dot symbol "·" denotes the Euclidean scalar product of two vectors
aand b:a·b=b·a=aibi. The product of a matrix Aand a vector bis given by (Ab)i=Aijbj,
the product of a vector aand a matrix Bis given by (aB)i=ajBji and the product of the
matrices Aand Bis given by (AB)ij =AikBkj. Note that we used Einstein’s summation
convention.
9
2. HYDRODYNAMICS AT LOW REYNOLDS NUMBER
2.1.2 The Navier-Stokes equation for incompressible fluids
Again, we consider an arbitrary volume Vwith the boundary ∂V , which is fixed
in space. Applying Newtons’ second law to the fluid, we get
d
dtP=Z
V
∂
∂tρudr3=−I
∂V
Jn dr2+I
∂V
tdr2+Z
V
fdr3.(2.5)
Here, P=RVρudr3is the momentum of the fluid in the volume V. The
first term on the right hand sight denotes the surface integral of the momentum
current density J=ρu⊗u1, which is a tensorial quantity and describes the
momentum per area and time, ndenotes the surface normal vector. The second
term on the right hand side of equation (2.5) describes the surface forces t(r, t)
due to interactions with the surrounding fluid and the last term describe the
volume forces f(r, t), which represents external forces acting on the fluid, e.g.
gravitation. The surface forces are formally given in terms of the stress tensor σ
[24]
I
∂V
tdr2=I
∂V
σn dr2=Z
V∇σdr3,(2.6)
where we used Gauss’s integral theorem. The divergence of the stress tensor ∇σ
can be interpreted as an internal volume force density, whose origin is from the
interactions with the surrounding fluid. In general, there are two contributions to
this term. One contribution arises from gradients in the hydrodynamic pressure
∇p, which is the only contribution in the static case [24]. So, the static stress
tensor
σstatic =−p1,(2.7)
is proportional to the identity matrix 1. The other contribution to the stress
tensor arises from the relative motion of neighboring fluid elements, which can
be expanded in power series with respect to gradients in the flow velocity. If the
velocity gradient are not too large, it is sufficient to take only the linear first order
derivatives ∇⊗uand ∇·u. For incompressible fluids the divergence ∇·u= 0
vanishes and the only remaining term for an isotropic fluid is the symmetric part
(∇⊗u+ (∇⊗u)T). Here, the subscript T denotes the transpose of the matrix.
Finally the full stress tensor for an isotropic incompressible fluid is given by [24]
σ=−p1+η(∇⊗u+ (∇⊗u)T).(2.8)
1The dyadic product of two vectors aand b, denoted by a⊗b, is defined such that the
components of the second rank tensor are given by (a⊗b)ij =aibj.
10
2.1 Hydrodynamic equations
The constant ηis the (dynamic) viscosity, which is specific to each material.
Fluids which are described by this equation are referred to as Newtonian fluids.
Combining this equation together with equation (2.6), we can write equation (2.5)
as
Z
V
∂
∂tρud3+I
∂V
Jn dr2=Z
V
(f−∇p+η∇2u)dr3.(2.9)
Again, using Gauss’s integral theorem on the second term on the left hand side
of this equation together with the continuity equation (2.3), we finally obtain the
Navier-Stokes equation for incompressible fluids [24, 66, 67]
ρ ∂
∂tu+ (u·∇)u!=η∇2u−∇p+f.(2.10)
2.1.3 The Reynolds number
The magnitude of the various terms in the Navier-Stokes equation (2.10) can
vary differently depending on the physical system. In the present work we are
interested in the flow field around small particles suspended in water. Let a
denote the typical length scale of the physical system (e.g. the radius of the
particle), τthe typical time scale and uthe typical velocity. Introducing the
following rescaled variables
u0=u/u,
r0=r/a,
t0=t/τ,
p0=p/vη
a,
f0=f/vη
a2,(2.11)
we obtain the dimensionless Navier-Stokes equation
a2ρ
ητ
∂
∂t0u0+Re(u0·∇0)u0=∇02u0−∇0p0+f0.(2.12)
The dimensionless factor Re is the so called Reynolds number, which is defined
as [24, 67]
Re =ρva
η.(2.13)
This number can be interpreted as the ratio of the inertial force Finert. ∼ρa2u2
and the dissipative force Fdis. ∼ηau. If the Reynolds number is small (Re
11
2. HYDRODYNAMICS AT LOW REYNOLDS NUMBER
Reynolds number
A large whale swimming at 10m/s 300000000
A tuna swimming at the same speed 30000000
A duck flying at 20m/s 300000
A large dragonfly going 7m/s 30000
A copepod in a speed burst of 0.2m/s 300
Flapping wings of the smallest flying insects 30
An invertebrata larva, 0.3mm long at 1mm/s 0.3
A sea urchin sperm advancing the species at 0.2mm/s 0.03
A bacterium, swimming at 0.01mm/s 0.00001
Table 2.1: A spectrum of Reynolds numbers for different sized organisms (after
Vogel [68])
1) the friction terms on the right hand side of the Navier-Stokes equation are
important, whereas for high Reynolds number (Re 1) the nonlinear convective
term on the left hand side dominates, which causes turbulence. Table (2.1) shows
the Reynolds number for different sized organisms.
In the present work we are interested in spherical colloids (a= 0.025µm to
4µm) suspended in water (ρ≈1g/cm3,η≈1mPa ·s) and moving in a Poiseuille
flow with a flow velocity of u≈1mm/s. For this situation, the Reynolds number
is in the order of
Re = 2.5·10−5to 4·10−31.(2.14)
Therefore, we can neglect the convective term (u· ∇)uin the Navier-Stokes
equation (2.10) and obtain the time-dependent Stokes equation
ρ∂
∂tu=η∇2u−∇p+f.(2.15)
The flow described by this linear partial differential equation is laminar and no
turbulence occurs [24]. Also, due to the linearity of equation (2.15) the principle
of superposition is valid.
2.1.4 The Stokes equations
In section 2.1.3 we showed that for small Reynolds number the Navier-Stokes
equation simplifies to the time-dependent Stokes equation (2.15), but so far we did
12
2.2 Solutions of the Stokes equations
not specified the time scale τthat is relevant for the dynamics to be considered.
To do so, we consider the equation of motion for a free particle with mass mand
friction coefficient ζsuspended in a fluid
m˙
v+ζv= 0.(2.16)
The time that the particle takes to lose its initial momentum due to the friction
is given by the momentum relaxation time τr, which follows from equation (2.16)
as
τr=m
ζ.(2.17)
As we will see in section 2.3.1, the friction coefficient of a spherical particle is
ζ= 6πηa. So, we can write down the time-dependent Stokes equation (2.15) in
dimensionless form as
9
2
ρ
ρp
τr
τ
∂
∂t0u0=∇02u0−∇0p0+f0,(2.18)
where ρp=m/(4/3πa3)is the mass density of the suspended particle. For spher-
ical colloids the momentum relaxation time τrranges typically from 1ns to 100ns,
whereas the relevant time scale τin experimental observations, which is referred
to as the Brownian timescale, ranges from 1ms to 1s [24], hence ττr. Further-
more, we can assume that the mass densities are in the same order of magnitude
ρ≈ρp. Therefore, at low Reynolds number and on the Brownian time scale we
can neglect the inertial terms on the left hand side of the Navier-Stokes equation
(2.10). The over-damped fluid dynamics is then well described by the so called
Stokes equation [24, 64]
∇p−η∇2u=f,
∇·u= 0,(2.19)
where the last equation in (2.19) is just again the incompressibility condition
(2.4). These two equations are also referred to as creeping-flow equations.
2.2 Solutions of the Stokes equations
As we showed in section 2.1, in the creeping-flow regime the fluid dynamics is
described by the Stokes equation and the incompressibility condition (2.19). In
this section we discuss the solutions to these equations.
13
2. HYDRODYNAMICS AT LOW REYNOLDS NUMBER
2.2.1 Green’s function of the Stokes equations
We consider an external force acting only at a single point r0in the fluid. Math-
ematically, this point force is described by a delta distribution
f=δ(r−r0)F,(2.20)
where Fis the total force acting on the fluid. Due to the linearity of the Stokes
equations (2.19) the fluid flow u(r)and the pressure field p(r)at some point r
are directly proportional to the force
u(r) = T(r,r0)F,
p(r) = g(r,r0)·F.(2.21)
The matrix T(r,r)and the vector g(r,r0)are the Green’s functions of the Stokes
equations for the flow u(r)and the pressure field p(r). They connect the point
force at a point r0to the resulting flow velocity and pressure at a point r. The
velocity field of a point force is also referred to as Stokeslet [69].
For a continuous force density f(r)the flow and the pressure field at some
point rare just the superpositions of the fields resulting from the forces acting
in each point of the fluid [24]
u(r) = ZT(r,r0)f(r0)dr03,
p(r) = Zg(r,r0)·f(r0)dr03.(2.22)
Once these Green’s functions are known and the external force is specified, we
can calculate the resulting flow velocity and pressure using the equations above.
In the following we introduce expressions for Green’s functions for an unbounded
fluid. For this situation the matrix T(r,r0)is called Oseen tensor, which we will
denote by the superscript O.
2.2.2 The Oseen tensor
To calculate the Oseen tensor TO(r,r0)and the pressure vector gO(r,r0)1we
substitute the integral representation of the fluid flow and the pressure field (2.22)
into the Stokes equations (2.19). This leads to
Zh∇⊗gO(r,r0)−η∇2TO(r,r0)−1δ(r−r0)if(r0)dr03= 0,
Zh∇TO(r,r0)i·f(r0)dr03= 0.(2.23)
1Note that the superscript O in gO(r,r0)denotes the pressure vector for the unbounded
geometry.
14
2.2 Solutions of the Stokes equations
Figure 2.1: Flow field around a point force acting at the origin in the fluid. The
yand zcoordinates are rescaled using a typical length scale Wof the system.
Since this equation is valid for any arbitrary external force density f(r), the
expressions in the brackets have to be equal to zero. So, that the Green’s functions
have to satisfy the following differential equations
∇⊗gO(r,r0)−η∇2TO(r,r0) = 1δ(r−r0),
∇TO(r,r0) = 0.(2.24)
The Green’s functions for an unbounded fluid have to vanish at an infinite distance
from the point force (TO(r,r0)→0and gO(r,r0)→0as (r−r0)→ ∞). Together
with this boundary condition, equation 2.23 can be solved. Hence we obtain the
following expressions [24]
TO(r,r0) = TO(r−r0) = 1
8πη|r−r0| 1+(r−r0)⊗(r−r0)
|r−r0||r−r0|!,
gO(r,r0) = gO(r−r0) = 1
4π
r−r0
|r−r0|3.(2.25)
Note that, due to translational invariance the Oseen tensor TO(r−r0)and the
pressure vector gO(r−r0)are functions of the relative vector r−r0only. Or in
other words, for an unbounded fluid the position of the origin is not significant.
From equation (2.25), we see that the fluid flow around a point force decays with
1/(r−r0), which means it is a long ranged field. In figure 2.1, we show the flow
field around a point force (Stokeslet) acting at the origin in the fluid.
15
2. HYDRODYNAMICS AT LOW REYNOLDS NUMBER
2.3 Flow past a sphere
In section 2.2 we introduced the Oseen tensor TO(r−r0), which connects the
force acting at a point r0in an unbounded fluid to the flow velocity at some
point r. In this section we find expressions for the fluid flow around a translating
and rotating sphere. For this we consider a sphere with radius asuspended in
a quiescent fluid at a position rp. Due to attractive interactions between the
fluid particles and the core material of the colloid, the fluid is assumed to "stick"
onto the surface of the sphere. Therefore, we assume stick boundary conditions
throughout this work. Thus the fluid velocity and the particle velocity are the
same at the surface ∂Vpof the sphere
u(r) = vp+ωp×(r−rp),r∈∂Vp.(2.26)
Here, vpdenotes the translational velocity of the sphere and ωpthe angular
velocity. The particle exerts a force on the fluid, which is concentrated on its
surface. According to equation (2.22), the resulting flow and pressure fields are
given by
u(r) = I
∂Vp
TO(r−r0)fp(r0)dr02,
p(r) = I
∂Vp
gO(r−r0)·fp(r0)dr02.
(2.27)
Note that f(r0)is now a force density per unit area that a surface element of the
particle exerts on the fluid. For an unbounded fluid, with the boundary condition
at infinity
u(r)→0, r → ∞,(2.28)
the flow field around a moving sphere can be calculated exactly. In the follow-
ing, we discuss the resulting flow past a translating and rotating particle in an
unbounded fluid separately.
2.3.1 The flow past a translating sphere
We consider a sphere moving with the constant velocity vpthrough an unbounded
and quiescent fluid. The stick boundary condition at the particle’s surface is given
by
u(r) = vp,r∈∂Vp.(2.29)
16
2.3 Flow past a sphere
There are two equivalent methods to calculate the resulting flow field uaround
this sphere, via the Stokes equations (2.19) or via the integral Green’s function
representation (2.22). In the following we choose the second option. A successful
ansatz for the force density is given in [24]. There, fp(r)is chosen constant and
proportional to the velocity of the sphere
fp=c
4πa2vp,(2.30)
where cis a constant. Before inserting the force density in equation (2.27), we
first expand the Oseen tensor around the center of the sphere using Taylor’s series
TO(r−r0) =
∞
X
n=0
1
n![(r0−rp)·∇¯
r]nTO(r−¯
r)¯
r=rp.(2.31)
Here, ∇¯
rdenotes the Nabla operator for the position vector ¯
r. Putting this into
equation (2.27), we obtain for the fluid flow
u(r) = c
4πa2I
∂Vp
∞
X
n=0
1
n![(r0−rp)·∇¯
r]nTO(r−¯
r)¯
r=rp
vpdr02
=c 1 + a2
6∇2
¯
r!TO(r−¯
r)¯
r=rp
vp.(2.32)
Note that higher order terms are proportional to ∇2
¯
r∇2
¯
r, which are identical zero
due to the biharmonic equations:∇2∇2u= 0 [24], such that we obtain for an
arbitrary point force Facting at r0
∇2∇2u(r)(2.21)
=∇2∇2TO(r−r0)F= 0 ⇒ ∇2
¯
r∇2
¯
rTO(¯
r)¯
r=r−r0=0,(2.33)
for r6=r0. Defining now the translational differential operator
Lt=1+a2
6∇2,(2.34)
we can express the resulting flow field as
u(r) = cLt
¯
rTO(r−¯
r)¯
r=rp
vp
=c
6πηa "3
4
a
|r−rp| 1+(r−rp)⊗(r−rp)
|r−rp|2!
+1
4 a
|r−rp|!3 1−3(r−rp)⊗(r−rp)
|r−rp|2!
vp.(2.35)
17
2. HYDRODYNAMICS AT LOW REYNOLDS NUMBER
To satisfy the no-slip condition, the constant chas to be chosen as c= 6πηa.
Since cis known, we can now easily calculate the total force that the particle
exerts on the fluid
Fp=I
∂Vp
fpdr2=ζtvp,(2.36)
with the translational friction coefficient
ζt=1
µt= 6πηa, (2.37)
where µtis the translational mobility of the sphere. This is Stokes friction law
for translational motion of a sphere.
2.3.2 The flow past a rotating sphere
For a sphere, rotating with the constant angular velocity ωp, the stick boundary
condition on its surface reads
u(r) = ωp×(r−rp),r∈∂Vp.(2.38)
Again, the force density is chosen to be proportional to the velocity of the surface
element
fp(r) = c
4πa2ωp×(r−rp),(2.39)
where cis again a constant, which has to be chosen such that the boundary
condition (2.38) is fulfilled. Putting this into equation (2.27) yields [24]
u(r) = c
12πηa a
|r−rp|!3
ωp×(r−rp),(2.40)
Obviously, we have to choose cas c= 12πηa. Therefore, we obtain for the total
torque that the rotating sphere exerts on the fluid
Tp=I
∂Vp
(r0−rp)×fp(r0)dr02=ζrωp,(2.41)
with the rotational friction coefficient
ζr=1
µr= 8πηa3,(2.42)
where µris the rotational mobility of the sphere. This is Stokes friction law for
rotational motion of a sphere.
18
2.3 Flow past a sphere
Figure 2.2: A spherical colloid is immersed in a fluid with the flow velocity u0(r).
This fluid flow influences the translational and rotational motion of the particle.
2.3.3 Faxén’s theorem
So far, we introduced expressions for the fluid flow past a moving sphere. Here,
we are interested how a given flow field u0(r)influences the translational and
rotational motion of a suspended particle. This flow is a solution of the Stokes
equations (2.19) and may represents an imposed flow or the flow field induced
by a moving sphere. At the particle’s surface, centered at rpthe stick boundary
condition has to be fulfilled
vp+ωp×(r−rp) = u0(r) + I
∂Vp
TO(r−r0)fp(r0)dr02,r∈∂Vp.(2.43)
The left hand side of this equation represents the particle’s surface velocity and
the right hand side the resulting flow field u(r), which is the superposition of the
incident flow u0(r)in the absence of the particle and the disturbance flow due
to the present of the sphere H∂VpTO(r−r0)fp(r0)dr02. Multiplying this equation
with 1/(4πa2)and integrating over the sphere’s surface ∂Vp, we obtain [24]
vp=µtFp+1
4πa2I
∂Vp
u0(r)dr2.(2.44)
Again, Fp=H∂Vpfp(r0)dr02is the total force that the particle exerts on the fluid.
Note that, due to rotational symmetry the integration over the angular velocity
in equation (2.43) vanishes. We expand the imposed flow field around the center
of the sphere
u0(r) =
∞
X
n=0
1
n![(r−rp)·∇¯
r]nu0(¯
r)|¯
r=rp(2.45)
19
2. HYDRODYNAMICS AT LOW REYNOLDS NUMBER
and insert this into equation (2.44). As we already discussed in section 2.3.1,
higher orders terms are proportional to ∇2
¯
r∇2
¯
r, which vanish due to the biharmonic
equations [24]. So finally, we obtain Faxén’s theorem for the translational motion
of a sphere [24, 64, 70]
vp=µtFp+u0(rp) + a2
6∇2
¯
ru0(¯
r)|¯
r=rp
=µtFp+Lt
¯
ru0(¯
r)¯
r=rp,(2.46)
where Lt
¯
ris the translational differential operator (2.34) with respect to ¯
r. Note
that in the absence of the incident flow u0(r) = 0this equation reproduces Stoke’s
friction law for the translational motion (2.36).
Multiplying equation (2.43) with (r−rp)×and integrating again over the
surface of the particle, we obtain Faxén’s law for the rotational motion [24, 64, 70]
ωp=µrTp+1
2∇¯
r×u0(¯
r)|¯
r=rp
=µrTp+Lr
¯
ru0(¯
r)|¯
r=rp.(2.47)
Here, we have introduced the rotational differential operator
Lr=1
2∇×,(2.48)
with Cartesian components (Lr)ik =1
2ijk ∂
∂xj, where ijk is the Levi-Civita tensor.
In the absence of an external flow u0(r) = 0, equation (2.47) reproduces Stoke’s
friction law for the rotational motion (2.41).
20
2.4 Hydrodynamic interactions
Figure 2.3: Sketch of hydrodynamic interactions. The motion of one sphere, due
to external forces and torques acting on it induces a flow field. This fluid flow
influences the motion of the other particles and is reflected at its surfaces, so that
the translational and rotational motion of all particles are all mutually coupled via
the surrounding flow field
2.4 Hydrodynamic interactions
In section 2.3 we calculated the flow field past a moving sphere and described
how the particle’s motion is influenced by the presence of a given fluid flow.
We consider now Nspheres suspended in an unbounded and quiescent fluid.
External forces Fiand torques Tiare acting on these spheres such that they
move with translational and angular velocity viand ωi. The motion of each
sphere induces a flow field, which on the other hand influences the motion of the
other particles. In this section we introduce the theoretical framework to describe
these hydrodynamic interactions.
2.4.1 Hydrodynamic interactions at the Oseen level
We consider Nspherical particles suspended in an unbounded and quiescent
fluid. Each sphere is subjected to external forces Fisuch that they exert the
21
2. HYDRODYNAMICS AT LOW REYNOLDS NUMBER
force density fi(r)on the fluid. Due to the linearity of the Stokes equations
(2.19) the resulting flow field reads
u(r) =
N
X
j=1 I
∂Vj
TO(r−r0)fj(r0)dr02.(2.49)
If we neglect the rotational motion of the spheres, the stick boundary condition
at the surface of particle iis given by
vi=I
∂Vi
TO(r−r0)fi(r0)dr02
+
N
X
j6=iI
∂Vj
TO(r,r0)fj(r0)dr02,r∈∂Vi,(2.50)
where viis the translational velocity of particle i. Multiplying this equation by
1/(4πa2)and integrating both sides over the surface of particle i, we obtain [24]
vi=µtFi+1
4πa2I
∂Vi
N
X
j6=iI
∂Vj
TO(r−r0)fj(r0)dr02dr2,r∈∂Vi,(2.51)
If the separation of the particles is very large (|ri−rj| 2a), we can in good
approximation evaluate the Oseen tensor at center of the spheres: TO(r−r0)≈
TO(ri−rj)for r∈∂Viand r0∈∂Vj. Therefore, equation (2.51) simplifies to
vi=µtFi+
N
X
j6=i
TO(ri−rj)Fj.(2.52)
Note, in this far-field approximation the particles are treated as point like parti-
cles (Stokeslets). The motion of each particle consists of their own motion (self
mobility) and the influence of the induced flow fields by the other spheres (cross
mobility). Introducing the mobility tensor for the translational motion
µtt
ii =µt1,
µtt
ij =TO(ri−rj)(2.53)
we can describe the many particle interaction completely by
v1
.
.
.
vN
=
µtt
11 ··· µtt
1N
.
.
.....
.
.
µtt
N1··· µtt
NN
F1
.
.
.
FN
.(2.54)
As we already mentioned, the above expression for the mobility tensor is only
valid for far separated particles. They are the leading terms in an expansion
with respect to the inverse particle separation a/rij. To describe the interactions
more accurately, higher orders have to be considered. They can be calculated
iteratively by the method of reflections.
22
2.4 Hydrodynamic interactions
2.4.2 Method of reflections
The idea behind this method is to determine the flow field u0(r)by iteration,
where the flow field, induced by the motion of one particle is reflected by the
other spheres and then bounces from particle to particle. Consider an external
force Fiand torque Tiacting on sphere iat position ri. Then in zeroth order
the sphere moves with the velocity (2.36,2.41)
v0
i=µtFi,
ω0
i=µrTi(2.55)
and generates the flow field u0(r), which is the superposition of the velocity field
caused by translation (2.35) and rotation of the sphere (2.40). Therefore, the
stick boundary condition is fulfilled at this particle’s surface ∂Vi. A neighboring
particle jat position rjundergoes this fluid flow. With the help of Faxén’s
theorem we can calculate its translational and rotational velocities
v1
j=µtFj+Ltu0(r)|r=rj,
ω1
j=µrTj+Lru0(r)|r=rj.(2.56)
It is obvious, that the presence of this particle jdisturbs the flow u0(r)and an
additional contribution u1(r)has to be considered. One says the incident flow
field u0(r)is reflected on sphere j. The reflected flow can be calculated by solving
the Stokes equations (2.19) considering the boundary condition u0(r) + u1(r) =
v1
j+ω1
j×(r−rj)on ∂Vj[24].
On the other hand, the reflected flow influences the motion of particle i. The
corrected velocities are given by Faxén’s law
v2
i=Ltu1(r)|r=ri,
ω2
i=Lru1(r)|r=ri.(2.57)
Note that here the force Fiand the torque Tiare not considered, since we already
took them into account in equation (2.55). The fluid flow u0(r) + u1(r)is now
reflected on particle iand we have to consider the further contribution u2(r),
which has to satisfy the boundary condition u1(r) + u2(r) = v2
i+ω2
i×(r−ri)
on ∂Vi. Once, u2(r)is obtained, we can calculate it’s influence on particle jusing
Faxén’s theorem for Fj= 0 and Tj= 0. Continuing this procedure, we obtain
the series expansion of the exact flow field
u(r) = u0(r) + u1(r) + u2(r) + ··· .(2.58)
The corresponding expansion for the velocities of the particles are given by
vi=v0
i+v2
i+··· and vj=v1
i+v3
i+··· ,
ωi=ω0
i+ω2
i+··· and ωj=ω1
i+ω3
i+··· .(2.59)
23
2. HYDRODYNAMICS AT LOW REYNOLDS NUMBER
Note that, strictly speaking, the above representation of the hydrodynamic in-
teractions are restricted to two particles only. For more than two particles, we
additionally have to take into account that the flow field reflected by particle j
does not only influences the motion of the original particle i, but also the mo-
tion of all other particles. Nevertheless, these results can be used to describe
suspension, where the simultaneous hydrodynamic interaction of three or more
particles is improbable as compared to pair interactions, which is the case for
dilute solutions [24].
The first iteration step (2.56) yields expressions for the particles velocity up
to the order of (a/rij)3[24], which is often referred to as the Rotne-Prager ap-
proxiamtion1. At this level, only pair interactions play a role. Here, the induced
flow due to the motion of particle icaused by external forces and torques acting
on it is calculated according to equations (2.35) and (2.40). The influence of this
flow field on the motion of the other spheres is then simply determined using
Faxén’s theorem (see section 2.3.3). Therefore, we can write down the N-particle
velocities in terms of mobilities
v1
.
.
.
vN
ω1
.
.
.
ωN
=
µtt
11 ··· µtt
1Nµtr
11 ··· µtr
1N
.
.
.....
.
..
.
.....
.
.
µtt
N1··· µtt
NN µtr
N1··· µtr
NN
µrt
11 ··· µrt
1Nµrr
11 ··· µrr
1N
.
.
.....
.
..
.
.....
.
.
µrt
N1··· µrt
NN µrr
N1··· µrr
NN
F1
.
.
.
FN
T1
.
.
.
TN
(2.60)
with
µtt
ii =µt1,
µrr
ii =µr1,
µrt
ii =µtr
ii =0,
µtt
ij =Lt
¯
r0Lt
¯
rTO(¯
r0−¯
r)¯
r0=ri,¯
r=rj
,
µrt
ij = (µtr
ji)T=Lr
¯
r0Lt
¯
rTO(¯
r0−¯
r)¯
r0=ri,¯
r=rj
.(2.61)
An alternative approach to obtain expressions for the many particle mobility
matrix µis by a multipole expansion, which we will discuss in section 2.5.3.
1Usually only the translational part µtt is referred to as the Rotne-Prager tensor [71, 72]
(also known as Yamakawa tensor [73]). However, in the following we mean the mobility tensor
µas a whole (equation (2.4.2)), when we consider the Rotne-Prager approximation.
24
2.5 Hydrodynamic interactions between two parallel planar walls in Poiseuille
flow
Figure 2.4: Two infinity extended parallel plates with the distance 2Ware lo-
calized at z=−Wand z=W. The coordinate system is chosen such that an
arbitrary force acting on the fluid can be decomposed in components parallel and
perpendicular to the walls: F=Fk+F⊥=Fyey+Fzez.
2.5 Hydrodynamic interactions between two parallel planar
walls in Poiseuille flow
So far we described the hydrodynamic interactions of colloidal particles suspended
in an infinite and otherwise quiescent fluid. But in general, the situation differs.
Biological and experimental systems are often confined between boundaries and
are driven out of equilibrium by an imposed flow. Examples are single actin
filaments confined by the dense cytoskeleton network inside cells [5, 74], red blood
cells in vessels [75, 76] or colloidal particles in microchannels [77, 78, 79].
In this section we describe the hydrodynamic interactions for a system con-
fined between two parallel planar walls under the influence of an imposed Poiseuille
flow. We first briefly derive the corresponding Green’s function for this situation
and then describe the interactions on the level of point-particles and for extended
particles.
2.5.1 The two-wall Green tensor
Let us consider a fluid which is confined between two infinity extended parallel
plates with the distance 2W(see figure 2.4). The plates are laying in the x-y
plane and the origin is located at the center between both plates, so that the
walls are localized at z=−Wand z=W. Consider now a point force acting on
the fluid between the two plates at the position r0. According to equation (2.21)
25
2. HYDRODYNAMICS AT LOW REYNOLDS NUMBER
the flow field around this point force is given by
u(r) = T2W(r,r0)F.(2.62)
Here, the superset 2W denotes the Green’s tensor of the Stokes equations (2.19)
for the bounded geometry. To determine T2W(r,r0), we decompose the point
force in components parallel and perpendicular to the walls F=Fk+F⊥and
choose the coordinate system such that Fk=Fyey(see figure 2.4). Here, eyis
the unit vector along the y-axis. So, we can write down the Stokes equations for
the parallel component
η∇2uk(r)−∇pk(r) = −Fyeyδ(r−r0),
∇·uk= 0.(2.63)
Compared to an unbounded fluid (see section 2.2.2), the presence of the walls
violates the translational symmetry along the zdirection. Therefore, we expect
the two-wall Green’s tensor to depend only on the differences x−x0and y−y0.
Introducing the two dimensional space vector s= (x−x0)ex+ (y−y0)eyand the
wave vector q=qxex+qyey, we perform a Fourier transform in the in the x-y
plane [55]
T2W(s, z, z0) = Zeiq·sˆ
T2W(q, z, z0)dq2,
uk(s, z) = Zeiq·sˆ
uk(q, z)dq2,
pk(s, z) = Zeiq·sˆpk(q, z)dq2,
δ(r−r0) = 1
(2π)2δ(z−z0)Zeiq·sdq2.(2.64)
Inserting these terms into the Stokes equations (2.63), we obtain a coupled system
of differential equations for the parallel component [55]
η ∂2
∂z2−q2!ˆukx(q, z)−iqxˆpk(q, z) = 0,
η ∂2
∂z2−q2!ˆuky(q, z)−iqyˆpk(q, z) = −Fy
(2π)2δ(z−z0),
η ∂2
∂z2−q2!ˆukz(q, z)−∂
∂z ˆpk(q, z) = 0,
iqxˆukx(q, z) + iqyˆuky(q, z) + ∂
∂z ˆukz(q, z) = 0,(2.65)
where q2=q·q. Note that the last equation represents the incompressibility con-
dition. The stick boundary condition has to be satisfied at both walls such that all
26
2.5 Hydrodynamic interactions between two parallel planar walls in Poiseuille
flow
velocity components ˆ
uk(s, z)have to vanish at z=±W. From the incompress-
ibility condition we see that in addition ∂zˆukz(s, z)must also vanish at z=±W.
Together with these conditions equation (2.65) can be solved completely [55].
Once ˆ
ukzis determined, we can directly read off the components ˆ
T2W
xy ,ˆ
T2W
yy and
ˆ
T2W
zy from equation (2.62). Using reciprocity ˆ
T2W
ij (r1,r2) = ˆ
T2W
ji (r2,r1)[55, 64]
and rotational symmetry about the z-axis, we additional obtain the components
ˆ
T2W
yx ,ˆ
T2W
xx ,ˆ
T2W
zx ,ˆ
T2W
zx and ˆ
T2W
yz . To obtain ˆ
T2W
zz the whole procedure has to be
repeated for the perpendicular force component (F⊥=Fzez) [55]. Since all com-
ponents are now known, we can write the two-wall Green tensor in Fourier space
as [55]
ˆ
T2W(q, z, z0) = 1
16π2ηq (tnn(q, z, z0)ezez+itnp(q, z, z0)ezˆ
q+itpn(q, z, z0)ˆ
qez
−tpp(q, z, z0)ˆ
qˆ
q+rpp(q, z, z0)(1−ezez)) ,(2.66)
with ˆ
q=q/q and the scalar functions tnn, . . . , rpp which are explicitly given in
Ref. [55] and Appendix A. Similar to the Blake tensor1[80], this solution splits
up naturally into the Oseen tensor and a reflectional part
T2W(r,r0) = TO(r−r0) + TR(r,r0).(2.67)
Using the polar coordinates qx=qcos φand qy=qsin φ, we transform the two-
wall Green tensor back into real space. So, we finally obtain for the reflectional
part2[55]
TR
xx =−1
16πη
∞
Z
0 J0(qs) + (y−y0)2−(x−x0)2
s2J2(qs)!t1pp(q, z, z0)dq,
+1
8πη
∞
Z
0
J0(qs)r1pp(q, z, z0)dq,
TR
zz =1
8πη
∞
Z
0
J0(qs)t1nn(q, z, z0)dq,
1The Blake tensor is Green’s function of the Stokes equations for the fluid flow in the
presence of one bounding wall
2Note, expressions for the two-wall Green tensor were first derived by Liron and Mochon
in 1976 [81]. In analogy to electrostatics, they set up an infinite series of images in order
to obey the no-slip boundary condition at the two bounding walls. Since the above integral
representation derived by Jones [55] is more symmetric and easier to handle than the original
form, we use this representation.
27
2. HYDRODYNAMICS AT LOW REYNOLDS NUMBER
Figure 2.5: Flow field around a point force acting a) perpendicular and b) parallel
to the walls. The yand z-axis are rescaled by the distance from the channel center
to the walls W.
TR
xy =1
8πη
(x−x0)(y−y0)
s2
∞
Z
0
J2(qs)t1pp(q, z, z0)dq,
TR
xz =1
8πη
x−x0
s
∞
Z
0
J1(qs)t1pn(q, z, z0)dq. (2.68)
Here, Jndenotes the Bessel function and s=q(x−x0)2+ (y−y0)2the distance
in the x−yplane. Due to the rotational symmetry abot the z-axis the component
TR
yy is the same as TR
xx with the change (x−x0)→(y−y0)and (y−y0)→(x−x0),
TR
yz is the same as TR
xz with the change (x−x0)→(y−y0). The other components
follow from the reciprocal theorem [55, 64]
T2W
ij (r,r0) = T2W
ji (r0,r).(2.69)
Figure 2.5 shows the flow field around a point force acting on the fluid between
the bounding walls.
28
2.5 Hydrodynamic interactions between two parallel planar walls in Poiseuille
flow
2.5.2 Hydrodynamic interactions of point-like particles
Let us consider now Nspherical particles suspended in a fluid which is confined
between two parallel planar walls and undergoing a pressure driven Poiseuille flow
u0(r) = v0"1−z
W2#ey,(2.70)
where v0is the maximum velocity at the centerline. As we already discussed in
section 2.3.3, the resulting flow field is given by the superposition of the imposed
flow u0(r)in the absence of the particles and the disturbance flow due to the
present of the spheres
u(r) = u0(r) +
N
X
j=1 I
∂Vj
T2W(r,r0)fj(r0)dr02.(2.71)
Since we are here interested only in the hydrodynamic interactions of point-like
particles, we neglect again the rotational motion. So, the stick boundary condition
on each particle reads
vi=u0(r) + I
∂Vi
TO(r,r0)fi(r0)dr02+I
∂Vi
TR(r,r0)fi(r0)dr02
+
N
X
j6=iI
∂Vj
T2W(r,r0)fj(r0)dr02,r∈∂Vi.(2.72)
Here, we explicitly split the two-wall Green tensor into the Oseen and reflectional
part (2.67). Multiplying this equation by 1/(4πa2)and integrating both sides
over the surface of particle i, we obtain
vi=Ltu0(r) + µtFi+1
4πa2I
∂ViI
∂Vi
TR(r,r0)fi(r0)dr02dr2
+1
4πa2I
∂Vi
N
X
j6=iI
∂Vj
T2W(r,r0)fj(r0)dr02dr2,r∈∂Vi,(2.73)
where we used Faxén’s formula for the translational motion of a sphere (2.46).
Assuming now point particles (a→0), so that the diameter of the spheres is
much smaller than the particles separation, we can approximate TR(r,r0)for
r,r0∈∂Viand T2W(r,r0)for r∈∂Vi,r0∈∂Vjby TR(ri,ri)and T2W(ri,rj),
respectively. So, we finally obtain
vi=u0(ri) + hµt1+TR(ri,ri)iFi+
N
X
j6=i
T2W(ri,rj)Fj.(2.74)
29
2. HYDRODYNAMICS AT LOW REYNOLDS NUMBER
We can write down this equation in terms of mobilities µtt
ij, where the self- and
cross mobilities for this situation are given by
µtt
ii =µt1+TR(ri,ri),
µtt
ij =T2W(ri,rj).(2.75)
Therefore, the many particle interaction is completely described by the mobility
tensor
v1
.
.
.
vN
=
u0(r1)
.
.
.
u0(rN)
+
µtt
11 ··· µtt
1N
.
.
.....
.
.
µtt
N1··· µtt
NN
F1
.
.
.
FN
.(2.76)
Note that this expression for the hydrodynamic interactions are only valid for
point like particles, which means in this case that not only the separation of
the Nparticles has to be much larger than the particle’s diameter but also the
channel width d= 2Whas to be much larger. The reflectional part of the two-
wall Green tensor can be interpreted as an infinite set of mirror images on both
walls supplemented by other singularities of the stokes equation to fulfill the stick
boundary condition [81]. Therefore, the distance to these mirror particles has to
be larger than the particle’s size.
In order to describe the hydrodynamic interactions for finite particle sizes,
we could use the methods of reflections (see section 2.4.2) and replace the Oseen
tensor by the two-wall Green tensor. But since the flow past a moving sphere
cannot be exactly calculated for this geometry and the imposed Poiseuille flow has
to be taken into account in this framework, we are using an alternative method
to obtain expressions for the mobility matrix.
2.5.3 Multipole expansion
An alternative way to calculate the mobility tensor for spherical particles is given
by a multipole expansion [55, 82, 83, 84, 85]. Let us consider again Nspherical
particles confined between two parallel planar walls under the influence of a pres-
sure driven Poiseuille flow u0(r). The presence of the particles disturb this flow,
such that the resulting flow field u(r)is given by1
u(r) = u0(r) +
N
X
j=1 ZT2W(r,r0)fj(r0)d3r0.(2.77)
1Note that here fi(r)is the force density per unit volume that the particle iexerts on the
fluid.
30
2.5 Hydrodynamic interactions between two parallel planar walls in Poiseuille
flow
Again, we assume stick boundary conditions, such that the flow velocity u(r)is
equal to the surface velocity of particle i
vi+ωi×(r−ri) = u0(r) +
N
X
j=1 ZT2W(r,r0)fj(r0)d3r0,r∈∂Vi.(2.78)
We project this integral equation on an appropriate set of complex basis functions
to obtain a system of linear algebraic equations, which connects velocity to force
multipoles [55, 56]. The basis functions are solutions of the homogeneous Stokes
equations (equation (2.19) with f= 0) [55, 56, 86], for which two sets exist:
v+
lmσ(r), which is regular at r= 0 and grow at infinity, and v−
lmσ(r), which is
singular at r=0and vanish at infinity. They are equivalent to those due to
Lamb [87] and explicitly given in Ref. [56]. The indices l,mand σtake the
respective values l= 1,2,3, . . . ,m=−l, −l+ 1, . . . , l −1, l and σ= 0,1,2. In
terms of these functions we expand the Oseen tensor (2.25) about an arbitrary
point R[88]
TO(r−r0) = 1
ηX
lmσ
v−
lmσ(r>)⊗v+∗
lmσ(r<),(2.79)
where r>,r<refer, respectively, to the larger and smaller of the relative position
vectors r−Rand r0−R, and the superscript ∗denotes the complex conjugate.
The reflectional part TR(r,r0)can be expanded similarly, but using only the
regular functions [55]
TR(r,r0) = X
lmσ X
l0m0σ0
TR(lmσ, l0m0σ0)v+
lmσ(r−R)⊗v+∗
l0m0σ0(r0−R),(2.80)
where TR(lmσ, l0m0σ0)denotes the expansion coefficients. Expanding now the
difference field u(r)−u0(r)in terms of regular solutions centered on the arbitrary
point R, we obtain
u(r)−u0(r) = X
lmσ
c(lmσ)v+
lmσ(r−R).(2.81)
Here, c(lmσ)denotes the expansion coefficients for the difference field. In order
to calculate the expansion coefficients in equation (2.80) and equation (2.81) one
needs an adjoint set of functions w±
lmσ [56, 89] which are bi-orthonormal with
respect to v±
lmσ on the surface of a sphere of arbitrary radius ain the sense that
Dw±
lmσ(r)a−1δ(r−a)v±
l0m0σ0(r)E=δll0δmm0δσσ0.(2.82)
31
2. HYDRODYNAMICS AT LOW REYNOLDS NUMBER
Explicit expressions for this adjoint set of functions are given in Ref. [56] and in
Appendix A. In analogy to quantum mechanics, we have introduced the scalar
product hg(r)|h(r)ibetween two complex functions g(r)and h(r)as
hg(r)|h(r)i=Zg∗(r)·h(r)dr3.(2.83)
Therefore, we can calculate the expansion coefficients for the difference field
u(r)−u0(r)expanded around the center of particle i
c(ilmσ) = Dw+
lmσ(i)δiu(r)−u0(r)E
(2.77,2.78)
=Dw+
lmσ(i)δivi+ωi×(r−ri)−u0(r)E.(2.84)
where c(ilmσ)denotes the expansion coefficient expanded around the center of
particle i,w+
lmσ(i)is short for w+
lmσ(r−ri),δiis short for a−1δ(|r−ri|−a)and
adenotes the radius of the spheres. Using the right hand side of equation (2.78)
and inserting the two-wall Green tensor expanded around the center of particle
jyields the linear relation [56]
c(ilmσ) =
N
X
j=1 Dw+
lmσ(i)δiT2W(r,r0)fj(r0)E
=X
jl0m0σ0
M(ilmσ, jl0m0σ0)f(jl0m0σ0)(2.85)
between the coefficients c(ilmσ), also called the velocity multipole moments and
the force multipole moments f(jl0m0σ0), which are defined as
f(jl0m0σ0) = Dv+
l0m0σ0(j)|fjE.(2.86)
The infinite-dimensional matrix
M(ilmσ, jl0m0σ0) = Dw+
lmσ(i)δi|T2W|w+
l0m0σ0(j)δjE(2.87)
represents the two-wall Green function. One writes equation (2.85) in a symbolic
notation,
c=Mf.(2.88)
where c and f are infinite dimensional vectors with respective components c(ilmσ)
and f(ilmσ)and matrix M contains the elements M(ilmσ, jl0m0σ0), which is
composed of M =M0+TR
0+M1+TR
1. Here, M0,M1arise from the Oseen tensor
and TR
0,TR
1from the reflectional part, M1and TR
1are off diagonal in particle
labels (i6=j), while M0and TR
0are diagonal and contain the single particle
contribution (i=j). Inverting equation (2.88) gives
f=Zc,Z=M−1,(2.89)
32
2.5 Hydrodynamic interactions between two parallel planar walls in Poiseuille
flow
where Z denotes the friction matrix for this situation. Note that the equations
(2.88) and (2.89) are exact and compared to equation (2.4.2) the knowledge and
construction of the flow field u(r)is not needed. The complex quantities in
equation (2.88) and (2.89) can be transformed to a real basis, where for example
the velocities viand ωiof particle iare linear combinations of the first multipole
orders c(i1m0) and c(i1m1). Likewise the forces Fiand torques Ticorrespond to
the first order multipoles f(i1m0) and f(i1m1) [55]. So, we can rewrite equation
(2.89) in an equivalent real Cartesian representation (details are given in Ref.
[56])
F
T
E
Hk
.
.
.
=
ζtt ζtr ζtd ζth
k···
ζrt ζrr ζrd ζrh
k···
ζdt ζdr ζdd ζdh
k···
ζht ζhr ζhd ζhh
k···
.
.
..
.
..
.
..
.
....
v−u0
ω−ω0
−g0
−hk
0
.
.
.
.(2.90)
The column vector on the left-hand side contains the known force multipoles
acting on the particles. Here, F= (F1,...,FN)is a N-particle column vector
of forces; similarly, one has torques T, stresslets E, and higher force multipoles
Hk. The column vector on the right-hand side consists of the velocity multipoles
of vi−u0, where v−u0= (v1−u0(r1),...,vN−u0(rN)) contains the relative
translational velocities and ω−ω0the relative angular velocities with the vorticity
ω0=1
2∇×u0of the imposed flow. Since the motion of rigid bodies is completely
described by Uand Ω, one only has the strain rate tensor g0=1
2(∇ ⊗ u0+
(∇⊗u0)T)and higher multipoles of the imposed flow. In case of the Poiseuille
flow, the third-rank Cartesian tensor h0αβγ =1
2∂β∂γv0α(r)decomposed in its
three irreducible moments hk
0(k= 1,2,3) is also relevant [55]. Finally, the grand
resistance matrix contains the N×Nmatrices ζab, the elements of which are
Cartesian tensors ζab
ij .
Typically, one controls forces and torques acting on particles as well as the im-
posed flow, while translational and rotational particle velocities are unknown. So,
a partial inversion of the grand resistance matrix introduces the grand mobility
matrix [64, 90],
v−u0
ω−ω0
−E
−Hk
.
.
.
=
µtt µtr µtd µth
k···
µrt µrr µrd µrh
k···
µdt µdr µdd µdh
k···
µht µhr µhd µhh
k···
.
.
..
.
..
.
..
.
....
F
T
g0
hk
0
.
.
.
,(2.91)
with
µtt µtr
µrt µrr != ζtt ζtr
ζrt ζrr !−1
(2.92)
33
2. HYDRODYNAMICS AT LOW REYNOLDS NUMBER
and the remaining mobilities are related to the resistant matrices ζab. We use
them to explicitely give the translational and rotational velocities of Nspherical
colloids subject to external forces, torques, and an imposed Poiseuille flow,
v=u0+µttF+µtrT+ (µttζtd +µtrζrd)g0
+
3
X
k=1
(µttζth
k+µtrζrh
k)hk
0,
ω=ω0+µrtF+µrrT+ (µrtζtd +µrrζrd)g0
+
3
X
k=1
(µrtζth
k+µrr
ζrh
k)hk
0.
(2.93)
2.5.4 L-multipole approximation and Lubrication correction
In principle, all the elements of the matrices M and Z can be calculated, and
therefore the exact Stokes velocities of the spheres can be obtained. However, in
practice the infinite-dimensional equations (2.88) and (2.89) need to be approxi-
mated to use them in simulations. In the L-multipole approximation the matrices
M(ilmσ, j0l0m0σ0),Z(ilmσ, j0l0m0σ0)and vectors c(ilmσ),f(ilmσ)are truncated
such that only elements with l, l0≤Lare taken into account [56, 85]. This defines
approximations for the involved matrices and vectors
cL=MLfL,fL=ZLcL,(2.94)
and gives rise to approximated grand resistance and mobility matrices, ζLand
µL, respectively. As long as the particles are sufficiently far from each other, the
L-multipole approximation converges rapidly to the exact mobility matrix with
increasing L[56]. However, for nearby particles the convergence is very slow and
an alternative method is needed. To account for the lubrication regime, we apply
a lubrication correction on the level of the friction matrix [56, 85, 91].
Consider two identical spherical particles in an unbounded fluid at positions
r1and r2, whose distance is very small. The lubrication forces between the
spheres are then very strong. For this situation the exact two particle friction
matrix ζ(1,2) can be evaluated using the Jeffrey-Onishi asymptotic formulas
supplemented by a series expansion in the inverse particle distance [92, 93]. Here,
ζ(1,2) =
ζtt(1,2) ζtr(1,2) ···
ζrt(1,2) ζrr(1,2) ···
.
.
..
.
....
(2.95)
denotes the exact two particle friction matrix between particle 1 and 2. In order
to describe the lubrication regime for a system of Nspheres, we introduce the
34
2.5 Hydrodynamic interactions between two parallel planar walls in Poiseuille
flow
lubrication correction in a pairwise additive correction [94]. For each element
of the N-particle friction matrix ζLwe know the exact two-sphere resistance
interaction of the particles iand jalone in an unbounded fluid, given by the
two-particle friction matrix ζ(i, j), which we simply add to ζL. This defines a
N-particle resistance matrix
˜
ζL=ζL+ζpair,(2.96)
where ζpair
Ldenotes the superposition of the exact two-particle friction matrix
ζ(i, j). However, part of the two-sphere resistance interactions, the far-field part,
has already been included by the friction matrix ζL. To avoid double counting of
the far-field part, we subtract the two-body interactions already contained in ζL.
To do so, we subtract from each element of ˜
ζLthe two-particle friction matrix in
unbounded fluid ζL(i, j)which one determines in the L-multipole approximation:
¯
ζL=ζL+ ∆pair
L.(2.97)
Here, the remainder is given by ∆pair
L=ζpair −ζpair
L, where ζpair
Ldenotes the
superposition of the approximated two-body friction matrices ζL(i, j). Note that
for L→ ∞,∆pair
Lbecomes zero and ζLconverges to the exact friction matrix.
The standard lubrication correction (2.97) affects not only the relative, but
also the the collective motions. In particular, it changes the total force and the
total torque exerted by a pair of particles on the fluid [91]. As a result, equation
(2.97) gives the incorrect large-distance asymptotic behavior for more than two
particles, for example, when nearby particles move together and interact with a
distant particle [91]. Since the lubrication phenomena arises due to the relative
motion of the particles only, one has to project out the forces due to the collective
motion. It was shown, that it is sufficient to project out the collective translational
motion only [91] such that we have to replace the two-particle friction matrices
by
s(i, j) = qTζ(i, j)q,sL(i, j) = qTζL(i, j)q,(2.98)
where projector qis given by [91]
qtt
ij =1
2(−1)i+j1,qa=b
ij =δij1,qa6=b
ij =0.(2.99)
In full analogy to the friction matrix ζpair composed of the exact two-body
friction matrices, we construct a matrix ζwall from the exact one-particle friction
matrices in the presence of the lower ζlow(i)and upper wall ζup(i). They can
be calculated from the Jeffrey-Onishi asymptotic formulas supplemented by a
series expansion in the inverse particle distance in the limit that one particle is
fixed in space and becomes infinitely large [95]. Finally, the complete L-multipole
35
2. HYDRODYNAMICS AT LOW REYNOLDS NUMBER
approximation for the friction matrix corrected by lubrication forces between two
particles as well as between a particle and the bounding walls, is given by
¯
ζL=ζL+ ∆pair
L+ ∆wall
L.(2.100)
In analogy to above, we have introduced ∆wall
L=ζwall −ζwall
L, where ζwall
Lis com-
posed from the L-multipole approximation for the one-particle friction matrix in
the presence of the lower and the upper wall ζlow
L(i) + ζup
L(i). Once ¯
ζLin equa-
tion (2.100) is calculated, one determines the lubrication-corrected L-multipole
approximation for the grand mobility matrix ¯
µLby partial inversion.
2.6 Stokesian dynamics simulations
In the previous sections we derived the equations of motions for Nspherical
particles suspended in a Newtonian fluid, where in the most general case we
considered, the particles are confined between two planar parallel walls and are
subjected to external forces, torques and an imposed Poiseuille flow (see section
2.5.3). As we showed, the particle velocities are linearly related to these forces,
torques and the imposed flow via the grand mobility matrix (2.91). Since the
dependence of the grand mobility matrix on the particle coordinates is highly
nonlinear, we describe in this section the integration scheme for the numerical
solution of the equations of motion.
2.6.1 Equations of motion
According to equation (2.91), the translational and angular velocities of the
spheres are given by
v=u0+¯
µtt
LF+¯
µtr
LT+¯
µtd
Lg0+
3
X
k=1
¯
µth
L,khk
0,
ω=ω0+¯
µrt
LF+¯
µrr
LT+¯
µrd
Lg0+
3
X
k=1
¯
µrh
L,khk
0.(2.101)
Note that, in the case of point-particles, the rotational motion is neglected and
equation (2.101) reduces to (2.76). This evolution equation that we use in the
Stokesian-dynamics simulations can then be written in compact differential form
as
dx=vdt, v=u0(x) + m(x)h(x) + µ(x)F(x).(2.102)
36
2.6 Stokesian dynamics simulations
Here, we have introduced the 6N-dimensional differential dx= (dx,dφ) =
(dr1,...,drN,dφ1,...,dφN), where the vector dφidefines an infinitesimal ro-
tation of particle i, the generalized N-particle velocities v= (v,ω)and u0=
(u0,ω0), the generalized forces F= (F,T), the higher order velocity multipole
moments h= (g0,h1
0,h2
0,h3
0), the mobility matrix
µ= ¯
µtt
L¯
µtr
L
¯
µrt
L¯
µrr
L!,(2.103)
and the matrix mcontaining the other elements of the grand mobility matrix
m= ¯
µtd
L¯
µth
L,1¯
µth
L,2¯
µth
L,3
¯
µrd
L¯
µrh
L,1¯
µrh
L,2¯
µrh
L,3!.(2.104)
As we already mentioned in the introduction to this section, the matrices mand
µare highly nonlinear in the particle coordinates. Therefore, we use a higher-
order integration scheme for the numerical solution of this equation, which is
presented in the following.
2.6.2 Numerical integration scheme
The simplest method to obtain an numerical solution of equation (2.102) is the
well known Euler method [96]. Here, the equation of motion is discretized with
respect to time t
∆xt=v(xt)∆t, (2.105)
and then the configuration update xt+∆t=xt+ ∆xtis calculated. However,
this procedure is not recommended for practical use, because it is less accurate
in comparison to other methods and it is not very stable [97]. Therefore, one
should use integration schemes of higher order in the time step ∆tto avoid the
accumulation of numerical errors. Higher order algorithm typically combine the
informations of several Euler-steps such that they match the Taylor expansion
of this function up to the order of ∆tn. By far the most often used integration
scheme is the so called fourth-order Runge-Kutta algorithm (also referred to as the
classical Runge-Kutta method) [98, 99]. This algorithm consists of the following
intermediate Euler steps
∆x1=v(xt),
∆x2=v(xt+1
2x1),
∆x3=v(xt+1
2x2),
∆x4=v(xt+x3).(2.106)
37
3
Brownian motion
«Je planmäßiger die Menschen
vorgehen, desto wirksamer vermag
sie der Zufall zu treffen.»
Friedrich Dürrenmatt
39
3. BROWNIAN MOTION
In the previous chapter we derived the equations of motion for Nspherical
objects suspended in a Newtonian fluid, where the hydrodynamic interactions are
completely described by the mobility matrix. The components of this matrix are
highly nonlinear in the particle coordinates and can be calculated for different
geometries (compare unbounded fluid [see section 2.4] with bounded fluid [see
section 2.5]). Next to this deterministic motion (Stokesian dynamics) the particles
also undergo erratic movements, the well known Brownian motion. The origin of
this motion arise due to collisions with the surrounding fluid molecules [24, 100].
On the molecular level the dynamics of the collisions are completely deter-
ministic. However, even if we would know all initial conditions, it is impossible
to predict the future of the system, since the computational effort is exorbitant
(around 1022 water molecules per cm3). Moreover, chaos theory shows that even
if we had all informations on the interactions and thus the complete equations of
motions, any slight error in the initial conditions would be magnified exponen-
tially [101]. So, in fact, there is practically no predictability.
Instead of including all degrees of freedom of the water molecules, we describe
the influence of the collision between them and the suspended objects by a sta-
tistical theory. In this description the fluid is first considered to be continuous,
which corresponds to averaging over a huge number of collisions of the suspended
particles with the fluid molecules. This averaging leads to the deterministic hy-
drodynamic forces and torques acting on the suspended particles as described in
detail in the previous chapter. Secondly, the solvent can be treated as a source
of noise that causes a random motion due to the continuous collisions with the
suspended particles.
This picture can be formulated by two different but equivalent concepts of the-
oretical mechanics. In the first approach the deterministic equations of motion are
supplemented by random forces and torques acting on the particles. This results
in differential stochastic equation (SDE), which is referred to as the Langevin
equation. In the second approach, one is not interested in the single particle
trajectories but instead, in ensemble averages. These can be calculated by the
probability distribution of the system in phase space, where the time evolution
is governed by a diffusion-type differential equation, the so called Smoluchowski
equation.
40
3.1 The Langevin equation
3.1 The Langevin equation
As discussed in section 2.6.1, the equations of motion for Nspherical particles
suspended in a solvent can be expressed in the following compact equation
v=u0(x) + m(x)h(x) + µ(x)F(x),(3.1)
where v= (v,ω)and u0= (u0,ω0)are the generalized N-particle velocities,
F= (F,T)the generalized forces, h= (g0,h1
0,h2
0,h3
0)the higher order velocity
multipole moments, µthe mobility matrix and the matrix mcontains the other
elements of the grand mobility matrix. Note that in the absence of the imposed
flow u0, the first two terms on the right hand side of this equation vanish and
the components of the mobility matrix µdepend on the chosen system geometry
(see section 2.4 and section 2.5). In addition to this deterministic motion, every
particle experiences a rapidly fluctuating random force and torque Γ(t), which
arises from collisions with the surrounding fluid particles [24, 100]. These forces
fluctuate on a timescale of about 10−14s, which is the typical relaxation time of
the surrounding fluid. Considering these stochastic forces the equation of motion
reads
v=u0(x) + m(x)h(x) + µ(x)(F(x) + Γ(t)).(3.2)
This equation is a first order stochastic differential equation, the so called Langevin
equation. In 1908, Paul Langevin was able to describe the diffusional process of
the Brownian motion with the help of such an equation [102]. In general the
Langevin equation reads
d
dtx(t) = h(x(t), t) + gΓ(t),(3.3)
where the last term on the right hand side is the so called noise term. If the
quantity gis constant then the last term is referred to as additive noise. If
g=g(x(t)) is a function of the variable x(t)then the last term is referred to
as multiplicative noise, which leads to the Ito-Stratonovich dilemma (see section
3.1.2).
Since the stochastic force Γ(t)is a random number, the rules of ordinary
differential equations can no longer be applied to equation 3.3. Therefore this
equation is meaningless in a strict mathematical sense. With the help of the
stochastic calculus of Wiener process we will be able to treat the trajectories of
this random motion analytically. But before, we want to discuss the properties
of the random force Γ(t).
41
3. BROWNIAN MOTION
3.1.1 Fluctuation-dissipation theorem
In order to treat equation (3.2) analytically, we need more informations about
the random force Γ(t). For simplicity, we consider one particle in one dimension
with mass mand friction coefficient ζ
m˙v+ζv = Γ(t).(3.4)
Since the systematic interactions with the solvent molecules are made explicit,
all interactions which will not vanish in the ensemble average belong to the de-
terministic part. Therefore, the ensemble average of the fluctuating force is zero
hΓ(t)i= 0.(3.5)
Here, the brackets h. . . idenote the ensemble average. As mentioned above, the
random forces fluctuate on a timescale of about 10−14s, whereas the relaxation
time of the suspended particle is in the order of a few nanoseconds. Due to a
large separation in the time scales, it is sufficient to assume a delta correlated
random force in time [24]
hΓ(t)Γ(t0)i=qδ(t−t0).(3.6)
Here, δdenotes the delta distribution and qis a constant which is referred to as
fluctuation strength and has to be determined. In equilibrium the equipartition
theorem must hold, which implies (in one dimension)
1
2mDv2E=1
2kBT, (3.7)
where kBis the Boltzman constant and Tthe absolute temperature. Together
with equations (3.4),(3.5) and (3.6) we find [24]
1
2mDv2E=1
2mv2(0)e−2t
τr+1
2m
t
Z
0
t
Z
0
e−(2t−t0−t00 )
τrhΓ(t0)Γ(t00)idt0dt00
tτr
=1
4
τr
mq, (3.8)
where τr=m/ζ is the relaxation time of the suspended particle. If we compare
this result with the equipartition theorem (3.7), we find q= 2kBTζ. This leads
to
hΓ(t)Γ(t0)i= 2kBTζδ(t−t0).(3.9)
This relation is referred to as the fluctuation dissipation theorem. It connects the
fluctuation strength with the friction coefficient, which determines the dissipation
42
3.1 The Langevin equation
t
w(t)
Figure 3.1: Three simulated sample paths of the one-dimensional Wiener process,
illustrating their great variability.
into heat. In three dimensions and for Nparticles equations (3.5) and (3.9) read
[103, 104]
hΓ(t)i=0
hΓ(t)⊗Γ(t0)i= 2kBTζδ(t−t0),(3.10)
where ζ=µ−1is the 6N×6Nfriction matrix (2.92).
3.1.2 The Wiener process
In the previous section we discussed the properties of the random force Γ(t). In
this section we want to evaluate the Langevin equation (3.2). Therefore we first
rewrite this equation as
v=f(x) + µ(x)Γ(t),(3.11)
where f(x) = u0(x) + m(x)h(x) + µ(x)F(x)). This equation can be formally
solved by integration, which leads to
x(t)−x(0) =
t
Z
0
f(x(t0)) dt0+
t
Z
0
µ(x(t0))Γ(t0)dt0.(3.12)
43
3. BROWNIAN MOTION
The second integral on the right hand side of equation (3.12) cannot be solved in
the sense of a Riemann sum. In order to analyze this term we need the stochastic
calculus of Wiener processes [105, 106, 107, 108]. The integral equation for the
increment ∆x(t) = x(t+ ∆t)−x(t)is given by
∆x(t) =
t+∆t
Zt
f(x(t0)) dt0+
t+∆t
Zt
µ(x(t0))dw(t0),(3.13)
where dw(t) = w(t+dt)−w(t)is the increment of a Wiener process
w(t) =
t
Z
0
Γ(t0)dt0,(3.14)
which is a Gaussian variable with zero mean hw(t)i= 0 and time correlation
hw(t)⊗w(t0)i= 2kBTζmin(t, t0).(3.15)
The increments dw(t)are then Gaussian variables too, with the following prop-
erties
hdw(t)i= 0,
hdw(t)⊗dw(t)i= 2kBTζdt, (3.16)
which one can easily validate by the definition (3.14) and equation (3.10). The
increments at different times are independent of each other. Also, one easily sees
from equation (3.16) that dw(t)is proportional to √dt.
The second integral on the right hand side of equation (3.13) has the form of
aStieltjes-integral. This integral is defined on the interval [t, t + ∆t]subdivided
into the equidistant parts t=t0< t1<, ··· <, tn=t+ ∆tas the limes of the
partial sums [109]
t+∆t
Zt
µ(x(t0))dw(t0) = ms-lim
n→∞
n
X
i=1
µ(x(τi))[w(ti)−w(ti−1)],(3.17)
where τi∈[ti−1, ti)and ms-lim denotes the mean-square limit1. It is quite easy to
show, that in contrast to ordinary Riemann integrals this integral (3.17) depends
on the particular choice of the intermediate point τi[107, 110], which we write as
τi= (1 −α)ti−1+αti.(3.18)
1Let Xnbe a sequence of random variables and let Xbe a random variable. The sequence
Xnis said to converge to Xin mean-square (ms-limn→∞Xn=X) if limn→∞ (Xn−X)2= 0.
44
3.1 The Langevin equation
Here, αcan take any value between zero and one. There are mainly two different
types of interpretations of the integral (3.17):
α= 0 Itô interpretation,
α=1
2Stratonovich interpretation.(3.19)
The difference between these interpretations can be clearly seen: Itô evaluates
µ(x(τi)) at the beginning of the interval, whereas Stratonovich evaluates this
term at the mid-point. If we consider ∆tto be very small, ∆t→dt, we can write
equation (3.13) as the following Itô stochastic differential equation
dx=fdt+µdw,(3.20)
where µ(x(t)) is evaluated at the beginning of the interval [t, t+dt]. Respectively,
we can write down the following Stratonovich stochastic differential equation
dx=fdt+µ◦dw,(3.21)
where the symbol "◦" denotes that µ(x(t+1
2dt)) is evaluated at the intermediate
point t+1
2dt. Now, the question arises, which one of these two interpretations
one has to choose. In principle, the answer depends on the chosen system. For
physically systems one usually tends to choose the Stratonovich interpretation.
This is the natural choice for an interpretation which assumes Γis real noise in
equation (3.10) (not a white noise) with finite correlation time, which is then
allowed to become infinitesimally small after calculating measurable quantities
[107]. Whereas in a system which is intrinsically discontinuous, e.g., in the stock
exchange or in the evolution of biological populations the Itô interpretation is
usually applied [109]. Another advantage of the Stratonovich interpretation is
that it allows the use of the ordinary rules of differential calculus, whereas Itôs
formula shows that rule of changing variables is not given by ordinary calculus
[107, 110]. On the other hand, Stratonovich stochastic differential equations are
usually translated to equivalent Itô stochastic differential equations for mathe-
matical analysis of their properties. To do so, we expand µ(x)in terms of tin
equation (3.21)
dx=fdt+µ(x+1
2dx+. . . )dw,(3.22)
with x(t+1
2dt) = x(t) + 1
2dx(t) + . . . . Considering only terms up to the order
of dtyields1[111, 112]
dx=fdt+µ+1
2(dx·∇)µ+O(dt)dw
=fdt+1
2([(dw⊗dw)µ]∇)µ+µdw+O(dt)3
2.(3.23)
1Note that here ∇= (∇1,...,∇N)denotes the generalized 3N-dimensional Nabla operator.
45
3. BROWNIAN MOTION
Note that dwis proportional to √t (3.16) and that dxalso consists of a term
proportional to dw(3.22). To eliminate the stochastic variable dwwe take the
mean of equation (3.23). Together with the properties of the Wiener increment
(3.16), we obtain
hdxi= [f+kBT(µζ)∇µ]dt
= [f+∇D]dt,
hdx⊗dxi= 2kBT(µζ)µ
= 2Ddt. (3.24)
Here, we used equation (2.92) and the generalized Einstein relation D=kBTµ,
where Dis the diffusivity [113]. So, finally we can write Stratonovichs interpre-
tation as an equivalent Ito stochastic differential equation [107, 110]
dx= [f+∇D]dt+µdw.(3.25)
The additional term ∇Din equation (3.25) compared to equation (3.20) is re-
ferred to as spurious drift [105], which can be interpreted as an additional drift
velocity. Since this term is proportional to the temperature T, this effect is noise
induced.
3.2 The Smoluchowski equation
Instead of calculating the trajectories of a stochastic system explicitly, one can
also describe the stochastic dynamics with the help of a probability density
ψ(x, t), where ψ(x, t)dxis the probability to find the system in [x,x+dx]at time
t. The evolution of this probability density is governed by the so called Fokker-
Planck equation. It can be shown that one can assign a unique Fokker-Planck
equation to each Langevin equation [110]. Therefore, different Langevin equa-
tions that possess the same Fokker-Planck equation are stochastically equivalent.
The Fokker-Planck equation can be obtained by the Kramer-Moayl expansion
[110]
∂
∂tψ(x, t) =
∞
X
n=1
3N
X
i1,...,in
(−1)n∂n
∂xi1···∂xin
D(n)
i1...in,(3.26)
where the Kramer-Moyal coefficients D(n)are given by the following moments
[110]
D(n)=1
n!lim
τ→0
1
τD[x(t+τ)−x(t)]⊗nE.(3.27)
46
3.3 Brownian dynamics simulation
Here, the subscript ⊗ndenotes the n-th dyadic product. In the over-damped
regime, we can read off the first two Kramer-Moyal coefficients directly from
equations (3.24) and (3.25) together with the properties of the Wiener increment
(3.16):
D(1) =f+∇D,
D(2) =kBTM=D.(3.28)
It can be shown that for the Langevin equations with Gaussian white noise (3.11),
higher Kramer-Moyal coefficients are identical zero [110]
D(n)=0, n > 2.(3.29)
Therefore, together with these coefficients (3.28) and (3.29), we obtain the fol-
lowing Fokker-Planck equation
∂
∂tψ(x, t) = −∇·[fψ−D∇ψ].(3.30)
Note, in the over-damped limit, the Fokker-Planck equation is applied to particle
position distributions ψ(x, t), which is better known as the Smoluchowski equation
[110]. We can write down the Smoluchowski equation (3.30) in the form of a
continuity equation
∂
∂tψ=−∇·j,
=−∇·(jdrift +jdiff),(3.31)
where jis probability density current, which consists of a deterministic part
jdet =fψand a diffusional part jdiff =−D∇ψ.
3.3 Brownian dynamics simulation
In section 3.1.2, we derived the stochastic differential equations (3.20) and (3.21).
With their help we are able to describe the fluctuating forces due to collisions
between the molecules and the suspended particles mathematically. In this sec-
tion we will describe how we solve this kind of equations numerically using the
method of Brownian dynamics simulation.
47
3. BROWNIAN MOTION
3.3.1 Numerical integration scheme
The simplest discretization scheme of the Langevin equation in the Itô form (3.25)
is the Euler step (see also section 2.6.2)
∆xt= [f(xt) + ∇D(xt)] ∆t+H(xt)∆¯
wt.(3.32)
Here, we have introduced the matrix H(x), which is the solution of the following
relation:
2D(x) = H(x)HT(x).(3.33)
This relation indicates, to say it with sloppy words, that the matrix H(x)is
somehow the square root of 2D(x). We can determine this quantity by Cholesky
decomposition, which we will describe in section 3.3.2. The first two moments of
the the Wiener increment ∆¯
ware given by
h∆¯
wti= 0,
h∆¯
wt⊗∆¯
wti=1∆t, (3.34)
which are independent of each other at different times. To obtain the correct
Wiener increments, we generate a new set of 6Nindependent random numbers
niat each time step, obeying the Gaussian distribution
p(n) = 1
2πe−n2/2(3.35)
with zero mean and unit variance and multiply these numbers with √∆t
∆¯
wt=
n1
.
.
.
n6N
√∆t. (3.36)
The algorithm we use to calculate these Gaussian random numbers is presented
in section 3.3.3. The representation (3.32) of stochastic differential equation is
equivalent to equation (3.25) for ∆t→dt, which one can easily verify by calcu-
lating the first two moments:
hdxi= [f+∇D]dt+Hhd¯
wi
(3.34)
= [f+∇D]dt,
hdx⊗dxi=h(Hd¯
w)⊗(Hd¯
w)i(3.34)
=HHTdt
(3.33)
= 2Ddt. (3.37)
48
3.3 Brownian dynamics simulation
After performing the Euler step (3.32), we can calculate the configurational up-
date xt+∆t=xt+ ∆xt. The main problem of this direct implementation is the
calculation of the divergence of the diffusion tensor ∇Dat each time step, which
requires enormous computational time. To avoid the explicit calculation of this
term, we write down the discretized Itô stochastic differential equation (3.32) in
the corresponding equivalent Stratonovich form (see section 3.1.2)
∆xt=f(xt)∆t+H(xt)◦∆¯
wt.(3.38)
The advantage of this form is that the explicit calculation of the spurious drift is
not necessary. On the other hand, the quantities in equation (3.38) are evaluated
at the intermediate time step t+1
2∆t. Unfortunately, this requires the knowledge
of the configurational update xt+1
2∆xt, which is not available. To overcome
this problem, we use a predictor-corrector scheme [114, 115, 116, 117]. First, we
perform an Euler step without spurious drift
∆¯
xt=f(xt)∆t+H(xt)∆¯
wt.(3.39)
Then we use the intermediate configuration xt+ ∆¯
xto calculate the corrector
step
∆xt=1
2[f(xt) + f(xt+ ∆¯
xt)] ∆t
+1
2h1+D(xt+ ∆¯
xt)D−1(xt)iH(xt)∆¯
wt.(3.40)
Now, we can compute the final configuration xt+∆t=xt+ ∆xt. Note that the
expression H(xt)∆¯
wtin the corrector step (3.40) is the same as in the predictor
step (3.39), thus computation time can be saved by calculating this term only
once. Note that the above scheme is actually correct to order ∆t2for a constant
diffusion matrix. There are also higher order integration schemes of Runge-Kutta
type available for stochastic differential equations [118]. However, these methods
require the explicit calculation of the spurious drift ∇D.
3.3.2 The Cholesky decomposition
In section 3.3.1, we have introduced the matrix H, which is the solution of
2D=HHT. In general every symmetric and positive defined matrix A(as the
diffusion matrix Dis) can be decomposed in
A=LLT,(3.41)
49
3. BROWNIAN MOTION
where Lis a lower triangular matrix. The single elements of this matrix can be
determined recursively by [119, 120]
Lii = aii −
i−1
X
k=1
L2
ik!1
2
,
Lij =1
Ljj
Aij −
j−1
X
k=1
LikLjk
for i > j, (3.42)
which is called Cholesky decomposition.
3.3.3 Gaussian random numbers: the Box-Muller transform
In this section we will describe how we compute Gaussian random numbers.
Consider a random number xwith a given probability distribution p(x). The
distribution for a transformed variable y=y(x)is then easily given by
p(y) = p(x)
dx
dy
.(3.43)
More generally, let us consider the random numbers x1, x2, . . . with the joint prob-
ability distribution p(x1, x2, . . . ). If each of the transformed variables y1, y2, . . .
is a function of all the xi, then the joint probability distribution of the yiis given
by [97, 114]
p(y1, y2, . . . ) = p(x1, x2, . . . )
∂(x1, x2, . . . )
∂(y1, y2, . . . )
,(3.44)
where |∂(x1, x2, . . . )/∂(y1, y2, . . . )|denotes the Jacobian determinant of xiwith
respect to yi. Let us now consider two independent random numbers x1and x2,
which are uniformly distributed in the interval [0,1]. The transformation [121]
y1=q−2 ln x1cos 2πx2,
y2=q−2 ln x1sin 2πx2,(3.45)
with the inversion
x1= exp −1
2(y2
1+y2
2),
x2=1
2πarctan y2
y1
,(3.46)
yields the following probability distribution [121]
p(y1, y2) = 1
√2πe−1
2(y2
1+y2
2).(3.47)
50
3.3 Brownian dynamics simulation
The transformation (3.45) is referred to as the Box-Muller transform, which gen-
erates two independent random numbers satisfying a Gaussian distribution with
zero mean and unit variance.
In order to make this algorithm more efficient, one usually picks two ran-
dom numbers ¯x1and ¯x2, which are homogeneously distributed on the unit circle
(¯x2
1+ ¯x2
2<1) instead of the unit square [122]. The probability to find the radial
coordinate r=q¯x2
1+ ¯x2
2in the interval (r, r +dr)is then given by 2rdr. There-
fore, robeys the distribution p(r)=2r. Then, the variable ris transformed to
s=r2. Satisfying the condition p(r)dr=p(s)ds, where ds=d(r2) = 2rdryields
p(s)=1. Therefore, s= ¯x2
1+ ¯x2
2is uniformly distributed in the interval (0,1)
and can be used instead of x1. Equally so, the angle φ= arctan ¯x2/¯x1divided
by 2πis uniformly distributed in the interval (0,1) and independent of s. So,
we can replace x2by φ/(2π). Together with the relations cos φ= ¯x1/√sand
sin φ= ¯x2/√swe can write equation (3.45) as
y1=q−2 ln s/s¯x1,
y2=q−2 ln s/s¯x2.(3.48)
The advantage is that calculating the trigonometric functions directly can be
avoided. Following this polar form, the Box-Muller method then reads: First,
we pick two random numbers ¯x1and ¯x2, which are uniformly distributed in the
interval (−1,1). Then, we calculate s= ¯x2
1+ ¯x2
2and accept (¯x1,¯x2)only if s < 1
(otherwise, we pick new numbers). Finally, we obtain two independent Gaussian
random numbers y1and y2by yi=q−2 ln s/s¯xi.
51
3. BROWNIAN MOTION
52
4
Semiflexible polymer in confined Poiseuille flow
«Everyone is sure of this (that er-
rors are normally distributed), Mr.
Lippman told me one day, since the
experimentalists believe that it is
a mathematical theorem, and the
mathematicians that it is an exper-
imentally determined fact.»
Henri Poincaré
53
4. SEMIFLEXIBLE POLYMER IN CONFINED POISEUILLE FLOW
Figure 4.1: Adapted with permission from Ref. [37]. Copyright 2012 American
Chemical Society. (a) Schematic sketch of the actin filament and of the experimen-
tal setup of Ref. [37]. The filaments are under the influence of a pressure driven
Poiseuille flow and they are only considered within the focal plane. (b) Center-of-
mass probability distributions of Ref. [37] for different flow velocities. The plot
shows the distributions of filaments along the channel cross-section for a half chan-
nel with the channel center at ycm = 0µm and the channel wall at ycm = 5.5µm.
4.1 Introduction
Polymers are molecules made from repeating units of identical or nearly identical
compounds called monomers that are linked together by a series of covalent bonds
[123]. They are essential and omnipresent in everyday life, from plastics and
elastomers on the one hand to natural biopolymers on the other. Biopolymers
such as proteins and nucleic acids play fundamental roles in biological processes.
For example, the specific sequence of the four nucleotide bases adenine, guanine,
cytosine and thymine in a DNA molecule makes each of us unique. Besides
the building block of life, the DNA, actin is also an important component of
cells, which is part of the cytoskeleton and play an important role in cellular
motility and mechanical stability. Typically, actin filaments have a diameter
of a few nanometer, an averaged length of about 5−10µm and a persistence
length of about 13µm [5]. Due to this properties actin filaments can be used
as model systems for semiflexible polymers and, therefore, are central in recent
experimental studies [5, 37]. In the experiments in Ref. [37] single αactin
filaments with an averaged length of about 8µm are investigated in a squared
54
4.1 Introduction
microchannel [see figure 4.1(a)]. The channel has a width of 11µm, a height of
10µm and only filaments within the focal plane (≈0.5µm) are considered. The
filaments are driven out of equilibrium using a pressure driven Poiseuille flow
and analyzed at the end of the microchannel. Besides the conformation and
orientation dynamics of the filament the steady state center-of-mass probability
distribution was measured, which shows an interesting behavior for different flow
velocities [see figure 4.1(a)].
In this chapter we address the flow induced migration of a semiflexible polymer
across streamlines by Brownian dynamics simulations. This effect is commonly
called cross-streamline migration. We explain it by analyzing the Smoluchowski
equation (3.30) for the probability distribution of the polymer’s center-of-mass.
In the regime of low Reynolds numbers flexible polymers show cross-streamline
migration, which has been intensively studied. A depletion of polymers near
walls in planar shear flow as well as in Poiseuille flow was observed in computer
simulations [44, 45, 46, 47, 48, 124, 125, 126] and experiments [49, 50, 51]. Since a
polymer interacts hydrodynamically with bounding walls, migration towards the
centerline occurs, where the thickness of the depletion layer increases with the
strength of the flow. This explanation was confirmed by analytical arguments [53,
54] and it describes well repulsion from bounding walls. In addition, simulations
reported migration away from the channel’s centerline [45, 46, 125, 126], especially
under strong confinement [46]. As a result, the maximum concentration of the
polymer occurs at a finite distance from the centerline and the depletion in the
center increases with flow strength [46, 125, 126]. The spatially varying shear rate
of a Poiseuille flow changes orientation and conformation of a polymer within the
channel. Close to the walls polymers are stretched and at the center they are
coiled [44, 44, 45, 46, 47, 125]. This gives rise to a position-dependent diffusivity
and thereby a diffusion current away from the center [45, 46, 53, 125, 126], which
generates the observed bimodal concentration profile.
Cross-streamline migration of semiflexible polymers, such as α-actin filaments,
is much less studied. Since their bending rigidity significantly determines flow-
induced conformations in the Poiseuille flow, we also expect an influence on cross-
streamline migration. Indeed, recent experimental [37] and simulation [38, 52]
studies report a much more pronounced bimodal concentration profile across the
channel when compared to flexible polymers. This effect occurs even under less
strong confinement. As for flexible polymers, the concentration profile is also
explained by the competition between hydrodynamic polymer-wall interactions
and enhanced diffusion away from the centerline.
In this chapter we reexamine the problem and perform Brownian dynamics
simulations for a bead-spring chain with bending elasticity. Hydrodynamic in-
teractions between the point-like beads are taken into account by the two-wall
Green tensor of the Stokes equations. We indeed observe the bimodal distribu-
55
4. SEMIFLEXIBLE POLYMER IN CONFINED POISEUILLE FLOW
tion for the center-of-mass and present results on the orientational order of the
end-to-end vector. Generalizing the approach of Ma and Graham [53], we formu-
late and interpret our findings with the help of a Smoluchowski equation for the
center-of-mass probability distribution by carefully analyzing all contributions to
the probability current. In contrast to flexible polymers, we show that the de-
terministic drift current, where hydrodynamic interactions along the polymer are
essential, mainly determines the bimodal distribution across the channel, whereas
diffusional currents become less important with increasing polymer stiffness. This
section is based on publication [A].
56
4.2 Modeling
Figure 4.2: A semiflexible polymer modeled by a bead-spring chain with bending
rigidity is confined between two parallel planar plates with distance 2W. The space
vector ripoints from the centerline at z= 0 to bead i.
4.2 Modeling
In this section we describe how we model a semiflexible polymer in a pressure
driven Poiseuille flow. We first introduce the semiflexible polymer as a bead-
spring chain with bending elasticity, explain how it couples to the Poiseuille
flow, and finally summarize details of our Brownian dynamics simulations. The
semiflexible polymer is confined between two parallel planar walls at positions
z=−Wand z=W, which are infinitely extended in the x−yplane (see
figure 4.2). For simplicity, we assume that all parts of the chain only move in the
y−zplane at x= 0. This corresponds to the experiments reported in Ref. [37],
where only filaments were recorded and analyzed whose motion occurred in a
narrow region (±0.5µm) around the focal plane of the microscope. Furthermore,
this also reduces our computational efforts considerably. Between the walls a
pressure driven Poiseuille flow is created, which we express as
u0(r) = v0"1−z
W2#ey,(4.1)
where v0denotes the maximum velocity at the centerline z= 0,eyis the unit
vector along the yaxis, and Wthe distance between the centerline and the walls.
4.2.1 The worm-like chain model
A very popular model to describe a polymer is the so called worm-like chain
model, which was introduced by Kratky and Porod in 1949 [127]. In this model
the polymer is described by the curve r(s), where sis the arc length in the interval
[0, L]and Lthe contour length of the polymer (see figure 4.3). It can be shown,
57
4. SEMIFLEXIBLE POLYMER IN CONFINED POISEUILLE FLOW
Figure 4.3: Worm-like chain model: the polymer is described by the curve r(s)
parametrized by the arc length s.t(s)denotes the tangential vector, Lthe total
contour length and Rthe end-to-end vector
that the correlation of the tangent vector t(s) = ∂r(s)/∂s,|t(s)|= 1 is given by
[127]
ht(s)·t(s+ ∆s)i=e−∆s/Lp.(4.2)
Here, Lpis the persistence length of the polymer, which quantifies its stiffness.
In principle, the mechanical properties of the polymer can be divided into three
classes [128]. Is the contour length much larger than the persistence length (L
Lp), then the polymer is considered to be flexible, is the contour length much
smaller than the persistence length (LLb)the polymer is considered to be
stiff and is the contour length on the same order as the persistence length (L≈Lp)
the polymer is considered to be semiflexible. To illustrate this, we have a look
at the correlation of the tangent vectors at the beginning and at the end of the
filament
ht(0) ·t(L)i=e−L/Lp.(4.3)
One immediately sees that for flexible polymers (LLp)the correlation is zero
and, therefore, the configuration at the end does not depend on the configuration
at the beginning. For semiflexible polymers (L≈Lp)the correlation is in the
order of e−1. In this case the configuration at the end depends on the configuration
at the beginning, but the contour does not has to be straight. For stiff polymers
(LLb)the correlation is equal to one and the contour of the polymer is
completely straight.
58
4.2 Modeling
Figure 4.4: Bead-spring model: the polymer is described by Nbeads, which are
connect by (N−1) frictionless springs. Here, qn=rn+1 −rndenotes the n-th
chain segment.
Any deformation of a short segment dsof the filament has an energy cost as-
sociated with it. By analogy to Hooke’s law, the relationship should be quadratic
for small deformations. If we neglect a torsional deformation, the free energy
(associated with bending the chain) is given by [129]
H=UB=1
2
L
Z
0
κ ∂t(s)
∂s !2
ds. (4.4)
Here, κ=LpkBTdenotes the bending rigidity and ∂t/∂s the change of the
tangential vector along the segment ds.
4.2.2 The bead-spring model
Since the hydrodynamic interactions are well described for spherical or point-like
particles (see section 2.4), we are using a discretized representation of the worm-
like chain model, the so called bead-spring model. In this coarse grained model
the filament is composed of Nbeads, which are connected by (N−1) springs
(see figure 4.4). We consider that the connecting springs are very thin, so that we
can neglect their friction with the surrounding fluid. For this situation, we only
have to take the hydrodynamic interactions of the beads into account, for which
we know the friction coefficients. The main idea behind this discretisation is to
describe the stretching and bending forces such that they act on finite discrete
points on the filament, such that we can apply the methods we derived in section
2.5.2.
59
4. SEMIFLEXIBLE POLYMER IN CONFINED POISEUILLE FLOW
To account for the constant contour length of the polymer, we apply Hooke’s
law on each chain segment. So, the total energy stored in stretching the polymer
is determined by the sum over the squared deviations of each segment from its
equilibrium length l=hqni
US=1
2H
N−1
X
n=1
(qn−l)2,(4.5)
where qn=|qn|is the length of the n-th chain segment
qn=rn+1 −rn.(4.6)
A single chain segment cannot be bended, but neighbouring segments may be
inclined to each other. To obtain an expression for the discretized bending energy,
we simply substitute dtwith (qn+1/qn+1 −qn/qn)and dswith lin equation (4.4).
Therefore, we obtain [130]
UB=κ
l
N−2
X
n=1 1−qn+1
qn+1 ·qn
qn!.(4.7)
Finally, the force acting on bead ndue to bending and stretching the bead-spring
chain follows from
Fn=−∇n(UB+US),(4.8)
where ∇nis the nabla operator with respect to the position vector rn.
4.2.3 Equations of motion
The beads in the model polymer also experience stochastic forces due to collisions
with fluid molecules from the surrounding solvent. As usual, we take them into
account as Gaussian white noise that obeys the fluctuation-dissipation theorem
(see section 3.1). Using a compact notation, the resulting Langevin equation in
differential form reads [114, 131]
dx=u0+1
kBTDF +∇Ddt+Hn√∆t , (4.9)
with D=kBTµ=1
2HHT.(4.10)
Here we have introduced 3Ndimensional vectors to collect the positions of the
beads, x= (r1, ..., rN), the forces acting on them, F= (F1, ..., FN), and
the velocities from the external flow field at the positions of the beads, u0=
(u0(r1), ..., u0(rN)). The mobilities µtt
ij form the mobility tensor µand Dstands
60
4.2 Modeling
for the generalized diffusion tensor, where kBis Boltzmann’s constant and Tthe
absolute temperature. The tensor His determined such that it fulfills equation
(4.10) (see section 3.3.2). Finally, ndenotes a random vector with expectation
values for the mean and the second moment,
hni= 0 and hn⊗ni=1.(4.11)
The first term on the right hand side in equation (4.9) describes the drift motion
of the beads, including the noise-induced spurious drift (∇D), and the second
term their Brownian motion [131]. To numerically integrate Eqs. (4.9), we use a
predictor-corrector scheme summarized in section 3.3.1 to avoid the explicit calcu-
lation of ∇D. An Euler step without the spurious drift generates an intermediate
configuration ¯
x(t+ ∆t) = x(t)+∆¯
xwith
∆¯
x=u0+1
kBTDF ∆t+Hn√∆t. (4.12)
Then the corrector step calculates the final configuration x(t+ ∆t) = x(t) + ∆x
with
∆x=1
2u0+¯
u0+1
kBTDF +¯
D¯
F∆t
+1
2h1+¯
DD−1iHn√∆t. (4.13)
Here the symbol "¯" means that these quantities are calculated for the intermediate
configuration ¯
x.
4.2.4 System parameters
In the experiments in Ref. [37], α-actin filaments with a contour length of about
L= 8µmwere analyzed in a microchannel with a square cross section of width
2W= 10µm. The actin filament is a thin semiflexible biopolymer with a diameter
of 8nm and a persistence length of Lp= 13µm. This gives the experimental value
Lp/L = 1.6. The velocity of the imposed Poiseuille flow at the centerline varied
from v0= 0mm/sto v0= 2.4mm/s. To be close to the experimental values,
we choose the following parameters for our simulations. We set the width of the
channel to the experimental value, 2W= 10µm, and refer all our lengths to W.
The polymer consists of N= 16 monomers and the equilibrium bond length to
the nearest neighbor is l/W = 0.1so that the total contour length amounts to
L/(2W) = 0.75. The bead radius is set to a/W = 0.005 meaning l/a = 20 which
allows to treat the beads as point particles when hydrodynamic interactions are
calculated. The spring constant is set to H= 1.2·104(kBT/W2)and we will
61
4. SEMIFLEXIBLE POLYMER IN CONFINED POISEUILLE FLOW
explore how the model polymer behaves when the persistence length varies from
Lp/L = 1 to Lp/L = 16. We use the viscosity of water, ηH2O= 0.01P, and
perform the Brownian dynamics simulations at a temperature T= 300K. The
imposed flow velocity at the centerline varies from v0= 0mm/sto v0= 2.5mm/s.
The configurational update of the polymer is calculated after a time step of
∆t= 10−5s. Since the evaluation of the two-wall Green tensor at each time step is
very time consuming, we stored its values on a grid with 200×200 grid points and
with a mesh size of 0.01W. Tensor values at intermediate locations are determined
from the grid points by linear interpolation. The probability distribution and
all averages presented in the following sections are calculated after the polymer
reached its steady state for various independent initial conditions.
62
4.3 Results
Figure 4.5: Snapshots of the bead-spring chain for different persistence lengths
at the centerline (top) and near the wall (bottom) at v0= 2.5mm/s. (a) Lp/L = 1,
(b) Lp/L = 4, and (c) Lp/L = 16.
4.3 Results
The imposed parabolic Poiseuille flow of equation (4.1) drives the filament out
of equilibrium and determines its steady state, which we quantify, e.g., by the
center-of-mass probability distribution n(zC). In this section we will investigate
how flow strength, hydrodynamic interactions, polymer stiffness, and the ratio
between bond length and bead radius, l/a, influence the polymer conformation,
the center-of-mass distribution n(zC), and the orientational order of the filament’s
end-to-end vector.
Figure 4.5 shows different snapshots of the filament near the centerline and
near a bounding wall. The filament constantly moves up and down along the
lateral direction. Most of the time, it is aligned parallel to the flow direction
although it is not perfectly straight. At the centerline, we also observe that the
filament bends to a typical U shape that traces the Poiseuille flow profile (see
figure 4.5, top). This configuration is not stable when the filament moves away
from the centerline. Outside the center, the filament also tumbles when the front
end moves closer to the walls where it experiences a smaller flow velocity (see
figure 4.5, bottom). After one tumbling event, the filament aligns along the flow
field until thermal motion moves the front end again closer to the wall. Increasing
63
4. SEMIFLEXIBLE POLYMER IN CONFINED POISEUILLE FLOW
v0also increases the local shear rates and thereby the tumbling frequency.
4.3.1 Polymer conformation and orientation
To study how the Poiseuille flow influences the filament shape and orientation
within the channel, we calculate the end-to-end vector Rand determine its ori-
entational order around the channel axis by the order parameter S
S(zC) = hcos(2θ)i.(4.14)
Here θis the angle between Rand the yaxis and Sdescribes the orientational
order in a two-dimensional system [132]. A polymer with randomly oriented
end-to-end vector gives S= 0, whereas S= 1 means that the polymer is per-
fectly aligned along the flow direction and S=−1indicates perfect orientation
perpendicular to the flow direction.
Figures 4.6 (a) and (b) show, respectively, the orientational order parameter
and the average end-to-end distance hRifor different flow strengths. Without
flow, Sis nearly zero within the region |zC|/W ≤0.25, where the filament can
freely rotate and where it assumes all orientations with equal probability. Outside
this region, the filament’s orientation is more and more restricted due to the steric
polymer-wall interaction and Sgrows monotonically until it reaches S= 1 at the
walls. For a rigid rod, one readily derives the order parameter as a function of
|zC|/W:
S(zC) =
1,|zC|/W ≤1−L/2W
u√1−u2
arcsinu,1−L/2W < |zC|/W ≤1,(4.15)
where uis given in equation (4.17). The inset of figure 4.6(a) shows very good
agreement with our simulations. Deviations result from thermal fluctuations of
the filament. They also explain that the end-to-end distance in figure 4.6(b)
deviates slightly from one. In the presence of external flow both the order pa-
rameter and the end-to-end distance have a minimum in the center of the channel
which we attribute to the U-shaped conformations illustrated at the bottom of
figure 4.5. Sand hRithen increase monotonically until a nearly constant value
or plateau is reached around |z|/W = 0.2. Here the filament is aligned along the
flow (S > 0.7) and tumbling events cause the end-to-end distance to be smaller
than for zero flow field. The end-to-end distance shows a weak minimum near
the wall which we attribute to an increasing number of tumbling events since
the shear stress on the filament within the channel becomes largest. Finally,
very close to the wall, when tumbling events can no longer occur, the filament
is perfectly aligned with the flow. The order parameter Sreaches its maximum
64
4.3 Results
Figure 4.6: Order parameter Sand average end-to-end distance hRiplotted
versus the lateral center-of-mass position zC/W: (a),(b) for a fixed stiffness Lp/L =
1and different flow velocities; (c),(d) for a fixed flow velocity v0= 2.5mm/s and
different stiffnesses Lp/L. Inset in (a): For v0= 0mm/s, the simulated order
parameter Sis compared to the analytic result of equation (4.15).
65
4. SEMIFLEXIBLE POLYMER IN CONFINED POISEUILLE FLOW
value one and the end-to-end distance even shows values above Lfor large flow
strengths. Interestingly, the plateau value of Sincreases with the flow strength
while hRidecreases. Higher shear rates align the filament better along the flow
direction but they also initiate more tumbling events which reduce hRibut not
S. We repeated our simulations without hydrodynamic interactions acting across
the filament and from the wall. The results agree well with the behavior discussed
in Figs. 4.6 (a) and (b). The same behavior occurs when we reduce the ratio
l/a of bond length to bead radius from the value l/a = 20 used in most of our
simulations to l/a = 5, where hydrodynamic interactions between the beads are
stronger. This indicates that overall shape and orientation of the filament is de-
termined by its overall friction with the solvent and in particular by the applied
Poiseuille flow.
For stiffer filaments the end-to-end distance hRiincreases as indicated in figure
4.6(d). The plateau value of the order parameter in figure 4.6(c) does not change
strongly when the stiffness increases. Interestingly, the minimum value for Sat
the centerline first decreases to even negative values and then increases again.
The negative values are caused by the U-shaped conformation where the end-to-
end vector is perpendicular to the centerline. A persistence length LP/L between
4 and 8 seems to be the best choice for realizing the U-shaped conformation.
4.3.2 Steady state center-of-mass distribution
One quantity to characterize the filament within the microchannel is the steady-
state probability distribution n(zC)for its center-of-mass across the channel. We
always normalize it to one, RW
−Wn(z)dz= 1. The experiments of Ref. [37] clearly
revealed that the distribution depends on how strongly the filament is driven out
of equilibrium. In figure 4.7 we plot the distribution for different flow strengths
v0from the centerline to the bounding wall. Without flow the center-of-mass
location zCis equally distributed within the interval |zC|/W ≤1−L/2W= 0.25.
In thermal equilibrium, such a behavior is expected for the freely moving filament
with arbitrary orientation. For locations zCcloser to the walls, the orientation of
the stiff filament is more and more restricted due to the steric interaction with
the wall and the distribution decreases to zero at the walls. For a rigid rod, one
readily derives the distribution
n(zC) =
π
2(π−L/W),|zC|/W ≤1−L/2W
arcsin u
π−L/W,1−L/2W < |zC|/W ≤1
(4.16)
66
4.3 Results
Figure 4.7: Center-of-mass probability distribution plotted from the centerline
(zC= 0) to the wall (zC=W) for different flow velocities v0and a persistence
length of LP/L = 1. The inset compares the simulation results for zero flow velocity
v0to the hard-needle distribution of equation (4.16).
with
u=1−|zC|/W
L/2W,(4.17)
which is compared to the simulated profile in the inset of figure 4.7. Deviations
result from thermal fluctuations of the filament. When we turn on the Poiseuille
flow (see figure 4.7, v0= 0.5mm/s), the probability distribution close to the wall
first increases. We understand this behavior. The tumbling events illustrated
in figure 4.5 (bottom) and initiated by the Poiseuille flow allow the filament to
move closer to the wall. Increasing the flow velocity further, the filament is
depleted more and more at both walls. This is initiated by the hydrodynamic
repulsion of the filament from the wall, which pushes it away from the wall [45,
53, 54]. The repulsion becomes clear with the help of figure 2.5 on the left. Due
to the strong shear flow close to the wall, the filament is under tension which
initiates a flow field that drives the model polymer away from the wall [45, 53].
In particular, at both ends of the rod tensional forces point into the rod. Each of
these forces initiates a flow field similar to the one of figure 2.5(a) which drives
the other end of the rod away from the wall. This phenomenon is different from
the orientational lift forces reported in Ref. [54]. The most pronounced behavior
67
4. SEMIFLEXIBLE POLYMER IN CONFINED POISEUILLE FLOW
Figure 4.8: Center-of-mass probability distributions for different persistence
lengths LP/L at a fixed center flow velocity v0= 2.5mm/s.
of the system is the depletion of the distribution at the centerline. It leads to a
bimodal probability density with a maximum at a finite distance from the center.
The maximum increases with the flow strength. At the maximum, the migration
towards the centerline induced by the hydrodynamic filament-wall interaction
has to be cancelled by a probability current away from the centerline. Based on
literature [45, 48, 53], Refs. [37, 38, 52] attribute such a current to the spatially
varying diffusivity of the filament across the channel: the U-shaped conformation
at the centerline has a larger diffusivity than the straight filament outside the
center. In section 4.4 we will demonstrate that the current away from the center
is mainly due to deterministic drift motion which becomes dominant over the
diffusional current when the rigidity of the filament increases. Figure 4.8 shows
the center-of-mass distributions at a fixed flow velocity v0= 2.5mm/s for different
persistence lengths LP. Here, the depletion layer at the walls decreases slightly
with increasing bending rigidity, whereas the bimodal distribution becomes more
pronounced. For larger rigidity, the migration away from the centerline increases.
As a result, the depletion at zC= 0 becomes stronger and the position of the
maximum is shifted towards the walls which causes the decreasing depletion layer.
We will explain this behavior in section 4.4. For LP/L = 16, the local minimum
and the maximum of the distribution differ by a factor of two, which is comparable
to the results of Ref. [37].
Most of our simulations are done for a bead distance to radius of l/a = 20.
68
4.3 Results
Figure 4.9: Center-of-mass probability density for different ratios l/a at a fixed
persistence length LP/L = 2 and a fixed flow velocity v0= 2.5mm/s. In order to
change the ratio l/a, we keep the length Land the number of beads Nconstant
and vary the radius aof the beads. This corresponds to increasing the thickness
of our model polymer.
This might mimic the thin actin filament. However, it also underestimates the
overall friction between the filament and the solvent. In figure 4.9 the center-of-
mass distributions for different bead sizes aat a fixed flow velocity v0= 2.5mm/s
are shown. Increasing the overall friction of the model polymer with aenhances
tensional forces within the polymer and thereby the drift currents discussed in
section 4.4. As a result, the bimodal distribution and the depletion at the walls
become more pronounced.
Hydrodynamic interactions between different parts of the model polymer and
with the wall are crucial for the observed bimodal distribution and the depletion
at the walls. In figure 4.10 we show the resulting center-of-mass distributions at
different flow strengths when hydrodynamic interactions are switched off during
the simulations. The equilibrium profile at zero flow field is the same as in figure
4.7 since it should not depend on dynamic properties such as hydrodynamic inter-
actions. For nonzero flow velocity, the minimum in the center vanishes compared
to figure 4.7 whereas the depletion close to the wall is less pronounced. It even
decreases with increasing v0since the hydrodynamic repulsion from the wall is
missing and the tumbling polymer is stronger confined due to larger viscous shear
stresses [37].
69
4. SEMIFLEXIBLE POLYMER IN CONFINED POISEUILLE FLOW
Figure 4.10: Center-of-mass probability distribution simulated without hydrody-
namic interactions for different flow velocities v0and persistence length LP/L = 1.
4.4 Kinetic theory for a semiflexible polymer
In this section we derive and analyze a Smoluchowski equation for the center-
of-mass probability distribution of the bead-spring chain with bending elasticity.
We thereby generalize the approach of Ma and Graham who formulated a kinetic
theory for a bead-spring dumbbell [53]. We derive the Smoluchowski equation
in section 4.4.1 and study and discuss in detail the different contributions to the
lateral center-of-mass current in section 4.4.2 in order to identify the cause for
cross-streamline migration.
4.4.1 Smoluchowski equation for center-of-mass current
The probability distribution ψ(r1, ..., rN;t)for finding the bead-spring chain in a
state determined by the bead coordinates r1, ..., rNat time tis governed by the
Smoluchowski equation
∂
∂tψ=−∇i·ji,(4.18)
where we have introduced the probability density current of bead i
ji=˙riψ−Dij∇jψ. (4.19)
70
4.4 Kinetic theory for a semiflexible polymer
Here,
˙ri=u0(ri) + µtt
ijFj(4.20)
denotes the deterministic drift velocity and Dij =kBTµtt
ij is the diffusion ten-
sor connected to the mobility tensor by the Einstein relation. Note that we use
the convention where we sum over bead indices that occur twice in an expres-
sion. Since we are interested in the center-of-mass probability distribution, we
introduce the respective center-of-mass position and bond vectors,
rC=1
N
N
X
i=1
ri,
qn=rn+1 −rn, n = 1, . . . , N −1.(4.21)
Using the new coordinates and
∇i=1
N∇rC+∇qi−1−∇qi,(4.22)
we rewrite the Smoluchowski equation (4.18) as
∂ψ
∂t =−∇rC·(˙rCψ−DrC∇rCψ−DT
qi∇qiψ)
−∇qi·(˙qiψ−Dqi∇rCψ−¯
Dij∇qjψ).(4.23)
We have introduced the respective deterministic velocities for center-of-mass and
bond i,
˙rC=1
NX
i
(u0(ri) + µtt
ijFj)(4.24)
˙qi=u0(ri+1)−u0(ri) + µtt
i+1,jFj−µtt
ijFj(4.25)
(4.26)
and various diffusivities
DrC=1
N2
N
X
i,j
Dij,(4.27)
Dqi=1
N
N
X
j
Di+1,j −Dij,(4.28)
¯
Dij =Di+1,j+1 −Di,j+1 −Di+1,j +Dij.(4.29)
The average over all diffusion tensors Dij is the Kirkwood diffusivity DrC, whereas
Dqiand ¯
Dij are diffusivities related to the bond vectors. We now write the full
71
4. SEMIFLEXIBLE POLYMER IN CONFINED POISEUILLE FLOW
probability distribution ψin the new coordinates rC,q1,...,qN−1and introduce
the center-of-mass probability distribution
n(rC, t) = Z. . . Zψ(rC,q1,...,qN−1;t)dq1. . . dqN−1.(4.30)
Integrating equation (4.23) over all bond vectors qi, one is able to derive a Smolu-
chowski equation for the center-of-mass,
∂n
∂t =−∇rC·jC,(4.31)
with the center-of-mass probability current
jC=hh˙rCiq−∇rChDrCiq+h∇rCDrCiq
+D∇qiDT
qiEqn−hDrCiq∇rCn. (4.32)
The brackets h. . . iqdenote an ensemble average over all bond vectors qior poly-
mer conformations with fixed rC,
hAiq=Z. . . ZAˆ
ψdq1. . . dqN−1,(4.33)
where
ˆ
ψ(rC,q1,...,qN−1;t) = ψ(rC,q1,...,qN−1;t)
n(rC, t).(4.34)
Note that the second and third terms on the right-hand side of equation (4.32) are
not necessarily the same. In integrating equation (4.23) over all bond vectors qi,
we use the reasonable assumption that all surface terms vanish. We therefore can
immediately skip the terms in the second line of equation (4.23). In the remaining
terms, ψis replaced by nˆ
ψfrom equation (4.34) and partial integrations are
performed such that no gradient acts on ˆ
ψand equation (4.33) can be applied.
This procedure finally gives the center-of-mass probability current (4.32).
Equations (4.31) and (4.32) determine the center-of-mass probability distri-
bution. Whereas the last term on the right-hand side of equation (4.32) describes
conventional diffusion [53], the remaining terms proportional to the distribution
nformally are drift terms. However, only the first term is due to deterministic
motion of the center-of-mass. The second to fourth terms result from the diffu-
sional currents in equation (4.19). As we demonstrate explicitly in the following
section, all of these terms, in principle, can lead to cross-streamline migration.
Here we add some general remarks.
72
4.4 Kinetic theory for a semiflexible polymer
The first term on the right-hand side of equation (4.32) describes migration
due to hydrodynamic interactions and the applied external flow. Without hydro-
dynamic interactions the cross mobilities µtt
ij vanish and equation (4.24) reduces
to ˙rC=Pi[u0(ri)]/N. The frictional term Piµtt
ii Fi= (PiFi)/(6πηa) = 0does
not contribute to the center-of-mass motion since the total force acting on the
filament has to be zero. As a result, without hydrodynamic interactions cross-
streamline migration cannot occur. The second term on the right-hand side of
equation (4.32) is the divergence of the Kirkwood diffusivity and a natural can-
didate for migration away from the centerline as demonstrated in Refs. [45, 53].
Finally, the third and fourth term vanish whenever the divergence of the diffusion
tensors or mobilities is zero. In particular, this is the case in an unbounded fluid
when one treats hydrodynamic interactions on the level of the Oseen tensor or
in the next higher order via the Rotne-Prager tensor. In our case, the third and
fourth term will be non-zero due to wall-induced hydrodynamic interactions.
4.4.2 Analysis of the lateral center-of-mass current
Due to translational symmetry along the channel direction, the center-of-mass
distribution does not vary along the yaxis and the zcomponent of the center-
of-mass current, i.e., the current across the channel is
jC,z =hh˙zCiq−∂zChDrC,zziq+h∂zCDrC,zziq
+D∂qi,yDT
qi,zy +∂qi,zDT
qi,zzEqn
−hDrC,zziq∂zCn. (4.35)
In steady state the center-of-mass current jCis constant. Since the current at
the walls has to vanish, jC,z is zero everywhere across the channel. In figure 4.11,
we plot all contributions of the center-of-mass current jC,z proportional to nfor
different bending rigidities of the filament. The sum of these currents balances the
diffusional current −hDrC,zziq∂zCn. Concentrating on LP/L = 1, we recognize
that close to the wall the deterministic drift current h˙zCiqnis directed away
from the wall, e.g., it is positive at zc=−W, and it dominates all the other
currents. Hence, the hydrodynamic repulsion from the walls is responsible for
the depletion at the walls. On the other hand, close to the centerline, the current
is directed towards the wall and, therefore, causes depletion at the centerline.
For LP/L = 1, the diffusional current −∂zChDrC,zziqndue to the gradient of the
conformation-averaged Kirkwood diffusivity also points away from the centerline
and contributes to the observed depletion. We understand this since the U-shaped
conformation close to the centerline has a larger diffusivity than the straight
conformation occuring outside the centerline [45, 53]. However, already at a
73
4. SEMIFLEXIBLE POLYMER IN CONFINED POISEUILLE FLOW
Figure 4.11: Plots of the different contributions of the lateral center-of-mass
current across the channel at flow velocity v0= 2.5mm/s for different bending
rigidities: (a) LP/L = 1, (b) LP/L = 2, (c) LP/L = 4, and d) LP/L = 16.
stiffness of LP/L = 2, the deterministic drift current is clearly the dominant part
for causing centerline depletion and at LP/L = 16 all the diffusional currents are
negligible. The experiments of Ref. [37] were performed for Lp/L = 1.6. So, we
74
4.5 Summary and conclusion I
conclude that the observed centerline depletion is mainly due to the deterministic
drift current. The third term h∂zCDrC,zziqnis approximately zero for each Lp/L
and the fourth term D∂qi,yDT
qi,zy +∂qi,zDT
qi,zzEqnleads to a current towards the
centerline due to hydrodynamic repulsion from the wall.
We summarize our results here. Whereas the filament at the centerline dis-
plays more compact conformations like the U shape, close to the walls it is mainly
aligned along the flow lines due to the large shear rate. Therefore, the friction co-
efficient for motion across the channel increases from the centerline to the wall and
its inverse, the Kirkwood diffusivity, decreases monotonically from the centerline
towards both walls. This initiates migration away from the centerline [45, 53, 54],
which can cause centerline depletion [45, 53, 125, 126], and was therefore used to
interpret the recent experiments of Steinhauser et al. [37, 52]. Here we demon-
strate that for increasing bending rigidity, centerline depletion is mainly caused
by the deterministic drift current. In the U-shaped conformation of figure 4.12(a)
the filament experiences bending forces. When it relaxes, different parts of the
filament interact hydrodynamically via flow fields initiated along the filament.
These hydrodynamic interactions cause the center-of-mass to move away from
the centerline. Increasing the rigidity of the filament also increases the bending
forces which explains why the deterministic drift currents in figure 4.11 close to
the centerline increase with LP/L.
The shear induced hydrodynamic repulsion of a dumbbell from a wall has
been treated in Refs. [45, 53, 54]. Similarly, close to the wall the shear flow
stretches the filament and creates tensional forces which initiate flow fields as
sketched in figure 2.5 on the left. They drive the filament away from the wall [see
figure 4.12(b)]. These tensional forces do not depend on the filament’s bending
rigidity. Therefore, in figure 4.11 the deterministic drift currents close to the wall
do not change with LP/L.
4.5 Summary and conclusion I
In this chapter we treated a paradigmatic model for studying how the properties
of a system change when it is driven out of equilibrium. Motivated by recent
experiments [37], we analyzed a single semiflexible polymer confined between two
planar walls and under the influence of an imposed Poiseuille flow with special
emphasis on the observed cross-streamline migration. We performed Brownian
dynamics simulations for a bead-spring chain with bending rigidity and used the
two-wall Green tensor of the Stokes equations to take into account hydrodynamic
interactions along the polymer and with the wall. We carefully analyzed how
polymer conformations, center-of-mass distribution, and orientational order of
75
4. SEMIFLEXIBLE POLYMER IN CONFINED POISEUILLE FLOW
Figure 4.12: (a) At the centerline, the relaxing U-shaped filament initiates flow
fields relative to the applied Poiseuille flow. The resulting hydrodynamic interac-
tions drive the filament away from the centerline. The strength of the flow is given
by the color code in arbitrary units. (b) Close to the wall the filament is under
tension. This initiates flow fields illustrated in figure 2.5 on the left that drive the
filament away from the wall.
the end-to-end vector within the channel depend on parameters such as the flow
strength and the stiffness of the model polymer. Analytic expressions for hard
needles at zero flow reproduce the simulation results and demonstrate how the
behavior of the model polymer changes when the Poiseuille flow is turned on.
Our results are in agreement with experiments [37] and simulations [38, 52] that
employ a different, particle based method to simulate the viscous environment
called multi-particle collision dynamics.
In particular, we observed the characteristic bimodal probability distribu-
tion for the polymer’s center-of-mass and showed that hydrodynamic interactions
along the model polymer and with the wall are essential for this distribution to oc-
cur, whereas shape and orientation of the filament are mainly determined by the
applied Poiseuille flow. Based on a Smoluchowski equation for the center-of-mass
probability distribution, we investigated cross-streamline migration and the origin
of the bimodal distribution in detail. Whereas the migration away from the wall
is due to hydrodynamic interactions with the bounding walls, in agreement with
Refs. [37, 38, 52] and work on flexible polymers [44, 44, 45, 46, 47, 49, 50, 125, 126],
we clearly identified a deterministic drift current as the major cause for migration
76
4.5 Summary and conclusion I
away from the centerline, especially when the bending stiffness of the model poly-
mer increases. The current is set up when bent conformations of the polymer relax
towards the straight filament. Diffusional currents due to a position-dependent
diffusivity become completely irrelevant with increasing rigidity of the polymers.
This demonstrates that bending rigidity leads to a clear difference in the behavior
of flexible and semiflexible polymers in Poiseuille flow and in the explanation of
the observed cross-streamline migration.
77
4. SEMIFLEXIBLE POLYMER IN CONFINED POISEUILLE FLOW
78
5
Nonlinear dynamics of spherical particles in
confined Poiseuille flow
«The scientist does not study na-
ture because it is useful to do so.
He studies it because he takes plea-
sure in it, and he takes pleasure in
it because it is beautiful.»
Henri Poincaré
79
5. NONLINEAR DYNAMICS OF SPHERICAL PARTICLES IN CONFINED
POISEUILLE FLOW
Figure 5.1: Systematic sketch of the classes of pair trajectories already observed
in linear shear flow (displayed in the center-of-mass frame). (a) passing trajectory:
two spheres pass each other and return to their initial lateral positions after the
binary encounter is completed. (b) swapping trajectory: the spheres swap their
lateral positions and leave in the same direction from where they came.
5.1 Introduction
An essential aspect for understanding microfluidic devices is the flow induced
cross-streamline migration of suspended objects in confined geometries. There-
fore, recent investigations have focused on the nonlinear dynamics in such geome-
tries. For example, hydrodynamic coupling under confinement with and without
imposed flow has been discussed close to a planar wall[32, 56, 133], close to a
planar fluid-fluid interface [134], between two planar walls[55], or in a spherical
cell [135]. From a fundamental point of view cross-streamline migration gives
rise to interesting dynamical patterns in particulate flow. Examples are the
bimodal density distribution in dilute polymer suspensions, as we discussed in
chapter 4, characteristic phonon modes of one-dimensional microfluidic crystals
in flat microchannels [27, 28, 29] or collective dynamics in two-dimensional arrays
[30, 31, 32].
To gain a better understanding for the rheology of such suspensions of rigid
particles, their relative motion is central. In 1972 Batchelor and Green [34] showed
that two spheres pass each other in unbounded shear flow under creeping-flow
condition and return to their initial lateral positions after a binary encounter is
completed [34] [see figure 5.1(a)]. This behavior was also observed in experimen-
tal studies [19]. In addition, in bounded shear flow a different class of particle
80
5.1 Introduction
trajectories was observed [35]. Here, wall proximity gives rise to swapping tra-
jectories, where the particles do not pass each other but instead they swap their
lateral positions and leave in the same direction from where they came [see figure
5.1(b)]. At non-zero Reynolds numbers these trajectories have also been observed
in unconfined shear flow [136, 137, 138].
In this chapter we study the relative motion of two identical spherical col-
loids in a pressure driven Poiseuille flow. They are confined between two planar
parallel walls, whose distance is comparable to the particle diameter. Several
methods for evaluating hydrodynamic interactions between colloids, which are
confined between two parallel walls, are used [139, 140, 141, 142, 143]. Here, we
employ the grand mobility matrix which incorporates the multipole expansion of
the hydrodynamic force densities acting on the particle surfaces and the lubri-
cation correction which we introduced in chapter 2. The nonlinear dynamics of
the particle pair in bounded Poiseuille flow gives rise to new classes of trajecto-
ries. Depending on their initial positions, the particles either move on oscillatory
trajectories while staying bound together in the Poiseuille flow or they separate
after an encounter. While the bound states have not been identified before, we
also observe a new class of trajectories in the unbound state. In addition to the
swapping trajectories, we find, what we call, cross-swapping trajectories. Finally,
using results for the two-particle systems, we illustrate the stability of particle
trains. This chapter is based on publication [B].
81
5. NONLINEAR DYNAMICS OF SPHERICAL PARTICLES IN CONFINED
POISEUILLE FLOW
Figure 5.2: Force- and torque-free hard spheres, with radius aand positions ri,
are confined between two parallel planar plates with distance 2W. A Poiseuille
flow is applied.
5.2 Modeling
In this section we describe how we treat hydrodynamic interactions between force-
and torque-free spherical particles in Poiseuille flow by using the grand mobility
matrix for particle configuration bounded by two parallel plates. We first define
our system and then we summarize details of our Stokesian dynamics simulations.
5.2.1 System geometry
We consider the dynamics of, in general, Nforce- and torque-free hard spheres,
each with radius a. They are confined between two parallel planar walls at po-
sitions z=−Wand z=W, which are infinitely extended in the x−yplane
(see figure 5.2). For simplicity, we assume that the particles only move in the
y−zplane at x= 0, which we set initially. Since there are no forces and torques
acting on the particles and due to the translational symmetry along the xdirec-
tion the spheres stay in that plane [55]. This reduces our computational efforts
considerably. Between the walls a pressure driven Poiseuille flow is created
u0(r) = v0"1−z
W2#ey,(5.1)
where again v0denotes the maximum velocity at the centerline (z= 0),eyis the
unit vector along the yaxis and Wthe distance between the centerline and the
walls.
82
5.2 Modeling
5.2.2 Equations of motion
As we already discussed in section 2.5, the translational and rotational veloci-
ties of Nspherical colloids subject to external forces, torques, and an imposed
Poiseuille flow are given by
v=u0(x) + m(x)h(x) + µ(x)F(x),(5.2)
where v= (v,ω)and u0= (u0,ω0)denote the 6N-dimensional generalized N-
particle velocities, F= (F,T)the generalized forces, h= (g0,h1
0,h2
0,h3
0)the
higher order velocity multipole moments,
µ= ¯
µtt
L¯
µtr
L
¯
µrt
L¯
µrr
L!,(5.3)
the mobility matrix and
m= ¯
µtd
L¯
µth
L,1¯
µth
L,2¯
µth
L,3
¯
µrd
L¯
µrh
L,1¯
µrh
L,2¯
µrh
L,3!.(5.4)
the matrix containing the other elements of the grand mobility matrix. The
elements ¯
µab
Lof the grand mobility matrix (2.91) are obtained by a multipole
expansion, truncated after the multipolar order Land supplemented by a lu-
brication correction, as discussed in detail in section 2.5. Since in our case the
particles are force- and torque-free (F=T=0), the equation of motion (5.2)
simplifies to
v=u0(x) + m(x)h(x),(5.5)
which is the basis for our Stokesian dynamics simulation.
5.2.3 System parameters
In order to determine the matrices ML(ilmσ, jl0m0σ)and ZL(ilmσ, jl0m0σ0)(see
section 2.5.3), we closely followed the implementation described in Ref. [55],
where one particle in Poiseuille flow is treated, and generalized it to Npar-
ticles by appropriate modifications (see Appendix A). Since the evaluation of
ML(ilmσ, jl0m0σ)at each time step is very time consuming, we defined in the
microchannel a grid with 160 ×25 grid points and a mesh size of 0.05W. Then,
we stored the values of ML(ilmσ, jl0m0σ)on the grid for different configurations
riand rj, which we define by the axial distance and the two lateral positions
of the two particles i, j. Matrix values at intermediate locations on the grid are
determined from the grid points by linear interpolation. Matrices and vectors are
truncated after the multipolar order L= 5.
83
5. NONLINEAR DYNAMICS OF SPHERICAL PARTICLES IN CONFINED
POISEUILLE FLOW
After transforming to real basis functions and applying the lubrication cor-
rection, which we described in section 2.5.4, we obtain the grand mobility matrix
(2.91). According to equation (5.5) we compute the angular and translational
velocities for a given particle configuration x= (r1,...,rN)and, then, update
the configuration with a small timestep ∆t= 2 ×10−3in units of W/v0(see
section 2.6.2).
5.3 Results
We first investigate the translational and angular velocity of one suspended parti-
cle and compare our results with these of Staben et al. [144]. Then we investigate
the nonlinear dynamics of a system of two equal-sized particles and, finally, we
use our findings to interpret the dynamics of a linear array of colloids.
5.3.1 Translational and angular velocity of a sphere
As we already discussed in section 4.1, at low Reynolds numbers close to zero,
single spherical particles just follow the streamlines in a laminar flow as a result of
the kinematic reversibility of the Stokes equations. In this section we analyze the
translational and angular velocity of one spherical particle with radius a/W = 0.6.
Figure 5.3 shows the translational velocity along the y-direction vyof the
sphere at different lateral positions z. For a spherical particle suspended in a
Poiseuille flow without bounding walls, the translational velocity is simply given
by Faxén’s law (F=0) (see section 2.3.3)
v=Ltu0(r)with Lt=1+a2
6∇2.(5.6)
At the centerline (z= 0), the velocity of the unconfined particle is nearly identical
to the translational velocity of a particle suspended in a confined Poiseuille flow
(see figure 5.3). But due to the interaction with the bounding walls, the velocity
of the confined particle decreases more towards the walls. We observe, that
truncating the multipole expansion after the multipolar order L= 5 and applying
the lubrication correction, the results are in a very good agreement to those of
Staben et al. [144].
Figure 5.4 shows the angular velocity of the sphere. Without bounding walls,
the velocity grows linear with z, which one can easily verify by inserting (5.1)
into Faxén’s law for the rotational motion (2.47). In the presence of the walls
the angular velocity increases monotonically with zuntil it reaches a maximum,
then the velocity decreases towards the walls due to the hydrodynamic interaction
with the walls.
84
5.3 Results
0.00 0.08 0.16 0.24 0.32 0.40
z/W
0.00
0.25
0.50
0.75
1.00
vy/v0
no walls
L= 3
L= 5
L= 5+ LC
Staben et al.
Figure 5.3: Translational velocity along the y-direction at different lateral posi-
tions zof the sphere. Our results for the multipolar orders L= 3 and L= 5 (with
and without lubrication correction [LC]), are compared to those of Staben et al.
[144] and for a particle in Poiseuille flow without bounding walls.
0.00 0.08 0.16 0.24 0.32 0.40
z/W
0.0
0.1
0.2
0.3
0.4
0.5
ωx/[v0/W]
no walls
L= 3
L= 5
L= 5+ LC
Staben et al.
Figure 5.4: Angular velocity at different lateral positions zof the sphere. The
results are shown for the multipolar orders L= 3 and L= 5 (with and without
lubrication correction), for a sphere in unbounded Poiseuille flow and from Ref.
[144].
85
5. NONLINEAR DYNAMICS OF SPHERICAL PARTICLES IN CONFINED
POISEUILLE FLOW
Figure 5.5: Two identical spherical colloids with an initial axial distance y21 =
y2−y1and initial lateral positions z1and z2.
5.3.2 Relative motion of two identical spherical colloids
We analyze the dynamics of two identical spherical colloids for different initial
conditions, namely the axial distance y21 and the lateral positions z1and z2, as
indicated in figure 5.5. We always assume that the leading sphere 2 initially
moves slower than the trailing sphere meaning |z1|<|z2|. So there will always be
a binary encounter of the two colloids. Then, the other case, |z1|>|z2|, always
belongs to trajectories after a binary encounter, as can be directly understood
from the kinematic reversibility of the Stokes equations.
5.3.2.1 Unbound and oscillatory bound states
In the case that the two particles cannot pass each other (2a/W > 1), we identify
different types of trajectories. They are illustrated for a/W = 0.6in figure 5.6
in the center-of-mass frame. Depending on the initial separation (but also the
initial lateral positions), we observe two states of particle motion, an unbound
state and an oscillatory bound state. In the unbound state the spheres approach
each other, turn around, and separate [see figure 5.6 (a) and (b)]. They move
either on swapping trajectories (ST) or cross-swapping trajectories (CS). In 1972,
Batchelor and Green showed that binary encounters between two spherical par-
ticles in unbounded shear flow lead to passing trajectories, where each particle
returns to its initial lateral position after a binary encounter is completed [34].
The swapping trajectories were first observed by Zurita-Gotor et al. for a pair
of spherical colloids in confined shear flow [35]. Here the particles migrate across
the streamlines such that they exchange their lateral coordinates at the end of
the binary encounter [see figure 5.6(a)]. Hence, such particle collisions displace
the spheres across the streamline of the external flow. In addition to this type of
trajectories, we also observe particles moving on cross-swapping trajectories [see
86
5.3 Results
−3−2−10123
−1.0
−0.5
0.0
0.5
1.0
z/W
a)
ST
−3−2−10123
(y−yC)/W
−1.0
−0.5
0.0
0.5
1.0
z/W
c)
8P
−3−2−10123
−1.0
−0.5
0.0
0.5
1.0
b)
CS
−3−2−10123
(y−yC)/W
−1.0
−0.5
0.0
0.5
1.0
d)
OT
Figure 5.6: Trajectories for different initial conditions plotted in the center-of-
mass frame. In the unbound state the spheres move on either (a) swapping trajec-
tories (ST) or (b) cross-swapping trajectories (CS). In the bound state they move
on trajectories (c) forming the pattern of an eight (8P) or (d) with an oval shape
(OT).)
figure 5.6(b)]. In contrast to swapping motion, the colloids approach closer to
each other and they exchange their lateral positions crosswise meaning
z0
1=−z2and z0
2=−z1,(5.7)
where z0
1,z0
2are the lateral positions after the binary encounter.
In addition to these unbound trajectories, where the particles separate at long
times, we can also observe oscillatory bound states when the initial axial distance
y21 is decreased. Here the particles oscillate on their trajectories while staying
together in the Poiseuille flow [see figure 5.6 (c) and (d)]. The spheres move on
trajectories either forming the pattern of an eight (8P) or with an oval shape
(OT) when viewed in the center-of-mass frame. Figure 5.7 shows state diagrams
for particle trajectories when the initial axial distance y21 is fixed. We varied the
the initial lateral positions z1and z2and determined the type of particle trajec-
tories as discussed in figure 5.6. The frequencies of the oscillatory bound states
are color-coded and show interesting structural features. The region of bound
states concentrates around the diagonal z2=−z1. For an initial axial distance of
y21/W = 1.3, where the spheres with radius a/W = 0.6are nearly touching, we
observe a wide region of oscillatory bound states [see figure 5.7(a)]. The oscilla-
tory state of type 8P occurs around z1≈0and z2≈0with smaller oscillation
87
5. NONLINEAR DYNAMICS OF SPHERICAL PARTICLES IN CONFINED
POISEUILLE FLOW
−0.4−0.2 0.0 0.2 0.4
z1/W
0.4
0.2
0.0
−0.2
−0.4
z2/W
8P
ST
ST OT
OT
CS
CS
CS
CS
a)
y21/W = 1.3
−0.4−0.2 0.0 0.2 0.4
z1/W
0.4
0.2
0.0
−0.2
−0.4
z2/W
8P
ST
ST OT
OT
CS
CS
CS
CS
b)
y21/W = 2.0
−0.4−0.2 0.0 0.2 0.4
z1/W
0.4
0.2
0.0
−0.2
−0.4
z2/W
8P
ST
ST
OTOT
CS
CS
CS
CS
c)
y21/W = 3.0
−0.4−0.2 0.0 0.2 0.4
z1/W
0.4
0.2
0.0
−0.2
−0.4
z2/W
ST
CS
CS
d)
y21/W = 4.0
0.000 0.005 0.010 0.0.15 0.020 0.025
ωW/v0
Figure 5.7: State diagram for the different types of trajectories of a two-particle
system. It is plotted when the initial lateral positions z1and z2are varied, while the
initial axial distance y21 is fixed. The particle radius is a/W = 0.6. a) y21/W = 1.3,
b) y21/W = 2.0, c) y21/W = 3.0, and d) y21/W = 4.0. The frequencies ωof the
oscillatory states [oval trajectory (OT) and pattern of eight (8P)] are color-coded.
Outside the colored region, unbound states with swapping- (ST) or cross-swapping
(CS) trajectories occur. The trajectories in the gray shaded areas show binary
encounters of either ST or CS type when they are traced back into the past.
88
5.3 Results
−0.6−0.3 0.0 0.3 0.6
z1/W
0.6
0.3
0.0
−0.3
−0.6
z2/W
8P
ST
ST OT
OT
PT
PT
PT
PT
CS
CS
a)
a/W = 0.4
−0.4−0.2 0.0 0.2 0.4
z1/W
0.4
0.2
0.0
−0.2
−0.4
z2/W
8P
ST
ST OT
OT
CS
CS
CS
CS
b)
a/W = 0.6
−0.2−0.1 0.0 0.1 0.2
z1/W
0.2
0.1
0.0
−0.1
−0.1
z2/W
8P
ST
ST OT
OT
CS
CS
CS
CS
c) a/W = 0.8
0.000 0.005 0.010 0.0.15 0.020 0.025
ωW/v0
Figure 5.8: State diagrams in z1, z2for an initial axial distance y21/W = 2.0and
different particle sizes: a) a/W = 0.4, b) a/W = 0.6, and c) a/W = 0.8.
frequencies than in the OT type. For initial lateral coordinates outside this re-
gion, the particles exhibit the unbound state. They either move on swapping
trajectories for initial coordinates situated around the diagonal z2=z1or cross-
89
5. NONLINEAR DYNAMICS OF SPHERICAL PARTICLES IN CONFINED
POISEUILLE FLOW
swapping trajectories. In the gray areas the leading particle 2 in the Poiseuille
flow has a larger velocity and a binary encounter and thus particle swapping does
not occur. However, tracing the particles into the past, the respective swapping
or cross-swapping trajectories are observed. This is clear either from the kine-
matic reversibility of the Stokes equations or when the direction of the Poiseuille
flow is reversed. Then the leading particle 2 becomes the trailing sphere. The
boundary between white and gray areas are the diagonals z2=z1and z2=−z1.
Interestingly, the region of bound 8P trajectories results from intersecting the
region with oscillatory motion and the region with swapping trajectories. This
makes sense since in the 8P trajectories as well as in the swapping trajectories
both particles cover the same zregion.
Now, when we increase the initial lateral distance y21 of the particles [see fig-
ure 5.7 (a)-(d)], the region of bound oscillatory states becomes narrower and the
frequencies decrease until ultimately the bound states vanish completely. Obvi-
ously, the hydrodynamic coupling between the two particles is no longer strong
enough to hold them together. Simultaneously, we also observe that the region
of swapping trajectories expands with increasing y21.
5.3.2.2 Dependence on particle size
So far we analyzed particle trajectories at a fixed radius of a/W = 0.6of the
two spheres. Figure 5.8 shows state diagrams in z1, z2for different particle sizes
at an initial axial distance of y21/W = 2. Small particles (a/W = 0.4) can pass
each other and in the state diagram of figure 5.8(a) we indeed observe passing
trajectories (PT) in addition to the trajectories reported so far.
Such passing trajectories were first observed in unbounded [34] and then
bounded [35] shear flow. The particles pass each other during approach, separate
again, and return to their initial cross-streamline positions. In figure 5.9 we show
several particle trajectories. Particle 1 starts at z1/W = 0.05 while the lateral
position z2of particle 2 is gradually decreased so that the transition from passing
to cross-swapping trajectories takes place. In contrast to the trajectories in linear
shear-flow [35], the trajectories of both particles are not symmetric since the shear
rate in Poiseuille flow varies across the streamlines. In the passing trajectories the
trailing sphere 1 has a larger velocity than the leading sphere 2. The transition to
cross-swapping takes place when during approach the trailing sphere migrates to
a position with smaller flow velocity than the leading sphere, v0y(z1)< v0y(z2), as
indicated in figure 5.9. As a consequence, the particles separate before the pass-
ing trajectories can be completed and they move on cross-swapping trajectories.
This behavior is a feature of the Poiseuille flow and cannot be observed in linear
shear flow.
For larger particle sizes the passing trajectories vanish and the oscillatory
90
5.3 Results
−1.0−0.5 0.0 0.5 1.0
(y−yC))/W
−1.0
−0.5
0.0
0.5
1.0
z/W
v0y(z1)< v0y(z2)
Particle 1
Particle 2
Figure 5.9: Passing trajectories of two identical spheres with radius a/W = 0.4
turn into cross-swapping trajectories (plotted in the center-of-mass frame): sphere
1 starts at a lateral position of z1/W = 0.005 while the lateral coordinate z2of
sphere 2 is gradually decreased. The transition to cross-swapping occurs when
v0y(z1)< v0y(z2).
state occupies a larger area in the accessible region of the z1, z2state diagram, as
figure 5.8 (b) and (c) demonstrate. So, when the confinement becomes stronger,
the oscillatory bound state becomes more important. In figure 5.10 we plot
the oscillation frequencies in the bound state against the particle radius afor
different initial axial distances y21. In plot a) we fix the initial lateral positions to
z1/W = 0.1and z2/W = 0. Small spheres then lie outside the narrow region of the
bound states [see figure 5.8(a)] and do not oscillate (ω= 0). When the particles
size increases, this region expands and at a critical radius acrit. the particle system
enters the regime of the oscillating 8P trajectories. The oscillation frequency
increases from zero with increasing a. Larger initial distances y21 reduce the
hydrodynamic coupling. Therefore the spheres have to be larger to initiate the
bound state and the frequency is smaller. Increasing the particle separation even
further, would require particle radii a/W > 1. Since this cannot be, oscillatory
states do not occur anymore.
91
5. NONLINEAR DYNAMICS OF SPHERICAL PARTICLES IN CONFINED
POISEUILLE FLOW
0.2 0.3 0.4 0.5 0.6
a/W
0.000
0.001
0.002
0.003
0.004
0.005
0.006
ω/[v0/W]
a) b)
z1/W = 0.1,z2/W = 0
y21/W
1.4
1.6
1.8
1.9
2.0
2.1
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
a/W
0.000
0.005
0.010
0.015
0.020
0.025
ω/[v0/W]
y21 4.2
>
≈
z1/W = 0.3,z2/W =−0.3
y21/W
1.4
1.8
2.2
2.6
3.0
3.4
Figure 5.10: Oscillation frequency ωin the bound states plotted versus particle
radius afor different initial axial distances y21. The lateral coordinates are set to
a) z1/W = 0,z2/W = 0.1and b) z1/W = 0.3,z2/W =−0.3.
In figure 5.10(b) we choose initial lateral coordinates z1/W = 0.3and z2/W =
−0.3. They lie on the diagonal in our z1, z2state diagrams where the spheres move
on oscillatory trajectories of type OT. Obviously, decreasing the particle radius
to zero at fixed y21, the width of the oscillatory region goes to zero. This explains
that the oscillation frequency also continuously approaches zero with a. With
increasing axial distance y21, the whole curve rises less and there is a limiting
axial distance y21 ≈4.2where the bound state vanishes completely.
5.3.3 Linear array of particles
Particles in bound states do not pass each other and do no separate as they do
when moving on passing, swapping, or cross-swapping trajectories. This may
have significant consequences for microfluidic systems, where precise manipula-
tion of linear arrays of, e.g, nearly aligned particles is very important [8, 145]
The bound two-particle states now suggest that particles in such particle trains
do not aggregate when they move with slightly different velocities on different
streamlines. Instead, they stay separated while moving on oscillatory trajecto-
ries. To investigate this question, we implemented linear arrays of particles as
shown in Figs. 5.11 and 5.12 using approximate periodic boundary conditions.
We extend the particle train to the right and to the left by a respective group
of three particles which we copy to the other side of the microchannel after each
time step in our simulations. In Figs. 5.11 and 5.12 we illustrate the dynamics
92
5.3 Results
of the particle trains.
In figure 5.11(a) we first investigate a particle train where the particles of
radii a/W = 0.6are approximately arranged along the centerline. Each sphere
starts with a small random lateral offset ∆zi/W ∈[−0.05,0.05] relative to the
centerline at z= 0 which may mimic experimental conditions. In addition, the
initial axial distance between two neighboring particles is set to yi+1,i/W = 1.6
to ensure that they belong to the two-particle oscillatory bound state of type
8P. Indeed, the snapshots in figure 5.11(a) show that the particle train is stable.
Between the snapshots the particles travel several times through the simulated
piece of microchannel. The particles do not aggregate and each particle moves
on an approximately oscillatory trajectory. It oscillates around the centerline
as illustrated at the bottom of figure 5.11(a), where we plot the zcoordinate
of the blue particle versus time. It also oscillates between its two neighbors as
the ycoordinate relative to the center-of-mass coordinate yCshows but with a
much smaller frequency compared to the lateral direction. A closer inspection
reveals slightly irregular oscillations which we attribute to the two neighboring
particles. In total, these oscillations demonstrate that the particle train does not
disintegrate.
In figure 5.11(b) we illustrate the dynamics of a particle train when we initially
arrange the spheres in a zig-zag pattern. The lateral positions of neighboring
particles approximately fulfill zi+1 =−zito which we also add a small random
offset ∆zi/W ∈[−0.05,0.05]. Together with the initial axial distance of two
neighboring particles, yi+1,i/W = 1.6, we expect that each particle moves on a
trajectory of type OT as in the two-particle bound state. Indeed, the snapshots
again reveal a stable particle train. The lateral and relative axial displacements
are less regular than in a). However, they are both bound so that each particle
in the center-of-mass frame moves on a restricted area similar to the two-particle
trajectories of type OT. As soon as we disturb the regular particle trains too
much, they are no longer stable. We show this in figure 5.12 for the two cases
of the linear (a) and the zig-zag (b) pattern. The blue particle starts with a
sufficiently large vertical offset so that with one of its neighbors it belongs to a
two-particle unbound state. Hence, all the particles no longer perform oscillatory
trajectories. They approach each other and migrate across streamlines until they
arrange in a zig-zag order. Here they constantly form transient clusters which
then dissolve in time at the downstream end of the cluster. The completely
irregular particle motion also becomes clear from the lateral position and the
relative axial displacement of the blue particle plotted versus time at the bottom
of Figs. 5.12a) and b).
93
5. NONLINEAR DYNAMICS OF SPHERICAL PARTICLES IN CONFINED
POISEUILLE FLOW
a)
t= 0W/v0
133W/v0
267W/v0
400W/v0
0 200 400 600 800 1000
−0.4
0.0
0.4
z/W
0 200 400 600 800 1000
t/[W/v0]
−2.3
−0.8
0.7
(y−yC)/W
b)
t= 0W/v0
133W/v0
267W/v0
400W/v0
0 200 400 600 800 1000
−0.4
0.0
0.4
z/W
0 200 400 600 800 1000
t/[W/v0]
−3.9
−2.4
−0.9
(y−yC)/W
Figure 5.11: Snapshots of a linear array of particles with radius a/W = 0.6
at different times t. In addition, the lateral coordinate zand the relative axial
displacement y−yCin the center-of-mass frame is plotted versus time for the blue
colored particle. a) Initially, the spheres are placed close to the centerline with a
small random offset ∆zirelative to z= 0. b) Initially, the spheres are arranged in
a zig-zag pattern, where the particles have a small random offset ∆zirelative to
zj+1 =−zj=±0.25W.
94
5.3 Results
a)
t= 0W/v0
133W/v0
267W/v0
400W/v0
0 200 400 600 800 1000
−0.5
0.0
0.5
z/W
0 200 400 600 800 1000
t/[W/v0]
−4.8
−0.8
3.2
(y−yC)/W
b)
t= 0W/v0
133W/v0
267W/v0
400W/v0
0 200 400 600 800 1000
−0.5
0.0
0.5
z/W
0 200 400 600 800 1000
t/[W/v0]
−6.4
−2.4
1.6
(y−yC)/W
Figure 5.12: Same setting as in figure 5.11. Now the offset of the blue colored
particle is larger so that neighboring particles belong to the unbound state.
95
5. NONLINEAR DYNAMICS OF SPHERICAL PARTICLES IN CONFINED
POISEUILLE FLOW
5.4 Summary and conclusion II
We identified and analyzed the different classes of trajectories for two identical
spherical colloids subject to Poiseuille flow under creeping-flow condition. The
distance of the bounding walls was comparable to the particle diameter. In ad-
dition, we investigated the stability of particle trains. To treat hydrodynamic in-
teractions between the particles, we determined the grand mobility matrix (2.91)
using a multipole expansion for the hydrodynamic force densities induced on the
particle surfaces. To be sufficiently accurate, we included multipoles up to the
fifth order and applied a lubrication correction to the friction matrix.
Two identical particles move either in an unbound or an oscillatory bound
state depending on their initial separation and lateral positions in the microchan-
nel. In the unbound state the spheres follow swapping-, cross-swapping or even
passing trajectories when the particles are sufficiently small. Whereas the passing
and swapping trajectories were already observed in bounded linear shear flow [35],
the cross-swapping trajectories are a distinct feature of the bounded Poiseuille
flow. Here, the particles move first on passing trajectories until the trailing sphere
migrates to a lateral position with lower velocity compared to the leading par-
ticle. As a consequence, the particles separate before they can pass each other.
For initial small axial distances y21 even bound states exist which do not occur in
linear shear flow. Here, the particles follow oscillatory trajectories while staying
together in the Poiseuille flow. Depending on the initial lateral positions, they
move on trajectories tracing either an eight or an oval shape. To categorize the
different states of motion, we provided state diagrams in the lateral positions for
different initial axial distances and also analyzed how the diagrams changed with
the particle radii. In addition, we determined the oscillation frequencies in the
bound states.
Particles moving on bound oscillatory trajectories do not separate which may
have significant consequences for microfluidic systems. The bound two-particles
states can explain why particle trains that either are nearly straight or have a
zig-zag pattern are stable. Defects in these ordered particle trains correspond
to particles with distinctly different flow velocities so that the particle trains
become unstable. They disintegrate and aggregate to transient clusters. Thus,
a proper understanding of the migration pattern in the two-particle bound state
may contribute towards finding better methods of controlling particle motion in
microfluidic devices.
96
6
Conclusions
«In science one tries to tell peo-
ple, in such a way as to be under-
stood by everyone, something that
no one ever knew before. But in
the case of poetry, it’s the exact
opposite!»
Paul Adrien Maurice Dirac
97
6. CONCLUSIONS
The aim of this work was to explore the potential of hydrodynamic interactions
to induce novel motional patterns in confined flow, which then provides a deeper
understanding for relevant physical and biological processes, such as the cross
streamline migration of semiflexible polymers [37] or the collective dynamics of
dense particle arrays [31].
We started by introducing the theoretical framework of this work. In chap-
ter 2, we presented the hydrodynamic equations, which are the basis to describe
the hydrodynamic interactions of suspended particles. Then we presented con-
cepts to describe these interactions between spherical particles suspended in an
unbounded fluid as well as in a bounded fluid under the influence of a pressure
driven Poiseuille flow. In chapter 3, we introduced the basic concept of Brownian
motion. First we presented the Langevin equation, a stochastic differential equa-
tion describing the time evolution of the suspended particles. Then we introduced
an equivalent approach to consider the probability distribution of the particles,
which is given by the Smoluchowski equation.
In chapter 4, we investigated the dynamics of a bead-spring chain confined
between two parallel walls under the influence of a pressure driven Poiseuille
flow. The hydrodynamic interactions between the point-like beads of the model
polymer were taken into account by the two-wall Green tensor of the Stokes
equations [55]. Performing Brownian dynamics simulations, we showed that the
characteristic bimodal probability distribution for the polymer’s center of mass
is governed by the hydrodynamic interactions along the bead-spring chain and
with the bounding walls. Without these interactions, the center of mass is dis-
tributed homogeneously across the channel. In addition, we presented numerical
results on the orientational order of the end-to-end vector of the model polymer
together with analytical hard-needle expressions at zero flow velocity. To gain
a better understanding of the migration mechanism of the polymer, we derived
a Smoluchowski equation for the center-of-mass distribution and carefully ana-
lyzed the different contributions to the probability current. In contrast to Refs.
[37, 38, 52], where the migration away from the centerline was explained by a
position-dependent diffusivity, we clearly identified a deterministic drift current
as the major cause for migration away from the centerline. We even showed that
diffusional currents due to a position-dependent diffusivity became less important
with increasing polymer stiffness.
Therefore, with our results, we provide a full understanding of the migration
mechanism of semiflexible polymers in confining geometries. In particular, we
clarified how the polymers migrate away from the centerline. Thereby we are
now able to describe the observed bimodal concentration profiles completely.
In chapter 5 we investigated the relative motion of a pair of force and torque-
free spherical particles in Poiseuille flow. Since the particle’s diameter was com-
parable to the channel width, we could no longer work with the simplified point
98
particles but had to take into account their finite extent. We determined the
grand mobility matrix using a multipole expansion of the force densities on the
spheres supplemented by a lubrication correction. Performing Stokesian dynam-
ics simulations, we showed, that the bounded Poiseuille flow gives rise to new
classes of trajectories resulting in cross-streamline migration of the two spheres.
Particles moving on these new trajectories exhibited either bound or unbound
states. In the unbound state the spheres followed swapping-, cross-swapping or
even passing trajectories. Whereas the passing and swapping trajectories were al-
ready observed in bounded linear shear flow [35], the cross-swapping trajectories
are a distinct feature of the bounded Poiseuille flow. Here, the particles moved
first on passing trajectories until the trailing sphere migrated to a lateral posi-
tion with lower velocity compared to the leading particle. As a consequence, the
particles separated before they can pass each other. In the bound state, which do
not occur in linear shear flow, the particles followed oscillatory trajectories while
staying together in the Poiseuille flow. Depending on the initial lateral positions,
they moved on trajectories tracing either an eight or an oval shape.
We analyzed the different classes of trajectories depending on the initial po-
sitions of the two particles and their size, and presented state diagrams, where
we categorized the trajectories and color coded the oscillation frequencies of the
bound states. Finally, we discussed how the results on the two-particle system
help to understand the stability of particle trains composed of several particles.
With help of our results we offer a deeper understanding of the cross-streamline
migration of colloidal particles in confined geometries under the influence of a
Poiseuille flow. This can have potential consequences for microfluidic devices.
For example, we showed that the center-of-mass distribution of semiflexible poly-
mers strongly depend on the stiffness of the filament, which may offer potential
applications in separation of different species of biopolymers. We also understood
the stability of linear particle trains, which is important for microfluidic systems,
where a precise manipulation of linear arrays of particles is necessary.
A possible extension of this work is to include external forces and torques
in our system. For example, one can think to add permanent dipole moments
to the colloids and use an electric or magnetic field to control their directions.
As a possible next step, it would be interesting to investigate the dynamics of
suspensions of different species (e.g. different in size) under confined flow. Also,
a very challenging task would be to change the system geometry, which may be
realized by fixed colloids. For example, one could add spherical obstacles inside
the channel or simulate a constriction of the channel by a linear array of fixed
particles.
To summarize: We investigated the nonlinear dynamics in two different col-
99
6. CONCLUSIONS
loidal model systems where the particles are confined between two planar walls
and driven out of equilibrium by an imposed Poiseuille flow. We showed, that
the bounded Poiseuille flow gives rise to interesting dynamic pattern formations.
On the one hand, our results reproduced the characteristic bimodal distribu-
tion of the the center of mass of a semiflexible polymer. For the first time, we
were able to identify the mechanism for the migration away from the centerline.
Therefore, our results provides a full understanding of the migration mechanism
of semiflexible polymers in confining geometries, which may offer potential appli-
cations in separation of biopolymers.
On the other hand, we showed that the bounded Poiseuille flow gives rise
to new and interesting classes of trajectories for a pair of force and torque-free
particles. With the help of our results, we were able to understand the stability
of particle trains composed of several particles, which may have significant con-
sequences for microfluidic systems, where precise manipulation of linear arrays of
particles is very important.
100
A
Appendix
In order to determine the matrices ML(ilmσ, jl0m0σ)and ZL(ilmσ, jl0m0σ0)(see
section 2.5.3), we closely follow the implementation described in Ref. [55], where
one particle in Poiseuille flow is treated, and generalize it to Nparticles by
appropriate modifications.
In section 2.5.1 we showed that the two-wall Green tensor splits up into the
Oseen tensor and a reflectional part. The reflectional part is given by
TR(r,r0) = Zeiq·sˆ
TR(q, z, z0)e−iq·s0d2q , (A.1)
where s= (x, y)and s0= (x0, y0)are the two dimensional space vectors and
q= (qx, qy)is the wave vector. In Fourier space, the reflectional part reads
ˆ
TR(q, z, z0) = 1
16πηq [t1nnezez+it1npezˆq+it1pn ˆqez
−t1pp ˆqˆq+r1pp(1−ezez)] ,(A.2)
where the scalar functions are given by
t1.. =t1..(q, z, ¯z) = F(w)TE..(u)F(v),
r1pp =r1pp(q, z, ¯z) = F(w)TGpp(u)F(v).(A.3)
Here, we have introduced the scalars u=qW,v=qz0and w=qz, the vec-
tor F(v)=(f1(v), f2(v), f3(v), f4(v)) = (cosh v, sinh v, v cosh v, v sinh v)and the
101
A. APPENDIX
matrices
Epp(u) =
E−C−−2utanh uE−0 0 −E−A−
0E+D−−2ucoth uE+−E+B−0
0−E+B−E+0
−E−A−0 0 E−
,
Gpp(u) =
−2e−u
cosh u0 0 0
0−2e−u
sinh u0 0
0 0 0 0
0 0 0 0
,
Enn(u) =
−E+C+0 0 E+A+
0−E−D+E−B+0
0E−B+−E−0
E+A+0 0 −E+
,
Enp(u) =
0u2E+−E+A+0
u2E−0 0 −E−B+
−E−A−0 0 E−
0−E+B−E+0
,
Epn(u) = −ET
np(u),(A.4)
with
A±=A±(u) = u±sinh ue−u,
B±=B±(u) = u±cosh ue−u,
C±=C±(u) = u(1 + u)±sinh ue−u,
D±=D±(u) = u(1 + u)±cosh ue−u,
E±=E±(u) = 1
sinh ucosh u±u.(A.5)
With the help of the adjoint set of functions
w+
lm0=1
l(2l+ 1)r−l Alm −2l+ 3
2Blm!,(A.6)
w+
lm1=i
l(l+ 1)r−l−1Clm,(A.7)
w+
lm2=1
(l+ 1)(2l+ 1)r−l−2Blm,(A.8)
102
with
Alm =lYlmer+∂Ylm
∂θ eθ+1
sin θ
∂Ylm
∂φ eφ,(A.9)
Blm =−(l+ 1)Ylmer+∂Ylm
∂θ eθ+1
sin θ
∂Ylm
∂φ eφ,(A.10)
Clm =1
sin θ
∂Ylm
∂φ eθ−∂Ylm
∂θ eφ(A.11)
where Ylm are the spherical harmonics, we can calculate the matrix elements (see
section 2.5.3)
TR(ilmσ, jl0m0σ0) = Dw+
lmσ(i)δiTR(r,r0)w+
l0m0σ0(j)δjE
=Zw∗+
lmσ(r−ri)1
aδ(|r−ri|−a)·. . .
. . . ·ZTR(r,r0)·w+
l0m0σ0(r0−rj)1
aδ(|r0−rj|−a)d3r0d3r
=Zw∗+
lmσ(¯
r)1
aδ(¯r−a)·. . .
. . . ·ZTR(¯
r+ri,¯
r0+rj)·w+
l0m0σ0(¯
r0)1
aδ(¯r0−a)d3¯r0d3¯r
(A.12)
with ¯
r=r−riand ¯
r0=r0−rj. Inserting equation A.1, we obtain
TR(ilmσ, jl0m0σ0) = Zw∗+
lmσ(¯
r)1
aδ(r0−a)·. . .
. . . ·ZZeiq·(¯
s+si)ˆ
TR(q,¯z+zi,¯z0+zj)e−iq·(¯
s0+sj)d2q·. . .
. . . ·w+
l0m0σ0(¯
r0)1
aδ(¯r0−a)d3¯r0d3¯r. (A.13)
Defining v=qz0=q¯z0+qzj=: ¯v0
j+ρjand w=qz =q¯z+qzi=: ¯wi+ρiwith
ρi/j =qzi/j, we can express fi(v)in terms of fi(¯v0
j)by a matrix operation
F(v) = M(ρj)F(¯v0
j),
F(w) = M(ρi)F( ¯wi)(A.14)
with
M(ρ) =
f1(ρ)f2(ρ) 0 0
f2(ρ)f1(ρ) 0 0
f3(ρ)f4(ρ)f1(ρ)f2(ρ)
f4(ρ)f3(ρ)f2(ρ)f1(ρ)
.(A.15)
103
A. APPENDIX
Thus from equations (A.2), (A.3) and (A.14) we find that ˆ
TRis a sum of terms,
each of which is separable in the variables ¯v0
jand ¯wi, while the ¯
sand ¯
s0depen-
dence in equation (A.13) is already separable, so the calculation of the matrix
element reduces to a product of integrals, one over ¯
rand one over ¯
r0. Explicit
results for the integrals over ¯
rand ¯
r0in equation (A.13) are given in Ref. [55].
Therefore, to determine the matrix elements TR(ilmσ, jl0m0σ), we have to solve
the remaining integral over qin equation (A.13), which is calculated by quadra-
ture. Note, that for i=jequation A.13 is identical to the equation given in Ref.
[55].
In order to obtain the whole matrix M (see section 2.5.3), we additional need
the matrix elements M0(ilmσ, il0m0σ)and M1(ilmσ, jl0m0σ)(i6=j), which are
given in Ref. [56] and [146], respectively.
104
List of Publications
[A] Sebastian Reddig & Holger Stark, Cross-streamline migration of a semiflexible
polymer in a pressure driven flow, Journal of Chemical Physics, 135, 165101
(2011)
[B] Sebastian Reddig & Holger Stark, Nonlinear dynamics of spherical particles
in Poiseuille flow under creeping-flow condition, Journal of Chemical Physics,
138, 234902 (2013)
105
106
References
[1] D. N. Breslauer, P. J. Lee, and L. P. Lee.Mol. BioSyst.,2:97, 2006.
[2] S.K. Sia and G.M. Whitesides.Electrophoresis,24:3563, 2003.
[3] M. B. Romanowski, A. R. Abate, A. Rotem, Holtze C., and D. A.
Weitz.Lab Chip,12:802, 2012.
[4] H. Stone, A. Stroock, and A. Ajdari.Annu. Rev. Fluid Mech.,
36:381, 2004.
[5] S. Köster, D. Steinhauser, and T. Pfohl.J. Phys.: Condens. Mat-
ter,17:4091, 2005.
[6] H. Yun, K. Kim, and W. G. Lee.Biofabrication,5(2):022001, 2013.
[7] R. D. Sochol, B. P. Casavant, M. E. Dueck, L. P. Lee, and
L. Lin.J. Micromech. Microeng.,21:054019, 2011.
[8] G. F. Christopher and Anna S. L. J. Phys. D: Appl. Phys.,40:R319,
2007.
[9] J. L. M. Poiseuille.Ann. Sci. Nat.,5:111, 11836.
[10] G. Segré and A. Silberberg.J. Fluid Mech.,14:115, 1962.
[11] R. Fåhræus and T. Lindqvist.Amer. J. Physiol.,96:562, 1931.
[12] S. S. Kuntaegowdanahalli, A. A. S. Bhagat, and I. Papautsky.
Lab on a Chip,9:2973, 2009.
107
REFERENCES
[13] D. Di Carlo, J. F. Edd, D Irimia, and M. Tompkins.Anal. Chem.,
80:2204, 2008.
[14] Z. Wu, B. Willing, J. Bjerketorp, J. Jansson, and K. Hjort.
Lab on a Chip,9:1193, 2009.
[15] R. C. Jeffrey and J. R. A. Pearson.J. Fluid Mech.,22:721, 1965.
[16] T. W. Pan and R. Glowinski.J. Comput. Phys.,181:260, 2002.
[17] A. Karnis, H. L. Goldsmith, and S. G. Mason.Can. J. Chem. Eng.,
44:181, 1966.
[18] D. Di Carlo, D. Irimia, R. G. Tompkins, and M. Toner.Proc.
Natl. Acad. Sci.,104:18892, 2007.
[19] C. L. Darabaner and S. G. Mason.Rheol. Acta,6:273, 1967.
[20] E. S. Asmolov.J. Fluid Mech.,381:63, 1999.
[21] J. A. Schonberg and E. J. Hinch.J. Fluid Mech.,203:517, 1989.
[22] L. G. Leal.Ann. Rev. Fluid Mech.,12:435, 1980.
[23] C. Prohm, M. Gierlak, and H. Stark.Phys. J. E,35:80, 2012.
[24] J. K. G. Dhont.An introduction to dynamics of colloids. Elsevier Science,
Amsterdam, 1996.
[25] H. L. Goldsmith and S. G. Mason.The Microrheology of Dispersions
in Rheology Theory and Applications,4. Academic Press, New York, 1967.
[26] D. Semwogerere, J. F. Morris, and E. R. Weeks.J. Fluid. Mech.,
581:437, 2007.
[27] T. Beatus, T. Tlusty, and R. Bar-Ziv.Nature Physics,2:743, 2006.
[28] P. Garstecki.Nature Phys.,2:733, 2006.
[29] T. Beatus, R. Bar-Ziv, and T. Tlusty.Phys. Rev. Lett.,99:124502,
2007.
[30] M. Baron, J. Bĺawzdziewicz, and E. Wajnryb.Phys. Rev. Lett.,
100:174502, 2008.
[31] J. Bławzdziewicz and E. Wajnryb.J. Phys.: Conf.Ser.,392:012008,
2012.
108
REFERENCES
[32] S. Bhattacharya.J. Chem. Phys.,128:074709, 2008.
[33] M. Zurita-Gotor, J. Bławzdziewicz, and E. Wajnryb.Phys. Rev.
Lett.,108:068301, 2012.
[34] G. K. Batchelor and J. T. Green.J. Fluid Mech.,56:375, 1972.
[35] M. Zurita-Gotor, J. Bławzdziewicz, and E. Wajnryb.J. Fluid
Mech.,592:447, 2007.
[36] C. W. Hsu and Y. L. Chen.J. Chem. Phys.,133:034906, 2010.
[37] D. Steinhauser, S. Köster, and T. Pfohl.ACS Macro Letters,1:541,
2012.
[38] R. Chelakkot, R. G. Winkler, and G. Gompper.Europhys. Lett.,
91:14001, 2010.
[39] A. M. Slowicka, E. Wajnryb, and M. L. Ekiel-Jezewska.Eur.
Phys. J. E,36:31, 2013.
[40] K. Sadlej, E. Wajnryb, and M. L. Ekiel-Jezewska.J. Chem. Phys.,
133:054901, 2010.
[41] O. B. Usta, D. Perchak, A. Clarke, J. M. Yeomans, and A. C.
Balazs.J. Chem. Phys.,130:234905, 2010.
[42] H. Noguchi and G. Gompper.Proc. Natl. Acad. Sci.,102:14159, 2005.
[43] G. Danker, P. Vlahovska, and C. Misbah.Phys. Rev. Lett.,
102:148102, 2009.
[44] R. M. Jendrejack, E. T. Dimalanta, D. C. Schwartz, M. D. Gra-
ham, and J. J. De Pablo.Phys. Rev. Lett.,91:38102, 2003.
[45] R. M. Jendrejack, D. C. Schwartz, J. J. De Pablo, and M. D.
Graham.J. Chem. Phys.,120:2513, 2004.
[46] O. B. Usta, J. E. Butler, and A. J. C. Ladd.Phys. Fluids,
18:031703, 2006.
[47] O. B. Usta, J. E. Butler, and A. J. C. Ladd.Phys. Rev. Lett.,
98:98301, 2007.
[48] J. P. Hernández-Ortiz, H. Ma, J. J. De Pablo, and M. D. Gra-
ham.Phys. Fluids,18:123101, 2006.
109
REFERENCES
[49] D. Stein, F. H. J. Van Der Heyden, W. J. A. Koopmans, and
C. Dekker.Proc. Natl. Acad. Sci.,103:15853, 2006.
[50] L. Fang, H. Hu, and R. G. Larson.J. Rheol.,49:127, 2005.
[51] A. Balducci, C. C. Hsieh, and P. S. Doyle.Phys. Rev. Lett.,
99:238102, 2007.
[52] R. Chelakkot, R. G. Winkler, and G. Gompper.J. Phys.: Con-
dens. Matter,23:184117, 2011.
[53] H. Ma and M. D. Graham.Phys. Fluids,17:083103, 2005.
[54] C. Sendner and R. R. Netz.Europhys. Lett.,79:58004, 2007.
[55] R. B. Jones.J. Chem. Phys.,121:483, 2004.
[56] B. Cichocki, R. B. Jones, R. Kutteh, and E. Wajnryb.J. Chem.
Phys.,112:2548, 2000.
[57] A. J. C. Ladd.Phys. Rev. Lett.,70:1339, 1993.
[58] S. Chen and G. D. Doolen.Annu. Rev. Fluid Mech.,30:329, 1998.
[59] P. J. Hoogerbrugge and J. M. V. A. Koelman.Europhys. Lett.,
19:155, 1992.
[60] R. D. Groot and P. D. Warren.J. Chem. Phys.,107:4423, 1997.
[61] A. Malevanets and R. Kapral.J. Chem. Phys.,110:8605, 1999.
[62] R. Kapral.Advances in Chemical Physics,140:89, 2008.
[63] T. Ihle and D. M. Kroll.Phys. Rev. E,67:066705, 2003.
[64] J. Happel and H. Brenner.Low Reynolds Number Hydrodynamics.
Kluwer, Dordrecht, 1983.
[65] J. F. Brady and G. Bossis.Annu. Rev. Fluid Mech.,20:111, 1988.
[66] G. K. Batchelor.An Introduction to Fluid Dynamics. Cambridge Uni-
versity Press, Cambridge, 2000.
[67] L. D. Landau and E. M. Lifshitz.Fluid Mechanics. Pergamon Press,
London, 1959.
110
REFERENCES
[68] S. Vogel.Life in Moving Fluids: The Physical Biology of Flow. Princeton
University Pres, Princeton, 1994.
[69] S. Kim and S. J. Karrila.Microhydrodynamics: Principles and Selected
Applications. Butterworth-Heinemann, Boston, 1991.
[70] H. Faxén.Annalen der Physik,373:89, 1922.
[71] J. Rotne and S. Prager.J. Chem. Phys.,50:4831, 1969.
[72] J.M.J. den Toondor and P.R. Onck.Artificial Cilia. RSC Publishing,
Cambridge, 1994.
[73] H. Yamakawa.J. Chem. Phys.,53:436, 1970.
[74] S. Köster.Biological matter in microfluidic environment - from single
molecules to self-assembly. PhD thesis, Göttingen, 2006.
[75] M. Sharan and A. S. Popel.Biorheology,38:415, 2001.
[76] R. Skalak and P. I. Branemark.Science,164:717, 1969.
[77] H. Míguez, S. M. Yang, and G. A. Ozin.Langmuir,19:3479, 2003.
[78] H. N. Unni and C. Yang.J. Colloid Interface Sci.,291:28, 2005.
[79] E. Kumacheva, P. Garstecki, H. Wu, and G. M. Whitesides.
Phys. Rev. Lett.,91:128301, 2003.
[80] J. R. Blake and A. T. Chwang.J. Eng. Math.,8:23, 1974.
[81] N. Liron and S. Mochon.J. Eng. Math.,10:287, 1976.
[82] C. Chang and R. L. Powell.Part. Sci. Technol.,12:189, 1994.
[83] S. Weinbaum, P. Ganatos, and Z. Yan.Annu. Rev. Fluid Mech.,
22:275, 1990.
[84] B. Cichocki, B. U. Felderhof, K. Hinsen, E. Wajnryb, and
J. Bławzdziewicz.J. Chem. Phys.,100:3780, 1994.
[85] S. Bhattacharya, J. Bławzdziewicz, and E. Wajnryb.Phys. Flu-
ids.,18:053301, 2006.
[86] R. Schmitz and B. U. Felderhof.Physica A,113:90, 1982.
[87] H. Lamb.Hydrodynamics. Cambridge University Press, Cambridge, 1932.
111
REFERENCES
[88] G. Perkins and R. B. Jones.Physica A,171:575, 1991.
[89] B. Cichocki, B. U. Felderhof, and R. Schmitz.PhysicoChem. Hy-
drodyn.,10:383, 1988.
[90] R. Schmitz and B. U. Felderhof.Physica A,116:163, 1982.
[91] B. Cichocki, M. L. Ekiel-Jezewska, and J. Wajnryb.J. Chem.
Phys.,111:3265, 1999.
[92] D. J. Jeffrey and Y. Onishi.J. Fluid Mech.,139:261, 1984.
[93] D. J. Jeffrey.Phys. Fluids A,4:16, 1992.
[94] G. Bossis and J. F. Brady.J. Chem. Phys.,80:5141, 1984.
[95] B. Cichocki and R. B. Jones.Physica A,258:273, 1998.
[96] G. Wanner and E. Hairer.Solving Ordinary Differential Equations II.
Springer-Verlag, Berlin, 1991.
[97] W. H. Press, S. A. Flannery, and W. T. Vetterling.Numerical
Recipes in C: The Art of Scientific Computing. Cambridge University Press,
Cambridge, 1992.
[98] G. Sewell.The Numerical Solution of Ordinary and Partial Differential
Equations. Academic Press, San Diego, 1988.
[99] D. F. Griffiths and D. J. Higham.Numerical Methods for Ordinary
Differential Equations. Springer-Verlag, Berlin, 2010.
[100] I. Karatzas and S. E. Shreve.Brownian Motion and Stochastic Cal-
culus. Springer-Verlag, Berlin, 1991.
[101] J. Argyris, G. Faust, M. Haase, and R. Friedrich.Die Erforschung
des Chaos: Eine Einführung in die Theorie nichtlinearer Systeme. Springer-
Verlag, Berlin, 2010.
[102] D. S. Lemons.An introduction to stochastic processes in physics: contain-
ing "On the theory of Brownian motion" by Paul Langevin. Johns Hopkins
University Press, Baltimore, 2002.
[103] R. F. Fox and G. E. Uhlenbeck.Phys. Fluids,13:1893, 1970.
[104] P. Mazur.Physica A,110:128, 1982.
112
REFERENCES
[105] H. C. Öttinger.Stochastic Processes in Polymeric Fluids. Springer-
Verlag, Berlin, 1996.
[106] W. Paul and J. Baschnagel.Stochastic Processes: From Physiscs to
Finance. Springer-Verlag, Berlin, 1999.
[107] C. W. Gardiner.Handbook of Stochastic Methods: For Physics, Chem-
istry, and the Natural Sciences. Springer-Verlag, Berlin, 1986.
[108] W. T. Coffey, Y. P. Kalmykov, and J. T. Waldron.The Langevin
Equation With Applications to Stochastic Problems in Physics, Chemistry
and Electrical Engineering. World Scientific Publishing Co. Pte. Ltd., Sin-
gapore, 2004.
[109] U. Hassler.Stochastische Integration Und Zeitreihenmodellierung:
Eine Einführung mit Anwendungen aus Finanzierung und Ökonometrie.
Springer-Verlag, Berlin, 2007.
[110] H. Risken.The Fokker-Planck Equation: Methods of Solution and Appli-
cations. Springer-Verlag, Berlin, 1996.
[111] W. Hess and R. Klein.Physica A,94:71, 1978.
[112] M. Fixman.J. Chem. Phys.,69:1527, 1978.
[113] T. J. Murphy and J. L. Aguirre.J. Chem. Phys.,57:2098, 1972.
[114] M. Reichert.Hydrodynamic Interactions in Colloidal and Biological Sys-
tems. PhD thesis, Konstanz, 2006.
[115] M. Hütter and H. C. Öttinger.J. Chem. Soc., Faraday Trans.,
94:1403, 1998.
[116] J. C. Butcher.Numerical Methods for Ordinary Differential Equations.
John Wiley & Sons, New York, 2003.
[117] P. Grassia, E. J. Hinch, and L. C. Nitsche.J. Fluid Mech.,282:373,
1995.
[118] J. Wilkie.Phys. Rev. E,70:017701, 2004.
[119] G. H. Golub and C. F. Van Loan.Matrix computations. Johns Hopkins
University Press, Baltimore, 1996.
[120] L. N. Trefethen and D. Bau.Numerical linear algebra. Society for
Industrial and Applied Mathematics, Philadelphia, 1997.
113
REFERENCES
[121] G. E. P. Box and M. E. Muller.Ann. Stat. Mech.,29:610, 1958.
[122] J. Bell.Communications of the ACM,11:498, 1968.
[123] M. P. Stevens.Polymer chemistry: An introduction. Oxford University
Press, New York, 1999.
[124] R. M. Jendrejack, D. C. Schwartz, M. D. Graham, and J.J.
De Pablo.J. Chem. Phys.,119:1165, 2003.
[125] R. Khare, M. D. Graham, and J. J. De Pablo.Phys. Rev. Lett.,
96:224505, 2006.
[126] L. Cannavacciuolo, R. G. Winkler, and G. Gompper.Europhys.
Lett.,83:34007, 2008.
[127] O. Kratky and G. Porod.Rec. Trav. Chim. Pays-Bas.,68:1106, 1949.
[128] P. S. Doyle and P. T. Underhill.Handbook of Materials Modeling,
page 2619. Springer-Verlag, Berlin, 2005.
[129] L. D. Landau and E. M. Lifshitz.Statistical Physics. Pergamon Press,
London, 1958.
[130] P. A. Wiggins and P. C. Nelson.Phys. Rev. E,73:031906, 2006.
[131] D. L. Ermak and J. A. McCammon.J. Chem. Phys.,69:1352, 1978.
[132] D. Frenkel and R. Eppenga.Phys. Rev. E,31:1776, 1985.
[133] B. U. Felderhof.J. Chem. Phys.,137:084906, 2012.
[134] J. Bławzdziewicz, M. L. Ekiel-Jezewska, and E. Wajnryb.J.
Chem. Phys.,133:114703, 2010.
[135] B. U. Felderhof and A. Sellier.J. Chem. Phys.,136:054703, 2012.
[136] D. R. Mikulencak and J. F. Morris.J. Fluid Mech.,520:215, 2004.
[137] G. Subramanian and D. L. Koch.Phys. Fluids,18:074402, 2006.
[138] P. M. Kulkarni and J. F. Morris.J. Fluid Mech.,596:413, 2008.
[139] P. J. Mucha, S. -Y. Tee, D. A. Weitz, B. I. Shraiman, and M. P.
Brenner.J. Fluid Mech.,501:71, 2004.
114
REFERENCES
[140] S. Bhattacharya, J. Bławzdziewicz, and E. Wajnryb.Physica A,
356:294, 2005.
[141] S. Bhattacharya, J. Bławzdziewicz, and E. Wajnryb.J. Comput.
Phys.,212:718, 2006.
[142] J. P. Hernández-Ortiz, J. J. De Pablo, and M. D. Graham.J.
Chem. Phys.,125:164906, 2006.
[143] J. W. Swan and J. F. Brady.J. Fluid Mech.,687:254, 2011.
[144] M. E. Staben, A. Z. Zinchenkom, and H. D. Robert.Phys. Fluids,
15:1711, 2003.
[145] T. M: Squires and S. R. Quake.Rev. Mod. Phys.,77:977, 2005.
[146] B. U. Felderhof and R. B. Jones.J. Math. Phys.,30:339, 1989.
115
REFERENCES
116
Zusammenfassung
In den letzten Jahrzehnten wurden nicht nur die Computerchips immer kleiner,
auch die Größe von mikrofluidischen Geräten hat sich drastisch geändert. Diese
sogenannten Hosentaschenlabore haben sich enorm weiterentwickelt und finden
mittlerweile Anwendung in vielen physikalischen, biologischen oder chemischen
Systemen. Ein wesentlicher Bestandteil all dieser Geräte ist die präzise Mani-
pulation von suspendierten Teilchen auf der Mikrometerskala. Nicht nur für die
Entwicklung solcher Geräte, sondern auch für ein tieferes Verständnis biologi-
scher Vorgänge ist es von fundamentaler Bedeutung die nichtlineare Dynamik
von kolloidalen Systemen in begrenzten Geometrien zu verstehen.
In dieser Arbeit präsentieren wir zwei verschiedene kolloidale Modellsysteme,
bei denen die suspendierten Teilchen zwischen zwei parallelen Platten begrenzt
sind und durch eine druckgetriebene Poiseuille Strömung aus dem thermodyna-
mischen Gleichgewicht herausgetrieben werden.
Im ersten Modell untersuchen wir die Cross-Streamline-Migration von se-
miflexiblen Polymeren. Dazu beschreiben wir das Polymer mit dem Bead-Spring-
Modell, eine diskrete Realisierung des bekannten Worm-Like-Chain Modells. Die
hydrodynamischen Wechselwirkungen zwischen den Punktteilchen dieses Modells
beschreiben wir mit dem Zweiwand-Green-Tensor der Stokesschen Gleichungen.
Mit der Hilfe von Brownsche Dynamik Simulationen untersuchen wir die Schwer-
punktverteilung des Polymers quer des Kanals. Unsere Ergebnisse reproduzieren
die typische bimodale Verteilung, welche in früheren Experimenten und Simula-
tionen beobachtet wurde. Um die Ursache dieser Verteilung besser zu verstehen,
leiten wir die Smoluchowski Gleichung für die Schwerpunktverteilung her und
analysieren die verschiedenen Anteile der Wahrscheinlichkeitsstromdichte, welche
diese bimodale Verteilung verursachen. Im Gegensatz zu früheren Arbeiten, bei
denen die ortsabhängige Diffusionskonstante des Schwerpunktes als Verursacher
117
der Migrationsbewegung weg von der Kanalmitte galt, identifizieren wir eindeu-
tig den deterministischen Drift als Hauptverursacher. Wir zeigen sogar, dass der
diffusive Anteil der Wahrscheinlichkeitsstromdichte immer mehr an Bedeutung
verliert, je mehr wir die Biegesteifigkeit des Modellpolymers erhöhen. Als zweites
System untersuchen wir die nichtlineare Dynamik kräftefreier kugelförmiger Kol-
loide, deren Durchmesser in der Größenordnung der Kanalbreite liegt. Im Gegen-
satz zum vorherigen System können wir hier nicht mehr das vereinfachte Modell
von Punktteilchen annehemen. Demnach bestimmen wir die hydrodynamischen
Wechselwirkungen über eine Multipolentwicklung um die endliche Ausdehnung
der Teilchen zu berücksichtigen. Wir zeigen, dass zusätzlich zu den Trajektorien
eines Teilchenpaares, welche bereits in begrenzter und unbegrenzter Scheerströ-
mung beobachtet wurden, die begrenzte Poiseuille Strömung neue Klassen von
Trajektorien entstehen lässt. Je nach Anfangsbedingungen und Teilchengröße,
bewegen sich zwei Kugeln auf geschlossenen oder offenen Bahnen. Wir katego-
risieren diese Trajektorien in einem Zustandsdiagramm und diskutieren wie die
Ergebnisse des Zweiteilchensystems dazu beitragen die Stabilität in einer Reihe
von vielen Teilchen zu verstehen.
118
Abstract
In the recent years, not only the physical size of computer chips decreased, also
the size of microfluidic devices changed drastically. These so called Lab on a chip
devices have seen great progress and find applications in many physical, biological
or chemical systems. An essential part of all these devices is the precise manipu-
lation of suspended particles on the micron scale. Therefore, it is of fundamental
importance to understand the nonlinear dynamics of colloidal systems in confined
geometries, not only for the development of such devices, but also to gain a better
understanding of biological processes.
In this work we present two different colloidal model systems, where the sus-
pended particles are confined between two planar walls and are driven out of
equilibrium using the pressure driven Poiseuille flow.
In the first model, we investigate the cross-streamline migration of semiflexible
polymers. We introduce the semiflexible polymer as a bead-spring chain, which
is a discrete representation of the well known worm-like chain. We take the
hydrodynamic interactions between the pointlike beads into account by the two-
wall Green tensor of the Stokes equations. With the help of Brownian dynamics
simulations we investigate the probability distribution of the center-of-mass of
the polymer across the channel. Our simulations reproduce the typical bimodal
distribution, as observed in previous experiments and simulations. To gain a
better understanding of the origin of this distribution, we derive a Smoluchowski
equation for the center-of-mass distribution and carefully analyze the different
contributions to the probability current that causes the bimodal distribution. In
contrast to previous studies, where the migration away from the centerline was
explained by a position-dependent diffusivity, we clearly identify a deterministic
drift current as the major cause for migration away from the centerline. We even
show that diffusional currents due to a position-dependent diffusivity become less
119
important with increasing polymer stiffness.
In the second model, we investigate the dynamics of force- and torque-free
spherical colloids, whose diameter is comparable to the distance of the two walls.
In contrast to the previous model, we can no longer work with the simplified point
particles. Instead we have to take the finite size of the spheres into account and
determine the hydrodynamic interactions by a multipole expansion. In addition
to the trajectories of a pair of particles, which were already observed in bounded
and unbounded shear flow, we show that the bounded Poiseuille flow gives rise
to new classes of trajectories resulting in cross-streamline migration. Depending
on their initial positions and their size, two particles move along either closed or
open trajectories (thus exhibiting bound or unbound states). We categorize these
trajectories in a state diagram and discuss how the results on the two-particle
system help to understand the stability of particle trains composed of several
particles.
120
Danksagung
Es gibt eine Reihe von Menschen, bei denen ich mich bedanken möchte. Als erstes
möchte ich mich bei meinem Doktorvater Prof. Dr. Holger Stark bedanken, der
mir überhaupt erst meine Promotion ermöglicht hat. Er stand mir immer mit Rat
und Tat zur Seite und gab mir stets die Freiheit meinen eigenen Weg zu gehen.
Die konstruktiven Gespräche und Kritiken stellen einen wertvollen Beitrag dieser
Arbeit dar.
Auch möchte ich mich für die freundlichen und konstruktiven Gespräche bei
meinem Zweitbetreuer Prof. Dr. Alexander Mikhailov bedanken.
Mein Dank richtet sich auch an die ganze Arbeitsgruppe der AG Stark. Her-
vorheben möchte ich insbesondere Reinhard Vogel, Andreas Zöttl, Christopher
Prohm, Katrin Wolff und Oliver Pohl, die mich immer wieder mit Ratschlägen
unterstützt haben und mir beim Korrekturlesen geholfen haben.
Desweiteren bedanke ich mich bei der Deutsche Forschungsgemeinschaft, wel-
che mich im Rahmen des Graduiertenkollegs GRK 1558 finanziell unterstützt
hat.
Zu guter Letzt möchte ich mich bei meiner Familie bedanken, die mich immer
unterstützt haben und mir motivierend beiseite stand.
121
