Efficient simulation of filament elastohydrodynamics in three dimensions

TL;DR

This paper introduces an efficient computational method for simulating the complex three-dimensional motion of slender, flexible filaments in viscous fluids, overcoming previous limitations to planar configurations.

Contribution

It extends coarse-graining simulation frameworks to three dimensions, addressing coordinate singularities and enabling rapid, accurate modeling of filament dynamics in complex environments.

Findings

Enables efficient 3D filament simulations

Circumvents coordinate-system singularities

Applicable to various biological and physical systems

Abstract

Fluid-structure simulations of slender inextensible filaments in a viscous fluid are often plagued by numerical stiffness. Recent coarse-graining studies have reduced the computational requirements of simulating such systems, though have thus far been limited to the motion of planar filaments. In this work we extend such frameworks to filament motion in three dimensions, identifying and circumventing coordinate-system singularities introduced by filament parameterisation via repeated changes of basis. The resulting methodology enables efficient and rapid study of the motion of flexible filaments in three dimensions, and is readily extensible to a wide range of problems, including filament motion in confined geometries, large-scale active matter simulations, and the motility of mammalian spermatozoa.

Equations

$${\mathbf{n}}_s-\mathbf{f}$$

$${\mathbf{m}}_s+{\mathbf{x}}_s\times \mathbf{n}-\mathbf{\tau }$$

$$\mathbf{n}(0)=\mathbf{n}(L)=\mathbf{m}(0)=\mathbf{m}(L)=0.$$

$$\mathbf{n}(s)=-\underset{s}{\overset{L}{\int }}\mathbf{f}(\tilde{s})\mathrm{d}\tilde{s}.$$

$$0=\underset{0}{\overset{L}{\int }}\mathbf{f}(\tilde{s})\mathrm{d}\tilde{s},$$

$$-\underset{s}{\overset{L}{\int }}\left[(\mathbf{x}(\tilde{s})-\mathbf{x}(s))\times \mathbf{f}(\tilde{s})+\mathbf{\tau }(\tilde{s})\right]\mathrm{d}\tilde{s}=\mathbf{m}(s).$$

$$\frac{\partial {\mathbf{d}}_{\alpha }}{\partial s}=\mathbf{\kappa }\times {\mathbf{d}}_{\alpha }$$

$$\mathbf{m}=EI\left({\kappa }_1{\mathbf{d}}_1+{\kappa }_2{\mathbf{d}}_2+\frac{1}{1+\sigma }{\kappa }_3{\mathbf{d}}_3\right),$$

$$-{\mathbf{d}}_{\alpha }(s)\cdot \underset{s}{\overset{L}{\int }}\left[(\mathbf{x}(\tilde{s})-\mathbf{x}(s))\times \mathbf{f}(\tilde{s})+\mathbf{\tau }(\tilde{s})\right]\mathrm{d}\tilde{s}=\frac{EI}{1+{\delta }_{\alpha ,3}\sigma }{\kappa }_{\alpha }(s),$$

$$-{\mathbf{d}}_{\alpha }^i\cdot \sum _{j=i}^N\underset{s_j}{\overset{s_{j+1}}{\int }}\left[(\mathbf{x}(\tilde{s})-{\mathbf{x}}_i)\times \mathbf{f}(\tilde{s})+\mathbf{\tau }(\tilde{s})\right]\mathrm{d}\tilde{s}=\frac{EI}{1+{\delta }_{\alpha ,3}\sigma }{\kappa }_{\alpha }(s_i),$$

$$-{\mathbf{d}}_{\alpha }^i\cdot \left({\mathbf{I}}_i^f+{\mathbf{I}}_i^{\tau }\right)=\frac{EI}{1+{\delta }_{\alpha ,3}\sigma }{\kappa }_{\alpha }(s_i),$$

$${\mathbf{I}}_i^f=\sum _{j=i}^N\left\{\left[\frac{\Delta s}{2}\left({\mathbf{x}}_j-{\mathbf{x}}_i\right)+\frac{\Delta s^2}{6}{\mathbf{d}}_3^j\right]\times {\mathbf{f}}_j+\left[\frac{\Delta s}{2}\left({\mathbf{x}}_j-{\mathbf{x}}_i\right)+\frac{\Delta s^2}{3}{\mathbf{d}}_3^j\right]\times {\mathbf{f}}_{j+1}\right\},$$

$${\mathbf{I}}_i^{\tau }=\sum _{j=i}^N\Delta s{\mathbf{\tau }}_j.$$

$$-\frac{\Delta s}{2}\sum _{j=1}^N({\mathbf{f}}_j+{\mathbf{f}}_{j+1})=\mathbf{n}(0).$$

$$-\mathcal{B}\left[\begin{array}{c}\mathbf{F}\\ \mathbf{T}\end{array}\right]=\mathbf{R},$$

$${\mathcal{B}}_1$$

$${\mathcal{B}}_2$$

$${\mathcal{B}}_3$$

$${\mathbf{d}}_{\alpha }^i\cdot {\mathbf{I}}_i^f$$

$$=\sum _{j=i}^N\left\{{\mathbf{d}}_{\alpha }^i\times \left[\frac{\Delta s}{2}\left({\mathbf{x}}_j-{\mathbf{x}}_i\right)+\frac{\Delta s^2}{6}{\mathbf{d}}_3^j\right]\cdot {\mathbf{f}}_j+{\mathbf{d}}_{\alpha }^i\times \left[\frac{\Delta s}{2}\left({\mathbf{x}}_j-{\mathbf{x}}_i\right)+\frac{\Delta s^2}{3}{\mathbf{d}}_3^j\right]\cdot {\mathbf{f}}_{j+1}\right\},$$

$${\mathbf{d}}_{\alpha }^i\cdot {\mathbf{I}}_i^{\tau }=\Delta s{\mathbf{d}}_{\alpha }^i\cdot \sum _{j=i}^N{\mathbf{\tau }}_j,$$

$$\mathbf{R}=\frac{EI}{1+{\delta }_{\alpha ,3}\sigma }[0,0,0,{\kappa }_1(s_1),{\kappa }_2(s_1),{\kappa }_3(s_1),{\kappa }_1(s_2),\dots ,{\kappa }_3(s_N){]}^{\top },$$

$$f_t=-C_t{u}_t,f_n=-C_n{u}_n.$$

$$C_t=\frac{2\pi \mu }{\log \left(2L/ϵ\right)-0.5},C_n=\frac{4\pi \mu }{\log \left(2L/ϵ\right)-0.5},$$

$$\dot{\mathbf{X}}=A\mathbf{F},$$

$$\bm{\tau}_{i}=4\pi\mu\epsilon^{2}\omega_{i}\bm{d}_{3}^{i}\$$

