## Abstract

Rod–rod contact is critical in simulating knots and tangles. To simulate contact, typically a contact force is applied to enforce nonpenetration condition. This force is often applied explicitly (Euler forward). At every time-step in a dynamic simulation, the equations of motions are solved over and over again until the right amount of contact force successfully imposes the nonpenetration condition. There are two drawbacks: (1) Explicit implementation brings numerical convergence issues. (2) Solving equations of motion iteratively to find this right contact force slows down the simulation. In this article, we propose a simple, efficient, and fully implicit contact model with high convergence properties. This model is shown to be capable of taking large time-steps without forfeiting accuracy during knot tying simulations when compared to previous methods. We introduce a new contact potential, based on a smoothed segment–segment distance function, that is an analytic function of the four endpoints of the two contacting edges. Since this contact potential is differentiable, we can incorporate its force (negative gradient of the energy) and Jacobian (negative Hessian of the energy) into the elastic rod simulation.

## 1 Introduction

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.

## 2 Methodology

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.

### 2.1 Discrete Elastic Rods.

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 *θ*^{i} 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 4*N* − 1 total degrees-of-freedom (DOF), which can be represented by the following vector **q** = [**x**_{1}, *θ*_{1}, **x**_{2}, *θ*_{2}, …, **x**_{N−1}, *θ*_{N−1}, **x**_{N}]^{T}; here, ^{T} denotes transposition operator. In this representation, the location of the *i*th node, **x**_{i}, corresponds to the 4*i* − 3, 4*i* − 2, and 4*i* − 1th entries of **q**.

*E*

^{elastic}, is composed of the stretching energy

*E*

^{s}, bending energy

*E*

^{b}, and twisting energy

*E*

^{t}as follows:

**e**

_{i}, $Eib$ is the bending energy at node

**x**

_{i}, and $Eit$ is the twisting energy at node

**x**

_{i}. For a rod with Young’s modulus

*E*, shear modulus

*G*, area moment of inertia

*I*, polar moment of inertia

*J*, and cross-sectional area

*A*, the energies are formulated as follows:

**e**

_{i}in undeformed state,

*κ*

_{i}is the curvature vector at node

**x**

_{i}, $\kappa i0$ is the undeformed curvature for the same node,

*τ*

_{i}is the integrated twist at node

**x**

_{i}, and $dLi=(|e\xafi\u22121|+|e\xafi|)/2$ is the Voronoi length in the undeformed state. Typically, the rod is uniformly discretized and the length of each edge

*dL*is the same throughout the rod. All these energies are functions of the configuration of the rod represented by

**q**[4].

*q*

_{i}(i.e.,

*i*th element of

**q**), the elastic forces (for nodal positions) and elastic moments (for twist angles) can be written as follows:

*i*th element of (4

*N*− 1) sized elastic force vector,

**F**

^{int}.

**F**

^{ext}is the external force vector (e.g., contact forces, gravity),

**M**is the diagonal mass matrix, and $q\xa8$ is the second derivative of the DOFs with respect to time. In DER, backward Euler method is used to solve the 4

*N*− 1 equations of motion to update the DOF vector

**q**. For the march from time

*t*

_{i}to

*t*

_{i+1}=

*t*

_{i}+

*dt*, where

*dt*is the time-step, Eq. (4) can be rewritten as follows:

**q**(

*t*

_{i}) are the DOFs at time-step

*t*

_{i}and $q\u02d9(ti)$ are the velocities at time-step

*t*

_{i}. In the setup studied in this article, the three external forces are (1) the contact force

**F**

_{c}, (2) the friction force

**F**

_{fr}, and (3) the viscous force

**F**

_{v}; and therefore, we have

**F**

^{ext}=

**F**

_{c}+

**F**

_{fr}+

**F**

_{v}. The formulation of these three forces will be discussed throughout the remainder of this article.

**q**(

*t*

_{i+1}) iteratively. Each element of the Jacobian matrix

**J**at row

*i*and column

*j*is expressed as follows:

**J**

_{c},

**J**

_{fr}, and

**J**

_{v}are square matrices representing the gradient of the three external forces—contact, friction, and viscous, respectively—with respect to the DOFs. Once

**q**(

*t*

_{i+1}) is known, the new velocity is simply $q\u02d9(ti+1)=(q(ti+1)\u2212q(ti))/dt$.

Normally, if the gradient of an external force, $\u2202Fiext/\u2202qj$, 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.

### 2.2 Contact Model.

