Dynamics of flexible fibers in viscous flows and fluids

TL;DR

This paper reviews recent advances in understanding the behavior of flexible fibers in viscous flows at low Reynolds numbers, highlighting experimental and modeling developments relevant to industrial and biological applications.

Contribution

It provides a comprehensive overview of experimental techniques and mathematical models that have advanced the study of flexible fiber dynamics in viscous fluids.

Findings

New fiber synthesis methods enable precise experimental studies.

Enhanced microfluidic and imaging techniques improve measurement accuracy.

Progress in mathematical modeling captures complex fiber-fluid interactions.

Abstract

The dynamics and deformations of immersed flexible fibers are at the heart of important industrial and biological processes, induce peculiar mechanical and transport properties in the fluids that contain them, and are the basis for novel methods of flow control. Here we focus on the low Reynolds number regime where advances in studying these fiber-fluid systems have been especially rapid. On the experimental side this is due to new methods of fiber synthesis, microfluidic flow control, and of microscope based tracking measurement techniques. Likewise, there have been continuous improvements in the specialized mathematical modeling and numerical methods needed to capture the interactions of slender flexible fibers with flows, boundaries, and each other.

Equations

$$\mu \Delta \mathbf{u}-\nabla p=0\text{mbox}\&\nabla \cdot \mathbf{u}=0,$$

$$\mathbf{u}(\mathbf{x})=\bar{\mathbf{u}}(\mathbf{x})+\sum _{m=1}^N{\mathbf{u}}_F^m(\mathbf{x})\text{mbox}where$$

$${\mathbf{V}}^n(s)=\bar{\mathbf{u}}({\mathbf{X}}^n)+\sum _{m=1,m\ne n}^N{\mathbf{u}}_F^m({\mathbf{X}}^n)+c_m\left(\mathbf{I}+{\mathbf{X}}_s^n{\mathbf{X}}_s^n\right){\mathbf{f}}^m+{\mathbf{K}}_n[{\mathbf{f}}^n]$$

$${\mathbf{f}}^n={\mathbf{f}}_{\text{E}}^n=-E{\mathbf{X}}_{ssss}^n+(T^n{\mathbf{X}}_s^n{)}_s,0\le s\le L_n.$$

$$\mathbf{R}(s,t)\cdot \left[\mathbf{V}(s,t)-\bar{\mathbf{u}}(\mathbf{X}(s,t),t)\right]=\mathbf{f}(s,t),$$

$${\mathbf{X}}_t=\mathbf{V}\equiv \mathbf{u}(\mathbf{X},t)+{\eta }^{-1}(\mathbf{I}+{\mathbf{X}}_s{\mathbf{X}}_s)\left[-E{\mathbf{X}}_{ssss}+(T{\mathbf{X}}_s{)}_s\right]\&{\mathbf{X}}_s\cdot {\mathbf{V}}_s=0$$

$$2T_{ss}-|{\mathbf{X}}_{ss}{|}^{2}T=S,$$

$$\mathbf{V}={\alpha }^{-1}\mathbf{U}(\alpha \mathbf{X},t)-{\tilde{\eta }}^{-1}(\mathbf{I}+{\mathbf{X}}_s{\mathbf{X}}_s)\left[{\mathbf{X}}_{ssss}-(T{\mathbf{X}}_s{)}_s\right]$$

$$\mathbf{R}\cdot \left({\mathbf{X}}_t-\bar{\mathbf{u}}\right)=-E{\mathbf{X}}_{ssss}+{\left(T{\mathbf{X}}_s\right)}_s+{\mathbf{f}}^B(s)$$

$$⟨{\mathbf{f}}^B(s,t)⟩=\text{b}0\&⟨{\mathbf{f}}^B(s,t){\mathbf{f}}^B(s^{\prime },t^{\prime })⟩=2kT\mathbf{R}(s,t)\delta (s-s^{\prime })\delta (t-t^{\prime })$$

$${\dot{\mathbf{X}}}_c=\mathbf{A}{\mathbf{X}}_c,\dot{\mathbf{e}}=\left(\mathbf{I}-\mathbf{e}{\mathbf{e}}^T\right)\mathbf{A}\mathbf{e},\bar{T}=-\frac{\tilde{\eta }}{4}({\mathbf{e}}^T\mathbf{E}\mathbf{e})(s^2-1/4)$$

$$\sigma w=w+s{w}_s+\frac{1}{4}\left(s^2-1/4\right)w_{ss}-{\tilde{\eta }}^{-1}{w}_{ssss}$$

$$\tilde{\eta }\left(w_t-w\cos \tilde{\omega }t\right)=-w_{ssss}+2{\bar{T}}_s{w}_s+\bar{T}{w}_{ss},\&\bar{T}=-\frac{\tilde{\eta }}{4}(s^2-\frac{1}{4})\cos \tilde{\omega }t,$$

Full text

I INTRODUCTION AND BACKGROUND

An important class of micro-scale fluid-structure interactions involves the interactions and deformations of flexible fibers with fluid flows. This is evident in the many biological transport processes, such as microorganismal swimming (Lauga & Powers (2009)) or nuclear positioning in eukaryotic cells (Shelley (2016)), that involve flexible fibers either actuated or passive. Fibers are the microstructure of many complex fluids – biological, industrial, and synthetic – studied both for their scientific and industrial importance and for their peculiar mechanical responses (Broedersz & MacKintosh (2014)). Such suspensions are particularly challenging to study because the suspended fibers have many degrees of freedom in deformation, and can exhibit microscopic instabilities. This makes the interaction of fibers with background flows surprisingly complex. Flexible fibers, both anchored and freely suspended, have been studied of late for the purposes in microfluidic flow control, and when actuated can exhibit complex collective and transport dynamics.

This review focuses on the low Reynolds number regime, as progress on fiber-fluid interactions in this regime has been especially rapid. On the experimental side, this is due to improvements in fiber synthesis and characterization, microfluidic flow control and improved microscope-based measurement and particle tracking techniques. By taking advantage of the relative simplicity of the Stokes equations (as opposed to Navier-Stokes), combined with adaptive resolution and fast summation approaches, there has been a likewise rapid improvement in numerical methods for simulating the deformations and interactions of fibers with flows, as well as with each other and other immersed structures. Coarse-grained descriptions of fiber assemblies and suspensions remain in their early stages, but their development is being sped by these new numerical methodologies.

In outline, we first discuss the current state-of-art in experimental synthesis and measurement techniques, and in numerical methods for dynamical simulation fiber-fluid interactions. We then describe the results of both experiment and theory for the dynamics of free and anchored fibers, followed by a review of the current state of research in many-fiber/fluid systems. We close with a discussion of future directions.

I.0.1 Setting the stage

Here we briefly introduce a few important physical parameters and relations. Consider a slender elastic fiber of length $L$, of circular cross-section with radius $a$ (hence $\epsilon=a/L\ll 1$), and flexural rigidity $E=YI$ with $Y$ the material Youngs modulus and $I$ the areal moment of inertia ($I=\pi a^{4}/4$). This fiber is immersed in a Newtonian fluid of shear viscosity $\mu$ with the fluid motion characterized by a strain-rate ${\dot{\gamma}}$. Neglecting inertial forces, three important forces are at play: Brownian forces $\sim k_{B}T/L$ ($k_{B}$ is the Boltzmann constant and $T$ the temperature), drag forces $\sim\mu{\dot{\gamma}}L^{2}$, and elasticity forces $\sim E/L^{2}$. For most of the work reviewed here, though not all, viscous drag and elasticity forces dominate Brownian forces. That predominance requires that $l_{p}/L\gg 1$ and $Pe=8\pi\mu\dot{\gamma}L^{3}/k_{b}T\gg 1$, where $l_{p}=E/k_{b}T$ is the persistence length of the fiber against thermal fluctuations, and the Péclet number, $Pe$, is the ratio of viscous to Brownian forces. Taking water as the solvent, a fluid strain-rate of ${\dot{\gamma}}=1~{}s^{-1}$ and a fiber of $L=4\,\mu m$ we find $Pe\sim 400$. For fibers of a length greater than a few microns $Pe\gg 1$ and center of mass diffusion can thus always be neglected compared to advection by viscous flow. Now for a material modulus of $Y=1~{}GPa$ and an aspect ratio $\epsilon=10^{-2}$ we find $l_{p}/L\sim 10^{5}$ and Brownian forces can be neglected over elastic forces. This is the case of most of the synthetic fibers treated in this review. Decreasing the aspect ratio to $\epsilon=10^{-3}$ and $Y$ by an order of magnitude reduces $l_{p}/L$ to $l_{p}/L\sim 1$ and shape fluctuations resulting from Brownian forces will become important. This is the case of semi-flexible polymers as for example actin filaments.

Finally, while this dimensionless parameter will appear naturally later in the review, we introduce $\tilde{\eta}=8\pi\mu{\dot{\gamma}}L^{4}/Ec$, where $c=-\ln(e\epsilon^{2})$, as the ration between viscous and elastic forces, a control parameter in many fiber-fluid problems.

