Efficient Implementation of Elastohydrodynamics via Integral Operators

TL;DR

This paper introduces an efficient integral operator-based method for simulating the complex dynamics and interactions of flexible filaments in fluid flows, addressing computational challenges and enabling new biological and physical insights.

Contribution

The manuscript develops a computationally efficient integral equation framework using regularized stokeslets for simulating multiple interacting filaments, including active and passive cases.

Findings

Filament interactions can promote buckling instabilities.

The method accurately models filament dynamics in various flow conditions.

Open-source MATLAB implementation facilitates further research.

Abstract

The dynamics of geometrically non-linear flexible filaments play an important role in a host of biological processes, from flagella-driven cell transport to the polymeric structure of complex fluids. Such problems have historically been computationally expensive due to numerical stiffness associated with the inextensibility constraint, as well as the often non-trivial boundary conditions on the governing high-order PDEs. Formulating the problem for the evolving shape of a filament via an integral equation in the tangent angle has recently been found to greatly alleviate this numerical stiffness. The contribution of the present manuscript is to enable the simulation of non-local interactions of multiple filaments in a computationally efficient manner using the method of regularized stokeslets within this framework. The proposed method is benchmarked against a non-local bead and link model, and recent code utilizing a local drag velocity law. Systems of multiple filaments (1) in a background fluid flow, (2) under a constant body force, and (3) undergoing active self-motility are modeled efficiently. Buckling instabilities are analyzed by examining the evolving filament curvature, as well as by coarse-graining the body frame tangent angles using a Chebyshev approximation for various choices of the relevant non-dimensional parameters. From these experiments, insight is gained into how filament-filament interactions can promote buckling, and further reveal the complex fluid dynamics resulting from arrays of these interacting fibers. By examining active moment-driven filaments, we investigate the speed of worm- and sperm-like swimmers for different governing parameters. The MATLAB(R) implementation is made available as an open-source library, enabling flexible extension for alternate discretizations and different surrounding flows.

Equations

$$\mathbf{X}(s,t)={\mathbf{X}}_0(t)+{\int }_0^s[\cos \theta (s^{\prime },t),\sin \theta (s^{\prime },t),0{]}^{T}d{s}^{\prime },$$

$$0=-\nabla p+{\nabla }^2\mathbf{u}+\mathbf{F},\nabla \cdot \mathbf{u}=0,$$

$$u_j(\mathbf{x},t)=\frac{1}{8\pi }{\int }_0^1{S}_{jk}^{\epsilon }\left(\mathbf{x},\mathbf{X}(s^{\prime },t)\right)f_k(s^{\prime },t)d{s}^{\prime }+\mathcal{O}\left({\epsilon }^r\right),$$

$$S_{jk}^{\epsilon }\left(\mathbf{x},\mathbf{X}\right)={\delta }_{jk}\frac{|\mathbf{x}-\mathbf{X}|+2{\epsilon }^2}{{\left(|\mathbf{x}-\mathbf{X}{|}^2+{\epsilon }^2\right)}^{3/2}}+\frac{\left(x_j-X_j\right)\left(x_k-X_k\right)}{{\left(|\mathbf{x}-\mathbf{X}{|}^2+{\epsilon }^2\right)}^{3/2}}.$$

$${\partial }_t{X}_j(s,t)=\frac{1}{8\pi }{\int }_0^1{S}_{jk}^{\epsilon }\left(\mathbf{X}(s,t),\mathbf{X}(s^{\prime },t)\right)f_k(s^{\prime },t)d{s}^{\prime },0⩽s⩽1,t⩾0.$$

$${\int }_0^1-\mathbf{f}(s,t)ds$$

$${\int }_0^1\mathbf{X}(s,t)\times \left(-\mathbf{f}(s,t)\right)ds$$

$$M(s,t)=EI{\partial }_s\theta (s,t),$$

$${\partial }_s\mathbf{N}(s,t)=\mathbf{f}(s,t).$$

$$0={\partial }_{s}M(s,t)+{\mathbf{e}}_3\cdot {\partial }_s\mathbf{X}(s,t)\times \mathbf{N}(s,t).$$

$$0=-M(s,t)+\bm{e}_{3}\cdot\left[\big{(}\bm{X}(s^{\prime},t)-\bm{X}(s,t)\big{)}\Big{|}_{s}^{1}-\int_{s}^{1}(\bm{X}(s^{\prime},t)-\bm{X}(s,t))\times\bm{f}(s^{\prime},t)\,ds^{\prime}\right].$$

$${\partial }_s\theta (s,t)=-{\mathbf{e}}_3\cdot {\int }_s^1(\mathbf{X}(s^{\prime },t)-\mathbf{X}(s,t))\times \mathbf{f}(s^{\prime },t)d{s}^{\prime },0⩽s⩽1,t⩾0.$$

$${\partial }_t{X}_j(s,t)$$

$$0$$

$${\partial }_s\theta (s,t)$$

$$\text{mbox}where\mathbf{X}(s,t)$$

$$\text{mbox}and{\partial }_t\mathbf{X}(s,t)$$

$${\mathbf{X}}^{[n]}(t):=\mathbf{X}\left((n-1)\Delta s,t\right),n=1,\dots ,Q+1.$$

$${\mathbf{X}}^{[m+1]}(t)={\mathbf{X}}^{[1]}(t)+\sum _{n=1}^m\Delta s{\left[\cos {\tilde{\theta }}^{[n]}(t),\sin {\tilde{\theta }}^{[n]}(t),0\right]}^T,m=1,\dots ,Q,$$

$${\tilde{\mathbf{X}}}^{[m]}(t)={\mathbf{X}}^{[1]}(t)+\sum _{n=1}^{m-1}\Delta s{\left[\cos {\tilde{\theta }}^{[n]}(t),\sin {\tilde{\theta }}^{[n]}(t),0\right]}^T+\frac{\Delta s}{2}{\left[\cos {\tilde{\theta }}^{[m]}(t),\sin {\tilde{\theta }}^{[m]}(t),0\right]}^T,$$

$$\mathbf{f}(s,t)\approx {\tilde{\mathbf{f}}}^{[m]}(t)=\mathbf{f}({\tilde{s}}^{[m]},t),(m-1)\Delta s⩽s<m\Delta s,$$

