Efﬁcient simulation of ﬁlament elastohydrodynamics in three dimensions

Fluid-structure simulations of slender inextensible ﬁlaments in a viscous ﬂuid are often plagued by numerical stiffness. Recent coarse-graining studies have reduced the computational requirements of simulating such systems, though they have thus far been limited to the motion of planar ﬁlaments. In this paper we extend such frameworks to ﬁlament motion in three dimensions, identifying and circumventing coordinate-system singularities introduced by ﬁlament parametrization via repeated changes of basis. The resulting methodology enables efﬁcient and rapid study of the motion of ﬂexible ﬁlaments in three dimensions, and is readily extensible to a wide range of problems, including ﬁlament motion in conﬁned geometries, large-scale active matter simulations, and the motility of mammalian spermatozoa. which we explore and resolve numerically in the subsequent sections.


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, utilized for driving flows and even locomotion, while 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, utilizing 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 [6,7,12]. 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 realizing 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 discretization of which yielded a relatively simple system of ordinary differential equations (ODEs) 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 nonlocal hydrodynamics in both infinite and semi-infinite domains, and applied to a variety of single and multifilament problems [15,16]. However, being confined to two dimensions limits the potential scope and applicability of these approaches, with three-dimensional (3D) 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 [3,17].
However, for nonplanar 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 explorations 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 utilizes a quarternion representation of filament orientation to parametrize 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], and Bouzarth et al. [22], each of which are equipped with nonlocal 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 were afforded by the latter in a range of biological and biophysical modeling 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 nonspecialized 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 paper 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 generalization of the Frenet triad and integration of the governing equations of elasticity. We will overcome fundamental issues with simple single parametrizations of a filament in three dimensions, presenting an effective computational approach utilizing adaptive reparametrization and basis selection. We will then validate the presented framework by consideration of three candidate test problems, simulating well-documented behaviors 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.

A. Equations of elasticity
We consider a slender, inextensible, unshearable filament in a viscous Newtonian fluid, with its centreline described by x(s), parametrized by arclength s ∈ [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 for contact force and torque denoted n, m, respectively, and where a subscript of s denotes differentiation with respect to arclength. The quantity 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ẋ, where here a dot denotes a time derivative. Similarly, τ is the torque per unit length applied on the fluid medium by the filament. While any standard boundary condition may be considered in the formalism below, throughout we impose zero force and torque at the filament ends so that 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 n(s) at the earliest opportunity; using the boundary condition n(L) = 0 and Eq. (1) we have This relation is also subject to the constraint that the boundary condition n(0) = 0 is satisfied, generating a global force balance constraint for the applied force per unit length, which we carry forward into the formalism below together with the elimination of 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 m(L) = 0 reveals the integrated moment balance: Given a right-handed orthonormal director basis {d 1 (s), d 2 (s), d 3 (s)}, generalizing the Frenet triad such that d 3 corresponds to the local filament tangent, following Nizette and Goriely [25] we define the twist vector κ by for α = 1, 2, 3. Writing κ = α κ α d α , for bending stiffness EI we use the Euler-Bernoulli constitutive relation of the Kirchhoff formalism to relate the contact torque m and the twist vector κ [24], via where σ is the Poisson ratio [25], assumed to be constant. With this constitutive relation the integrated moment balance equations in the d α directions are simply for α = 1, 2, 3 and where δ a,b denotes the Kronecker delta. 123103-3

B. Filament discretization
In discretizing 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 s, with segment end points having positions denoted by x 1 , . . . , x N+1 and the constraints of inextensibility and the absence of cross-section shear being satisfied inherently. The end points of the ith segment correspond to x i and x i+1 for i = 1, . . . , N, with the local tangent d 3 being constant on each segment and denoted d i 3 . In what follows we will consider a discretization of d 1 , d 2 such that they are also constant on each segment, and we denote these constants similarly as d i 1 , d i 2 . Writing s i for the constant arclength associated with each material point x i , we apply Eq. (9) at each of the s i for i = 1, . . . , N, splitting the integral at the segment end points to give for α = 1, 2, 3. On the jth segment, x may be written as is given by η = (s − s j )/ s. Additionally discretizing the force per unit length as a continuous piecewise-linear function, with η as above we have where we write f j = f (s j ). Substitution of these parametrizations into Eq. (10) and subsequent integration yields, after simplification, where the integral contribution of the force and torque densities are denoted I f i and I τ i , respectively. With this discretization I f i has reduced to 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 discretization for I τ i . Hence, taking the piecewise constant discretization τ = τ j on the jth segment, we have the simple expression 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 f j and the τ j . Similarly, with this piecewise-linear force discretization the integrated force balance of Eq. (1) simply reads We ] for components f j,x , f j,y, f j,z of f j with respect to some fixed laboratory frame with basis {e x , e y , e z }, and similarly 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 and correspond to the force balance of Eq. (14), with the remaining 3N columns zero. The remaining rows of B encode the moment balance of Eq. (11) as expanded in Eq. (12), organized in triples such that B 3(i−1)+3+α projects the ith moment balance equation 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 with the latter expression readily transcribed as a linear operator acting on the f j for j = i, . . . , N. Analogously, we have from which a linear operator acting on the τ j for j = i, . . . , N can be constructed. Accordingly, the (3N + 3)-vector R is given by 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 B and R are well defined for a general filament in three dimensions, given the local directors d 1 and d 2 and computing the components of the twist vector as In terms of the discretized 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 κ 1 (0) = κ 2 (0) = κ 3 (0) = 0.