$$-\mathcal{B}\left[\begin{array}[]{cc}A^{-1}&0\\ 0&\tilde{A}\end{array}\right]\left[\begin{array}[]{c}\dot{\bm{X}}\\ \bm{\omega}\end{array}\right]=-\mathcal{B}\mathcal{A}\left[\begin{array}[]{c}\dot{\bm{X}}\\ \bm{\omega}\end{array}\right]=\bm{R}\,,$$

$${\mathbf{d}}_1^i$$

$${\mathbf{d}}_2^i$$

$${\mathbf{d}}_3^i$$

$${\theta }_i=\arccos \left({\mathbf{d}}_3^i\cdot {\mathbf{e}}_z\right),{\varphi }_i=\arctan \left(\frac{{\mathbf{d}}_3^i\cdot {\mathbf{e}}_y}{{\mathbf{d}}_3^i\cdot {\mathbf{e}}_x}\right),{\psi }_i=\arctan \left(\frac{{\mathbf{d}}_1^i\cdot {\mathbf{e}}_z}{{\mathbf{d}}_2^i\cdot {\mathbf{e}}_z}\right).$$

$${\mathbf{x}}_j={\mathbf{x}}_1+\Delta s\sum _{i=1}^{j-1}{\mathbf{d}}_3^i,{\dot{\mathbf{x}}}_j={\dot{\mathbf{x}}}_1+\Delta s\sum _{i=1}^{j-1}{\dot{\mathbf{d}}}_3^i.$$

$$Q\dot{\mathbf{\Theta }}=\dot{\mathbf{X}},$$

Full text

I Introduction

The coupled elasticity and hydrodynamics of flexible inextensible filaments on the microscale are of significance to much of biology, biophysics and soft matter physics. For example, many organisms possess slender flagella or cilia, utilised for driving flows and even locomotion, whilst investigation into the role of synthetic filaments as both soft deformable sensors and methods of propulsion has been the subject of recent enquiry [1, 2, 3, 4, 5, 6, 7, 8, 9]. As a result, the complex mechanics of fluid-structure interaction has been well-studied, utilising methods such as the slender body and resistive force theories of Hancock [2], Gray and Hancock [10], Johnson [11], through to the exact representations of boundary integral methods as used by Pozrikidis [12, 6, 7]. A fundamental barrier to much numerical investigation has been the severe stiffness associated with the equations of filament elasticity when coupled to viscous fluid dynamics. Hence, as remarked in the recent and extensive review of du Roure et al. [13], an appropriate framework, capable of realising efficient simulation of filament elastohydrodynamics, is crucial for the numerical study of filament mechanics.

Recently, significant progress has been made in resolving the dynamics of planar filaments, with the work of Moreau et al. [14] presenting a coarse-grained model of filament elasticity that overcame much of the stiffness previously associated with slender elastohydrodynamics. Key to this approach was the integration of the pointwise force and moment balance equations, the spatial discretisation of which yielded a relatively simple system of ordinary differential equations to solve in order to describe filament motion. Additionally demonstrated to be flexible in the original publication of Moreau et al., this framework has been extended to include non-local hydrodynamics in both infinite and semi-infinite domains, and applied to a variety of single and multi-filament problems [15, 16]. However, being confined to two dimensions limits the potential scope and applicability of these approaches, with three dimensional filament motion being readily and frequently observed in a plethora of biophysical systems, such as the complex flagellar beating found in spermatozoa or the helically-driven monotrichous bacterium Escherichia coli [17, 3].

However, for non-planar filaments in three dimensions there is currently no methodology analogous to that of Moreau et al., with state-of-the-art frameworks still plagued by extensive numerical stiffness, necessitating costly computation to the extent that practical simulation studies have been limited and parameter space studies are largely prohibited. With three dimensions inherently more challenging than lower dimensional settings, this field has seen developments such as the recent work of Schoeller et al. [18], which utilises a quarternion representation of filament orientation to parameterise the three dimensional shape of the slender body. However, in this framework, numerical care is required to satisfy the inextensibility condition, with similar such consideration necessary in the earlier methodologies of Olson et al. [19], Simons et al. [20], Ishimoto and Gaffney [21], Bouzarth et al. [22], each of which are equipped with non-local slender-body hydrodynamics and consider nearly inextensible filaments. Consequently, these existing approaches often require the use of sophisticated computing hardware in order to simulate filament motion, with typical simulations of Ishimoto and Gaffney having a runtime of multiple hours on high performance computing clusters. The recent work of Jabbarzadeh and Fu [23] compared and contrasted these nearly inextensible approaches with a truly inextensible scheme, concluding that both accuracy and efficiency was afforded by the latter in a range of biological and biophysical modelling scenarios. Despite their improved efficiency, typical walltimes for these filament simulations are measured on a timescale of hours on typical hardware. Thus, there remains significant scope for the development of an efficient framework for the simulation of inextensible elastic filaments in three dimensions, one in which filament dynamics can be rapidly computed on non-specialised hardware on timescales of seconds or minutes, thus facilitating a wealth of future studies and explorations into complex and previously intractable biological and physical systems.

Hence, the fundamental objective of this study is to develop and describe an efficient framework for the numerical simulation of filament mechanics in three dimensions. We will build upon the recent and significant work of Moreau et al. [14], extending their approach to include an additional spatial dimension via a generalisation of the Frenet triad and integration of the governing equations of elasticity. We will overcome fundamental issues with simple single parameterisations of a filament in three dimensions, presenting an effective computational approach utilising adaptive reparameterisation and basis selection. We will then validate the presented framework by consideration of three candidate test problems, simulating well-documented behaviours of filaments in a viscous fluid and including a side-by-side comparison against the existing and recent methodology of Ishimoto and Gaffney [21]. Finally, we will showcase the flexibility and general applicability of the presented approach by describing and exemplifying a number of methodological extensions.

II Methods

II.1 Equations of elasticity

We consider a slender, inextensible, unshearable filament in a viscous Newtonian fluid, with its centreline described by $\bm{x}(s)$, parameterised by arclength $s\in[0,L]$ for dimensional filament length $L$. We model the filament as a Kirchhoff rod with arclength-independent material parameters, circular cross-sections, and, in the first instance, no intrinsic curvature or intrinsic torsion. Both the filament and fluid inertia are negligible for the physical scales associated with many applications, especially those associated with cellular flagella and cilia; hence there is no inertia here and throughout. Along the filament we have the pointwise conditions of force and moment balance, given explicitly by

[TABLE]

