1 Introduction
Throughout human history, flexible filamentary structures have been essential to human society, serving various purposes such as fastening, sailing, climbing, weaving, and hunting. As our understanding of material properties improved, so did our ability to engineer rods with enhanced material properties (e.g., flexibility, strength, and resilience). This in turn has resulted in the need to study and better understand the complicated mechanics of filaments. Thus, several previous works have sought out to understand the various mechanics of rodlike structures including the deployment of rods jawed2014coiling ; jawed2014pattern ; jawed2015geometric , elastic gridshells baek2018form ; panetta2019x ; baek2019rigidity , plant growth goriely2006mechanics , knots jawed2015untangling ; audoly2007elastic ; patil2020topological ; patil2020discharging ; grandgeorge2021mechanics , and propulsion of bacterial flagella jawed2015propulsion ; jawed2016deformation ; jawed2017dynamics .
With realworld experiments being costly and tedious to implement, the need for accurate physicsbased numerical simulations is prevalent. Such simulations not only allow for advanced mechanicsbased study, they also open up the avenue to challenging problems in robotics ranging from simulating soft robot dynamics to closing the sim2real gap for deformable material manipulation. In recent times, discrete differential geometrybased (DDGbased) simulations have shown surprisingly successful performance in capturing the nonlinear mechanical behaviors of rod structures bergou2008discrete ; bergou2010discrete . However, frictional contact handling still lacks a descriptive understanding.
Contact handling methods can generally be divided into three distinct categories: impulse methods, constraintbased optimization methods, and penalty energy methods. As the name suggests, impulse methods compute contact forces based on the required impulse to keep rod segments from penetrating, with an example being the impulse force model by Spillmann et al. spillmann2008adaptive . Although computationally efficient and straightforward to implement, unrealistic visual jittering often occurs when using too large time steps during simulations as the generated forces are handled explicitly choi2021imc . Therefore, impulse methods often must deal with insufficient physical accuracy or use sufficiently small time steps.
Constraintbased optimization methods treat frictional contact as a constrained optimization problem. For example, in Ref kaufman2014adaptive , frictional contact is formatted as constraints according to discrete SignoriniFishera conditions kikuchi1988contact and the maximal dissipation principle goyal1991planar . Other approaches treat frictional contact as a conic optimization problem alart1991mixed ; daviet2011hybrid ; bertails2011nonsmooth ; li2018implicit . Constraintbased methods can often produce physically realistic results but are usually difficult to implement. Furthermore, constraintbased methods require additional handling for the generated constraints, often producing expensive computational costs.
The final contact method type, penalty energy methods, utilize a formulated “contact energy” whose gradient is treated as the contact force. Due to the requirement of a smooth differentiable gradient (and Hessian for implicit formulations), such methods utilize smooth differentiable functions to best approximate the behavior of frictional contact. Such methods have become popular in recent times as they have been shown capable of generating accurate frictional contact li2020incremental ; choi2021imc ; patil2020topological while remaining simple to implement (relative to constraintbased methods) and computationally efficient. Building upon this, we propose Implicit Contact Model (IMC), a fullyimplicit penaltybased contact model for frictional contact based on our previous work in Ref. choi2021imc . We improve upon this iteration by (1) reformulating frictional contact to be fullyimplicit for enhanced physical accuracy, (2) squaring our contact energy term for more rigid contact, (3) changing our smoothly approximated distance formulation to a more stable piecewise analytical formulation, and (4) adding a line search method for increased robustness. Our proposed numerical framework can generate contact for any rodrod contact scenario and can also generate contact for 3D meshes with proper alterations.
In this paper, we choose to showcase our frictional contact model by simulating the novel and difficult contact scenario of flagella bundling flores2005study ; cisneros2008unexpected ; reigh2012synchronization ; maniyeri2014numerical ; hintsche2017polar ; nguyen2018impacts ; lee2018bacterial ; huang2021numerical ; powers2002role , a significant natural phenomena that occurs when microorganisms with multiple flagella swim in fluid (e.g., Escherichia coli and Salmonella typhimurium kim2003macroscopic ). Each flagellum consists of a rotary “head”, a short flexible hook, and a helical filament. By rotating their filaments, these microorganisms can navigate their environments through sophisticated manipulation of the solidfluid interaction between their flexible structures and the surrounding flows. This has led to biomimicry, where flagella have inspired the design of several soft robot locomotion strategies in viscous fluids magdanz2013development ; ye2013rotating ; beyrand2015multi ; son2013bacteria ; rusconi2014bacterial . However, studying the mechanics of flagella is exceptionally challenging due to flagella possessing radii smaller than their optical wavelengths as well as having high rotational speeds jawed2015propulsion .
This highlights the necessity for accurate simulators, the development of which is nontrivial due to the multifaceted problem of having to deal with hydrodynamic interactions, geometrically nonlinear deformations, and realistic contact handling. To the best of our knowledge, this paper is the first to design a fast and accurate simulator that can capture the bundling phenomena of multiple flagella rotating in low Reynolds number fluids. We believe our simulator has promise for aiding the study of flagella bundling as well being an efficient data generator for training bioinspired flagellar robots for datadriven control approaches.
The primary contributions of our work are outlined below.

