A simulation platform for slender, semiflexible, and inextensible fibers with Brownian hydrodynamics and steric repulsion

Abstract

The last few years have witnessed an explosion of new numerical methods for filament hydrodynamics. Aside from their ubiquity in biology, physics, and engineering, filaments present unique challenges from an applied-mathematical point of view. Their slenderness, inextensibility, semiflexibility, and meso-scale nature all require numerical methods that can handle multiple lengthscales in the presence of constraints. Accounting for Brownian motion while keeping the dynamics in detailed balance and on the constraint is difficult, as is including a background solvent, which couples the dynamics of multiple filaments together in a suspension. In this paper, we present a simulation platform for deterministic and Brownian inextensible filament dynamics which includes nonlocal fluid dynamics and steric repulsion. For nonlocal hydrodynamics, we define the mobility on a single filament using line integrals of Rotne-Prager-Yamakawa regularized singularities, and numerically preserve the symmetric positive definite property by using a thicker regularization width for the nonlocal integrals than for the self term. For steric repulsion, we introduce a soft local repulsive potential defined as a double-integral over two filaments, then present a scheme to identify and evaluate the nonzero components of the integrand. Using a temporal integrator developed in previous work, we demonstrate that Langevin dynamics sample from the equilibrium distribution of free filament shapes, and that the modeling error in using the thicker regularization is small. We conclude with two examples, sedimenting filaments and cross-linked fiber networks, in which nonlocal hydrodynamics does and does not generate long-range flow fields, respectively. In the latter case, we show that the effect of hydrodynamics can be accounted for through steric repulsion.

Figures (18)

• Figure 1: Spectral discretization. The discrete degrees of freedom $\boldsymbol{\tau}$ and $\boldsymbol{X}_\text{MP}$ define a set of Chebyshev points $\boldsymbol{X}$ and continuum interpolant $\mathbb{X}(s)$.
• Figure 2: Eigenvalues of possible corrections to the nonlocal mobility for $\hat{\epsilon}=10^{-3}$. Left panel: correcting the mobility via \ref{['eq:MFirst']}, by subtracting oversampled quadrature (with number of oversampled points $N_u$) and adding the special quadrature, gives a correction matrix with negative eigenvalues. Right panel: using a fatter radius ${\hat{a}}^*$ for the nonlocal mobility via \ref{['eq:MNew']}, then using the correction matrix $\widetilde{\boldsymbol{M}}_\text{SQS}-\widetilde{\boldsymbol{M}}_\text{SQS}^{(*)}$, always gives a symmetric positive definite correction. For these plots, $\widetilde{\boldsymbol{M}}_\text{SQS}$ is symmetrized, but the eigenvalues are not truncated (truncating them does not change the plots by eye). These plots are generated with a specific fiber ($N_x=25$ and $\hat{\epsilon}=10^{-3}$), but are unchanged when we vary the shape.
• Figure 3: Schematic of segment-based steric interaction algorithm. (a) We divide the fiber into pieces and use the midpoint of each piece (black circles) to identify possible close segments. (b) We then approximate each piece using a straight segment (dashed lines), find the minimum distance between the two segments, and use that as an initial guess for Newton's method. (c) Newton's method finds the local minimizer of the distance on the two fibers, around which we draw quadrature intervals (thick lines).
• Figure 4: End-to-end distance distribution with nonlocal hydrodynamics. In all cases we compare the results with Langevin dynamics with various time step sizes $\Delta t_f$ (in units of the slowest bending timescale) to MCMC, thereby determining the time step required to simulate the equilibrium distribution, and how it scales with $\ell_m$ (no change), $\hat{\epsilon}$ (no change), $\ell_p^*=\ell_p/L$ (relative time scale scales with $\ell_p^*$), and $N$ (time scale goes roughly as $N^{-4}$). The insets in each case show a probability histogram of the number of GMRES iterations required for convergence (tolerance $10^{-3}$).
• Figure 5: Deterministic dynamics of a single falling filament. The left panels show a series of snapshots of the sedimenting filament for a variety of elasto-gravitational numbers $\beta=gL^3/\kappa$. The right panel shows a summary of the vertical extension of the filament $h$ as a function of $\beta$, compared to a previous bead-spring model with slightly different model geometry cunha2022settling. Colored dots on the right plot correspond to the colored fiber snapshots shown in the left panels.