II EXPERIMENTAL TECHNIQUES

A limiting factor in experimentation of fiber-fluid interactions is achieving good control on fibers in terms of their shape, dimensions, and mechanical properties as well as good flow control. Even if low Reynolds number experiments can be carried out at macroscopic scales by using highly viscous suspending fluids, it is much more comfortable and advantageous to work at microscopic scales where microscopy and microfluidics can be combined to perform reliable experiments. The microscopic approach requires microfabrication techniques, smart choice of the microchannel geometry, and ways to track the particle and its deformation while transported. All these aspects will be detailed in the following discussion.

II.1 Fiber synthesis and properties

Reliable experiments require fibers that have controlled dimensions and high aspect ratios as well as known and reproducible mechanical properties. Recent fabrication techniques that can be directly implemented into the channels are specially well-suited for this goal.

II.1.1 Fiber fabrication

The fabrication of microparticles with non-spherical shapes has received much attention in recent years, but elongated objects require specific approaches. Electro-spinning can be an efficient technique to fabricate fibers of typical diameters $\sim 200~{}nm$ and lengths of $1-100~{}\mu m$ (Fig. 1i) with controlled mechanical properties (Li & Xia (2004); Jun et al. (2014); Nakielski et al. (2015)). Microfluidics has opened new routes to tailor elongated objects. Using the control over channel geometry given by soft lithography, it is possible to form a jet of a solution that will solidify either by photo-polymerization or by diffusion of small molecules from the surrounding fluid (Figs. 1a&b) (Jeong et al. (2004); Mercader et al. (2010); Choi et al. (2011); Nunes et al. (2012); Perazzo et al. (2017)). Another option is to photo-polymerize a photo-sensitive solution, directly in the microchannel, through a mask which determines the object shape (Figs. 1c-e). The presence of oxygen close to the PDMS walls prevents polymerization in these regions keeping the fabricated object free to flow in the microchannel (Helgeson et al. (2011); Berthet et al. (2016)). Finally, self-assembly of colloids is a reliable technique that produces very flexible, high aspect-ratio particles but their sedimentation due to their high density can be limiting (Fig. 1h) (Goubault et al. (2005); Zhang et al. (2011); Berthet et al. (2016)).

Another attractive strategy is to rely on biology for controlled and monodisperse elongated objects. Actin filaments, semi-flexible polymers resulting from the polymerization of globular protein, are models for flexible Brownian fibers (Harasim et al. (2013); Kirchenbuechler et al. (2014); Liu et al. (2018)). Microtubules are another candidate. Bacterial colonies can also form elongated objects that allow fluid-structure interaction problems to be studied (Amir et al. (2014); Rusconi et al. (2010, 2011)).

More complex shapes can be produced. In 2D, approaches based on self-assembly or photo-lithography can be implemented (Figs. 1c-e and Li & Xia (2004)). 3D fabrication is more challenging, and attention has focused on microhelices given their similarity to bacterial flagella. Electrospinning is a widespread technology (reviewed in Silva et al. (2017)); The spontaneous formation of helices from flat ribbons is an elegant alternative (Pham et al. (2013)) (Fig. 1f). Recent developments in 3D printing using two-photon excitation allows sub-micron resolution but is currently limited to rigid objects (Fig. 1g).

II.1.2 Mechanical Characterization

The mechanical properties of macro-scale fibers can be measured by rheometry or by the bending of the fiber under its own weight (see Fig. 1j). When working with elastomers these two measurements are in good agreement (Quennouz et al. (2015)). Beam-bending measurements can be implemented in a microfluidic channel to perform in situ measurements (see fig.1k) (Duprat et al. (2015)).

When the smallest dimension of the fiber is a micron or less, Brownian fluctuations can be large enough to induce bending deformations. In this case the flexural rigidity $E$ can be expressed as $l_{p}k_{b}T$, with $l_{p}$ the persistence length, and an elegant way to measure $E$ is to extract $l_{p}$ from the analysis of many filament configurations (fig. 1l). This approach is common for biological semi-flexible polymers like actin or microtubules (Gittes et al. (1993)) and has also been used for tiny electrospun filaments (Nakielski et al. (2015)).

II.2 Flow control

At large spatial scales, working at low Reynolds number requires long-time experiments with low-velocity flows of highly viscous liquids. However, at small scales low Reynolds numbers are easily achievable even with low viscosity fluids. For this regime the channel geometry determines completely the flow and thus different flows can be obtained depending on the device design.

To avoid flow complexity, studies are often done using linear extensional or shear flows (Fig. 2a&b). Experimentally, working with linear shear flow is only possible in a Couette geometry if the gap is small compared to the cylinder diameters (Fig. 2a). In channel flows, pure shear can be approximated by Poiseuille flow if the object size is much smaller than the channel height (Fig. 2a) (Liu et al. (2018)). Cross-slot geometries generate extensional flows but with small spatial extension. An elegant but difficult way of overcoming this limitation is to keep the object at the stagnation point by adjusting the pressures in the different entries as was implemented by Schroeder et al. (2003) and later by Kantsler & Goldstein (2012) (see Fig. 2b). A more robust way is to optimize numerically the exact geometry to obtain a purely extensional flow over a larger distance even in the presence of the walls (Fig.2b) (Haward et al. (2012); Zografos et al. (2016)).

Numerical simulations are sometimes done in 2D due to their relative simplicity, while in experiments 2D flows are not always easy to achieve. One option is to use a Hele-Shaw geometry, in which lateral dimensions are much larger than the cross-channel dimension (Fig. 2c). On a microscope, the vertical configuration allows a nearly 2D flow to be generated in the observation plane above the boundary (Haward et al. (2012)).

Flows can be controlled by syringe pumps imposing flow rates, or by pressure control devices. Pressure control is well adapted to low velocity and shear rate. Particle velocimetry techniques are used to quantitatively determine local velocities and thus shear rates.

A new development in the field are microfluidic rheometers that access the rheological properties of complex suspensions of various types (Lindner & Arratia (2016)) (Fig. 2d). Their main advantages are their sensitivity to low variations in viscosity and normal stresses and the possibility to directly observe the flowing particles within the device.

II.3 Particle observations and tracking

Typically, fiber deformations and trajectories need to be measured simultaneously even though they occur at different spatial scales: Images with high-magnification are necessary to obtain good resolution of fiber deformations whereas trajectories are inherently large-scale. In a microchannel, bright field, phase contrast, and fluorescence microscopy are usually best suited to acquire sufficient quality images of fiber deformations. The use of high-magnification and a high-numerical-aperture objective is in general required. In addition, a long working distance can be very useful to observe filaments far from the bottom of the chamber where the coupling to the wall is negligible. At macro-scale, the 3D-shape may be reconstructed from images taken by two cameras or by the smart use of mirrors.

To track an object along its trajectory at the macro-scale, automated motion of the camera can be used. Keeping a fiber in the field of view in a microfluidic channel requires that the microscope stage moves at or near the fiber velocity. In combination with a fast camera this avoids image blurring. Following trajectories in 3D becomes necessary when objects do not follow streamlines as in cross-stream migration, or for active particles. In this case, one can use 3D tracking algorithms like those developed to follow swimming bacteria (Darnige et al. (2017); Qu et al. (2018)).

III MATHEMATICAL MODELING AND SIMULATION METHODS

The interaction of flexible fibers with Stokesian fluids has a special structure for which specialized mathematical and computational methods have been developed. The most basic is local slender-body theory (SBT), which gives a leading-order local relation between elastic and drag forces and is a mainstay for analysis. Nonlocal interactions arising from fluid incompressibility can be captured through use of higher-order, and more complex, slender-body formulations, or through other approaches such as immersed boundary methods, bead-rod models, or regularized Stokeslet methods. Here we first discuss the nonlocal SBT formulation, as it is compact and amenable to fast solution methodologies.

III.1 Nonlocal Slender Body Theory

The appeal of nonlocal SBT is its reduction of filament-fluid and filament-filament interactions to relatively simple dynamics equations for filament centerlines. For this we assume that the fluid is described by the incompressible Newtonian Stokes equation:

[TABLE]

with $\mu$ the shear viscosity, $\mathbf{u}$ the fluid velocity, and $p$ the pressure. Generally, solutions to the Stokes equations can be found using boundary integral methods (Pozrikidis (1992)), where the fluid velocity is represented as a distribution of fundamental solutions to the Stokes equations on all immersed and bounding surfaces. The fundamental solution corresponding to a $\mathbf{\delta}$-function of force is given by the Stokeslet tensor: $\mathbf{G}({\bf x})=\frac{1}{8\pi\mu}\frac{\mathbf{I}+\hat{\bf x}\hat{\bf x}}{|{\bf x}|}$, where $\hat{\bf x}={\bf x}/|{\bf x}|$.