C. Coupling hydrodynamics
We now relate the force density f acting on the fluid to the velocity of each segment end point, utilizing the commonly applied method of resistive force theory as introduced by Hancock [2] and 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 = 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 123103-5 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 utilize the expressions of Gray and Hancock [10], with where μ is the medium viscosity, noting the relation C n = 2C t . We approximate the local filament tangent at the segment end point x i as the average of d i−1 3 and d i 3 for i = 2, . . . , N, with the tangent for i = 1 and N + 1 simply being taken as d 1 3 and d N 3 , respectively. By linearity, and again assuming a piecewise-linear force density along segments, we may write the coupling of translational kinematics to hydrodynamics asẊ where A is a square matrix of dimension 3(N + 1) × 3(N + 1) and is a function only of the segment end points x i . Of dimension 3(N + 1), the vectorẊ corresponds to the linear velocities of the segment end points, and is constructed analogously to F with respect to the laboratory frame. This relation results from the application of the no-slip condition at the segment end points, coupling the filament to the surrounding fluid.
In order to relate the rate of rotation of each segment to the viscous torque τ 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 ith segment with the rotation ω i about its local tangent d i 3 via the relation of Chwang and Wu [26]: (24) and in particular the 2 scaling entails that the torque per unit length contributions are relatively small. Here we recall that μ is the viscosity of the fluid medium, and is the radius of the filament. We may write this relation as a linear operator on ω = [ω 1 , . . . , ω N ] , written simply as T =Ãω. This crude approximation may readily be substituted for nonlocal hydrodynamics via the method of regularized Stokeslet segments, which will likely be a topic of future work. Similarly, nonlocal hydrodynamics may be utilized in place of Eq. (23), as used for two-dimensional (2D) 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 where A is assumed to be invertible and A is defined to be a block matrix of dimension (6N + 3) × (4N + 3) with nonzero blocks A −1 andÃ.

