A fluctuating boundary integral method for Brownian suspensions Yuanxun Bill Bao 1 Manas Rachh 2 Eric Keaveny 3 Leslie Greengard 1 Aleksandar Donev 1 1 Courant Institute, NYU 2 Yale University 3 Imperial College London Complex Creeping Fluids: Numerical Methods and Theory BIRS-CMO, October 5, 2017 Y. Bao (Courant) Fluctuating BIE
Brownian Particles of Complex Shape Figure: Brownian motion of passive boomerang colloidal particles from Chakrabarty et al. 2013 Y. Bao (Courant) Fluctuating BIE
Brownian Particles of Complex Shape Figure: Brownian suspension of 40 starfish-shaped particles. Y. Bao (Courant) Fluctuating BIE
Brownian Dynamics with Hydrodynamic Interactions (BD-HI) ⊲ Consider a suspension of N rigid bodies with configuration Q = { q β , θ β } N β =1 consisting of positions and orientations immersed in a Stokes fluid. ⊲ The Ito stochastic equation of Brownian Dynamics (BD) is d Q 1 2 W ( t ) + ( k B T ) ∂ Q · N , dt = N F + (2 k B T N ) where N ( Q ) is the body mobility matrix , F = { f β , τ β } N β =1 is the applied forces and torques , k B T is the temperature, and W ( t ) is a vector of independent white noise processes. ⊲ The stochastic noise amplitude satisfies the � � ∗ 1 1 fluctuation-dissipation balance : N N = N . 2 2 ⊲ The stochastic drift term ∂ Q · N = � j ∂ j N ij is related to the Ito interpretation of the noise. Y. Bao (Courant) Fluctuating BIE
Hydrodynamic Body Mobility Matrix ⊲ The body mobility matrix N ( Q ) ≻ 0 is a symmetric positive definite (SPD) and it includes hydrodynamic interactions and (periodic) boundary conditions . ⊲ For viscous-dominated flows ( Re → 0), we can assume steady Stokes flow and solve the deterministic Stokes mobility problem , U = N F , where U = { u β , ω β } N β =1 collects the linear and angular velocities. ⊲ At every time step of BD simulation, we need to generate particle velocity in the form of (dropping k B T ), 1 � 2 W . U = N F + N ⊲ This talk: How to accurately and efficiently compute the action of 1 2 ? N and N Y. Bao (Courant) Fluctuating BIE
Motivations The majority of BD-HI methods in chemical engineering: ⊲ are limited to spherical particles; ⇒ low ⊲ use an uncontrolled truncation of multipole expansion = accuracy in the near field; 1 2 W . ⊲ scale super-linearly for generating N We combine ideas from Positively Split Ewald (Fiore et al. 2017) and boundary integral methods to develop methods that ⊲ handle particles of nonspherical, complex shape ; ⊲ have controlled accuracy for dense suspensions; ⊲ achieve linear-scaling with the number of particles. = ⇒ 1st boundary integral method that accounts for Brownian motion of nonspherical particles immersed in a viscous incompressible fluid. Bao et al. 2017, submitted to JCP [arXiv:1709.01480]. Y. Bao (Courant) Fluctuating BIE
First-Kind Boundary Integral Formulation 1 2 W and solve a mobility ⊲ Let us first ignore Brownian terms N problem to compute N F . ⊲ For simplicity, consider only a single body Ω. The first-kind boundary integral equation for the mobility problem, � u + ω × ( x − q ) = v ( x ∈ ∂ Ω) = G ( x − y ) ψ ( y ) dS y , (1) ∂ Ω along with force and torque balance conditions, � � ψ ( x ) dS x = f and ( x − q ) × ψ ( x ) dS x = τ , (2) ∂ Ω ∂ Ω where ψ ( x ∈ ∂ Ω) is the traction and G is the (periodic) Green’s function for the Stokes flow (Stokeslet). ⊲ Note that one can alternatively use a completed second-kind or a mixed first-second kind formulation for improved conditioning. We only know how to generate Brownian displacements efficiently in the first-kind formulation. Y. Bao (Courant) Fluctuating BIE
First-Kind Boundary Integral Formulation ⊲ Assume that the surface of the body ∂ Ω is discretized in some manner, and the single-layer operator M is approximated by some quadrature, � G ( x − y ) ψ ( y ) dS y ≡ M ψ → M µ , ∂ Ω where M is a SPD operator with kernel G with r − 1 singularity in 3D (log r in 2D), discretized as a SPD single-layer matrix M . ⊲ In matrix notation the mobility problem can be written as a saddle-point linear system for the surface forces µ and rigid-body motion U = { u , ω } , � � � µ � � 0 � − K M = − , (3) − K ⊤ 0 U F where K ⊤ µ is the discrete force and torque balances. ⊲ Using Schur complement to eliminate µ , we get U = N F = ( K ⊤ M − 1 K ) − 1 F . Y. Bao (Courant) Fluctuating BIE
Brownian Displacements 1 2 ? ⊲ How do we compute the action of N More precisely, how to generate a Gaussian random vector with covariance N ? ⊲ Assume for now we knew how to generate a random surface velocity ˘ v ( x ∈ ∂ Ω) with covariance M (periodic Stokeslet). ⊲ In the continuum setting, ˘ v is a distribution , not a function. Formally, � ˘ v � = M . v ˘ Y. Bao (Courant) Fluctuating BIE
Brownian Displacements 1 2 W , ⊲ Key idea 1 : solve the mobility problem with ˘ v = M � � � � � ˘ � M − K v µ = − , (4) − K ⊤ 0 U F 1 1 ⇒ U = N F + N K ⊤ M − 1 M 2 W , 2 W = N F + N = 1 2 with the correct covariance: which defines a N � � ∗ � � ∗ 1 1 1 1 = N K ⊤ M − 1 M M − 1 K N N N M 2 2 2 2 = N ( K ⊤ M − 1 K ) N = N ( N ) − 1 N = N . (5) ⊲ ˘ v comes from the stochastic stress tensor in fluctuating hydrodynamics after eliminating the fluid velocity and pressure from the fluid+body equations of motion. Y. Bao (Courant) Fluctuating BIE
Ewald Splitting of M 1 2 W . ⊲ We need accurate and fast algorithms for M µ and M ⊲ Key idea 2 : split the Stokeslet into SPD kernels: ) ≡ G ( r ) + G ( w ) G = H ∗ G + ( G − H ∗ G , (6) � �� � � �� � ξ ξ singular smooth far-field near-field � � 1 + k 2 e − k 2 / 4 ξ 2 . where the Hasimoto function is � H ( k ; ξ ) = 4 ξ 2 This is the key idea of Positively Split Ewald (Fiore et al., J. Chem. Phys., 2017). ⊲ Recall that � G ( x − y ) ψ ( y ) dS y ≡ M ψ ≈ M µ . ∂ Ω The splitting of G naturally induces the splitting of M , and subsequently the splitting of M into SPD matrices , M = M ( r ) + M ( w ) . (7) Y. Bao (Courant) Fluctuating BIE
Fast Algorithm for M µ ⊲ M ( r ) µ : To handle the singularity of G ( r ) ξ , we need to employ singular quadrature (e.g., Alpert quadrature in 2D). Since G ( r ) decays as exp( − ξ 2 r 2 ), the mat-vec product M ( r ) µ ξ can be computed rapidly with linear-scaling in the real space. ⊲ M ( w ) µ : The far-field sum can be accelerated by the Spectral Ewald method (Lindbo/Tornberg) in Fouier space, M ( w ) = D † BD , (8) where D is the non-uniform FFT (Greengard/Lee), and B is a SPD block-diagonal matrix (in Fourier space), � H ( k , ξ ) ( I − ˆ k ˆ B ( k , ξ ) = k ) . k 2 Y. Bao (Courant) Fluctuating BIE
Splitting of M 1 / 2 W ⊲ Key Idea 3 : generate random surface velocity with covariance M by � M ( w ) � 1 � M ( r ) � 1 2 W ( w ) + 2 W ( r ) , 1 d . 2 W M = (9) if both M ( w ) and M ( r ) are SPD and � W ( w ) W ( r ) � = 0 . ⊲ The far-field matrix M ( w ) is SPD by construction, and we can write � � � � † M ( w ) = D † BD = 1 1 D † B D † B , (10) 2 2 so that the wave-space random surface velocity can be generated with a single call to the NUFFT, � M ( w ) � 1 2 W ( w ) = D † B 1 2 W ( w ) . (11) Y. Bao (Courant) Fluctuating BIE
Near-Field Contribution of M 1 / 2 W � M ( r ) � 1 2 W ( r ) can be ⊲ The near-field random surface velocity generated by a Krylov Lanczos method of Chow/Saad. ⊲ Because of the short-ranged nature of M ( r ) , the Lanczos method converges in a reasonable amount of iterations. ⊲ However, in general, M ( r ) is not symmetric, so M ( r ) is not SPD strictly speaking, because the Alpert quadrature does not enforce the SPD of M . � M ( r ) � ⊤ � � M ( r ) + ⊲ Nevertheless, we find that symmetrizing 1 2 preserves the order of accuracy of Alpert quadrature, and the Krylov Lanczos iteration is rather insensitive to any small negative eigenvalues. Y. Bao (Courant) Fluctuating BIE
Block-Diagonal Preconditioners � � � µ � � � 1 − K 2 W M M = − . − K ⊤ 0 U F ⊲ To mitigate the inherent ill-conditioning of M due to the use of a first-kind boundary integral formulation, we apply a block-diagonal preconditioner , i.e. , we simply neglect all hydrodynamic interactions between distinct bodies in the preconditioner, both when solving the saddle-point mobility problem using GMRES, and � M ( r ) � 1 2 W ( r ) . in the Lanczos iteration for generating ⊲ Both preconditioners can be precomputed using LAPACK for a single body, and then applied to many bodies via two fast vector rotations per body. ⊲ GMRES and Lanczos converge in a constant number of iterations, growing only weakly with packing fraction. Y. Bao (Courant) Fluctuating BIE
Recommend
More recommend