The flows produced by the motion of slender fibers are treated especially as their surface integrals can be reduced, through asymptotics, to integrals of Stokeslets along their center-lines (Keller & Rubinow (1976); Johnson (1980); Götz (2000)). Specifically, consider $N$ fibers, each of radius $a$ and lengths $L_{n}$, with center-line positions ${\bf X}^{n}(s,t)$ for $n=1,...,N$, with $s$ the arclength ($0\leq s\leq L_{n}$). Each fiber is assumed to be slender, i.e., $\epsilon_{n}=a/L_{n}\ll 1$, and to exert a force per unit length ${\bf f}^{n}(s,t)$ on the surrounding fluid. Given a background flow $\bar{{\bf u}}$ the fluid velocity at a point ${\bf x}$ within the fluid is, to leading order, given by distributions of Stokeslets along fiber center-lines:

[TABLE]

Götz (2000) derives higher-order corrections in terms of integrals of Stokes doublets.

The self-induced velocity of a fiber is not found by direct evaluation of Eqs. (2) upon that fiber. This yields a logarithmically divergent integral. Instead, fiber self-induction is a matter of careful asymptotic matching (Keller & Rubinow (1976); Johnson (1980); Götz (2000)). Such analyses give the velocity of the centerline ${\bf X}^{n}(s)$ as

[TABLE]

where $\mathbf{X}_{s}$ is the tangent vector and subscript $s$ denotes a partial derivative. Here $c_{n}=-\ln(\epsilon_{n}e^{2})$ is an asymptotically large parameter, and ${\bf K}_{n}$ is an $O(\epsilon^{0})$ linear operator that contains a ”finite-part” integral of Stokeslets that captures fiber self-interactions. The use of this integral in computation is discussed by Tornberg & Shelley (2004) who show that it must be appropriately regularized both to maintain asymptotic accuracy and to control the growth of short wavelength modes treated inaccurately by slender-body theory.

Equation (3) relates the fiber velocity to the fiber forces acting upon the fluid. As a dynamics problem, this system is closed by specifying the fiber elastic forces and identifying a Lagrangian or material variable for the fiber. These are done simultaneously by modeling each fiber as an inextensible Euler-Bernoulli beam:

[TABLE]

with $-E{\bf X}_{ssss}$ the bending force per unit length, and $T$ the axial tension with $(T\mathbf{X}_{s})_{s}$ the tension force per unit length. Inextensibility identifies $s$ as a material variable, and hence ${\bf X}_{t}^{n}(s,t)={\bf V}^{n}$. Inextensibility also generates an auxiliary integro-differential equation for $T$ by imposing $\mathbf{V}^{n}_{s}\cdot{\bf X}^{n}_{s}=0$ which follows from differentiating the identity ${\bf X}^{n}_{s}\cdot{\bf X}^{n}_{s}=1$ (Tornberg & Shelley (2004)).

III.2 The special case of local SBT

Local SBT has been the mainstay for theoretical analyses of the dynamics of flexible fibers, and sits at the core of many numerical approaches. In particular, retaining from Eq. (3) only the logarithmic leading order in $\epsilon$ we have a local balance of fluid drag forces, and fiber forces exerted on the fluid:

[TABLE]

where ${\bf R}=\eta\left({\bf I}-\mathbf{X}_{s}\mathbf{X}_{s}/2\right)$ and $\eta=8\pi\mu/c$. For an elastic fiber modeled by Eq. (4)) we have a dynamical equation, plus a constraint:

[TABLE]

The inextensibility constraint generates an elliptic equation for $T$:

[TABLE]

where $S$ depends upon background flow and the bending force. Given appropriate boundary conditions, $T$ has a unique solution, and Eq. (6) can be used to evolve the fiber’s shape, position, and orientation.

Basic scaling: Consider a background flow with a characteristic length-scale $W$ and time-scale ${\dot{\gamma}}^{-1}$, expressed as $\bar{{\bf u}}={\dot{\gamma}}W{\bf U}({\bf x}/W,{\dot{\gamma}}t)$. By scaling space on $L$, time on ${\dot{\gamma}}^{-1}$, and tension $T$ on $E/L^{2}$, Eq. (6) can be given in nondimensional form as

[TABLE]

with $\alpha=L/W$, and where $\tilde{\eta}=8\pi\mu{\dot{\gamma}}L^{4}/Ec$ is the effective strength of flow forcing. Note that for linear background flows $\alpha$ disappears from the equation. If the tension is negative, reflecting compressive fluid stresses, then in Eq. (8) is seen as a competition of fourth-order diffusion and a second-order anti-diffusion.

III.2.1 Boundary conditions

Finally, boundary conditions are needed on the tensile and bending forces. For freely suspended fibers, these are the “natural” boundary conditions of having zero applied forces and torques, which simplify to $\mathbf{X}_{sss}=\mathbf{X}_{ss}=\mathbf{0}$, and $T=0$ on the each end of the fiber. These boundary conditions suit the fiber force containing ${\bf X}_{ssss}$ and $T$ satisfying essentially an elliptic boundary value problem. Having filaments attached (clamped or hinged) to surfaces, mobile or immobile, generates other boundary conditions and explores other dynamics, which are detailed in Nazockdast et al. (2017b).

III.3 Brownian semi-flexible fibers

Thus far we have assumed the fibers to be rigid enough that thermal fluctuations have negligible effect upon their conformations. For this to hold, the length of the fibers needs to be significantly smaller than their persistence length $l_{p}$. If this assumption is not met (as can be for carbon nanotubes, actin filaments, and wormlike polymers), thermal fluctuations can have important consequences on the overall mechanics of the suspension, and has been the subject of many studies. A review of this topic is beyond the scope here and we refer the interested reader to Broedersz & MacKintosh (2014). Here we only give a limited summary of previous work and the standing challenges.

In the presence of thermal fluctuations Eq. (5) is modified to

[TABLE]

where $\mathbf{f}^{B}$ is the thermal fluctuating force with zero mean and finite variance given by fluctuation dissipation theorem:

[TABLE]

The brackets denote ensemble averages. Munk et al. (2006) used this formulation to study a 2D Brownian semi-flexible fiber in shear flow. Manikantan & Saintillan (2013) also used this formulation in studies of Brownian flexible fibers in two-dimensional arrays of counter-rotating Taylor-Green vortices and extensional flows.

It follows from Eq. (10) that computing the displacements from thermal fluctuations requires computing ${\bf R}^{1/2}$. For local SBT ${\bf R}^{1/2}$ can be computed analytically. For non-local SBT, the resistance tensor, ${\bf R}$, is nonlocal and depends on the position of all fibers. Numerically, this means computing the square root of a square matrix with dimensions $N_{T}^{2}=(N\times N_{F})^{2}$, where $N_{F}$ is the number of fibers and $N$ is the number of discretization points per fiber. The direct computational of ${\bf R}^{1/2}$ requires $\mathcal{O}(N_{T}^{3})$ operations, which makes simulation of systems beyond a few fibers very expensive. More efficient techniques have been proposed to reduce the computational cost to $N_{T}^{\alpha}$ with $\alpha$ varying between $1<\alpha<2$, depending on the conditions. In a recent development Fiore et al. (2017) use Ewald summation to reduce the computational cost to $O(N_{T})$.

Several other works have modeled Brownian filaments by using bead-rod/spring models for inextensible polymeric filaments. The main complexity involves imposing the inextensibility constraint. Specifically, it turns out that very stiff springs and inextensible rods differ in their degrees of freedom and hence in the resulting equilibrium statistical mechanics (Fixman, 1978). Hinch (1994) showed that a psuedo-potential force needs be applied to a connection with a very stiff spring in order to exactly recover the behavior and statistical physics of an inextensible connection.

III.4 Fast Numerical Methods for SBT

Tornberg & Shelley (2004) developed a stable and numerically tractable version of nonlocal SBT for flexible filaments with free-ends, and devised specialized quadrature schemes for the regularized finite-part integral in ${\bf K}[{\bf f}]$ of Eq. (3). The fibers’ centerlines were discretized uniformly and derivatives evaluated using a second-order accurate finite difference scheme, while trapezoidal rule was used for evaluating the remaining integrals. For time-stepping, a backward differentiation finite difference scheme was used and the bending forces were treated implicitly to remove temporal stiffness, while the tension equation ? which imposed inextensibility of the fibers ? was treated explicitly.

Nazockdast et al. (2017b) extended the work of Tornberg & Shelley (2004) to a platform suitable for simulating a substantially larger number of fiber, to $O(1000)$ or more. They used a Chebyshev basis for the representation of fiber position and to compute high-order derivatives with spectral accuracy. They used implicit-explicit time discretization and treated both bending and tensile forces implicitly. As a result, they observed a three order of magnitude improvement in the maximum stable time-step in their numerical experiments, compared to the explicit treatment of tension by Tornberg & Shelley. They solved the linearized system of equations using GMRES (Saad & Schultz (1986)) with a Jacobi or block Jacobi preconditioner. A kernel independent fast multipole method (FMM) is utilized for fast computation of nonlocal hydrodynamic interactions via fast matrix?vector products, yielding a computational cost of $O(M)$ operations per time-step, where $M$ is the number of unknowns. The entire computational scheme is parallelized and scalable to many computational cores, which allowed simulation of $O(1000)$ semi-flexible filaments.

