Brownian Dynamics of Fully Confined Suspensions of Rigid Particles Without Green's Functions

TL;DR

This paper introduces a scalable RB-FIB method for simulating Brownian suspensions of rigid particles in confined domains without Green's functions, accurately capturing equilibrium and dynamic behaviors.

Contribution

The authors develop a novel RB-FIB approach that avoids explicit Green's functions, enabling efficient large-scale Brownian dynamics simulations of confined rigid particle suspensions.

Findings

Accurately reproduces Gibbs-Boltzmann equilibrium distribution.

Demonstrates particle dynamics in confined geometries and depinning transitions.

Shows boundary particles are displaced first in confined flows.

Abstract

We introduce a Rigid-Body Fluctuating Immersed Boundary (RB-FIB) method to perform large-scale Brownian dynamics simulations of suspensions of rigid particles in fully confined domains, without any need to explicitly construct Green's functions or mobility operators. In the RB-FIB approach, discretized fluctuating Stokes equations are solved with prescribed boundary conditions in conjunction with a rigid-body immersed boundary method to discretize arbitrarily-shaped colloidal particles with no-slip or active-slip prescribed on their surface. We design a specialized Split--Euler--Maruyama temporal integrator that uses a combination of random finite differences to capture the stochastic drift appearing in the overdamped Langevin equation. The RB-FIB method presented in this work only solves mobility problems in each time step using a preconditioned iterative solver, and has a computational complexity that scales linearly in the number of particles and fluid grid cells. We demonstrate that the RB-FIB method correctly reproduces the Gibbs-Boltzmann equilibrium distribution, and use the method to examine the time correlation functions for two spheres tightly confined in a cuboid. We model a quasi--two-dimensional colloidal crystal confined in a narrow microchannel and hydrodynamically driven across a commensurate periodic substrate potential mimicking the effect of a corrugated wall. We observe partial and full depinning of the colloidal monolayer from the substrate potential above a certain wall speed, consistent with a transition from static to kinetic friction through propagating kink solitons. Unexpectedly, we find that particles nearest the boundaries of the domain are the first to be displaced, followed by particles in the middle of the domain.

Figures

Equations

$$\rho {\partial }_t\mathbf{v}$$

$$\nabla \cdot \mathbf{v}$$

$$⟨{\mathcal{Z}}_{ij}(\mathbf{x},t){\mathcal{Z}}_{kl}({\mathbf{x}}^{\prime },t^{\prime })⟩=\left({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk}\right)\delta (t-t^{\prime })\delta (\mathbf{x}-{\mathbf{x}}^{\prime }).$$

$$\mathbf{v}\left(\mathbf{r}\right)={\mathbf{u}}_p+\left(\mathbf{r}-{\mathbf{q}}_p\right)\times {\mathbf{\omega }}_p+{\breve{\mathbf{u}}}_p.$$

$$m_p\frac{\partial {\mathbf{u}}_p}{\partial t}$$

$${\mathbf{I}}_p\frac{\partial {\mathbf{\omega }}_p}{\partial t}$$

$$\displaystyle-\eta\boldsymbol{\nabla}^{2}\boldsymbol{v}+\boldsymbol{\nabla}\pi=\boldsymbol{g},\$$

$$\mathbf{v}(\mathbf{r})$$

$${\mathbf{f}}_p={\int }_{\partial {\mathcal{B}}_p}\mathbf{\sigma }\cdot \mathbf{n}\left(\mathbf{r}\right)dA(\mathbf{r}),{\mathbf{\tau }}_p$$

$$\mathbf{U}$$

$$=\mathcal{N}\mathbf{F}+k_{B}T\left({\partial }_{\mathbf{Q}}\cdot \mathcal{N}\right)+\sqrt{2k_{B}T}{\mathcal{N}}^{1/2}\mathcal{W}$$

$${\mathbf{U}}_{\text{EM}}=\mathcal{N}\mathbf{F}+\sqrt{\frac{2k_{B}T}{\Delta t}}{\mathcal{N}}^{1/2}\mathbf{W},$$

$$-\eta {\nabla }^2\mathbf{v}+\nabla \pi =\sqrt{\frac{2k_{B}T\eta }{\Delta t}}$$

$$\nabla \cdot \mathbf{v}$$

$${\int }_{\Omega }\delta \left(\mathbf{x}-\mathbf{r}\right)\mathbf{v}\left(\mathbf{x}\right)dV(\mathbf{x})$$

$${\mathbf{f}}_p={\int }_{\partial {\mathcal{B}}_p}\mathbf{\lambda }\left(\mathbf{r}\right)dA(\mathbf{r}),{\mathbf{\tau }}_p$$

$$-\eta \mathbb{L}\mathbf{v}+\mathbf{G}\pi$$

$$\mathbf{D}\mathbf{v}$$

$$-{\mathcal{L}}^{-1}=\frac{1}{\eta }\left({\mathbb{L}}^{-1}-{\mathbb{L}}^{-1}\mathbf{G}{\left(\mathbf{D}{\mathbb{L}}^{-1}\mathbf{G}\right)}^{-1}\mathbf{D}{\mathbb{L}}^{-1}\right),$$

$${\left(\mathcal{J}\mathbf{v}\right)}_i^p$$

$${\left(\mathcal{S}\mathbf{\lambda }\right)}_{\alpha }$$

$${\left(\mathcal{K}\mathbf{U}\right)}_i^p$$

$${\left({\mathcal{K}}^T\mathbf{\lambda }\right)}_p$$

$$-\eta \mathbb{L}\mathbf{v}+\mathbf{G}\pi$$

$$\mathbf{D}\mathbf{v}$$

$$\mathcal{J}\mathbf{v}$$

$${\mathcal{K}}^T\mathbf{\lambda }$$

$${\mathcal{M}}_{ij}\approx {\iint }_{\Omega \times \Omega }{\delta }_h\left(\mathbf{x}-{\mathbf{r}}_i\right)\mathbb{G}\left(\mathbf{x},\mathbf{y}\right){\delta }_h\left(\mathbf{y}-{\mathbf{r}}_j\right)dV(\mathbf{y})dV(\mathbf{x})$$

$$\mathcal{M}\mathbf{\lambda }$$

$${\mathcal{K}}^T\mathbf{\lambda }$$

$${\mathbf{U}}_{\text{EM}}=\mathcal{N}\mathbf{F}-\mathcal{N}{\mathcal{K}}^T{\mathcal{M}}^{-1}\breve{\mathbf{u}}+\sqrt{\frac{2k_{B}T}{\Delta t}}{\mathcal{N}}^{1/2}\mathbf{W},$$

$${\partial }_{\mathbf{Q}}\cdot \mathbf{R}\left(\mathbf{Q}\right)$$

$$=\underset{\delta \to 0}{\lim }\frac{1}{\delta }⟨\mathbf{R}\left(\mathbf{Q}+\frac{\delta }{2}\Delta \mathbf{Q}\right)-\mathbf{R}\left(\mathbf{Q}-\frac{\delta }{2}\Delta \mathbf{Q}\right)⟩\Delta \mathbf{P},$$

$${\partial }_{\mathbf{Q}}\cdot \mathcal{N}={\partial }_{\mathbf{Q}}\cdot {\left({\mathcal{K}}^T{\mathcal{M}}^{-1}\mathcal{K}\right)}^{-1}=$$

$$-\mathcal{N}\left({\partial }_{\mathbf{Q}}{\mathcal{K}}^T\right):{\mathcal{M}}^{-1}\mathcal{K}\mathcal{N}-\mathcal{N}{\mathcal{K}}^T{\mathcal{M}}^{-1}\left({\partial }_{\mathbf{Q}}\mathcal{M}\right):{\mathcal{M}}^{-1}\mathcal{K}\mathcal{N}+\mathcal{N}{\mathcal{K}}^T{\mathcal{M}}^{-1}\left({\partial }_{\mathbf{Q}}\mathcal{K}\right):\mathcal{N}.$$

$${\partial }_{\mathbf{Q}}\cdot \mathcal{M}={\partial }_{\mathbf{Q}}\cdot \left(\mathcal{J}{\mathcal{L}}^{-1}\mathcal{S}\right)=\left({\partial }_{\mathbf{Q}}\mathcal{J}\right):{\mathcal{L}}^{-1}\mathcal{S}+\mathcal{J}{\mathcal{L}}^{-1}\left({\partial }_{\mathbf{Q}}\cdot \mathcal{S}\right).$$

$${\mathbf{W}}^{\text{FT}}=k_{B}T\left[\begin{array}{c}\frac{1}{L_p}{\mathbf{W}}_p^{\text{f}}\\ {\mathbf{W}}_p^{\tau }\end{array}\right],$$

$${\mathbf{U}}^{RFD}$$

$${\mathbf{\lambda }}^{RFD}$$

Full text

I Introduction

