1. No category

Brownian Suspensions of Rigid Particles
Aleksandar Donev &
Steven Delong, Bakytzhan Kallemov and Florencio Balboa
and others
Courant Institute, New York University
Hydrodynamics at Small Scales, SIAM CSE15
Salt Lake City, Utah,
March 18th 2015
A. Donev (CIMS)
R-BD
3/2015
1 / 22
Introduction
Non-Spherical Colloids near Boundaries
Figure: (Left) Cross-linked spheres; Kraft et al. [1]. (Right) Lithographed
boomerangs; Chakrabarty et al. [2].
A. Donev (CIMS)
R-BD
3/2015
2 / 22
Introduction
Bent Active Nanorods
Figure: From the Courant Applied Math Lab of Zhang and Shelley [3]
A. Donev (CIMS)
R-BD
3/2015
3 / 22
Introduction
Thermal Fluctuation Flips
QuickTime
A. Donev (CIMS)
R-BD
3/2015
4 / 22
Introduction
Steady Stokes Flow (Re → 0)
Consider
of Nb rigid bodies with positions
a suspension
Q = %1 , . . . , %Nb and orientations Θ = {θ 1 , . . . , θ Nb }.
We describe orientations using quaternions.
For viscous-dominated flows we can assume steady Stokes flow and
define the body mobility matrix N (Q, Θ),
[U , Ω]T = N [F , T ]T ,
where the left-hand side collects the linear U = {υ 1 , . . . , υ Nb } and
angular Ω = {ω 1 , . . . , ω Nb } velocities,
and the right hand side collects the applied forces
F (Q, Θ) = {F1 , . . . , FNb } and torques T (Q, Θ) = {τ 1 , . . . , τ Nb }.
A. Donev (CIMS)
R-BD
3/2015
5 / 22
Introduction
Brownian Motion
The Brownian motion of the rigid bodies is described by the
overdamped Langevin equation, symbolically:
1
dQ/dt
U
F
=
=N
+ (2kB T N ) 2 W (t) .
dΘ/dt
Ω
T
How to represent orientations using normalized quaternions and
handle the constraint kΘk k = 1?
What is the correct thermal drift (i.e., what does mean)?
1
How to compute (the action of) N and N 2 and simulate the
Brownian motion of the bodies?
A. Donev (CIMS)
R-BD
3/2015
6 / 22
Introduction
Difficulties/Goals
Stochastic drift It is crucial to handle stochastic calculus issues carefully
for overdamped Langevin dynamics. Since diffusion is slow
we also want to be able to take large time step sizes.
Complex shapes We want to stay away from analytical approximations
that only work for spherical particles.
Boundary conditions Whenever observed experimentally there are
microscope slips (glass plates) that modify the
hydrodynamics strongly. It is preferred to use no Green’s
functions but rather work in complex geometry.
Gravity Observe that in all of the examples above there is gravity and
the particles sediment toward the bottom wall, often very
close to the wall (∼ 100nm). This is a general feature of all
active suspensions but this is almost always neglected in
theoretical models.
Many-body Want to be able to scale the algorithms to suspensions of
many particles–nontrivial numerical linear algebra.
A. Donev (CIMS)
R-BD
3/2015
7 / 22
Blob models of complex particles
Blob/Bead Models
Figure: Blob or “raspberry” models of: a spherical colloid, and a lysozyme [4].
The rigid body is discretized through a number of “beads” or “blobs”
with positions Q = {q1 , . . . , qN }.
Describe the fluid-blob interaction using a localized smooth kernel
δa (r ) with compact support of size a giving the effective
hydrodynamic radius of the blob (diffuse sphere).
Standard in fluctuating/stochastic immersed boundary methods but
with stiff springs instead of truly rigid agglomerates.
A. Donev (CIMS)
R-BD
3/2015
8 / 22
Blob models of complex particles
Rigidly-Constrained Blobs
∇π − η∇2 v =
N
X
λi δa (qi − r) +
p
2ηkB T ∇ · W
i=1
∇ · v = 0 (Lagrange multiplier is π)
N
X
λi = F (Lagrange multiplier is υ)
(1)
i=1
N
X
qi − %0 × λi = τ (Lagrange multiplier is ω),
i=1
Z
∀i :
δa (qi − r) v (r, t) dr = υ + ω × qi − %0 + slip (Multiplier is λi )
A. Donev (CIMS)
R-BD
3/2015
9 / 22
Blob models of complex particles
Notation
Composite velocity U = {u1 , . . . , uN } and
rigidity forces Λ = {λ1 , . . . , λN }.
Define the composite local averaging linear operator J (Q) operator,
and the composite spreading linear operator, S (Q) = J ? (Q),
Z
ui = (J v)i = δa (qi − r) v (r, t) dr
λ (r) = (SΛ) (r) =
N
X
λi δa (qi − r) .
i=1
Denote the (potentially discrete) operators scalar gradient G ≡ ∇,
vector divergence D = −G? ≡ ∇·, tensor divergence Dv , and vector
Laplacian L = −Dv D?v ≡ ∇2 .
A. Donev (CIMS)
R-BD
3/2015
10 / 22
Blob models of complex particles
Saddle-Point Problem
Define the geometric matrix K that converts body kinematics to blob
kinematics,
U = KY = K [U , Ω]T = U + Ω × Q − Q0 .
We get the symmetric constrained Stokes saddle-point problem,