III.5 The Regularized Stokeslet Method

A popular method for solving the Stokes equations is the method of regularized Stokeslets of Cortez (2001) which uses superpositions of regularized fundamental solutions. Flores et al. (2005) used a superposition of regularized Stokeslets and Rotlets to simulate the dynamics of driven flagellae. In their study, a flagellum is a network of flexible springs, and a helical shape so composed is driven by a torque at its base. In a study also modeling flagella (and cilia) Smith (2009) developed a version of regularized Stokeslets that utilizes a boundary-element approach to the discretization and incorporates the presence of walls (see also Jiang et al. (2014)). Bouzarth et al. (2011) use a regularized representation of a one-dimensional curve of two-dimensional Stokeslets to simulate the nonlocal dynamics of flexible, slightly extensible fibers. Olson et al. (2013) have combined the regularized Stokeslet method for evolving slender rods with the internal mechanics of an elastica with intrinsic twist and curvature. And, in theoretical work, Cortez et al. (2012) have developed a asymptotically consistent slender-body theory based on using regularized center-line forces.

III.6 The Immersed Boundary Method

The immersed boundary method (IBM) of (Peskin, 2002) has also been applied to this class of problems. In this method, a filament is discretized by Lagrangian markers connected by spring elements, and their relative displacements by fluid motions are used to calculate the filament’s elastic response. These elastic forces are then distributed onto a background grid covering the fluid volume, and are used as forces acting upon the fluid, thus modifying the fluid flow that displaces the markers (Lim & Peskin, 2004). The advantage of the IBM is that detailed immersed mechanical structures can be simulated, but at the cost of having to solve the flow equations in the entire fluid volume. Other volume-based methods have also been used in this context, including finite element/volume methods (Mitran, 2007), and Lattice Boltzmann (Wu & Aidun, 2010).

Within the immersed boundary framework, for a fiber to have a physical width the volume grid size needs to be several times smaller, making the method computationally expensive for simulating slender bodies such as fibers and thin disks. To overcome this, Nguyen & Fauci (2014) employed an adaptive-mesh version of the IBM (Griffith et al., 2007) to study the dynamics of flexible fibers as models, in part, for diatom chains (diatoms are unicellular phytoplankton). For this, they constructed fibers that are composite structures made of alternating segments that mimic the structure of diatom chains (see Fig. 4c). While a simple beam model, like Eq. (8), describes its bending deformations well, the results of compressive strain are not. Wiens & Stockie (2015) used their parallel implementation of the IBM to simulate a suspension of up to 256 semi-flexible fibers in shear flow. Their analysis of parallel performance shows their algorithm is weakly scalable. Recently, Stein et al. (2017) implemented a fourth-order accurate IBM for studying flows over arbitrary smooth domains. A combination of this with adaptive-grid immersed boundary (Griffith et al., 2007) might improve the computational efficiency and accuracy of immersed boundary in simulating slender objects in flows.

III.7 Bead-and-Rod Models

While closely related to methods based upon SBT or the immersed boundary method, bead-and-rod models have a somewhat independent lineage. In bead-and-rod models, a flexible fiber is represented as a one-dimensional chain of linked rigid bodies (e.g. spheres, spheroids, rods) that experience local drag and interact with each other through short-range forces (e.g. repulsion, lubrication, friction), while sometimes neglecting long-ranged hydrodynamic interactions or only including a subset of them. A recent review of these methods is found in Hämäläinen et al. (2011). In recent work Delmotte et al. (2015) have elaborated upon basic bead models by introducing a new Lagrange multiplier method to impose constraints, and consider several flow problems – Jeffery orbits, buckling in shear, actuated swimming filaments – by using an approximate accounting of Stokesian hydrodynamics.

IV FREE FIBERS IN FLOW

How fibers are buckled by flow is central to much of the interesting nonlinear dynamics observed in simulations and experiments of fiber motion. Here we first discuss a prototypical situation – a straight fiber moving in a linear straining background flow – where buckling arises as an instability. We then discuss other prototypical problems such as morphological transitions of fibers transported in linear shear flows and more complex flow geometries where the coupling of fiber deformation and transport can be observed.

IV.1 Fiber morphology in simple flows

Here we discuss the shape dynamics of fibers in linear shearing and straining flows.

IV.1.1 Buckling instabilities in compressive flows

Because of its simplicity, local SBT is usually the preferred tool for studying linear stability of immersed fibers and is used here to illustrate the buckling instability of free fibers in the simplest of cases.

A perfectly straight isolated fiber will remain straight in any incompressible and linear flow, even as the fiber is translated and rotated. Hence, consider the background flow ${\bf U}({\bf x},t)={\bf A}(t)\cdot{\bf x}$ with $tr({\bf A})=0$. A straight fiber can be represented as ${\bf X}(s,t)={\bf X}_{c}(t)+s{\bf e}(t)$ with ${\bf X}_{c}$ its center point, and ${\bf e}$ a unit orientation vector. Inserting this form into Eqs. (8) & (7), and applying the “natural” conditions of zero total force and torque, yields

[TABLE]

where ${\bf E}=({\bf A}+{\bf A}^{T})/2$ is the symmetric rate-of-strain tensor. Hence, the rod center is carried with the local flow, the orientation vector obeys Jeffrey’s equation (Jeffrey (1922)), and the tension ${\bar{T}}$ is quadratic in $s^{2}$ with its sign determined by the orientation of ${\bf e}$ relative to the principle axes of ${\bf E}$. Thus, if the fiber is aligned with compressive straining of the flow then the tension is negative and hence compressive. This is the condition required for a buckling instability.

The most straightforward case is given by a simple straining flow $\bar{{\bf u}}=(x,-y)$ where the straight fiber is aligned with the $y$-axis (i.e., ${\bf e}\equiv{\hat{\bf y}}$). This choice maximizes flow compression (Young & Shelley (2007)) and is free of explicit time-dependence from the background flow. This leads to a variable coefficient eigenvalue problem for the growth or decay of an in-plane perturbation $w$ to the straight fiber:

[TABLE]

where $w_{ss}=w_{sss}=0$ at $s=\pm 1/2$. With an eigenvalue solver we can track the system’s eigenvalues and eigenfunctions as $\tilde{\eta}$, the effective viscosity or strain-rate, is increased. For small $\tilde{\eta}$ the straight fiber is stable to perturbations. With increase in $\tilde{\eta}$ we find the successive crossing to the right half-plane of eigenvalues coupled to eigenfunctions associated with increasingly higher order bending modes. The first three crossings occur at $\tilde{\eta}_{1}=153.2$, $\tilde{\eta}_{2}=774.3,$ and $\tilde{\eta}_{3}=1930$ (Quennouz et al. (2015)), and their associated eigenfunctions are shown in Fig. 3a (left panel). The respective eigen-shapes are the classical U-, S-, and W-shaped buckling modes.

In experiment it is nontrivial to achieve a linear straining flow and reasonably long residence times for a fiber moving in it. One possibility is to use the approach of a fiber to a stagnation point. Fiber buckling has been investigated in a macroscopic realization of such a system by Wandersman et al. (2010) who used centimetric soft elastomer fibers. These fibers move in a viscous quasi-2D cellular flow, an array of counter-rotating vortices surrounding hyperbolic stagnation points. Above a critical value of $\tilde{\eta}$, fibers are observed to buckle (Fig. 3b (right)) at or near stagnation points. With increasing $\tilde{\eta}$ more complex fiber shapes were observed (Fig. 3a (right)), corresponding to the eigen-shapes predicted by the linear stability analysis. Fiber trajectories from experiments and simulations (3b) show fiber transport across the cellular array. Due to the complexity of the flow, buckling does not take place at each passage of a stagnation point, but depends on the exact conformation and position of the approaching fiber (Quennouz et al. (2015)).

Buckling instabilities have also been investigated in the presence of Brownian fluctuations. Manikantan & Saintillan (2015) have shown numerically that Brownian fluctuations play a minor role on the buckling instability and serve mainly to broaden the threshold of the transitions towards the different buckling modes. Kantsler & Goldstein (2012) investigated the deformation of a micrometric actin fiber held at a stagnation point created in a microfluidic cross-slot device. Their experiments as well as numerical simulations by Manikantan & Saintillan (2015) show modes qualitatively similar to the deterministic predictions, albeit with shape fluctuations owing to Brownian fluctuations (see Figs. 3c, d).