$$\bm{X}(s,t)\approx\tilde{\bm{X}}^{[m]}(t)+\left(s-\tilde{s}^{\,[m]}\right)\Big{[}\cos\theta\left(\tilde{s}^{\,[m]},t\right),\,\sin\theta\left(\tilde{s}^{\,[m]},t\right),\,0\Big{]}^{T},$$

$${\partial }_t{\tilde{X}}_j^{[m]}(t)=\frac{1}{8\pi }\sum _{n=1}^Q{\tilde{f}}_k^{[n]}(t)$$

$$\cdot {\int }_{(n-1)\Delta s}^{n\Delta s}{S}_{jk}^{\epsilon }\left({\tilde{\mathbf{X}}}^{[m]},{\tilde{\mathbf{X}}}^{[n]}(t)+\left(s^{\prime }-{\tilde{s}}^{[n]}\right){\left[\cos \theta ({\tilde{s}}^{[n]},t),\sin \theta ({\tilde{s}}^{[n]},t),0\right]}^T\right)d{s}^{\prime }.$$

$${\partial }_t{\tilde{X}}_j^{[m]}(t)=\frac{1}{8\pi }\sum _{n=1}^Q{I}_{jk}^{[m,n]}\left(t;\Delta s,\epsilon \right){\tilde{f}}_k^{[n]}(t),m=1,\dots ,Q,$$

$$0$$

$$\frac{{\tilde{\theta }}^{[m+1]}(t)-{\tilde{\theta }}^{[m]}(t)}{\Delta s}$$

$${\partial }_t{\tilde{X}}_j^{[m]}(t)$$

$$\text{mbox}where{\mathbf{X}}^{[m]}(t)$$

$${\tilde{\mathbf{X}}}^{[m]}(t)$$

$$\text{mbox}and{\partial }_t{\tilde{\mathbf{X}}}^{[m]}(t)$$

$$+\frac{\Delta s}{2}{\partial }_t{\tilde{\theta }}^{[m]}(t){\left[-\sin {\tilde{\theta }}^{[m]}(t),\cos {\tilde{\theta }}^{[m]}(t),0\right]}^T,$$

$$\dot{\mathbf{Y}}=\mathcal{F}(\mathbf{Y}),\mathbf{Y}(0)={\mathbf{Y}}_0,$$

$$\displaystyle\textup{A}\left[\begin{array}[]{c}\dot{\bm{Y}}\\ \tilde{\bm{f}}\end{array}\right]$$

Full text

I Introduction

Flexible filaments are ubiquitous in the natural world, and thus a clear understanding of their behavior is paramount in many biological problems. Models for simulating the dynamics of these filaments, while plentiful, have historically been mathematically complex and numerically expensive, with even simple computational experiments taking in the order of hours or even days to solve (see Moreau et al. for detailed benchmarking moreau2018asymptotic ).

Micro-scale filament problems have been previously tackled using other modeling approaches which can be broadly separated into those based upon (a) a nonlinear PDE in the filament position (such as the method by Schoeller et al. schoeller2019method ), or (b) a discretization into simpler elements, such as beads with connecting springs jayaraman2012autonomous or interlocking gears delmotte2015general . For category (a), inextensibility is enforced using Lagrange multipliers of tension, which are often costly to compute. For the discrete approaches in category (b), other ways of enforcing this condition are used. For example, the bead model of Jayaraman et al. jayaraman2012autonomous prescribes large spring constants between each bead, contributing to the numerical stiffness of the system. Equivalently, the gears model of Delmotte et al. delmotte2015general imposes a non-holonomic constraint to ensure non-penetrability between adjacent beads, but as a result requires large numbers of points to represent a single filament.

A recent promising development via Moreau et al. moreau2018asymptotic , referred to as coarse graining, is based on reformulating the problem via an integral equation with the filament tangent angle as the dependent variable. The method, initially developed using a local hydrodynamic drag law, provides an efficient framework for simulating non-interacting filament dynamics. This approach builds upon the early studies of Brokaw brokaw1971bend ; brokaw1972computer and Hines & Blum hines1978bend , and contrasts with Cartesian formulations gadelha2010nonlinear ; tornberg2004simulating .

The contribution of the present manuscript is to enable efficient and accurate simulation of multiple, non-locally interacting, passive and active filaments in ambient flows by incorporating recent developments in the regularized stokeslet method cortez2018regularized ; gallagher2018meshfree with the integral formulation in terms of the tangent angle of Moreau et al. moreau2018asymptotic .

The potential applications for a fast and accurate filament modeling framework are numerous. There has long been interest in understanding the mechanics and regulation of sperm flagellar movement, in particular problems relating to: understanding the mechanical structure and motor regulation brokaw1971bend ; lindemann1994model ; guo2018bistability ; oriola2017nonlinear , investigating the response of the flagellar beat to its rheological environment huang2018hydrodynamic ; gadelha2013counterbend ; smith2009bend , understanding the dynamics of sperm due to surrounding solid walls montenegro2015spermatozoa ; denissenko2012human , and studying the effect of viscosity on sperm swimming woolley2001study . For a detailed review surrounding the importance of the sperm flagellum see Gallagher et al. gallagher2018casa . Furthermore, such a method could be used to investigate phenomena associated with epithelial cilia driven flows such as: cilia waveform modulation by length brokaw1971bend , the effects of flow induced by cilia on embryonic development montenegro2012modelling , studying the physical limits of flow sensing ferreira2017physical , and investigating the mechanical structure of the axoneme in cilia omori2017nodal . Another avenue of active-filament research to which the proposed framework could be applied is magnetic swimmers gadelha2013optimal . These models have wider relevance in the field of synthetic biology, with particular application to microscopic bacteriophage-based fibre sensors pacheco2011detection ; lobo2015direct ; gallagher2017model and flexible filament microbots montenegro2018microtransformers . The proposed framework could be used to further investigate the dynamics of bundles of filaments coy2017counterbend , and additionally has applications in the multi-scale studies of complex polymeric fluids, and flagellar movement through them yang2017dynamics ; wrobel2016enhanced .