Referring to Fig. 2(b), denote $xi,xi+1,xj,xj+1\u2208R3$ as the Cartesian nodal coordinates of the *i*th and *j*th 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\u2208X$. 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\u2192R1$. 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 $\u2212\u2207E(xij)$ as well as the negative Hessian $\u2212\u22072E(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 2*h* 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.

*h*and then sharply increases as Δ approaches 2

*h*. Although the gradients in the region Δ > 2

*h*are nonzero, the inclusion of these “cushioning” forces greatly aid convergence and reduce any unwanted oscillating behavior that can often occur in penalty methods. The effect of these nonzero gradients are explained in detail in Secs. 3.3 and 4. Moving on, the key component of the contact model involves obtaining a differentiable analytical expression for Δ, which is difficult as computing the minimum distance between two line segments is a highly nonlinear and noncontinuous process. Next, we briefly describe Lumelsky’s min-distance algorithm and a new modified smooth approximation that will be used as Δ(

**x**).

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 $\u2219$ is the dot product.

**e**

_{ij}=

**x**

_{i}−

**x**

_{j}to the ones already previously formulated:

**e**

_{i}and

**e**

_{j}. With these vectors, we can then calculate the necessary intermediary values as follows, where

*i*and

*j*subscripts are left out for clarity:

*F*(

*x*) as a fix bound function, where all values above 1 are 1 and all values below 0 are 0. Anything between is outputted identically. This function is the piecewise function shown below and is first of the two remaining noncontinuous components.

*t*,

*u*∈ [0, 1] are the ratios that determine at which point along the length of each edge the connecting min-distance vector lies as shown in Fig. 2(b). With this in mind, the fix bound function

*F*(

*x*) is used to ensure that these values do not go outside the appropriate range. As two edges become increasingly closer to parallel,

*λ*approaches 0 and becomes 0 when perfectly parallel. To prevent division by zero, a piecewise function is used to describe the assignments of

*t*.

*u*

_{f}≠

*u*(i.e.

*u*< 0 or

*u*> 1), the

*t*value is reassigned as follows.

### 2.3 Contact Energy, Its Gradient, and Hessian.

*E*(Δ(

**x**)), Δ(

**x**) must be differentiable. In Sec. 2.2, we introduced an algorithm that can compute Δ. Now, we modify the min-distance algorithm into a differentiable analytical expression. As Eq. (8) is analytical, the only necessary modifications lie in Eqs. (10) and (11). First, the fixbound function

*F*(

*x*) can be modeled with the following smooth approximation, which is denoted by

*H*(

*x*).

*k*is a hyperparameter, which determines how stiff the curves are. A larger

*k*value will result in a more accurate approximation but will result in “stiff” first and second derivatives leading to reduced convergence; thus, this value should be determined empirically. Next is the conditional reassignment in Eq. (11). As the reassignment only depends on whether

*u*

_{f}≠

*u*, this is equivalent to the reassignment only occurring when

*u*< 0 or

*u*> 1. To model this, we can use a boxcar function denoted

*B*(

*x*), which consists of two compounded logistic functions.

*H*(

*x*) and

*B*(

*x*) are plotted in Fig. 3 for a value of

*k*= 50, which is the value used to produce the simulation results. This

*k*value is a good tradeoff between accuracy and reliable convergence.

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*_{1} · *D*_{2} − *S*_{2} · *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.

*E*(

**x**). It should be noted that an end-to-end differentiation of

*E*(

**x**) ends up with an extremely large and complex equation when using symbolic differentiation [35]. Therefore, to greatly simplify the expression and improve computational efficiency, the gradient and Hessian for several of the intermediary algorithmic values are taken and then chain ruled together. Effectively, we can define

*E*(

**x**) by the following equivalent functional:

*f*(·) are all functions of

**x**, chain rule tells us that we can obtain the gradient of the contact energy by the following:

Here, we see that for any arbitrary edge combination **x**, the produced force $\u2212\u2207E(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 (4*N* − 1)-sized **F**_{c} vector located at the following positions: 4*i* − 3, 4*i* − 2, 4*i* − 1, 4(*i* + 1) − 3, 4(*i* + 1) − 2, 4(*i* + 1) − 1, 4*j* − 3, 4*j* − 2, 4*j* − 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 $\u2212\u22072E(x)$. Once obtained, the Hessian is added to the (4*N* − 1 × 4*N* − 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.

### 2.4 Adding Friction.

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\u2208R12$ 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)).

Frictional force is independent of velocity, and

‖

**F**_{fr}‖ =*μ*_{k}·*F*_{n}during sliding,

*μ*

_{k}is the dynamic friction coefficient and

*F*

_{n}is the normal force (more details later in this section). From the contact model, we were able to derive $\u2212\u2207E(x),$ which is equivalently the normal contact forces

**F**

_{c}. As mentioned earlier, the gradient of the contact energy will always produce forces that are along the contact normal. Therefore, we can obtain

*F*

_{n}at the

*i*th edge of the contact pair

**x**

_{ij}by simply summing up the contact forces on the

*i*th and

*i*+ 1th nodes: $Fn=\Vert Fci+Fci+1\Vert $. We can also use these contact forces to obtain the contact norm $n=(Fci+Fci+1)/Fn$. The direction of friction is then determined by the tangential relative velocity $vrelT$ of edge

*i*with respect to edge

*j*. This can be obtained using the nodal velocities

**v**

_{i},

**v**

_{i+1},

**v**

_{j}, and

**v**

_{j+1}as shown.

*i*th edge using Coulomb’s friction equation as shown below. Here, we add a weight

*γ*∈ [0, 1] to get rid of frictional forces between edges with extremely small relative velocities as this can cause unwanted behavior. To maintain differentiability, this weight is obtained using a smooth Heaviside step function with the tangential relative velocity as input and the stiffness

*k*set to 50. Here,

*c*determines the limit for the step transition and was set to 0.15 for our experiments. Note that this limit

*c*must take into consideration the scaling of the model, which is explained in Sec. 2.5.

*i*th edge $Ffrie$, we can do the same for the

*j*th edge as well to obtain $Ffrje$. The computed frictional forces are then equally distributed to each node and then concatenated four times to form the final friction vector $Ffrij\u2208R12$.

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.

#### Implicit contact method.

#### Discrete elastic rods.

### 2.5 Full Algorithm.

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:

*δ*, the collision limit,*k*_{c}, the contact stiffness,*ce*_{k}, the contact energy stiffness, and*ω*, the number of iterations before the hybrid algorithm computes the Jacobian.

*S*= 1/

*h*, so that the adjusted rod radius equals a unit value of 1. To ensure that the distance at which two edges experience a force is very close to the rod surface, the following energy function is used.

*δ*is the threshold value used to determine when two edges are “in contact.” This value is fed into a collision detection algorithm, which returns all edge combinations falling into this threshold, which is denoted by the following set:

*k*

_{c}accordingly.

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 Δ > 2*h* and shoot to ∞ as soon as Δ = 2*h*. 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* = 2*h*. 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 < 2*h* − *ε*_{1} or > 2*h* + *ε*_{2}, where *ε*_{1} and *ε*_{2} are limits indicating an acceptable contact range. This algorithm updates *k*_{c} only when *md* is deviating from the region defined by [2*h* − *ε*_{1}, 2*h* + *ε*_{2}] 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.

## 3 Results

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.

### 3.1 Damping for Stability.

**F**

_{v}on a node as denoted in Eq. (23), where

*η*(Pa · s) is a viscosity coefficient and

*dL*is the node’s Voronoi length.

*N*− 1). This damping force was also added to SPT for fair comparison. For all simulations, a viscosity value of