The optimal flow geometry for the investigation of buckling instabilities is a hyperbolic channel, where constant extension rates can be achieved over long residence times. This was realized recently using actin filaments transported in an optimized microfluidic flow geometry (sec. II.2) and preliminary results can be seen on Fig. 3e. The filament is first stretched in the converging part of the channel where Brownian fluctuations are suppressed by the viscous forces. Then the filament enters the diverging part of the channel where it experiences strong compression leading to higher buckling modes. A similar succession of stretching and coiling has been observed in microchannels with constrictions for actin filaments (Strelnikova et al. (2017)) and polymeric microfibers (Nunes et al. (2012)).

IV.1.2 Fiber morphology in shear flows

A linear shear flow is the superimposition of a rotational and a straining flow and, as a consequence, fiber dynamics are a combination of tumbling and periodic deformation. A rigid ellipsoid tumbles periodically in a linear shear flow, while resisting stretching or compression as it rotates through the extensive and compressive quadrants. A flexible fiber on the other hand will react to the straining flow, leading to a transition to a buckling instability as $\tilde{\eta}$ is increased (Becker & Shelley (2001); Liu et al. (2018); Forgacs & Mason (1959a)). When $\tilde{\eta}$ is increased further, fibers become too flexible to sustain rotation of the fiber as a whole and a transition to ”snaking” motion is observed. Here the fiber remains mainly aligned with the flow direction, but performs a tank trading motion (Forgacs & Mason (1959b); Harasim et al. (2013); Liu et al. (2018); Nguyen & Fauci (2014); Delmotte et al. (2015)).

Early experiments by Forgacs & Mason (1959b) investigated the deformation of millimetric elastic fibers in corn syrup under linear shear in a Couette cell. They were the first to address fiber buckling (Forgacs & Mason (1959a)) and observed what they baptized ”springy rotations” (Forgacs & Mason (1959b)). For longer fibers they identified ”snake” turns as shown on Fig. 4a. The transition to fiber buckling was clearly identified only later in a numerical investigation by Becker & Shelley (2001) for a fiber rotating in a linear shear flow. The first transition to buckling takes place at $\tilde{\eta}=306.4$. This is twice the value $\tilde{\eta}_{1}=153.2$ found for straining flow, as the magnitude of the straining part in shear flow is proportional to only half of the magnitude of the shear flow. The two threshold values are thus effectively identical. Using the local SBT formulation, their numerical simulations studied the very nonlinear shape dynamics of the fiber above the buckling transition (see Fig. 4d). Experiments by Harasim et al. (2013) using actin filaments in microscopic channel flows have investigated in detail the ”snake” turn motions and were the first to describe and model this configuration as a tank trading motion where a narrow bend travels along the filament whose two ends remain mainly aligned with the flow. They also addressed the role of rotational diffusion on the period of Jeffery orbits. Numerous investigations have observed the ”snaking” dynamics numerically (Nguyen & Fauci (2014); Stockie & Green (1998); Delmotte et al. (2015)) (see Fig. 4c). Only recently have all subsequent regimes been observed (Fig. 4e) and classified (Fig. 4f) as a function of $\tilde{\eta}$ and $\l_{p}/L$, by experiments using actin filaments and a numerical investigation using non-local SBT in the presence of Brownian fluctuations (Liu et al. (2018)). Special care was taken when designing the experimental flow geometry to assure linear shear conditions (sec. II.2). The authors also addressed the transition to the ”snaking” mode which is initiated by a buckling instability occurring, under sufficiently strong shear, in a filament that has not yet aligned with the compressive axis. Rotation of the deformed fiber will then align one end with the flow direction whereas the other can bend. The theoretically predicted threshold is in good agreement with experimental observations (Fig. 4f). Well above the onset of the ”snaking” motion more complex dynamics such as coiling (Nguyen & Fauci (2014)), knot formation (Kuei et al. (2015)), or helical motion (Forgacs & Mason (1959b)) were also reported.

IV.2 Fiber morphology and transport in more complex flows

Here we discuss the shape dynamics of fibers in more complex background flows, such as those arising from confinement or imposition of a background flow through external forcing.

IV.2.1 Poiseuille flow and interaction with bounding walls

Channel (or Poiseuille) flows represent the majority of experimental systems and, unlike linear shear flows, exhibit a shear gradient across the channel width. This can affect fiber morphology when fiber lengths are comparable to the channel dimensions. A first consequence is that buckling instabilities (as described above) are not observed and bending modes dominate. Typical fiber morphologies in Poiseuille flows are reminiscent of ”snaking” motions (Harasim et al. (2013)) (Fig. 4b) and can vary depending on the exact positioning within the flow (Reddig & Stark (2011); Chelakkot et al. (2010)) (Fig. 5a,b). Secondly, several groups have demonstrated, using bead numerical models (Reddig & Stark (2011); Chelakkot et al. (2010)) or experiments with actin filaments (Steinhauser et al. (2012),Fig. 5c), that semi-flexible polymers show strong cross-stream migration in shear gradients and accumulate at preferred distances from channel walls. Similar predictions have been made for flexible fibers (Slowicka et al. (2013)). Very recently Nakielski et al. (2015) observed cross-stream migration using Brownian electrospun hydrogel nanofilaments in an oscillatory flow (Fig. 5d) and highlighted the dependance of migration dynamics on the initial fiber positioning.

IV.2.2 Confined flows

Fiber transport in a Hele-Shaw cell (Fig.  2c) represents a very different situation. Fibers with a diameter comparable to the cell height are confined by the top and bottom walls and now evolve in 2D in a plug flow. Such confined microfabricated fibers (sec. II.1) have been observed to bend into a C-shape facing backwards relative to the flow direction (Fig. 5e). At first sight it is surprising that a fiber transported at zero force in a plug-flow deforms, but this can be explained as follows: friction with the top and bottom walls decreases the fiber velocity compared to the average flow velocity and the fiber can be seen as a moving obstacle. The fiber thus creates a flow perturbation inducing an inhomogeneous force distribution along the fiber. Cappello et al. (2018) have evaluated the force distribution as a function of fiber confinement using modified Brinkman equations. The resulting fiber shapes are found to be in good agreement with experimental findings (Cappello et al. (2018)). The coupling between fiber deformation, drift, and interaction with the side walls (Nagel et al. (2018)) leads to complex trajectories and finally a preferred alignment with the flow direction (Fig. 5e).

IV.2.3 Cellular flows

From their numerical study of fiber-flow interactions Young & Shelley (2007) predicted that flexible fibers could move as random walkers across a closed stream-line cellular flow (Fig. 5h). Fiber buckling, occuring close to the stagation points, yields due to its many degrees of freedom an effective randomness, allowing the fiber to cross streamlines and to be transported across the flow. The precise transport properties depend in a non-trivial way on fiber flexibility and length, as has also been shown in experiments (Wandersman et al. (2010)). Numerical simulations by Manikantan & Saintillan (2013) addressed the role of Brownian fluctuations which were shown to increase trapping of filaments in vortices and thus decrease transport across the vortex array. Deng et al. (2015) have numerically confirmed these findings. See Bouzarth et al. (2011) and Young (2009) for a numerical study of two fibers interacting in such cellular flows and Yang & Fauci (2017) for transport in a polymeric cellular flow (Fig. 5f.)

V SINGLE FIBERS MOVING UNDER EXTERNAL FORCES AND CONSTRAINTS

V.1 Sedimenting Fibers

When flexible fibers are subject to external forces such as gravity, the coupling between the elastic and viscous forces can induce complex deformation and reorientation that strongly affects their macroscopic transport. Sedimentation is the classical example (though centrifugation would have similar consequences). The sedimentation of flexible fibers was first studied from the theoretical point of view by Xu & Nadim (1994). Numerical studies have since been conducted using bead models (Cosentino Lagomarsino et al. (2005); Schlagberger & Netz (2005); Delmotte et al. (2015); Saggiorato et al. (2015)) or slender-body theory (Li et al. (2013)). Only very recently have experimental studies been performed (Marchetti et al. (2018); Raspa et al. (2018)).

During sedimentation, fiber weight and viscous drag balance, and a rigid fiber settles at constant velocity and orientation. Hydrodynamic interactions along the fiber yield a viscous drag that is smaller in the middle than at the edges. For a flexible fiber, this inhomogenous drag deforms the fiber and leads to rotation. The theoretical work of Li et al. (2013) first explored how this affects settling trajectories. As the drag depends on fiber orientation, the rotation leads to a nontrivial coupling between horizontal and vertical translations (Fig. 6a). This occurs until the fiber has reached an equilibrium configuration characterized by a constant settling velocity and a bent shape, the middle being perpendicular to gravity independently of the release conditions (Fig. 6b). Saggiorato et al. (2015) have explored numerically the behavior of very flexible and very heavy filaments whose shape can become non-planar resulting in drift or helical trajectories (fig.6c).

