Implicit Contact Model for Discrete Elastic Rods in Knot Tying

Discrete differential geometry-based simulations have shown surprisingly successful performance in simulating slender structures, e.g., rods [1–4], viscous threads [2,5], ribbons [6], and plates/shells [7–9]. Discrete elastic rod (DER) algorithm [2,4]—originally developed in the computer graphics community to simulate hair, fur, and other filamentary structures in movies—has been borrowed by the engineering community to solve a variety of problems involving deployment of rods [10–12], propulsion of bacterial flagella [13–15], elastic gridshells [16–18], and self-assembly of carbon nanotubes [19]. One of the most interesting and challenging problems in these tasks is the simulation of knots, which introduces complex and intricate patterns of self-contact within an elastic rod. The mechanics of knots [20–26] tied in elastic rods is an intricate interplay of the elastic forces and friction. Therefore, when simulating knots, a fast, stable, and physically accurate contact model is desired.

Contact handling methods can generally be divided into three categories: impulse methods [27], constraint-based methods [28], and penalty methods [22,29]. DER, which implements a variation of Kirchhoff’s rod model [30], has been used to handle contact during knot tying using a contact method proposed by Spillmann and Teschner [27]. This model resolves contact by computing the contact forces that will exactly lead to the desired, collision-free state. Although computationally efficient, unrealistic visual jittering during knot tying occurs for sufficiently large time-steps due to its explicit nature. Kaufman et al. [28] proposed a method capable of simulating large assemblies of elastic rods by adaptively adjusting its degree of nonlinearity. This method formulates frictional contact using discrete Signorini–Fischera conditions [31] and the maximal dissipation principle [32]. By adaptively incorporating sufficient nonlinearity into the collision response, more physically realistic results are produced compared to impulse methods. More recently, Li et al. [29] proposed an implicit time-stepping method that utilizes smooth barrier functions to simulate contact between deformable objects and induce deformations. Here, approximated functions are introduced to smooth the contact energy and frictional responses between two elements. Finally, Patil et al. [22] formulated contact forces in knots by using strain potential, while friction was formulated as a damping force. Although robust, many of these state-of-the-art methods [28,29] are difficult to implement, computationally expensive, and can be often overkill for simulating knots. This can be attributed to the fact that they are formulated in a way that robustly handles difficult contact scenarios such as high velocity impacts, which are often absent in knot tying. To combat this, we formulate an algorithm that is efficient, intuitive, easy to implement, and is specially catered toward knot tying. We then compare simulation results with the contact method proposed by Spillmann and Teschner (SPT) [27] as it is most comparable in implementation complexity and computation cost.

In this article, we introduce implicit contact model (IMC), a contact handling method for DER simulations [1,2] where the contact forces are handled implicitly. DER discretizes the elastic rod into a number of “nodes.” Two consecutive nodes are connected by an “edge” To deal with contact, we define a contact energy, which will be used to derive the normal contact forces (responsible for enforcing nonpenetration) and Coulombic frictional forces. Instead of formulating this contact energy as a function of the distance between nodes, we formulate it as a function of the minimum distance between two edges, which results in more visually and physically realistic results. Figure 1 presents snapshots from our simulations of the knot tying process, where the two free ends of the “tails” of the knots are pulled. Throughout this article, we consider open overhand knots [20]; such knots can be described by the unknotting number n, related to the number of turns in the “braid” of the knot. Figures 1(a)–1(d) show knots with n = 1, …, 4. Interestingly, when the knots with n = 3 and n = 4 are sufficiently tight, they undergo snap-through buckling and the “loop” of the knot suddenly transitions from a near-circular shape to a distorted configuration. The simulation can reliably capture this behavior.

This article is organized as follows. In Sec. 2, we discuss the methodology of the proposed version of DER with IMC as the contact model. Next, in Sec. 3, we undertake both theoretical validation of the IMC contact model and compare it to a well-established preexisting contact model in terms of pull force accuracy and runtime. Finally, conclusive remarks as well as potential future research directions are discussed in Sec. 4.