for contact force and torque denoted $\bm{n},\bm{m}$ respectively and where a subscript of $s$ denotes differentiation with respect to arclength. The quantity $\bm{f}$ is the force per unit length applied on the fluid medium by the filament, which we will later express in terms of the filament velocity $\dot{\bm{x}}$, where here a dot denotes a time derivative. Similarly, $\bm{\tau}$ is the torque per unit length applied on the fluid medium by the filament. Whilst any standard boundary condition may be considered in the formalism below, throughout we impose zero force and torque at the filament ends so that

[TABLE]

In the Kirchhoff framework note that the contact force is not constitutive, but simply an undetermined Lagrange multiplier for the intrinsic constraints of inextensibility and no shearing of the filament cross section in any direction [24]. Thus we may eliminate $\bm{n}(s)$ at the earliest opportunity; using the boundary condition $\bm{n}(L)=0$ and Eq. 1 we have

[TABLE]

This relation is also subject to the constraint that the boundary condition $\bm{n}(0)=\bm{0}$ is satisfied, generating a global force balance constraint for the applied force per unit length,

[TABLE]

which we carry forward into the formalism below together with the elimination of $\bm{n}(s)$ via Eq. 4. Thus, and as originally considered in Moreau et al. [14], integration of Eq. 2 and use of the boundary condition $\bm{m}(L)=\bm{0}$, reveals the integrated moment balance

[TABLE]

Given a right-handed orthonormal director basis $\{\bm{d}_{1}(s),\bm{d}_{2}(s),\bm{d}_{3}(s)\}$, generalising the Frenet triad such that $\bm{d}_{3}$ corresponds to the local filament tangent, following Nizette and Goriely [25] we define the twist vector $\bm{\kappa}$ by

[TABLE]

for $\alpha=1,2,3$. Writing $\bm{\kappa}=\sum_{\alpha}\kappa_{\alpha}\bm{d}_{\alpha}$, for bending stiffness $EI$ we use the Euler-Bernouilli constitutive relation of the Kirchhoff formalism to relate the contact torque $\bm{m}$ and the twist vector $\bm{\kappa}$ [24], via

[TABLE]

where $\sigma$ is the Poisson ratio [25], assumed to be constant. With this constitutive relation the integrated moment balance equations in the $\bm{d}_{\alpha}$ directions are simply

[TABLE]

for $\alpha=1,2,3$ and where $\delta_{a,b}$ denotes the Kronecker delta.

II.2 Filament discretisation

In discretising the filament we follow the approach of Walker et al. [16], as previously applied to planar filaments and itself building upon the earlier work of Moreau et al. [14]. We approximate the filament shape with $N$ piecewise-linear segments, each of constant length $\Delta s$, with segment endpoints having positions denoted by $\bm{x}_{1},\ldots,\bm{x}_{N+1}$ and the constraints of inextensibility and the absence of cross section shear are satisfied inherently. The endpoints of the $i$th segment correspond to $\bm{x}_{i}$ and $\bm{x}_{i+1}$ for $i=1,\ldots,N$, with the local tangent $\bm{d}_{3}$ being constant on each segment and denoted $\bm{d}_{3}^{i}$. In what follows we will consider a discretisation of $\bm{d}_{1},\bm{d}_{2}$ such that they are also constant on each segment, and we denote these constants similarly as $\bm{d}_{1}^{i},\bm{d}_{2}^{i}$. Writing $s_{i}$ for the constant arclength associated with each material point $\bm{x}_{i}$, we apply Eq. 9 at each of the $s_{i}$ for $i=1,\ldots,N$, splitting the integral at the segment endpoints to give

[TABLE]

for $\alpha=1,2,3$. On the $j$th segment, $\bm{x}$ may be written as $\bm{x}(s)=\bm{x}_{j}+\eta(\bm{x}_{j+1}-\bm{x}_{j})$, where $\eta\in[0,1]$ is given by $\eta=(s-s_{j})/\Delta s$. Additionally discretising the force per unit length as a continuous piecewise-linear function, with $\eta$ as above we have $\bm{f}(s)=\bm{f}_{j}+\eta(\bm{f}_{j+1}-\bm{f}_{j})$ on the segment, where we write $\bm{f}_{j}=\bm{f}(s_{j})$. Substitution of these parameterisations into Eq. 10 and subsequent integration yields, after simplification,

[TABLE]

where the integral contribution of the force and torque densities are denoted $\bm{I}_{i}^{f}$ and $\bm{I}_{i}^{\tau}$ respectively. With this discretisation $\bm{I}_{i}^{f}$ has reduced to

[TABLE]

in agreement with expressions for planar filaments found in Moreau et al. [14], Walker et al. [16]. As we will highlight below, the contributions of the applied torque per unit length are relatively small given the slenderness of the filament, motivating a less refined discretisation for $\bm{I}_{i}^{\tau}$. Hence, taking the piecewise constant discretisation $\bm{\tau}=\bm{\tau}_{j}$ on the $j$th segment, we have the simple expression

[TABLE]

From the above we see explicitly that the integral component of each moment balance equation may be written as a linear operator acting on the $\bm{f}_{j}$ and the $\bm{\tau}_{j}$. Similarly, with this piecewise-linear force discretisation the integrated force balance of Eq. 1 simply reads

[TABLE]

We write $\bm{F}=[f_{1,x},f_{1,y},f_{1,z},\ldots,f_{N+1,x},f_{N+1,y},f_{N+1,z}]^{\top}$ for components $f_{j,x},f_{j,y,}f_{j,z}$ of $\bm{f}_{j}$ with respect to some fixed laboratory frame with basis $\{\bm{e}_{x},\bm{e}_{y},\bm{e}_{z}\}$, and similarly $\bm{T}$ for the vector of components of applied torque per unit length. Here and throughout, forces and torques are written with respect to the laboratory reference frame. With this notation, we may write the equations of force and moment balance as

[TABLE]

where $\mathcal{B}$ is a matrix of dimension $(3N+3)\times(6N+3)$ with rows $\mathcal{B}_{k}$. For $k=1,2,3$ the first $3N+3$ columns are given by

[TABLE]

and correspond to the force balance of Eq. 14, with the remaining $3N$ columns zero. The remaining rows of $\mathcal{B}$ encode the moment balance of Eq. 11 as expanded in Eq. 12, organised in triples such that $\mathcal{B}_{3(i-1)+3+\alpha}$ projects the $i$th moment balance equation onto $\bm{d}_{\alpha}^{i}$, in that this $(3i+\alpha)$th row of $\mathcal{B}$ captures the $\bm{d}_{\alpha}^{i}$ component of $-(\bm{I}_{i}^{f}+\bm{I}_{i}^{\tau})$. The cross products inherited from Eq. 12 may now be notationally simplified by use of the cyclic property of the scalar triple product, explicitly giving

