Complex dynamics of long, flexible fibers in shear

TL;DR

This paper investigates the complex behaviors of long, flexible fibers in shear flow, revealing transitions from rotation to buckling and snaking, using advanced simulations to capture multiple buckling sites and coiling phenomena.

Contribution

The study introduces a simulation framework based on regularized Stokeslets and fast multipole methods to model complex fiber dynamics, including multiple buckling and coiling, extending previous experimental and computational work.

Findings

Fibers exhibit buckling and snaking behaviors at high elasto-viscous numbers.

Simulations successfully reproduce multiple buckling sites and coiling.

Transitions depend on fiber slenderness and elasto-viscous number.

Abstract

The macroscopic properties of polymeric fluids are inherited from the material properties of the fibers embedded in the solvent. The behavior of such passive fibers in flow has been of interest in a wide range of systems, including cellular mechanics, nutrient aquisition by diatom chains in the ocean, and industrial applications such as paper manufacturing. The rotational dynamics and shape evolution of fibers in shear depends upon the slenderness of the fiber and the non-dimensional "elasto-viscous" number that measures the ratio of the fluid's viscous forces to the fiber's elastic forces. For a small elasto-viscous number, the nearly-rigid fiber rotates in the shear, but when the elasto-viscous number reaches a threshhold, buckling occurs. For even larger elasto-viscous numbers, there is a transition to a "snaking behavior" where the fiber remains aligned with the shear axis, but its ends curl in, in opposite directions. These experimentally-observed behaviors have recently been characterized computationally using slender-body theory and immersed boundary computations. However, classical experiments with nylon fibers and recent experiments with actin filaments have demonstrated that for even larger elasto-viscous numbers, multiple buckling sites and coiling can occur. Using a regularized Stokeslet framework coupled with a kernel independent fast multipole method, we present simulations that capture these complex fiber dynamics.

Tables

Table 1. Table 1: Geometric and material parameters of fibers and numerical parameters (all non-dimensional).

Quantities	Dimensionless
	value
Fiber length, $L$	$0.139-2.28$
Fiber radius, $R_f$	$0.005$
Slenderness ratio, $\rho =R_f/L$	$0.0022-0.036$
Spring stiffness, $k$	$0.0112$
Bending rigidity, $EI$	$4.8\times {10}^{-6}$
Shear scale, $\alpha$	$2.5$
Elasto-viscous number, $\overline{\mu }$	$5.35\times {10}^2-2.37\times {10}^7$
Numerical Parameters
Cross sections along fiber, $N_f$	$80-1313$
Points per fiber cross section, $N_c$	$18$
Total number of points on fiber $N$	$1440-23634$
Spacing between fiber nodes, $\Delta s$	$0.00174$
Blob size, $\epsilon$	0.00225
Time step, $\Delta t$	$1.0\times {10}^{-4}$

Table 2. Table 2: Nomenclature for KIFMM to compute Eq. ( 6 ).

$S_{\epsilon }$	regularized Stokeslet kernel
$S_0$	singular Stokeslet kernel
${𝐰}_j={𝐟}_j{A}_j$	the source strength at source point $j$
${𝐰}_k^{equiv}$	the equivalent source strength at equivalent point $k$
$p$	the number of equivalent points along each cubic box edge
$N_{leaf}$	the maximum number of points per leaf box

Equations

$$\begin{array}{rl}0& =-\nabla \hat{P}+\mu \Delta \hat{\mathbf{u}}+\hat{\mathbf{F}},\\ 0& =\nabla \cdot \hat{\mathbf{u}},\end{array}$$

$${\varphi }_{\epsilon }(\mathbf{x}-\mathbf{y})=\frac{15{\epsilon }^4}{8\pi (r^2+{\epsilon }^2{)}^{7/2}},$$

