Coarse graining the dynamics of immersed and driven fiber assemblies

An important class of fluid-structure problems involves the dynamics of ordered arrays of immersed, flexible fibers. While specialized numerical methods have been developed to study fluid-fiber systems, they become infeasible when there are many, rather than a few, fibers present, and these methods do not lend themselves to analytical calculation. Here, we introduce a coarse-grained continuum model, based on local-slender-body theory, for elastic fibers immersed in a viscous Newtonian fluid. It takes the form of an anisotropic Brinkman equation whose skeletal drag is coupled to elastic forces. This model has two significant benefits: (1) the density effects of the fibers in a suspension become analytically manifest and (2) it allows for the rapid simulation of dense suspensions of fibers in regimes inaccessible to standard methods. As a first validation, without fitting parameters, we achieve very reasonable agreement with three-dimensional (3D) immersed boundary simulations of a bed of anchored fibers bent by a shear flow. Second, we characterize the effect of density on the relaxation time of fiber beds under oscillatory shear and find close agreement to results from full numerical simulations. We then study buckling instabilities in beds of fibers, using our model both numerically and analytically to understand the role of fiber density and the structure of buckling transitions. We next apply our model to study the flow-induced bending of inclined fibers in a channel, as has been recently studied as a flow rectifier, examining the nature of the internal flows within the bed, and the emergence of inhomogeneous permeability. Finally, we extend the method to study a simple model of metachronal waves on beds of actuated fibers, as a model for ciliary beds. Our simulations reproduce qualitatively the pumping action of coordinated waves of compression through the bed.

Figure 1

Many biological and industrial problems involve a regular bed of fibers that are clamped to a boundary, as shown in panel (a) (which depicts the boundary to which the fibers are fixed as flat, though it need not be). When fiber orientation and motion are coherent, it is viable to coarse grain the fluid-fiber system. The effect of the fibers on the fluid is coarse grained by averaging the forces exerted by the fibers on the fluid in a disk of radius $r$ and thickness $d$ that lies in the plane normal to the unit tangent $\mathbf{n}$ at $\mathbf{x}$, denoted by ${\mathcal{D}}_r\left(\mathbf{x}\right)$, shown as the gray disk in panel (b).

Figure 2

Panel (a) shows a dense fiber bed, with all physical fibers represented. Because the bed is extremely dense, a direct discretization (e.g., via slender-body theory or the immersed boundary method) would be computationally expensive. However, the orientation field of the fibers is relatively smooth and can be accurately approximated using the Brinkman-elasticae model given in Eq. (26) and a sparser discretization, such as the one shown in panel (b).

Figure 3

Panel (a) shows a possible configuration of fibers. The smooth force density $\mathcal{F}$ is easily computable on the discretized fibers, shown in purple. In order to facilitate rapid solution of the Stokes equation, these forces are interpolated to the nodes of a regular grid. Panel (b) shows an enlarged view of the gray rectangle in panel (a). Forces are known at the discretization points of the fiber field, shown as purple squares. Nodes of the regular grid are shown as circles. A line separating the fiber region from the purely fluid region is formed by connecting the tips of the fibers. Points are determined to be inside of the fiber region if they fall within a polygon enclosing the fiber region (in this geometry, simply under the dashed line). Points within the fiber region are shown in black. Forces are interpolated from the known values at the purple squares to the black circles.

Figure 4

Panel (a) shows the shear-flow setup used for immersed boundary simulations when ${\rho }_0=0.25$. The initial and final positions of the fiber are shown in gray and black, respectively, and the gray curve shows the path of the fiber tip. The three contour plots show the $x$ component of the velocity at $z=Z\left(L\right)/3, z=2Z\left(L\right)/3$, and $z=Z\left(L\right)$. Panels (b)–(d) show a comparison of results for immersed boundary simulations (markers) and the Brinkman-elasticae model (lines), for a bed of fibers subjected to a uniform shear, across a range of densities. Panel (b) shows the total fiber deflection $\left[X\right(L)-X(0\left)\right]$ at steady state, panel (c) shows the time at which the fiber was deflected by $95\%$ of its steady state deflection, and panel (d) shows the flow occlusion at steady state (i.e., the ratio of flow through the channel with the deformed fibers present to the flow through the channel with no fibers present).

Figure 5

Panel (a) shows the deflection of a fiber bed subject to a shear flow across a range of effective densities $\xi$. The blue line denotes $\xi =1$, where the bed transitions from a regime with behavior dominated by single-fiber dynamics to a regime where the behavior is dominated by many fiber dynamics. The value $\xi =1$ at which this transition occurs corresponds to when the average interfiber spacing is approximately half of the fiber length. Panel (b) shows the maximum number of Newton iterations required per time step and the maximum number of GMRES iterations required to invert the Jacobian in each Newton iteration.