We will extend the framework introduced by Moreau et al. moreau2018asymptotic , augmenting and reformatting their formulation with the method of regularized stokeslets of Cortez et al. cortez2001method ; cortez2005method . These methods have proven to be accurate and effective in modeling the hydrodynamics in various multiple-fiber scenarios stein2019coarse ; olson2015hydrodynamic . The method of regularized stokeslets enables the modeling of non-local hydrodynamics within and between filaments, and between filaments and surrounding structures. The method is implemented in a numerically efficient manner, retaining the computational economy and low hardware requirements inherited from the Moreau et al. formulation.

The structure of this manuscript is as follows: in Sec. II, the Elastohydrodynamic Integral Formulation (EIF) for a single filament is proposed. In Sec. III, alterations to the EIF for various single- and multi-filament scenarios are presented. Verification and benchmarking of the method is given in Sec. IV. Simulation results of the problems formulated in Sec. III are then presented in Sec. V, followed by discussion of the results and of further possible applications in Sec. VI. The MATLAB® code for the methods described within this report are provided in the associated GitLab repository, available at gitlab.com/atticushm/eif.

II Model formulation

The dynamics of elastic filaments in Stokes flow will be modeled by constructing integral operator formulations of the governing fluid and elastodynamic equations. Each filament is described by the time-evolving tangent angle $\theta(s,t)$ along its arclength. Taking $\bm{X}_{0}\left(t\right)$ to be the leading point at time $t$, the filament geometry is then given by

[TABLE]

where $s\in[0,1]$ and $\theta\in[-\pi,\pi)$ are the dimensionless arclength (scaled by the filament length $L$) and tangent angle respectively. The velocity of a point on the filament is given by differentiating Eq. (1) with respect to dimensionless time, scaled by $\tau$. The behavior of planar filaments in a Newtonian fluid is considered, which can be described by the three-dimensional dimensionless Stokes flow equations,

[TABLE]

where $\bm{u}(\bm{x},t)$ is fluid velocity, $p(\bm{x},t)$ is pressure, and $\bm{F}(\bm{x},t)$ is a force exerted by the body onto the fluid. As shown by Cortez cortez2005method , an approximate solution to Eq. (2) is given by the regularized stokeslet integral,

[TABLE]

where $\bm{f}(s,t)$ is the force per unit length exerted by the filament on the fluid, non-dimensionalized with the scaling $\mu L/\tau$ for a given fluid dynamic viscosity $\mu$, and $\bm{X}(s,t)$ denotes the filament position as a function of arclength and time. The error arises from the regularization of the stokeslet, for a chosen regularization parameter $0<\varepsilon\ll 1$, where $r=1$ or $2$ in the near- and far-field respectively. The combined process of solving Eq. (3) for the unknown force densities $\bm{f}(s,t)$ will result in errors of order $\varepsilon$ everywhere martin2019use . Summation convention dictates that repeated indices are summed over and unrepeated indices range over $\{1,2,3\}$. The kernel of the integral in Eq. (3) is known as the regularized stokeslet, defined

[TABLE]

Applying the no-slip condition $\bm{u}(\bm{X}(s,t),t)=\partial_{t}\,\bm{X}(s,t)$ on the boundary of the filament yields the regularized stokeslet integral equation

[TABLE]

Inertialess dynamics requires that the filament is force and moment free at each instant, giving

[TABLE]

The integral formulation for the hydrodynamic problem thus comprises Eqs. (5), (6) and (7). Next, the elastic behavior of the filaments is considered.

The integral operator for the elastodynamic behavior is formulated by considering the force and moment balance over an infinitesimal segment of a filament. We denote by $\bm{N}(s,t)$ and $M(s,t)$ the contact force and contact moment respectively exerted by a distal segment of filament $[s,1]$ on a proximal segment $[0,s)$. Constitutively linear elasticity and bending in the $xy$-plane implies that the dimensional contact moment relates to the curvature via

[TABLE]

where $EI$ is the bending modulus of the filament. The fluid dynamic force density $\bm{f}(s,t)$ and $\bm{N}(s,t)$ are related by

[TABLE]

Under the assumption of free boundary conditions, the contact force and moment are zero at each end of a filament. Moment balance over the infinitesimal segment reveals that the contact force and moment are related as

[TABLE]

Integration of Eq. (10) over a distal segment $[s,1]$ yields, on application of the moment-free boundary condition and Eq. (9),

[TABLE]

Application of the the force-free condition $\bm{N}(1,t)=0$, the constitutive Eq. (8), and the choice of time scale $\tau=\mu L^{4}/EI$, gives the non-dimensionalized integral formulation for the filament tangent angle evolution

[TABLE]

Combining Eqs. (1), (5), (6), (7), and (12), we obtain the elastohydrodynamic integral formulation (EIF) for the evolution of a filament in Stokes flow

[TABLE]

with unknowns $\bm{X}_{0}(t)$, $\theta(s,t)$ and $\bm{f}(s,t)$. In solving the EIF for a particular physical problem, we must specify the initial position $\bm{X}_{0}(0)$ and tangent angle $\theta(s,0)$.

II.1 Spatial discretization of the EIF

To solve the EIF, we spatially discretize filaments to obtain a set of integro-differential equations which can be numerically evaluated. Dividing the arclength domain into $Q$ segments of equal length $\Delta s=1/Q$, the positions of the resulting segment endpoints are denoted as

[TABLE]

The angle connecting $\bm{X}^{[n]}(t)$ to $\bm{X}^{[n+1]}(t)$ approximating \theta\big{(}(n-1)\Delta s,t\big{)} is denoted $\tilde{\theta}^{\,[n]}(t)$, for $n=1,\ldots,Q$. The positions of the endpoints are given in terms of the initial point $\bm{X}^{[1]}(t)$ and discretized tangent angle $\tilde{\theta}^{\,[n]}(t)$ as

[TABLE]

with the segment midpoints

[TABLE]

for each $m=1,\dots,Q$. An illustration of this discretization is displayed in Fig. 1. Differentiating Eq. (20) with respect to time yields the kinematic equation for the segment velocities.

For the fluid dynamics, rather than using a conventional (and potentially expensive), quadrature rule for evaluating the rapidly-varying kernel $S_{jk}^{\varepsilon}$, we employ the method of Smith smith2009boundary . By approximating the force density in Eq. (13) as piecewise constant along each segment, the kernel can be analytically integrated, reducing the level of quadrature needed to evaluate the slowly-varying force density (for higher-order force discretizations see Cortez cortez2018regularized ). Writing