Suspensions of micron sized rigid colloidal particles can be investigated by an ever expanding set of experimental tools such as optical tweezers and externally applied gravitational, magnetic and electric fields. Passive and active colloidal suspensions studied in the lab or in nature are almost always confined to be near one or more physical boundaries such as microscope slips ConfinedSphere_Sedimented ; BoomerangDiffusion . The confinement geometry strongly influences the underlying hydrodynamics of micro-particle suspensions due to the long-ranged nature of hydrodynamic interactions in the steady Stokes regime AutophoreticSpheres_Adhikari , and boundaries may also play a key role in propulsion mechanisms in active suspensions ActiveDimers_EHD . In recent years, precise and detailed measurements have been performed on quasi–two-dimensional (quasi–2D) colloidal crystal monolayers in order to interrogate kinetic friction at the microscopic scale Bohlein2011 ; ColloidSheetFriction , but numerical simulations of this kind of system with proper accounting for hydrodynamic interactions are lacking. In this paper we introduce a numerical method that can simulate Brownian suspensions of rigid colloidal particles of complex shape in fully confined rectangular domains, with arbitrary combinations of periodic, no-slip, or free-slip boundary conditions along different dimensions.

Designing scalable simulation techniques for Brownian suspensions of many passive or active non-spherical colloids in confined domains is still an outstanding challenge in the field. While a number of existing methods can efficiently handle triply-periodic domains using Ewald techniques BrownianDynamics_OrderNlogN ; StokesianDynamics_Brownian ; FluctuatingFCM_DC ; SpectralSD , many experiments are carried out in some form of tight confinement, such as the dynamics of particles pressed between two microscope slides ColloidsInTightSlit ; BoomerangDiffusion or flowing through a microchannel. Stokesian Dynamics (SD) has become an industry standard in chemical engineering circles for Brownian suspensions BrownianDynamics_OrderNlogN ; StokesianDynamics_Brownian ; SD_SpectralEwald ; SpectralSD , and the method has been adapted to non-spherical rigid particles RigidBody_SD ; StokesianDynamics_Rigid ; StokesianDynamics_Confined , as well as particles in partial or full confinement StokesianDynamics_Confined ; StokesianDynamics_Wall ; StokesianDynamics_Slit . However, traditional SD and recent closely-related boundary-integral AutophoreticSpheres_Adhikari ; BoundaryIntegralWall_Adhikari ‘implicit-fluid’ methods remain limited to spherical particles in specific geometries, for which a (grand) mobility tensor can be constructed explicitly. Furthermore, existing methods only have computational complexity that scales linearly in the number of particles for periodic boundary conditions, for which Fast Fourier Transforms (FFTs) accelerate the many-body computations BrownianDynamics_OrderNlogN ; StokesianDynamics_Brownian ; FluctuatingFCM_DC ; SpectralSD . Fluctuating Lattice Boltzmann simulations have also been used for suspensions LB_SoftMatter_Review . These techniques, however, require introducing artificial fluid compressibility and fluid inertia, which imposes a severe restriction on the time step size in order to achieve a physically-realistic Schmidt number.

The need to explicitly construct and compute Green’s functions for nontrivial boundary conditions can be avoided by using a grid-based solver for the fluctuating Stokes equations in order to compute the action of the Green’s function and generate Brownian forces. The authors of SELM_FEM presented a Finite Element Method (FEM) capable of simulating Brownian suspensions in very general, confined domains. This method however, was limited to minimally resolved (point-like) spherical particles and did not correctly account for the stochastic drift term that appears in the overdamped Langevin equations when the mobility is configuration dependent. The authors of DECORATO_FEM proposed a FEM scheme with body-fitted grids that is capable of simulating Brownian particles with arbitrary shape in confined domains. However, the method requires complex remeshing every time step and the temporal integrator used in this work requires solving expensive resistance problems. Because of this, the body-fitted FEM approach does not scale well in the number of particles and is in practice limited to one or a few individual particles.

The Rigid-Body Fluctuating Immersed Boundary Method (RB-FIB) method presented in this work is, to our knowledge, the first method that can simulate Brownian suspensions of rigid particles with arbitrary shape in fully confined domains with controllable accuracy and in computational time that scales linearly in the number of particles (at finite packing fractions). The method is built on contributions from a number of past works by some of us. In BrownianBlobs , the authors presented a Fluctuating Immersed Boundary (FIB) method which could simulate fluctuating suspensions of minimally resolved spheres (or blobs) in general physical domains. In RigidMultiblobs , the authors presented a rigid multiblob method to simulate deterministic suspensions in general domains by constructing complex particle shapes out of agglomerates of minimally resolved spheres/blobs, termed rigid multiblobs. Both BrownianBlobs and RigidMultiblobs employ the Immersed Boundary (IB) method to handle the fluid-particle coupling. IB methods provide a low-accuracy but inexpensive and flexible alternative to body-fitted grid-based methods since no remeshing is required as particles move around in the domain. In BrownianMultiblobSuspensions , some of us presented an efficient temporal integration scheme to simulate the dynamics of Brownian suspensions with many rigid particles of arbitrary shape confined above a single no-slip wall. This method relies on the simple geometry of the physical domain for which an explicit form for the hydrodynamic mobility operator is available. Here we present an amalgamation of these past approaches: we develop a generalization of the temporal integration of BrownianMultiblobSuspensions that fits the IB framework used in RigidMultiblobs to handle the hydrodynamic interactions including with boundaries, and use the fluctuating hydrodynamics approach proposed in BrownianBlobs to account for Brownian motion. We develop a novel Split Euler–Maruyama (SEM) temporal integration scheme to capture the stochastic drift which strongly affects even a single particle. The SEM scheme modifies the preconditioned Krylov method of RigidMultiblobs to maintain its linear scaling but includes the necessary stochastic contributions to the rigid-body dynamics.

In Bohlein2011 the authors experimentally observed soliton wave patterns in a driven colloidal monolayer moving above a bottom substrate. The bottom wall is patterned with a periodic potential meant to mimic surface roughness. The monolayer is forced into quasi-2D confinement by laser-induced forces, and driven by the flow generated by moving the bottom wall. Brownian motion is crucial in activating the transitions of the colloids between the minima of the patterned potential, and must be captured accurately to resolve the dynamics of the monolayer. In section V.2, we use the RB-FIB method to numerically investigate a modified version of the experiment performed in Bohlein2011 where we confine the monolayer in a thin microchannel. We observe novel wave patters in the colloidal monolayer which emerge due to the physical confinement. While simulations of this type have been performed ColloidSheetFriction ; LatticeNoHydro , our work is, to our knowledge, the first which includes an accurate treatment of the hydrodynamics.

The rest of this work is organized as follows. In section II we give the continuous formulation of the problem and introduce some relevant notation. In section III we formulate the spatial discretization of the continuum equations. In section IV we introduce the SEM scheme as an efficient temporal discretization that maintains discrete fluctuation dissipation balance. To numerically validate our scheme, we consider several test cases. Appendix A considers a boomerang shaped particle in a slit channel, and confirms that the RB-FIB method is first order weakly accurate for expectations with respect to the Gibbs-Boltzmann equilibrium distribution. Section V.1 considers two spherical particles trapped in a tight cuboidal box to examine the effect of spatial resolution on dynamic statistics. In section V.2 we study the transition from static to dynamic friction for a suspension of many spherical colloids confined to a quasi–two dimensional slit channel, and hydrodynamically driven across a periodic substrate potential by translating the microchannel with constant velocity. We conclude with a summary and discussion of future directions in section VI

II Continuum Formulation

The fluctuating Stokes equations in a physical domain $\Omega$ can be written as FluctHydroNonEq_Book

[TABLE]

where the fluid has density $\rho$, shear viscosity $\eta$, and temperature $T$, all of which we take to be constant in this work. Denoting Cartesian coordinates with $\boldsymbol{x}\in\Omega$ and time with $t$, $\boldsymbol{g}\left(\boldsymbol{x},t\right)$ represents a fluid body force and $\sigma=-\pi\boldsymbol{I}+\eta\left(\boldsymbol{\nabla}\boldsymbol{v}+\boldsymbol{\nabla}\boldsymbol{v}^{T}\right)$ is the dissipative component of the fluid stress tensor, where $\boldsymbol{v}\left(\boldsymbol{x},t\right)$ is the fluid velocity, and $\pi\left(\boldsymbol{x},t\right)$ is the pressure. The stochastic stress tensor $\sqrt{2k_{B}T\eta}\boldsymbol{\mathcal{Z}}$ accounts for the fluctuating contribution to the fluid stress and ensures fluctuation dissipation balance, where $\boldsymbol{\mathcal{Z}}\left(\boldsymbol{x},t\right)$ is a symmetric random Gaussian tensor whose components are delta correlated in space and time,

[TABLE]

We consider $N_{b}$ arbitrarily shaped rigid particles (bodies) $\mathcal{B}_{p}$, $1\leq p\leq N_{b}$, suspended in the fluctuating Stokesian fluid. Because the particles are treated as rigid, their position in space can be completely described by the Cartesian location $\boldsymbol{q}_{p}(t)$ of a representative tracking point, and an orientation $\boldsymbol{\theta}_{p}(t)$ relative to a reference configuration. The bulk of our discussion will be independent of how one wishes to describe the orientation of the particles, however, in practice, we use unit quaternions, as described in BrownianMultiBlobs . The composite configuration of particle $p$ will be denoted by $\boldsymbol{Q}_{p}=\left[\boldsymbol{q}_{p},\boldsymbol{\theta}_{p}\right]$.