D. Parametrization
We may parametrize the tangents d i 3 on each linear segment by the Euler angles θ i ∈ [0, π], φ i ∈ (−π, π], ψ i ∈ (−π, π] for i = 1, . . . , N [24]. With this parametrization we may make a choice of d 1 and d 2 , taking here the three orthonormal vectors to be written with respect to the laboratory frame and where s θ ≡ sin θ i , c θ ≡ cos θ i , and analogously for s φ , c φ , s ψ , and c ψ . From the directors we recover As the discretized filament is piecewise linear, for j = 1, . . . , N + 1 we may write With d i 3 parametrized as above, we can thus expressẋ j as a linear combination of the derivatives of θ i and φ i for i = 1, . . . , j − 1, in addition to including the time derivative of the base point x 1 . Hence we may write where Q is a 3(N + 1) × 3(N + 1) matrix and x 1,x , x 1,y, x 1,z are the components of x j in the basis {e x , e y , e z }. Explicitly, Q may be constructed viã where the matrices Q k1 are of dimension (N + 1) × 3, with Q k2 and Q k3 being of dimension (N + 1) × N, for k = 1, 2, 3. In the definition of Q, the subscript P denotes that the ith row ofQ is permuted to the P(i)th row of Q, where This permutation ofQ allows us to define the sub-blocks simply, given explicitly as Further, in this parametrization we may readily relate the local rate of rotation about d i 3 , denoted ω i as in Eq. (24), to θ, φ, ψ and their time derivatives. Explicitly, this relationship is ω = cos(θ )φ +ψ, and is notably linear in the derivatives of the Euler angles. Thus, we form the composite matrix where I N is the N × N identity matrix and the N × N matrix C has diagonal elements C i = cos(θ i ) for i = 1, . . . , N, with all other elements zero. The upper block, Q, maps the parametrization into the laboratory frame, while the lower blocks convert between the parametrization and the local rate of rotation about d 3 as written in director basis. The (4N + 3) × 3(N + 1) matrix Q now encodes the expressions of velocities and rotation rates in terms of the parametrization, via noting that the representation of ω is relative to the director basis, while the representation ofẊ is relative to the basis of the laboratory frame.
Having constructed Q, we now combine Eqs. (25) and (36) to give noting in particular that the matrix BAQ is square and of dimension (3N + 3) × (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 parametrization 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.

E. Coordinate singularities
Consider a straight filament aligned with the e z axis, with each of the d i 3 = [0, 0, 1] written in the laboratory frame. For this filament θ i = 0 for all i, while the φ 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 φ 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 θ = {0, π} of this coordinate system, φ would be undetermined on the segment and arbitrary, with attempts to solve our parametrized system of ordinary differential equations failing. Further, were a segment to pass close to but not through a pole, time derivatives of φ would necessarily become large, with φ 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 parametrization and not the underlying physics. This problem is well known for Euler angle parametrizations, and is commonly referred to as the "gimbal lock." Analogous issues with arclength derivatives occur when considering neighboring segments, with the value of φ 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 φ s , instead considering only derivatives of the smooth quantities d α , 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 parametrization, we exploit the finiteness of the set of angles θ i along with the independence of the underlying elastohydrodynamical problem from the parametrization. Throughout this paper we have assumed a fixed laboratory frame with basis {e x , e y , e z }, present only so that vector quantities may be written componentwise 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 θ, φ, ψ, from which the aforementioned coordinate singularities appear if any of the θ i approach zero or π . Thus, if one makes a choice of basis {e x , e y , e z } such that the corresponding Euler angles θ i are some δ neighborhood away from the poles of the new parametrization, 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 θ = 0, π, 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 δ > 0 such a choice of basis {e x , e y , e z } necessarily exists due to the finiteness of the set of θ i , with δ 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 x, such a process of repeatedly changing basis when necessary will prevent issues associated with the parametrization described above, and will in practice enable the efficient simulation of filament motion without introducing significant artificial stiffness or singularities.

A. Selecting a new basis
Initially choosing an arbitrary laboratory basis {e x , e y , e z }, the above formulation is implemented in MATLAB, with the system of 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, variableorder 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 discretization. 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 time step, the values of θ i are checked to determine if they are within δ of a coordinate singularity, typically with δ = π/50. Should the parametrization 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 (θ i , φ i ) uniformly and disjointly covering the surface of the unit sphere, which together θ and φ parametrize, the probability that any random basis results in a scenario with min i {θ i , π − θ i } < δ is given by 2N sin 2 (δ/2), a consequence of elementary geometry. With this quantity being significantly less than unity for a wide range of N with δ 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 δ = π/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 paper, one may instead proceed in a deterministic manner, selecting a near-optimal basis from knowledge of the existing parametrization. Given the set of parameters θ i and φ i , we may choose aθ ∈ [0, π] andφ ∈ [0, 2π ) so as to maximize the distance of (θ,φ) from each of the (θ i , φ i ) and their antipodes. In practice, an approximate solution to this problem is attained by selecting (θ,φ) from a selection of preset test points in order to maximize the distance from the (θ i , φ i ), where distance is measured on the surface of the unit sphere that θ and φ naturally parametrize, 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θ andφ, we (a) Before a change of basis and subsequent reparametrization, we see that these points are located near to the θ = 0 and π axes, which are shown as black vertical lines. A new potential location for the θ = 0 axis is shown in orange, selected so as to maximize the distance from the (θ i , φ i ) and their antipodes. (b) Following reparametrization, the points and antipodes are located maximally away from the new axis. This example scenario corresponds to a helical filament initially with θ i = π/6 for i = 1, . . . , 50.
Choosing the other members of the orthonormal basis e x , e y arbitrarily, expressed in this new basis the accompanying filament parametrization will be removed from any coordinate singularities by construction, as exemplified in Fig. 1(b).

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

Relaxation of a planar filament
Noting that 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 utilize invariance of the center of mass as a gold standard below, where applicable. First, 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 e x 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 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. that this paper generalizes. Present throughout the computed motion is the left-right symmetry of the  [14], Hall-McNair et al. [15], Walker et al. [16]. (b) Translation of the center of mass throughout the motion, analytically zero, is captured numerically with errors on the order of 10 −3 L by the proposed methodology, notably the same order of magnitude as that attained with the two-dimensional methodology of Walker et al. [16]. Axes x x and x y correspond to the unit vectors e x and e y , respectively. initial condition, with the filament shape evolving smoothly even with a small number of segments, as shown in Fig. 2. Further, the center 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. 2(b) 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].