$${\mathbf{u}}_{st}(\mathbf{x})={\int }_{\Sigma }{S}_{\epsilon }(\mathbf{x},\mathbf{y})\hat{\mathbf{f}}(\mathbf{y})d{S}_y=\frac{1}{8\pi \mu }{\int }_{\Sigma }\frac{(r^2+2{\epsilon }^2)\hat{\mathbf{f}}(\mathbf{y})+(\hat{\mathbf{f}}(\mathbf{y})\cdot (\mathbf{x}-\mathbf{y}))(\mathbf{x}-\mathbf{y})}{(r^2+{\epsilon }^2{)}^{3/2}}d{S}_y,$$

$$\mathcal{E}=\frac{1}{2}\sum _j{k}_j{l}_j{\left(\frac{\Vert {\mathbf{x}}_{j_1}-{\mathbf{x}}_{j_2}\Vert }{l_j}-1\right)}^2$$

$${\text{mbox}\mathbf{f}}_{j_1}{A}_{j_1}=-\frac{\partial \mathcal{E}}{\partial {\mathbf{x}}_{j_1}}$$

$${\mathbf{u}}_{st}({\mathbf{x}}_i)=\sum _{j=1}^N{S}_{\epsilon }({\mathbf{x}}_i,{\mathbf{x}}_j){\text{mbox}\mathbf{f}}_j{A}_j.$$

$$EI$$

$$\bar{\mu }=\frac{8\pi \mu \dot{\gamma }{\hat{L}}^4}{c\hat{E}\hat{I}}=\frac{8\pi \alpha L^4}{cEI},$$

$$\mathbf{u}({\mathbf{x}}_i)$$

Full text

1 Introduction

The motion of flexible fibers in flow is central to many biological systems at the microscale. Mammalian sperm flagella propel these cells through the female reproductive tract [4], while microtubule fibers are ingredients of the mitotic spindle in cell division [21]. While the dynamics of these fiber-fluid systems are actuated by molecular motors, other biological systems contain passive fibers that are transported and undergo shape deformations due to the flow. Examples of these passive fibers include microtubules transported by cytoplasmic streaming in fungal hyphae [18] and chains of diatom cells that move with water currents through the ocean [9].

Early experiments by Forgacs and Mason [6] on synthetic fibers in shear demonstrated a spectrum of orbits and shape deformations. Shorter, stiffer fibers experienced a signature Jeffery orbit, where periodic tumbling was accompanied by little to no deformation. Longer fibers exhibited periodic orbits with shape deformations that were qualitatively catalogued as S-turns (buckling) and snake turns. For the longest fibers, Forgacs and Mason [6] observed complex shape deformations that they described as coiled orbits with entanglement. Decades later, with the availabilty of microfluidic technology, these periodic shape deformations in shear have been observed in DNA strands and actin filaments [8, 7, 13]. Most recently, the complex coiling and entanglement of long actin fibers was measured in [12], where an actin filament of length of more than sixty microns subjected to a shear flow. A fiber that is initially straight and aligned with the shear direction develops the shape of a hook at the ends of the fiber during a snake turn. Later in time, the fiber exhibits a more complex behavior, including multiple buckling sites in the middle of the fiber and three-dimensional entanglement.

In a review article on the dynamics of flexible fibers in flow, du Roure et. al. [3] describe recent technological advances in experimentation and recent algorithmic advances in computational modeling that have given rise to deeper understanding and probing of fiber-fluid systems. Exploiting the inertia-free environment at the microscale, and the slender geometry of the fibers, much progress has been made in using slender body theory [19] to describe the orbits and buckling of fibers in flow e.g. [22, 20, 15, 13]. In this manuscript, we present a mathematical model and numerical method that captures the complex shape deformations of the longest fibers, that is not easily achieved with a slender body formulation. The fibers we consider need not be thin and are not represented purely by their centerlines, but by a discretization of their surface. We use a similar fiber model as that used to examine the dynamics of diatom chains in a non-zero Reynolds number environment [16]. However, here we assume that the length and velocity scales are small enough so that the fluid dynamics is well-described by the Stokes equations, and use a regularized Stokeslet formulation [2, 1] of the fluid-fiber system.