*η*= 0.01 Pa · s was used unless otherwise specified.

### 3.2 Theoretical Validation.

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^{3}, and Young’s Modulus *E* = 1.8*e*5 Pa is used and is discretized into 301 nodes. We plot a regularized pull force *Fh*^{2}/*EI* against $h/R$, where $EI=\pi 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.

*n*= 1) are compared with the predictive model by Audoly et al. [23]. This model states the following theoretical equivalence, where

*σ*is a numerical constant (

*σ*= 0.492 for trefoil knots) and $\u03f5=h/R$. The ± term is the frictional component from tightening (+) and loosening (−).

*η*to 0.0005 as loosening can be sensitive to small forces. As shown in Fig. 5(a), the recorded pull forces roughly follow the curves of the predictive model albeit with some displacements when tightening increases sufficiently. This can be attributed to imperfections in the predictive model as Eq. (25) is an elegant lightweight solution that does not perfectly model friction when the knot becomes sufficiently tight, i.e., $h/R$ becomes large. Still, aside from the displacement, the pull forces follow the rate of increase/decrease of the predictive model well, which is a good indicator of correctness.

*ε*

_{0}value, we obtain

*μ*

_{theory}= 0.1306, which is reasonably close to the friction coefficient used in simulation,

*μ*

_{k}= 0.10.

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}.

### 3.3 Pull Force Accuracy.

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.

### 3.4 Runtime.

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.

Model | n | AIPT | ATPI (ms) | Iters | Time (s) |
---|---|---|---|---|---|

IMC | 1 | 2.23 | 2.24 | 245,824 | 551 |

2 | 2.67 | 2.20 | 245,830 | 542 | |

3 | 3.04 | 2.18 | 261,707 | 570 | |

4 | 3.24 | 2.12 | 239,987 | 509 | |

SPT | 1 | 8.25 | 2.30 | 907,541 | 2085 |

2 | 9.28 | 2.26 | 853,935 | 1931 | |

3 | 10.25 | 2.26 | 881,907 | 1990 | |

4 | 10.94 | 2.21 | 810,160 | 1790 |

Model | n | AIPT | ATPI (ms) | Iters | Time (s) |
---|---|---|---|---|---|

IMC | 1 | 2.23 | 2.24 | 245,824 | 551 |

2 | 2.67 | 2.20 | 245,830 | 542 | |

3 | 3.04 | 2.18 | 261,707 | 570 | |

4 | 3.24 | 2.12 | 239,987 | 509 | |

SPT | 1 | 8.25 | 2.30 | 907,541 | 2085 |

2 | 9.28 | 2.26 | 853,935 | 1931 | |

3 | 10.25 | 2.26 | 881,907 | 1990 | |

4 | 10.94 | 2.21 | 810,160 | 1790 |

Note: AIPT, average iterations per time-step; ATPI, average time per iteration; Iters, the total number of Newton’s iterations that were necessary to complete the simulation. The time is the total computational time to completion.

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++.

## 4 Concluding Remarks

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 Δ > 2*h*, whereas the usage of Eq. (7) produces a force *F*_{c} ≠ 0 when Δ > 2*h*. 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 Δ > 2*h*, 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.

## 5 Source Code

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.

## Funding Data

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

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

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