Rigid body $p$ moves with translational velocity $\boldsymbol{u}_{p}=d\boldsymbol{q}_{p}/dt$ and rotates with angular velocity $\boldsymbol{\omega}_{p}$ around the tracking point. The composite velocity $\boldsymbol{U}_{p}=\left[\boldsymbol{u}_{p},\boldsymbol{\omega}_{p}\right]$ can be used to express the velocity of an arbitrary point $\boldsymbol{r}\in\partial\mathcal{B}_{p}$ on the surface of the particle to give the no-slip condition

[TABLE]

Here $\breve{\boldsymbol{u}}_{p}$ is an apparent slip between the fluid and the particle surface which may be freely prescribed up to the condition that the integral of the slip velocity vanish over the surface of the particle FBIM . The free slip velocity can be used to account for active layers formed by, say, electrokinetic flows or beating flagella.

Newton’s second law gives the acceleration of particle $p$ in terms of the applied force $\boldsymbol{f}_{p}$ and the applied torque $\boldsymbol{\tau}_{p}$,

[TABLE]

where $\boldsymbol{n}\left(\boldsymbol{r}\right)$ is the outward pointing surface normal vector, and $m_{p}$ and $\boldsymbol{I}_{p}$ are the particles mass and angular inertia tensor respectively.

In this work we are interested in the overdamped or steady Stokes limit of the dynamics, which is the one relevant for colloidal suspensions due to the very large Schmidt numbers and very small Reynolds numbers. In the absence of thermal fluctuations, the steady Stokes or inertia-less limit of equations (3)-(4), (2), and (1) is taken by simply deleting the inertial terms to obtain

[TABLE]

The solution to the linear system of equations (5)–(7) can be expressed using a symmetric, positive semi-definite body mobility matrix $\boldsymbol{\mathcal{N}}(\boldsymbol{Q})$ which is a function of the composite configuration $\boldsymbol{Q}=\left[\boldsymbol{q}_{1},\boldsymbol{\theta}_{1},\ldots,\boldsymbol{q}_{N_{b}},\boldsymbol{\theta}_{N_{b}}\right]$. The body mobility matrix acts on the external forces and torques applied to each particle $\boldsymbol{F}=\left[\boldsymbol{f}_{1},\boldsymbol{\tau}_{1},\ldots,\boldsymbol{f}_{N_{b}},\boldsymbol{\tau}_{N_{b}}\right]$ to produce a composite vector $\boldsymbol{U}=\left[\boldsymbol{u}_{1},\boldsymbol{\omega}_{1},\ldots,\boldsymbol{u}_{N_{b}},\boldsymbol{\omega}_{N_{b}}\right]$ of the rigid body velocities of the particles, $\boldsymbol{U}=\boldsymbol{\mathcal{N}}\boldsymbol{F}$.

If we include the thermal fluctuations, the overdamped limit of equations (3)-(4), (2), and (1) is the overdamped Langevin equation LangevinDynamics_Theory ; FBIM

[TABLE]

where $\diamond$ is the kinetic stochastic product KineticStochasticIntegral_Ottinger , and $\boldsymbol{\mathcal{W}}(t)$ is a collection of independent Wiener processes. Equation (9) is the conversion of (8) to Ito form, by introducing the term $k_{B}T\left(\partial_{\boldsymbol{Q}}\cdot\boldsymbol{\mathcal{N}}\right)$ which we call the ‘thermal’ or ‘stochastic’ drift 111The precise mathematical interpretation of this notation when $\boldsymbol{Q}$ includes particle orientations expressed in terms of quaternions is explained in BrownianMultiBlobs .. Fluctuation dissipation balance is maintained in (9) through the random vector $\boldsymbol{\mathcal{N}}^{1/2}\boldsymbol{\mathcal{W}}$ such that $\boldsymbol{\mathcal{N}}^{1/2}\left(\boldsymbol{\mathcal{N}}^{1/2}\right)^{\star}=\boldsymbol{\mathcal{N}}$, where star denotes an $L_{2}$ adjoint (i.e., a conjugate transpose for matrices). Note that the dynamics of the particle configurations $d\boldsymbol{Q}/dt$ can be directly expressed in terms of $\boldsymbol{U}$ using the quaternion representation of the particle orientations, as discussed in BrownianMultiBlobs .

Numerical methods for temporal integration of (9) in time are discussed in BrownianMultiblobSuspensions . These methods require an efficient method to generate both the deterministic and Brownian velocities of the particles. Specifically, efficient Brownian dynamics for rigid bodies requires efficiently computing an Euler-Maruyama approximation of the apparent linear and angular velocities over a time step of duration $\Delta t$,

[TABLE]

where $\boldsymbol{W}$ is a collection of independent standard Gaussian random variables. Much of the strength and flexibility of the RB-FIB method we introduce in this work comes from the fact that we compute $\boldsymbol{U}_{\text{EM}}$ by discretizing a semi–continuum formulation of the overdamped particle dynamics (temporarily neglecting the thermal drift) instead of relying on explicit representations of Greens functions as in BrownianMultiblobSuspensions ; FBIM . The missing stochastic drift term in (9) involving $\partial_{\boldsymbol{Q}}\cdot\boldsymbol{\mathcal{N}}$ can be obtained in expectation by adding a correction to the Euler-Maruyuama method, as we explain in detail in Section IV.

Temporarily neglecting the terms that contribute to the stochastic drift, the overdamped limit of equations (3)-(4), (2), and (1) can be obtained quite simply by deleting all of the inertial terms and replacing the space-time white noise field $\boldsymbol{\mathcal{Z}}(\boldsymbol{x},t)$ with a white-in-space random Gaussian field $\boldsymbol{Z}(\boldsymbol{x})/\sqrt{\Delta t}$. In the immersed boundary approach we use here, we extend the fluid equation over the whole domain, including inside the bodies, since (2) is satisfied for all $\boldsymbol{r}\in\mathcal{B}_{p}$ including the rigidly moving interior RigidMultiblobs . This gives a system of semi–continuum linear equations for $\boldsymbol{U}_{\text{EM}}=\left[\boldsymbol{u}^{\text{EM}}_{1},\boldsymbol{\omega}^{\text{EM}}_{1},\ldots\boldsymbol{u}^{\text{EM}}_{N_{b}},\boldsymbol{\omega}^{\text{EM}}_{N_{b}}\right]$,

[TABLE]

We do not add a superscript ‘EM’ to the velocity and pressure here with the understanding that in the overdamped limit they are just auxiliary variables used to obtain the motion of the particles. Here $\boldsymbol{\lambda}$ is the jump in the fluid stress across the boundary of the particles, $\boldsymbol{\lambda}\equiv\llbracket\boldsymbol{\Sigma}\cdot\boldsymbol{n}\rrbracket$, which can be simply identified as the traction force when $\breve{\boldsymbol{u}}=\boldsymbol{0}$; see Appendix A in RigidMultiblobs for an explanation why the same formulation works even when the apparent slip $\breve{\boldsymbol{u}}$ is nonzero. The use of a Dirac–delta distribution to restrict quantities to their values on the surface of the bodies $\boldsymbol{r}\in\partial\mathcal{B}_{p}$ is the basis for the immersed boundary spatial discretization described next.

III Discrete Formulation

In this section we describe an efficient and robust means of discretizing the problem formulated in section II. We will use standard finite difference and immersed boundary methods to construct matrix discretizations of the differential and integral operators appearing in (11)–(14). We will briefly review efficient methods to solve the large-scale linear system which arises in the discretized equations; details can be found in RigidIBM ; RigidMultiblobs .

III.1 Discrete Fluctuating Stokes Equations

To begin discretizing equations (11)–(14) we consider the discretization of the fluctuating Stokes equations without any suspended particles. Specifically, we will temporarily ignore the effect of the rigid bodies on the fluid in equations (11)-(12) by setting $\boldsymbol{\lambda}=0$. Importantly, all of the discussion in this section will remain agnostic to the choice of physical boundary conditions on $\Omega$.