While the diatom chain model of [16] was able to capture complex buckling behavior of long fibers, it was based on an adaptive mesh immersed boundary method that needed fine grid resolution on the fluid domain near the fiber. In contrast, the regularized Stokeslet formulation described here, while requiring fine resolution of the fiber surface to capture the complex behavior, is based on fundamental solutions of the Stokes equations and does not require a spatial grid on the surrounding fluid. Although one of the most attractive features of the regularized Stokeslet framework is the ease of implementation – the velocities at $N$ nodes are computed based upon the forces at each of the $N$ nodes - the direct $N^{2}$ evaluation becomes costly for $N$ large. In order to resolve the surface of a long, thin fiber using a discretization such that the distance between nodes around a cross section is on the order of the distance between cross sections, the number of nodes $N$ becomes necessarily large. Here we demonstrate that the incorporation of a kernel independent fast multipole method [26] to compute velocities rather than a direct evaluation allows for faster simulations of the longest fibers, and will be a promising tool for multi-fiber investigations.

In the following sections we will discuss the construct of the model fiber and its coupling to a Stokes fluid using the method of regularized Stokeslets. We will demonstrate the shape deformations of fibers of increasing length in shear, and discuss how these results compare with recent studies. Moreover, we will present simulations of long fibers (at large elasto-viscous numbers) that capture the complex coiling and entanglement observed in experiments.

2 Methodology

2.1 Stokes equations

Assuming that length and time scales are small, we model a flexible fiber coupled to a viscous fluid in three-space by the incompressible Stokes equations:

[TABLE]

where $\hat{P}$ is the pressure, $\hat{\mathbf{u}}$ is the fluid velocity, $\mu$ is the fluid viscosity, and $\hat{\mathbf{F}}$ is the external force per volume exerted by the fiber on the fluid. The forces in Eqn. 1 are localized at the fiber surface, and will be described below. The motion of the passive fiber will be driven by an imposed background shear flow ${\hat{\bf u}}_{b}(\hat{x},\hat{y},\hat{z})=(\dot{\gamma}\hat{y},0,0)$.

We use a regularized Stokeslet framework [1] to model the elastohydrodynamic system. Rather than using a surface integral of Dirac delta function forces $A(\mathbf{y})\hat{\mbox{\bf f}}(\mathbf{y})\delta(\mathbf{x}-\mathbf{y})$, we consider regularized force density $A(\mathbf{y})\hat{\mbox{\bf f}}(\mathbf{y})\phi_{\varepsilon}(\mathbf{x}-\mathbf{y})$ supported on a patch of area $A(\mathbf{y})$ on the surface of the fiber. The regularization (or blob) function is chosen to be:

[TABLE]

where $r=\|\mathbf{x}-\mathbf{y}\|$, and $\varepsilon$ is the regularizaiton parameter [1]. This leads to the velocities due to the regularized Stokeslets as follows:

[TABLE]

where ${\mathbf{\hat{f}}}(\mathbf{y})$ is force per unit area and $\Sigma$ denotes the surface of the fiber. As the regularization parameter $\varepsilon$ approaches zero, the kernel $S_{\varepsilon}$ approaches the classical singular Stokeslet $S_{0}$.

We nondimensionalize this coupled fluid-fiber problem by choosing a characteristic length scale $\mathcal{L}=\hat{L}_{0}$ on the order of a fiber length, a time scale $\mathcal{T}=\alpha\dot{\gamma}^{-1}$, a velocity scale $\mathcal{U}=\hat{L}_{0}/\mathcal{T}$, and a force scale ${{\mathcal{F}}}=\mu\hat{L}_{0}^{2}/\mathcal{T}$. Here $\alpha$ is a non-dimensional tuning parameter for background shear flow that is chosen to be $\alpha=2.5$ in all simulations shown below. We will use these non-dimensionalized quantities throughout the manuscript.

2.2 Representation of fiber