[TABLE]

where $\tilde{s}^{\,[m]}$ denotes the arclength at the midpoint of the $m$th segment, with a piecewise linear discretization for the filament,

[TABLE]

we obtain the spatially discrete equation,

[TABLE]

The integral in Eq. (24) can be calculated exactly by transforming to a local coordinate system in which one axis is aligned with the segment tangent (details are given in App. B, Eqs. B1–B3 of smith2009boundary , although note that the simplified form for the near-field integral in Eq. B4 contains a typographical error in the $\delta s/\varepsilon$ fraction, which is upside-down). For brevity we denote the integral in Eq. (24) as $I_{jk}^{\,[m,n]}(t;\Delta s,\varepsilon)$, yielding the system of linear equations,

[TABLE]

with the force and moment balance equations given in Eqs. (6) and (7) discretized via the midpoint rule as

[TABLE]

The semi-discrete form of the EIF is then

[TABLE]

where $m=1,\dots,Q-1$ in Eq. (30) and $m=1,\dots,Q$ in Eq. (31) and (33).

By using Eqs. (32) and (33), the variables $\tilde{\bm{X}}^{[1]}(t),\ldots,\tilde{\bm{X}}^{[Q+1]}(t)$ and $\bm{X}^{[2]}(t),\ldots,\bm{X}^{[Q]}(t)$ can be eliminated from Eqs. (29), (30) and (31). The resulting system is then linear in $\partial_{t}\,\bm{X}^{[1]}(t),\,\partial_{t}\,\tilde{\theta}^{\,[1]},\ \dots,\,\partial_{t}\,\tilde{\theta}^{\,[Q]},\,\tilde{\bm{f}}^{\,[1]}(t),\dots,\,\tilde{\bm{f}}^{\,[Q]}(t)$. Hence, given the discrete configuration $\tilde{\bm{X}}^{[1]}(t),\,\tilde{\theta}^{\,[1]}(t),\dots,\,\tilde{\theta}^{\,[Q]}(t)$, the rate of change of position and angle can be found by solving a dense $(3Q+2)\times(3Q+2)$ system of linear equations. Thus, the semi-discrete system can be expressed concisely as an autonomous non-linear initial value problem,

[TABLE]

where $\bm{Y}(t):=[\tilde{\bm{X}}^{[1]}(t),\ \tilde{\theta}^{\,[1]}(t),\dots,\tilde{\theta}^{\,[Q]}(t)]$ is a $Q+2$ column vector and $\dot{\bm{Y}}(t)=\partial_{t}\,\bm{Y}(t)$. By augmenting the problem with the unknown force densities, this can be written as the matrix system,

[TABLE]

where A is a $(3Q+2)\times(3Q+2)$ block matrix, $\bm{b}$ is a $3Q+2$ column vector, and $\tilde{\bm{f}}$ is a $2Q$ column vector, so that the concatenation $[\,\dot{\bm{Y}},\tilde{\bm{f}}\,]^{T}$ is a $3Q+2$ column vector. The matrix blocks of A $\left(\textup{A}_{E},\ \textup{A}_{K},\text{ and }\textup{A}_{H}\right)$ encode the elastodynamic, kinematic, and hydrodynamic equations given by Eqs. (30), (33) and (31) respectively. In the vector $\bm{b}$, the first entry corresponds to the moment balance on the whole filament (Eq. (29)), the subsequent $Q-1$ rows are the moment balances about each interior joint (Eq.(30)), and the next 2 rows correspond to the total force balance (Eq. (28)). The remaining zero entries correspond to the equivalence between the hydrodynamic and kinematic velocities (Eqs. (31) and (32)). The matrix system given in Eq. (37) is solved for $\dot{\bm{Y}}$ and $\tilde{\bm{f}}$ at each time step using the MATLAB® backslash command, and the resulting rates vector $\dot{\bm{Y}}$ is integrated using the built-in variable-step, variable-order ODE solver ode15s shampine1997matlab . To demonstrate the ease of application of the EIF framework, all simulations are performed using the MATLAB® R2019a default settings for ode15s. In particular, the absolute and relative error tolerances are $10^{-6}$ and $10^{-3}$ respectively. While this IVP exhibits some stiffness, it is less stiff than the systems produced by other methods, as the integral formulation avoid the need of additional Lagrange multipliers to ensure filament inextensibility. At each time step, the filament is constructed according to Eq. (32) with $\Delta s=1/Q$, ensuring that filament arclength is preserved over the course of the simulation.

III Single and multi-filament problems

In the following section, we apply the EIF framework described in Sec. II to problems involving single or multiple filaments in the presence of body forces, surrounding flow, or undergoing self-powered propulsion.

III.1 A single passive filament in shear flow

We investigate the dynamics of a single passive filament in a linear shear flow, ${\bm{u}^{*}_{s}(t)=\dot{\gamma}\,x_{2}^{*}(s^{*},t^{*})\,\bm{e}_{1}}$, with shear rate $\dot{\gamma}$, where ∗ denotes a dimensionful variable. Non-dimensionalizing with respect to the length of the filament $L$, time scale $\tau=\dot{\gamma}^{-1}$, and force density scaling $\mu L\dot{\gamma}$ (where $\mu$ is the dynamic viscosity of the surrounding fluid), yields the non-dimensionalized equation for the hydrodynamic velocity

[TABLE]

Additionally, from the dimensional version of Eq. (16), together with the scalings defined above, we obtain the non-dimensionalized elastohydrodynamic integral equation

[TABLE]

where the dimensionless viscous-elastic parameter

[TABLE]

quantifies the ratio of viscous to elastic forces on a shear timescale, where $EI$ is the bending rigidity of the filament. The apparent flexibility of a filament is completely characterized using $\mathcal{V}$, with large values describing floppy fibers, and small values stiff fibers. The similarity between Eqs. (16) and (44) results in a similarly discretized form as in Eq. (30).

The hydrodynamic shear flow equation (Eq. (43)) is semi-discretized following the steps in Sec. II to obtain

[TABLE]

The system of equations for a single passive filament in shear flow is thus given by Eq. (37) but with the alterations

[TABLE]

