TL;DR
This paper introduces a computational framework for simulating suspensions of passive and active filaments, addressing numerical stiffness and complex fluid-structure interactions with quaternion-based kinematics and implicit integration.
Contribution
The authors develop a novel simulation method using quaternions, implicit time-integration, and quasi-Newton methods to efficiently model filament suspensions with complex fluid interactions.
Findings
Framework effectively handles filament bending and twisting.
Method demonstrates good performance in test simulations.
Implementation available in MATLAB/Octave for further use.
Abstract
Flexible filaments and fibres are essential components of important complex fluids that appear in many biological and industrial settings. Direct simulations of these systems that capture the motion and deformation of many immersed filaments in suspension remain a formidable computational challenge due to the complex, coupled fluid--structure interactions of all filaments, the numerical stiffness associated with filament bending, and the various constraints that must be maintained as the filaments deform. In this paper, we address these challenges by describing filament kinematics using quaternions to resolve both bending and twisting, applying implicit time-integration to alleviate numerical stiffness, and using quasi-Newton methods to obtain solutions to the resulting system of nonlinear equations. In particular, we employ geometric time integration to ensure that the quaternions…
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Code & Models
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
\cleanlookdateon\xpretocmd
Methods for suspensions of passive and active filaments
Simon F Schoeller Contributed equally. Department of Mathematics, Imperial College London, London SW7 2AZ, UK
Adam K Townsend
Department of Mathematics, Imperial College London, London SW7 2AZ, UK
∗ Timothy A Westwood
Department of Mathematics, Imperial College London, London SW7 2AZ, UK
∗ Eric E Keaveny Email addresses: [email protected], [email protected],
[email protected], [email protected] Department of Mathematics, Imperial College London, London SW7 2AZ, UK
Abstract
Flexible filaments and fibres are essential components of important complex fluids that appear in many biological and industrial settings. Direct simulations of these systems that capture the motion and deformation of many immersed filaments in suspension remain a formidable computational challenge due to the complex, coupled fluid–structure interactions of all filaments, the numerical stiffness associated with filament bending, and the various constraints that must be maintained as the filaments deform. In this paper, we address these challenges by describing filament kinematics using quaternions to resolve both bending and twisting, applying implicit time-integration to alleviate numerical stiffness, and using quasi-Newton methods to obtain solutions to the resulting system of nonlinear equations. In particular, we employ geometric time integration to ensure that the quaternions remain unit as the filaments move. We also show that our framework can be used with a variety of models and methods, including matrix-free fast methods, that resolve low Reynolds number hydrodynamic interactions. We provide a series of tests and example simulations to demonstrate the performance and possible applications of our method. Finally, we provide a link to a MATLAB/Octave implementation of our framework that can be used to learn more about our approach and as a tool for filament simulation.
Contents
-
2.8.2 Differential-algebraic system and time integration in two dimensions
-
2.9.3 Differential-algebraic system and time integration in three dimensions
1 Introduction
*Microscopic flexible filaments and rigid fibres suspended in fluid arise in many industrial and biological applications, ranging from the wood pulp fibres found in paper production **(**Derakhshandeh et al., 2011; Stockie & Green, 1998; Ross & Klingenberg, *1997; *Pettersson et al., 2017; du Roure et al., 2017) to the flagella and cilia used by cells to swim, pump, or mix fluid (Lauga & Powers, *2009; *Elgeti et al., 2015; Brennen & Winet, 1977; Supatto & Vermot, *2011; *Smith et al., *2019; *Faubel et al., 2016; Elgeti & Gompper, 2013). In suspension, such as in polymeric fluids, their presence alters the bulk rheological properties of the medium leading to non-Newtonian responses such as shear thinning or viscoelasticity (Larson, 1999). When entangled or connected in networks, filaments and fibres form gels and disordered solids as is the case in important biological materials such as mucus **(**Hwang et al., 1969; Sheehan & Carlstedt, *1984; *Quraishi et al., 1998; Lai et al., 2009) and the extra-cellular matrix *(Baker et al., *2009; *Heck et al., 2017). Filaments can also be active biological elements that facilitate transport through self-deformation (Lauga & Powers, 2009; Brennen & Winet, 1977), or through growth and the motion of motor proteins along their lengths (Shelley, 2016).
*Simulation of filament or fibre motion requires solving a complex low Reynolds number fluid–structure interaction problem where the motion and deformation of all filaments are coupled through the surrounding fluid. As a result, computations involving filaments use many existing models and methods developed over the years to resolve the hydrodynamic aspects of the problem at various levels of approximation, ranging from drag-based resistive force theory **(Moreau et al., 2018; Yamamoto & Matsuoka, 1995) and point and regularised singularity methods **(Cortez, *2001; *Cortez et al., *2005; *Cosentino Lagomarsino et al., *2005; *Delmotte et al., *2015; *Olson et al., 2013; Smith, 2009; Smith et al., 2007) for single filament problems, to the immersed boundary method **(Peskin, 2002; Fauci & Peskin, 1988; Stockie & Green, 1998; Lim et al., 2008; Lim, 2010; Wiens & Stockie, 2015) and discretisations of the integral equations given by slender-body theory coupled with fast-summation techniques **(*Saintillan et al., 2005; Tornberg & Shelley, 2004; Gustavsson & Tornberg, 2009; Nazockdast et al., 2017; Tornberg & Gustavsson, 2006) for larger collections of flexible filaments or rigid fibres.
*While there are a variety of models and methods available to treat filament hydrodynamics, a similar catalogue of methods and algorithms for filament elasticity is less developed. Due to the higher-order derivatives associated with bending and twisting, the equations of filament motion are numerically stiff. In addition, conditions such as inextensibility give rise to constraints that must be satisfied as filaments move and deform. Despite this, explicit time integration is often used **(Peskin, *2002; *Li et al., *2013; *Olson et al., 2013; Simons et al., 2015) to advance the filaments in time, while constraints of inextensibility are mimicked through stiff springs **(Chelakkot et al., 2010; Stockie & Green, 1998), both of which tend to limit stable timestep sizes. There have been recent efforts to use implicit time integration schemes **(Moreau et al., *2018; *Hall-McNair et al., 2019), as well as treat the constraints directly using Lagrange multipliers **(Delmotte et al., 2015) or through a suitable choice of the filament degrees of freedom *(Moreau et al., *2018; *Hall-McNair et al., 2019). These studies, however, are limited in the number of filaments that they can treat as they rely on explicit knowledge of the filament mobility matrix and require dense linear algebra to solve for the unknown Lagrange multipliers, or to implicitly integrate the equations of motion. For specific hydrodynamic approaches intended for larger-scale simulation of many interacting filaments, similar techniques have been used to reduce computational times. For example, predictor–corrector schemes have been successfully used with the immersed boundary method (Wiens & Stockie, 2015), while recent work (Nazockdast et al., 2017), building from Tornberg & Shelley (2004), presents a complete approach utilising Lagrange multipliers and implicit time integration along with matrix-free hydrodynamics using the fast multipole method.
*In this study, we provide a comprehensive, computationally-scalable methodology for simulating filament dynamics that is not dependent on the specific method or model used to describe the hydrodynamic forces experienced by the filaments. Our approach only assumes that hydrodynamics forces and torques experienced by the discretised beam segments are linearly related to their translational and rotational velocities, consistent with low Reynolds number hydrodynamics. The methods that we present in this paper expand on those we have already used for planar simulations of thousands of interacting sperm cells (Schoeller & Keaveny, 2018) and the motion of undulatory swimmers through structured **(Majmudar et al., 2012) and unstructured (Kamal & Keaveny, 2018) complex environments. This paper presents the fully three-dimensional version of the methodology that takes advantage of unit quaternions to describe filament kinematics and geometric time integration schemes **(*Iserles et al., 2000; Faltinsen et al., 2001; Park & Chung, 2005) to advance the unit quaternions in time.
*Beginning from Kirchhoff rod theory (Kirchhoff, *1859; *Lim et al., 2008), we provide a description of the filament model and highlight the differences between deformations restricted to a plane and those that are fully three-dimensional. In particular, for fully three-dimensional deformations, we show how both filament bending and twisting can be captured through the use of unit quaternions to describe how the filament’s material frame varies with arclength. We employ Lagrange multipliers to obtain the constraint forces and moments necessary to link the material frame with the positions of points along the filament, and implicit time integration to handle numerical stiffness. By comparing Jacobian-free Newton–Krylov (Knoll & Keyes, 2004) and Broyden’s methods (Broyden, 1965; Kvaalen, 1991), we show how the resulting system of nonlinear equations, whose solution provides the updated generalised positions and Lagrange multipliers, can be solved most effectively using Broyden’s method with a suitable choice of approximate Jacobian based on a diagonal mobility matrix. Also novel to our approach is our application of geometric time integration to ensure the quaternions remain unit. We find these schemes help to guarantee robustness of the iterative solver for the nonlinear system of equations. We provide a number of tests of our method and explore several examples using both a hydrodynamic approach based on direct summation of the Rotne–Prager–Yamakawa (RPY) tensor **(Wajnryb et al., 2013) and matrix-free computations based on the force-coupling method (FCM) (Maxey & Patel, 2001; Lomholt & Maxey, 2003). In doing so, we show that our method provides an effective tool for studying a variety of applications involving filament dynamics including tethered filaments and cilia arrays, interacting undulatory swimmers, and sedimenting suspensions and clouds of filaments. We provide a basic MATLAB/Octave implementation *(Schoeller et al., 2019) of the method that can be used to reproduce our numerical results, applied to other problems of interest, or altered to couple with the user’s preferred hydrodynamic solver.
2 Filament model
We begin by describing the filament model for a single filament. In this work, a filament refers to an inextensible, slender (i.e. high-aspect ratio) body, whose shape is given by a curve in space and whose deformations are related to the rotations of a right-handed frame defined along the length of the curve. The filament has circular cross-section with radius and length . Its physical size is assumed to be sufficiently small such that its inertia, and that of the surrounding fluid, can be ignored. It can be subject to external forces and moments along its length, such as those due to the surrounding fluid, in addition to the internal moments due to bending and twist, and the internal force that arises due to constraints on filament motion. In addition, we consider the important case where, in the absence of external forces and torques, the filament may relax to a curve with nonzero curvature and twist.
2.1 Kinematics
The positions of the points along the filament’s centreline at time are given by , where is the filament arclength. Along with the position, there is a right-handed, orthonormal material frame at each and , see Fig. 1(a). At this stage, the unit vector and position, , are independent of each other, but later we will introduce the kinematic constraint, , that links these quantities and will lead to forces within the filament.
Filament deformation is determined by how the material frame varies with . Following Landau & Lifshitz (1986, §18); Powers (2010), this deformation can be described by the vector that describes the spatial-rate of rotation of the local frame with such that for any vector in the material frame, . Using and , can be expressed as
[TABLE]
where the first term is related to bending based on a curve with tangent , while the second term describes twisting about . By combining the derivatives of the frame vectors with the fact that , we can write the bending contribution as
[TABLE]
while the twist contribution is given by
[TABLE]
2.2 Constitutive law and internal moments
The internal moments, , are linearly related to filament bending, Eq. 2, and twisting, Eq. 3, through the bending, , and twist, , moduli. These moduli are related to the shear modulus, , and the Young’s modulus, , of the underlying filament material via , , where is the moment of inertia about any radial axis lying in the circular filament cross-section (Landau & Lifshitz, 1986, §16, 17). Specifically, we substitute Eqs. 2 and 3 into (19.7) of Landau & Lifshitz (1986, §19)* to yield*
[TABLE]
With this form of , the filament equilibrium configuration is straight and untwisted. Following **Lim et al. (2008); Lim (2010); Olson et al. (2013)**, to incorporate a nontrivial equilibrium shape into the model, we can introduce the preferred curvatures, and , and preferred twist, , such that
[TABLE]
In the absence of applied forces and torques acting on the filament, its shape is determined completely by , and . Additionally, by allowing , and to be time dependent, we can introduce shape changes that allow the filament to propel itself in a surrounding fluid without experiencing a net force or net torque, in accordance with the constraints on swimming at low Reynolds number (Purcell, 1977; Lauga & Powers, 2009).
2.3 Kinematic constraint and internal stress
Along with the internal moments, filament deformation will give rise to an internal force, , or the force on the filament cross-section. In **Lim et al. (2008); Lim (2010); Olson et al. (2013)**, this force is captured through a constitutive law that relates and the local frame such that in equilibrium,
[TABLE]
We instead impose Eq. 6 as a kinematic constraint, requiring that be the unit tangent to the curve defined by . This constraint gives rise to the internal force, and as such, it can be viewed as the Lagrange multiplier111Rather than considering force and moment balances along the filament, the filament model can alternatively be derived through the introduction of a Lagrangian (*Lim *et al., 2008, appendix). The internal force arises from the constraint term,
(7)
in the Lagrangian, which will also contain the bending and twist energies, as well as any external potential.* that enforces (6). The internal force provides both the tension in the direction that keeps the filament inextensible, but also the necessary normal force that couples the force and moment balances as we describe below.*
2.4 Force and moment balances
In general, the internal force and moment described above will balance external forces, , and torques, , per unit length that act along the filament, see Fig. 1(c). The external forces and torques can arise due to, for example, an external field such as gravity or short-ranged filament–filament interactions. In particular, we insist that they include the force and torque per unit length exerted on the filament by the surrounding fluid. We assume that the interaction with the fluid is governed by low Reynolds number hydrodynamics, and hence these forces and torques are linearly related to the filament’s motion. As described in Landau & Lifshitz (1986, §19); Powers (2010), the resulting force and moment balances are
[TABLE]
In addition to these equations, we will have conditions at the filament ends. For example, for filaments with free ends, we have222The boundary conditions require that . In our case, this will be satisfied since includes forces due to the surrounding fluid which will automatically balance any net force applied to the filament by other means (e.g. gravity, filament–filament repulsion).* and . Filaments that are tethered will need to satisfy conditions imposed on and at the tethered end. At this stage, it is important to note that Eqs. 8 and 9 establish the low Reynolds number mobility problem that couples the motion of all filaments through the fluid flows that they generate. This mobility problem, which we describe and solve in the discrete setting (see Section 2.7), yields equations for the translational and angular velocities of the discrete filament segments. The values of these velocities, however, will depend on which remains unknown. The values of are determined later by imposing the kinematic constraint, Eq. 6, when integrating the differential equations for the filament segment positions and orientations.*
2.5 Discrete force and moment balances
To solve the equations numerically, the filament is first discretised into segments of length as shown in Fig. 1(b). The segments have positions, , and orientations, , for . The segment orientations are the discrete representation of the frame vector, , at the segment position. Applying central differencing to Eqs. 8 and 9, the force and moment balances become
[TABLE]
where and are the internal moment and stress, respectively, between segments and ; is the external force per unit length on segment ; and is the external torque per unit length on ; see Fig. 1(d). Free end conditions are established by requiring that and . Similarly, for the kinematic constraint, Eq. 6, we have
[TABLE]
Multiplying Eqs. 10 and 11 by , we convert the force and moment per unit length balances to the force and torque balances such that
[TABLE]
for segment . The elastic torque on is given by
[TABLE]
Since the internal moment is directly linked to how the local frame changes with arclength, the computation of will depend on whether filament deformation is strictly planar, or fully three-dimensional. We present both cases below in Sections 2.8.1 and 2.9.2.
The constraint force and torque, respectively, on segment , are given by
[TABLE]
The values of are the Lagrange multipliers associated with the discrete version of the kinematic constraint, Eq. 12, that relates the translational and rotational motions of segments and . The Lagrange multipliers are unknown and must be computed as the filament moves and deforms.
Finally, the external forces and torques acting on a segment are defined as and , respectively. The external forces and torques can arise due to an external field, such as gravity or an external flow, but also through interactions with neighbouring filaments as mediated by the fluid, or through direct collisions and steric interactions. We describe here how we capture these interactions and, in particular, how hydrodynamic interactions lead to segment motion. In the case where the external forces and torques are due only to steric interactions and the surrounding fluid, we have
[TABLE]
where are the steric forces experienced by segment and and are the hydrodynamic force and torque on segment . We have chosen to introduce the negative sign here so that and denote, respectively, the force and torque segment exerts on the fluid.
2.6 Steric interactions
In our simulations, segments from different filaments and non-neighbouring segments from the same filament repel each other via a short-ranged, pairwise barrier force. For simplicity, we take this force to be a function of the centre-to-centre separation of the segments, though a more detailed representation of the segment shape could also be taken.
Specifically, in our simulations, the total barrier force on segment is
[TABLE]
*where the force on segment due to is *(Dance et al., 2004)
[TABLE]
and if . We also have that for any separation distance if and are neighbouring segments from the same filament. In Eq. 21, the displacement between segments is , while . The parameter controls the strength of repulsion at contact, while controls the range over which the barrier force acts. Recall that is the filament thickness. In practice, to prevent the segments from overlapping, we choose the reference force to be large compared to the typical bending force and set to ensure that the force only acts for particle separations that are 10% of their contact distance.
2.7 Hydrodynamics interactions and segment motion
The hydrodynamic forces and moments on the segments balance the internal stresses and moments and the inter-segment repulsive forces. In the absence of fluid inertia, this leads to a low Reynolds number mobility problem which provides a linear relationship between the forces and torques on the segments and their velocities and angular velocities,
[TABLE]
where is the vector of all segment velocity components. Similarly contains the angular velocity components for all segments, while and are the vectors that hold the forces and torques, respectively, that all the filament segments exert on the surrounding fluid. The methods presented in this paper can be coupled with a variety of hydrodynamic models and associated computational methodologies, including those that are mobility matrix-free, provided that they retain this linear relationship between the forces that the segments exert on the fluid and segment motion. Such approaches include Stokesian Dynamics (Brady & Bossis, 1988; Sierou & Brady, 2001), boundary integral methods (Pozrikidis, 1992), the method of regularised Stokeslets (Cortez, 2001) and the immersed boundary method applied to the Stokes equations (Peskin, 2002).
*In this work, we employ two approaches to perform the mobility matrix multiplication and resolve segment motion: a direct computation based on the Rotne–Prager–Yamakawa (RPY) tensors *(Wajnryb et al., 2013), and a matrix-free approach based on the force-coupling method (FCM) (Maxey & Patel, 2001; Lomholt & Maxey, 2003). Both approaches provide a similar resolution of segment hydrodynamic interactions based on the Stokes flows generated by singular and regularised force and torque distributions. In addition, RPY and FCM are closely related to the methods mentioned above and are chosen to demonstrate how our filament model can be used in conjunction with a wide range of methods, models, and solvers for filament motion in viscous fluids. We note that both RPY and FCM link the radius of the filament cross-section to the segment length. Thus, increasing the number of segments changes the aspect ratio of the filament, and hence its hydrodynamic response. While we will see in Section 4.3 that this hydrodynamic model provides accurate values for the hydrodynamic force on the filament over a range of aspect ratios, it could be important to be able to increase (or decrease) the spatial resolution of filament deformation without affecting hydrodynamic response. This decoupling is achieved, for example, with slender body theory (Tornberg & Shelley, 2004), as well as through a modified version of the method of regularised Stokeslets (Cortez, 2018). With this drawback in mind, we designed the overarching computational framework to be flexible and readily accommodate changes to the hydrodynamic model provided it satisfies Eq. 22.
2.7.1 RPY mobility matrix
*In one approach, we perform a direct pairwise evaluation of the RPY tensor **(Wajnryb et al., 2013) to compute the mobility matrix–vector multiplication. The operation count for this computation scales like the number of segments squared, though more sophisticated, faster techniques *(Liang et al., 2013) could also be used. With RPY, and given a collection of non-overlapping particles of radius and separations of , in a fluid with viscosity , the velocities of particle are given by
[TABLE]
For an unbounded fluid, the sub-tensors are defined as
[TABLE]
*where , is the identity matrix, and is the three-dimensional Levi-Civita symbol. We note that there are extensions and variants of these tensors which allow for particles of different sizes, overlapping particles, the inclusion of background shear flows, and the presence of no-slip boundaries *(Wajnryb et al., *2013; *Zuk et al., 2014; Swan & Brady, 2007).
2.7.2 Force-coupling method
We also present computations where the mobility matrix–vector multiplication is performed using the matrix-free FCM (Maxey & Patel, 2001; Lomholt & Maxey, *2003; *Liu et al., 2009). With FCM, the forces and torques the segments exert on the fluid are transferred to the fluid through a truncated and regularised force multipole expansion in the Stokes equations,
[TABLE]
where the sums are over all segments, is the fluid viscosity, is the fluid velocity and is the pressure. The Dirac delta functions in the multipole expansions are replaced by the Gaussians,
[TABLE]
Taking advantage of the ratios established for spherical particles (Maxey & Patel, 2001; Lomholt & Maxey, 2003), we take the Gaussian envelope sizes to be related to the filament thickness through and . In Eq. 33, we have ignored the contribution of the stresslets that would enforce a vanishing rate of strain within the segments (Lomholt & Maxey, 2003). The stresslets can also be included at the cost of several conjugate gradient iterations (Yeo & Maxey, 2010) each time the segment velocities and angular velocities are determined.
After the fluid velocity is found, it is spatially averaged using the same Gaussian functions to obtain the velocity and angular velocity of each segment . Specifically, we have that
[TABLE]
Following Yeo & Maxey (2010); Keaveny (2014)**, in our computations, the regularised forcing is first evaluated on a regular grid. The Stokes equations subject to periodic boundary conditions are then solved using a Fourier spectral method. Finally, the trapezoidal rule is used to integrate numerically Eqs. 37 and 38 to obtain the translational and angular velocity for each segment.
2.8 Planar filament deformation and motion
In order to keep track of filament deformation and motion, as well as compute the internal moments and forces, we need to describe how the local frame, , at each point along the filament rotates with time.
In many important problems, such as sperm locomotion through planar flagellar beats (Schoeller & Keaveny, *2018; *Yang et al., 2008), filament deformation and motion are restricted to a single plane and only bending occurs. Due to the absence of both elastic and preferred twist, the local frame at each point along the filament will rotate about a single, fixed direction. Without loss of generality, we can take this direction to be aligned with the fixed vector , as well as for all and . This direction is always normal to the plane of motion.
2.8.1 Computing the internal moments from the angles
As a result of taking the out-of-plane direction as , we have that in Eq. 5. The internal moment is then purely the result of bending and is given by
[TABLE]
with the understanding that since lies in the same plane for all , will be in the direction for all and . In the discretised model, the elastic moments are then given by
[TABLE]
where again it is understood that is in the direction for all and . The restriction to a single plane also allows for to be conveniently expressed in terms of an angle such that .
With this expression for found, we now have expressions for all the forces and torques on the segments appearing in the discrete force and torque balances, Eqs. 13 and 14.
2.8.2 Differential-algebraic system and time integration in two dimensions
With expressions for forces and , and torques and , we now move on to integrating the system in time. The mobility problem in Eq. 22 expresses the velocities, , and angular velocities, , for each segment , in terms of these aforementioned forces and torques through
[TABLE]
Due to the motion being restricted to the plane, we have and for each . We emphasise that we cannot simply evaluate Eq. 41 to obtain the segment motion since the values of appearing in (see Eq. 16) are unknown and are obtained by applying the kinematic constraint, Eq. 12, for each . Thus, the segment positions and angles can be updated by integrating numerically in time the differential-algebraic system of equations,
[TABLE]
Due to the numerical stiffness associated with the overdamped dynamics of elastic beams (Powers, *2010; *Nazockdast et al., 2017), as well as the need to ensure that the constraints are satisfied after advancing in time, we discretise in time using the implicit second-order backwards differential formula (BDF) (Ascher & Petzold, 1998) and impose the constraint at the updated positions and orientations (time ), such that
[TABLE]
where is the timestep size and and are given by the mobility problem, Eq. 41, with the right-hand side evaluated at using , , and . This yields a nonlinear system of equations for the positions, , and angles, , as well as the internal force Lagrange multipliers, , associated with each constraint. If is the number of filaments and is the number of segments per filament, the position and orientation updates, Eqs. 45a and 45b, and the constraints, Eq. 45c, constitute a system of nonlinear equations.
Before solving this system numerically, we first reduce the system size by substituting the constraint, Eq. 45c, into the position updates, Eq. 45a, as to reduce the number of degrees of freedom to those of a so-called robot arm whose motion is completely described by a single position and the orientations of each link comprising the arm. Specifically, in the case of a single filament, we replace for by
[TABLE]
to obtain the new, reduced system of equations,
[TABLE]
This system may be expressed as , where the solution contains the updated position for the first segment, all orientation angles and the Lagrange multipliers. For filaments each with segments, as a result of this substitution, the dimension of the nonlinear system is reduced from to .
2.9 Three-dimensional filament deformation and motion
In general, there is no restriction about how the filament can bend or twist, and the local frame along the filament can rotate about any axis. Rather than considering the local frame vectors explicitly, we keep track of these rotations using the unit quaternions that map the standard basis to the local frame at each point along the filament. Unit quaternions allow for successive rotations to be computed easily while requiring less storage and fewer floating point operations than rotation matrices. Quaternions also avoid gimbal lock that is typically experienced with Euler angles (Allen & Tildesley, 2017). Before describing how we employ quaternions in our methodology, we provide a brief overview of representing rotations using quaternions.
2.9.1 Representing rotations as quaternions
Quaternions (Allen & Tildesley, 2017; Dunn & Parberry, 2011; Vince, 2011) can be viewed as an extension of the complex numbers in that they have one real part and three imaginary components. As such, the quaternions inherit the notion of conjugation from the complex numbers, , and the norm of a quaternion is the Euclidean norm . Quaternions are also frequently identified with elements of as
[TABLE]
We use the operation to extract the vector in constructed from the last three entries of a quaternion.
Quaternions are subject to standard element-wise addition, real multiplication by a scalar, and the associative, non-commutative product,
[TABLE]
Under this product, the unit quaternions, i.e. quaternions with , form a group with identity and inverse .
As unit quaternions can be written as
[TABLE]
they can be identified with spatial rotations of an angle, , anticlockwise about the unit vector, . Specifically, the rotation applied to a vector is given by
[TABLE]
where is the image of after its rotation. This rotation can also be written as
[TABLE]
where is the rotation matrix whose entries are related to those of the quaternions through
[TABLE]
For successive rotations first by , say, followed by , we have that .
2.9.2 Computing the internal moments from the quaternions
As we use quaternions to describe how the local frame rotates as the filament bends and twists, it is convenient to have an expression for the internal moments in terms of the quaternions themselves.
Specifically, the unit quaternions provide the rotation of the standard basis to the local frame vectors, and we have that at each point along the filament and at each time,
[TABLE]
In Section A.1, we describe how to obtain the quaternion that satisfies this condition at .
Using Eq. 54 and the expression for the internal moments in terms of the frame vectors, Eq. 5, the internal moments can be expressed as
[TABLE]
where . The details of this derivation are presented in Section A.2. This expression is a special case of more general constitutive laws discussed in Zupan et al. (2009).
For the discretised system, the internal moments are then given by
[TABLE]
where the interpolated quaternion, , is constructed by performing half of the rotation from to . Specifically, we have where the square root is defined in Section A.1.
2.9.3 Differential-algebraic system and time integration in three dimensions
With the translational and angular velocities of the segments given by the mobility problem, Eq. 41, but with Eq. 56 used to compute , the positions and quaternions are updated by integrating the differential-algebraic system,
[TABLE]
The differential equation Eq. 58 for the quaternions describes the time evolution of the entire frame of segment such that taking for a vector is equivalent to (see Section A.3) integrating in time the perhaps more familiar expression
[TABLE]
*As when filament motion and deformation were restricted to a plane, we can update the positions using the second-order BDF scheme. While we would like to update the quaternions using a similar scheme, we must also ensure that after each update they continue to have unit norm as to continue to represent rotations. This could be accomplished by simply applying the second-order BDF scheme to Eq. 58 and subsequently normalising the result. We have found in practice, however, that this approach interfaces poorly with the quasi-Newton methods we use to solve the nonlinear system for the segment positions, orientations and Lagrange multipliers as it introduces timestep restrictions for numerical stability. Instead, we apply a geometric multi-step method *(Faltinsen et al., 2001) to perform the multiplicative update
[TABLE]
for some appropriate unit quaternion . This approach preserves the norm of the quaternion to machine precision and, in practice, leads to better numerical stability properties of the larger numerical method.
Applying the general framework set out in Faltinsen et al. (2001) to our specific case, we have that for sufficiently small , Eq. 58 has solution
[TABLE]
*where the vector is an element of the Lie algebra (in this case ) and the exponential map *(Iserles et al., 2000) is given by
[TABLE]
The Lie algebra element itself satisfies the differential equation
[TABLE]
*with . The function is the inverse of the differential of the exponential mapping and is given by *(Iserles et al., 2000)
[TABLE]
In a general sense, the Lie algebra element can be interpreted as the integral of the angular velocity over time. This holds true if the rotation is about a single axis. For example, if for a constant , then Eq. 64 with initial condition gives and we obtain .
*By introducing Eq. 64, the problem of updating the unit quaternions is transferred to one of updating the Lie algebra elements. This can be done using a standard, additive scheme with the only caveat being that since the previously updated Lie algebra elements are in the tangent spaces of the previous quaternions, we must also re-centre our coordinate system at the current quaternion when using a multi-step method *(Faltinsen et al., 2001). We show in Section A.4 that re-centring this yields the simpler second-order BDF scheme
[TABLE]
for the Lie algebra element of segment .
Combining the Lie algebra element update with that for the positions, as well as the constraints linking the segment orientations and positions, we obtain the system of equations
[TABLE]
for the updated segment positions, Lie algebra elements, and Lagrange multipliers. As in 2D, and are given by the mobility problem, Eq. 41, with the right-hand side evaluated at . For filaments, each discretised into segments, this is an system of nonlinear equations.
As done in the case of planar filament deformations, we reduce the system size before seeking a numerical solution by substituting
[TABLE]
into the position updates. This yields the new system of equations
[TABLE]
As done in the planar case, we may again express this system as , where the solution contains the updated position for the first segment, all the Lie algebra elements and the Lagrange multipliers. This substitution reduces the dimension of the nonlinear system from to .
3 Numerical solvers for the nonlinear systems
As described in the previous section, updating the filament positions, orientations, and the Lagrange multipliers requires finding a solution of the nonlinear system of equations of the form . This system is given by Eq. 47 for planar filament deformations and Eq. 69 in the fully three-dimensional case.
To solve the system, one would typically apply Newton’s method. Starting with as an initial guess of the solution, the solution is computed iteratively from
[TABLE]
where is the Jacobian of the system, (i.e. ). The process is continued until is smaller than a given tolerance.
For the systems given by Eqs. 47 and 69, the Jacobian will contain the mobility matrix, (see Eq. 22), as well as its derivative with respect to segment positions. In general, when one uses fast, matrix-free methods, the mobility matrix is not known explicitly. Even if the mobility matrix was known, the Jacobian would still be very complicated and challenging to compute, especially given the dependence of the mobility matrix on particle positions. Additionally, storing the Jacobian would incur significant memory costs when the number of filament segments is large.
To avoid the costs and complications associated with standard Newton’s method, we explore how solutions can be obtained using Jacobian-free Newton–Krylov and Broyden’s methods. The remainder of this section provides a description of these methods and, in particular, highlights the advantages and drawbacks associated with these approaches when applied to our specific computation. We also provide results from a series of tests comparing these methods and show that Broyden’s method with a suitable approximate Jacobian provides the most effective approach to finding the solution.
3.1 Jacobian-free Newton Krylov (JFNK)
The Jacobian-free Newton–Krylov (JFNK) method (Knoll & Keyes, 2004) uses Newton’s method to find the solution of the nonlinear system, but avoids an explicit computation of the Jacobian through judicious evaluations of the function, . With this method, the linear system of the form, for unknown , that arises during each Newton iteration is solved using a Krylov subspace method, typically GMRES (Saad & Schultz, 1986). Within GMRES, the Jacobian matrix–vector multiplication required at each iteration is approximated by the finite difference formula,
[TABLE]
Since in JFNK multiplication by the Jacobian is replaced by Eq. 71, each GMRES iteration requires evaluating .
In our case, this involves evaluating Eq. 47 or Eq. 69 and computing the translational and angular velocities of each segment, which itself requires performing a multiplication by the mobility matrix, the most computational costly aspect of the simulation. Thus, although JFNK does not require explicit knowledge of the Jacobian, it could still incur high computational costs if GMRES is slow to converge and many GMRES iterations are required for each Newton iteration. As a result, and as is often the case with GMRES, preconditioning may be required to limit iteration counts and to ensure good convergence rates.
3.2 Broyden’s method
Another approach, and the one that we will ultimately adopt, to obtain the solution of the nonlinear system is Broyden’s quasi-Newton method (Broyden, 1965). Like the traditional Newton’s method, Broyden’s method finds the solution iteratively, but instead uses an initial approximate Jacobian, , that is improved at each iteration through rank-one updates based on the secant equation. The choice of update, however, is not unique and the so-called “good” Broyden’s method updates the Jacobian itself (see Algorithm 1 in Algorithm 1), while the so-called “bad” Broyden’s method updates its inverse (Algorithm 2 in Algorithm 2). In both cases, we can employ a limited-memory version (van de Rotten, 2003) of Broyden’s method where the updates are stored as vectors rather than altering and writing to memory the full Jacobian at each iteration. Thus, the only matrix that may require storage is , which can be chosen with this memory requirement in mind.
For Broyden’s method to converge, not only must the convergence conditions for Newton’s method be satisfied, but also the initial guess for the Jacobian, , must be sufficiently close to . The process of choosing is similar to that of finding a suitable preconditioner in that one seeks with the condition that systems involving should be easily solved. Below we present two approaches that we have explored to obtain approximate Jacobians based on simplified, yet related, physical systems.
3.2.1 Jacobian-free Broyden’s method (JFB)
Our first approach involves applying the JFNK-style finite differencing, Eq. 71, to a function, , to perform the matrix–vector multiplications by and solve the systems appearing in Algorithm 1 or Algorithm 2 using GMRES. As this combines a Jacobian-free approach with Broyden’s method, we refer to this approach as Jacobian-free Broyden’s (JFB) method.
In the tests presented below, the function involves all the steps of the model presented in Section 2, however, we omit both the hydrodynamic interactions between filament segments by replacing the mobility matrix by and remove the steric interactions between filament segments. In this way, the computationally expensive mobility matrix multiplication is only performed once per Broyden iteration.
3.2.2 A block diagonal, explicit (EJB)
Along with JFB, we also construct explicitly the approximate Jacobian based on with a diagonal mobility matrix and no steric interactions between filaments. By differentiating , we obtain expressions for the entries of which we evaluate at each timestep. These expressions are provided in Appendix B.
In general, for simulations involving filaments, each with segments, will be for filament motion confined to a plane, and for the fully 3D case. When interactions between the filaments are ignored, will be block diagonal and consist of blocks of size for 2D motion and for fully 3D simulations. As a result, solutions to systems involving can be found by solving the smaller, independent linear systems corresponding to each block. Additionally, only the blocks, rather than the entire matrix, need to be stored in memory.
3.3 A comparison of the nonlinear system solvers
Here, we test and compare the different approaches for solving . Based on the tests presented below, we find that bad Broyden’s method using the explicit expressions for (see Appendix B) provides the most effective approach for solving the system arising from the discretisation of the differential-algebraic system, Eqs. 47 and 69.
*In our tests, we consider a monolayer of filaments, each formed of segments and subject to a constant force per unit length of magnitude in the plane of the monolayer. The simulations are performed using the MATLAB with FCM implementation described in the next section. The filaments are initially straight and are distributed uniformly and isotropically in the centre-plane of a periodic domain of size . Varying between and , we explore a range of filament lengths . FCM is used to resolve the hydrodynamic interactions between the filaments, and the Stokes equations are solved on a grid of points. With the domain size fixed, the computational cost of FCM scales linearly with number of filament segments. The filaments are allowed to settle for two characteristic sedimentation times defined as . The parameter is set such that the dimensionless elasto-gravitational number *(Cosentino Lagomarsino et al., 2005) is and filament deformation is significant. The timestep size is set to be .
In the following tests, the GMRES residual tolerance for both JFNK and JFB is and the JFNK finite difference parameter is set to . The Newton and Broyden iteration tolerances for the infinity norm of the residual relative to are both . These parameters have been selected to minimise wall times while still ensuring sufficient accuracy. To quantify solver performance, we measure average wall times, the average number of mobility matrix multiplications per timestep, and the average iteration count per timestep. These averages of are computed between simulation times and once the suspension has evolved away from its initial configuration. All simulations were performed on a 16-core 2.5GHz AMD Opteron 6380 processor.
3.3.1 JFNK versus Broyden’s method (JFNK v EJB-b)
We first compare wall times per timestep for the simulation run with JFNK with those using “bad” Broyden’s method with the explicit approximate Jacobian (EJB-b), Algorithm 2. A key aspect of using JFNK effectively is the choosing the correct preconditioner to accelerate GMRES convergence. Here, we explore two preconditioners. In the first instance, employ a diagonal right preconditioner that mitigates ill-conditioning due to the factor of appearing in entries of the Jacobian associated with the system Eqs. 47 and 69. Since other entries in the Jacobian are , small timesteps can increase the Jacobian condition number, resulting in slow GMRES convergence. As demonstrated explicitly in Appendix B for a drag-only mobility matrix, this factor of appears in the Jacobian as a result of variations with respect to the Lagrange multipliers.
Accordingly, the relevant diagonal entries of the right preconditioner are set to . Specifically, for a single, segment filament in 2D, the right preconditioner is
[TABLE]
with the system ordered as in Eq. 47.
To provide a balanced comparison with Broyden’s method, we also use as a left preconditioner the block-diagonal, explicit Jacobian described in Section 3.2.2 and derived in Appendix B. Based on the filament model with a diagonal mobility matrix and no steric iteractions, this preconditioner captures the local physics of the problem and its block structure allows for rapid system solves and reduced storage costs.
In Fig. 2, we compare the performance of JFNK with these preconditioners against that of EJB-b. The left panel of Fig. 2 shows the average wall time per timestep, while the right panel indicates the average number of mobility solves required per timestep. By “mobility solve”, we mean multiplying by the mobility matrix using the force-coupling method (FCM) as described in Section 2.7.2. We see that left preconditioning with the approximate Jacobian improves the performance of JFNK much more so than . Based on our wall time measurements, however, EJB-b is still approximately five times faster than JFNK. The additional costs per timestep associated with JFNK are due to the higher number of mobility solves that it requires. While EJB-b needs only one mobility solve per Broyden iteration, each Newton iteration in JFNK requires multiple GMRES iterations that each cost one mobility solve. Due to these additional costs, JFNK is not competitive for solving the nonlinear system.
3.3.2 Choice of Broyden’s method (EJB-b v JFB-g)
In this section, we compare the performance of a JFB implementation of good Broyden’s method (JFB-g) with EJB-b. Implementing good Broyden’s method with the explicit approximate Jacobian is computationally expensive since it would require solving a dense linear system that changes with each Broyden iteration (see Algorithm 1 in Algorithm 1). In contrast, for bad Broyden’s method, we avoid the JFB implementation since we update directly the inverse of the explicit approximate Jacobian. As a result, the system we need to solve at each Broyden iteration remains the same.
As discussed in Section 3.3.1 for JFNK, preconditioning is also required to accelerate GMRES convergence in JFB-g. In the tests presented below, we again trial the right preconditioner , as well as the approximate Jacobian, , as a left preconditioner.
The performance measurements for JFB-g with and without preconditioners and EJB-b simulations are presented in Fig. 3. For all approaches, the number of Broyden iterations per timestep increases similarly with both the number of segments per filament (Fig. 3(b)) and the number of filaments (Fig. 3(a)). This increase is expected, as increasing the number of filaments or segments increases the concentration of the suspension. This in turn increases filament interactions and the number of iterations required.
When we examine the wall time per Broyden iteration, however, differences between the approaches begin to emerge. Fig. 3(c) shows the wall time per iteration as a function of the number of filaments. We find that for all methods, the wall time per iteration scales linearly with the number of filaments . This is a result of the FCM mobility solve being the most expensive part of the iteration and its cost scaling linearly with the total number of segments.
The wall time per iteration as the number of segments is varied is indicated in Fig. 3(d). For JFB-g with as the preconditioner, or no preconditioner at all, the computational time increases at a more rapid rate of approximately . This rapid increase in cost is a result of the approximate Jacobian system size growing with and poor GMRES convergence. For example, we find that for a JFB-g simulation of filaments each with segments, without preconditioning, 94% of the wall time per iteration is spent in the GMRES routine to solve the Jacobian system. We do see, however, using an effective preconditioner, in this case , can bring down this cost. For EJB-b and JFB-g with preconditioner , we observe a slightly faster than linear growth in computational time with the number of segments. Despite the reduction due to preconditioning, EJB-b still outperforms JFB by approximately a factor of 2.
In addition to being the more computationally efficient method, EJB-b is demonstrably more robust. In practice, for JFB-g, we found that the convergence of Broyden’s method relied heavily on the quality of the solution given by GMRES. Consequently, a low residual tolerance was typically required. We also observed that GMRES could be prone to stagnation and found that GMRES convergence is sensitive to the choice of in Eq. 71.
3.3.3 Effects of tolerance and timestep size on performance
The Broyden iteration tolerance and the timestep size not only affect the accuracy of the simulation but also the computation time. In particular, while increasing the timestep size reduces the number of timesteps, larger timestep sizes require more Broyden iterations in order to reach the prescribed tolerance as initial guesses are extrapolated from solutions at past timesteps.
Measurements of the wall time for a full settling time, , at different timestep sizes and different Broyden’s tolerances are presented in Fig. 4 for EJB-b simulations of a concentrated suspension (area fraction 12.5%). In general, increasing the Broyden’s tolerance decreases the wall time, though with reduced effect at the smallest timestep sizes. Due to differences in iteration numbers, we find that decreasing the timestep from to increases the computation time by only a factor of . We also note that for , Broyden’s method fails to converge in the case where the tolerance was . Reducing the timestep further to , however, does not lead to additional reductions in the number of Broyden iterations, and consequently, the computation times for these simulations are approximately times greater than those with .
3.3.4 Effects of segment and filament number on performance
In general, the wall time impact of changing the number of segments, , or filaments, , can be complicated. In practice, we found that increasing local concentrations of filaments can lead to higher numbers of Boyden iterations. For each iteration, the predominant computations under EJB-b are the hydrodynamic mobility solve and the inversion of the approximate Jacobian, . Given that is block diagonal, the complexity of this operation is at worst . For filaments with , we have found that the hydrodynamic mobility solve is the more costly aspect of the computation.
The total number of segments impacts the computational work required to apply the mobility matrix. In Section 2.7, we discussed the two hydrodynamic models used in this work: the force-coupling method (FCM) for periodic domains and RPY for unbounded domains. As we use a direct pairwise summation, the computational work with our implementation of RPY scales like . For FCM, the computational work scales as , and is linked to the cost of constructing the body-force in the Stokes equations used to represent the forces the segments exert on the fluid. This scaling is seen in Fig. 3(c, d). It is important to note that for FCM this cost is in addition to the fixed cost of the FFT-based Stokes solver that scales like , where is the total number of grid points. Thus, the computational work associated with FCM will scale linearly with the number of segments only once the filament density is sufficiently high.
4 Accuracy and convergence
In this section, we present results from numerical experiments that demonstrate the finite difference scheme’s second-order convergence in space, as well as the second-order convergence in time for the 2D and 3D implementations of the BDF schemes. In addition, we provide a comparison between the RPY version of the hydrodynamic model and numerical solutions to the boundary integral equations for the hydrodynamic forces experienced by a straight filament.
4.1 Convergence in space
To demonstrate the second-order convergence of the finite difference scheme used to discretise the force and moment balances, Eqs. 10, 11 and 12, we consider an initially horizontal filament that is clamped at one end and subject to a vertical force of magnitude at the free end . The filament bends in the direction of the force, and after some time, reaches an equilibrium shape described by , the angle the filament makes with the vertical. After integrating the continuous force and moment balances, one finds that the equilibrium angle satisfies (Landau & Lifshitz, 1986, §19)
[TABLE]
Fig. 5* shows the error,*
[TABLE]
where is the angle for segment given by the finite difference scheme, as a function of . We observe second-order convergence to the analytical solution as we decrease .
4.2 Convergence in time
Here, we confirm the second-order convergence of the BDF2 schemes for two- and three-dimensional filament motion. To test the 2D scheme, we consider a pair of filaments that are initially straight, oriented vertically, and separated horizontally. The filaments are then allowed to settle as a result of a constant vertical force per unit length applied to each. For the 3D scheme, we consider the same settling problem, but for four initially straight filaments that are aligned vertically and arranged such that when viewed from above, they are at the corners of a square (cf. Fig. 19). In all 2D and 3D simulations, we use pairwise hydrodynamics via RPY and the dimensionless elasto-gravitational number is , where is the force per unit length acting on each filament. Each filament consists of segments and the simulations are run for . The settling time is denoted as . With these parameters, the filaments bend as they settle and move apart from one another until they each reach a horseshoe-like shape (cf. Fig. 15).
Fig. 6* shows the error based on the maximum Euclidean distance of the positions to the very accurate solution that was computed using timesteps per settling time and a tolerance of . We remind the reader that the positions along the length of the filament are computed using the values of the angles (in 2D) or quaternions (in 3D), and thus their accuracy depends on that of the BDF scheme. For both the 2D and 3D schemes, we do observe second-order accuracy provided that the tolerance for the Broyden’s iterations is sufficiently low. Thus, a key parameter in ensuring accuracy of the solution is not only the timestep, but the tolerance to which the resulting system of equations is solved.*
4.3 Accuracy of the hydrodynamic model
In this section, we examine the accuracy of the hydrodynamic models that we employ in this work to solve the mobility problem and obtain the translational and angular velocity of the segments. In this test, we consider the simplest case of a translating rigid, straight filament with a circular cross-section of radius . We compute the force on each segment using the RPY model described above and compare these values with those given by a second-order discretisation of the boundary integral (BI) equations for rigid body motion (Keaveny & Shelley, 2011). For the boundary integral solver, we take the surface of the filament to be that of a cylinder with radius that has spherical caps. The spacing between the segments with the RPY model is . We only compare the RPY model with the BI values since both solvers are configured for an unbounded fluid, while our implementation of FCM requires the imposition of periodic boundary conditions. We do expect, though, that for an unbounded fluid, FCM will yield force values similar to those of RPY as the expressions (Maxey & Patel, 2001; Lomholt & Maxey, 2003) for the FCM fluid velocity and particle mobility are asymptotic to their RPY counterparts.
Fig. 7*(a) shows the force on the filament segments as a function of given by the RPY model and BI computations for the case where the filament translates with speed in a direction perpendicular to its axis. The segment force values are shown for filaments with different aspect ratios, . For the BI solutions, the segment forces are computed by integrating numerically the surface tractions over sections of the cylinder of length , corresponding directly to those in the RPY model. The numerical integration is performed using the trapezoidal rule. Similarly, Fig. 7(b) shows the force on the segments given by the two methods when the filament translates with speed in the direction of its axis, while Fig. 7(c) shows the total force as a function of aspect ratio. We see that while the RPY model slightly overestimates the segment and total force values with a maximum relative error of approximately 4%, it does reproduce the general dependence of the segment force on given by the BI computations. We also observe that the error is largely independent of the filament aspect ratio.*
5 Simulations
In this section, we present results from several simulations that range from a single tethered filament whose position is fixed at one end, to many interacting filaments moving freely in the surrounding fluid.
5.1 Implementations
In the simulations that follow, we used one of three implementations of the method:
MATLAB with RPY*: For the single filament simulations, we have used a MATLAB implementation of the method in which the hydrodynamic interactions are accounted for through the RPY mobility matrices. We have made this implementation available on GitHub *(Schoeller et al., 2019).
For our simulations involving tethered filament arrays, we have accelerated the RPY evaluation by using MATLAB’s MEX functionality and performing the pairwise computation in C. 2. 2.
MATLAB with FCM*: For the two-dimensional simulations in periodic domains involving many sedimenting filaments or active swimmers, we use a MATLAB–MEX implementation of the method with FCM resolving the hydrodynamic interactions between the segments. Specifically, since with FCM we must solve the Stokes equations, we use C functions to loop through the grid points when assembling the FCM fluid forcing, and to loop through the wave numbers when inverting the Stokes operator. The remaining aspects of the computation, including the FFTs, are performed in MATLAB.* 3. 3.
C++ with RPY:* The final simulations of sedimenting filament clouds were performed with a more powerful C++ implementation of the method that takes advantage of the Armadillo linear algebra library (Sanderson & Curtin, 2016, 2018). In addition, the RPY computation, though still performed pairwise, was parallelised using OpenMP.*
5.2 Tethered filaments
5.2.1 Rotational dynamics of a single tethered filament
In Coq et al. (2008, 2009), an elastic filament is immersed in a tank of pure glycerine. The filament is connected to a motor at its base such that at rest, the filament forms an angle of [math] with the motor axis. The base of the filament is displaced from the motor axis by 2\text{,}\mathrm{m}\mathrm{m}* (see **Fig. 8(a)**). The filament has radius *$a=$435\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}, and its length is varied between and . A dimensionless Sperm number, Sp, that characterises the ratio of viscous to elastic forces is introduced as
[TABLE]
where is the transverse drag coefficient for a straight rod, is the rotation rate of the motor, and is the bending modulus of the elastic filament. By varying , a range of Sp can be explored and the steady-state distance, , of the free end from the rotation axis is recorded. At lower Sp the filament barely deforms from its straight configuration and undergoes a nearly rigid body rotation such that d\approx\delta_{0}+L\sin($$). At larger Sp, the filament slowly collapses onto the rotation axis and .
In order to simulate this experiment, the system of equations Eq. 67 is supplemented with the additional constraints such that for timestep ,
[TABLE]
to tether the filament base to the motor, and
[TABLE]
to rotate the filament base with the motor. These constraints are associated with two new vector Lagrange multipliers and , and consequently the first segment experiences additional constraint forces and torques given by
[TABLE]
where is the matrix such that . The total differential algebraic system is modified by substituting these new constraints, Eqs. 76 and 77, into the update equations for the first segment.
To compare our simulations to the experimental results, we match the ratio and the length range by varying the number of segments, . For each length considered, the simulation values of or are varied to explore the range of . In addition, the simulations are run using the RPY mobility matrices (see Section 2.7.1) to capture segment hydrodynamic interactions.
In Fig. 8(b), we present the filament end distance, , as a function of from our simulations for the range of filament lengths. Along with this data, we show the values measured in the experiments in Coq et al. (2008). Though the exact filament lengths associated with the data points were not indicated in Coq et al. (2008), the values of measured in the experiments lie broadly within the range given by our simulations. For comparison, we have also included similar results given by the gears model (GM) from Delmotte et al. (2015). While these results largely coincide with those from our simulations and the experiments, there are some differences that may be attributed to the interpretation of the offset distance, , as discussed in Delmotte et al. (2015).
5.2.2 Coordination in cilia arrays
In this section, we extend our study of a single tethered elastic filament to the simulation of a model ciliary array. Cilia are biologically occurring elastic filaments that undergo time-periodic motion and tend to exhibit metachronal waves and coordination (Brennen & Winet, *1977; *Guo et al., 2018). A simplified model of cilia dynamics involves the “geometric switch” (Guo et al., 2018; Kim & Netz, 2006), wherein a driving torque is applied to the base of an elastic filament until some fixed angle is reached with respect to the vertical. At this point, the direction of the torque is reversed until the filament reaches an angle of , when it switches again and so on. In Kim & Netz (2006), both the direction and the magnitude of the driving torque change once the critical angle is reached (i.e. with ), allowing for distinct “effective” and “recovery” strokes. The resulting filament motion, which is non-reciprocal due to its elasticity, allows for the net pumping of the fluid.
Using our methodology, we first simulate a single model cilium comprised of segments. The position of the base segment is constrained to remain fixed to where it is attached to a no-slip wall. The hydrodynamic effects of the wall are included through the RPY tensor for a half-space given in Swan & Brady (2007)**. In addition, for numerical stability of our implicit scheme, we allow for a smooth transition between the torque values once or are reached. This is accomplished by introducing a transition period of several timesteps during which the driving torque is taken from a sigmoid curve connecting to . Fig. 9(a) shows the shape of the cilium over a period of oscillation given by our simulations. The differences in filament shape and speed due to the different values of and can be clearly observed.
Since the switching between and is determined by the tangent vector at the model cilium’s base, Kim & Netz (2006)* observed that hydrodynamically interacting pairs of model cilia can phase-lock regardless of their initial phase difference. Using our model, we performed similar simulations of two interacting model cilia that are separated by a distance at the base. Each model cilium is again comprised of segments and has torque ratio and . The torque applied at the base of the model cilia is perpendicular to both the base-to-base separation vector and the surface normal. Figure 9(b) shows , the angle between the surface normal and base tangent for the second filament, when the first filament switches from the fast to the slow stroke. Similar to the kind of coordination observed in Kim & Netz (2006)**, we find that regardless of the initial phase difference, reaches the same value of for each simulation.*
Expanding from the pair simulations, we now consider arrays of these cilia and investigate the effect of inter-filament spacing on coordination. We arranged filaments on an grid at a base-to-base separation of in the plane of motion and in the orthogonal plane (see Fig. 10). Fixing , we vary the ratio and explore its effect on cilia self-organisation. Building from the coordination parameter of Guo et al. (2018), we introduce the time-dependent parameter
[TABLE]
to measure the coordination of cilia, where is the sign of the driving torque on cilium at time . We set the moving window size, , to be , where is the period of an isolated cilium. When , there is perfect synchrony. Values of close to zero indicate a lack of coordination amongst the cilia.
In Fig. 11, we show over time for arrays of cilia with different values of ranging from to . Each array was initialised with all cilia in-phase and the simulations are run to a final time of approximately . Overall, we find that denser arrays provide a greater level of coordination, as indicated by the higher values of . The decrease in coordination with is also evident in Fig. 12, where we show snapshots of the array near the final time for the different cases of . At the extremes of the range, quickly reaches a value about which it exhibits fairly small fluctuations. For intermediate , however, we see more complex behaviour where fluctuates over a wide range with intermittent periods of coordinated and uncoordinated motion.
5.3 Bidisperse suspension of undulatory swimmers
In this section, we employ the 2D version of our computational model to simulate a suspension of slender undulatory swimmers moving in the -plane. The suspension consists of two populations (indexed by and ) that have different swimming gaits. Similar simulations (Agrawal & Babu, 2018) have been performed using multiparticle collision dynamics for two swimmer populations with different frequencies, but the same waveform.
All swimmers have length and are formed of segments. The motion of the swimmers is driven by a time-dependent preferred curvature,
[TABLE]
in Eq. 5, where is the amplitude, is the wavenumber for population , is the undulation frequency with as the undulation period, and is the phase. The two different swimming gaits are prescribed by specifying different and .
In our simulations, the curvature amplitude is , while the undulation frequency is set such that the ratio of viscous to elastic forces is . The phase, , for each swimmer is drawn randomly from a uniform distribution. To achieve sufficiently different waveforms, we set and . Figure 13 shows waveforms for isolated and swimmers with these parameters. With the swimming speeds and , the swimmers from population travel at less than a quarter of the speed of those from .
In our simulations, we consider a suspension of 100 swimmers with 50 in population and 50 in . The swimmers are all initially straight and distributed uniformly and isotropically in the centre-plane of a periodic domain of size . FCM is used to resolve the hydrodynamic interactions with the fluid flow being solved on a grid. Steric interactions are implemented as described in Section 2.6.
*The configuration of the suspension at three times is shown in Fig. 14. After 150 undulation periods, we see that all swimmers have clustered, with swimmers tending to cluster with those from their own population, i.e. with a similar waveform. Clustering is also observed in planar monodisperse suspensions of model sperm cells **(Schoeller & Keaveny, 2018; Yang et al., 2008) and slender undulatory swimmers *(Yang et al., 2010), as well as in planar bidisperse suspensions (Agrawal & Babu, 2018) of undulatory swimmers. We first observe swimmers from forming clusters with each other, while the slower swimmers from remain dispersed in the surrounding fluid. Once clusters are established, we begin to see swimmers joining these clusters, but tending to remain close to other swimmers. Interestingly, filament softness has contributed to the development of these slightly mixed clusters. We sometimes observed high-wavenumber swimmers deforming to match the shape of neighbouring low-wavenumber swimmers.
5.4 Sedimentation
In this section, we use our methodology to study filament sedimentation. We both revisit the case of a single filament falling under gravity, as well as explore cases of multiple filaments, ranging from small ensembles to clouds and suspensions. We address these situations in cases where filament motion is restricted to a plane, as well as those where their motion is completely three-dimensional.
5.4.1 Settling filament
We first revisit the case of a single filament settling under gravity and measure its deformation and effective drag coefficient. This problem, first investigated numerically by Cosentino Lagomarsino et al. (2005) and more recently in experiments by Marchetti et al. (2018), has also been used previously as a test problem for filament models (Keaveny, *2008; *Delmotte et al., 2015).
In our simulations, we consider an initially straight and horizontal filament formed of segments in an unbounded fluid with the hydrodynamic interactions resolved using the RPY tensor with hydrodynamic radius, . As in Section 3.3, a constant force per unit length of magnitude is applied vertically to the filament. Using the parameters arising in the problem, we can define the settling time as , as well as the elasto-gravitational number . Running the simulation with timestep size , the filament is allowed to sediment until it reaches a steady-state shape and steady-state settling speed, . The filament at various points in time is shown in Fig. 15 for the case where . Once the steady state is reached, we measure the normal deflection , defined as the distance between the highest and lowest points of the filament (see Fig. 15), and the effective drag coefficient , where and is the value of in the limit corresponding to a straight filament.
The values of and as a function of are presented in Fig. 16. We compare our results with those presented in Delmotte et al. (2015) for both their gears model, where filaments are formed of 34 particles with no separation between the particles, and a joint model, formed of 31 particles with gaps of between the particles. As we have , our spacing is equivalent to that used in the joint model. In all cases, the RPY mobility matrices are used to capture hydrodynamics interactions.
All models show excellent agreement in their respective measurements of normal deflection, . In all cases, the models predict the same initial linear increase of with , followed by a plateau to the same value at high . For the drag coefficient, , the differences, though still small, are more pronounced. In particular, our measurements find good agreement with the gears model for as the sudden decay in given by both models coincide. For higher , however, we observe better agreement with the joint model, which yield a similar decay with increasing . We suspect that the agreement with the gears model for lower is due to the better resolution of the bending moments given by these models. At higher , we attribute the agreement between our results and those given by the joint model to the similar viscous force on the filament that results from setting .
5.4.2 Sedimentation of large 2D suspensions
Here, we investigate the sedimentation of a large, planar suspension of filaments in a periodic domain. This set of simulations was inspired by Gustavsson & Tornberg (2009); Tornberg & Gustavsson (2006)**, where a similar simulation was also performed in a periodic domain, but for a suspension of rigid filaments that can move in all three dimensions. We note as they do, that sedimentation in periodic domains, while a canonical problem, is quite different from that in unbounded or bounded systems and can exhibit dependencies on domain sizes. Nevertheless, these simulations do provide a context with which to explore the effects of filament flexibility on sedimentation.
The simulations presented here are larger versions of the test problem described in Section 3.3. We also run these simulations for longer times. Specifically, we consider monolayers of 100, 500, and 1000 filaments, each formed of segments. The dimensionless elasto-gravitational number is set at two values, and , to examine the effect of deformation on suspension dynamics as the filaments settle for as long as . The timestep is chosen such that . The filaments are initially straight and are distributed uniformly and isotropically in the centre-plane of a periodic domain of size . FCM is used to resolve the hydrodynamic interactions between the filaments, with the Stokes equations being solved on a grid of points. Steric interactions are implemented as in Section 2.6.
Snapshots of the sedimenting suspensions for with and are shown in Fig. 17. In general, we find that both stiff () and flexible () filaments over time form a single large cluster surrounded by isolated filaments that are found to tumble, but are most often oriented in the direction of gravity. We do see, however, that for the flexible filaments, the cluster formed is of much higher density, and the filaments surrounding it can buckle and flex, as well as tumble, due to the shear induced by the falling cluster.
Following Gustavsson & Tornberg (2009)**, for each of the simulations that we have performed, we compute the mean sedimentation velocity,
[TABLE]
where is the vertical centre-of-mass velocity of filament . The resulting values over time are shown in Fig. 18. We observe that for each concentration, the settling velocity and its growth rate increase with filament flexibility (higher ). We suspect that this is a result of two factors. First, as we observed in Section 5.4.1, flexibility leads to a reduction in the filaments’ effective drag coefficient. Second, we also see in Fig. 17 that deformable filaments form denser clusters. This increase in local filament density in turn enhances interactions, yielding higher sedimentation velocities.
In previous work (Mackaplow & Shaqfeh, 1998; Gustavsson & Tornberg, 2009) on rigid fibre sedimentation in periodic domains, the onset of a plateau in the mean settling velocity coincided with the formation of a single cluster in the centre of the domain that is surrounded by an upwelling of clear fluid. Additionally, observed fluctuations about the plateau value corresponded to cluster break-up and reformation events that occur periodically. In our case, we find that even after , the mean settling velocity continues grow with time. For each case of , the setting velocity growth rate is higher for the more flexible filaments, but perhaps more surprising is that the growth rate appears to be independent of . Indeed, increasing the filament density by setting , we see that over the that for both , the mean sedimentation velocity parallels that for . The continued growth may be attributed to the observation that there many filaments have yet to join the cluster. Additionally, we suspect that the continued growth in sedimentation velocity as compared to the fully three-dimensional cases (Mackaplow & Shaqfeh, 1998; Gustavsson & Tornberg, 2009) is due to both the larger lateral distances in our simulation that the filaments must traverse, as well as their restriction to move only in two dimensions. We suspect also that the reduced out-of-plane thickness enhances the hydrodynamic interactions between the filaments, allowing for the higher growth rates in our simulations.
5.4.3 Tumbling, settling and buckling of a filament square
Gustavsson & Tornberg (2009*)** have shown that rigid fibres evenly distributed along a circle can periodically tumble as they settle. Similar periodic orbits have been seen both numerically (Claeys & Brady, 1993) and experimentally *(Jung et al., 2006) for ellipsoids and other rigid structures released from an initially symmetric configuration. Here, we explore how the introduction of filament flexibility can affect these dynamics.
In our simulations, we consider four filaments whose centres coincide with the corners of a horizontal square of side length . The filaments are initially aligned vertically in the direction of gravity. We discretise the filaments into segments and, to allow for an unbounded fluid, we resolve their hydrodynamic interactions using the RPY mobility matrix. Steric interactions are implemented as in Section 2.6. The timestep is chosen such that . The elasto-gravitational number, , is varied from (very stiff filaments) to (very soft filaments). The long-term behaviour reveals significant differences between these cases, including symmetry breaking (Fig. 19). We note that similar observations have been made previously **Llopis et al. (2007); Saggiorato et al. (2015)**.
At the lowest value of , we observe periodic tumbling (Fig. 19a) where the filaments first separate and then reapproach when they return to their initial configuration. This matches the periodic orbits seen in Gustavsson & Tornberg (2009); Jung et al. (2006)**. For and , while we again observe the filaments returning very close to their initial configuration, we see that the time it takes to return to this configuration decreases with , as does the distance between the filaments. For the case , when returning after tumbling, the filaments, whose deformation is now apparent, approach each other so closely that they collide and experience steric interactions (Fig. 19b). At this point, the filaments cease tumbling and instead begin to settle as a group, rearranging their relative positions as they move downwards. At , we see that filaments do not exhibit any tumbling and instead adopt a horseshoe-like shape and a steady configuration as they settle, see Fig. 19c. Finally, for , while we see that at early times the dynamics are similar to those for , we observe also that this state is unstable. The accumulation of round-off errors are sufficient to trigger the instability (Fig. 19d) and the onset of more complex filament behaviour.
5.4.4 Sedimentation of flexible filament clouds
*As a final demonstration of our methodology, we consider the dynamics of a cloud of flexible filaments settling under gravity. For clouds of rigid, spherical particles (Nitsche & Batchelor, *1997; *Metzger et al., 2007), as the cloud settles, the particles within the cloud move along trajectories reminiscent of the streamlines within a settling spherical drop (Batchelor, 1967). As the cloud settles, it sheds particles at the cloud axis, until it becomes a torus, leading to the cloud dividing. While similar studies have been performed for clouds of rigid fibres *(Park et al., *2010; *Nazockdast et al., 2017), showing how particle anisotropy accelerates cloud break-up, the role of flexibility, to the best of our knowledge, has not been explored.
In our fully three-dimensional cloud simulations, the filaments are distributed uniformly and isotropically in cubes (Figs. 20 and 21) of side length and using the initialisation method describe in Section A.5. The filaments are discretised into segments and the timestep size is . The hydrodynamic interactions in the unbounded domain are resolved using the RPY motility matrix and steric interactions are implemented as in Section 2.6. To highlight the effect of filament flexibility, simulations were performed with elasto-gravitational numbers and .
Fig. 20* shows snapshots from the cloud simulations with filaments for both and . The initial conditions for these simulations are identical. At short times, we observe that both clouds exhibit very similar behaviour. Later in the simulations, however, marked differences arise. The cloud of stiff filaments () splits into two clusters as a result of filament depletion at its centre. This behaviour matches results from Park et al. (2010). The cloud consisting of more flexible filaments () breaks up later, and in a much more complicated way. We see that it disintegrates into several clusters, rather than splitting in two. We suspect that this could be connected to individual filaments’ orbits within clouds, as well as their sedimentation velocities. We have also performed a larger, though less dense, cloud simulation with filaments and . As shown in Fig. 21, we find that the features of the cloud’s initial shape persist as the cloud overturns and only starts to disappear once the cloud begins to break up.*
6 Discussion and conclusions
*In this paper, we have presented a scalable, efficient and flexible method for performing large-scale simulations of passive and active filaments and their suspensions. In particular, we have shown how to describe the fully 3D filament deformations due to bending and twisting through the use of quaternions, and in doing so, demonstrated how quasi-Newton methods coupled with an implicit, geometric multi-step scheme can be used to overcome the numerical stiffness associated with deformation while ensuring constraints on filament motion are satisfied. Our framework can be easily integrated with many existing models and methods for resolving hydrodynamic interactions. In addition, the second-order spatial discretisation could be substituted for a higher-order finite difference method on the equispaced grid, though new expressions for the approximate Jacobian would need to be established given the change in stencil. Our usage of iterative schemes for the nonlinear systems that arise ensures that efficient, matrix-free methods for the hydrodynamics are also applicable, which we have also demonstrated. We have shown that our method can be used to study many relevant problems in biofluid dynamics, as well as suspension mechanics, potentially at scales previously inaccessible. We encourage the reader to experiment with the MATLAB/Octave implementation *(Schoeller et al., 2019) of the method to study similar problems or expand upon it to explore new directions.
*While we have developed and explored the method exclusively in the limit where fluid inertia is negligible, it would be of interest to understand how our framework could be adapted for use with fluid solvers where inertia is included (Wiens & Stockie, 2015), to study, for example, filament suspensions in turbulent flows **(Rosti et al., 2018). Additionally, it could be of interest to explore filament motion in non-Newtonian fluids, such as viscoelastic fluids, especially with respect to modelling flagellar motion and cell swimming *(Teran et al., 2010). Our method could be used in the context where the filaments are components in more complex physical systems, such as networks, or, by incorporating Brownian motion (Keaveny, *2014; *Delong et al., 2015), suspended polymers in solution, or entangled polymers in melts.
Acknowledgements
EEK and AKT gratefully acknowledge support from EPSRC Grant EP/P013651/1. SFS gratefully acknowledges funding by an Imperial College President’s PhD scholarship. TAW is thankful for funding through an EPSRC Studentship (Ref: 1832024). All authors also thank Smitha Maretvadakethope for a detailed reading of the manuscript.
Appendix A Quaternions
A.1 Quaternion initialisation
In this section we detail how to initialise quaternion encoding the rotation from the standard basis to the material frame . Since , we need only provide initial conditions for the tangent and one of normal vectors. The main idea is to determine the quaternion which rotates a given unit vector to another and then apply this twice in succession by first mapping to and then mapping the image of under this first rotation to .
Let and be unit vectors. Recall that , , where is normal to the plane spanned by and , and is the angle between the two vectors as measured anticlockwise about . Recalling Eq. 50, we see that the quaternion represents the rotation of about the vector . Since we need only rotate by to align with , we seek the quaternion for which . By expanding and rearranging we find that a square root of a quaternion may be defined as
[TABLE]
It is readily checked that for all : (a) , and (b) . With this defined, we produce the quaternion initial conditions as follows:
Let and be the initial tangent and normal vectors to the filament centreline. 2. 2.
Evaluate the quaternion which maps as . 3. 3.
Let be the image of the -axis direction under this first rotation. 4. 4.
Evaluate the rotation mapping as 5. 5.
Assign the quaternion an initial value of .
A.2 Internal moments in terms of quaternions
In this section, we derive an expression for the internal moments, Eq. 5, expressed solely in terms of the unit orientation quaternion. Recall Eq. 5:
[TABLE]
Since , we have that
[TABLE]
Next, considering , a function of arclength and time , we observe that
[TABLE]
where we have used , a consequence of . Now, since , we have that is a pure quaternion (i.e. one with zero scalar part) and thus may be unambiguously identified with its vector part . Observing that for pure quaternions we have , the above yields
[TABLE]
Analogous expressions hold for and . Using the invariance under circular shift of the scalar triple product, we find that
[TABLE]
Thus, we simplify,
[TABLE]
which yields
[TABLE]
Defining we find that Eq. 85 reduces to
[TABLE]
A.3 The quaternion ODE
In this section, we show that if a frame vector satisfies , where is defined through a unit quaternion and an initial vector (i.e. ), then the quaternion must satisfy . Employing the same simplifications as in Eq. 86, we observe that
[TABLE]
This implies
[TABLE]
Given that this holds true for all the vectors of the frame, we obtain , and since is a pure quaternion we have
[TABLE]
Dividing by 2 and right-multiplying by yields the desired result.
A.4 BDF2 in the Lie algebra
In this section, we derive the form of the BDF2 update equation which applies to the Lie algebra elements. This derivation follows the more general framework laid out in Faltinsen et al. (2001). Let denote the Lie algebra element corresponding to the quaternion at timestep expressed in the coordinate system from timestep , that is,
[TABLE]
To integrate Eq. 64 using BDF2 we need to solve
[TABLE]
given and . The coordinate transform between tangent spaces needs to be such that a given Lie algebra element represents the same quaternion independent of the current coordinate system, hence
[TABLE]
needs to hold.
To this end, we use the Baker–Campbell–Hausdorff formula which is defined to satisfy
[TABLE]
While in our case (* with BDF2), there is a closed-form expression for the Baker–Campbell–Hausdorff formula (Engø, 2001), it will only be necessary to use the identity .*
For Eq. 99 to hold, the Lie algebra elements need to be transferred from one tangent space to the next. It follows from the condition that this can be achieved via
[TABLE]
Substituting this into Eq. 98 we obtain
[TABLE]
Observing in the right-hand side of this expression that
[TABLE]
and, since from Eq. 97, that
[TABLE]
we obtain
[TABLE]
Noticing that the Lie algebra elements only appear here expressed according to the coordinate system in which they were first produced, we may unambiguously write the BDF2 update equation as
[TABLE]
A.5 Filament initialisation in cloud simulations
The isotropic initialisation of straight filaments works as follows. We generate uniform random points in , where is the linear box size. These points set the positions of the filament segments at the ends. We then generate four-dimensional random vectors drawn from the standard normal distribution. After normalising, these vectors provide the initial quaternions for all segments comprising each filament. We then initialise the remaining segment positions for each filament in straight lines along directions based on the initial quaternions. If there is any overlap between filaments, one is discarded and we reinitialise its position and quaternion. This process generates a uniform distribution of filament orientations (Karney, 2007).
Appendix B Form of the approximate Jacobian
In this section, we detail how we generate an analytic approximation to the Jacobian, , by including only local physics for each particle – a local Stokes drag model, without hydrodynamic or steric interactions. The resulting approximate Jacobian retains the essence of the full Jacobian, but can be rapidly generated and posesses a block diagonal structure.
*The supplementary code *(Schoeller et al., 2019) implements the generation of the Jacobian. Here, we demonstrate how the expressions are derived.
B.1 In 2D
The system as described in Eq. 47 is formed of the position update for particle in each filament, the orientation updates for all , , and the modified constraints for :
[TABLE]
This system is repeated for each filament.
In general, the velocities and angular velocities of the filament segments are found by solving the mobility problem,
[TABLE]
where is the vector of all segments’ velocity components, similarly of all angular velocity components, of all force components and of all torque components acting on the filaments’ segments.
In the approximate Jacobian, we implement a local Stokes drag model, where the mobility matrix is diagonal, giving and . The only force we apply is the constraint force, Eq. 16, and the torques applied are the elastic torque, Eq. 15, and the constraint torque, Eq. 17.
Recall that for planar motion, without loss of generality, we take the third frame vector to be aligned with the fixed vector . This direction is always normal to the plane of motion.
Together, we have
[TABLE]
repeated for each filament. Since the approximate Jacobian assumes no interaction between filaments, we can build the approximate Jacobian for a multifilament system from independent blocks, placed along the diagonal, with each block holding a single-filament Jacobian. In the following we discuss the derivation of a single filament’s approximate Jacobian without loss of generality.
Having written out the system in Eq. 109, it can be differentiated with respect to the state variable, , to form the Jacobian
[TABLE]
The remaining expressions for the submatrices relating segments and are given by, after some algebraic manipulation and dropping the time (superscript) indices for clarity,
[TABLE]
where and is the Kronecker delta.
B.2 In 3D using Lie algebra elements
Here we discuss the differences in the construction of the approximate Jacobian in 3D compared to 2D. Again we consider the system , this time represented by Eq. 69. The standard BDF2 update for the first segment in the filament is handled analogously to the 2D case, Eq. 109a, so we restrict our attention to the Lie algebra update equations and the modified constraints. Substituting in the Stokes drag approximation, each filament-level sub-block of the Jacobian is based on the equations
[TABLE]
where the Lie algebra update equations are for all segments , and the modified constraints are for segments .
For the sake of legibility, we discard any terms with a time index other than , as they will have zero derivative in all cases, and drop the time superscript. Hence, to produce the Jacobian it suffices to consider the equations
[TABLE]
The Jacobian then corresponds to the system in Eqs. 109a and 113, differentiated with respect to the state variable, , forming
[TABLE]
The expressions for and are the same as in the 2D case, Eqs. 111a and 111e. We therefore consider the remaining three blocks.
We first observe that
[TABLE]
where is the Kronecker delta. By taking Eq. 64 and replacing by the first term of its Taylor series expansion,
[TABLE]
we can approximate . Combining this with from Eq. 60, we have
[TABLE]
where we have introduced the notation for the skew-symmetric matrix satisfying (and satisfying ). Hence,
[TABLE]
For the final blocks, we employ all three terms written in the Taylor series expansion in Eq. 116 to reduce Eq. 112a to
[TABLE]
From this we can see that the matrices of interest take the forms
[TABLE]
The problem of constructing these matrices thus reduces to that of evaluating the derivatives of the torque on segment . Segment experiences elastic and constraint torques
[TABLE]
where is the elastic moment between segments and . Since the elastic moments do not depend on the Lagrange multipliers, we have
[TABLE]
Now considering the derivatives with respect to the Lie algebra elements, we have
[TABLE]
If we use a frame vector approximation for the elastic moment,
[TABLE]
where , then we can write
[TABLE]
where the last term is a dyadic product. Similarly,
[TABLE]
The derivative with respect to any other Lie algebra element is zero, and thus we have provided all of the elements necessary to form an approximate Jacobian. We find that level of approximation employed here works well in our implementation.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Agrawal & Babu (2018) Agrawal, A. & Babu, S. B. 2018 Self-organization in a bimotility mixture of model microswimmers. Physical Review E 97 (2), 020401.
- 2Allen & Tildesley (2017) Allen, M. P. & Tildesley, D. J. 2017 Computer Simulation of Liquids , 2nd edn. Oxford University Press, Oxford.
- 3Ascher & Petzold (1998) Ascher, U. M. & Petzold, L. R. 1998 Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations , vol. 61. SIAM.
- 4Baker et al. (2009) Baker, E. L., Bonnecaze, R. T. & Zaman, M. H. 2009 Extracellular matrix stiffness and architecture govern intracellular rheology in cancer. Biophysical Journal 97 (4), 1013–1021.
- 5Batchelor (1967) Batchelor, G. K. 1967 An Introduction to Fluid Dynamics . Cambridge University Press.
- 6Brady & Bossis (1988) Brady, J. F. & Bossis, G. 1988 Stokesian Dynamics. Annual Review of Fluid Mechanics 20 (1), 111–157.
- 7Brennen & Winet (1977) Brennen, C. E. & Winet, H. 1977 Fluid mechanics of propulsion by cilia and flagella. Annual Review of Fluid Mechanics 9 , 339–398.
- 8Broyden (1965) Broyden, C. G. 1965 A class of methods for solving nonlinear simultaneous equations. Mathematics of Computation 19 (92), 577–593.