√


 
−ηL G −S 0
v
2ηkB T W
∇·
 −D 0

 

0
0 
0

 π  = 
,
 −J 0

0 K  Λ  
0
?
0
0 K
0
Y
R
where Y = [U , Ω]T and R = [F , T ]T , and recall that S = J ? .
A. Donev (CIMS)
R-BD
3/2015
11 / 22
Blob models of complex particles
Mobility Matrix
Eliminate velocity and pressure using the Schur complement
#
"
slip M −K
Λ
,
=
e
−K?
0
Y
− R+R
e = √2ηkB T K? M−1 J L−1 Dv W are the random
where R
(stochastic) forces and torques.
Here the all-important 3N × 3N blob mobility matrix M is
M = J L−1 S,
−1
where L−1 = −L−1 + L−1 G DL−1 G
DL−1 denotes the Stokes
solution operator.
A. Donev (CIMS)
R-BD
3/2015
12 / 22
Blob models of complex particles
Rigidly-Constrained Blobs
The physical interpretation is simple:
MΛ = KY + slip
e
K? Λ = R + R,
where the unknown Y = [U , Ω]T are the body kinematics,
e are the
R = [F , T ]T are the applied forces and torques and R
random (stochastic) forces and torques.
Here Λ are the unknown rigidity forces (Lagrange multipliers) acting
on the blobs that needs to be solved for.
The 3N × 3N block mobility matrix M has a simple pairwise
physical interpretation:
The 3 × 3 block Mij maps a force on blob j to a velocity of blob i,
Z
−1
Mij ≈ η
δa (qi − r)G(r, r0 )δa (qj − r0 ) drdr0
(2)
where G is the Green’s function (Oseen tensor for unbounded).
A. Donev (CIMS)
R-BD
3/2015
13 / 22
Blob models of complex particles
Suspensions of Rigid Bodies
Taking yet one more Schur complement we get
1
U
F
=N
+ (2kB T N ) 2 W.
Ω
T
The many-body mobility matrix N takes into account rigidity and
higher-order hydrodynamic interactions,
−1
N = K? M−1 K
.
If a fluctuating fluid solver is used it gives an explicit square root of
p
1
N 2 = 2kB T N K? M−1 J L−1 Dv W.
Observe that discrete fluctuation-dissipation balance is
guaranteed,
1 ?
1
N 2 N 2 = N K? M−1 J L−1 LL−1 S M−1 KN =
N K? M−1 MM−1 KN = N K? M−1 K N = N N −1 N = N .
A. Donev (CIMS)
R-BD
3/2015
14 / 22
Blob models of complex particles
How to Approximate the Mobility
If we have a way to approximate the (action of) the mobility matrix
M we can also do this without invoking a fluid solver.
We need to be able to solve
e
N −1 Y = K? M−1 K Y = R + R,
which we can either do using direct or iterative solvers.
There are different ways to obtain M:
In unbounded domains we can just use the Rotne-Prager-Yamakawa
tensor (RPY) (always SPD!).
In simple geometries such as a single wall we can use a generalization
of RPY [5].
For periodic domains we can use Ewald-type summations or
non-uniform FFTs with a fluctuating spectral fluid solver.
In more general cases we can use a fluctuating FEM/FVM fluid
Stokes solver [6, 7].
A. Donev (CIMS)
R-BD
3/2015
15 / 22
Results
Brownian motion under gravity
We consider the Brownian motion of a single rigid body near a no-slip
boundary.
Temporal integration of the overdamped equations is done using a
random finite different (RFD) approach as described by Steven
Delong.
Number of blobs is small and we have a simple geometry so we use
approximate Blake-Rotne-Prager tensor (Brady & Swan [8])
For this test we use direct linear algebra to compute N and
1
Cholesky factorization to compute N 2 .
We add gravity which makes the equilibrium Gibbs-Boltzmann
distribution be
mgh + Usteric
PGB (Q, Θ) ∼ exp −
,
kB T
where h is the center-of-mass height and Usteric is a Yukawa-type
repulsion with the wall.
A. Donev (CIMS)
R-BD
3/2015
16 / 22
Results
Quasi-2D Diffusion
Brownian motion is confined near the bottom wall so it quasi-two
dimensional.
Without external forcing the Brownian motion along the wall should
be isotropic diffusive at long time scales.
A naive guess for the effective 2D diffusion coefficient would be
the Gibbs-Boltzmann average of the parallel translational mobility:
Dk = kB T µk GB .
This is in fact a theorem for a sphere because rotational Brownian
motion does not change the mobility.
Is it true for non-spherical particles?
A. Donev (CIMS)
R-BD
3/2015
17 / 22
Results
MSD for a sphere
MSD(t) for Icosahedron with Hydrodynamic Radius = 0.5
12
10
RFD Icosahedron Parallel MSD
RFD Blob Parallel MSD (a = 0.5)
Blob Parallel Mobility
RFD Icosahedron Perpendicular MSD
RFD Blob perpendicular MSD
Icosahedron Asymptotic Perp MSD
MSD
8
6
4
2
00
50
100
time
150
200
Figure: Mean square displacement (MSD) for a non-uniform spherical particle of
unit diameter discretized as an icosahedron of 12 blobs or just a single blob.
A. Donev (CIMS)
R-BD
3/2015
18 / 22
Results
MSD for a tetrahedron
MSD(t) for Tetrahedron
MSD
15
Fixman Parallel MSD
RFD Parallel MSD
Parallel Mobility
Fixman Perpendicular MSD
RFD Perpendicular MSD
Asymptotic Perpendicular MSD
Average Perpendicular Mobility
10
5
00
100
200
300
time
400
500
600
700
Figure: MSD for a non-spherical particle (tetrahedron/tetramer).
A. Donev (CIMS)
R-BD
3/2015
19 / 22
Results
The choice of tracking point matters
MSD(t) for Free Tetrahedron
RFD Parallel MSD Vertex
RFD Parallel MSD CoM
Mu Parallel Vertex
Mu Parallel CoM
16
14
12
MSD
10
8
6
4
2
00
50
100
150
time
200
250
300
350
Figure: MSD for a non-spherical particle (tetrahedron/tetramer).
A. Donev (CIMS)
R-BD
3/2015
20 / 22
Results
Resolving lubrication forces
41−cyl
221−cyl
121−cyl
834−cyl
39−shell
Theory(dense)
6
10
F
µU
4
10
2
10
0
10 −2
10
−1
10
ε
0
10
Figure: The drag coefficient for a periodic array of cylinders in steady Stokes flow
for close-packed arrays with inter-particle gap ε, showing the correct asymptotic
5
ε− 2 lubrication force divergence.
A. Donev (CIMS)
R-BD
3/2015
21 / 22
Results
References
Daniela J. Kraft, Raphael Wittkowski, Borge ten Hagen, Kazem V. Edmond, David J. Pine, and Hartmut L¨
owen.
Brownian motion and the hydrodynamic friction tensor for colloidal particles of complex shape.
Phys. Rev. E, 88:050301, 2013.
Ayan Chakrabarty, Andrew Konya, Feng Wang, Jonathan V Selinger, Kai Sun, and Qi-Huo Wei.
Brownian motion of boomerang colloidal particles.
Physical review letters, 111(16):160603, 2013.
Daisuke Takagi, Adam B Braunschweig, Jun Zhang, and Michael J Shelley.
Dispersion of self-propelled rods undergoing fluctuation-driven flips.
Phys. Rev. Lett., 110(3):038301, 2013.
Jos´
e Garc´ıa de la Torre, Mar´ıa L Huertas, and Beatriz Carrasco.
Calculation of hydrodynamic properties of globular proteins from their atomic-level structure.
Biophysical Journal, 78(2):719–730, 2000.
Eligiusz Wajnryb, Krzysztof A Mizerski, Pawel J Zuk, and Piotr Szymczak.
Generalization of the rotne–prager–yamakawa mobility and shear disturbance tensors.
Journal of Fluid Mechanics, 731:R3, 2013.
S. Delong, F. Balboa Usabiaga, R. Delgado-Buscalioni, B. E. Griffith, and A. Donev.
Brownian Dynamics without Green’s Functions.
J. Chem. Phys., 140(13):134110, 2014.
Pat Plunkett, Jonathan Hu, Christopher Siefert, and Paul J Atzberger.
Spatially adaptive stochastic methods for fluid–structure interactions subject to thermal fluctuations in domains with
complex geometries.
Journal of Computational Physics, 277:121–137, 2014.
James W. Swan and John F. Brady.
Simulation of hydrodynamically interacting particles near a no-slip boundary.
Physics of Fluids, 19(11):113306, 2007.
A. Donev (CIMS)
R-BD
3/2015
22 / 22