and where the vector of unknowns $[\,\dot{\bm{Y}},\tilde{\bm{f}}\,]^{T}$ remains unchanged from Eq. (37).

III.2 A single passive filament sedimenting under gravity

The EIF method also allows for implementation of a body force to the system. In this section, the simulation of passive filaments sedimenting under gravity is considered. Assuming uniform mass per unit length $\rho$, the force per unit length due to gravity acting upon the filament is $-\rho g\bm{e}_{2}$. Non-dimensionalizing with respect to the filament length $L$, time scaling $\tau=\mu L^{4}/EI$, and force density $\mu L/\tau$, addition of the gravitational force density produces the force and moment balance equations

[TABLE]

where $\bm{X}_{c}(t)$ is the center of gravity of the filament, $\bm{f}(s,t)$ is the force per unit length the filament exerts on the fluid, and $\mathcal{G}$ is the elasto-gravitational parameter

[TABLE]

Following a derivation similar to that presented in Sec. II, the elastohydrodynamic equation is found as

[TABLE]

which, along with Eqs. (51) and (52), has spatially discretized form

[TABLE]

As in Sec. II.1, we form a matrix system encoding Eqs. (55)–(57) along with Eqs. (31) and (33)

[TABLE]

which can be solved to find the filament velocities $\partial_{t}\,\bm{X}_{0}(t),\,\partial_{t}\,\tilde{\theta}^{\,[m]}(t)$ and force densities $\tilde{\bm{f}}^{\,[m]}(t)$ for $m=1,\dots,Q$. The matrix $\textup{A}_{g}$ has the same block form as A, but with an alteration in the first row of the elasticity block $\textup{A}_{E}$ due to do the inclusion of the center of gravity in the total moment balance equation (Eq. 56). The right hand side vector is constructed as $\bm{b}_{g}=\bm{b}+\hat{\bm{b}}$, where $\hat{\bm{b}}$ encodes all the changes required to the right hand side vector $\bm{b}$ (given in Eq. (41)) after expansion and rearrangement of Eqs. (55), (56), and (57).

III.3 A single active filament

For the problems considered in Secs. III.1 and III.2, we can consider active filaments by including a time-dependent moment density term added to the elastodynamic formulation given in Eq. (12). The moment balance in Eq. (10) is extended

[TABLE]

where $\mathfrak{m}(s,t)$ is the moment per unit length that drives the actuation of the filament. Continuing the non-dimensionalization as in Sec. II, with time scaling $\tau=\omega^{-1}$, yields the non-dimensionalized elastohydrodynamic equation for a single active filament

[TABLE]

where $\mathcal{S}$ is the dimensionless swimming number, defined

[TABLE]

with $\omega$ being the angular velocity of the swimming beat prescribed to the filament. Note that the swimming parameter is different from the commonly-used sperm number Sp, which has dependence on a chosen resistance coefficient. Following work by Gadêlha et al. gadelha2010nonlinear and Montenegro-Johnson et al. montenegro2015spermatozoa , the active moment $\mathfrak{m}(s,t)$ is set as a traveling wave with amplitude $\mathfrak{m}_{0}$, wave number $k$ and angular frequency $\omega$. The semi-discretized form of Eq. (60) is

[TABLE]

for $m=1,\dots,Q$. In the following section, we consider the how the presented framework can be extended to model systems of multiple filaments.

III.4 Systems of multiple filaments

The EIF can be used to simulate the dynamics of large groups of filaments, accounting for the non-local hydrodynamic interactions between them. For each of the problems presented in Secs. III.1, III.2, and III.3, the kinematic and elastodynamic equations apply to each filament in the system individually. The hydrodynamic equations are extended so that interactions between all filaments are considered. For a system of $N$ passively relaxing filaments, this equation reads

[TABLE]

where the superscript $\{\alpha\}$ in this continuous equation refers to the $\{\alpha\}$th filament in the system. Whilst three-dimensional hydrodynamic effects are computed, the kinematics and hence elasticity of the multiple filament problem remains two-dimensional, ensuring planarity. Modification of Eq. (63) to consider external flows or body forces follows from the single filament derivations presented in Secs. III.1, III.2, and III.3, and non-dimensionalizing yields the same governing dimensionless parameters $\mathcal{V}$, $\mathcal{G}$, and $\mathcal{S}$ respectively. Discretization of the resulting equations is performed in the same manner, producing an ${N(3Q+2)\times N(3Q+2)}$ linear system similar in structure to those presented in each of the single filament problem descriptions. Full details of the equations as well as the associated linear systems for each of the multiple-filament problems are given in the Supplemental Material Sec. SII.

IV Model verification

To verify the accuracy and efficacy of the EIF, we compare the computed dynamics of a single relaxing filament between the proposed method and a high resolution bead and link method (BLM) formulation. Based upon the work of Jayaraman et al. jayaraman2012autonomous , the BLM accounts for non-local hydrodynamic interactions, and when highly resolved, provides accurate solutions (details provided in the Supplemental Material Sec. SIII). For this reason, the BLM is used to both verify the proposed EIF method, and as a reliable point of comparison between other existing methods. We consider the case of a filament, bent into a parabola, with initial condition $\bm{Y}$ constructed by sampling symmetrically from the curve $y=0.5x^{2}$, and ensuring unit arclength. This filament is then allowed to relax with no external forcing. Motion of the filament in this scenario is due to the constitutive bending moments along the arclength, given in Eq. (30). This experiment is simulated using each of four methods:

EIF-RSM: the proposed EIF method, which uses regularized stokeslets to account for non-local hydrodynamic interactions, 2. 2.

EIF-RFT: the EIF method, with resistive force theory in place of the method of regularized stokeslets to model local hydrodynamic interactions (refer to the associated Supplemental Material Sec. SIV for full details), 3. 3.

MGG: the original angle formulation method by Moreau, Giraldi and Gadêlha moreau2018asymptotic , which uses resistive force theory to model local hydrodynamic interactions, 4. 4.

BLM: the bead and link method, which accounts for non-local hydrodynamic interactions.

We include the EIF-RFT method to verify the equivalence of the present implementation with that of Moreau et al. moreau2018asymptotic under the reduction to local hydrodynamics. In all simulations in this paper the regularization parameter for use in the method of regularized stokeslets is chosen as $\varepsilon=0.01$.