[TABLE]

with the latter expression readily transcribed as a linear operator acting on the $\bm{f}_{j}$ for $j=i,\ldots,N$. Analogously, we have

[TABLE]

from which a linear operator acting on the $\bm{\tau}_{j}$ for $j=i,\ldots,N$ can be constructed. Accordingly, the $(3N+3)$-vector $\bm{R}$ is given by

[TABLE]

so that the local moment balance is expressed relative to the local director basis. We remark that each of the quantities involved in the construction of $\mathcal{B}$ and $\bm{R}$ are well-defined for a general filament in three dimensions, given the local directors $\bm{d}_{1}$ and $\bm{d}_{2}$ and computing the components of the twist vector as $\kappa_{1}=\bm{d}_{3}\cdot\partial_{s}\bm{d}_{2}$, $\kappa_{2}=\bm{d}_{1}\cdot\partial_{s}\bm{d}_{3}$, and $\kappa_{3}=\bm{d}_{2}\cdot\partial_{s}\bm{d}_{1}$. In terms of the discretised filament, these arclength derivatives are approximated via finite differences in practice. Additionally, we will proceed assuming that the filament is moment-free at the base, which additionally enforces $\kappa_{1}(0)=\kappa_{2}(0)=\kappa_{3}(0)=0$.

II.3 Coupling hydrodynamics

We now relate the force density $\bm{f}$ acting on the fluid to the velocity of each segment endpoint, utilising the commonly-applied method of resistive force theory as introduced by Hancock [2], Gray and Hancock [10] and adopted by Moreau et al. [14] for planar filaments, incurring typical errors logarithmic in the aspect ratio of the filament. Here taking the radius of the filament to be $\epsilon=10^{-2}L$, which more generally is assumed to be small in comparison to the filament length, simple resistive force theory gives the leading order relation between filament velocity and force density as

[TABLE]

Here $f_{t}$ and $f_{n}$ denote the components of the force density tangential and normal to the filament, with analogous definitions of $u_{t}$ and $u_{n}$. We will utilise the expression of Gray and Hancock [10], with

[TABLE]

where $\mu$ is the medium viscosity, noting the relation $C_{n}=2C_{t}$. We approximate the local filament tangent at the segment endpoint $\bm{x}_{i}$ as the average of $\bm{d}_{3}^{i-1}$ and $\bm{d}_{3}^{i}$ for $i=2,\ldots,N$, with the tangent for $i=1$ and $i=N+1$ simply being taken as $\bm{d}_{3}^{1}$ and $\bm{d}_{3}^{N}$ respectively. By linearity, and again assuming a piecewise-linear force density along segments, we may write the coupling of translational kinematics to hydrodynamics as

[TABLE]

where $A$ is a square matrix of dimension $3(N+1)\times 3(N+1)$ and is a function only of the segment endpoints $\bm{x}_{i}$. Of dimension $3(N+1)$, the vector $\dot{\bm{X}}$ corresponds to the linear velocities of the segment endpoints, and is constructed analogously to $\bm{F}$ with respect to the laboratory frame. This relation results from the application of the no-slip condition at the segment endpoints, coupling the filament to the surrounding fluid.

In order to relate the rate of rotation of each segment to the viscous torque $\bm{\tau}_{i}$ acting on it, we here consider an approximation of the finite segment as an infinite rotating cylinder, associating the torque per unit length on the $i$th segment with the rotation $\omega_{i}$ about its local tangent $\bm{d}_{3}^{i}$ via the relation of Chwang and Wu [26]:

[TABLE]

and in particular the $\epsilon^{2}$ scaling entails the torque per unit length contributions are relatively small. Here we recall that $\mu$ is the viscosity of the fluid medium, and $\epsilon$ is the radius of the filament. We may write this relation as a linear operator on $\bm{\omega}=[\omega_{1},\ldots,\omega_{N}]^{\top}$, written simply as $\bm{T}=\tilde{A}\bm{\omega}$. This crude approximation may readily be substituted for non-local hydrodynamics via the method of regularised Stokeslet segments, which will likely be a topic of future work. Similarly, non-local hydrodynamics may be utilised in place of Eq. 23, as used for two-dimensional filament studies by Hall-Mcnair et al. [15] and Walker et al. [16], the latter incorporating a planar no-slip boundary and still yielding an explicit linear relation analogous to Eq. 23.

Combining Eqs. 15, 23 and 24 yields the linear system

[TABLE]

where $A$ is invertible and $\mathcal{A}$ is defined to be a block matrix of dimension $(6N+3)\times(4N+3)$ with non-zero blocks $A^{-1}$ and $\tilde{A}$.

II.4 Parameterisation

We may parameterise the tangents $\bm{d}_{3}^{i}$ on each linear segment by the Euler angles $\theta_{i}\in[0,\pi]$, $\phi_{i}\in(-\pi,\pi]$, $\psi_{i}\in(-\pi,\pi]$ for $i=1,\ldots,N$ [24]. With this parameterisation we may make a choice of $\bm{d}_{1}$ and $\bm{d}_{2}$, taking here the three orthonormal vectors to be

[TABLE]

written with respect to the laboratory frame and where $s_{\theta}{}\equiv\sin{\theta_{i}}$, $c_{\theta}{}\equiv\cos{\theta_{i}}$, and analogously for $s_{\phi}{},c_{\phi}{},s_{\psi}$ and $c_{\psi}$. From the directors we recover

[TABLE]

As the discretised filament is piecewise linear, for $j=1,\ldots,N+1$ we may write

[TABLE]

With $\bm{d}_{3}^{i}$ parameterised as above, we can thus express $\dot{\bm{x}}_{j}$ as a linear combination of the derivatives of $\theta_{i}$ and $\phi_{i}$ for $i=1,\ldots,j-1$, in addition to including the time derivative of the base point $\bm{x}_{1}$. Hence we may write

[TABLE]

where $Q$ is a $3(N+1)\times 3(N+1)$ matrix and $x_{1,x},x_{1,y,}x_{1,z}$ are the components of $\bm{x}_{j}$ in the basis $\{\bm{e}_{x},\bm{e}_{y},\bm{e}_{z}\}$. Explicitly, $Q$ may be constructed via

[TABLE]

where the matrices $Q_{k1}$ are of dimension $(N+1)\times 3$, with $Q_{k2}$ and $Q_{k3}$ being of dimension $(N+1)\times N$, for $k=1,2,3$. In the definition of $Q$, the subscript $P$ denotes that the $i$th row of $\tilde{Q}$ is permuted to the $P(i)$th row of $Q$, where

[TABLE]

This permutation of $\tilde{Q}$ allows us to define the sub-blocks simply, given explicitly as