We discretize the Stokes equations on a regular Cartesian ‘Eulerian’ grid of cells with volume $\Delta V=h^{d}$, where $h$ is the grid spacing and $d$ is the space dimension. Velocity variables are staggered on the faces of the grid relative to the cell-centered pressure variables. The infinite dimensional white noise field $\boldsymbol{Z}$ is spatially discretized as $\boldsymbol{W}$, a collection of random Gaussian variables generated on the faces and nodes of the fluid grid, as described in LLNS_Staggered . We use staggered discretizations of the vector divergence $\boldsymbol{D}$ and scalar gradient $\boldsymbol{G}$ operators that obey the adjoint relation $\boldsymbol{G}=-{\boldsymbol{D}}^{\star}$. We also define discrete tensor divergence $\boldsymbol{\mathds{D}}$ and vector gradient $\boldsymbol{\mathds{G}}=-\boldsymbol{\mathds{D}}^{\star}$ operators, which account for the imposed boundary conditions on $\partial\Omega$. The discrete scalar Laplacian is $\boldsymbol{L}=\boldsymbol{D}\boldsymbol{G}$ and the discrete vector Laplacian is $\boldsymbol{\mathbb{L}}=\boldsymbol{\mathds{D}}\boldsymbol{\mathds{G}}$ in order to satisfy a discrete fluctuation dissipation balance principle LLNS_Staggered .

With $\boldsymbol{\lambda}=0$, we discretize equations (11)-(12) as

[TABLE]

This maintains discrete fluctuation dissipation balance even in the presence of physical boundaries. It is convenient at this point to introduce the symmetric positive-semidefinite discrete Stokes solution operator

[TABLE]

such that the solution to (15) can be written as $\boldsymbol{v}=\boldsymbol{\mathcal{L}}^{-1}\boldsymbol{f}.$

III.2 Immersed Boundary Method for Rigid Bodies

To mediate the fluid-structure interaction we use a discrete approximation $\delta_{h}$ to the Dirac delta distribution appearing in (11)–(14). In the results presented here $\delta_{h}$ is the 6-point kernel developed in New6ptKernel . While the discretization on the grid introduces numerical artifacts in general IBM_PeskinReview , the 6-point kernel we use here was specifically designed with grid invariance and smoothness in mind New6ptKernel .

Throughout the rest of this work, we will represent a body $\mathcal{B}_{p}$ as a rigid agglomerate of markers or blobs with positions $\boldsymbol{r}^{p}_{i}\in\partial\mathcal{B}_{p}$ and we refer to the collection of these points as the ‘Lagrangian’ grid. The ideal spacing between the blobs $s\sim h$ is related to the meshwidth $h$ used in discretizing the fluid equations (15), as discussed in detail in section IV of RigidMultiblobs . We discretize $\boldsymbol{\lambda}$ on the Lagrangian grid as a collection of force vectors $\boldsymbol{\lambda}_{i}^{p}\approx\boldsymbol{\lambda}\left(\boldsymbol{r}^{p}_{i}\right)\Delta A\left(\boldsymbol{r}^{p}_{i}\right)$. It is important to note that the discrete $\boldsymbol{\lambda}_{i}^{p}$ has units of force rather than force density as the continuum $\boldsymbol{\lambda}\left(\boldsymbol{r}\right)$. Recall that the fluid velocity $\boldsymbol{v}$ is defined on the centers of the faces $\boldsymbol{x}_{\alpha}$ of the Eulerian grid.

To discretize equations (11)–(14) we use simple trapezoidal rule quadratures to approximate the integrals as appropriate sums over Eulerian or Lagrangian grid points. This leads us to define the spreading operator $\boldsymbol{\mathcal{S}}$ and the interpolation operator $\boldsymbol{\mathcal{J}}$ as

[TABLE]

It is important to note that $\boldsymbol{\mathcal{J}}$ and $\boldsymbol{\mathcal{S}}$ satisfy the adjoint relation $\boldsymbol{\mathcal{J}}=\Delta V\boldsymbol{\mathcal{S}}^{\star}$. This property ensures discrete conservation of energy IBM_PeskinReview as well as discrete fluctuation dissipation balance. The definitions of $\boldsymbol{\mathcal{S}}$ and $\boldsymbol{\mathcal{J}}$ are modified when the support of the kernel overlaps with a physical boundary of $\Omega$. Appropriate ghost points across the physical boundary are used which ensures the effects of the boundary are incorporated into the spreading and interpolation operations while preserving the adjoint property, as described in Appendix D of RigidIBM .

We discretize the integrals giving the total force and torque in (14) using simple trapezoidal quadrature to define the geometric matrix $\boldsymbol{\mathcal{K}}(\boldsymbol{Q})$ RigidMultiblobs_Swan ,

[TABLE]

III.3 The Discrete System

We can now compactly state the spatially discretized system (11)–(14) as

[TABLE]

For simplicity in the following discussion, we will take the fluid body force $\boldsymbol{g}=\boldsymbol{0}$.

Using $\boldsymbol{\mathcal{S}}$ and $\boldsymbol{\mathcal{J}}$ we may define a regularized, symmetric, positive semi-definite blob-blob mobility matrix $\boldsymbol{\mathcal{M}}=\boldsymbol{\mathcal{J}}\boldsymbol{\mathcal{L}}^{-1}\boldsymbol{\mathcal{S}}$, where we recall that $\boldsymbol{\mathcal{L}}^{-1}$ denotes the discrete Stokes solution operator. The block $\boldsymbol{\mathcal{M}}_{ij}$ gives the pairwise mobility matrix between two blobs $\boldsymbol{r}_{i}$ and $\boldsymbol{r}_{j}$ BrownianBlobs ; RigidIBM ; RigidMultiblobs ,

[TABLE]

where $\mathbb{G}$ is the Green’s function for the Stokes equations in $\Omega$ with the specified boundary conditions. For two markers/blobs that are sufficiently far apart in a sufficiently large domain, $\boldsymbol{\mathcal{M}}_{ij}$ approximates the mobility for a pair of spheres of radius 222Specifically, $a=1.47h$ for the 6-point kernel used in this work RigidIBM ; RigidMultiblobs . $a\sim h$. When the Green’s function is available analytically, $\boldsymbol{\mathcal{M}}$ can be computed explicitly; here we handle more general boundary conditions by solving the discretized steady Stokes equations numerically to compute the action of $\mathbb{G}$.

We may use the definition of $\boldsymbol{\mathcal{M}}$ to eliminate the velocity and pressure from (20)–(23) to obtain the reduced linear system

[TABLE]

Here we can identify the matrix $\boldsymbol{\mathcal{M}}^{1/2}=\sqrt{\eta/\Delta V}\boldsymbol{\mathcal{J}}\boldsymbol{\mathcal{L}}^{-1}\boldsymbol{\mathds{D}}$ as was done in BrownianBlobs , such that $\boldsymbol{\mathcal{M}}^{1/2}\left(\boldsymbol{\mathcal{M}}^{1/2}\right)^{\star}=\boldsymbol{\mathcal{M}}$. Equation (25)–(26) is identical to equation (9) in BrownianMultiblobSuspensions – the only difference here is that we do not explicitly have access to $\boldsymbol{\mathcal{M}}$ because we do not necessarily know the Green’s function $\mathbb{G}$. Following BrownianMultiblobSuspensions , the solution of (25)–(26) can be written in terms of the body mobility matrix $\boldsymbol{\mathcal{N}}=\left(\boldsymbol{\mathcal{K}}^{T}\boldsymbol{\mathcal{M}}^{-1}\boldsymbol{\mathcal{K}}\right)^{-1}$,

[TABLE]

where we have identified $\boldsymbol{\mathcal{N}}^{1/2}\equiv\boldsymbol{\mathcal{N}}\boldsymbol{\mathcal{K}}^{T}\boldsymbol{\mathcal{M}}^{-1}\boldsymbol{\mathcal{M}}^{1/2}$, such that $\boldsymbol{\mathcal{N}}^{1/2}\left(\boldsymbol{\mathcal{N}}^{1/2}\right)^{\star}=\boldsymbol{\mathcal{N}}$.

We have now shown that we can compute $\boldsymbol{U}_{\text{EM}}$ given by (27) efficiently by solving the linear system (11)–(14), which can be done with a complexity linear in the number of particles thanks to the preconditioned GMRES solver developed in RigidMultiblobs . All the remains to perform Brownian Dynamics is an efficient means of computing the stochastic drift term $k_{B}T\left(\partial_{\boldsymbol{Q}}\cdot\boldsymbol{\mathcal{N}}\right)$ in (9), as we discuss next.

IV Temporal Integration

In this section, we develop a temporal integration scheme for equation (9) which efficiently captures the contribution from the stochastic drift term $k_{B}T\left(\partial_{\boldsymbol{Q}}\cdot\boldsymbol{\mathcal{N}}\right)$. Specifically, we use random finite differences (RFDs) BrownianBlobs ; BrownianMultiBlobs ; MagneticRollers to compute terms that will be included on the right hand side of equations (20)–(23) to account for the stochastic drift in expectation. The algorithm we develop here requires the solution of an additional linear system similar to (20)–(23) each time step in order to capture the stochastic drift. This is still more efficient than the classical Fixman scheme which requires solving a resistance problem, which is a lot more expensive than solving a mobility problem when iterative methods are used RigidMultiblobs ; libStokes .

IV.1 Random Finite Differences