In this section, first, we very briefly discuss the DER model (for a more in-depth introduction, please refer to Ref. [4]). Then, we formulate the IMC contact model and its intuitive integration into the existing DER algorithm. Finally, hyperparameters are explained in terms of their effect on algorithmic performance.

Under the DER framework, the centerline of a rod is discretized into N nodes, x_i (1 ≤ i ≤ N) which corresponds to N − 1 edges, e_i (1 ≤ i ≤ N − 1), where e_i = x_i+1 − x_i. Each edge has an orthonormal adapted reference frame, ${d1i,d2i,ti}$⁠, and a material frame, ${m1i,m2i,ti}$ as shown in Fig. 2(a). The reference frame is updated through parallel transport in time [1,2], while the material frame can be obtained from the scalar twist angle θⁱ with respect to the reference frame. Following this, for a given rod, the Cartesian coordinates of each node as well as the twist angle for each edge results in 4N − 1 total degrees-of-freedom (DOF), which can be represented by the following vector q = [x₁, θ₁, x₂, θ₂, …, x_N−1, θ_N−1, x_N]^T; here, ^T denotes transposition operator. In this representation, the location of the ith node, x_i, corresponds to the 4i − 3, 4i − 2, and 4i − 1th entries of q.

Normally, if the gradient of an external force, $∂Fiext/∂qj$⁠, cannot be analytically evaluated, this term is omitted during Newton iterations and that external force is considered “explicitly” (Euler forward). This generally requires a smaller time-step size dt, leading to larger computation time. Later, we will show that for IMC, the contact and friction Jacobians are analytically obtainable. Implicit treatment of contact and friction is a key contribution of this article.

Referring to Fig. 2(b), denote $xi,xi+1,xj,xj+1∈R3$ as the Cartesian nodal coordinates of the ith and jth edges in a rod configuration. Next, denote an edge “combination” as the following vector concatenation: x_ij := (x_i, x_i+1, x_j, x_j+1). We denote the set of all valid edge combinations as $X$⁠; two consecutive edges are always in contact, and those combinations are not included in this valid combination $X$⁠.

From here, an arbitrary edge combination is denoted as simply x and all edge combinations are assumed to be valid: $x∈X$⁠. The contact energy E is then expressed as a differentiable analytical expression, which takes the four nodes of the two contacting edges as inputs, $E(x):R12→R1$⁠. Under this formulation, we can see that the proposed contact energy is only dependent on the nodal coordinates of the discretized rod and not on the twist angles of the edges, θ, as the contact forces are computed based on the minimum distance, Δ, between two contacting edges.

Following this, to calculate Δ, we utilize an efficient and accurate algorithm for computing the minimum distance between two finite line segments in N dimensions proposed by Lumelsky [33]. Originally a piecewise function, the algorithm is then modified to become a twice differentiable smooth approximation. By using this completed expression, we can then obtain the negative gradient of the energy $−∇E(xij)$ as well as the negative Hessian $−∇2E(xij)$⁠, which are then used to evaluate F_c and J_c in Eqs. (5) and (6). The gradient of E(x) produces contact forces that act in the direction of the contact normal and whose magnitude varies with Δ, which results in physically realistic forces when dealing with rod–rod contact. Finally, as this method is essentially a penalty method, a stiffness parameter k_c is then used to scale F_c and J_c appropriately.

By producing the contact forces in this way, dynamic friction can be calculated according to Coulomb’s friction law. In the past, previous methods [27,34] have been unable to simulate Coulomb friction due to the inability of obtaining its Jacobian. As we have access to the normal contact force Jacobian, we can calculate Coulombic friction forces as well as the friction Jacobian. By having access to the Jacobians of the normal and friction force, our contact model reliably converges even in complex contact states. As the cost of producing the Jacobian is relatively high and leads to having to solve a nonbanded matrix, we introduce a hybrid approach that ensures computational speed by only calculating the Jacobian when necessary.

To model contact energy, we compute the minimum distance Δ between two edges and feed this value into a smooth inverse ReLU function whose origin is based on the contact distance 2h as shown below where h is the cross-sectional radius of the rod and ce_k is a stiffness term that determines how sharp the curve is.