[TABLE]

Further, in this parameterisation we may readily relate the local rate of rotation about $\bm{d}_{3}^{i}$, denoted $\omega_{i}$ as in Eq. 24, to $\theta,\phi,\psi$ and their time derivatives. Explicitly, this relationship is $\omega=\cos({\theta})\dot{\phi}+\dot{\psi}$, and is notably linear in the derivatives of the Euler angles. Thus, we form the composite matrix

[TABLE]

where $I_{N}$ is the $N\times N$ identity matrix and the $N\times N$ matrix $C$ has diagonal elements $C_{i}=\cos({\theta_{i}})$ for $i=1,\ldots,N$, with all other elements zero. The upper block, $Q$, maps the parameterisation into the laboratory frame, whilst the lower blocks convert between the parameterisation and the local rate of rotation about $\bm{d}_{3}$ as written in director basis. The $(4N+3)\times 3(N+1)$ matrix $\mathcal{Q}$ now encodes the expressions of velocities and rotation rates in terms of the parameterisation, via

[TABLE]

noting that the representation of $\bm{\omega}$ is relative to the director basis, whilst the representation of $\dot{\bm{X}}$ is relative to the basis of the laboratory frame.

Having constructed $\mathcal{Q}$, we now combine Eqs. 25 and 36 to give

[TABLE]

noting in particular that the matrix $\mathcal{B}\mathcal{A}\mathcal{Q}$ is square and of dimension $(3N+3)\times(3N+3)$. Naively, this system of ordinary differential equations can be readily solved numerically to give the evolution of the filament in the surrounding fluid. However, the use of a single parameterisation to describe the filament will in general lead to degeneracy of the linear system and ill-defined derivatives in both space and time, issues which we explore and resolve numerically in the subsequent sections.

II.5 Coordinate singularities

Consider a straight filament aligned with the $\bm{e}_{z}$ axis, with each of the $\bm{d}_{3}^{i}=[0,0,1]^{\top}$ written in the laboratory frame. For this filament $\theta_{i}=0$ for all $i$, whilst the $\phi_{i}$ are undetermined, arbitrary and notably need not be the same on each segment. Were we to attempt to formulate and solve the linear system of Eq. 37, both $\phi$ and its derivatives would be ill-defined, and correspondingly we would be unable to solve the system for the filament dynamics, which are physically trivial in this particular setup. In more generality, if a filament were to have any segment pass through one of the poles $\theta=\{0,\pi\}$ of this coordinate system, $\phi$ would be undetermined on the segment and arbitrary, with attempts to solve our parameterised system of ordinary differential equations failing. Further, were a segment to pass close to but not through a pole, time derivatives of $\phi$ would necessarily become large, with $\phi$ well-defined but varyingly rapidly as the segment moves close to the pole of the coordinate system. These large derivatives would artificially introduce additional stiffness to the elastohydrodynamical problem, inherent only to the parameterisation and not the underlying physics. This problem is well known for Euler angle parameterisations, and is commonly referred to as the ‘gimbal lock’. Analogous issues with arclength derivatives occur when considering neighbouring segments, with the value of $\phi$ varying rapidly and artificially between segments that reside near the pole of the coordinate system. In this latter case however, our formulation of the elastohydrodynamical problem circumvents the need for evaluation of $\phi_{s}$, instead considering only derivatives of the smooth quantities $\bm{d}_{\alpha}$, though we are not able to resolve issues with temporal derivatives in the same way.

In order to avoid the numerical and theoretical problems associated with singular points in the filament parameterisation, we exploit the finiteness of the set of angles $\theta_{i}$ along with the independence of the underlying elastohydrodynamical problem from the parameterisation. Throughout this work we have assumed a fixed laboratory frame with basis $\{\bm{e}_{x},\bm{e}_{y},\bm{e}_{z}\}$, present only so that vector quantities may be written component-wise for convenience. Our choice of such a basis is arbitrary, with the physical problem of filament motion being independent of our selection of particular basis vectors. It is with respect to this basis that we have defined the Euler angles $\theta,\phi,\psi$, from which the aforementioned coordinate singularities appear if any of the $\theta_{i}$ approach zero or $\pi$. Thus, if one makes a choice of basis $\{\bm{e}_{x}^{\star},\bm{e}_{y}^{\star},\bm{e}_{z}^{\star}\}$ such that the corresponding Euler angles $\theta_{i}^{\star}$ are some $\delta$-neighbourhood away from the poles of the new parameterisation, the system of ordinary differential equations given in Eq. 37 may be readily solved, at least initially. Should the solution in the new coordinate system approach one of the new poles $\theta^{\star}=0,\pi$, a new basis can again be chosen, and this process iterated until the filament motion has been captured over a desired interval.

We note that for sufficiently small $\delta>0$ such a choice of basis $\{\bm{e}_{x}^{\star},\bm{e}_{y}^{\star},\bm{e}_{z}^{\star}\}$ necessarily exists due to the finiteness of the set of $\theta_{i}$, with $\delta$ in practice able to be sufficiently large so as to limit the effects of coordinate singularities. Thus, subject to reasonable assumptions of smoothness of the filament position $\bm{x}$, such a process of repeatedly changing basis when necessary will prevent issues associated with the parameterisation described above, and will in practice enable the efficient simulation of filament motion without introducing significant artificial stiffness or singularities.

III Implementation, verification and extensions

III.1 Selecting a new basis

Initially choosing an arbitrary laboratory basis $\{\bm{e}_{x},\bm{e}_{y},\bm{e}_{z}\}$, the above formulation is implemented in MATLAB®, with the system of ordinary differential equations (ODEs) of Eq. 37 being solved using the inbuilt stiff ODE solver ode15s, as described in detail by Shampine and Reichelt [27]. This standard variable-step, variable-order solver allows for configurable error tolerances, typically set here at $10^{-5}$ for absolute error and $10^{-4}$ for relative error, in general significantly below the error associated with the piecewise-linear filament discretisation. Derivatives with respect to arclength are approximated with fourth order finite differences, with the resulting dynamics insensitive to this choice of scheme. Initially and at each timestep, the values of $\theta_{i}$ are checked to determine if they are within $\delta$ of a coordinate singularity, typically with $\delta=\pi/50$. Should the parameterisation be approaching a singularity, a new basis is chosen and the problem recast in this basis.

