Brownian Dynamics without Green’s Functions Aleksandar Donev Courant Institute, New York University & Steven Delong, Courant Rafael Delgado-Buscalioni, UAM Florencio Balboa Usabiaga, UAM Boyce Griffith, Courant
Particle Methods for Micro- and Nano-flows ECFD VI, Barcelona, Spain July 21st 2014
A. Donev (CIMS)
FIB
7/2014
1 / 25
Outline
1
Fluid-Particle Coupling
2
Overdamped Limit
3
Results
4
Outlook
A. Donev (CIMS)
FIB
7/2014
2 / 25
Fluid-Particle Coupling
Levels of Coarse-Graining
Figure: From Pep Espa˜ nol, “Statistical Mechanics of Coarse-Graining” A. Donev (CIMS)
FIB
7/2014
4 / 25
Fluid-Particle Coupling
Incompressible Fluctuating Hydrodynamics The colloidal are immersed in an incompressible fluid that we assume can be described by the time-dependent fluctuating incompressible Stokes equations, p (1) ρ∂t v + ∇π = η∇2 v + f + 2ηkB T ∇ · Z ∇ · v = 0, along with appropriate boundary conditions. Here the stochastic momentum flux is modeled via a random Gaussian tensor field Z(r, t) whose components are white in space and time with mean zero and covariance hZij (r, t)Zkl (r0 , t 0 )i = (δik δjl + δil δjk ) δ(t − t 0 )δ(r − r0 ).
A. Donev (CIMS)
FIB
7/2014
(2)
5 / 25
Fluid-Particle Coupling
Brownian Particle Model Consider a Brownian “particle” of size a with position q(t) and ˙ and the velocity field for the fluid is v(r, t). velocity u = q, We do not care about the fine details of the flow around a particle, which is nothing like a hard sphere with stick boundaries in reality anyway. Take an Immersed Boundary Method (IBM) approach and describe the fluid-blob interaction using a localized smooth kernel δa (∆r) with compact support of size a (integrates to unity). Often presented as an interpolation function for point Lagrangian particles but here a is a physical size of the particle (as in the Force Coupling Method (FCM) of Maxey et al). We will call our particles “blobs” since they are not really point particles. A. Donev (CIMS)
FIB
7/2014
6 / 25
Fluid-Particle Coupling
Local Averaging and Spreading Operators Postulate a no-slip condition between the particle and local fluid velocities, Z q˙ = u = [J (q)] v = δa (q − r) v (r, t) dr, where the local averaging linear operator J(q) averages the fluid velocity inside the particle to estimate a local fluid velocity. The induced force density in the fluid because of the force F applied on particle is: f = −Fδa (q − r) = − [S (q)] F, where the local spreading linear operator S(q) is the reverse (adjoint) of J(q). The physical volume of the particle ∆V is related to the shape and width of the kernel function via Z −1 ∆V = (JS)−1 = δa2 (r) dr . (3) A. Donev (CIMS)
FIB
7/2014
7 / 25
Fluid-Particle Coupling
Fluctuation-Dissipation Balance One must ensure fluctuation-dissipation balance in the coupled fluid-particle system. The stationary (equilibrium) distribution must be the Gibbs distribution Peq (q) = Z −1 exp (−U(q)/kB T ) , (4) where F(q) = −∂U(q)/∂q with U(q) a conservative potential. No entropic contribution to the coarse-grained free energy because our formulation is isothermal and the particles do not have internal structure. In order to ensure that the dynamics is time reversible with respect to an appropriate Gibbs-Boltzmann distribution), the thermal or stochastic drift forcing f th = (kB T ) ∂q · S (q)
(5)
needs to be added to the fluid equation [1, 2, 3]. A. Donev (CIMS)
FIB
7/2014
8 / 25
Overdamped Limit
Viscous-Dominated Flows We consider n spherical neutrally-buoyant particles of radius a in d dimensions, having spatial positions q = {q1 , . . . , qN } with (1) (d) qi = (qi , . . . , qn ). Let script J and S denote composite fluid-particles interaction operators. Let us assume that the Schmidt number is very large, Sc = η/ (ρχ) 1, where χ ≈ kB T / (6πηa) is a typical value of the diffusion coefficient of the particles [4]. To obtain the asymptotic dynamics in the limit Sc → ∞ heuristically, we delete the inertial term ρ∂t v in (1), ∇ · v = 0 and p ∇π = η∇2 v + SF + 2ηkB T ∇ · Z ⇒ (6) p v = η −1 L−1 SF + 2ηkB T ∇ · Z , where L−1 0 is the Stokes solution operator. A. Donev (CIMS)
FIB
7/2014
10 / 25
Overdamped Limit
Overdamped Limit A rigorous adiabatic mode elimination procedure informs us that the correct interpretation of the noise term in this equation is the kinetic stochastic integral, s " # dq (t) 2kB T −1 1 = J (q)L S(q)F(q) + ∇ Z (r, t) . (7) dt η η This is equivalent to the standard equations of Brownian Dynamics (BD), 1 dq = MF + (2kB T M) 2 W(t)+kB T (∂q · M), dt
(8)
where M(q) 0 is the symmetric positive semidefinite (SPD) mobility matrix M = η −1 J L−1 S. A. Donev (CIMS)
FIB
7/2014
11 / 25
Overdamped Limit
Brownian Dynamics via Fluctuating Hydrodynamics It is not hard to show that M is very similar to the Rotne-Prager mobility used in BD, for particles i and j, Z −1 Mij = η δa (qi − r)K(r, r0 )δa (qj − r0 ) drdr0 (9) where K is the Green’s function for the Stokes problem (Oseen tensor for infinite domain). The self-mobility defines a consistent hydrodynamic radius of a blob, Mii = Mself =
1 I. 6πηa
For well-separated particles we get the correct Faxen expression, r=q a2 2 a2 2 −1 Mij ≈ η I + ∇r I + ∇r0 K(r − r0 ) r0 =qj . i 6 6 At smaller distances the mobility is regularized in a natural way and positive-semidefiniteness ensured automatically. A. Donev (CIMS)
FIB
7/2014
12 / 25
Overdamped Limit
Numerical Methods Both compressible and incompressible, inertial and overdamped, numerical methods have been implemented by Florencio Balboa (UAM) on GPUs for periodic BCs (public-domain!), and in the parallel IBAMR code of Boyce Griffith by Steven Delong for general boundary conditions (to be made public-domain next fall!). Spatial discretization is based on previously-developed staggered schemes for fluctuating hydro [5] and the immersed-boundary method kernel functions of Charles Peskin. Temporal discretization follows a second-order splitting algorithm (move particle + update momenta), and is limited in stability only by advective CFL. We have constructed specialized temporal integrators that ensure discrete fluctuation-dissipation balance, including for the overdamped case. A. Donev (CIMS)
FIB
7/2014
13 / 25
Overdamped Limit
(Simple) Midpoint Scheme Fluctuating Immersed Boundary Method (FIBM) method: Solve a steady-state Stokes problem (here δ 1) p ∇π n = η∇2 vn + 2ηkB T ∇ · Z n + Sn F (qn ) kB T δ fn δ fn n n fn + S q + W −S q − W W δ 2 2 ∇ · vn = 0. Predict particle position: 1
qn+ 2 = qn +
∆t n J v 2
Correct particle position, 1
qn+1 . = qn + ∆tJ n+ 2 v. A. Donev (CIMS)
FIB
7/2014
14 / 25
Results
Slit Channel
PDF × a
0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.000
10
Euler-Maruyama (biased) Midpoint Biased Distribution Gibbs Boltzmann
20
30 40 H/a
50
60
Figure: Probability distribution of the distance H to one of the walls for a freely-diffusing blob in a two dimensional slit channel. A. Donev (CIMS)
FIB
7/2014
16 / 25
Results
Equilibrium Radial Correlation Function
2.50
dt=10 dt=50 MC (dense) dt=10 dt=100 MC (dilute)
RDF
2.00 1.50 1.00 0.50 0.00 0.5
1
2 r/σ
1.5
2.5
3
3.5
Figure: Equilibrium radial distribution function g2 (r) for a suspension of blobs interacting with a repulsive LJ (WCA) potential. A. Donev (CIMS)
FIB
7/2014
17 / 25
Results
Hydrodynamic Interactions 20 me=mf, α=0.01 me=mf, α=0.1 me=mf, α=0.25 Rotne-Prager (RP) RP + Lubrication RPY
F / FStokes
16 12 8 4 0.5
1
2 Distance d/RH
4
8
Figure: Effective hydrodynamic force between two approaching blobs at small 2F0 Reynolds numbers, FFSt = − 6πηR . H vr A. Donev (CIMS)
FIB
7/2014
18 / 25
Results
Diffusive Dynamics At long times, the motion of the particle is diffusive with a diffusion R∞ coefficient χ = limt→∞ χ(t) = t=0 C (t)dt, where E ∆q 2 (t) 1 D = [q(t) − q(0)]2 . 2t 2dt The Stokes-Einstein relation predicts χ(t) =
χ=
kB T kB T (Einstein) and χSE = (Stokes), µ 6πηa
(10)
where for our blob a is on the order of the fluid solver grid spacing. The dimensionless Schmidt number Sc = ν/χSE controls the separation of time scales. Self-consistent theory predicts a correction to Stokes-Einstein’s relation for small Sc , χ kB T χ ν+ = . 2 6πρa A. Donev (CIMS)
FIB
7/2014
19 / 25
Results
Stokes-Einstein Corrections
χ / χSE
1
0.9 Self-consistent theory From VACF integral From mobility
0.8 1
2
4
8 16 32 64 128 256 Schmidt number Sc
Figure: Corrections to Stokes-Einstein with changing viscosity ν = η/ρ, me = mf = ρ∆V . A. Donev (CIMS)
FIB
7/2014
20 / 25
Results
Colloidal Gellation: Cluster collapse 7.4 BD with HI BD without HI FIBM (4pt, IBAMR) FIBM (3pt, fluam)
7.2
Rg
7.0 6.8 6.6 6.4 6.2 6.0 -2 10
10
-1
0
10 10 3 t / tB = t kBT / ηa
1
10
2
Figure: Relaxation of the radius of gyration of a colloidal cluster of 13 spheres toward equilibrium, taken from Furukawa+Tanaka. A. Donev (CIMS)
FIB
7/2014
21 / 25
Outlook
Immersed Rigid Blobs
Unlike a rigid sphere, a blob particle would not perturb a pure shear flow. In the far field our blob particle looks like a force monopole (stokeset), and does not exert a force dipole (stresslet) on the fluid. Similarly, since here we do not include angular velocity degrees of freedom, our blob particle does not exert a torque on the fluid (rotlet). It is possible to include rotlet and stresslet terms, as done in the fluctuating force coupling method [6] and Stokesian Dynamics in the deterministic setting. Maintaining fluctuation-dissipation balance more challenging.
A. Donev (CIMS)
FIB
7/2014
23 / 25
Outlook
Conclusions Fluctuating hydrodynamics seems to be a very good coarse-grained model for fluids, and coupled to immersed particles to model Brownian suspensions. The minimally-resolved blob approach provides a low-cost but reasonably-accurate representation of rigid particles in flow. We have recently successfully extended the blob approach to reaction-diffusion problems (with Amneet Bhalla and Neelesh Patankar). Particle and fluid inertia can be included in the description, or, an overdamped limit can be taken if Sc 1. More complex particle shapes can be built out of a collection of blobs to form a rigid body.
A. Donev (CIMS)
FIB
7/2014
24 / 25
Outlook
References P. J. Atzberger. Stochastic Eulerian-Lagrangian Methods for Fluid-Structure Interactions with Thermal Fluctuations. J. Comp. Phys., 230:2821–2837, 2011. F. Balboa Usabiaga, R. Delgado-Buscalioni, B. E. Griffith, and A. Donev. Inertial Coupling Method for particles in an incompressible fluctuating fluid. Comput. Methods Appl. Mech. Engrg., 269:139–172, 2014. Code available at https://code.google.com/p/fluam. 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. F. Balboa Usabiaga, X. Xie, R. Delgado-Buscalioni, and A. Donev. The Stokes-Einstein Relation at Moderate Schmidt Number. J. Chem. Phys., 139(21):214113, 2013. F. Balboa Usabiaga, J. B. Bell, R. Delgado-Buscalioni, A. Donev, T. G. Fai, B. E. Griffith, and C. S. Peskin. Staggered Schemes for Fluctuating Hydrodynamics. SIAM J. Multiscale Modeling and Simulation, 10(4):1369–1408, 2012. Eric E. Keaveny. Fluctuating force-coupling method for simulations of colloidal suspensions. J. Comp. Phys., 269(0):61 – 79, 2014.
A. Donev (CIMS)
FIB
7/2014
25 / 25