Solving equations (20)–(23) without any modifications yields an efficient means of computing $\boldsymbol{U}_{\text{EM}}$ and all that remains to simulate equation (9) is a means of computing $k_{B}T\left(\partial_{\boldsymbol{Q}}\cdot\boldsymbol{\mathcal{N}}\right)$. In past work BrownianBlobs ; BrownianMultiBlobs ; MagneticRollers , some of us have proposed to use random finite differences (RFD) to generate this term as follows. Consider two Gaussian random vectors $\Delta\boldsymbol{P}$ and $\Delta\boldsymbol{Q}$, such that $\left\langle\Delta\boldsymbol{P}\right\rangle=\left\langle\Delta\boldsymbol{Q}\right\rangle=0$ and $\left\langle\Delta\boldsymbol{P}\Delta\boldsymbol{Q}^{T}\right\rangle=\boldsymbol{I}$. For an arbitrary matrix $\boldsymbol{\mathcal{R}}\left(\boldsymbol{Q}\right)$, it holds that

[TABLE]

where $\left\langle\cdot\right\rangle$ denotes an average over realizations of the random vectors.

The limits in equations (28) may be discretely approximated by simply choosing a small value for $\delta$ at the cost of introducing a truncation error of $\mathcal{O}\left(\delta\right)$ if the one–sided difference is used (first line) or $\mathcal{O}\left(\delta^{2}\right)$ if the centered difference is used (second line). The value of $\delta$ should be chosen to balance the magnitude of the truncation errors in the RFD with any numerical error associated with the application of (multiplication with) $\boldsymbol{R}$. Specifically, we take $\delta=\epsilon^{1/2}$ for the one-sided difference and $\delta=\epsilon^{1/3}$ for the centered difference, where $\epsilon$ is the relative error with which matrix-vector products with $\boldsymbol{R}$ are computed.

A simple modification of the Euler–Maruyama scheme is to use a RFD with $\boldsymbol{\mathcal{R}}\equiv\boldsymbol{\mathcal{N}}$ to account for the stochastic drift BrownianMultiBlobs . Since application of $\boldsymbol{\mathcal{N}}$ requires solving (20)–(23) iteratively with some loose relative tolerance $\epsilon$, the required value of $\delta$ would be relatively large especially for the one-sided difference. Using the two-sided difference requires solving two additional mobility problems per time step, which is quite expensive. We now propose an alternative approach.

IV.2 The Split–Euler–Maruyama (SEM) Scheme

Our goal to design a means of computing the stochastic drift with as few linear solves as possible. In our prior work BrownianMultiblobSuspensions , we accomplished this by expanding $\partial_{\boldsymbol{Q}}\cdot\boldsymbol{\mathcal{N}}$ using the chain rule,

[TABLE]

In this work we do not have explicit access to $\boldsymbol{\mathcal{M}}$ so we carry this expansion one step further as done in BrownianBlobs ,

[TABLE]

Our aim is to generate each term in (IV.2) and (31) through a separate RFD on $\boldsymbol{\mathcal{J}}$, $\boldsymbol{\mathcal{S}}$, $\boldsymbol{\mathcal{K}}$, or $\boldsymbol{\mathcal{K}}^{T}$. This is particularly advantageous as these operators can all be applied efficiently to within roundoff tolerance in linear time, without requiring linear solvers.

The Euler–Maruyama–Traction (EMT) scheme proposed in BrownianMultiblobSuspensions can be adapted to the present context to give the Split–Euler–Maruyama (SEM) scheme outlined in Algorithm 1. In the first step of this algorithm we generate random forces and torques on each body

[TABLE]

where $\boldsymbol{W}_{p}^{\text{f}/\tau}$ are standard Gaussian random variables generated independently for each body $p$. Here $L_{p}$ is a length scale for body $p$, which we take as the maximum pairwise distance between blobs on a body. In step 2 of Algorithm 1 we solve a mobility problem with $\boldsymbol{W}^{\text{FT}}$ as the applied force/torque on each body, to obtain the random variables

[TABLE]

Defining a random translational and rotational displacement for body $p$

[TABLE]

gives two randomly diplaced positions 333Here for simplicity of notation we use addition to denote a random rotation of the body by an oriented angle $(\delta/2)\boldsymbol{W}_{p}^{\tau}$ even though in practice this is realized as a quaternion multiplication in three dimensions BrownianMultiBlobs . of body $p$,

[TABLE]

We can produce the desired drift term through the following random finite differences on the matrices $\boldsymbol{\mathcal{K}}^{T}$ and $\boldsymbol{\mathcal{K}}$:

[TABLE]

which are analogous to the quantities $\boldsymbol{D}^{F}$ and $\boldsymbol{D}^{S}$ computed in Algorithm 1 of BrownianMultiblobSuspensions .

However, the RFD performed directly on $\boldsymbol{\mathcal{M}}$ to compute $\boldsymbol{D}^{S}$ in the EMT scheme (see Algorithm 1 of BrownianMultiblobSuspensions ) is computed in the SEM scheme using (31) as a sum of RFDs on the interpolation and spreading operators,

[TABLE]

Note that the computation of $\boldsymbol{D}^{\boldsymbol{\mathcal{S}}}$ in step 4 of Algorithm 1 requires an additional application of $\boldsymbol{\mathcal{L}}^{-1}$ and therefore an additional unconstrained Stokes solve. This does not add much additional complexity to the computation because unconstrained (fluid only) Stokes systems of the form (15) have fewer degrees of freedom and are far better conditioned than constrained (fluid + rigid bodies) systems of the form (20)–(23) RigidMultiblobs .

To produce the correct drift term, we add $\boldsymbol{D}^{\boldsymbol{\mathcal{K}}}-\boldsymbol{D}^{\boldsymbol{\mathcal{J}}}-\boldsymbol{D}^{\boldsymbol{\mathcal{S}}}$ as a random slip and add $\boldsymbol{D}^{\boldsymbol{\mathcal{K}}^{T}}$ as a random force on the right hand side of the linear system in step 6 of Algorithm 1. This generates an additional contribution to the velocity of the rigid particles, $\boldsymbol{U}^{n}=\boldsymbol{U}_{\text{EM}}+\boldsymbol{U}_{\text{Drift}}$, where

[TABLE]

We used equations (38)–(41) as well as equation (IV.2) to simplify from the first to the second line in (IV.2). On average (i.e. in expectation) this will produce the desired drift term,

[TABLE]

as shown in more detail in Appendix B.

V Numerical Results

To numerically investigate the RB-FIB algorithm, we have implemented it in the IBAMR code IBAMR , freely available at https://github.com/IBAMR, by modifying existing codes developed for deterministic Stokesian suspensions in RigidMultiblobs . In particular, we have reused the existing linear solvers in steps 2 and 6 of Algorithm 1. The scaling and convergence of the numerical linear algebra routines, as well as optimal parameters for the outer and inner Krylov and multigrid iterative solvers are discussed in RigidMultiblobs . In all of the following numerical examples (including the Appendix), we use an absolute tolerance proportional to the time step size in the outer FGMRES solver required by steps 2 and 6 of the SEM scheme, as well as in the unconstrained GMRES solve required by step 4. As recommended for the 6-point immersed boundary kernel in RigidMultiblobs , we take the spacing of the Lagrangian blobs $s\approx 3h$, where $h$ is the Eulerian grid spacing. Time steps which generate unphysical configurations (such as a blob overlapping the wall) are rejected and the step repeated. However, these instances are very rare because we use repulsive potentials to prevent particle-particle and particle-wall overlaps, and because we employ the modifications of $\boldsymbol{\mathcal{S}}$ and $\boldsymbol{\mathcal{J}}$ near the boundaries introduced in Appendix D in RigidIBM . In all of the following examples, the fluid is water at room temperature $T=300$ K and viscosity $\eta=1$ mPa s.

Data for the examples studied in this section was gathered on Northwestern University’s QUEST computing cluster. Multiple independent trajectories are run for each case considered in both of the following examples as well as for Appendix A, where we take as many time steps as would complete in a fixed amount of computation time (typically on the order of one week). We accounted for the variable lengths of the runs by using means and variances weighted by the trajectory length when computing relevant statistics. The error bars included in the figures of this section show $95\%$ confidence intervals, computed using the weighted means and variances.

In Appendix A we examine a single colloidal boomerang suspended in a slit channel to demonstrate the first order weak accuracy of the SEM scheme and validate our implementation. This appendix also demonstrates the viability of the RB-FIB method for arbitrary particle shapes and shows the importance of numerically capturing the stochastic drift term in (9).

In Sections V.1 and V.2 we numerically investigate how tight physical confinement affects hydrodynamic coupling between particles. In section V.1 we use a simple example of two spheres confined in a cuboid in order to establish the spatial resolution required for the RB-FIB method to capture dynamic statistics with sufficient accuracy. In section V.2 we numerically investigate a variant of an experiment reported in Bohlein2011 that measured the friction forces and resulting wave patterns in a hydrodynamically driven colloidal monolayer. In our setup the monolayer is confined in a narrow slit channel and is also confined in the lateral directions by walls, mimicking an experiment performed in a microfluidic channel. Using the RB-FIB method we find some novel behavior in the propagation of the waves through the mololayer.

V.1 Two Spheres in a Tight Cavity