We propose a fullyimplicit penaltybased frictional contact model that has improved computational efficiency and accuracy compared with our previous work in Ref. choi2021imc .

We formulate a full endtoend framework for the novel and difficult contact scenario of flagella bundling in low Reynolds number fluids, which incorporates a DDGbased simulation framework, our frictional contact framework, and fluidsolid interaction.

We conduct an indepth sidebyside comparison between our proposed method (IMC) with the stateoftheart (IPC) li2020incremental as well as showcase convincing results for the sticking slipping phenomena of friction.
The remainder of the paper is as follows. In Sec. 2, we introduce the DDGbased framework for simulating a flexible rod. Next, in Sec. 3, we formulate our frictional contact model. Simulation results for IMC are then shown in Sec. 4 along with sidebyside comparisons with IPC and analysis of friction performance. Finally, in Sec. 5, we provide concluding remarks as well as potential future research directions. Further details on the fluidsolid interaction model, IMC algorithmic components, and miscellaneous results can be found in Supplementary Information Appendices A, B, and C, respectively. All source code used in this paper is released publicly and can be found at https://github.com/StructuresComp/rodcontactsim.
2 Discrete Elastic Rods
In order to simulate the geometric nonlinear behaviors of flagella in viscous fluids, we utilize the DDGbased framework Discrete Elastic Rods (DER) bergou2008discrete ; bergou2010discrete . As shown in Fig. 2(a), DER expresses the centerline of an elastic rod with discrete nodes: . This results in a total of edges where . Note that for DER, we use subscripts to denote indices for quantities associated with nodes and superscripts for indices for quantities associated with edges. Following this, each edge is described using two orthogonal frames: a reference frame and a material frame as shown in Fig. 2(b). The reference frame is predefined at initial time s. The material frame shares the same director as the reference frame and is obtainable through a twist angle with respect to the reference frame. A total of nodes, each represented by a Cartesian coordinate , and twist angles constitute a total of degrees of freedom: .
To simulate the elastic properties of a rod, we must compute elastic energies as defined by strain. Based off of Kirchhoff’s rod model kirchhoff1859uber , strains can be divided into three categories: stretching, bending, and twisting. Starting off, the stretching strain of an edge is described by
(1) 
where is the undeformed length of edge . From hereafter, quantities with a represent their values in their undeformed state.
Moving forwards, the bending strain for a node is evaluated by a curvature binormal which captures the misalignment between two consecutive edges:
(2) 
Here, , where is the turning angle shown in Fig. 2(b). The material curvatures are the components of the curvature binormal via the directors of the material frame:
(3)  
Finally, the twisting strain for a node is computed as
(4) 
where is the difference between the two consecutive reference frames of the and th edges.
With all strains defined, we can now formulate the stretching, bending, and twisting energy of a discretized elastic rod:
(5)  
where is Young’s Modulus; is the crosssectional area; is the shear modulus; is the polar second moment of area along tangent ; and are moments of inertia along material directors and ; and is the Voronoi length.
Next, the internal forces (for each nodal degree of freedom) and moments (for each twist degree of freedom) can be obtained via the partial derivative of the sum of elastic energies given in Eq. 5 as shown:
(6) 
These quantities result in a sized force vector . Following this, we can write the system of equations of motions as the sum of inertial terms, internal forces, and external forces (e.g. contact, friction, gravity). This results in the equation
(7) 
where is the diagonal mass matrix, is the second derivative of the DOFs with respect to time, is the external force vector, and is the total force. This external force vector will be made up of our contact forces described in Section 3 as well as viscous drag forces from solidfluid interactions described in Supplementary Information A. As Eq. 7 is a rootfinding problem, we use Newton’s method to solve the system of equations to march through time. Note that the Jacobian of the viscous drag forces is unobtainable. Therefore, viscous drag forces are treated explicitly, i.e. their Jacobian is simply ignored. Still, as our contact model is fully implicit, we are able to robustly simulate flagella bundling regardless of the explicitly added forces, as we will soon show.
3 Contact Model Methodology
Moving forwards, we denote the following vector concatenation describing an edgetoedge contact pair: , where to exclude consecutive edges from consideration when enforcing contact. We describe the set of all valid edge combinations as . In future equations, we simply denote the subscriptless as an arbitrary edge combination for clarity. We design a contact energy to increase as the minimum distance between two bodies approaches a contact threshold ( for our application where is the radius of the flagella). With this, the contact energy gradient is used as the contact forces while the contact energy Hessian is used as the contact force Jacobian, where is the contact stiffness which scales the contact forces appropriately to enforce nonpenetration. In the upcoming sections, we will now formulate contact energy , minimum distance between two edges , as well as friction.
3.1 Contact Energy
In the ideal setting, contact energy must satisfy two properties: (1) it is zero for any distance and (2) it is nonzero at exactly distance . A Heaviside step function can essentially describe these properties. Such a function is nonsmooth with a very sudden discontinuous change in value, and therefore, cannot be solved reliably by rootfinding algorithms such as Newton’s method. To remedy this, IPC uses the following energy formulation to smoothly approximate contact:
(8) 
where is the distance tolerance that defines the region for which nonzero forces are experienced. This contact energy approaches when decreasingly approaches and is therefore undefined for the region . Although this barrier formulation allows IPC to strictly enforce nonpenetration, the solver must be careful never to allow any contact pairs in the penetration zone and/or venture into this undefined region during the optimization process. This is ensured by the inclusion of a custom line search method which conservatively sets an upper bound for the Newton update coefficient .
In contrast to this, we design our energy formulation to allow for optimization into the penetrated region, thus expanding the range contact forces are experienced from to . This in turn allows us to take advantage of more aggressive line search methods, which leads to faster convergence for the flagella contact problem. Although this in theory allows our model to be susceptible to penetration, a sufficient contact stiffness remedies this issue. We provide a method that adaptively sets an appropriate stiffness value in Supplementary Information B.3. In addition, to further ensure nonpenetration, we take our previous energy formulation from choi2021imc and square it so that our gradient grows exponentially instead of linearly. In the end, we use the smooth approximation
(9) 
where indicates the stiffness of the energy curve.
We incorporate the piecewise term for two reasons. First, this term equivalently models our energy formulation for the region and has a simpler gradient and Hessian, resulting in computational efficiency. Second, and more importantly, the piecewise term also ensures numeric stability by preventing the exponential term in Eq. 9 from exploding. We show our plotted energy term in Fig. 3(a) for various values. As shown, the energy starts to increase at an exponential rate as decreases towards the contact limit which is shown as 0 here. As decreases, more realistic contact is achieved (enhancing accuracy) in exchange for a stiffer equation (more difficult to converge).
3.2 Computing Distance
As mentioned in li2020incremental , the minimum distance between two edges and can be formulated as the constrained optimization problem
(10) 
where and
represent the contact point ratios along the respective edges. Minimum distance between two edges can be classified into three distinct categories: pointtopoint, pointtoedge, and edgetoedge. As the names suggest, these classifications depend on which points of the edges the minimum distance vector
lies as described by and shown in Fig. 2(c).In our previous work choi2021imc , we altered Lumelsky’s edgetoedge minimum distance algorithm lumelskey1985min_dist (which implicitly computes the values) to be fully differentiable through smooth approximations. In this work, we now change the distance formulation to use piecewise analytical functions as shown below in Eqs. 11, 12, and 13, similar to li2020incremental , as we found more stable performance compared to our smooth formulation despite the nonsmooth Hessian when changing between contact categories.
We now describe the conditions for each contact type classification. First, if lies on the ends of both edges (i.e. both constraints are active), then the distance formulation degenerates to the pointtopoint case which can easily be solved using the Euclidean distance formula,
(11) 
where and are the nodes for first and second edges in contact, respectively.
If only lies on one end of one rod (i.e. only one constraint is active), then the contact type degenerates to pointtoedge. This can be solved as
(12) 
where and are the nodes of the edge for which the minimum distance vector does not lie on an end and is the node of the edge which the minimum distance vector does lie on. Finally, edgetoedge distance (i.e. no active constraints) for the th and th edges can be solved as
(13)  
where indicates a unit vector. With fully defined, this concludes our contact energy formulation. To correctly classify contact pairs, we use Lumelsky’s algorithm to compute values.
3.3 Adding Friction
Similar to before, we model friction according to Coulomb’s friction law, which describes the conditions necessary for two solids to transition between sticking and sliding. This law states that the frictional force is (1) equal to during sliding, (2) is in the region of when sticking, and (3) is independent of the magnitude of velocity. Here, is the friction coefficient and is the normal force experienced by the body.
Let us denote the following equivalencies for clarity: and . Following this, for a contact pair , we can obtain the normal force on the th and th nodes as and , respectively. This in turn allows us to obtain the contact norm vector
(14) 
The direction of friction is then the tangential relative velocity between edges and . To compute this, we must first compute the relative velocities of the edges at the point of contact, which can be done using as shown below:
(15)  
where , and are the velocities of the th, th, th, and th nodes, respectively. The tangential relative velocity of edge with respect to edge can then be computed as
(16) 
where is our friction direction.
Now, we must also make our contact model capable of simulating the transition between sticking and sliding. Coulomb’s law tells us that during static friction and that for sliding friction. Sticking occurs up until the tangential force threshold is surpassed, after which sliding begins. This relation (similar to ideal contact energy) can also be described by a modified Heaviside step function. For the same reasons as before, we replace this step function for another smooth approximation described by
(17) 
where (m/s) is our desired slipping tolerance and is the stiffness parameter. As shown in Fig. 3(b), smoothly scales the friction force magnitude from zero to one as increases from zero. The slipping tolerance describes the range of velocities for which a friction force is experienced. In other words, we consider velocities within this range to be “sticking”.
Finally, the friction experienced by a node for a contact pair can be described as
(18) 
With friction fully defined, we can now formulate the friction Jacobian . Note that due to Eq. 15, our formulation depends on , which means that the gradient is required. We can avoid this computation through the realization that the magnitudes of the contact forces and have an underlying linear relationship with where
(19)  
Therefore, we can obtain by simply solving
(20) 
We can now treat as a function of
, resulting in a simplified chain ruling procedure. Let us denote Eq.
18 as the functional ). The friction Jacobian can then be computed through chain rule as(21) 
This concludes our fully implicit friction scheme. Full psuedocode for the IMC algorithm can be found in Alg. 1.
4 Simulation Results
In this section, we showcase extensive quantitative and qualitative results for IMC. First, we discuss all our simulation parameters. We then conduct a detailed comparison between IMC and the stateoftheart contact handling method: Incremental Potential Contact (IPC) li2020incremental . Afterwards, we showcase comprehensive results concerning friction and display IMC’s ability to simulate the sticking sliding transition.
Model  AIPTS  ATPTS [ms]  Total Iters  Total Run Time [hr]  RTI  