The model fiber that we consider has a native straight shape and equilibrium length $L$. We construct the discretization of the surface of the cylindrical fiber by placing cross-sections of radius $R_{f}$ along the centerline, perpendicular to the centerline (see Figure 1). Each cross-section is discretized by $N_{c}$ = 18 points, and we take $N_{f}$ cross-sections along the fiber so that the spacing between neighboring cross-sections is approximately equal to the spacing between adjacent points on a cross section.

Each of the $N=N_{c}\times N_{f}$ discrete points on the surface of the fiber is connected to a subset of the other surface points by a Hookean spring, giving elasticity to the structure. We define the elastic energy in the system as

[TABLE]

where $k_{j}$ is the stiffness of a spring with resting length $l_{j}$ that connects points $j_{1}$ and $j_{2}$. The sum is over all springs. The force at $\mathbf{x}_{j_{1}}$ is ${\mbox{\bf f}}_{j_{1}}A_{j_{1}}$ where $A_{j_{1}}$ is the area of a patch of surface centered at $\mathbf{x}_{j_{1}}$ in the discretization. The elastic forces are derived from the energy function:

[TABLE]

This network of nodes and elastic linkages will impart tensile stiffness and bending rigidity to the fiber, calibrated by the connectivity of the nodes and the stiffness constants of individual linkages. Similar constructs of semi-flexible filaments coupled to an incompressible fluid have been used to model bacterial flagella [10, 11, 5] and diatom chains [17]. In all simulations shown here, we choose a network connectivity so that each point on a given cross-section is connected to every other point on that cross-section, as well as to every other point on the two cross-sections adjacent to it. This means that each node is connected to $17+2\times 18=53$ other nodes. In addition, in all simulations shown, the stiffness constant $k_{j}=k$ in Eqn. 4 is taken to be the same for all springs. The resting lengths of the springs, $l_{j}$ in Eqn. 4 do vary with $j$, and are determined by the straight fiber configuration.

Due to the imposed background shear, the flexible fiber will depart from its equilibrium shape as the network springs become stretched or compressed, causing forces at the nodes to develop. The fluid velocity due to these elastic forces is evaluated at each material point $\mathbf{x}_{i},i=1,...,N$ of the fiber surface, using a discrete version of (Eq. (3)):

[TABLE]

Using the forward Euler method, this velocity, added to the background shear velocity $\mathbf{u}_{b}$, is used to update the positions of the nodes of the flexible fiber. The range of the fiber’s (non-dimensional) geometrical and elastic parameters used in simulations are shown in Table 1, along with the numerical parameters used.

2.3 Calculation of fiber bending rigidity and non-dimensionalization

The macroscopic bending rigidity $EI$ of the node-spring structure depends upon the individual spring constants $k_{j}$ and the topology of the spring network. Following [11, 17], we precompute the bending rigidity $EI$ for the model fiber by bending it into a circular arc with a prescribed radius of curvature $\kappa$. Because this curved shape stretches and compresses the network of springs, a non-zero elastic energy $\mathcal{E}$ from Eqn. 4 results, arriving at:

[TABLE]

As in other elastohydrodynamic systems where flexible fibers are coupled to a Stokesian incompressible fluid, the dynamics are governed by two non-dimensional parameters: the fiber aspect ratio (slenderness parameter) $\rho=R_{f}/L$, and the elasto-viscous number that measures the relative importance of flow forces to elastic forces (e.g. [22, 17, 25, 23, 13]). We define the elasto-viscous number:

[TABLE]

where we have indicated the ratio both in dimensional parameters and using our non-dimensional scaling. For a given slenderness ratio $\rho$, the tangential drag on the filament is given by the geomentric parameter $c=\ln({4e/{\rho}^{2}})$. Moreover, in this work, we do not consider the effect of thermal fluctuations.

3 Using mixed kernels in Kernel Independent FMM for regularized Stokeslet

3.1 Kernel Independent FMM steps

During the solution of the fiber dynamics, the velocity of each node on the fiber must be evaluated at every time step. The direct evaluation of the equation

[TABLE]

