A Fully Implicit Method for Robust Frictional Contact Handling in Elastic Rods

by   Dezhong Tong, et al.

Accurate frictional contact is critical in simulating the assembly of rod-like structures in the practical world, such as knots, hairs, flagella, and more. Due to their high geometric nonlinearity and elasticity, rod-on-rod contact remains a challenging problem tackled by researchers in both computational mechanics and computer graphics. Typically, frictional contact is regarded as constraints for the equations of motions of a system. Such constraints are often computed independently at every time step in a dynamic simulation, thus slowing down the simulation and possibly introducing numerical convergence issues. This paper proposes a fully implicit penalty-based frictional contact method, Implicit Contact Model (IMC), that efficiently and robustly captures accurate frictional contact responses. We showcase our algorithm's performance for the challenging and novel contact scenario of flagella bundling in fluid medium, a significant phenomenon in biology that motivates novel engineering applications in soft robotics. In addition to this, we offer a side-by-side comparison with Incremental Potential Contact (IPC), a state-of-the-art contact handling algorithm. We show that IMC possesses comparable accuracy to IPC while converging at a faster rate for the flagella bundling case.


page 1

page 2

page 3

page 4


Factory: Fast Contact for Robotic Assembly

Robotic assembly is one of the oldest and most challenging applications ...

Linear-frictional contact model for 3D discrete element simulations of granular systems

The linear-frictional contact model is the most commonly used contact me...

Codimensional Incremental Potential Contact

We extend the incremental potential contact (IPC) model for contacting e...

Simulation of Parallel-Jaw Grasping using Incremental Potential Contact Models

Soft compliant jaw tips are almost universally used with parallel-jaw ro...

Reduced-order modeling of a sliding ring on an elastic rod with incremental potential formulation

Mechanical interactions between rigid rings and flexible cables are wide...

Critically fast pick-and-place with suction cups

Fast robotic pick-and-place with suction cups is a crucial component in ...

An Efficient Model Order Reduction Scheme for Dynamic Contact in Linear Elasticity

The paper proposes an approach for the efficient model order reduction o...

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 rod-like 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 .

Figure 1: Rendered snapshots of flagella bundling with varying amounts of flagella. Rows contain (a) , (b) , (c) , and (d)

flagella. Each column indicates the flagella configuration at the moment of time indicated in the top row.

With real-world experiments being costly and tedious to implement, the need for accurate physics-based numerical simulations is prevalent. Such simulations not only allow for advanced mechanics-based 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 geometry-based (DDG-based) 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, constraint-based 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.

Constraint-based optimization methods treat frictional contact as a constrained optimization problem. For example, in Ref kaufman2014adaptive , frictional contact is formatted as constraints according to discrete Signorini-Fishera conditions kikuchi1988contact and the maximal dissipation principle goyal1991planar . Other approaches treat frictional contact as a conic optimization problem alart1991mixed ; daviet2011hybrid ; bertails2011nonsmooth ; li2018implicit . Constraint-based methods can often produce physically realistic results but are usually difficult to implement. Furthermore, constraint-based 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 constraint-based methods) and computationally efficient. Building upon this, we propose Implicit Contact Model (IMC), a fully-implicit penalty-based 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 fully-implicit 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 rod-rod 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 micro-organisms 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 micro-organisms can navigate their environments through sophisticated manipulation of the solid-fluid 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 bio-inspired flagellar robots for data-driven control approaches.