The fiber shape at equilibrium depends on the relative magnitudes of the gravitational force and the elastic restoring force. The relevant control parameter is $B=F_{g}L^{2}/E$, with $F_{g}$ the gravitational force (Fig. 6d) (Cosentino Lagomarsino et al. (2005); Delmotte et al. (2015); Saggiorato et al. (2015); Marchetti et al. (2018)). For small $B$ the steady velocity is close to the horizontal settling velocity of a rigid fiber, while for large $B$, the fiber experiences large deformations and the velocity is close to the vertical settling velocity of a rigid fiber (Fig. 6e) (Delmotte et al. (2015)). Marchetti et al. (2018) recently confirmed this experimentally and identified the intermediate regime as one of continuous drag reduction through fiber reconfiguration. Defining an elasto-viscous number $\tilde{\eta}$ is not straightforward for this problem as the viscous forces are function of fiber shape.

The case of a fiber settling while aligned with gravity is very different: a inhomogenous fiber tension results from the inhomogenous viscous drag, and is compressive in the leading-half and extensile in the trailing part (Li et al. (2013)). Above a certain threshold of the control parameter, buckling may occur in the leading half (Fig. 6f).

V.2 Anchored fibers

Passive anchored fibers, such as the primary cilium (Young et al. (2012)), are found in biological systems, and can also form spontaneously under flow conditions as in the formation of biofilm streamers (Drescher et al., 2013). In engineered micro-fluidic flow geometries, fabricated anchored fibers can be used as flow sensors, or conversely the micro-fluidic flows can be used to measure the bending properties of unknown materials. Depending on the anchoring and flow orientation either buckling or bending can occur.

V.2.1 Fiber buckling

Using local SBT, Guglielmini et al. (2012) investigated the stability of elastic fibers when held against an impinging linear or quadratic stagnation point flow (Fig. 7a). One fiber end is either clamped or hinged at the wall (i.e. free to rotate with zero applied torque). To identify critical values in $\tilde{\eta}$ for buckling transitions they discretized the linearized dynamics equations and its boundary conditions and posed it as a generalized finite-dimensional eigenvalue problem for growth rates. For the clamped filament a first unstable bending mode was identified for both flow fields. The next unstable mode (in $\tilde{\eta}$) corresponds to buckling (Fig. 7a) with slightly different thresholds for the two flows. The hinged fiber is always unstable to rotation around its base, and its higher modes correspond to buckling, having significantly lower thresholds as compared to free fibers.

Experimental studies of fiber buckling in impinging flows are difficult as even a slight misalignment between the fiber base and the stagnation point leads to fiber bending (or rotation) instead of buckling. However, situations where long fibers remain either temporarily pinned, as observed for fibers flowing in rough fractures (D’Angelo et al. (2009)), or which experience a compressive flow only locally, as observed during continuous fiber fabrication in a widening microfluidic device (Nunes et al. (2013)), show buckling instabilities similar to those of anchored fibers. In particular, Nunes et al. (2013) show that such buckling can be used to form ”crimped” microfibers with a well defined morphology (Fig. 7b).

V.2.2 Fiber bending

The deformation of anchored fibers by viscous flows is, in a number of simple flow geometries, akin to a bending beam where viscous forces play the role of a gravitational load. Microscopic beam bending has been used by Duprat et al. (2015) to determine the bending modulus of an elastomeric fiber fabricated using the UV projection method introduced in Sect. II.1. Here the fiber was confined between the top and bottom channel walls and held on both ends in a straight microchannel (Fig. 1l). In this configuration flow takes place through the small gap above and below the fiber and the viscous force applied can be determined using a simple lubrication analysis.

The flow geometry is more complex when a confined fiber is attached to one wall and so partially blocks a confined channel. The resulting flow profile has been visualized (Fig. 7c) in experiment by Wexler et al. (2013) using fibers attached to the wall of a micro-channel and perpendicular to the flow direction. A theoretical model (Wexler et al. (2013)) gave insight into the competition between bending and leakage flow, showing good agreement with their experimental results, and was also used to measure fiber rigidity. By knowing its mechanical properties, such a fiber could be used as a flow sensor in micro-fluidic devices.

Amir et al. (2014) studied the flow-induced bending of single-cell Escherichia coli growing from slits in the side-walls of a micro-fluidic channel (Fig. 7d). By applying a flow perpendicularly to the cells in their unconfined flow geometry, the experimental set-up corresponds again to a simple bending beam experiment where the force applied to the cell results from the viscous friction of the fluid. By estimating this viscous force and using linear elasticity theory, the authors were able to estimate the bending stiffness of the E. coli cells from the measured deflection. The simple hydrodynamic set-up gave values similar to those obtained using much more costly techniques such as AFM measurements.

Another biological example is the primary cilium. The primary cilium is a non-motile and flexible hair-like protrusion from a cell into the extracellular space and is found in a variety of vertebrate cells. The primary cilium is believed to act as a mechanoreceptor by bending in response to flow as is observed in kidney tubule cells. Its dynamics has been investigated in a combined experimental and theoretical study by Young et al. (2012).

More complex flow geometries are experienced by biofilm streamers that have been shown to grow into long, filamentous elastic structures (Rusconi et al., 2010, 2011) under flow. The shape of elastic fibers within flows with curved streamlines has been determined numerically by Autrusson et al. (2011) for fibers anchored, either hinged or clamped, at given positions near a two-dimensional corner. This work shows that, due to tension and bending forces within the fibers, the fibers do not align with the flow but rather cross flow streamlines. This is in agreement with the experimental observations from Rusconi et al. (2010, 2011) shown on Fig. 7e.

V.2.3 Stretching of helical fibers

The stretching by flow of very flexible anchored microhelices has been studied by Pham et al. (2015). Here the viscous drag represents a distributed loading leading to a nonuniform pitch distribution (Fig. 7f). The overall helix extension is well described by the FENE spring model, with the linear part simply given by the ratio of viscous to bending forces. The micro-scale helical ribbons with nano-scale thickness serve both as model systems for bacteria flagella and as very sensitive flow sensors.

V.3 Actuated fibers

Biological cilia and flagella are of central importance to biophysical processes like propulsion or embryogenesis (Bray (2001); Hirokawa et al. (2006); Babu & Roy (2013)). Both are elongated fibers anchored to cellular membranes and which beat or rotate through the action of internal molecular motors (like the cilia covering the surface of micro-organisms; see Camalet et al. (1999)) or by rotating motors anchored at their base (like bacterial flagella; see Bray (2001)). Their interactions with the surrounding fluid medium lead to shape deformations to useful work like transport processes like propulsion or mucus clearance in the lung. The study of anchored filaments submitted to forces is also very important to understand the coupling between cytoskeletal filaments and molecular motors that occurs in the cell (e.g. Trong et al. (2015); Shelley (2016)) during intracellular trafficking or to design micro-robots capable of motion in a fluid.

In the large literature on this topic, largely beyond the scope of this review, fibers are usually anchored to a wall but in some studies fibers are attached to an object small enough to be moved by the fiber motions (i.e. flagella and bacteria). This was studied at the micron-scale in two different experiments (Figs. 8a-b). Early to the field, Wiggins et al. (1998) attached an actin filament to a bead that was periodically driven by an optical trap. The bead oscillations induced traveling waves on the fiber that generated propulsive force. In another realization of this phenomenon Dreyfus et al. (2005) built an artificial micro-swimmer from a chain of magnetic beads attached to a red-blood cell - only here to break symmetry - and actuated the chain with a time-dependent magnetic field.

In another class of studies, fiber motion is induced by the rotation of the anchoring point. Wolgemuth et al. (2000) took inspiration from biology to model the deformations induced by rotation of a fiber around its axis in a viscous fluid. Rotation, twisting and bending combine to yield different regimes of rotation, one of them having the free end repelled from the axis. The interaction of two rotating fibers has been recently studied numerically and experimentally by Man et al. (2017) (Fig. 8d), who show that the induced rotational flows tend to bend the fibers around each other. If a fiber is tilted, at the anchoring point, relative to its rotation axis, the interaction with the fluid leads to the collapse of the fiber onto the rotation axis (Fig. 8c) while assuming a helicoidal shape that generates a propulsive force (Manghi et al. (2006); Qian et al. (2008); Coq et al. (2008)). The rotation of a flexible helix can result in a propulsive force sufficiently large to induce buckling (Fig. 8e, Jawed et al. (2015)). The flexibility and buckling of the flagellar hook at the base of the flagellum plays an important role in the swimming of mono-flagellated bacteria where its buckling yields changes in swimming direction (Son et al. (2013)). Interactions between two weakly flexible helices result in bundling as has been shown experimentally by Kim et al. (2003). As a model for cilia, Guo et al. (2018) studied the hydrodynamic coupling of immersed elastic fibers driven at their base by an applied torque. For two fibers their simulations and theory showed that bistability of in-phase and out-of-phase oscillations can arise from a coupling of elasticity and drive to flow.