The geometric configuration of the relaxing filament is simulated between $t=0$ and $t=0.02$, using both a high-resolution ($Q=200$) BLM formulation and the EIF-RSM method with $Q=100$. The initial and final filament shapes from each simulation are displayed in Fig. 2a, which show excellent agreement between these formulations. Quantitative comparison is shown in Fig. 2b, where we plot the root mean squared error (RMSE) between the position of the filament described by each model against that of the high-resolution BLM. Here excellent convergence is evident, with the RMSE being minimal for even small values of $Q$. Comparisons with MGG and EIF-RFT are also presented in this figure. While it might be counter-intuitive that the local hydrodynamic models initially increase in error with $Q$, this is due to drift in the center of mass that the non-local methods correctly capture. Increasing $Q$ in the local methods leads to increasingly resolved shape of the relaxing filament, leading to the convergence of this error. These comparisons highlight the change in dynamics when considering the inclusion of non-local hydrodynamics for even a simple single filament problem.

Moreau et al. demonstrated the significant reduction in computational cost achieved when formulating elastohydrodynamic problems as an integro-differential equation. With the added complexity of including non-local hydrodynamic modeling, the EIF-RSM still performs very well. In Fig. 2c we plot logarithmic comparisons of the simulation runtime for each of EIF-RSM, MGG, and BLM formulations. The walltime recorded is the total computational time for the method including setup time. For the tangent angle formulation methods, the majority of this time is accounted for by the linear solve at each time step.

It it unsurprising that the local hydrodynamic formulation (MGG) outperforms the methods with non-local interactions due to reduced complexity of the problem. The EIF-RSM and BLM method perform on par for increasing $Q$. However, the BLM method requires large numbers of beads in order to accurately capture the correct filament dynamics (see Sec. SIII of the provided Supplemental Material), whereas the EIF-RSM can achieve similar accuracy with fewer than half of the degrees of freedom required by the BLM (see Fig. 2b).

V Results

In the following section, results from the numerical experiments outlined in Sec. III are presented. We begin by examining the dynamics of single passive filaments in shear flow and sedimenting under gravity, and active filaments with prescribed internal moments. Additionally, we investigate larger arrays of passive filaments, again in shear flow and undergoing sedimentation, and systems of multiple swimming filaments. The regularization parameter is $\varepsilon=0.01$ and, unless otherwise stated, $Q=40$, with segment lengths $\Delta s=1/Q$. This choice of $Q$ is motivated by the convergence results given in the Supplemental Material Sec. SI. Simulations are run on a computer equipped with an Intel i7-8750H processor and 16GB of RAM.

V.1 Results: a single passive filament in shear flow

It is well known that for critical values of the characteristic parameter, a filament in shear flow exhibits shape buckling liu2018morphological ; tornberg2004simulating due to a stress difference across the fiber while under compression by the fluid. Dynamics of a filament in shear flow are simulated by solving Eq. (37) with A and $\bm{b}$ given by Eqs. (49) and (50), and characterized by parameter $\mathcal{V}$ (Eq. (45)).

We initialize our filament as a straight rod with a small perturbation following the method of Young young2009hydrodynamic , writing the initial angular configuration as

[TABLE]

for initial angle $\theta_{0}$ and small perturbation parameter $\Delta\theta_{0}$. As discussed by Tornberg & Shelley tornberg2004simulating , prescribing a small perturbation to a straight filament can drastically change the dynamics from rotational Jeffery orbits to interesting buckling phenomena.

In Fig. 3, we demonstrate how changing the value of $\mathcal{V}$ affects the amount of buckling a filament experiences. The fiber is perturbed with $\theta_{0}=0.9\pi$ and $\Delta\theta_{0}=0.1$. As the filament rotates, buckling modes form, which are directly linked to the size of $\mathcal{V}$ and the choice of perturbation (Eq. (64)). Large values produce high-order buckling modes, which can be seen in the fourth row of Fig. 3. Conversely, comparatively small values produce no buckling and the filament rotates through a standard Jeffery orbit. For $\mathcal{V}=5\times 10^{3}$, a first-order buckling mode begins to form, visible in the first row of Fig. 3. The amplitude of the buckling increases as $\mathcal{V}$ increases until higher-order modes are induced. Tornberg & Shelley tornberg2004simulating examined buckling governed by an effective viscosity parameter $\bar{\mu}=8\pi\mathcal{V}$, producing filament shapes akin to the second row of Fig. 3.

The buckling problem of a flexible filament in shear flow has multiple solutions becker2001instability ; liu2018morphological . As noted by Tornberg & Shelley tornberg2004simulating , the choice of initial perturbation to the filament shape can preferentially lead to one of the solutions. For example, the particular initial condition considered in tornberg2004simulating produces a reflected filament shape when changing the sign of the perturbation. Additionally, for the high-order buckling modes, increasing values of $Q$ must be chosen in order to resolve the filament dynamics. However, given an initial condition that does not uniquely determine a specific solution branch (as for this problem with large values of $\mathcal{V}$), the buckled solution becomes sensitive to the choice of discretization. Details of the convergence of the method for this problem is contained within Sec. SI.1 of the Supplemental Material.

At critical values of $\mathcal{V}$, unstable buckling can occur. For example, when ${\mathcal{V}=3.4\times 10^{4}}$, a higher-order buckling mode initially forms (row three of Fig. 3) but collapses into a lower-order mode. In Fig. 4, the filament evolution highlighted in Fig. 3 is displayed as a curvature plot in arclength and time. The unstable higher-order mode can be clearly seen collapsing into a lower-order configuration, indicated by the two dotted lines in Fig. 4a. Increasing $\mathcal{V}$ further ensures that a higher-order mode is stable throughout the compression phase (Fig. 4b).

The perturbation to the filament shape induced by buckling can be investigated by considering the evolving body-frame tangent angles ${\tilde{\theta}^{\,[m]}_{\text{body}}=\tilde{\theta}^{[m]}(t)-\overline{\theta}(t)}$, where

[TABLE]