IMC  2  3.00  10.2  0.57  1.82  
3  3.01  21.3  1.19  1.82  
5  3.02  67.5  3.34  1.40  
10  3.12  389.4  22.77  1.22  
IPC  2  4.00  18.75  1.04  N/A  
3  4.00  39.5  2.17  N/A  
5  4.01  95.3  4.68  N/A  
10  4.02  477.47  27.88  N/A 
4.1 Parameters and Setup
In the simulation, we design the flagella as righthanded helical rods manufactured with linear elastic material. We set the material properties as follows: Young’s modulus was set to MPa; Poisson’s ratio was set to 0.5; density of the rod was set to kg/m; the crosssectional radius was set to mm, and the fluid viscosity was set to 0.1 Pas. Here, a Poisson’s ratio of 0.5 was chosen to enforce the flagella to be an incompressible material. The topologies of the flagella are helices with helical radius m, helical pitch m, and axial length m. These parameters were chosen as they best mimic the geometries of biological flagella found in nature jawed2015propulsion ; jawed2017dynamics ; rodenborn2013propulsion ; huang2021numerical .
We explore the bundling phenomena with flagella () where the rotating ends of each flagella is fixed along the axis as shown in Fig. 2(a). These rotating ends are treated as boundary conditions and are spaced out equidistantly so as to form a regular polygon with angles with side length m as shown in Fig 4(c). We set the rotation speed of the flagella ends to rad/s which keeps the Reynolds number in our numerical simulation to be always smaller than , thus satisfying the Stokes flow.
Finally, we discretize each flagella into 68 nodes for a total of 67 edges. We found this discretization to have the best tradeoff between computational efficiency and accuracy. Furthermore, we set the time step size to ms. As the forces generated from our fluid model are handled explicitly, we found 1 ms to be the largest stable time step size before convergence performance became hampered. A distance tolerance of was used for all simulations.
4.2 Comparison between IMC and IPC
Both IMC and IPC were used to simulate 250 seconds of rotation for scenarios with 2, 3, 5, and 10 flagella as shown in Fig 1. As the friction coefficient between structures is usually trivial in viscous fluids, we consider purely contact without friction (). First, a sidebyside visual comparison for is shown for IMC and IPC in Fig. 4(a). For various time steps, we can see that the configurations of the flagella are near identical, indicating that both methods have comparable accuracy. To further study this similarity, we define normalized average difference to measure the difference in flagella nodal configurations between IMC and IPC:
(22) 
The relationship between normalized average difference and time is shown in Fig. 4(d). Here, we can find that the difference between the configurations is quite minimal, further cementing the notion that IMC and IPC possess comparable accuracy.
Where IMC starts to improve upon IPC is in terms of computational efficiency. Detailed metrics for all runs can be seen in Table 1 which showcase the average iteration per time step (AIPTS), average time per time step (ATPTS), total iterations, and total run time. All metrics were recorded using time steps with at least one contact. Here, we can see that IMC was able to converge with less average iterations than IPC for all flagella cases resulting in significant reductions in total run time. These run time improvements are most drastic for and and start to decrease as increases further as the RSS force computation starts to become a bottleneck. Regardless, a clear monotonic decrease can still be seen.
4.3 Friction Example
Although friction is usually negligible in a viscous fluid medium, influence of friction on flagella bundling is still intriguing since the effect of friction can become significant as the environment changes (e.g., granular medium). We assume an imaginary viscous environment where the friction coefficient between structures is nonnegligible. We present simulation data for two flagella () with friction coefficients . For all simulations, a slipping tolerance of was used. All other parameters are kept the same as before.
We first showcase the sticking slipping phenomena with snapshots for in Fig 5. Intuitively, as increases, we also see the amount of sticking increase as well. Convergence results for all friction examples can be seen in Table 2 where average iterations per time step and simulation length are reported. Here, we notice two trends. First, for , the time at which the simulation ends starts to decrease from 250 seconds. This is because is the point at which the flagella become completely tangled as shown in the bottom right frame of Fig 5. As increases past 0.7, the tangling happens earlier and earlier as shown. Furthermore, we observe that the number of average iterations starts to increase as increases. This is in line with our expectations as larger values result in greater sticking.
AIPTS  Total Iterations  Sim End [sim s]  