requires $O(N^{2})$ operations. The cost can be reduced to $O(N)$ operations with the Kernel Independent Fast Multipole Method (KIFMM) [26], which builds an adaptive octree for a given set of source (S) points $\mathbf{x}_{j}$ and target (T) points $\mathbf{x}_{k},1\leq k\leq N$, by recursively refining the octree with no more than $N_{leaf}$ points in each leaf box.

The interactions between points in a leaf box and other points are divided into near-field and far-field. Near-field interactions refer to points in all adjacent leaf boxes and are summed directly. Far-field interactions refer to points in non-adjacent leaf boxes and are approximated using a set of equivalent points with source strength $\mathbf{w}_{k}^{equiv}$ on each equivalent point $k$. In three dimensions, the equivalent points are chosen to be on a uniform cubic surface mesh surrounding each octree box, with $p$ points along each cubic edge. The total number of equivalent points is $6(p-1)^{2}+2$. This approximation converges exponentially with increasing number $p$. The contribution from near- and far-fields are added together to form $\mathbf{u}$ at each target point.

The tradeoff between cost and accuracy is controlled by $p$. In general, $p=10$ gives single-precision accuracy and $p=16$ gives double precision accuracy. $N_{leaf}$ controls the depth of the octree, and thus affects the total computation time. In practice, $N_{leaf}$ is set to about 2000 to fit the current CPU architecture.

There are two ways to invoke the KIFMM for $S_{\varepsilon}$. The first choice is to directly use $S_{\varepsilon}\in\mathcal{R}^{3\times 3}$ as the kernel for all operations throughout the octree, for both near-field and far-field interactions. This approach is straightforward to implement, but is limited to a common value of $\varepsilon$ for all points. Further, if periodic boundary condition are necessary, the periodizing operator $\mathcal{T}_{M2L}$ must be recalculated for each different $\varepsilon$ [24]. Another approach is to regard each source point in $S_{\varepsilon}$ as a four dimensional vector $(w_{1},w_{2},w_{3},\varepsilon)$. In this case, $S_{\varepsilon}$ becomes a nonlinear kernel because $\varepsilon$ appears nonlinearly in $S_{\varepsilon}$. The near-field interactions between source-target pairs are computed directly. For far-field interactions, the singular kernel $S_{0}$ is used for equivalent points, and the strength $\mathbf{w}_{k}^{equiv}$ is found by matching the equivalent flow field with the regularized flow field generated by source points. It is straightforward to find the equivalent $\mathbf{w}_{k}^{equiv}$ and the exponential convergence of KIFMM is maintained because the regularized kernel $S_{\varepsilon}$ is a solution to the Stokes equation. After $\mathbf{w}_{k}^{equiv}$ is found, the singular kernel $S_{0}$ is used throughout the tree traversal. This approach allows the regularization parameter $\varepsilon$ to vary for different source points, and allows the direct reuse of $\mathcal{T}_{M2L}$ computed for $S_{0}$ [24] to implement various types of periodic boundary conditions.

The computation code is implemented based on the parallel KIFMM library PVFMM [14]. Both MPI and OpenMP parallelism are implemented to improve parallelization efficiency. SIMD instructions are also utilized to improve efficiency on modern CPU architectures.

4 Results and Discussion

In recent laboratory experiments that investigate the dynamics of actin filaments in shear [13], the diameter and bending rigidity of the fibers are fixed by nature. The elasto-viscous number of the experiments is varied by either adjusting the background shear or observing actin filaments of different lengths. Motivated by these experiments, here we choose to examine the dynamics of fibers in shear by keeping their bending rigidity and their diameter fixed, but vary their length. The shear rate remains fixed in these simulations. As the fiber length increases, its slenderness ratio decreases, and the corresponding elasto-viscous number of the system increases.

