Accurate close interactions of Stokes spheres using lubrication-adapted image systems

Abstract

Stokes flows with near-touching rigid particles induce near-singular lubrication forces under relative motion, making their accurate numerical treatment challenging. With the aim of controlling the accuracy with a computationally cheap method, we present a new technique that combines the method of fundamental solutions (MFS) with the method of images. For rigid spheres, we propose to represent the flow using Stokeslet proxy sources on interior spheres, augmented by lines of image sources adapted to each near-contact to resolve lubrication. Source strengths are found by a least-squares solve at contact-adapted boundary collocation nodes. We include extensive numerical tests, and validate against reference solutions from a well-resolved boundary integral formulation. With less than 60 additional image sources per particle per contact, we show controlled uniform accuracy to three relative digits in surface velocities, and up to five digits in particle forces and torques, for all separations down to a thousandth of the radius. In the special case of flows around fixed particles, the proxy sphere alone gives controlled accuracy. A one-body preconditioning strategy allows acceleration with the fast multipole method, hence close to linear scaling in the number of particles. This is demonstrated by solving problems of up to 2000 spheres on a workstation using only 700 proxy sources per particle.

Figures (18)

• Figure 1: By using an adaptive number of additional source points per closely interacting particle pair, lubrication forces between 20 closely interacting spheres in a cluster (many separated by $10^{-3}$ radii) are resolved with our new technique. Each sphere travels with a random translational and angular velocity. The target accuracy of $10^{-3}$ in the relative residual is reached everywhere on the particle surfaces.
• Figure 2: Panel (a) displays the basic discretization, with the original source points $\boldsymbol{y}_j$, $j = 1,\dots,N$, on the radius $R_p$ proxy-surface in blue, and collocation points $\boldsymbol{x}_i$, $i = 1,\dots,M$, in black on the unit sphere surface. Both source and collocation points are obtained from spherical design nodes. Panels (b) and (c) display the image-enhanced grid of source and collocation points for a pair of close-to-touching spheres. Additional image points are distributed along a line segment in the interior of each sphere, while additional collocation points are distributed "above" each line of images. The collocation points are located on two spherical caps determined by angles $\beta_1$ and $\beta_2$ surrounding the point of closest approach to another particle, as in panel (b). If a particle is closely interacting with more than one neighbor, one region of extra source and collocation points of this type is added for every near-contact. In panel (c), collocation points are indicated in black and image points in red, generated via \ref{['line_points']}-\ref{['tj']}.
• Figure 3: Lines of image points obtained by repeated reflections of each sphere center in the other sphere via \ref{['reflect']}. Accumulation points for the lines of images depend on the separation distance $\delta$ between the two spheres.
• Figure 4: Convergence test for two unit spheres $\delta$ apart, discretized with $N$ source points at inner proxy-surfaces of radius $R_p = 0.6$, and affected by two types of boundary conditions, translation or rotation, as described in main text. Panel (a) displays rate dependence on $\delta$. In panel (b), dashed curves show convergence as predicted by $\mathcal{O}(R_{\text{acc}}^{\sqrt{N}})$ for the max relative residual $\epsilon_{\text{res}}^{\max}$, while the max force/torque error $\epsilon_{FT}$, relative to a reference solution determined with $N=3529$, converges at a rate of $\mathcal{O}(R_{\text{acc}}^{1.5\sqrt{N}})$.
• Figure 5: A pair of unit spheres discretized with $N$ sources on proxy-surfaces of radius $R_p$ are separated by a small distance $\delta$. A sweep over $R_p$ and $\delta$ is made to investigate accuracy and MFS coefficient magnitudes. The test is done for two settings: case 1 (squeezing motion), and case 2 (shearing motion). Squeezing motion is harder to resolve, confirmed by comparing panels (a) and (b). The error in the computed force relative to the results of Brenner in \ref{['Brenner']}, denoted by $\epsilon_{\text{FT}}$, is a more forgiving measure of the error than the maximum relative residual $\epsilon_{\text{res}}^{\max}$; compare panels (c) and (b) for $N=686$ and (f) and (e) for $N=1353$. For a given $R_p$, additional image points can only be added for $\delta<\delta^*$ (to the left of the red curve representing $\delta^*(R_p)$ in each panel). For large $\delta$, $\epsilon_{\text{FT}}$ and $\epsilon_{\text{res}}^{\max}$ are capped at $10^{-6}$ which is the tolerance chosen for GMRES. Vertical black lines are drawn for reference at $\delta = 5\cdot 10^{-2}$, $\delta = 1\cdot 10^{-1}$, $\delta = 5\cdot 10^{-1}$ and $\delta = 1$.

Theorems & Definitions (13)

• Remark 1: Recent Stokes MFS work
• Remark 2: Conditioning
• Remark 3: Stresslet strengths
• Remark 4: Extracting net forces and torques
• example 1
• example 2
• Remark 5: Left-preconditioning
• example 3: A self-convergence test for the exterior flow field
• example 4: Accuracy gain with adaptive images
• example 5: Comparison of accuracy and efficiency relative to a BIE scheme