Lumelsky’s algorithm produces the minimum distance between two line segments in $RN$ and contains three noncontinuous components; we can eliminate one of them by simply taking the assumption that no edge can be reduced to a point. In other words, each edge must have a finite length greater than zero. With this condition eliminated, we now briefly layout the simplified min-distance algorithm for an arbitrary edge combination x_ij. Note that here, we simply go over the steps of the algorithm. For further intuition on how exactly the algorithm is computing the min-distance, please refer to the original paper [33]. For notation, · is scalar multiplication and $∙$ is the dot product.

The last noncontinuous component of the algorithm lies in the piecewise function in Eq. (10), which is actually left noncontinuous. Although this introduces a piecewise function into the expression, this does not hurt convergence for the following reasons. First, it should be noted that λ will almost never equal exactly 0 due to floating point arithmetic. Therefore, the piecewise function is only required to prevent simulation crashes during simulation starts with perfectly parallel rod configurations. Furthermore, the numeric stability of this algorithm is maintained as whenever λ approaches zero, the numerator S₁ · D₂ − S₂ · R also approaches zero by a similar magnitude, effectively avoiding numeric overflow problems [33]. In terms of performance, we have found that the produced Hessian is an excellent indicator of gradient direction when approaching or passing through parallel configurations when validating against finite difference for a wide variety of edge configurations.

Here, we see that for any arbitrary edge combination x, the produced force $−∇E(x)$ will be a vector of size 12 consisting of four concatenated three-dimensional contact force vectors for every node making up x. These 12 elements contribute to the 12 entries of the (4N − 1)-sized F_c vector located at the following positions: 4i − 3, 4i − 2, 4i − 1, 4(i + 1) − 3, 4(i + 1) − 2, 4(i + 1) − 1, 4j − 3, 4j − 2, 4j − 1, 4(j + 1) − 3, 4(j + 1) − 2, and 4(j + 1) − 1. Once the contact forces are computed for every contacting edge combination during a time-step, the force values are added to F_c and then incorporated into DER.

To obtain the Hessian, we simply take Eq. (17) and differentiate once again to obtain $−∇2E(x)$⁠. Once obtained, the Hessian is added to the (4N − 1 × 4N − 1)-sized Jacobian matrix J_c in a similar manner. The derivation of the Hessian can be done by using the product rule and its derivation is left out for brevity. See Sec. 5 for source code implementing these expressions.

Just as we obtained contact force vectors of size $R12$⁠, we produce frictional forces in a similar manner where given an edge combination x_ij, we compute a friction force vector $Ffrij∈R12$ according to Coulomb’s dynamic friction law. For each edge combination, this 12-sized vector is added to the appropriate entries of F_fr (Eq. (5)).

Once these friction force vectors are computed for every contacting edge combination during a time-step, we can then compute the friction force Jacobian matrix $Jfrij$ as well. These are then added to F_fr and J_fr in exactly the same way as the contact energy gradient and Hessian as described in Sec. 2.3.

It should be noted that several simplifications were made for the friction model. First, the relative velocities were computed using the midpoint of the edges rather than the contact points. Likewise, the friction forces were evenly distributed rather than being dependent on the contact points. This was done as it greatly simplifies the friction force Jacobian, leads to improved convergence, and does not have any noticeable effects so long as the rod is sufficiently discretized.

Furthermore, we treat friction semi-explicitly by using the known velocities from the previous time-step. This allows the friction direction to remain constant during Newton iterations, which improves convergence considerably. Although a fully implicit scheme is possible, computing the necessary contact Hessian on every iteration is costly and the overall speed of the algorithm greatly benefits from this formulation.

In terms of limitations, this method clearly does not enforce static friction and so should only be used for continuous sliding scenarios such as knot tying. Furthermore, friction occurring due to the rod twist θ is not modeled. When an edge undergoes enough twist and is in contact with a receiving edge, friction forces occur slightly off the centerline of this receiving edge. As our model assumes that all friction occurs precisely on the centerline and formulates all contact only using x, such friction–twist coupling is neglected. Finally, the produced friction forces can possibly overtake the pull forces when the pull speed is very low leading to unrealistic sliding in the opposite direction. This must be remedied by pulling at a sufficiently high speed.