is the mean filament angle. By fitting Chebyshev polynomials to $\tilde{\theta}^{\,[m]}_{\text{body}}$ along the filament at a time $t$, allowing for $5\%$ interpolation error, we can assess the evolution of both the order of Chebyshev polynomials required and their associated magnitude. In Figs. 5a and 5b, we show the results of this fitting process for two choices of $\mathcal{V}$. An increase in $\mathcal{V}$ requires a commensurate increase in polynomial order required, illustrated by examining the Chebyshev coefficients in Figs. 5c and 5d.

V.2 Results: a single passive filament sedimenting under gravity

Following Sec. III.2, simulations of a single passive filament sedimenting under gravity (with no background flow) are considered. Filament dynamics in this case are simulated by solving Eq. (58) and characterized by the dimensionless elasto-gravitational parameter $\mathcal{G}$ (Eq. (53)).

By sampling $\theta$ from the very-low amplitude parabola ${y=1\times 10^{-7}x^{2}}$, the filament geometry is initialized with unit arclength, and pre-solved with a coarse discretization ($Q=10$) until the shape is sufficiently curved so that a higher-resolution representation can be employed (full details are provided in the Supplemental Material Sec. SI.2). In the following results, $Q=40$ is used for the upscaled initial condition.

For different choices of $\mathcal{G}$, various sedimentary buckling modes can be observed. In Fig. 6, a stable “U” filament configuration emerges over time. For large values of $\mathcal{G}>3000$, a metastable “W” configuration forms, which transitions into a stable “U” horseshoe shape (see Fig. 7, in which $\mathcal{G}=3500$). Such behavior has previously been observed by Delmotte et al. delmotte2015general and Cosentino Lagomarsino et al. consentino2005hydrodynamic , who witnessed buckling at the same threshold value for their identical elasto-gravitational parameter. This transition shifts the filament’s center of gravity to the left, creating an asymmetry which is then partially resolved as the horseshoe equilibrium configuration is reached.

V.3 Results: a single active filament

In the following section, we consider swimming in a stationary fluid caused by two types of traveling-wave moment densities, (1) sperm-like sinusoidal motion $\mathfrak{m}_{1}(s,t)$, and (2) a worm-like beating pattern $\mathfrak{m}_{2}(s,t)$, given respectively by

[TABLE]

where $\mathfrak{m}_{0}$ and $k$ are the dimensionless amplitude and wave number respectively. These swimmers are initialized by sampling $\theta$ from a low-amplitude parabola $y=1\times 10^{-3}x^{2}$ and constructed as before, ensuring unit arclength.

For appropriately small choices of $\mathfrak{m}_{0}$, substitution of Eq. (66) or (67) into Eq. (62) can induce swimming in a filament in stationary surrounding flow. Fixing the wave number $k$, the relationship between filament elasticity and choice of driving force (governed by $\mathfrak{m}_{0}$) is investigated. The Velocity Along a Line (VAL) is a measure of the swimming speed of a filament for a chosen $\mathcal{S}$ and $\mathfrak{m}_{0}$ pair, calculated via

[TABLE]

in which $T=2\pi$ is the period of the driving wave and $\bm{X}_{0}^{(j)}$ represents the position of the leading point of the filament after it has traveled $j$ wavelengths, with $j$ chosen such that the filament has established a regular motion after beginning to swim.

Swimming speed for different choices of parameter pairs $\left(\mathcal{S},\mathfrak{m}_{0}\right)$ are presented in Fig. 8. For critical values, filaments self-intersect, in which case the EIF is inapplicable. The shape of such filaments are shown in Fig. 9. For a sperm-like swimmer (left panel of Fig. 8), swimming speed increases as $\mathfrak{m}_{0}$ and $\mathcal{S}$ are increased in tandem for a fixed wave-number $k=4\pi$. This is in contrast to worm-like swimmers (right panel of Fig. 8), in which there is a clear optimal choice of $\mathfrak{m}_{0}$ for a given $\mathcal{S}$ to induce fastest swimming. Shape profiles for the fastest swimmers of each swimming type are also displayed in Fig. 9.

V.4 Results: multiple passive filaments in shear flow

Motivated by the work of Young young2009hydrodynamic , two filaments of equal length are placed into a linear shear flow so that their initial midpoints ${\bm{X}(s=1/2,t=0)}$ intersect the line $y=0$, and are separated by a distance $\Delta X_{0}$. The initial shape profiles for each filament are the same as in Sec. V.1. For two filaments characterized by the same viscous-elastic parameter, non-local hydrodynamic interactions lead to geometric asymmetry, as seen in Fig. 10. The value of $\mathcal{V}$ determines the level of dissimilarity between filaments with identical initial shape profiles, with higher values leading to larger deviations in shape throughout the rotation.

The initial separation distance also has a significant influence on the dynamics of each filament. The shapes of initially close filaments evolve in tandem, assuming similar geometries at a given moment in time (row one of Fig. 10). Increasing the separation distance results in a decoupling of the filament shape profiles (row two of Fig. 10). The EIF allows for each filament’s characteristic parameter to be independently chosen, and in Fig. 11, we highlight the variation in shape this can induce. In each of the panels in Figs. 10 and 11, color denotes the magnitude of the relative fluid velocity disturbance (found by subtracting the shear fluid velocity from the resultant fluid velocity), scaled with respect to the size of fluid velocity in the bottom left corner of each frame. Fluid lines show the instantaneous direction of the flow disturbance.

The geometric similarity between filaments can be assessed via their Procrustes score gower1975generalized , calculated using the procrustes MATLAB® command. We compare pairs of filaments initially separated by $\Delta X_{0}=0.5$ and $1$ for a range of $\mathcal{V}$, with results displayed in Fig. 12. As $\mathcal{V}$ is increased, hydrodynamic interactions through the fluid cause larger levels of asymmetry between the filaments during the period of maximum buckling (when they are perpendicular to the direction of shear at $t\approx 3$). Increasing the initial filament separation results in a higher baseline Procrustes score, but with reduced variability, demonstrating how non-local hydrodynamic interactions decay as the filaments are moved apart.

V.5 Results: multiple passive filaments sedimenting under gravity