A natural method of selecting a new basis is perhaps to choose one uniformly at random. Indeed, by considering the worst-case scenario of the $N$ tuples $(\theta_{i},\phi_{i})$ uniformly and disjointly covering the surface of the unit sphere, which together $\theta$ and $\phi$ parameterise, the probability that any random basis results in a scenario with $\min_{i}{\{\theta_{i},\pi-\theta_{i}\}}<\delta$ is given by $2N\sin^{2}(\delta/2)$, a consequence of elementary geometry. With this quantity being significantly less than unity for a wide range of $N$ with $\delta$ large enough to avoid severe artificial numerical stiffness, as discussed above, a practical implementation for the simulation of filament elastohydrodynamics as formulated above may simply select a new basis randomly, repeating until a suitable basis is found. With $\delta=\pi/50$ and $N=50$, the probability of rejecting a candidate new basis is bounded above by $10\%$, thus in practice one should expect to find an appropriate basis within few iterations of the proposed procedure.

Alternatively, and as we will do throughout this work, one may instead proceed in a deterministic manner, selecting a near-optimal basis from knowledge of the existing parameterisation. Given the set of parameters $\theta_{i}$ and $\phi_{i}$, we may choose a $\hat{\theta}\in[0,\pi]$ and $\hat{\phi}\in[0,2\pi)$ so as to maximise the distance of $(\hat{\theta},\hat{\phi})$ from each of the $(\theta_{i},\phi_{i})$ and their antipodes. In practice, an approximate solution to this problem is attained by selecting $(\hat{\theta},\hat{\phi})$ from a selection of preset test points in order to maximise the distance from the $(\theta_{i},\phi_{i})$, where distance is measured on the surface of the unit sphere that $\theta$ and $\phi$ naturally parameterise, as shown in Fig. 1. It should be noted that this process of selection impacts negligibly on computational efficiency with 10,000 test points. With these choices of $\hat{\theta}$ and $\hat{\phi}$, we form a new basis by mapping the original basis vector $\bm{e}_{z}$ to the vector $\bm{e}_{z}^{\star}$, given explicitly by

[TABLE]

Choosing the other members of the orthonormal basis $\bm{e}_{x}^{\star},\bm{e}_{y}^{\star}$ arbitrarily, expressed in this new basis the accompanying filament parameterisation will be removed from any coordinate singularities by construction, as exemplified in Fig. 1(b).

III.2 Validation

In what follows we validate the presented methodology against known filament behaviours and a sample three dimensional simulation via an existing methodology. Initial configurations and parameter values for each can be found in the Appendix A, with behaviours qualitatively independent of these parameter choices and filament setups.

III.2.1 Relaxation of a planar filament

Noting there is no analytical test solution for the dynamics of a fully 3D Kirchhoff rod in a viscous fluid, to the best of our knowledge, we consider validations by comparison with numerical studies in the literature, though we additionally utilise invariance of the centre of mass as a gold standard below, where applicable. Firstly, we validate the presented approach by considering the problem of filament relaxation in two dimensions, a natural and well-studied subset of the three dimensional framework. We consider the simple case of a symmetric curved filament in the $\bm{e}_{x}\bm{e}_{y}$ plane, which will provide a test of symmetry preservation, integrated moment balance, and hydrodynamics via qualitative comparisons with the earlier works of Moreau et al. [14] and Hall-Mcnair et al. [15]. Simulating the relaxation of such a filament to a straight equilibrium condition with $N=40$ segments takes less than $2\text{\,}\mathrm{s}$ on modest hardware (Intel® Core™ i7-6920HQ CPU), from which we immediately see retained the computational efficiency of the framework of Moreau et al. which this work generalises. Present throughout the computed motion is the left-right symmetry of the initial condition, with the filament shape evolving smoothly even with a small number of segments, as shown in Fig. 2. Further, the centre of mass, which in exact calculation would be fixed in space due to the overall force-free condition on the filament, is captured numerically with errors on the order of $10^{-3}L$, demonstrating very good quantitative satisfaction of this condition. Of note, in Fig. 2b we have verified that this error is of the same order of magnitude as that generated by the two-dimensional methodology of Walker et al. [16].

III.2.2 Planar bending of a filament in shear flow

Further, whilst the above is reassuring and serves to validate a subset of the implementation, we note from Eqs. 26 and 27 that motion cast in the $\bm{e}_{x}\bm{e}_{y}$ plane of the laboratory frame may often render the evolution of one of the directors $\bm{d}_{1},\bm{d}_{2}$ trivial. In order to avoid this we may consider planar problems in slightly more generality, posing a problem that is planar though not aligned with the $\bm{e}_{x}\bm{e}_{y}$ plane. As an exemplar such problem we attempt to recreate a typical but complex behaviour of a flexible filament in a shear flow, that of the J-shape and U-turn [28], aligning both the filament and the background flow in a plane spanned by $\bm{e}_{y}$ and $\bm{e}_{x}+\bm{e}_{z}$. In order to ensure the absence of alignment of the parameterisation with the $\bm{e}_{x}\bm{e}_{y}$ plane, we disable the adaptive system of basis selection for the purpose of this example, and simulate the motion of a filament in a shear flow. The background flow with velocity $\bm{u}_{b}$ and vorticity $\bm{\Omega}$ is incorporated into the framework via the transformation $\bm{u}\mapsto\bm{u}-\bm{u}_{b}$, $\omega_{i}\mapsto\omega_{i}-\bm{\Omega}\cdot\bm{d}_{3}^{i}/2$, yielding the modified system

[TABLE]

for a vector $\hat{\bm{U}_{b}}$ of background flows and vorticities evaluated at segment endpoints and midpoints, respectively. Details of the flow field and initial setup are given in the Appendix A.

Adopting the timescale $T$ to be the inverse of the shearing rate of the flow, we consider a parameter regime in which one would expect to see to formation of J-shape and subsequently a U-turn, defined by their characteristic morphologies in Liu et al. [28, Figure 1, Movies S5, S6] and from which an appropriate parameter choice is obtained. Notably, the impact of thermal noise perturbations is not considered here, in contrast to Liu et al. [28], preventing a quantitative comparison. In Fig. 3 we present the initial, J-shape, and U-turn configurations of the filament as simulated via the proposed methodology, with computation requiring approximately $20\text{\,}\mathrm{s}$ with $N=50$ segments. The simulated filament shapes are in qualitative agreement with those shown in Liu et al. [28, Figure 1], and serve as further validation of the coarse-grained methodology. Notably, enabling the described method of basis selection effectively casts this problem in the $\bm{e}_{x}\bm{e}_{y}$ plane, affording a twofold increase in computational efficiency and highlighting the benefits of adaptive reparameterisation.

III.2.3 Relaxation of a non-planar filament