In addition to the force and Jacobian generation, there are several additional steps to the IMC algorithm that are explained in this section as well as several hyperparameters that must be properly tuned for optimal performance and convergence which are listed as follows:

1. δ, the collision limit,

2. k_c, the contact stiffness,

3. ce_k, the contact energy stiffness, and

4. ω, the number of iterations before the hybrid algorithm computes the Jacobian.

The collision limit δ must be chosen carefully as a higher δ value results in additional computation due to more qualifying edge combinations, whereas a δ value that is too low will produce nonsmooth gradients that hamper convergence.

A good way to determine a proper δ value is to observe the plotted contact energy function from Eq. (21) for a chosen ce_k value. By observing the contact energy curve, the point at which the generated gradients are ≈ 0 can be found when the slope of the curve is nearly flattened out. Choosing a δ value that encompasses this region ensures that the generated gradients are sufficiently smooth. An example of this process can be found in Fig. 4.

The stiffness of the contact forces is determined by the contact energy stiffness value ce_k. As this value becomes higher, the region above the contact surface at which two edges experience a force decreases. This leads to stiff contact, which is more physically accurate as realistically the contact energy should be zero whenever Δ > 2h and shoot to ∞ as soon as Δ = 2h. Conversely, a ce_k value that is too large will result in convergence issues, while a ce_k value that is too low will have excessive oscillations between the contact bodies. A value that was found to be a good compromise between physical accuracy and convergence was ce_k = 50.

The next hyperparameter that must be specified is the contact stiffness k_c. Not to be confused with the approximation stiffness k in Eqs. (13) and (14), the contact stiffness k_c is a scalar value that is used to scale the contact force and Jacobian. This value is adaptively readjusted every time-step to ensure that excessive hovering or penetration is minimized, and so, only the initial value must be specified, which can be found empirically. This initial k_c should be reasonably close to the value that would result in md = 2h. In our experiments, we employed a simple algorithm that decreased or increased k_c by a fraction of a percent depending on whether md was < 2h − ε₁ or > 2h + ε₂, where ε₁ and ε₂ are limits indicating an acceptable contact range. This algorithm updates k_c only when md is deviating from the region defined by [2h − ε₁, 2h + ε₂] and is otherwise left constant to prevent overshooting.

As mentioned earlier, this algorithm applies a hybrid approach in which the Jacobian of the contact and friction forces are only computed once the number of iterations passes a limit ω. This is done because often convergence can be quickly obtained even without the contact Jacobian, and with the absence of the contact Jacobian, the overall Jacobian matrix remains a banded matrix, which can be solved significantly faster. Although the Jacobian results in a decrease in iteration count, the consequent increase in computational time outweighs this benefit as IMC is able to reliably converge rapidly without it for a majority of time-steps. With this in mind, the Jacobian is crucial for completing volatile contact states with high velocities and impacts that could otherwise end the simulation prematurely. Therefore, this hybrid approach maximizes computational speed while ensuring that the simulation can consistently reach the next time-step during especially difficult contact scenarios such as inversion and the initialization phase where the rod rapidly reverts to its lowest energy state. The limit ω should be chosen empirically, so that the Jacobian is only generated when necessary. In our experiments, we used an ω of 20.

Finally, although not a hyperparameter, a damping coefficient α is used to reduce the step size of the Newton solver as the number of iterations increase for a particular time-step. For this, a simple decaying algorithm is used which reduces α by a factor of 2 every other iteration. Overall, aside from the collision detection algorithm (which is only performed on the first iteration of every time-step), the time complexity of IMC with and without Jacobian generation is $O(n)$ and $O(n2)$⁠, respectively, where n is the number of collisions detected. The full contact algorithm and its implementation in DER can be seen in Algorithms 1 and 2, respectively. In Algorithm 2, the term “free” are the indices that correspond to the free degrees-of-freedom of the elastic rod. The remaining degrees-of-freedom are “fixed” and depends on the user defined boundary conditions.