The primary contributions of our work are outlined below.

  • We propose a fully-implicit penalty-based frictional contact model that has improved computational efficiency and accuracy compared with our previous work in Ref. choi2021imc .

  • We formulate a full end-to-end framework for the novel and difficult contact scenario of flagella bundling in low Reynolds number fluids, which incorporates a DDG-based simulation framework, our frictional contact framework, and fluid-solid interaction.

  • We conduct an in-depth side-by-side comparison between our proposed method (IMC) with the state-of-the-art (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 DDG-based 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 side-by-side 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 fluid-solid 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/rod-contact-sim.

2 Discrete Elastic Rods

Figure 2: (a) Discrete schematic of a flagellar elastic rod. Nodes , , and the edge between and is clamped along the dashed centerline and rotated with an angular velocity

. The rest of the nodes constitute the helical flagellum which revolves around the centerline. (b) A zoomed in snapshot of two edges showcasing their reference frame, material frame, turning angles, and twist angles. (c) Illustration of two edges approaching contact. The green dots showcase the nodes of the edges while the green dashed lines denote the centerlines of the edges. The red dashed line denotes the vector

whose norm is the minimum distance between the edges. is connected to edges and by and where . As approaches the contact threshold , repulsive forces increase at an exponential rate, thus enforcing non-penetration.

In order to simulate the geometric nonlinear behaviors of flagella in viscous fluids, we utilize the DDG-based 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


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:


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:


Finally, the twisting strain for a node is computed as


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:


where is Young’s Modulus; is the cross-sectional 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:


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


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 solid-fluid interactions described in Supplementary Information A. As Eq. 7 is a root-finding 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 edge-to-edge 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 non-penetration. In the upcoming sections, we will now formulate contact energy , minimum distance between two edges , as well as friction.

Figure 3: Plots for the approximation functions in (a) Eq. 9 and (b) Eq. 17 with varying tolerance values. Note that some of the tolerances displayed are unrealistically large for clarity.

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 non-zero at exactly distance . A Heaviside step function can essentially describe these properties. Such a function is non-smooth with a very sudden discontinuous change in value, and therefore, cannot be solved reliably by root-finding algorithms such as Newton’s method. To remedy this, IPC uses the following energy formulation to smoothly approximate contact:


where is the distance tolerance that defines the region for which non-zero 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 non-penetration, 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 non-penetration, 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


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


where and

represent the contact point ratios along the respective edges. Minimum distance between two edges can be classified into three distinct categories: point-to-point, point-to-edge, and edge-to-edge. 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 edge-to-edge 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 non-smooth 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 point-to-point case which can easily be solved using the Euclidean distance formula,


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 point-to-edge. This can be solved as


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, edge-to-edge distance (i.e. no active constraints) for the -th and -th edges can be solved as


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.

1 Function IMC():
         // compute velocity
         // Eq.9
         // scale by contact stiffness
         // Eq.18
         // Eq.21
2       return
Algorithm 1 Implicit Contact Model

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


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:


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


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


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


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


Therefore, we can obtain by simply solving


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


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 state-of-the-art 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.

Figure 4: Rendered snapshots for flagella simulated by (a) IMC and (b) IPC. We can observe that there is great qualitative agreement between both methods at the shown time steps. (c) A top down visualization of boundary conditions applied to the highest nodes (filled in red circles) of each flagella as well as the angular rotation applied to them. The larger hollow red circles represent the rest of the helical flagella. (d) The norm of the average difference in the nodal positions for the flagella simulated by IMC and IPC with respect to time.
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
Table 1: IMC vs. IPC li2020incremental run time data. Simulations are run for a total of 250 seconds with a time step size of ms and a rotation speed of rad/s. The contact model used can be seen in the far left column. Next to this, indicates the number of flagella. AIPTS stands for average iterations per time step. ATPTS stands for average time per time step. Total Iters indicates the total number of Newton’s iterations that were necessary to complete the simulation. The Total Run Time is the total computational time to completion. Finally, RTI stands for run time improvement and is the ratio of improvement between IMC’s and IPC’s total run time.

4.1 Parameters and Setup

In the simulation, we design the flagella as right-handed 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 cross-sectional 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 trade-off 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 side-by-side 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:


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.

Figure 5: Rendered snapshots for with varying friction coefficients. Each column indicates a moment in time as indicated by the time stamp in the top row. The first row shows the frictionless case as a baseline. The second row has where minor sticking can be observed as the point at where the flagella no longer contact is higher than the frictionless case. Still, still has plenty of slipping allowing the flagella to not become coiled. As we increase to 0.7 in the third row, we can see the amount of sticking increase, ultimately resulting in the flagella becoming completing coiled.

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 non-negligible. 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
Table 2: Friction results for varying friction coefficients. AIPTS stands for average iterations per time step. Total Iterations indicate the total number of Newton’s method iterations that were necessary to complete the simulation. Sim End indicates the total simulated time. All simulations were set to run for 250 seconds. As can be seen, simulations with end earlier due to excessive tangling of the flagella.

5 Conclusion

In this paper, we introduced an improved version of our fully-implicit and penalty-based frictional contact method, Implicit Contact Model. To test the performance of our contact model, we formulated an end-to-end 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 state-of-the-art while maintaining faster convergence. Furthermore, we showcased convincing frictional results in an imaginary viscous environment where friction is non-negligible.

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 physics-based 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 CMMI-2101751 and IIS-1925360. M.K.J. is grateful for support from NSF (CAREER-2047663, CMMI-2053971).

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 self-propelled 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:


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


where is the length of the edge . Through integration, we can then rewrite Eq. 24 as


where the coefficients are


and the sequences of are


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


where coefficient matrices and are


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


where (the velocity of the liquid at nodal position has the same speed of the rod due to non-slip 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


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


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


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


The maximum force magnitude of these forces can then be used to determine the contact stiffness


where is a constant scaling factor. In all our experiments, we set to be . The intuition behind this contact stiffness formulation is to achieve non-penetration through force equilibrium. Furthermore, by using the maximum of , if non-penetration 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.

Parameters :    // Initial interval for
Input:   // DOFs from DER
Output:   // Newton search magnitude
4 Function LineSearch():
5       iter
7       success False
9       while success is False do
10             if   then
11                   success True
12             else if   then
14             else
16             if  small value or iter iterMax then
17                   success True
19             iter iter + 1
20       end while
21      return
Algorithm 2 LineSearch

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 Goldstein-Price 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

Parameters : , tolerance
Require : boundary conditions
1 Function FlagellaSim():
2       Guess
         // Eq. 30
6       while  tolerance do
               // Eq. 6
7             if  then // run only on first iter
                     // Sec. B.2
                     // Sec. B.3
10             for  do
                     // Alg. 1
14             end for
               // Downsize to only include free DOFs
               // Solve
               // Alg. 2
               // update error
20       end while
23       return
Algorithm 3 Flagella Simulation

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.

Figure 6: Cumulative iteration and computation time with respect to simulation time.

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


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 .

Figure 7: Relationship between normalized propulsive force and friction coefficient . Normalized propulsive force for one specific friction coefficient is obtained as the average force from initial time s to final time (a) s and (b) s.

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 open-source C++ library SymEngine 

Symengine . All experiments were run single-threaded using an AMD Ryzen Threadripper 2990WX 32-core processor. Although parts of our code can be improved by multi-threading, we decided to forgo this as solving the dense matrix in Eq. 30 soon becomes a bottleneck.


  • (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. Sageman-Furnas, 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, X-shells: 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.