In this section we investigate how spatial resolution effects the accuracy of the RB-FIB method. We simulate two neutrally buoyant spheres of hydrodynamic radius $R_{h}=0.656\ \mu$m, tightly confined in a $3.478\times 1.739\times 1.739\ \mu\text{m}^{3}$ rectangular box. These physical dimensions ensure that the spheres are almost always in near contact with each other and/or a physical boundary, allowing us to highlight that the RB-FIB method can tackle problems in which confinement plays an important role in the dynamics.

To reduce the severe time step size restriction required to ensure that there are no particle-particle or particle-wall overlaps, we introduce a soft repulsive potential between the two particles as well as between the particles and the wall of the form

[TABLE]

For the interparticle repulsion, $d=2R_{H}$ and $r$ is the distance between the particle centers. For the repulsion between a particle and a wall, $d=R_{H}$ and $r$ is the distance between a particle center and a wall. We take $b=0.1R_{H}$ and $\Phi_{0}=4k_{B}T$ for both the particle-particle and particle-wall potentials as in MagneticRollers , as this choice ensures that the time scale associated with the steric repulsion isn’t much smaller than the diffusive time scale, while also maintaining a low probability of unphysical configurations. We use a dimensionless time step size $\Delta\tau=\frac{k_{B}T}{6\pi\eta R_{h}^{3}}\Delta t=0.0044$ as this was found to be small enough to ensure that the temporal integration errors are smaller than the statistical errors.

We discretize the domain using a grid spacing $h=0.0543\times(1,2,4)\mu$m and discretize the particles with $162,42,12$ approximately equally spaced blobs respectively, as shown in figure 1. The geometric radius of the particles is determined according to table I in RigidMultiblobs in order to maintain a constant hydrodynamic radius $R_{h}=0.656\mu$m as the resolution is refined. The physical domain is taken to be of dimensions $2L\times L\times L$ where $L=8\times 0.0543\mu$m.

In appendix A we study the temporal integration errors by examining the marginals of the equilibrium (static) Gibbs–Boltzmann (GB) distribution for a boomerang confined in a slit channel; we performed similar tests for the two-sphere example studied here and found negligible errors in the static distributions for all resolutions. This is expected for any equilibrium average for sufficiently small time step sizes because only dynamic quantities are affected by the resolution of the hydrodynamics. Here we investigate dynamic statistics in the form of the equilibrium mean squared displacement (MSD) block matrix with tensor blocks

[TABLE]

where the average is taken over equilibrium trajectories and the subscript $p,r=1,2$ denotes the particle. To compute the components of the MSD we use the SEM scheme to simulate several independent equilibrium trajectories for each of the three resolutions. In Fig. 2 we compare the components of the MSD for different spatial resolutions. Because of the symmetry in the problem, $\boldsymbol{\text{MSD}}_{11}=\boldsymbol{\text{MSD}}_{22}$ and $\boldsymbol{\text{MSD}}_{pp}^{yy}=\boldsymbol{\text{MSD}}_{pp}^{zz}\equiv\boldsymbol{\text{MSD}}_{pp}^{\perp}$, so we show the average of the equivalent components of the MSD tensor.

Because the particles are completely confined, every component of the translational MSD approaches an asymptote

[TABLE]

where $\boldsymbol{z}^{p}_{1},\boldsymbol{z}^{p}_{2}$ are two independent samples from the equilibrium distribution of particle $p$, generated from an MCMC method. The appropriate asymptotes are plotted in Figs. 2(a),(b) and we can see that the correct asymptotic MSD is approached using all three resolutions considered. This is again expected because the asymptotic MSD is controlled by the Gibbs-Boltzmann equilibrium distribution and not by the (hydro)dynamics.

The Stokes–Einstein relation gives the short time diffusion tensor

[TABLE]

where the superscript in $\left\langle\boldsymbol{\mathcal{N}}_{pr}^{(tt)}\right\rangle_{GB}$ refers to the translation–translation block of the body mobility matrix $\boldsymbol{\mathcal{N}}$ and $\left\langle\cdot\right\rangle_{GB}$ denotes an average with respect to the Gibbs–Boltzmann (GB) distribution. To compute $\left\langle\boldsymbol{\mathcal{N}}\right\rangle_{GB}$ we use $642$ blobs to discretize each sphere ($h=0.5\times 0.0543\mu$m) and compute a sample mean of $\boldsymbol{\mathcal{N}}$ over equilibrium configurations sampled using a Markov Chain Monte Carlo (MCMC) method 444We generate enough samples to measure each component of $\left\langle\boldsymbol{\mathcal{N}}_{pr}^{(tt)}\right\rangle_{GB}$ to within $1\%$ statistical error with $95\%$ confidence.. The insets of figure 2 show that the short-time Stokes–Einstein relation (47) is accurately maintained for both of the finer resolutions but not for the coarsest resolution. Theoretical results are unavailable for the self-diffusion coefficient at intermediate times. We see in Figs. 2(a),(b) close agreement between the two higher resolutions ($162$ and $42$ blobs), but with some visible deviations for the lowest resolution case ($12$ blobs), indicating insufficient resolution.

The authors of TwoSphereCouple found analytically and experimentally that two nearby colloidal spheres were strongly hydrodynamically coupled but the coupling weakens significantly near a confining boundary. In our example we see a competition of influence: the two spheres are always close to each other and hence their motion should be coupled, but they are also always close to a physical boundary which would decouple their motion. To get a clearer picture of this competition, Fig. 3 shows the non-vanishing components (the $xx$ or $\parallel$ and $yy,zz$ or $\perp$ components) of $\boldsymbol{\text{MSD}}_{1,2}$. The inset of Fig. 3(a) shows that the short time diffusive motion in the $\parallel$ direction between the two spheres (given by (47)) is very weakly coupled. However, the full panel of Fig. 3(a) shows a much stronger coupling in the parallel motion of the spheres for longer times. Therefore, the influence of another nearby particle eventually dominates over the influence of the walls, which initially decouples the particles’ motion. In Fig. 3(b) we see that the short time diffusion in the $\perp$ ($y,z$) directions between the spheres is somewhat strongly anti-correlated. After $t\approx 4$s the behavior of the MSD inflects and decays to nearly zero for larger times and the perpendicular motion of the spheres effectively decouples due to the confinement. The results in Fig. 3 indicate that 12 blobs per sphere is not sufficient to accurately predict the time correlation functions for more than one particle, and at least 42 blobs per sphere are required for particles this close to each other and the walls.

V.2 Friction in a Colloidal Monolayer

In Bohlein2011 , the authors performed an experiment in which a colloidal monolayer was hydrodynamically driven across a bottom wall on which a substrate potential was generated by optical traps. This potential mimics the effect of corrugation of the wall, which is a key contributor to the effective friction with the wall. In Bohlein2011 the system was kept quasi–two–dimensional by forcing the monolayer to remain near the wall using a vertically incident laser to form a confining potential that dramatically reduced out of plane motion. The colloids used in the experiment were negatively charged polystyrene spheres suspended in water. Due to their negative charge, the particles spontaneously formed a stable 2D triangular crystal YukawaLattice20KT . The corrugation potential used in the experiments reported in Bohlein2011 took the form of a periodic lattice with 3-fold symmetry around its minima. The colloidal crystal and the corrugation potential are called commensurate if the lattice constants agree and incommensurate otherwise; both cases were considered in Bohlein2011 .

In the experiments reported in Bohlein2011 , the sample cell was translated to generate a flow field and hence a fluid drag force on the colloids. This lateral drag force on the crystal served as a control parameter and, under commensurate conditions, the authors observed a critical translation velocity (and hence lateral force) at which the colloidal crystal became unpinned from the corrugation potential. This critical value represents the transition of the crystal from static friction to kinetic as it becomes unpinned. Just above this critical velocity, they observed localized density variations in the colloidal crystal taking the form of traveling kink solitons.

In this section, we will use the RB-FIB method to numerically investigate a novel variant of this experiment wherein the quasi–two–dimensionality of the problem (essential for the existence of a stable 2D monolayer) is achieved by confining the monolayer in a slit channel between two walls, see illustration in figure 4. This type of confinement is easy to realize in experiments using a microchannel. Additionally, this example serves to demonstrate an important feature of the RB-FIB method: not only can it account for tight confinement, but it can also easily account for external flow fields generated by the imposed boundary conditions.

V.2.1 Physical Parameters

The electrostatic repulsion between the colloids is accounted for by a pairwise Yukawa potential of the form

[TABLE]