In this section, we first validate the correctness of our contact model against theory. Afterward, to observe the benefits of using IMC, we compare simulation results with the contact model (SPT) proposed by Spillmann and Teschner [27] for knots of unknotting numbers n ∈ [1, 4], which are shown in Fig. 1. For both methods, the contact model is used to simulate knot tying until a mutual termination state, which is shown in the rightmost column of Fig. 1. Finally, we compare the computational efficiency and convergence properties between the two methods. All simulations for IMC used a contact energy stiffness ce_k = 50 and a collision limit δ = 0.15, which was obtained using the method mentioned in Fig. 4.

Coulomb’s friction law states that friction is independent of velocity. To show that our model abides by this, we plot the pull forces F in Newton (N) when tightening a knot with unknotting number n = 2 with friction coefficient μ_k = 0.10 for pull speeds u of 3, 6, and 9 mm/s (pulled from both ends). A 1 m long rod with cross-sectional radius h = 1.6 mm, density ρ = 1180 kg/m³, and Young’s Modulus E = 1.8e5 Pa is used and is discretized into 301 nodes. We plot a regularized pull force Fh²/EI against $h/R$⁠, where $EI=πEh4/4$ is the flexural modulus and R is the radius of the knot loop. The radius R can be computed using the knot circumference from Fig. 1 as R = e/(2π). As Fig. 5(a) shows, the magnitude of the n = 2 pull forces is approximately the same for all pull speeds. Dynamic friction force in Coulomb’s model is independent of velocity, and therefore, this observation supports the physical correctness of the generated friction forces.

Finally, we show that the pull forces monotonically increase with μ_k in Fig. 5(b). Here, we consider the n = 4 case. When recording the pull forces for μ_k of 0.1, 0.2, and 0.3, we can see that the rate of increase is held constant, while the magnitudes monotonically increase. Furthermore, we can see that as friction increases, the inversion point occurs sooner, which indicates that the point of inversion is highly dependent on μ_k.

Simulations of knot pulling are performed for both IMC and SPT using a pull speed u = 6 mm/s, time-step dt = 0.5 ms, friction coefficient μ_k = 0.10, and the same rod properties from Sec. 3.2. With these settings, knot tying was simulated for each method with the experienced pull forces being recorded at each time-step. Figure 6(a) plots these pull forces with respect to the end-to-end shortening. Here, we see that for both IMC and SPT, as the end-to-end shortening decreases (knot is pulled tight), the pull forces increase identically as expected. They also increase in magnitude as the unknotting number n goes up, and for n = 3 and n = 4, inversion occurs, which is indicated by the sudden drop in pull force. This is shown visually in Fig. 6(b) for n = 4.

The largest difference between the methods can be seen in the considerable amount of force jittering by SPT, which occurs due to the time-step being too large. This leads to visually unrealistic results, where the knot continuously “trembles” while being tied. Conversely, IMC produces much smoother pull forces, which directly translate to visually smooth simulations for the same time-step size. One caveat that should be noted is that SPT ensures exact nonpenetration during contact, whereas IMC allows small penetrations and hovering to occur, which is physically unrealistic. Still, through the adaptive contact stiffness and the inclusion of damping, IMC contains any hovering to stay within 20 μm above the contact surface, while penetrations rarely exceed 5 μm (for comparison, cross-sectional radius is h = 1.6 mm). Thus, this minor error is largely indiscernible both visually and physically.

Overall, the enforcement of nonpenetration at every time-step limits the maximal time-step that SPT can take without experiencing significant force jittering, whereas IMC produces smooth results in exchange for contact varying within a small region.

Next, we discuss runtime comparisons and convergence characteristics. Here, we show that in addition to IMC producing smoother results than SPT at sufficiently large time-steps, IMC is also computationally competitive. Both models use the same DER implementation, which is written in c++. One discrepancy is that while SPT is directly implemented into DER in c++, IMC is written entirely in python for ease of prototyping. To minimize any performance differences arising from using different languages, we employed an LLVM-based python JIT compiler [36] for certain computational intensive portions of the code such as the chain ruling procedure from Eq. (17). The computational time for each iteration for both contact methods is recorded using the ctime library. This timing was done so that any computational time arising from IO usage recording the data and rendering the rod graphically were excluded. One thing to note is that the timings for IMC include all of the shared memory overhead between the c++ and python programs. Therefore, a significant performance increase for IMC can be expected when fully implemented in c++ and compiled without this overhead. Finally, both methods used identical Newton tolerances for all simulations, and all simulations were run on a single thread on an Intel Core i7-9700K 3.60GHz CPU.