Although thermal fluctuations play a role in the dynamics of fiber-flow systems, we choose not to include any Brownian forces in the simulations presented here. Moreover, because we track the surface of the fiber rather than just its centerline as in slender body models, a straight fiber that is initialized with its centerline along the $x$-axis needs no perturbation to begin its orbit – spatial gradients in velocity on each cross section are immediately formed. In the absence of Brownian fluctuations, the dynamics of a perfectly cylindrical fiber whose centerline initially coincided with the $x$-axis in this unbounded shear flow must necessarily obey some symmetry constraints. In particular, the centerline must remain in the plane that it was initialized in and the shape deformations with respect to the centroid of the fiber must be odd. In the simulations presented below, we will see the effects of the small fluctuations that occur due to the numerical perturbations that arise from time integration and from the finite discretization of the fiber surface.

Figure 2 shows the periodic orbits of three fibers of increasing length that display the signature tumble, $S$-turn and snaking behavior reported by [6, 13, 22, 16]. The first column of Figure 2 depicts the rotational orbit of a fiber of length $L=0.139\ (\bar{\mu}=5.35\times 10^{2}$) that exhibits little deformation from its straight shape. Figure 3 (a) shows a surface plot of the evolving (absolute value of) curvature along the arclength of the fiber as a function of time during approximately five orbits. We see that two slight bends occur during each rotation when the fiber is aligned with the maximal compressive region of the shear (a forty-five degree angle with the negative $x$-axis), as in the $T=4.95$ snapshot in the first column of Figure 2. The positions of maximal curvature do not travel along the arclength of the fiber, but occur as standing waves that appear periodically in each orbit. The second column of Figure 2 depicts the rotational orbit of a fiber of length $L=0.196\ (\bar{\mu}=1.97\times 10^{3}$) that clearly buckles into an $S$-shape during its orbit. Figure 3 (b) shows the corresponding evolution of curvature. Again we see standing waves of curvature. In contrast, the third column of Figure 2 shows the snaking behavior of a fiber of length $L=0.523\ (\bar{\mu}=8.22\times 10^{4})$. Here the ends of the fiber curl in towards the middle of the fiber in an antisymmetric manner, and the points of maximal curvature travel along its arclength. These traveling waves of maximal curvature are evident in Figure 3 (c) during the times at which the fiber is bent from its straight configuration. In each of the simulations shown in Figure 2, the fibers exhibit periodic orbits, and during most of the orbit the fibers remain in a nearly straight configuration aligned with the $x$-axis. Note that the period of rotation increases with the length of the fiber.

We compare the above with the coordinated experiments and simulations of actin fibers with Brownian fluctuations by Liu et. al. [13] where the values of elasto-viscous number $\bar{\mu}$ at which transitions from tumbling to buckling to snaking occurred were quantified. The elasto-viscous numbers of the simulations presented in Figure 2 fall squarely within the ranges of these different dynamical regimes. We note that because of the imposed symmetry in our model, in the absence of a perturbed initial position or Brownian fluctuations, we do not observe the $C$-buckling or $U$-turns reported in [13], but rather their counterparts with odd symmetry. If perturbations were added to the initial placement of the fiber, we do observe such asymmetric shapes (not shown).

We now examine the dynamics of longer fibers with values of elasto-viscous number $\bar{\mu}$ that transition from the snaking behavior shown in the third column of Figure 2 to more complex dynamics and shape evolution. Figure 4 shows some time snapshots of a fiber of length $L=0.931$ in shear with a corresponding elasto-viscous number ($\bar{\mu}=7.51\times 10^{5}$) that is beyond those considered in [13]. At time $T=4.95$ we observe the emergence of the hooks at each end that are evident in snaking behavior. However, it is worthwhile to compare the snapshot of the $L=0.523$ fiber in Figure 2 at $T=10.5$ to the longer fiber in Figure 4 at $T=16.25$. In the longer fiber, we see that additional hooks have formed at the ends, giving multiple local extrema in curvature along the fiber. The evolution of the curvature is shown in Figure 3 (d), where we can see the propagation of traveling waves of curvature. We remark that this fiber never did regain its straight shape. Figure 5 shows the dynamics of an even longer fiber of length $L=2.28$ that exhibits even richer shape dynamics. Again, the hooks appear at the fiber ends at $T=10$, but the fiber is long enough to support multiple coils as time evolves, as well as additional buckling sites in the middle. Figure 3 (e) shows the evolution of the curvature along this fiber with a cross-section at time $T=36$ indicated. Figure 3 (f) shows the curvature plotted at $T=36$ as a function of arclength. The multiple peaks in curvature can be identified in the fiber configuration also shown. In Figure 5 we can observe that at about $T=50$, the perturbations due to numerical fluctuations cause the centerline of the fiber to move out of the plane, giving rise to the entanglement seen by Forgacs and Mason [6] and in [12]. Figure 6 shows a zoomed-in image of this entanglement at time $T=110$.