where $U_{0}$ is the repulsion strength, $\lambda_{D}$ is the Debye screening length, and $r$ is the distance between particle centers. We set the screening length $\lambda_{D}=0.16\mu$m to the value reported in Bohlein2011 . The repulsion strength $U_{0}$ is a free parameter which must be chosen to ensure that a stable colloidal crystal is formed in the absence of a trapping potential, and that this crystal remains stable in the presence of the thermal fluctuations and background flow. In agreement with YukawaLattice20KT ; YukawaLattice80KT , we find that choosing $20k_{B}T\leq U_{0}\leq 80k_{B}T$ is reasonable as the results presented later in this section were not strongly effected by taking $U_{0}$ to be either of these extremes; we will fix $U_{0}=20k_{B}T$ henceforth. To reduce the probability of blob–wall overlaps we include a soft repulsive potential between each blob and all of the physical boundaries of the domain. This potential takes the form of (44), where $d=a$, $b=0.1a$, and $\Phi_{0}=4k_{B}T$, and $r$ is the distance between a blob and a wall.555For spherical particles we could have put a wall-repulsive potential on each sphere’s center rather than each blob and avoided small spurious torques on the particles. However this would not generalize easily to arbitrary particle shapes. We use a dimensionless steric time step size $\Delta\tau=\Delta t\left(\Phi_{0}/(6\pi\eta a^{2}b)\right)=0.138$ which we found to be sufficiently small to make particle–wall overlaps infrequent.

To capture the effects of corrugation, each colloid feels the potential ColloidSheetFriction

[TABLE]

where a particle’s $y$–symmetrized 2D position is $\boldsymbol{X}_{2D}=[x,y-L_{y}/2]^{T}$. Here the scaled lattice directions are

[TABLE]

Following Bohlein2011 , we take the average particle separation to be $s_{\text{Lattice}}=5.7\mu$m and the particle radius $R_{H}=1.95$. In the absence of a substrate potential the spacing of the colloidal crystal is controlled by the number of particles and the domain dimensions. To ensure that the spacing of the substrate potential $s_{\text{Lattice}}$ is commensurate with the spacing in the colloidal crystal, we use 272 particles and take $L_{x}=16s_{\text{Lattice}}$, $L_{y}=15s_{\text{Lattice}}$. The domain is taken to be periodic in the direction of the applied flow field ($x$ direction) and no–slip boundaries are used in every other direction. While the periodicity in the $x$ direction introduces some unphysical artifacts, $L_{x}$ is large enough to produce kink solitons, which were found to have a support of $\approx 8s_{\text{Lattice}}$ Bohlein2011 . The width of the domain in the $z$ direction is taken to be $L_{z}=1.28D_{H}$, where $D_{H}=2R_{H}$ is the particle diameter. This is to ensure the fairly strict quasi–two dimensionality required for colloidal crystals to form.

To drive the colloidal monolayer, we move the top, bottom and side walls with velocity $v_{\text{wall}}$ along the $x$ axes, i.e., we impose a fluid velocity $\left(v_{\text{wall}},0,0\right)$ on all walls as a boundary condition for the Stokes equations 666In this simple case, the flow created by the wall slip is simple constant plug flow, so one could move/shift the corrugation potential with velocity $\left(-v_{\text{wall}},0,0\right)$ instead of imposing a velocity on the walls. However, in more general situations, e.g., (time-dependent) pressure-driven flow in square channels, it is much simpler and more flexible to let the Stokes solver compute the fluid flow generated by the boundaries.. We non–dimensionalize the control parameter $v_{\text{wall}}$ using the work required to move one colloid one lattice site in an unbounded domain,

[TABLE]

In what follows, we vary $W_{\text{wall}}=\left\{1.63,1.86,2.1,2.33,2.56\right\}$ and investigate the dynamics of the colloidal monolayer shown in figure 4.

V.2.2 Discretization Parameters

We discretize the colloids using 42 blobs, as this resolution was found in section V.1 to provide sufficient accuracy for the dynamics. The fluid grid spacing $h=s_{\text{Lattice}}/16$ so that substrate potential aligns with the periodic domain size in the $x$ direction. The blob spacing $s=0.95$ was chosen so that the particles’ hydrodynamic radius would be $R_{H}=1.95$ (geometric radius $R_{G}=1.7380$) and $s\approx 3h$ as was recommended in RigidMultiblobs .

V.2.3 Numerical Observations of Kinks in a Colloidal Monolayer

In this section we vary the dimensionless wall velocity $W_{\text{wall}}$ to investigate the critical value at which static friction is broken and the colloidal monolayer begins to slide through the thin ’corrugated’ channel via kink solitons. To measure the point at which the colloidal monolayer breaks from the corrugated substrate potential, we simply determine whether any particle displaces more than $s_{\text{Lattice}}$ from its initial position. Partial depinning of the colloidal crystal is observed at $W_{\text{wall}}=1.86$ wherein only some of the particles are sufficiently displaced. Full depinning is observed for all $W_{\text{wall}}>1.86$ wherein all of the particles are sufficiently displaced after some amount of time. Brownian motion is crucial to the formation and propagation of kinks. While it is certainly possible to see depinning of the monolayer in the absence of thermal fluctuations, the tiny ’kicks’ provided by the fluctuating fluid activate a particle’s transition between potential wells. So much so, in fact, that no depinning was observed in complimentary deterministic simulations for all values of $W_{\text{wall}}$ considered.

As noted in Bohlein2011 , a kink is formed when one particle escapes the potential well it is confined to through a combination of thermal forces as well as hydrodynamic drag from the background flow. Once a particle escapes, it enters the neighboring potential well in the direction of the flow. Once confined to this new well, which also typically has a particle trapped in it, a combination of Yukawa and steric repulsion forces the original occupant of the well into its neighboring well, and the process repeats. We identify as a kink the propagation of particles escaping their well along the direction of flow.

Figure 5 shows the propagation of kinks in three nearby rows of particles along the direction of flow; see the SI for a movie. Figure 5(b) shows the $x$ position of each particle from the ‘blue’ row over time. The transverse black contour lines show local extrema in the velocities of the particles (averaged over 50 time steps to filter out the Brownian velocities), which are seen to correspond to the times at which a particle jumps to another lattice site. The bends in the black velocity contours show a finite speed of propagation for the kink and the ‘S’ shaped profile is due to the periodicity in the $x$ direction. Figure 5(c) shows the contours of the maximum velocity for the particle rows highlighted in panel (a). The prevailing ‘S’ shape in all of the contours in panel (c) demonstrates that there are propagating kinks in all of the rows of particles in the monolayer. However, it is difficult to appreciate the correlations in the motion of the kinks in nearby rows in these results.

In figure 6 we investigate how kinks influence other kinks in the direction transverse to their motion. The authors of Bohlein2011 observed that kinks extend also perpendicular to the direction of the force with a small lag in time. Their system however, was large enough to be considered unconfined in the $y$ direction. Here we examine the correlation of kinks between rows of particles, where a row is defined by binning the particles’ $y$ coordinates with bin width $2s_{\text{Lattice}}$. By selecting a representative particle from each row in the monolayer, as seen in the top panels of figure 6, we can track the hopping or ‘displacement’ times for each row. That is, we track the time at which the representative particle of a row has been displaced by an integer multiple of $s_{\text{Lattice}}$. The bottom row of panels in Fig. 6 shows the hopping times for representative particles shown in the top row of panels, where the color of the markers corresponds to the color of the particle. The displacement times of the particles are grouped based on how many lattice sites the particles have displaced. That is, every particle’s first displacement time is connected, as is every particle’s second, and so on.

The bottom row of panels in Fig. 6 shows that for every value of $W_{\text{wall}}>1.63$, the first rows to be displaced are the rows closest to the walls which bound the domain in the $y$ direction. For $W_{\text{wall}}=1.87$, the only particles which become displaced are those in the rows nearest these walls and their immediate neighbors. This is likely because of the additional hydrodynamic screening provided by the walls, as we investigated in section V.1. By analogy with what was observed in Bohlein2011 , one might expect that the next rows to be displaced are the immediate neighbors of the rows nearest the walls, and then their neighboring rows, and so forth, with the middle rows displacing last. To the contrary, we see that one of the middle rows displaces soon after the displacement of the rows nearest the walls. This is a surprising result and may stem from the lateral ($y$) confinement of the system. A more thorough investigation of the unusual collective dynamics of kinks in this system is deferred to future work.

VI Conclusion

In this work, we described the Rigid Body Fluctuating Immersed Boundary (RB-FIB) method to simulate the Brownian dynamics of arbitrarily shaped rigid particles in fully confined domains. The fluctuating solvent was treated explicitly, allowing for well-known staggered finite difference schemes to be used and general boundary conditions to be applied on the boundaries of the computational domain. We designed an efficient Split–Euler–Maruyama temporal integrator that uses a nontrivial combination of random finite differences to capture the stochastic drift appearing in the overdamped Langevin equation. We implemented this method in the IBAMR software infrastructure, which is freely available at https://github.com/IBAMR.

We studied the dynamical correlation functions of two tightly confined spheres in close proximity to each other and physical walls, and examined the effect of spatial resolution on the accuracy of both the short- and long-time equilibrium mean-square displacement. We also used the RB-FIB method to model a quasi–2D colloidal crystal confined in a narrow slit channel. The layer was hydrodynamically driven across a commensurate periodic substrate potential mimicking the effect of a corrugated wall. We observed partial and full depinning of the colloidal monolayer from the substrate potential above a certain wall speed, consistent with the transition from static to kinetic friction observed in Bohlein2011 . Further, we observed the propagation of kink solitons parallel to the direction of flow. These kinks extended along the colloidal monolayer in the direction transverse to the flow. We observed a curious ‘M’ pattern in the particle displacements across the domain wherein particles nearest the boundaries of the domain and particles in the middle of the domain are the first to be displaced.