Planar bending of a filament in shear flow
Further, while 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 e x e y plane of the laboratory frame may often render the evolution of one of the directors d 1 , 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 e x e y plane. As an exemplar such problem we attempt to recreate a typical but complex behavior 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 e y and e x + e z . In order to ensure the absence of alignment of the parametrization with the e x 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 u b and vorticity is incorporated into the framework via the transformation for a vectorÛ b of background flows and vorticities evaluated at segment end points and midpoints, respectively. Details of the flow field and initial setup are given in the appendix. 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 formation of J shape and subsequently a U turn, defined by their characteristic morphologies in Liu et al. (see Fig. 1 and Movies S5 and S6 in [28]) 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 Fig. 1 and Movies S5 and S6 in Liu et al. [28]). We note that the choosing of an improved basis for computation has been disabled for this example, and yields a twofold increase in computational efficiency if enabled. Arrows indicate the direction of the background shear flow. Axes x x , x y , x z correspond to the unit vectors e x , e y , e z , respectively. The plane containing the filament and the shear flow is shown in gray.
Notably, enabling the described method of basis selection effectively casts this problem in the e x e y plane, affording a twofold increase in computational efficiency and highlighting the benefits of adaptive reparametrization.

Relaxation of a nonplanar filament
Finally, we consider truly nonplanar relaxation of filaments in three dimensions. Typical simulations of such a relaxation with N = 50 segments have a runtime of approximately 10 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 4(e) shows a quantitative evaluation of the computed solutions, with the deviation of the filament center of mass from the initial condition shown as solid black curves for each methodology. With the force-free condition implying that the filament center 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. 4(e) is a dimensional measure of the difference between the two computed solutions, here denoted E (t ) and defined as the non-negative root  [21], in turn heavily based on the work of Olson et al. [19], shown from multiple perspectives at multiple time points. In both cases we observe relaxation from a nonplanar, helical configuration to a straight filament, and good agreement between the two computed solutions, in particular given the logarithmic accuracy of resistive force theories. (e) A quantitative comparison between the two frameworks, with solid lines showing the Euclidean distance of the filament center of mass from its initial location, analytically zero by the force-free condition though here on the order of 10 −2 L. The dashed curve quantifies the error between the computed solutions, denoted E , with the square of this error defined as E 2 = ||x P − x IG || 2 2 ds/L, where x P and x IG are the locations of the filament centerline as computed by the proposed methodology and that used by Ishimoto and Gaffney, respectively. Computation with the presented methodology took approximately 30 s on a modest laptop computer, in comparison to the multiple hours required on sophisticated cluster hardware for the methodology of Ishimoto and Gaffney. Here we have simulated filament motion with N = 100 segments. Axes x x , x y , x z correspond to the unit vectors e x , e y , e z , respectively. of where x P (s, t ) and 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 ∈ [0, L] as a material parameter, with the differences in position x P (s, t ) − 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 s on personal computing hardware with ODE error tolerances of 10 −5 , while the computations utilizing the implementation of Ishimoto and Gaffney required 2.5 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 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 discretization parameter N.