Finally, we consider truly non-planar relaxation of filaments in three dimensions. Typical simulations of such a relaxation with $N=50$ segments have a runtime of approximately 10$\mathrm{s}$ on the modest hardware described above, often requiring at most one choice of basis though naturally problem dependent, and provide reasonable accuracy. Thus, even when considering inherently three-dimensional problems we see retained in this methodology the low computational cost of the formulation of Moreau et al. [14], representing significant improvements in computational efficiency over recent studies in three dimensions [19, 20, 21]. This is particularly evident when directly comparing the presented coarse-grained methodology with the results of Ishimoto and Gaffney [21], considering in this case the relaxation of a helical configuration to a straight equilibrium. A side-by-side comparison of the relaxation dynamics as computed by the proposed methodology and that presented in Ishimoto and Gaffney is shown in Fig. 4, noting that the work of Ishimoto and Gaffney [21] considers an actively driven nearly inextensible filament, of which relaxation dynamics are a natural subset. Figure 4 highlights good agreement between methodologies that is in line with the level of accuracy typically afforded by resistive force theories used here, recalling errors logarithmic in the filament aspect ratio. Figure 4e shows a quantitative evaluation of the computed solutions, with the deviation of the filament centre of mass from the initial condition shown as solid black curves for each methodology. With the force-free condition implying that the filament centre of mass should not move throughout the relaxation dynamics, this measured deviation serves as an assessment of the accuracy of both frameworks, with each exhibiting variation only on the order of $10^{-2}L$. Also shown in Fig. 4e is a dimensional measure of the difference between the two computed solutions, here denoted $E(t)$ and defined as the non-negative root of

[TABLE]

where $\bm{x}_{P}(s,t)$ and $\bm{x}_{IG}(s,t)$ denote the filament centreline as computed by the proposed methodology and that used by Ishimoto and Gaffney, respectively. Numerically approximating this integral with quadrature and noting that this error is consistently on the order of $10^{-2}L$ throughout the relaxation, we see evidenced good quantitative agreement between the two frameworks, thus validating the proposed methodology. We also remark that the solution of Ishimoto and Gaffney [21] does not perfectly satisfy filament inextensibility, with variation in total length of approximately 1%, which may have some impact on the computed dynamics. Thus, when computing $E(t)$ as defined above, we treat $s\in[0,L]$ as a material parameter, with the differences in position $\bm{x}_{P}(s,t)-\bm{x}_{IG}(s,t)$ therefore capturing discrepancies between the simulated locations of the material point with undeformed arclength $s$ at time $t$, with $s$ not necessarily equal to the deformed arclength in the solution of Ishimoto and Gaffney [21].

In contrast to the agreement between solutions, there is a stark difference between the associated time required for computation. Taking $N=100$ segments and computing until relaxation, the coarse-grained framework calculated the solution in approximately $30\text{\,}\mathrm{s}$ on personal computing hardware with ODE error tolerances of $10^{-5}$, whilst the computations utilising the implementation of Ishimoto and Gaffney required $2.5\text{\,}\mathrm{h}$ on a high performance computing cluster. A thorough investigation of the time required for computation with the presented methodology for various choices of the parameters $\epsilon$ and $N$ is showcased in Fig. 5, from which we note the remarkable performance of this simple implementation across parameter regimes, with the walltime naturally dependent on the discretisation parameter $N$.

III.3 Model extensions

III.3.1 Intrinsic curvature

In order to showcase the versatility of the presented framework, we demonstrate its simple extension to filaments with non-zero intrinsic or reference curvature, which can exhibit complex buckling behaviours [29]. Recalling the constitutive relation of Eq. 8, the effect of an intrinsic curvature $\bm{\kappa}^{0}$ is to alter the bending moments, yielding the modified constitutive relation

[TABLE]

where we have written $\bm{\kappa}^{0}=\sum_{\alpha}\kappa_{\alpha}^{0}\bm{d}_{\alpha}$ in the local director basis. Notably, the reference curvature can plausibly depend on a variety of quantities, including time, arclength, and spatial position, with an example being the modelling of an internally driven filament by a time-dependent intrinsic curvature [18, 19].

Practically, the inclusion of such an intrinsic curvature amounts to a simple subtraction of the reference curvature from the computed components of $\bm{\kappa}$ at each instant, with $\bm{R}$ as given in Eq. 20 being modified accordingly. Doing so, we simulate the relaxation of a straight filament with a non-zero intrinsic curvature, with the reference curvature explicitly given by $\bm{\kappa}^{0}=\pi\bm{d}_{2}+2\pi\bm{d}_{3}$ corresponding to a helical configuration. As shown in Fig. 6, the filament indeed relaxes to a helix, with the walltime of this simulation being $15\text{\,}\mathrm{s}$ on a laptop computer, having taken $\epsilon=10^{-2}L$ and $N=70$ segments, noting that the filament shape has been sufficiently resolved with this discretisation.

III.3.2 Clamped and internally driven filaments

Finally, we consider the simple extensions to both clamped filaments and to those with actively generated internal moments, the latter being akin to the active beating of biological cilia and flagella [30]. For time and arclength-dependent active moment density $\bm{m}^{a}$, its inclusion into the presented framework acts to modify the moment balance of Eq. 2 to

[TABLE]

Repeating the integration of the pointwise moment balance from $s=s_{i}$ to $s=L$ as in Section II leads to a modified form of Eq. 11, explicitly given as

[TABLE]

where the integrated active moment density is written as

[TABLE]

For a given active moment density, assumed to be integrable, $\bm{I}_{i}^{a}$ may be readily computed either analytically or numerically, with its components in the local $\bm{d}_{\alpha}^{i}$ direction then modifying $\bm{R}$ from Eq. 20 accordingly.

Clamping the filament at the base is somewhat simpler, in that the overall force and moment balance conditions on the filament are merely replaced by enforcing no motion or rotation at the base. These conditions may be stated concisely as

[TABLE]

supplanting the force and moment-free conditions of Eqs. 5 and 6, having taken $s=0$ in the latter. Implementing these minor modifications, as an example we specify a travelling wave of internal moment given by $\bm{m}^{a}=5\sin(s-t)\bm{d}_{1}+5\cos(s-t)\bm{d}_{2}$ and simulate the active motion of a clamped filament, taking $N=50$ and $\epsilon=10^{-2}L$. Snapshots of this eventually periodic motion are shown in Fig. 7, with the motion simulated up until $t=8\pi$ from a straight initial configuration and with a walltime of around $8\text{\,}\mathrm{s}$.

IV Discussion