0.1  3.01  250  
0.2  3.01  250  
0.3  3.61  250  
0.4  4.89  250  
0.5  6.67  250  
0.6  8.71  250  
0.7  14.47  235.89  
0.8  14.16  180.98  
0.9  11.1  142.72  
1.0  11.65  135.32 
5 Conclusion
In this paper, we introduced an improved version of our fullyimplicit and penaltybased frictional contact method, Implicit Contact Model. To test the performance of our contact model, we formulated an endtoend simulation framework for the novel and difficult contact scenario of flagella bundling in viscous fluid. For this contact problem, we showed that IMC has comparable accuracy to the stateoftheart while maintaining faster convergence. Furthermore, we showcased convincing frictional results in an imaginary viscous environment where friction is nonnegligible.
For future work, we wish to improve upon the stability and robustness of our friction model. Despite the implicit formulation, the number of iterations necessary to converge starts to increase as
increases. Another interesting avenue of research is the use of deep learning to learn physicsbased dynamics for simulation. Neural networks, when properly trained, have been known to be able to generate nearly identical outputs as numerical simulations while achieving orders of magnitude reduction in computation. Thus, utilizing the computational efficiency and differentiability of neural networks while maintaining physical realism is a promising direction.
6 Acknowledgements
We are grateful for financial support from the National Science Foundation (NSF) under award number CMMI2101751 and IIS1925360. M.K.J. is grateful for support from NSF (CAREER2047663, CMMI2053971).
Appendix A Fluid Model
In this section, we formulate the relationship between velocity and hydrodynamic forces for each node in a discretized rod. We utilize Regularized Stokeslet Segments (RSS) cortez2018regularized method to capture the viscous drag forces exerted on a slender rod moving in a viscous fluid as shown in Fig. 2(c). The Stokeslet, the primary Green’s function of Stokes flow, represents solutions for the Stokes equations as a superposition of fundamental solutions. The methods of regularized Stokeslets have been implemented in numerous studies such as selfpropelled microorganisms cisneros2010fluid spagnolie2012hydrodynamics , hyperactivated sperm motility olson2011coupling simons2015fully , and flagella bundling flores2005study huang2021numerical .
For a particular choice of regularization cortez2018regularized , the relationship of the velocity at an evaluation point and the corresponding regularized force exerted on a point is the regularized Stokeslet:
(23) 
where is the fluid viscosity; ; , and is the regularization parameter.
For an edge in a viscous fluid, a point along this edge can be defined by , where . We assume a linear density of the force applied on the edge: . With this, we obtain the relationship between the velocity at and the linear force density as
(24) 
where is the length of the edge . Through integration, we can then rewrite Eq. 24 as
(25) 
where the coefficients are
(26)  
and the sequences of are
(27)  
Recall that in DER a flexible rod constitutes nodes and edges each with length . Therefore, we can formulate hydrodynamic forces for the whole rod through the summation of Eq. 25 for each segment, resulting in
(28) 
where coefficient matrices and are
(29)  
Finally, by using Eqs. 28 and 29, we can describe the balance between the velocity and the force density along a discretized rod as the linear system
(30)  
where (the velocity of the liquid at nodal position has the same speed of the rod due to nonslip boundary); is the force density; is the hydrodynamic force vector, and is the Voronoi length at node . We assume the density of the fluid to be equal to the density of the rod so that buoyant forces are negligible. As stated in Ref cortez2018regularized , we choose the regularization parameter based on the value of where is the Voronoi length.
Appendix B IMC Algorithmic Components
b.1 Scaling
Since the radius of the rod may be very thin (in the range of a few mm), collision detection and gradient and Hessian computations may be susceptible to floating point inaccuracies. To mitigate this, we scale all nodes by a factor of for numeric stability. We denote and and as a result, normalize Eq. 9 to become
(31) 
Accounting for this scaling in the gradient and Hessian computation can easily be taken care of through the product and chain rule which ends in us multiplying the gradient by and the Hessian by . With this normalization in mind, we choose a stiffness for a desired through the relation
(32) 
b.2 Collision Detection
For collision detection, we simply obtain the set of all edge combinations whose minimum distance is less than , resulting in the contact set
(33) 
As this operation can be quite computationally expensive, we instead compute a candidate set at the beginning of each time step where by a large enough margin. We then compute the actual contact set from the candidate set at the start of each iteration, significantly reducing computational cost.
Note that if is not set large enough, certain edge combinations not belonging to the initial set may enter the contact zone or even penetrate by the end of the optimization. Our energy formulation in Eq. 9 is capable of dealing with this as minor penetrations do not lead to simulation failure and will be remedied in the next time step. This is in contrast to IPC, which may require a significantly larger and/or more robust collision checking during each iteration of the optimization process.
b.3 Contact Stiffness
As a penalty method, a contact stiffness is used to scale the contact force and Jacobian. An appropriate value must be used to ensure that penetration (caused by too low of a value) or excessive hovering (caused by too large of a value) does not occur. First, let us denote the set of all and node indices of the contact set edge combinations as the set . Then, to generate an appropriate scaling for the contact stiffness, we compute the norm of the sum of forces (minus contact and friction) experienced by each node in the contact set . We can denote these forces as the set
(34) 
The maximum force magnitude of these forces can then be used to determine the contact stiffness
(35) 
where is a constant scaling factor. In all our experiments, we set to be . The intuition behind this contact stiffness formulation is to achieve nonpenetration through force equilibrium. Furthermore, by using the maximum of , if nonpenetration can be achieved for the edge with the largest value in , then this value should be large enough to prevent penetration for all other contact pairs as well.
b.4 Line Search
Once the internal forces (e.g., bending force, twisting force, and stretching force), external forces (e.g., viscous dragging force and contact forces), and their respective Jacobians are computed, we can simply use Newton method to find the solution of the equations of motion. However, due to the high nonlinearity of the governing equations, convergence for Newton’s method may suffer without a line search method. To rectify this, we perform GoldsteinPrice line search in the Newton direction to ensure that the square of the Euclidean norm of total force in Eq. 7 decreases.
We design an inner loop for each Newton iteration where we deploy the line search algorithm. This inner loop returns an optimal search step size until is smaller than a certain tolerance or until a maximum number of iterations is reached. As mentioned in Sec. 3.1, our energy formulation allows us to use a more aggressive line search strategy compared to Ref li2020incremental , resulting in larger search step sizes and faster convergence. Pseudocode for the line search method can be found in detail in Algorithm 2.
b.5 Full Pseudocode
Appendix C Miscellaenous Results
c.1 Stability of IMC and IPC
We plot the cumulative iteration and computation time for both IMC and IPC for all flagella cases with respect to simulation time in Fig. 6. When observing the plots, we can see that both methods have near linear growth rates (a sign of stability) for all metrics with IMC consistently possessing higher computational efficiency as indicated by the lower slopes.
c.2 Normalized Propulsive Force vs. Friction Coefficient
We also look into the influence of on the propulsive forces generated by the flagella bundling. With the ends of flagella clamped, we can refer to the propulsive force as the external forces experienced by the clamped nodes. With this, we formulate a normalized propulsive force as
(36) 
where is the bending stiffness. We plot normalized propulsive force with respect to for two scenarios. In Fig. 7(a), the relationship between normalized propulsive force and is given for averaged over 250 seconds. We do this as simulations for high frictional cases end before the final time s due to the extreme tightening of rotating ends caused by sticking forces. For this plot, we observe that the propulsive force is nearly identical for and then starts to exponentially increase as increases past 0.3. Let us explain this behavior through physics. First, as grows, the sticking surface between flagella grows, eventually resulting in “one” flagellum. This essentially results in a single flagellum with a higher effective stiffness, which translates to higher thrust forces as a single flagellum has been shown to have higher propulsive efficiency compared to multiple flagella with equal cumulative volume riley2018swimming .
In Fig. 7(b), the same is done for averaged over 135 seconds as this is the earliest sim end time (). Here, we also observe that propulsive force increases as increases. The abnormal change between and can be attributed to insufficient data as seconds is nowhere near the simulation end for these two scenarios.
c.3 Solver and Computational Tools
On the computational end, we use the sparse matrix solver oneMKL PARDISO from Intel’s oneAPI Math Kernel Library MKL
and optimized all operations by compiling Eigen to use MKL as the BLAS and LAPACK backend. All necessary gradients and Hessians were computed symbolically and created into callable functions optimized through LLVM using the opensource C++ library SymEngine
Symengine . All experiments were run singlethreaded using an AMD Ryzen Threadripper 2990WX 32core processor. Although parts of our code can be improved by multithreading, we decided to forgo this as solving the dense matrix in Eq. 30 soon becomes a bottleneck.References
 (1) M. K. Jawed, F. Da, J. Joo, E. Grinspun, P. M. Reis, Coiling of elastic rods on rigid substrates, Proceedings of the National Academy of Sciences 111 (41) (2014) 14663–14668.
 (2) M. K. Jawed, P. M. Reis, Pattern morphology in the elastic sewing machine, Extreme Mechanics Letters 1 (2014) 76–82.
 (3) M. K. Jawed, P.T. Brun, P. M. Reis, A geometric model for the coiling of an elastic rod deployed onto a moving substrate, Journal of Applied Mechanics 82 (12) (2015) 121007.
 (4) C. Baek, A. O. SagemanFurnas, M. K. Jawed, P. M. Reis, Form finding in elastic gridshells, Proceedings of the National Academy of Sciences 115 (1) (2018) 75–80.
 (5) J. Panetta, M. KonakovićLuković, F. Isvoranu, E. Bouleau, M. Pauly, Xshells: A new class of deployable beam structures, ACM Transactions on Graphics (TOG) 38 (4) (2019) 83.
 (6) C. Baek, P. M. Reis, Rigidity of hemispherical elastic gridshells under point load indentation, Journal of the Mechanics and Physics of Solids 124 (2019) 411–426.
 (7) A. Goriely, S. Neukirch, Mechanics of climbing and attachment in twining plants, Physical review letters 97 (18) (2006) 184302.
 (8) M. Jawed, P. Dieleman, B. Audoly, P. M. Reis, Untangling the mechanics and topology in the frictional response of long overhand elastic knots, Physical review letters 115 (11) (2015) 118302.
 (9) B. Audoly, N. Clauvelin, S. Neukirch, Elastic knots, Physical Review Letters 99 (16) (2007) 164301.
 (10) V. P. Patil, J. D. Sandt, M. Kolle, J. Dunkel, Topological mechanics of knots and tangles, Science 367 (6473) (2020) 71–75.
 (11) V. P. Patil, Ž. Kos, M. Ravnik, J. Dunkel, Discharging dynamics of topological batteries, arXiv preprint arXiv:2002.11822 (2020).
 (12) P. Grandgeorge, C. Baek, H. Singh, P. Johanns, T. G. Sano, A. Flynn, J. H. Maddocks, P. M. Reis, Mechanics of two filaments in tight orthogonal contact, Proceedings of the National Academy of Sciences 118 (15) (2021).
 (13) M. K. Jawed, N. K. Khouri, F. Da, E. Grinspun, P. M. Reis, Propulsion and instability of a flexible helical rod rotating in a viscous fluid, Physical review letters 115 (16) (2015) 168101.
 (14) M. K. Jawed, P. M. Reis, Deformation of a soft helical filament in an axial flow at low reynolds number, Soft Matter 12 (6) (2016) 1898–1905.
 (15) M. Jawed, P. M. Reis, Dynamics of a flexible helical filament rotating in a viscous fluid near a rigid boundary, Physical Review Fluids 2 (3) (2017) 034101.
 (16) M. Bergou, M. Wardetzky, S. Robinson, B. Audoly, E. Grinspun, Discrete elastic rods, ACM transactions on graphics (TOG) 27 (3) (2008) 63.
 (17) M. Bergou, B. Audoly, E. Vouga, M. Wardetzky, E. Grinspun, Discrete viscous threads, ACM Transactions on Graphics (TOG) 29 (4) (2010) 116.
 (18) J. Spillmann, M. Teschner, An adaptive contact model for the robust simulation of knots, in: Computer Graphics Forum, Vol. 27, Wiley Online Library, 2008, pp. 497–506.
 (19) A. Choi, D. Tong, M. K. Jawed, J. Joo, Implicit Contact Model for Discrete Elastic Rods in Knot Tying, Journal of Applied Mechanics 88 (5) (03 2021).
 (20) D. M. Kaufman, R. Tamstorf, B. Smith, J.M. Aubry, E. Grinspun, Adaptive nonlinearity for collisions in complex rod assemblies, ACM Transactions on Graphics (TOG) 33 (4) (2014) 1–12.
 (21) N. Kikuchi, J. T. Oden, Contact problems in elasticity: a study of variational inequalities and finite element methods, SIAM, 1988.
 (22) S. Goyal, A. Ruina, J. Papadopoulos, Planar sliding with dry friction part 2. dynamics of motion, Wear 143 (2) (1991) 331–352.
 (23) P. Alart, A. Curnier, A mixed formulation for frictional contact problems prone to newton like solution methods, Computer methods in applied mechanics and engineering 92 (3) (1991) 353–375.