The SEM scheme presented here is based on the Euler–Maruyama method and as such is first order weakly accurate (as shown both numerically and theoretically in Appendix A and B), as well as only first order deterministically accurate. Higher order deterministic accuracy has been shown to be very beneficial in designing temporal integrators for Brownian dynamics of rigid particles in half-space domains BrownianMultiblobSuspensions . The SEM scheme can be extended to deterministically second-order accurate Adams–Bashforth, trapezoidal and midpoint variants, but whether or not the additional computational cost incurred by these methods is justified should be investigated in future work.

The RB-FIB method we present here has notable advantages over other methods. First, it can handle a variety of combinations of boundary conditions in different directions seamlessly. Second, it can handle colloidal particles of complex shapes with varying levels of spatial resolution (fidelity), i.e., with controllable accuracy. Third, the method scales linearly in the number of particles and fluid grid cells. Fourth, the explicit solvent approach also facilitates coupling to additional physics, including elastic bodies handled using the immersed boundary method, non-Newtonian or multicomponent solvents, and electrohydrodynamics. However, there are also some disadvantages compared to other methods. First, because the method uses an explicit representation of the fluid domain, infinite domains cannot be considered in the present formulation. Second, because the method scales linearly in the number of fluid and particle degrees of freedom, phenomena involving only a few particles and large fluid domains are particularly inefficient to simulate using the RB-FIB method. Third, lubrication forces between nearly touching particles are not accurately handled for realistic number of blobs per particle, and the method becomes less efficient (in terms of both memory use and computing time) when there are many blobs per particle. Some of these shortcomings can be addressed by using explicit Green’s functions BrownianMultiblobSuspensions , or specializing to spherical particles Galerkin_Wall_Spheres ; AutophoreticSpheres_Adhikari ; BrownianDynamics_OrderNlogN .

One major advantage of the RB-FIB method not exploited in this work is it’s ability to include fluid body forces in the momentum equation. This allows for the inclusion of a range of multiphysics phenomena whose coupling with the fluid is via a fluid body force. An example are electrohydrodynamic phenomena for which the body force is the divergence of the Maxwell stress tensor Bhalla2014FullyRI . Since electrostatic fields are non–dissipative, no additional effort is needed to be taken to account for thermal fluctuations. However, charge separation creates a thin Debye layer near solid surfaces which is difficult to resolve numerically, but can be approximated using asymptotics ElectroneutralAsymptotics_Yariv . The ability of the RB-FIB method scheme to prescribe active slip velocities on the surfaces of the particles and the boundaries can be used to account for electroosmotic slip flows without resolving the Debye layers.

Acknowledgements.

This work was supported by the MRSEC Program of the National Science Foundation under Award Number DMR-1420073. This work was also partially supported by the National Science Foundation under collaborative award DMS-1418706 and by DMS-1418672. We thank Northwestern University’s Quest high performance computing service for the resources used to perform the simulations in this work. Brennan Sprinkle and Aleksandar Donev were supported in part by the Research Training Group in Modeling and Simulation funded by the National Science Foundation via grant RTG/DMS-1646339. Amneet Bhalla acknowledges research support provided by the San Diego State University.

Appendix A Temporal Accuracy for a Confined Colloidal Boomerang

To investigate the weak accuracy of the RB-FIB method scheme, we consider a single colloidal ‘boomerang’ suspended in a periodic slit channel, as depicted in figure 7. In addition to examining the temporal accuracy of the SEM scheme, this example demonstrates the ability of the method to handle arbitrary particle shapes as well as fully confined physical domains.

The fluid is water at room temperature $T=300$ K and viscosity $\eta=1$ mPa s. We model the boomerang using 15 equally spaced blobs with 7 of them forming each of the $2.1\mu$m arms of the ‘L’ shaped body. We take the blob radius $a$ to be related to the inter–particle spacing $s=0.3\mu$m as $s=0.925a$, as was recommended in RigidMultiblobs . As we employ the 6–point delta kernel New6ptKernel here, we have that $a=1.47h$ and hence the Eulerian mesh width $h=0.2208\mu$m.

The physical domain has extent $64h\times 64h\times 16h$ with no slip conditions on the top and bottom boundaries and periodic conditions in the $x,y$ directions. These dimensions were chosen to have some similar qualitative features to the quasi-2D confinement of colloidal boomerangs in the experiments reported in BoomerangDiffusion ; AsymmetricBoomerangs ; the boomerang barely has enough room to rotate fully in any direction if it is centered between the walls. We include a soft repulsive potential between the blobs and the wall of the form (44) but with $d=a$ and $r$ interpreted as the distance between the blobs and the wall. We use $b=0.1R_{H}$ and $\Phi_{0}=4k_{B}T$ as in MagneticRollers , as this choice ensures that the steric time scale associated with the potential isn’t much smaller than the diffusive time scale, while also maintaining a low probability of blob-wall overlaps. As in BrownianMultiblobSuspensions , we give each blob a buoyant mass of $m_{e}=1.57\times 10^{-11}$mg and thus apply a gravitation force $-m_{e}g\hat{z}$ on each of the 15 blobs, where $g$ is the Earth’s gravitational acceleration.

Figure 8(a) shows the probability distribution of the height of the boomerangs cross–point (the blob at the bend in the ‘L’ shape) computed using the SEM scheme with 3 different non–dimensionalized time step sizes, $\Delta\tau=\Delta t\,(k_{B}T)/(6\pi\eta a^{3})=\left\{0.14,0.072,0.036\right\}$. The largest time step size considered was chosen to be just under the empirically measured stability limit of $\Delta\tau\approx 0.22$. The black curve represents the correct marginal distribution which we obtain from the equilibrium Gibbs–Boltzmann distribution using an MCMC method. Comparing the marginal distributions generated by the SEM scheme to the correct distribution, we see clear improvement as the time step size is decreased. The inset of figure 8(b) shows that this decrease in the error is consistent with the first order weak accuracy proven in Appendix B. To underscore the importance of accounting for the stochastic drift term in (9), figure 8 also shows the biased equilibrium distribution obtained by using the SEM scheme without any RFD terms to account for the drift. We can see that there is a marked error in the biased distribution emphasizing the importance of the stochastic drift.

The cross–point height distribution shown in figure 8(a) is clearly bimodal. The peaks arise because of the shape of the particle and how it interacts with the confinement of the parallel walls. The asymmetry in the two peaks of distribution is partly due to the gravitational force. To interrogate this further, we track a unit vector $\hat{\boldsymbol{p}}$ (shown in figure 7) that is in the plane of the boomerang and intersects the boomerang’s cross–point at $45^{\circ}$. Figure 8(b) shows the distribution of the angle $\theta$ that $\hat{\boldsymbol{p}}$ makes with the $z$–axis, sampled using a Markov Chain Monte Carlo (MCMC) method as well as the SEM scheme. Further, we can see that while the mean of this distribution is fairly close to $\pi/2$ (for which $\hat{\boldsymbol{p}}$ is parallel with the walls) the variance is fairly wide. This shows that the cross-point frequently points up or down, which accounts for the bimodality of the cross–point height distribution. The inset in figure 8(b) demonstrates that the SEM scheme is first order weakly accurate in both position and orientation statistics as shown in Appendix B.

Appendix B Proof of Weak Accuracy for the SEM Scheme

As the SEM scheme is an adaptation of the Traction–corrected Euler–Maruyama scheme developed in BrownianMultiblobSuspensions , much of proof of first order accuracy from Appendix B2 in BrownianMultiblobSuspensions remains valid for the SEM scheme. We may identify $\boldsymbol{D}^{\boldsymbol{\mathcal{K}}^{T}}$ from section IV of this work with $\boldsymbol{D}^{F}$ of BrownianMultiblobSuspensions . With this identification, all that remains to show is that $\left\langle\boldsymbol{D}^{\boldsymbol{\mathcal{S}}}+\boldsymbol{D}^{\boldsymbol{\mathcal{J}}}-\boldsymbol{D}^{\boldsymbol{\mathcal{K}}}\right\rangle$ is equal to $\left\langle\boldsymbol{D}^{S}\right\rangle$ from BrownianMultiblobSuspensions and the proof of weak accuracy from Appendix B2 will suffice for the SEM scheme. In the interest of brevity, we will reuse the notation and results already established in Appendix B of BrownianMultiblobSuspensions .

From (38)–(39) and (40)–(41), we may write

[TABLE]

Expanding this using (33), we have in expectation

[TABLE]

Here we used the definition of $\boldsymbol{\mathcal{M}}$ from equation (24) in the second to last equality. The last equality simply identified the result with that obtained in equation (B14) of BrownianMultiblobSuspensions .