Figure 6

Oscillatory rheology of a fiber bed. Panel (a) shows the storage ($G^{\prime }$) and loss ($G^{\prime \prime }$) moduli for the fiber bed, computed at $\xi =10$. The frequency at which $G^{\prime }$ and $G^{\prime \prime }$ intersect, ${\omega }_0^{*}$, is shown with a dashed gray line; this frequency gives a relaxation timescale $t^{*}=2\pi /{\omega }_0^{*}$ that depends on $\xi$. Panel (b) shows how the relaxation time $t^{*}$ depends on the effective density $\xi$. The $\xi =0$ limit, computed analytically, is shown as the horizontal gray line, and is in close agreement with the numerically computed values for low $\xi$. The blue line denotes $\xi =1$ and separates the single-fiber regime, which has a timescale independent of $\xi$, and the many-fiber regime, which has a timescale that scales with the inverse of $\xi$.

Figure 7

Stability of fiber beds under gravitational load. Panel (a) shows the growth rate of a small perturbation $f$, as a function of the effective buckling force $\tilde{g}$ and the effective density $\xi$. Panel (b) shows the growth rate computed from both the eigendecomposition (lines) and numerical simulation of the full nonlinear system (markers), as a function of the effective density $\xi$, for several values of the effective gravitational force $\tilde{g}$.

Figure 8

Physical setup for fiber-bed flow rectifier. Fibers of length $L$ are placed at an angle $\theta$ relative to the base of a channel of height $H$. Flow in different directions induces an asymmetric response in the fiber bed, allowing more fluid to pass in the forward direction than the reverse direction (see Fig. 9).

Figure 9

Panels (a) and (b) show the fiber deformation and flow fields for a flow rectifier with $\theta ={45}^{\circ }, \tilde{L}=1, \xi =100$, and $\tilde{E}=0.01$. Blue arrows show the flow field and the grayscale shows the fiber density, scaled consistently. Panel (c) shows the velocity field as a function of $z$ corresponding to panel (a), in black, and panel (b), in blue; $u$ for the “backward” direction has been multiplied by $-1$ to make comparison easier.

Figure 10

Panel (a) shows the impedance ratio (${\mathcal{I}}_r$), and the forward impedance (${\mathcal{I}}_f$), for $\theta ={45}^{\circ }$ and $\tilde{L}=1$, across a wide range of $\xi$ with $\xi \tilde{E}=1$. Panels (b) and (c) show the velocity field and fiber deformation, respectively, with solid lines denoting results in the forward direction and dashed lines denoting results in the backward direction. Colors denote density: Purple is the least dense ($\xi =10$); red is the most dense ($\xi =2000$). The colors are equispaced on a log scale in $\xi$; the two most purple lines are much closer together in $\xi$ than the two most red lines.

Figure 11

The rectification ratio ${\mathcal{I}}_r$, as a function of the mount angle $\theta$ and $\lambda =\xi \tilde{E}$ with $\xi =1000$. The white line shows the value of $\lambda$ that produces the largest rectification ratio, as a function of $\theta$.

Figure 12

Numerical performance and accuracy in two dimensions. Panel (a) shows the wall clock time per time step, in seconds, measured as the average over the first 10 time steps after $t=1$, for a range of densities and discretizations. Panels (b) through (d) show fiber deformation (in white) and velocity magnitude (color scale, larger magnitudes are darker) at time $t=5$, for coarser (b) to finer (d) discretizations. All discrete fibers used in the computation are shown.

Figure 13

Panel (a) shows the transient behavior and decay of oscillations in the net flux due to small amplitude ($\gamma =0.01$) actuation. The density and effective rigidity for the main figure are $\xi =100$ and $\tilde{E}={10}^{-4}$, respectively, and thus $\lambda =\xi \tilde{E}=0.01$. Transient behavior relaxes with a timescale ${\lambda }^{-1}=100$, while a persistent oscillation at the driving frequency damps with the timescale $2\pi {\lambda }^{-1}$. The inset shows the net flux at large amplitude ($\gamma =\pi /4$) and in the pulsatile regime ($\lambda =2$). In this case, the persistent oscillations are large and do not decay. Panel (b) shows the rate of decay of the persistent oscillations for large amplitude ($\gamma =0.5$) actuation with $\xi =1000$ over a range of $\tilde{E}$. Although both the decay rate and flux scale with $\lambda$ for large $\lambda$, a pulsatile regime exists in a range about $\lambda \approx 2$ in which these oscillations do not decay.