Forced fibers can either be actuated locally or all along the body. The latter case has been experimentally and numerically modeled by magnetic colloid chains (Dreyfus et al. (2005); Gauger & Stark (2006); Babataheri et al. (2011)). The former has been addressed in different experimental and numerical studies mainly with the forcing being applied to the anchoring point (Wolgemuth et al. (2000); Wiggins et al. (1998); Manghi et al. (2006); Qian et al. (2008); Coq et al. (2008); Guo et al. (2018)), as already discussed, and typically yield oscillations. De Canio et al. (2017) recently studied tangential forcing at the fiber free-end as a model for the action of motor proteins on microtubules, and demonstrated a bifurcation to fiber beating via a buckling instability (Fig. 8f-g).

VI MICRO-MACRO COUPLING AND COLLECTIVE BEHAVIOR

We have thus far focused on single fibers (or a few) interacting with flows. However, in most industrial and natural settings flexible fibers are mostly found collectively and in suspension. Understanding the transport and rheology of fibrous suspensions at macroscopic length scales is of significant importance. In this chapter we review studies on the relationship between microscopic dynamics of suspended fibers and macroscopic rheology, as well as situations where fibers act collectively in their interactions with flow as in ciliary beds or in new flow control devices.

The main difficulty in studying fibrous assemblies or suspensions is the complex and long-ranged nature of hydrodynamic interactions (HIs). The importance of fiber interactions is defined by their average separation distance. In the dilute regime the fibers separation distance is much larger than their length: $nL^{3}\ll 1$ where $n$ is the fiber number density. In this limit HIs can be neglected and for a suspension the ensemble behavior can be estimated from the dynamics of a single fiber interacting with the flow. Increasing the concentration causes HIs to be stronger. In the semi-dilute regime of $nL^{3}>\approx 1$, the distance between fibers becomes less than their length and both HIs and possibly steric interactions become important. With further increase of concentration the separation distance becomes comparable to the fiber radius and rotational motions are hindered. For suspensions this can result in formation of liquid crystalline phases (Doi & Edwards (1978)), which is beyond the scope of this review.

VI.1 Rheology of fiber suspensions

VI.1.1 Linear rheology for dilute non-Brownian suspensions

The analysis of rheology for non-Brownian suspensions in oscillatory background shear flows is very similar to the buckling instability analysis in simple flows in §IV.1.1. For small perturbations of fiber shape from a straight rod, whose dynamics are described by Eq. 11, Eq. (6) reduces to

[TABLE]

where $w$ is the perturbation transverse to the fiber orientation in the plane of shear, and $\tilde{\omega}$ is the oscillation frequency normalized by the maximal shear rate $\dot{\gamma}$.

Perturbations will decay so long as $\tilde{\eta}<306.4$, the first bifurcation to buckling for the constant shear case. The long-time behavior will then be that of an oscillating rigid fiber which, because of time-reversibility, will yield only a viscous response under time average. For values of $\tilde{\eta}$ beyond the buckling transition fibers under compressive stress will buckle, with that buckling relieved in the reverse part of the cycle. Since elasticity is relaxational throughout the oscillation it remains unclear whether the long-time behavior supports persistent buckling. If not, the rheology will remain purely viscous.

VI.1.2 Linear rheology of dilute Brownian suspensions

The linear viscoelastic behavior of semi-flexible polymers has been widely studied over the past two decades. Examples of these systems include worm-like polymers, cytoskeletal filaments and nanotubes. A detailed discussion of the literature in this limit is beyond the scope of this review, and readers are referred to Shankar et al. (2002). Here we give a very brief summary relevant to this review.

Rotational diffusion is the slowest time-scale for Brownian semi-flexible polymers, and is given by $\tau_{\text{rot}}\sim\mu L^{3}/(k_{b}Tc)$. Another time-scale arises from curvature relaxation, given by $\tau_{\text{curv}}=8\pi\mu L^{4}/Ec$. The ratio of bending to rotational diffusion time-scales is $\tau_{\text{curv}}/\tau_{\text{rot}}=L/l_{p}\ll 1$. The fastest relaxation time is associated with tensile forces. In a rigid fiber the tension is set instantaneously. The presence of thermal forces induces shape fluctuations in semi-flexible polymers. Thus, the average end-to-end length of the polymer will be less than its full contour length and will change under external flows/forces. The tension time-scale is roughly defined as the time $\tau_{T}$ over which the end-to-end distance of the polymer reaches mechanical equilibrium under small external forces Shankar et al. (2002). The tension dynamics is more involved, but has been studied in depth (Morse (1998); Everaers et al. (1999); Pasquali et al. (2001)). These studies find that $\tau_{\text{T}}/\tau_{\text{curv}}=(L/l_{P})^{4}\ll 1$. Thus, $\tau_{\text{T}}\ll\tau_{\text{curv}}\ll\tau_{\text{rot}}$. This separation of time-scales, and the linear response in small deformation allows for analytical solutions of Eq. (9) in §III.3, which is outlined in Shankar et al. (2002).

Having calculated the shape of the fiber, the time-dependent particle stress tensor is found by the force moment integral along the fiber: $\mathbf{S}=\int_{0}^{L}\mathbf{X}(s)\mathbf{f}(s)ds$, recalling that $\mathbf{f}=\mathbf{f}^{E}+\mathbf{f}^{B}$ where $\mathbf{f}^{E}=(T\mathbf{X}_{s})_{s}-E\mathbf{X}_{ssss}$, and $\mathbf{f}^{B}$ is given by Eq. (10). In dilute suspensions the interactions between the fibers are ignored, and the ensemble average stress reduces to $\mathbf{\Sigma}_{F}(t)=n\langle\mathbf{S}\rangle$. Shankar et al. (2002) used this expression to compute elastic and loss moduli in oscillatory shear flow for different values of $L/l_{p}$. Fig. 9a gives a schematic presentation of the tensile, bending, and thermal contributions to the time-dependent shear modulus relaxation in a step strain deformation (see Fig. 2 in Shankar et al. (2002)). At very short times where $t\ll\tau_{\text{curv}}$, the shear modules and its relaxation is determined by tensile forces, while for $t\sim\tau_{\text{curv}}$ all forces contribute equally to the rheology. Finally, at long times the behavior is dominated by thermal forces; thus the long-time behavior is viscous.

VI.1.3 Nonlinear rheology in shear flows

Introducing flexibility to the fiber breaks time reversal symmetry of Stokes equation and changes the behavior qualitatively. Becker & Shelley (2001) used local SBT to show that fiber buckling (see Fig. 4d) gave rise to positive first normal stress differences, and that the predicted threshold to buckling agreed well with the onset of positive first normal stress differences in shearing experiments of nylon fibers in glycerin (Goto et al., 1986).

Tornberg & Shelley (2004) confirmed this observation in their nonlocal SBT simulations. Figure 9a shows the numerical predictions of Tornberg & Shelley for first normal stress difference induced by a single flexible fiber in shear flow, $N_{1}=\mathbf{S}_{xx}-\mathbf{S}_{yy}$, during one rotation of the fiber (solid line). The results for a rigid rod are shown for comparison. To understand the physical basis for this behavior it helps to study the fiber’s shape in three instances: (i) the beginning of deviation from rigid rod, (ii) when the deviations are maximized, and (iii) the return to a rigid rod behavior. These instances are labeled with open circles and are shown with their corresponding shapes in Fig. 9a. We see clearly that the deviation coincides with the buckling of fibers under compressional stress ($\mathbf{p}\cdot\mathbf{E}\cdot\mathbf{p}<0$). As the fiber rotates to the extensional quadrant ($\mathbf{p}\cdot\mathbf{E}\cdot\mathbf{p}>0$), it is stretched back by the straining component of the shear flow to its initial straight shape, and the corresponding stress approaches that of a rigid rod. Intuitively, buckling reduces the effective hydrodynamic length (the end-to-end distance) of the fiber, decreasing the stress magnitude in comparison with a rigid rod (recall that $|\mathbf{S}|\sim L^{3}$). Once the fiber moves to the extensional flow quadrant, it straightens and its effective length increase causes the stress magnitude to approach that of a rigid rod.

VI.1.4 Shear rheology of semi-dilute suspension

Several aspects of the rheology of flexible fibers in semi-dilute regime remains poorly understood, due to the complexity of long-range HIs. To our knowledge, there is no large-scale simulation of flexible fibers in shear flow, with the exception of Tornberg & Shelley (2004) who studied the suspension of $N_{F}=25$ flexible fibers in oscillatory shear flow. They observed that, at large enough shear rates, the fibers deform during the period they are in compressional quadrants of the flow. The fibers are then stretched back once they move to the extensional quadrants. However, they don’t entirely relax, perhaps because of fiber-fiber interactions, and maintain strong buckling into the ensuing periods. They also observe that the total stored elastic energy is a strong function of the initial configuration of the filaments, which the authors relate to the small number of fibers in their simulations. For the same reason they did not report on the rheology of the suspension. To our knowledge, there is currently no prediction of shear viscosity and normal stresses vs shear rate.