Intrinsic curvature
In order to showcase the versatility of the presented framework, we demonstrate its simple extension to filaments with nonzero intrinsic or reference curvature, which can exhibit complex buckling behaviors [29]. Recalling the constitutive relation of Eq. (8), the effect of an intrinsic curvature κ 0 is to alter the bending moments, yielding the modified constitutive relation where we have written κ 0 = α κ 0 α d α 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 modeling 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 κ at each instant, with R as given in Eq. (20) being modified accordingly. Doing so, we simulate the relaxation of a straight filament with a FIG. 6. The relaxation of a straight filament with nonzero intrinsic curvature to a helical configuration, with reference curvature specified as κ 0 = π d 2 + 2π d 3 . Having taken N = 70 segments and = 10 −2 L, we observe the smooth relaxation of the filament away from its initial straight configuration, computed in 15 s on a modest laptop computer. Axes x x , x y , x z correspond to the unit vectors e x , e y , e z , respectively. nonzero intrinsic curvature, with the reference curvature explicitly given by κ 0 = π d 2 + 2π 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 s on a laptop computer, having taken = 10 −2 L and N = 70 segments, noting that the filament shape has been sufficiently resolved with this discretization.

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 m a , its inclusion into the presented framework acts to modify the moment balance of Eq. (2) to Repeating the integration of the pointwise moment balance from s = s i to L as in Sec. II leads to a modified form of Eq. (11), explicitly given as where the integrated active moment density is written as For a given active moment density, assumed to be integrable, I a i may be readily computed either analytically or numerically, with its components in the local d i α direction then modifying 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 aṡ 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 m a = 5 sin(s − t )d 1 + 5 cos(s − t )d 2 and simulate the active motion of a 123103-15 FIG. 7. The 3D beating of a clamped filament driven by prescribed active internal moments. Having specified a traveling wave of out-of-phase sinusoidal active moments m a = 5 sin(s − t )d 1 + 5 cos(s − t )d 2 , an initially straight filament deforms to a periodic driven motion, with the tip of the filament following a circular path. Here we have taken = 10 −2 L and N = 50, noting that the filament shape has been resolved smoothly with this level of discretization. Computation required approximately 8 s on modest hardware, simulating up until t = 8π . Axes x x , x y , x z correspond to the unit vectors e x , e y , e z , respectively. clamped filament, taking N = 50 and = 10 −2 L. Snapshots of this eventually periodic motion are shown in Fig. 7, with the motion simulated up until t = 8π from a straight initial configuration and with a walltime of around 8 s.

IV. DISCUSSION
In this paper 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 minimize numerical stiffness associated with the governing equations of elastohydrodynamics. This representation was readily implemented via an Euler angle parametrization, with adaptive basis selection and reparametrization 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 nonplanar filament dynamics, which is able to afford reductions in numerical stiffness even when separated from singularities of the parametrization. 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 [14,[30][31][32][33][34][35]. 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 nonlocal hydrodynamics, perhaps via the regularized 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 Tornberg and Shelley [4], Johnson [37], Walker et al. [38]. Such improvements may also include the consideration of confined geometries, with motion in a half space being 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 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 Hall-McNair et al. [15] and Walker et al. [16] for their nonlocal refinements of the 2D theory of Moreau et al. [14].
In summary, we have presented, verified, and exemplified a 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 unrealizable biological and biophysical studies into filament dynamics on the microscale.
The computer code used and generated in this paper is freely available from [40].
with respect to standard laboratory Euler angles. Filament motion is simulated with E h ≈ 1.6 × 10 −5 , 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).

Planar bending of a filament in shear flow
In order to generate the characteristic behaviors of the J shape and U turn we initialize the filament via We align a background shear flow in the same plane as the filament, proportional in strength to the coordinate in the e y direction, denoted y. Explicitly, this flow u b is given in the laboratory frame in dimensionless form by having taken the timescale T to be the inverse of the dimensional shear rate. Simulation proceeds with E h ≈ 4.7 × 10 5 , 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.

Relaxation of a nonplanar filament
Taking N = 100 segments, we impose the helical initial condition We simulate filament motion with E h ≈ 3.1 × 10 4 , 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, while 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.

Relaxation of a nonstraight filament
Taking N = 70 segments, we impose the straight initial condition The intrinsic curvature is specified as κ 0 = π d 2 + 2π d 3 . We simulate filament motion with E h ≈ 1.5 × 10 5 , with results insensitive to this choice.
The active moment density is specified as m a = 5 sin(s − t )d 1 + 5 cos(s − t )d 2 . We simulate filament motion with E h = 10 3 up until t = 8π . 123103-18