In Table 1, the runtime, total number of iterations, average iterations per time-step (AIPT), and average time per iteration (ATPI) are reported. For all knots, we can immediately see that IMC converges with a noticeably smaller amount of iterations than SPT. While the ATPI between both methods are about equal, the simulations for IMC finish approximately 4 × faster than SPT for all knots. For both methods, AIPT increases as the knot complexity n increases as expected.

One observation was that the number of iterations to complete a time-step increased for IMC as the knot loop became extremely small and/or tightly inverted. Dynamically reducing the time-step solves this issue, but this was left out so that all reported results were for a constant time-step. For the same time-step size, SPT can more reliably converge when the knot loop becomes very small albeit with the large force jittering still present. Therefore, it may be worth investigating performance differences between the two methods when SPT uses a time-step size that is small enough to compete with the smoothness of IMC, while IMC starts at a larger time-step and dynamically reduces as convergence becomes an issue.

Still, even for a constant time-step size of 0.5 ms, IMC is seen to be more computationally efficient when pulling the knots close to taut and for the majority of the knot tying procedure, takes far less iterations to converge. This difference should only increase with IMC being implemented directly into DER in c++.

In this article, we introduced a novel contact model for DER simulations in which the contact forces are handled implicitly using a smooth penalty force formulation. This model was shown to be able to model dynamic friction and simulate knot tying accurately. We showcased comparisons with previous methods [27] and concluded that our method was capable of producing more visually smooth realistic results, while also producing physically accurate data. Furthermore, our method was stable, computationally competitive, and took minimal iterations to converge.

Although this method has shown promising results for knot tying simulations, it is not without its drawbacks. First, IMC is largely unsuitable for contact scenarios that frequently involve sudden excessively large velocities and impacts as these will result in excessive penetration and possible overshoot of the contacting body without appropriate damping. Where IMC shines is scenarios with constant sliding contact such as knot tying where contact forces may or may not gradually rise.

In addition, to be absolutely physically realistic, contact forces should equal zero when Δ > 2h, whereas the usage of Eq. (7) produces a force F_c ≠ 0 when Δ > 2h. It has been shown that employing a smooth approximation such as this can greatly improve convergence in penalty methods [25], which is one of the goals for this algorithm. As the contact forces approach zero as ce_k becomes sufficiently large while Δ > 2h, this is not a significant problem as discussed in Sec. 3.3. Still, the fact that F_c is not exactly zero is something to consider. Finally, the large amount of hyperparameters leaves something to be desired as tuning may be necessary when switching between rod properties, which can be time consuming.

Some possible future research directions involve modifications that further improve the realism of the contact model. One of these pertains to the stiffness parameter k_c. A simplification that was employed is the usage of a global stiffness parameter. More realistic contact can be simulated using local stiffness parameters as shown previously in Ref. [25] in exchange for more computation. This addition may also fix the problem where friction forces overtake the pull force of the knot if the pull speed is too low as mentioned in Sec. 2.4.

Another research direction is to extend IMC to noncircular cross sections. Currently, IMC is only capable of handling circular cross sections due to computing the minimum distance between edges using purely the centerlines. Extending IMC to noncircular cross sections while maintaining differentiability is nontrivial but will be necessary to simulate more complexly shaped rods. Furthermore, handling more complex knots will require IMC to take into consideration cross section deformation as IMC currently assumes that the cross-sectional shape remains constant.

Finally, another challenging problem that remains is the proper modeling of static friction. Although dynamic friction is adequately modeled, ultimately, subtle frictional threshold events such as the transition from sticking to sliding and vice-versa are necessary to simulate realistic contact outside the realm of constant sliding.

The source code for IMC as well as the initial knot configurations for all conducted simulation tests can be found at https://github.com/QuantuMope/imc-der.

• National Science Foundation (Grant No. IIS-1925360).

There are no conflicts of interest.

The data and information that support the findings of this article are freely available at https://github.com/QuantuMope/imc-der.