Both Harasim et. al. [7] and Liu et. al [13] presented theoretical analysis of the dynamics of the $J$ shape at the fiber ends that evolve during a snaking orbit. In particular, it was assumed that the hook developed is well-approximated by a semicircle of a fixed radius. For fibers in our simulations that did cross the threshhold of $\bar{\mu}$ at which snaking occurs, we measured the radius of such a circle using a least squares fit (see Figure 7 (a) ). Even for large values of $\bar{\mu}$, the initiation of bending exhibits the emergence of hooks. Figure 7 (b) shows that the emergent radii is nearly independent of the length of the fiber and, hence, independent of $\bar{\mu}$. Because our computations are nearly at the limit of long filaments, this agrees with the predictions of [7]. During the initiation of the snaking in each of the simulations (even those that go on to complex shape deformations with no periodic orbits), we can measure the speed of propagation of the maximal curvature along the fiber arclengths. This is equivalent to the slope of the traveling waves in the curvature surface plots in Figure 3 (c-e). Although not immediately apparent, because the range of the spatial axis is the arclength of the fiber which is different in each panel, these slopes are nearly equal. As with the hook radii, we find that these “snaking velocities” are also nearly independent of $\bar{\mu}$ (see Figure 7 (b)).

4.1 Computational considerations

To capture the complex coiling and entanglement dynamics of slender fibers, fine resolution of their surface and along their length is required. While the regularized Stokeslet formulation relies on fundamental solutions of the Stokes equations and not a finite difference or finite element discretation of the surrounding three-dimensional fluid domain, a direct $N^{2}$ evaluation of $N$ velocities at $N$ nodes becomes prohibitive for large $N$. Figure 8 presents the timings for 1000 time steps of the fiber-fluid system using direct summation and KIFMM summation. We see that for small fiber lengths, the overhead for KIFMM outweighs its benefits, but that a cross-over occurs at about $N=1.7\times 10^{4}$ nodes, corresponding to a fiber length of $L=1.6$. For the longest fiber of length $L=2.2$ we already see a factor of 1.6 speed-up. This speed-up will be significant for simulations of multiple fiber dynamics.

Note that in the simulations presented above, and in the entangled state shown in Figure 6, there is no self-intersection of the fiber. Because of finite time steps and finite regularization parameters, we cannot rule out self-intersection. While we have not implemented any repulsive forces as fiber points get very close, such adjustments to the algorithm can be included.

5 Conclusions

In summary, we have presented a computational method that captures complex shape deformations of long flexible fibers in shear. Using the method of regularized Stokeslets in conjunction with a kernel independent fast multipole method we are able to simulate the dynamics of fibers at elasto-viscous numbers that are a couple of orders of magnitudes larger than those reported using slender body formulations. Because the fiber is represented by a spring-node system with individual spring elements, it is straightforward to model fibers with inhomogenous material properties by altering connectivities and stiffness constants of the individual elements. Moreover, the implementation of KIFMM will allow us to probe the interaction of multiple fibers in an array of background flows beyond a simple linear shear.

6 Acknowledgements

This research was supported by the National Science Foundation grant DMS-1043626 and the Gulf of Mexico Research Initiative. The authors would like to thank Olivia du Roure, Yanan Liu, Anke Lindner and Michael Shelley for helpful discussions.