In this work we have seen that the motion of inextensible unshearable filaments in three dimensions can be concisely described by a coarse-grained framework, building upon the principles of Moreau et al. [14] in order to minimise numerical stiffness associated with the governing equations of elastohydrodynamics. This representation was readily implemented via an Euler angle parameterisation, with adaptive basis selection and reparameterisation overcoming coordinate singularities associated with a fixed representation of the dynamics. We have further demonstrated the efficacy of the proposed deterministic method for basis selection by explicit simulation of non-planar filament dynamics, which is able to afford reductions in numerical stiffness even when separated from singularities of the parameterisation. The presented framework retains the flexibility and extensibility of the formulations of Moreau et al. [14], Hall-Mcnair et al. [15], Walker et al. [16], with background flows, active moment generation, body forces and other effects or constraints being simple to include in this representation. The simplicity of these possible extensions speaks to the broad utility of the proposed approach, with potential for use in the simulation of both single and multiple filamentous bodies in fluid under a wide variety of circumstances and conditions.

In formulating our methodology we have made the simplifying assumption of coupling fluid dynamics to forces via resistive force theory, which is well-known to incur errors logarithmic in the filament aspect ratio, though variations remain in widespread use [31, 30, 32, 33, 34, 35, 14]. Resistive force theories additionally suffer from locality, in that portions of the filament do not directly interact with one another through the fluid. A natural development of the presented approach would therefore be the inclusion of non-local hydrodynamics, perhaps via the regularised Stokeslet segment methodology of Cortez [36] as included in the work of Walker et al. [16], or lightweight singular slender body theories such as those of Johnson [37], Tornberg and Shelley [4], Walker et al. [38]. Such improvements may also include the consideration of confined geometries, with motion in a half space of particular pertinence to typical microscopy of flagellated organisms. Despite the many possible directions for hydrodynamic refinement, we note that the linearity of Stokes flows necessitates that any relations between forces and velocities be linear, with explicit formulations simply giving rise to modified linear operators $\mathcal{A}$ that may be readily inserted into the described framework.

In the wider context of methods for soft filament simulation, the scope of which is illustrated by the generality of the approach of Gazzola et al. [39], the proposed framework enables rapid simulation of the subclass of purely inextensible filaments. This modeling assumption has been demonstrated to be a valid approximation of nearly inextensible filaments in a variety of contexts [23], affording greatly improved computational efficiency over previous approaches that we have compounded here, albeit with reduced hydrodynamic fidelity. However, we expect that the addition of improved hydrodynamics will have minimal impact on the computational efficiency of the presented methodology, with this efficiency not being derived from our use of simple resistive force theory, as noted by Walker et al. [16], Hall-Mcnair et al. [15] for their non-local refinements of the 2D theory of Moreau et al. [14].

In summary, we have presented, verified and exemplified a novel framework for the rapid simulation of inextensible, unshearable filaments in a viscous fluid at zero Reynolds number. Despite the improved generality of this methodology over existing two-dimensional approaches, we have retained the computational efficiency and simplicity of the work of Moreau et al. [14], affording significant extensibility and thus facilitating a vast range of previously unrealisable biological and biophysical studies into filament dynamics on the microscale.

Acknowledgements

We are grateful to Prof. Derek Moulton for discussions on elastic filaments, and to Prof. David Smith for discussions on basis rotation. B.J.W. is supported by the UK Engineering and Physical Sciences Research Council (EPSRC), grant EP/N509711/1. K.I. acknowledges MEXT Leading Initiative for Excellent Young Researchers (LEADER), JSPS KAKENHI for Young Researcher (JP18K13456), and JST, PRESTO Grant Number JPMJPR1921, Japan. Elements of the simulations were performed using the cluster computing system within the Research Institute for Mathematical Sciences (RIMS), Kyoto University.

The computer code used and generated in this work is freely available from https://gitlab.com/bjwalker/3d-filaments.git.

Appendix A Parameters and initial conditions

We nondimensionalise the system Eq. 37 as in Walker et al. [16], resulting in a dimensionless system of the form

[TABLE]

for timescale $T$ and elastohydrodynamic number $E_{h}$, where dimensionless counterparts are given by

[TABLE]

having multiplied the force balance equations by $\Delta s$ to unify dimensions. All examples begin with the filament base, $\bm{x}_{1}$, coincident with the origin of the laboratory frame, and we set Poisson’s ratio to zero, i.e.  $\sigma=0$ throughout though note a lack of sensitivity of results to this choice.

A.1 Relaxation of a planar filament

In simulating the relaxation of a filament in the $\bm{e}_{x}\bm{e}_{y}$ plane, we impose the initial filament shape as

[TABLE]

with respect to standard laboratory Euler angles. Filament motion is simulated with $E_{h}\approx$1.6\text{\times}{10}^{5}\text{,}$$, though this choice is arbitrary given the invariance of the dynamics to rescalings in time (assuming that the rescaling is not so extreme as to break the inertialess assumption).

A.2 Planar bending of a filament in shear flow

In order to generate the characteristic behaviours of the J-shape and U-turn we initialise the filament via

[TABLE]

We align a background shear flow in the same plane as the filament, proportional in strength to the coordinate in the $\bm{e}_{y}$ direction, denoted $y$. Explicitly, this flow $\bm{u}_{b}$ is given in the laboratory frame in dimensionless form by

[TABLE]

having taken the timescale $T$ to be the inverse of the dimensional shear rate. Simulation proceeds with $E_{h}\approx$4.7\text{\times}{10}^{5}\text{,}$$, consistent with the regime found in Liu et al. [28], and we note that the tip of the filament initially curves into the oncoming background flow in $y<0$.

A.3 Relaxation of a non-planar filament

Taking $N=100$ segments, we impose the helical initial condition

[TABLE]

We simulate filament motion with $E_{h}\approx$3.1\text{\times}{10}^{4}\text{,}$$, with results insensitive to this choice. The parameters used in the implementation of Ishimoto and Gaffney [21] are as in their publication, with the image system for a plane wall accordingly removed and the actively generated torques set to zero to allow for filament relaxation in free space. In particular, whilst the filaments of Ishimoto and Gaffney [21] are extensible, the extensional modulus of these filaments is sufficiently high so as to provide near inextensibility in their results, enabling meaningful comparison.

A.4 Relaxation of a non-straight filament

Taking $N=70$ segments, we impose the straight initial condition

[TABLE]

The intrinsic curvature is specified as $\bm{\kappa}^{0}=\pi\bm{d}_{2}+2\pi\bm{d}_{3}$. We simulate filament motion with $E_{h}\approx$1.5\text{\times}{10}^{5}\text{,}$$, with results insensitive to this choice.

A.5 Active beating of a clamped filament

Taking $N=50$ segments, we impose the straight initial condition

[TABLE]

The active moment density is specified as $\bm{m}^{a}=5\sin(s-t)\bm{d}_{1}+5\cos(s-t)\bm{d}_{2}$. We simulate filament motion with $E_{h}=10^{3}$ up until $t=8\pi$.