On the experimental front, Huber et al. (2014) studied the dynamics of individual actin filaments in semi-dilute solutions of mixtures of labelled and non-labelled filaments flowing in microchannels. They showed that steric hindrance reduces the tumbling rate significantly in semi-dilute solutions as compared to isolated filaments. Entanglements between filaments can also lead to anomalous or aborted tumbling events. The authors link the microscopic filament dynamics to the shear thinning viscosity of the semi-dilute solution via an effective diffusion coefficient.

The effect of entanglements on the three-dimensional deformations of actin filaments in a semi-dilute solution was investigated experimentally by Kirchenbuechler et al. (2014). They used confocal microscopy on a cone-and-plate counter-rotating shear cell and were able to observe simultaneously microscopic filament dynamics and macroscopic flow properties (Fig. 9b). They showed that shear flow causes stretching and alignment of hair pin shaped filaments explaining the observed shear thinning.

A very different behavior is observed by Perazzo et al. (2017) for concentrated suspensions of very long fibers. Using a microfluidic fabrication method (Fig. 1i) they produced concentrated suspensions of fibers with an aspect ratio greater than 100 and showed that gelation, resulting from the formation of mechanical interlocking of chains, is triggered by flow. This is also a very simple way to produce biocompatible hydrogels.

VI.2 Sedimentation of fiber assemblies.

The sedimentation of particles in fluids has been a major area of research in suspension mechanics over the past several decades (Guazzelli & Hinch, 2011). There have also been several experimental (Herzhaft & Guazzelli, 1999; Metzger et al., 2005) and numerical (Butler & Shaqfeh, 2002; Saintillan et al., 2005; Tornberg & Gustavsson, 2006) studies on sedimentation of suspensions of rigid fibers. The complex behavior in sedimenting suspensions is due to HIs (Ramaswamy, 2001; Butler & Snook, 2018). For instance, sedimenting fiber suspensions form inhomogeneous clusters resulting in enhanced sedimentation rates of those clusters, in both numerical and experimental studies (Metzger et al. (2005); Saintillan et al. (2005); Tornberg & Gustavsson (2006)). Manikantan et al. (2014) studied the effect of fiber flexibility under small deformations on sedimentation using a continuum theory, later followed by particle simulations (Manikantan & Saintillan, 2016). These studies showed that flexibility can have both a stabilizing (no clusters) and destabilizing (cluster formation) effect, depending on the relative magnitudes of fiber flexibility and rotational diffusion.

Another interesting problem is the sedimentation of particle clouds. Experiments and simulations have used rigid rods and spheres to show that an initially spherical cloud evolves into a descending toroidal shape reminiscent of a high Reynolds number vortex ring. The torus eventually breaks into smaller clouds which evolve into tori (Machu et al. (2001); Metzger et al. (2007, 2005); Park et al. (2010)). Nazockdast et al. (2017b) used nonlocal SBT to study clouds of $N_{F}=1024$ flexible fibers in a semi-dilute regime. Their observations were in general agreement with the previous experimental and computational studies on rigid rods (Fig. 9c). They observed no qualitative effect of fiber flexibility on the collective behavior of the cloud though flexibility does result in changes in the shape of the torus ring (manuscript in preparation). Fig. 9c shows the tori for two clouds with fibers of differing effective rigidities (65 in ratio). In short, increasing flexibility results in a smaller tori.

VI.3 Collective dynamics of anchored fibers.

In many biological processes and industrial applications fiber assemblies are attached to substrates and particles. Detailed theoretical descriptions are still few, and a full accounting is beyond the scope of our review, but we discuss several interesting recent examples.

Alvarado et al. (2017) have recently shown how to construct soft flow rectifiers by mounting dense beds of tilted flexible fibers onto the walls of a flow channel (Fig. 9d, right panel). When fluid is driven in the direction of the fiber tilt, the fiber bed depresses and facilitates flow. When fluid is pumped in the opposite direction the bed rises and occludes flow. The mechanical responses were quantified in a Taylor-Couette device, and rationalized theoretically using a homogenous Brinkman approximation for fluid flow through the bed.

In a related, but very different, work Nazockdast et al. (2017b) studied the rheology of centrosomal arrays of flexible microtubules in eukaryotic cells. Using nonlocal SBT, they simulated the dynamics of flexible fibers radially anchored to a sphere moving under an external oscillatory force, $F=F_{0}\cos\omega t$ (schematic inset in Fig. 9d, top left). At low frequencies the response was purely viscous, and the fiber array was well-modeled as a homogeneous Brinkman medium. At higher frequencies, the force and velocity response moved out of phase, reflecting viscoelasticity, and they used the phase lag to approximate the elastic, $G^{\prime}$, and loss moduli, $G^{"}$ (Fig. 9d). For a single fiber $G^{"}/G^{\prime}$ has a minimum at: $\omega_{1}=\tau_{\text{curv}}^{-1}$. A minimum was also observed in the many-fiber simulation but at $\omega\approx 25\tau_{\text{curv}}^{-1}$. In other words, HIs between the fibers decreased the relaxation time by a factor of $\approx 25$. An explanation is that the fibers create a poroelastic medium whose fluid penetration length decreases with increasing fiber number. Thus, bending modes with wavelengths above that of the penetration length are prohibited, resulting in shorter effective fiber length, and a much shorter relaxation time.

Nazockdast et al. (2017a) used the same computational framework to study pronuclear positioning, an important transport process in eukaryotic cells. They showed that ignoring HIs leads to order-of-magnitude mispredictions of the necessary positioning forces. They also showed that different proposed models for the active forces gave rise to different cellular flows, and proposed flow measurement as a tool to differentiate between active mechanisms.

Beds of actively beating cilia are common in biology, propelling ciliated organisms, pumping fluid in the brain (Faubel et al. (2016)), performing selective filtering (Nawroth et al. (2017)), and clearing mucus from the lungs. Mitran (2007) carried out a 3D simulation of propagating metachronal waves in rows of pulmonary cilia. They considered a system where cilia sit primarily in a Newtonian viscous fluid, while the tip moves within a viscoelastic fluid modeling a mucus layer. They found that minimizing the work done by molecular motors results in synchronized beating of cilia, and suggested that HIs between cilia can lead to synchronized beating. Ding et al. (2014) use a regularized Stokeslet method to study the dynamics of an infinite bed of driven cilia. They showed that metachronal waves not only enhance fluid transport above the bed, but increases mixing of fluid within it. Finally, in experiments Coq et al. (2011) investigated the dynamics of extended arrays of artificial cilia driven by a precessing magnetic field. Whereas the dynamics of an isolated cilium was a rigid body rotation, collective driven beating results in a symmetry breaking of the precession patterns. The trajectories of the cilia are anisotropic and experience a significant structural evolution as the actuation frequency increases (Fig. 9d, bottom, left).

VII SUMMARY AND OUTLOOK

As should be clear, flexible fibers interacting with flowing liquids present a rich source of problems in fluid/structure interaction. They also present complicated phenomena, which require sophisticated experimental techniques to observe and measure, and complicated theories through which to understand them.

There are many areas of open inquiry, among them the interactions of fibers with complex media. Viscoelastic responses are typical of many biological environments, such as the reproductive tract Fauci & Dillon (2006) or inside of the cell Wirtz (2009). There is a developing literature on swimming of microorganisms in viscoelastic fluids (Elfring & Lauga (2015); Sznitman & Arratia (2015)). Using a 2D slender-body actuated elastica model, Thomases & Guy (2017) recently studied the role of body flexibility for undulatory swimming in viscoelastic fluids. There have been few if any studies of flexible fibers interacting with complex flows of complex fluids; Recent simulations of Yang & Fauci (2017) (Fig. 5f) of fiber transport by cellular viscoelastic flow is the first of which we are aware. The fundamental theoretical difficulty is the necessity of evolving bulk elastic stresses via transport nonlinearities. This makes fluid-structure problems, much less those with multiple elastic bodies, very challenging for viscoelastic flows and morally equivalent to those for the Navier-Stokes equations. On the experimental side, it remains challenging to synthesize complex fluids with well-characterized (and simple!) rheological responses.

New kinds of mathematical coarse-grained descriptions need to be developed to describe collective behavior of flexible fibers in fluids, especially when hydrodynamic interactions are strong. One regime where progress is being made is when the fibers can be considered as well-aligned. In recent work, Stein & Shelley (in preparation) have developed a continuum Brinkman-type model that captures the anisotropic drag from elongated flexible structures to compute the flow feedback to the bending and tensile response of a porous elastic medium (Moeendarbary et al. (2013); Strychalski et al. (2015)). Figure 9d (top right) shows the model’s result in simulating a soft flow rectifiern by the bending of a bed of tilted and anchored elastic fibers (Fig. 9d (top middle) from Alvarado et al. (2017)).

ACKNOWLEDGMENTS

E.N. and M.J.S. acknowledge the support of the National Science Foundation, the National Institutes of Health, and the Department of Energy. A.L. acknowledges support from European Research Council Consolidator Grant 682367 under the European Union’s Horizon 2020 program.

LITERATURE CITED