A range of filament systems with multiple choices of the elasto-gravitational parameter $\mathcal{G}$ are presented in Figs. 13, 14 and 15. In each case, the initial filament shape configurations are as those in Sec. V.2, and are separated from each other by $\Delta\bm{X}_{0}=[\Delta X_{0},\Delta Y_{0}]$. A stopping criterion is implemented which halts the integration when filaments intersect or otherwise touch. For small values of $\mathcal{G}$ filaments slide towards each other as they sediment due to the anisotropy of Stokes drag. For larger values of $\mathcal{G}$ (rows two and three of Fig. 13) the metastable and stable buckling modes observed in the single filament experiments (Figs. 6 and 7) are replicated. In these systems, the interaction of multiple filaments leads to the steady-state profiles being reached sooner than in the single-filament cases. Furthermore, the onset of buckling occurs at reduced values of $\mathcal{G}$ when multiple filaments are interacting.

As was the case for multiple filaments in shear flow, the inter-filament spacing $\Delta\bm{X}_{0}$ also has a significant effect on the resulting group dynamics, and can lead to symmetry breaking in the arrangement of filaments. In Fig. 14, two initial filament placement configurations are tested. For each setup, we initialize four filaments arranged in a grid, with $\Delta X_{0}=1.5$ and vertical separations $\Delta Y_{0}=0.5$ (row one of Fig. 14) and $\Delta Y_{0}=1.5$ (row two of Fig. 14). Smaller vertical spacing leads to the filaments “nestling” in a horizontally-mirrored configuration; when placed further apart, horizontal mirroring still occurs, but the filaments do not approach each other.

In Fig. 15 nine identical filaments are arranged in a grid with initial spacing ${\Delta X_{0}=1.5}$ and $\Delta Y_{0}=1$. As with the smaller arrays, filaments can group and buckle according to the choice of characteristic parameter $\mathcal{G}$. Although identical in governing parameter $\mathcal{G}$, the non-local hydrodynamic interactions between filaments induce different buckling behaviors depending on location within the array.

For all choices of $\mathcal{G}$, vertical symmetry-breaking occurs between the top-most and middle rows of filaments, with the second and third rows of filaments on the left and right of the frame nestling into those below them. For $\mathcal{G}=3500$, more prominent buckling is apparent, with filaments in the central column approaching a “W” shape, whereas those on the flanks assume a “U” configuration. Competing interactions between filaments in the central column and their surrounding neighbors causes them to remain in a metastable “W” shape longer than would be expected in the single-filament case.

V.6 Results: multiple active filaments

Multiple self-propulsive active filaments can be simulated using the EIF method. To illustrate this, we consider two filaments swimming alongside each other at different separation distances $\Delta Y_{0}$, with shapes initialized as in Sec. III.3. For fixed $\mathcal{S}=8$ and $\mathfrak{m}_{0}=0.05$, the size of the separation determines the resulting average swimming velocity, as shown in Fig. 16. At low levels of separation, hydrodynamic interactions between the synchronous beats result in a higher average VAL (Eq. (68)), which decays as the filaments are placed further apart (Fig. 16a).

VI Discussion

In this manuscript, the elastohydrodynamic integral equation (EIF) framework is formulated, and applied to problems involving single and multiple filaments in shear flows, sedimenting under gravity, and swimming due to a prescribed active moment. From the simulations presented in Sec. V, it is apparent how inter-filament non-local hydrodynamic interactions have a role in governing filament shapes and buckling behavior. By examining active moment driven swimmers, an optimum pairing of moment-amplitude and characteristic swimming number is revealed.

One of the key benefits to the integral formulation framework presented is that the need for computing Lagrange multipliers of tension is removed; inextensibility of the modeled filament is already guaranteed by the use of the method of lines discretization. This results in a method that is more computationally efficient than many similar contemporary models. The EIF method proposed by Moreau and co-authors moreau2018asymptotic highlighted the reduction in computational runtime attainable when using the integral formulation applied to elastodynamics, due in part to the reduction in the required degrees of freedom the integral formulation affords. Alleviating this numerical stiffness facilitates the study of non-linear problems, such as those involving filament buckling investigated in Sec. V. Additional computational costs are incurred by introducing non-local hydrodynamics to the EIF framework, in line with the increase in the dimension and density of the matrix system solved at each time-step. The BLM and the EIF with regularized stokeslets complete a relaxing rod simulation in about equivalent walltime (as expected, since both approaches compute three-dimensional non-local hydrodynamics), shown in Fig. 2c. However, to ensure that the point-wise error is small, fine discretizations need to be used with the bead model ($Q>60$), whereas the EIF method performs equally well with a much coarser segmentation of the filament ($Q>20$). This allows for adequately accurate results to be obtained for reduced computational cost when using the EIF with regularized stokeslets over the bead and link model.

The proposed method has the potential to quickly and accurately simulate arrays of filaments in various flows and surroundings. The EIF method presented here will enable the solution of more challenging problems such as planar beating cilia sheets, or multiple sperm swimming in a narrow channel. The methods’ modular framework enables such problems to be setup and executed in a simple and obvious manner. Additionally, the EIF method fits into the family of regularized stokeslet methods cortez2005method , which are increasingly widely used for problems in biological fluid mechanics. Crucially, the computational efficiency exhibited in the single filament experiments is retained. The method could be extended by incorporating the treecode formulation of Wang et al. wang2018treecode , or equivalent methods, to coarse-grain the far-field flow, simplifying the computational cost and enabling simulation of increasingly large numbers of filaments. The inclusion of a repulsive force term, such as the Lennard-Jones potential employed by Jarayaman et al. jayaraman2012autonomous , could additionally enable simulation of filaments in close proximity and help avoid filament self-intersection.

The results of Sec. V.1 suggest that the shape of a buckling filament can be described to a high degree of accuracy by a relatively low order Chebyshev polynomial (Fig. 5). These results suggest that a modified discretization based on orthogonal polynomials, perhaps in concert with suitable quadrature techniques muldowney1995spectral might provide further improvements in efficiency and scalability. In summary we hope that the integral operator formulation of elastohydrodynamics will be valuable in biological fluid dynamics and beyond.

VII Acknowledgments

D.J.S. and M.T.G. acknowledge funding from the Engineering and Physical Sciences Research Council (EPSRC), Healthcare Technologies Challenge Award (EP/N021096/1). T.D.M-J. acknowledges funding from the EPSRC (EP/R041555/1). A.L.H-M. acknowledges support from the EPSRC for funding via a PhD scholarship (EP/N509590/1).