Abstract

When an overhand knot tied in an elastic rod is tightened, it can undergo a sudden change in shape through snap buckling. In this article, we use a combination of discrete differential geometry (DDG)-based simulations and tabletop experiments to explore the onset of buckling as a function of knot topology, rod geometry, and friction. In our setup, two open ends of an overhand knot are slowly pulled apart, which leads to snap buckling in the knot loop. We call this phenomenon “inversion” since the loop appears to move dramatically from one side of the knot to the other. This inversion occurs due to the coupling of elastic energy between the braid (the portion of the knot in self-contact) and the loop (the portion of the knot with two ends connected to the braid). A numerical framework is implemented that combines discrete elastic rods with a constraint-based method for frictional contact to explore inversion in overhand knots. The numerical simulation robustly captures inversion in the knot and is found to be in good agreement with experimental results. In order to gain physical insight into the inversion process, we also develop a simplified model of the knot that does not require simulation of self-contact, which allows us to visualize the bifurcation that results in snap buckling.

1 Introduction

Found all throughout the world, knots are complex geometrical patterns of slender elastic rods in self-contact. They exist for various purposes in daily life, ranging from practical (shoelaces) to decorative (Chinese knotting). Given their importance as well as their topological and mechanical complexity, knots remain a valuable yet challenging research topic. One particular application of importance lies in the medical sciences as knots have been shown to almost always be present in the chains of sufficiently long polymers [1]. For example, DNA—a particularly long polymer—is known to have various knots arise in its twists and coils. Some specific enzymes, such as topoisomerases, even help to unwind such entanglements [2]. The implications of these knots in DNA molecules are numerous: knots can block DNA replication (promoting replicon loss), increase antibiotic sensitivity, and even block gene transcription [3]. In addition, Marenduzzo et al. [4] stated that self-contact in DNA plays a vital role in changing the ejection speed of a virus from the capsid. Aside from DNA, many proteins also contain knots; as many as 750 knotted configurations are included in the protein data base [5]. Such knotted proteins have even been shown to be present in the photoreceptors of plants [6], and their complex topologies are identified to catalyze various enzymatic reactions [7]. Going further, knots also possess an essential effect on the mechanical properties of polymers, as shown by knotted DNA having more muscular rigidity than actin [8]. With all these effects of knots being discovered, researchers in the field of chemistry have started focusing on developing methods for the synthesis of molecular knots [9]. Even when limiting our scope to the medical sciences, the implications of knots in structures are endless.

In order to fully understand and utilize knots in these structures, many questions concerning the mathematical and mechanical properties of knots must be addressed. In this article, we analyze the mechanics of overhand knots that are being tightened by force applied at the knot’s ends. For some knots, this tightening can lead to a snap buckling instability causing a sudden change in the knot’s shape. An important contributing factor in this mechanical phenomenon is self-contact, which occurs along the knotted rod. From the standpoint of mathematics, a knot is a closed loop that is infinitely thin. However, in reality, the thickness of a knot is finite, and the two ends are often open. Mechanically, researchers are more interested in the strength, equilibrium configuration, and dynamic behavior of a knotted filament [10]. Hence, self-contact should be carefully considered when handling problems of knots in mechanics.

Self-contact is a problematic field in continuum mechanics due to the fact that an elastic structure’s contact zones and configurations are unknown in advance. Furthermore, there exist multiple models for the contact between elastic objects. In this article, we focus on the Coulomb frictional contact model for self-contact in knots since Coulomb frictional contact is an adequate approximation for contact dynamics at the macroscopic scale. Various approaches have been developed in the prior work to account for Coulomb frictional contact. These approaches can be categorized into three main methods: impulse-based methods, penalty energy-based methods, and constraint-based methods (nonsmooth mechanics). Here, we discuss the pros and cons of these methods and describe the method used in this article for studying the mechanics of knots.

Impulse-based methods take the displacements at contact zones into account to compute the required impulse for preventing these zones from penetrating. This method is widely used in simulations for engineering applications for its simplicity, for example, contact in filaments [11,12] and granular media simulations [13]. However, such methods tend to suffer from unrealistic “jittering” when large time-steps are used, as the contacted forces are usually computed explicitly [14]. In addition, using displacements at contact zones to compute friction cannot return a rigorous formulation of Coulomb friction.

Penalty energy-based methods are also prevalent. In the work of Li et al. [15], Patil et al. [16], Choi et al. [14], and Tong et al. [17], a penalty energy term is formulated in which the gradient of the penalty energy is used as the contact force promoting nonpenetration. Such methods are attractive due to their low computational cost since contact forces are computed based on the simulated objects’ own degrees-of-freedom (DOFs) and are visually realistic. However, except for the work in Ref. [16], few of those listed prior works compared the simulation results with experiments. Furthermore, the hyper-parameters of the penalty energy term typically require careful tuning to achieve robust numerical results.

The last class of methods are the constraint-based methods. Jean and Moreau [18,19] proposed unilateral constraints to solve dry friction in granular media. Alart and Curnier [20] developed the first approach to solve constraint-based contact dynamics with Newton’s method to find the root of a nonsmooth complementary function. Building upon the work of Alart and Curnier [20], Bertails-Descoubes et al. [21] implemented this approach to simulate fibers; however, convergence issues for high contact ratios were noted. Daviet et al. [22] combined an analytical solver with the complementary function for contacts to capture Coulomb friction in hair dynamics robustly. The algorithm in Davie et al. [22] is implemented in Ref. [23] to capture the frictional contact in assemblies of discrete elastic rods (DERs) [11,24]. Based on the previous nonsmooth optimization method, Daviet [25] proposed a general framework for simulating contact in thin nodal objects. Since constraint-based methods have clear mathematical definitions of contact and its accuracy has been supported in the prior work, for example, sound generated by objects in contact [26], we choose to use the constraint-based method to simulate the contact that occurs in overhand knots.

Snap buckling in a tightened overhand knot is just one of the countless instabilities that can be observed in slender filaments. More broadly, thin elastic structures can showcase intriguing mechanical properties when appropriate boundary conditions are applied. For example, buckling instabilities, which are an essential consideration in the design of stable, flexible structures, and intelligent systems, are highly nonlinear in geometry. Snapping and bifurcation, a process of rapid, sudden change between different stable configurations caused by external stimuli, can be seen extensively in nature and engineering. For example, Forterre et al. [27] discussed the snapping of Venus flytraps; Kebadze et al. [28] studied bistable configurations of slap bracelets; Pandey et al. [29] analyzed snap-through of a toy popper, Tong et al. [30] studied the buckling in helical rods with robotics, and Chen et al. [31] explored bistability in soft robots.

So far, many pioneering works in snap buckling have studied thin elastic structures without contact. Gomez et al. [32] and Sano et al. [33] discussed the asymmetric constraints of two ends for one-dimensional systems. Starostin and van der Heijden [34] studied the multistability in inextensible helical ribbons induced by stretching. Morigaki et al. [35] explored the snap of a thin looped paper induced by stretching. Sano and Wada [33] explored the snap buckling of a rod-like structure under the twisting of two ends. Yu and Hanna [36] investigated bifurcations of a ribbon under shearing. Finally, Zhang et al. [37] and Wan et al. [38] characterized the bistability of ribbon-like systems under loading out of the plane. In addition to these works, the bulking instability of two-dimensional curved surfaces has also been studied in prior works, e.g., snapping transitions in cylinders [3942] and snap buckling in spheres and hemispheres [39,4348]. Self-contact, which is usually expressed as an inequality constraint, also plays an important role in the buckling instabilities of slender structures, although the inclusion of self-contact typically makes buckling analysis significantly more challenging. Reference [49] studied the stability of supercoiled plasmids with self-contact to mimic DNA configurations. In the work of Ref. [50], the buckling analysis of twisted DNA was implemented to study the bounds of topoisomerase relaxation.

Our work, which focuses on snap buckling in overhand knots with self-contact, builds upon the previous work that focused on characterizing the equilibrium properties of knots in thin elastic rods. The iconic overhand knot is one of the most fundamental knots, ubiquitous in nature, and daily life. Previous works studying overhand knots have considered the effects of tightening forces and equilibrium shapes. For example, in the study by Audoly et al. [51], the relationship between traction forces and shortening of a trefoil knot (31) and cinquefoil knot (51) was given. Przybyl and Pieranski [52] analyzed the curvature and torsion of a tightening trefoil knot. Jawed et al. [53] formulated an analytical expression for the relationship between tightening forces, friction, and shortening of overhand knots with different values of crossing number. Moulton et al. [10] found the stable configuration of an open trefoil knot without self-contact. Most closely related to our work is the phenomenon shown in Fig. 1 of Ref. [11] and the subsequent work by Clauvelin et al. [54]. In Fig. 1 of Ref. [11], a trefoil knot is shown to have both continuous and discontinuous shape changes when rotating the ends. Building upon this observation, Clauvelin et al. [54] derived an asymptotic solution for trefoil and cinquefoil knots tied on an infinitely long rod, consisting of straight tails and a circular loop. They then showed that a combination of tension and twist can cause the asymptotic knot solution to buckle via either helical buckling of the knot’s tails or out-of-plane buckling of the circular loop.

Motivated by this phenomenon, we use a combination of simulations and experiments to show that overhand knots can experience a sudden change in shape simply due to tightening of the knot, and without any external twist applied at the two ends, as shown in Fig. 1. This change in the knot’s shape is referred to as “inversion” going forward in this article. Since the boundary conditions of the tightening process are straightforward, we argue that this process is topology induced, and we explore how and why inversion occurs during tightening. Our contributions are as follows: (1) we implement the constraint-based method from Refs. [20,21] with a well-known DDG model (DERs) from Refs. [11,24] in a simple way to simulate overhand knots; (2) we quantify the mechanical responses when tightening an overhand knot through simulations and experiments and discuss the contributing factors of inversion in detail, and (3) we develop a simplified knot model that quantitatively captures the inversion process, but does not require the simulation of self-contact, which we then use to gain further insight into the bifurcation underlying the inversion process.

Fig. 1
Snapshots of overhand knots in the tightening process from both experiments and simulations: (a1) n = 2, lc/L = 0.61; (a2) n = 2, lc/L = 0.75; (a3) n = 2, lc/L = 0.79; (b1) n = 3, lc/L = 0.68; (b2) n = 3, lc/L = 0.78; (b3) n = 3, lc/L = 0.79; (c1) n = 4, lc/L = 0.59; (c2) n = 4, lc/L = 0.66; (c2) n = 4, lc/L = 0.68
Fig. 1
Snapshots of overhand knots in the tightening process from both experiments and simulations: (a1) n = 2, lc/L = 0.61; (a2) n = 2, lc/L = 0.75; (a3) n = 2, lc/L = 0.79; (b1) n = 3, lc/L = 0.68; (b2) n = 3, lc/L = 0.78; (b3) n = 3, lc/L = 0.79; (c1) n = 4, lc/L = 0.59; (c2) n = 4, lc/L = 0.66; (c2) n = 4, lc/L = 0.68
Close modal

The remainder of this article is organized as follows. In Sec. 2, both the problem of inversion and the geometry of an overhand knot are formulated. In Sec. 3, the numerical framework for the DDG model and frictional responses are discussed. Next, results from experiments and simulations are compared in Sec. 4. The simplified model, which is based on the geometry of overhand knots, and its use in analyzing inversion are described in Sec. 5. Finally, conclusions and avenues for future work are laid out in Sec. 6.

2 Problem Statement

In Fig. 2, the topology of an overhand knot is shown. These knots are made up of a braid with arc length lb, a loop with arc length ll, and two tails where boundary conditions are applied. The topology of the braid can be characterized by the unknotting number n, which is the number of times one end must pass the loop to fully untie the knot. The cross-sectional radius of the rod forming the knot is h. The two ends of the overhand knot are clamped, and the clamped-to-clamped length is lc. When pulling two ends of an overhand knot, there is another intrinsic quantity not shown in Fig. 2, which is the friction coefficient μ.

Fig. 2
Configurations of overhand knots with different unknotting numbers
Fig. 2
Configurations of overhand knots with different unknotting numbers
Close modal

As shown in Fig. 1, when pulling two ends of an overhand knot with n > 1, the loop will suddenly invert and contact the braid region at a specific point. This is a snap buckling process caused by the geometry of the overhand knot. In this problem, the overhand knot is formed with a naturally straight elastic rod, and the boundary conditions consist of the two clamped ends, which are pulled apart at a steady velocity. Note that when applying the boundary conditions, the two manipulated ends of the rod have the same material frames, and the two ends have zero relative rotation. When the two manipulated ends have nonzero relative rotation, more complicated behavior of the knot can arise. For example, Bergou et al. [11] showed that applying relative rotations on these two material frames results in a change of topology of a trefoil knot. In this article, we focus on the case when the relative rotation is zero to show the influence of the knot’s topology on the nonlinear behavior of the structure rather than the influence of the boundary rotations. Indeed, inversion is a geometry-dependent phenomenon, and the influence of Young’s Modulus is ignorable, which we validate in the Supplemental Material on the ASME Digital Collection. The rods used in the experiments are made up of rubber, and Young’s Modulus is measured by first quantifying the gravito-bending length Lgb = (E h2/8ρg)1/3, where ρ is density and g is gravity [55]. When a rod is in contact with a rigid substrate, its shape is governed by the gravito-bending length Lgb. Young’s modulus E can be computed by measuring this shape. Young’s modulus E was set to 1.8 MPa, and Poisson’s ratio ν was set to 0.5 (incompressible) in the experiments and simulations, resulting in a shear modulus G = E/(2(1 + ν)) = 0.6 MPa. The inversion point is experimentally measured in terms of the clamped length lc, and the influence of the contributing factors n, h, and μ on the inversion point is fully analyzed in Sec. 4.

3 Numerical Framework for Knots

In this section, we discuss how to incorporate a constraint-based method motivated by the work of Refs. [1821] with DER [11,24] to simulate a knot with frictional contact. Constraint-based methods for simulating frictional contact responses usually require complex and careful numerical treatment. In this section, we present a simple two-step approach to combine DER and the constraint-based method in an explicit way. In the first step, we compute the friction contact forces explicitly, and in the second step, we introduce the computed friction forces into DER to update the status of the simulated rod. In this section, we first show the equations of motion (EOM) based on DER; then, we discuss how to compute the friction contact forces within the EOM. This two-step approach is able to accurately capture the nonlinear behaviors of the knot.

3.1 Discrete Elastic Rods Method.

DER is an algorithm developed by [11] to capture the nonlinear mechanical properties of a Kirchhoff elastic rod. Here, we discuss how to construct the EOM of a rod with frictional contact responses. As shown in Fig. 3(a), the centerline of the elastic rod is discretized into K nodes [x0, …, xK−1], resulting in K − 1 edges [e0, …, eK−2]. In this discrete model, each edge ei has an orthogonal reference frame [di1,di2,ti] and a material frame [mi1,mi2,ti]. The shared director ti in these two frames is the edge tangent between successive nodes so that ti = ei/‖ei‖. The reference frame [di1,di2,ti] is predefined at initial time, and a time marching scheme matches the orthogonal directors at time t to their status at time t + Δt, where Δt is the time-step size. In Fig. 3(b), the material frame at the ith edge is obtained from the reference frame at the ith edge with rotation θi. In other words, the twist of each edge can be obtained through the reference and material frames. This suggests that the rod centerline can be represented with a total of 4K − 1 degrees-of-freedom (3K for the nodal positions and K − 1 for the twist angles of each edge), which can be expressed as the following vector:
q=[x0,θ0,x1,θ1,x2,θ2,,xK2,θK2,xK1]T
(1)
where T is the transpose operator.
Fig. 3
(a) Discrete schematic diagram of an elastic rod. (b) Closeup of two connected edges showcasing its reference frame, material frame, turning angle, and twisting angle in the discrete model. (c) Discrete schematic of two contacted rod segments.
Fig. 3
(a) Discrete schematic diagram of an elastic rod. (b) Closeup of two connected edges showcasing its reference frame, material frame, turning angle, and twisting angle in the discrete model. (c) Discrete schematic of two contacted rod segments.
Close modal
With the discretization scheme from Eq. (1) defined, the formulation of elastic strains, energies and the time marching scheme will be discussed in the remaining part of this section. The strains of a deformed Kirchhoff rod are made up of stretching, bending, and twisting strains. The stretching strain of an edge is simply
εi=eie^i1
(2)
Hereafter, quantities with a ()^ indicate the undeformed status, e.g., e^i is the undeformed length of the ith edge. Bending strain can be captured with the help of a curvature binormal vector, which represents the misalignment between two consecutive edges:
(κb)i=2ei1×eiei1ei+ei1ei
(3)
The norm (κb)i=2tan(ϕi/2) is the value of the curvature, and ϕi is the turning angle between two consecutive edges, which can be observed in Fig. 3(b). The material curvatures are then given by the inner products between the curvature binormal vector and the material frame directors [mi1,mi2]:
κi1=12(mi12+mi2)(κb)iκi2=12(mi11+mi1)(κb)i
(4)
Finally, the twisting strain of an edge is
τi=θiθi1+Δτiref
(5)
where Δτiref is the discrete integrated reference twist, which can be obtained through the parallel transport between two neighboring reference frames.
Given the elastic strains, the expressions for the discrete elastic energies of an elastic rod can be formulated as follows:
Es=12i=0K2EA(εi)2e^i
(6a)
Eb=12i=1K21ei[EI1(κi1κ^i2)2+EI2(κi2κ^i2)2]
(6b)
Et=12i=1K21eiGJ(τiτ^i)2
(6c)
where A is the rod’s cross-sectional area; E is Young’s modulus; G is the shear modulus; J is the geometrical factor of the torsional rigidity, which is equal to the second polar moment of the area along the tangent ti when the cross section is circular; and I1 and I2 are the second moments of the area along the material frames mi1 and mi2, respectively. By using the potential elastic energies in Eq. (6), a first-order implicit Euler method is used to integrate the equations of motion numerically from told to tnew = told + Δt. The final equations of motion are given as follows:
MΔt(qnewqoldΔtq˙old)FnewintFnewext=0
(7a)
Fnewint=(Es+Eb+Et)qnew
(7b)
Fnewext=FG+Ffr
(7c)

In Eq. (7), Fnewint is the internal force; Fnewext is the external force, which constitute of the gravity FG and frictional contact responses Ffr; M is the diagonal mass matrix; and the subscript new refers to the state at time tnew. Note that the formulation of Fnewint is given and FG is constant during simulation. We next outline the formulation of the frictional contact responses Ffr so that Newton’s method can be implemented to solve the EOM stated in Eq. (7).

3.2 Frictional Contact.

In our numerical framework, friction is modeled with Coulomb’s friction. We have found that treating frictional contact responses as linear complementary constraints can describe Coulomb’s friction accurately. Since the linear complementary constraints describing frictional contact were first formulated by Jean and Moreau [18], many modern numerical frameworks have implemented these constraints to simulate frictional responses [2023]. In this article, we implemented the formulation from Refs. [20,21] to compute the frictional contact between two contacting edges, as shown in Fig. 3(c).

In Fig. 3(c), when two edges are in contact, linear complementary constraints exist between contact forces Ffr and the contact points. Each segment is treated as a cylinder, and contact points Ca and Cb can be computed with the algorithm presented by Ref. [56]. The relative tangent velocity (Ca relative to Cb) is ut and the contact force is Ffr=Ftfr+Fnfr. Note that n = (CaCb)/d is the normal direction, and t=ut/ut is the tangential direction. Based on the linear complementary constraints from Refs. [18,20], the relationship between contact distance d and the normal component of the frictional contact force is given by Signorini’s condition:
Fnfrn00d2h
(8)
When two segments contact (d − 2h = 0), the linear complementary condition between the friction force Ftfr and tangential relative velocity ut is expressed as follows:
μFnfr+Ftfrt00ut
(9)
The relative velocity u at the contact point, the minimal distance d, the contact normal direction n, and tangential direction t can be acquired from the DER formulation. Therefore, we can solve for the frictional contact responses Ffr with the linear complementary condition. In the work by Alart and Cunier [20], an augmented Lagrangian formulation is proposed so that the linear complementary constraints in Eqs. (8) and (9) can be rewritten as follows:
Fac(Ffr)=0=[FfrnPR+(Fnfrnρd)(Ffrt)tPD(μ(Ffrnρd)+)((Ffrt)tρut)]
(10)
Equation (10) is the linearized formulation of Eqs. (8) and (9) proposed by Alart and Cunier [20]. In this formulation, P(·) is the projection operator and ρ is the penalty coefficient defined by users. (Details can be found in the Supplemental Material.) Equation (10) can be solved with Newton’s method. However, Bertails-Descoubes et al. [21] found that when the contact ratio is high, solving Eq. (10) globally will cause convergence issues. Daviet et al. [22,25] later developed a relatively complicated Gauss–Sidel solver that treats different contact pairs locally to improve convergence performance. In this article, we use an explicit method to solve the linear complementary conditions in Eqs. (8) and (9) since we study the behavior of the knot under quasi-static conditions. We used the configuration of the knot at the previous time-step (which is known) to approximate the required quantities, e.g., contact distance d, contact normal n, and tangent direction t in the linear complementary constraints so that we can find the frictional contact responses with Eq. (10) in a robust way. The full algorithm of solving for the frictional contact responses and the validation of the frictional contact framework is presented in the Supplemental Material.

Once the frictional contact response Ffr are computed for all contacting segments with the explicit method, we have all the items in Eq. (7). The DOFs vector q for the knot can be updated via time marching with Newton’s method.

4 Results From Simulation and Experiments

In this section, we present the results from the numerical framework stated in Sec. 3 and compare them with experimental data. We discretize the elastic rod into K = 301 nodes in our simulations. A validation of our frictional model is provided in the Supplemental Material. In our rod model, 14 DOFs [x0, θ0, x1, xK−2, θK−2, xK−1] are constrained; recall that each node, e.g., x0, corresponds to three DOFs. These 14 DOFs correspond to the Cartesian coordinates and twist at the two ends. We applied this clamped boundary at the two ends for both our simulations and experiments. All other nodes are free to deform. In our simulation, we assume that (1) the effects of gravity are negligible and (2) the rod material is homogeneous. However, if desired, the effect of gravity can be easily accounted for in the discrete model as an additional external force. In the Supplemental Material, we show that gravity has minimal influence on the inversion point through simulation.

4.1 Experimental Setup and Illustration of Inversion Point.

We now introduce the setup for our experiments as well as provide a formal description of the inversion point. In our experiments, two fixtures are used to clamp and pull the two ends of a suspended overhand knot. These fixtures move collinearly in opposite directions, as shown in Fig. 4. Recall that the end-to-end length, also referred to as the clamped length, is lc. With this, we define two new quantities, H and W, to measure the geometry of the closed loop of the overhand knot. As shown in Fig. 4, H is the height of the knot as defined by the distance along the vertical symmetric division of the knot and W is the width of the knot as defined by the distance between the two ends of the braid region. As the knot tightens, lc increases while H and W smoothly decrease. However, for knots with unknotting numbers n > 1, there will be a specific point during the tightening process, where H suddenly drops to zero, which we define as the inversion point. A formal definition of the inversion point is shown in Fig. 5.

Fig. 4
Experimental setup used to observe the inversion point of overhand knots
Fig. 4
Experimental setup used to observe the inversion point of overhand knots
Close modal
Fig. 5
Illustration of the inversion point in plots of (a) l¯k versus F¯ and (b) l¯k versus H¯. A graphic of the rod state (c) before inversion and (d) post-inversion are also displayed.
Fig. 5
Illustration of the inversion point in plots of (a) l¯k versus F¯ and (b) l¯k versus H¯. A graphic of the rod state (c) before inversion and (d) post-inversion are also displayed.
Close modal
Mathematically, the geometry of the knot is most relevant with the closed loop. The size of the closed loop in Fig. 4 can be approximated by Llc, where lc is the end-to-end length and L is the total length of the rod used to tie the knot. In our experiments, we use lk = Llc to measure the inversion point. In the simulations, the inversion point is measured using two quantities: traction force and geometrical changes. We convert all quantities to be dimensionless and use ()¯ to denote this. We nondimensionalize lk, F, and H as follows:
l¯k=lkh=LlchF¯=Fh2EIH¯=HW
(11)
where EI = (π/4)E h4 (EI = EI1 = EI2) is the bending stiffness; h is the cross-sectional radius of the rod; E is Young’s modulus; and F is the traction force. In our simulations, Young’s modulus of the rod is set to E = 1.8 MPa, which is the same as the measured Young’s modulus of the rods used in experiments. Poisson’s ratio is set to 0.5, which means the material is assumed to be incompressible.

As shown in Fig. 5(a), it can be seen that the traction force F¯ increases initially as the knot is tightened and then proceeds to rapidly decrease and then rebound once l¯k reaches the inversion point. Similarly, in Fig. 5(b), H¯ drops to zero at the inversion point. The inversion points obtained from both simulation and experiments are the same. Therefore, in this article, we will use H¯ to compare the inversion points from the simulations and experiments.

Since the problem is geometry based, the contributing factors are the ones that have a direct effect on the geometry of the closed loop. With this in mind, there are three quantities that play a key role in inversion: rod radius h, unknotting number n, and friction coefficient μ. Going forward, we will see how these quantities impact the inversion point in the system. For both the simulations and experiments shown in Fig. 5, the parameters h = 1.6 mm, L = 1 m, n = 3, and μ = 0.1 were used.

4.2 Inversion Points and Contributing Factors.

Through the numerical framework and desktop experiments, we explore the influence of the three contributing factors, h, n, and μ, systematically. For simulations, we performed a parameter sweep for each parameter (n, h, and μ) where varying values were used for each parameter, while the other two remained constant. Due to the restrictions of the rods’ material in the real world, we perform experiments covering only parts of the full parameter sweep done in simulation. Regardless of this, all experimental results still adequately show the influence of the contributing factors (n, h, and μ) on the inversion point.

4.2.1 Effect of Rod Radius.

First, we quantitatively study the influence of h on inversion. When studying each contributing factor, we set the remaining contributing factors to be constant in order to most accurately discover the direct influence of the parameter in question. For each contributing factor, a total of four plots are created: simulation results of the normalized traction forces F¯ with respect to the normalized knot close-loop length l¯k, simulation results of the normalized height H¯ with respect to l¯k, comparison of simulation and experimental results, and a phase diagram.

When studying the effect of h, the unknotting number n was set to 3 and μ was set to 0.1. Furthermore, the pull speed at both ends was set to Δu = 1 mm/s in order to ensure quasi-static responses from the system. Figure 6(a) shows the evolution of F¯ with respect to l¯k for different rod radii values, h ∈ [1.6, 2.0, 2.4, 2.8, 3.2] mm. In this plot, the sudden jump in force magnitude appears at similar values of l¯k as h increases. In other word, the inversion point always happen at the similar nondimensionalized length of the closed knot loop no matter how the rod radius changes. Indeed, the closed knot loop length is nondimensionalized with h, and therefore, inversion happens more easily with larger rod radii. Likewise, Fig. 6(b) plots H¯ with respect to l¯k and the inversion points can also be located from the sudden drops in H¯. For the experimental comparison, two rods with different rod radii were used: one with h = 1.6 mm and another with h = 3.4 mm. As can be seen in Fig. 6(c), we find excellent agreement between both the simulation and experiments for both rod radii. Finally, in Fig. 7, the phase diagram of the overhand knot with n = 3 and μ = 0.1 is given. Overall, when the nondimensionalized length of the closed knot loop is the same, the instability behaviors are approximately the same for different rod radii.

Fig. 6
Effect of normalized rod radius h on inversion point: (a) the relationship between normalized traction force F¯ and normalized knot closed-loop length l¯k for different h, (b) the relationship between normalized height H¯ and normalized knot closed-loop length l¯k for different h, and (c) comparison between simulated and experimental results of normalized height H¯ versus normalized knot closed-loop length l¯k.
Fig. 6
Effect of normalized rod radius h on inversion point: (a) the relationship between normalized traction force F¯ and normalized knot closed-loop length l¯k for different h, (b) the relationship between normalized height H¯ and normalized knot closed-loop length l¯k for different h, and (c) comparison between simulated and experimental results of normalized height H¯ versus normalized knot closed-loop length l¯k.
Close modal
Fig. 7
Phase diagram of normalized rod radius h versus inversion point with n = 3 and μ = 0.1
Fig. 7
Phase diagram of normalized rod radius h versus inversion point with n = 3 and μ = 0.1
Close modal

4.2.2 Effect of Unknotting Number.

Moving forward, we now turn to the next contributing factor, the unknotting number n. Similarly to before, we set the remaining contributing factors to constant in order to observe the influence of n on the system. The remaining contributing factors are set to h = 1.6 mm and μ = 0.1.

In Figs. 8(a) and 8(b), the relationship between F¯ and H¯ with respect to l¯k is given, respectively, for several values of n. For both plots, the jump points correspond to the inversion points. For the experiments, a rod with h = 1.6 mm and μ = 0.1 was used to make overhand knots of n = 3 and n = 4. As shown in Fig. 8(c), we once again find excellent agreement between the simulated and experimental results. In Fig. 9, a phase diagram showcasing the influence of n is shown. Overall, it has been shown that the critical value of l¯k tends to increase as the unknotting number n increases (corresponding to a decreasing clamped length l¯c). We also note that, to observe the inversion phenomena, n should always be larger than 1. In other words, when n = 1, there will be no inversion in the system.

Fig. 8
Effect of unknotting number n on inversion point: (a) the relationship between normalized traction force F¯ and normalized knot closed-loop length l¯k with different n, (b) the relationship between normalized height of closed-loop H¯ and knot closed-loop length l¯k with different n, and (c) comparison between simulated and experimental results of normalized height H¯ versus normalized knot closed-loop length l¯k
Fig. 8
Effect of unknotting number n on inversion point: (a) the relationship between normalized traction force F¯ and normalized knot closed-loop length l¯k with different n, (b) the relationship between normalized height of closed-loop H¯ and knot closed-loop length l¯k with different n, and (c) comparison between simulated and experimental results of normalized height H¯ versus normalized knot closed-loop length l¯k
Close modal
Fig. 9
Phase diagram of unknotting number n versus inversion point with μ = 0.1 and h = 1.6 mm
Fig. 9
Phase diagram of unknotting number n versus inversion point with μ = 0.1 and h = 1.6 mm
Close modal

4.2.3 Effect of Friction Coefficient.

Finally, we look at the effect of the friction coefficient μ on the system. Once again, we make the remaining contributing factors constant: n is set to 3 and h is set to 3.4 mm. Figures 10(a) and 10(b) show the relationship between F¯ and H¯ with respect to l¯k, respectively, for several values of μ. For both plots, the jump points correspond to the inversion points. Here, we can see that as μ increases, the normalized clamped length at which inversion occurs decreases. For the experiments, a rod with normalized rod radius h = 3.4 mm and unknotting number n = 3 was used. In addition, chalk and glycerin were used to change the surface of the rod for the purpose of manually adjusting the frictional coefficient μ. We found that when we added chalk to the surface of rod, μ became 0.7, while glycerin produced μ = 0.1. By using this setup, we once again find excellent agreement between simulation results and experiments as shown in Fig. 10(c). Finally, in Fig. 11, the phase diagram is shown. Overall, we conclude that the threshold for inversion decreases as the friction coefficient increases.

Fig. 10
Effect of friction coefficient μ on inversion point: (a) the relationship between normalized traction force and normalized knot closed-loop length l¯k for different μ, (b) the relationship between normalized height of closed-loop H¯ and normalized knot closed-loop length l¯k for different μ, and (c) comparison between simulated and experimental results of normalized height H¯ versus normalized knot closed-loop length l¯k
Fig. 10
Effect of friction coefficient μ on inversion point: (a) the relationship between normalized traction force and normalized knot closed-loop length l¯k for different μ, (b) the relationship between normalized height of closed-loop H¯ and normalized knot closed-loop length l¯k for different μ, and (c) comparison between simulated and experimental results of normalized height H¯ versus normalized knot closed-loop length l¯k
Close modal
Fig. 11
Phase diagram of friction coefficient μ versus inversion point with n = 3 and h = 3.4 mm
Fig. 11
Phase diagram of friction coefficient μ versus inversion point with n = 3 and h = 3.4 mm
Close modal

Through desktop experiments and numerical simulations, we have discovered that h, μ, and n all play positive roles in the inversion process of overhand knots. So far, we have clarified the effects of h, μ, and n on the inversion points of overhand knots. However, exactly how these contributing factors influence the system and why inversion happens remains to be explained. In Sec. 5, we develop a simplified model of the knot, which does not require the simulation of self-contact, to explain why inversion happens and what role the contributing factors play in the system.

5 A Simplified Knot Model

In this section, we introduce a simplified model for the overhand knot with the goal of investigating why inversion occurs during the knot tightening process. The primary goal of this simplified model is to eliminate the need to simulate frictional self-contact in the numerical simulations of the knot. However, to capture the geometrically nonlinear behavior of the rod, the simplified model still depends on numerical simulation of the contact-free knot loop. As shown in Fig. 2, an overhand knot can be divided into three parts: the loop, the braid, and the tails. In the simulations and experiments described in the previous section, it was observed that the braided region loosens, and portions of the braid are fed into the loop region when inversion occurs. Furthermore, the length of the tail regions remain constant before and after inversion. We therefore conclude that inversion occurs due to a coupling of elastic energy within the braid and loop regions. In the simplified model, we ignore the tails and instead focus solely on the closed loop, which consists of the braided region and the loop region.

A schematic of the simplified knot model is shown in Fig. 12. Within the simplified model, rods with helical centerlines are used to approximate the braided region, friction within the braided region is neglected, and self-contact within the loop region is neglected. As shown in Fig. 12(a), and consistent with the formulation in Sec. 2, ll is the loop’s arc length and lb is arc length of the braided region. We also use t1 and t2 to denote the tangent directions at the ends of the loop. As the knot in Fig. 12(a) is tightened, it will deform into the knot shown in Fig. 12(b). Throughout this process, we assume that the knot’s shape remains antisymmetric about the YZ plane. We can therefore parameterize the rod’s boundary conditions using two angles, θ and ϕ, as shown in Figs. 12(c) and 12(d). The angle θ is the rotation of the tangent vectors t1 and t2 about the z-axis. Due to continuity in the closed loop, the corresponding end of the braid will also rotate about the z-axis, and this rotation angle will also be θ. As shown in Fig. 12(d), the angle ϕ is the rotation of the tangent vectors t1 and t2 about the x-axis.

Fig. 12
Simplified model of the overhand knot
Fig. 12
Simplified model of the overhand knot
Close modal
Since the tails of the rod are neglected in the simplified model, there are six geometric parameters to consider: lb, ll, h, n, θ, and ϕ. We note that for a sufficiently small radius h, the axial length of the rod within the braided portion is approximately 2 lb, and we have lklb + ll. From the three length-dependent parameters lb, ll, and h, we can obtain two nondimensional parameters, which we choose to be l¯k=lk/h and ϵ=lb/ll. We note that the nondimensional lengths l¯b=lb/h and l¯l=ll/h can be recovered from l¯k and ϵ via
l¯b=ϵl¯k1+ϵ,l¯l=l¯k1+ϵ
As mentioned previously, we have observed that inversion occurs due to coupling of elastic energy within the braid and loop regions. Therefore, with Eb denoting the elastic energy of the braid and El denoting the elastic energy of the loop, we define the following nondimensional energies:
E¯b=EbhEI,E¯l=ElhEI
(12)
To compute the nondimensional braid energy E¯b, we first parameterize the helical rods that comprise the braid region as follows:
x¯=cos(t)y¯=sin(t)z¯=b¯tb¯=l¯b2nπ2θt[nπ+θ,nπθ]
(13)
Neglecting the twisting energy of the helical rods within the braided region, we can use Eq. (13) to obtain the nondimensional bending energy of the braided region E¯b in terms of the nondimensional curvature κ¯:
κ¯=11+b¯2E¯b=2κ¯2(nπθ)1+b¯2
(14)
Next, we describe the energy of the loop E¯l. Since the loop is a suspended elastic rod in three-dimensional space, determining its configuration and elastic energy is nontrivial. However, we can use DER-based numerical simulations described in earlier sections to find the loop’s shape. Doing so, we can express E¯l as the function:
E¯l=f(ϕ,θ,l¯k,ε)
(15)
where the value of f(·) is found with DER.

We can now express the total knot energy E¯(ϕ,θ,l¯k,ε) of the closed loop as E¯b+E¯l. Equilibrium configurations of the simplified knot must satisfy E¯=0. However, rather than coupling the analytical expression for the braid energy E¯b with the numerically computed loop energy E¯l, we will instead focus on the influence of a single variable on the equilibrium configurations and their stability. Specifically, we will investigate how the rotation angle θ affects knot configurations. From the numerical simulations described in previous sections, it was observed that inversion typically occurs at the approximate parameter values ϕ = π/3 and ε=0.7. Therefore, these nondimensional parameter values will be used in the remainder of this section.

To observe the effects of θ on equilibrium configurations for the parameters ϕ = π/3 and ε=0.7, we used DER-based numerical simulations to solve for the total energy E¯ with a rod radius of h = 1.6 mm, unknotting number of n = 3, a clamped length corresponding to l¯k=137.5, and with fixed rotation angles in the range θ ∈ [0, 5] radians. Equilibrium configurations of the knot correspond to rotation angles at which E¯/θ=0. Figure 13(a) shows the value of E¯ for θ ∈ [0, 5]. Here, we see that there are three critical points, with the outer two critical points being stable and the intermediate critical point being unstable. We emphasize that the configurations corresponding to these critical points are equilibria with respect to variations in θ only, and the other parameters (specifically, ϕ and ε) are held fixed. Figure 13(b) shows the value of E¯/θ for the range of θ, where we again see that there are three critical points.

Fig. 13
Effects of rotation angle θ and normalized knot closed-loop length l¯k on the energy E¯ in the simplified knot model with ϕ = π/3, ε=0.7, and n = 3: (a) the total energy E¯ has three critical points with respect to variations in θ when l¯k=137.5, (b) the gradient ∂E¯/∂θ is zero at the three critical points when l¯k=137.5, and (c) the values of θ at the critical points undergo twofold bifurcations as the normalized closed-loop length l¯k decreases, with the second fold corresponding to the snap-through instability observed in previous sections
Fig. 13
Effects of rotation angle θ and normalized knot closed-loop length l¯k on the energy E¯ in the simplified knot model with ϕ = π/3, ε=0.7, and n = 3: (a) the total energy E¯ has three critical points with respect to variations in θ when l¯k=137.5, (b) the gradient ∂E¯/∂θ is zero at the three critical points when l¯k=137.5, and (c) the values of θ at the critical points undergo twofold bifurcations as the normalized closed-loop length l¯k decreases, with the second fold corresponding to the snap-through instability observed in previous sections
Close modal

We can now observe how the number of critical points changes as the knot closed-loop length l¯k changes. Figure 13(c) shows the values of θ at the critical points as the clamped length varies. For sufficiently large l¯k, there is a single stable critical point. As l¯k decreases, two additional critical points, one stable and one unstable, appear in a fold bifurcation. Further decreasing l¯k causes the unstable critical point to collide with the original stable critical point in a second fold bifurcation, leaving a single stable critical point for sufficiently small l¯k. The second fold bifurcation corresponds to the snap-through instability observed in previous sections. We note that although the simplified model exhibits a hysteresis loop between the twofold bifurcations, this hysteretic behavior is not observed in experiments since increasing l¯k causes the knot to sag due to gravity, which is not included in our model.

This simplified knot model can now be used to predict how changes in the rod’s radius and the unknotting number affect the critical clamped length at which inversion occurs. This in turn allows us to compare our simplified knot model’s output with the simulated and experimental results from previous sections. For values of the rod radius h in the range [1.5, 3.5] mm and an unknotting number of n = 3, we computed the critical closed-loop length l¯k at which the snap-through instability occurs. These results are shown by the dashed line in Fig. 14, along with the experimental results and the simulation results based on the full DER model that includes contact and friction, with a friction coefficient of μ = 0.1. A similar comparison is shown in Fig. 15, where the radius is fixed at h = 1.6 mm and the unknotting number is varied. We again emphasize that ϵ and ϕ are held fixed in the simplified model for the data shown in Figs. 14 and 15.

Fig. 14
Comparison between simplified model, the full DER simulation, and the experimental measurements for different values of h with n = 3 and μ = 0.1
Fig. 14
Comparison between simplified model, the full DER simulation, and the experimental measurements for different values of h with n = 3 and μ = 0.1
Close modal
Fig. 15
Comparison between simplified model, the full DER simulation, and the experimental measurements for different values of n with h = 1.6 mm and μ = 0.1
Fig. 15
Comparison between simplified model, the full DER simulation, and the experimental measurements for different values of n with h = 1.6 mm and μ = 0.1
Close modal

From Figs. 14 and 15, we can observe both qualitative and quantitative agreement between the simplified model, the full DER simulation, and the experiments. Particularly, we see that l¯k is invariant under changes in the rod radius h¯ and grows with the unknotting number n, corresponding to a decreasing clamped length l¯c.

6 Conclusion

In this article, we studied the snap-through buckling of overhand knots during the pulling of the knot’s clamped ends. For this purpose, a discrete differential geometry-based model, one-dimensional discrete elastic rods, with a constraint-based method for frictional contact was introduced in order to study the effect of overhand knot topology on physical behavior when tightened. Through both experiments and simulations, we found that complex topologies of knots can induce intriguing phenomena such as snap buckling and that parameters such as μ, h, and n have an influence on the geometry of overhand knots and thus, an influence on the snap buckling point itself. As μ, h, and n increase, inversion, which is the process of snap buckling, occurs earlier with an increasingly loose overhand knot. To make it clear that inversion is a result of the energy coupling in the system, we analyzed the topology of overhand knots and formulated a simplified model to predict the inversion point. Furthermore, we have observed that as the closed loop of the knot tightens, some parts of the braid (the self-contact zone) loosen suddenly, which directly causes inversion, and beyond this point, the loop will stay in contact with the braided region.

Although a few previous studies on snap buckling in elastic structures have considered the influence of self-contact, contact-dependent buckling is still a challenging topic both numerically and analytically. Therefore, we believe our work to be promising. The inversion in overhand knots shows that topologies caused by self-contact can introduce interesting behaviors of elastic structures. Because of the energy coupling between the contact zones and regions without contact, inversion occurs, resulting in overhand knots that self-tangle in complex ways. Such an understanding of inversion may prove helpful in the explanations of tangles for polymers like DNA and proteins or even be a motivation in the design of soft robotics.

Funding Data

  • NSF (Grant No. IIS-1925360) for D.T., A.C., J.J, and M.K.J.

  • NSF (Grant Nos. CAREER-2047663, CMMI-2101751, and CMMI-2053971) for M.K.J.

Conflict of Interest

There are no conflicts of interest.

Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

References

1.
Sumners
,
D.
, and
Whittington
,
S.
,
1988
, “
Knots in Self-Avoiding Walks
,”
J. Phys. A: Math. Gen.
,
21
(
7
), p.
1689
.
2.
Sumners
,
D. W.
,
2011
, “
DNA Knots: Theory and Experiments
,”
Prog. Theor. Phys. Suppl.
,
191
, pp.
1
19
.
3.
Mann
,
J. K.
,
2007
,
DNA Knotting: Occurrences, Consequences & Resolution
,
The Florida State University
,
Tallahassee, FL
.
4.
Marenduzzo
,
D.
,
Micheletti
,
C.
,
Orlandini
,
E.
, and
Sumners
,
D. W.
,
2013
, “
Topological Friction Strongly Affects Viral DNA Ejection
,”
Proc. Natl. Acad. Sci. USA
,
110
(
50
), pp.
20081
20086
.
5.
Ziegler
,
F.
,
Lim
,
N. C.
,
Mandal
,
S. S.
,
Pelz
,
B.
,
Ng
,
W.-P.
,
Schlierf
,
M.
,
Jackson
,
S. E.
, and
Rief
,
M.
,
2016
, “
Knotting and Unknotting of a Protein in Single Molecule Experiments
,”
Proc. Natl. Acad. Sci. USA
,
113
(
27
), pp.
7533
7538
.
6.
Wagner
,
J. R.
,
Brunzelle
,
J. S.
,
Forest
,
K. T.
, and
Vierstra
,
R. D.
,
2005
, “
A Light-Sensing Knot Revealed by the Structure of the Chromophore-Binding Domain of Phytochrome
,”
Nature
,
438
(
7066
), pp.
325
331
.
7.
Lim
,
N. C.
, and
Jackson
,
S. E.
,
2015
, “
Molecular Knots in Biology and Chemistry
,”
J. Phys.: Condens. Matter.
,
27
(
35
), p.
354101
.
8.
Arai
,
Y.
,
Yasuda
,
R.
,
Akashi
,
K.-I.
,
Harada
,
Y.
,
Miyata
,
H.
,
Kinosita
,
K.
, and
Itoh
,
H.
,
1999
, “
Tying a Molecular Knot With Optical Tweezers
,”
Nature
,
399
(
6735
), pp.
446
448
.
9.
Leigh
,
D. A.
,
Pirvu
,
L.
, and
Schaufelberger
,
F.
,
2019
, “
Stereoselective Synthesis of Molecular Square and Granny Knots
,”
J. Am. Chem. Soc.
,
141
(
14
), pp.
6054
6059
.
10.
Moulton
,
D. E.
,
Grandgeorge
,
P.
, and
Neukirch
,
S.
,
2018
, “
Stable Elastic Knots With No Self-Contact
,”
J. Mech. Phys. Solids
,
116
, pp.
33
53
.
11.
Bergou
,
M.
,
Wardetzky
,
M.
,
Robinson
,
S.
,
Audoly
,
B.
, and
Grinspun
,
E.
,
2008
, “
Discrete Elastic Rods
,”
ACM Trans. Graphics (TOG)
, pp.
1
12
.
12.
Spillmann
,
J.
, and
Teschner
,
M.
,
2008
, “
An Adaptive Contact Model for the Robust Simulation of Knots
,”
Comput. Graphics Forum
,
27
(
2
), pp.
497
506
. Wiley Online Library.
13.
Lee
,
S. J.
, and
Hashash
,
Y. M.
,
2015
, “
idem: An Impulse-Based Discrete Element Method for Fast Granular Dynamics
,”
Int. J. Numer. Methods Eng.
,
104
(
2
), pp.
79
103
.
14.
Choi
,
A.
,
Tong
,
D.
,
Jawed
,
M. K.
, and
Joo
,
J.
,
2021
, “
Implicit Contact Model for Discrete Elastic Rods in Knot Tying
,”
ASME J. Appl. Mech.
,
88
(
5
), p.
051010
.
15.
Li
,
M.
,
Ferguson
,
Z.
,
Schneider
,
T.
,
Langlois
,
T.
,
Zorin
,
D.
,
Panozzo
,
D.
,
Jiang
,
C.
, and
Kaufman
,
D. M.
,
2020
, “
Incremental Potential Contact: Intersection- and Inversion-Free, Large-Deformation Dynamics
,”
ACM Trans. Graphics (TOG)
,
39
(
4
), pp.
49:1
49:20
.
16.
Patil
,
V. P.
,
Sandt
,
J. D.
,
Kolle
,
M.
, and
Dunkel
,
J.
,
2020
, “
Topological Mechanics of Knots and Tangles
,”
Science
,
367
(
6473
), pp.
71
75
.
17.
Tong
,
D.
,
Choi
,
A.
,
Joo
,
J.
, and
Jawed
,
M. K.
,
2023
, “
A Fully Implicit Method for Robust Frictional Contact Handling in Elastic Rods
,”
Extreme Mech. Lett.
,
58
, p.
101924
.
18.
Jean
,
M.
, and
Moreau
,
J. J.
,
1987
, “Dynamics in the Presence of Unilateral Contacts and Dry Friction: A Numerical Approach,”
Unilateral Problems in Structural Analysis—2
,
G.
Del Piero
and
F.
Maceri
, eds.,
Springer
, Berlin, pp.
151
196
.
19.
Jean
,
M.
, and
Moreau
,
J. J.
,
1992
, “
Unilaterality and Dry Friction in the Dynamics of Rigid Body Collections
,”
1st Contact Mechanics International Symposium
,
Lausanne, Switzerland
.
20.
Alart
,
P.
, and
Curnier
,
A.
,
1991
, “
A Mixed Formulation for Frictional Contact Problems Prone to Newton Like Solution Methods
,”
Comput. Methods Appl. Mech. Eng.
,
92
(
3
), pp.
353
375
.
21.
Bertails-Descoubes
,
F.
,
Cadoux
,
F.
,
Daviet
,
G.
, and
Acary
,
V.
,
2011
, “
A Nonsmooth Newton Solver for Capturing Exact Coulomb Friction in Fiber Assemblies
,”
ACM Trans. Graphics (TOG)
,
30
(
1
), pp.
1
14
.
22.
Daviet
,
G.
,
Bertails-Descoubes
,
F.
, and
Boissieux
,
L.
,
2011
, “
A Hybrid Iterative Solver for Robustly Capturing Coulomb Friction in Hair Dynamics
,”
ACM Trans. Graphics (TOG)
,
30
(
6
), pp.
1
12
.
23.
Kaufman
,
D. M.
,
Tamstorf
,
R.
,
Smith
,
B.
,
Aubry
,
J.-M.
, and
Grinspun
,
E.
,
2014
, “
Adaptive Nonlinearity for Collisions in Complex Rod Assemblies
,”
ACM Trans. Graphics (TOG)
,
33
(
4
), pp.
1
12
.
24.
Bergou
,
M.
,
Audoly
,
B.
,
Vouga
,
E.
,
Wardetzky
,
M.
, and
Grinspun
,
E.
,
2010
, “
Discrete Viscous Threads
,”
ACM Trans. Graphics (TOG)
,
29
(
4
), pp.
1
10
.
25.
Daviet
,
G.
,
2020
, “
Simple and Scalable Frictional Contacts for Thin Nodal Objects
,”
ACM Trans. Graphics (TOG)
,
39
(
4
), p.
61–1
.
26.
Schweickart
,
E.
,
James
,
D. L.
, and
Marschner
,
S.
,
2017
, “
Animating Elastic Rods With Sound
,”
ACM Trans. Graphics (TOG)
,
36
(
4
), pp.
1
10
.
27.
Forterre
,
Y.
,
Skotheim
,
J. M.
,
Dumais
,
J.
, and
Mahadevan
,
L.
,
2005
, “
How the Venus Flytrap Snaps
,”
Nature
,
433
(
7024
), pp.
421
425
.
28.
Kebadze
,
E.
,
Guest
,
S.
, and
Pellegrino
,
S.
,
2004
, “
Bistable Prestressed Shell Structures
,”
Int. J. Solids Struct.
,
41
(
11–12
), pp.
2801
2820
.
29.
Pandey
,
A.
,
Moulton
,
D. E.
,
Vella
,
D.
, and
Holmes
,
D. P.
,
2014
, “
Dynamics of Snapping Beams and Jumping Poppers
,”
EPL (Europhys. Lett.)
,
105
(
2
), p.
24001
.
30.
Tong
,
D.
,
Borum
,
A.
, and
Jawed
,
M. K.
,
2021
, “
Automated Stability Testing of Elastic Rods With Helical Centerlines Using a Robotic System
,”
IEEE Rob. Autom. Lett.
,
7
(
2
), pp.
1126
1133
.
31.
Chen
,
T.
,
Bilal
,
O. R.
,
Shea
,
K.
, and
Daraio
,
C.
,
2018
, “
Harnessing Bistability for Directional Propulsion of Soft, Untethered Robots
,”
Proc. Natl. Acad. Sci.
,
115
(
22
), pp.
5698
5702
.
32.
Gomez
,
M.
,
Moulton
,
D. E.
, and
Vella
,
D.
,
2017
, “
Critical Slowing Down in Purely Elastic ‘Snap-Through Instabilities
,”
Nat. Phys.
,
13
(
2
), pp.
142
145
.
33.
Sano
,
T. G.
, and
Wada
,
H.
,
2019
, “
Twist-Induced Snapping in a Bent Elastic Rod and Ribbon
,”
Phys. Rev. Lett.
,
122
(
11
), p.
114301
.
34.
Starostin
,
E.
, and
van der Heijden
,
G.
,
2008
, “
Tension-Induced Multistability in Inextensible Helical Ribbons
,”
Phys. Rev. Lett.
,
101
(
8
), p.
084301
.
35.
Morigaki
,
Y.
,
Wada
,
H.
, and
Tanaka
,
Y.
,
2016
, “
Stretching an Elastic Loop: Crease, Helicoid, and Pop Out
,”
Phys. Rev. Lett.
,
117
(
19
), p.
198003
.
36.
Yu
,
T.
, and
Hanna
,
J.
,
2019
, “
Bifurcations of Buckled, Clamped Anisotropic Rods and Thin Bands Under Lateral End Translations
,”
J. Mech. Phys. Solids
,
122
, pp.
657
685
.
37.
Zhang
,
Y.
,
Jiao
,
Y.
,
Wu
,
J.
,
Ma
,
Y.
, and
Feng
,
X.
,
2020
, “
Configurations Evolution of a Buckled Ribbon in Response to Out-of-Plane Loading
,”
Extreme Mech. Lett.
,
34
, p.
100604
.
38.
Wan
,
G.
,
Liu
,
Y.
,
Xu
,
Z.
,
Jin
,
C.
,
Dong
,
L.
,
Han
,
X.
,
Zhang
,
J. X.
, and
Chen
,
Z.
,
2020
, “
Tunable Bistability of a Clamped Elastic Beam
,”
Extreme Mech. Lett.
,
34
, p.
100603
.
39.
Bende
,
N. P.
,
Evans
,
A. A.
,
Innes-Gold
,
S.
,
Marin
,
L. A.
,
Cohen
,
I.
,
Hayward
,
R. C.
, and
Santangelo
,
C. D.
,
2015
, “
Geometrically Controlled Snapping Transitions in Shells With Curved Creases
,”
Proc. Natl. Acad. Sci. USA
,
112
(
36
), pp.
11175
11180
.
40.
Pezzulla
,
M.
,
Stoop
,
N.
,
Jiang
,
X.
, and
Holmes
,
D. P.
,
2017
, “
Curvature-Driven Morphing of Non-Euclidean Shells
,”
Proc. R. Soc. A: Math. Phys. Eng. Sci.
,
473
(
2201
), p.
20170087
.
41.
Jiang
,
X.
,
Pezzulla
,
M.
,
Shao
,
H.
,
Ghosh
,
T. K.
, and
Holmes
,
D. P.
,
2018
, “
Snapping of Bistable, Prestressed Cylindrical Shells
,”
EPL (Europhys. Lett.)
,
122
(
6
), p.
64003
.
42.
Lavrenčič
,
M.
, and
Brank
,
B.
,
2018
, “
Simulation of Shell Buckling by Implicit Dynamics and Numerically Dissipative Schemes
,”
Thin-Walled Struct.
,
132
, pp.
682
699
.
43.
Huang
,
N.-C.
,
1964
, “
Unsymmetrical Buckling of Thin Shallow Spherical Shells
,”
ASME J. Appl. Mech.
,
31
(
3
), pp.
447
457
.
44.
Lazarus
,
A.
,
Florijn
,
H.
, and
Reis
,
P. M.
,
2012
, “
Geometry-Induced Rigidity in Nonspherical Pressurized Elastic Shells
,”
Phys. Rev. Lett.
,
109
(
14
), p.
144301
.
45.
Shim
,
J.
,
Perdigou
,
C.
,
Chen
,
E. R.
,
Bertoldi
,
K.
, and
Reis
,
P. M.
,
2012
, “
Buckling-Induced Encapsulation of Structured Elastic Shells Under Pressure
,”
Proc. Natl. Acad. Sci. USA
,
109
(
16
), pp.
5978
5983
.
46.
Marthelot
,
J.
,
López Jiménez
,
F.
,
Lee
,
A.
,
Hutchinson
,
J. W.
, and
Reis
,
P. M.
,
2017
, “
Buckling of a Pressurized Hemispherical Shell Subjected to a Probing Force
,”
ASME J. Appl. Mech.
,
84
(
12
), p.
121005
.
47.
Evkin
,
A.
,
Kolesnikov
,
M.
, and
Prikazchikov
,
D. A.
,
2017
, “
Buckling of a Spherical Shell Under External Pressure and Inward Concentrated Load: Asymptotic Solution
,”
Math. Mech. Solids
,
22
(
6
), pp.
1425
1437
.
48.
Hutchinson
,
J. W.
, and
Thompson
,
J. M. T.
,
2018
, “
Imperfections and Energy Barriers in Shell Buckling
,”
Int. J. Solids Struct.
,
148
, pp.
157
168
.
49.
Coleman
,
B. D.
,
Swigon
,
D.
, and
Tobias
,
I.
,
2000
, “
Elastic Stability of Dna Configurations. II. Supercoiled Plasmids With Self-Contact
,”
Phys. Rev. E
,
61
(
1
), p.
759
.
50.
Thompson
,
J. M. T.
,
2008
, “
Single-Molecule Magnetic Tweezer Tests on DNA: Bounds on Topoisomerase Relaxation
,”
Proc. R. Soc. A: Math. Phys. Eng. Sci.
,
464
(
2099
), pp.
2811
2829
.
51.
Audoly
,
B.
,
Clauvelin
,
N.
, and
Neukirch
,
S.
,
2007
, “
Elastic Knots
,”
Phys. Rev. Lett.
,
99
(
16
), p.
164301
.
52.
Przybyl
,
S.
, and
Pieranski
,
P.
,
2009
, “
Tightening of the Elastic Overhand Knot
,”
Phys. Rev. E
,
79
(
3
), p.
031801
.
53.
Jawed
,
M. K.
,
Dieleman
,
P.
,
Audoly
,
B.
, and
Reis
,
P. M.
,
2015
, “
Untangling the Mechanics and Topology in the Frictional Response of Long Overhand Elastic Knots
,”
Phys. Rev. Lett.
,
115
(
11
), p.
118302
.
54.
Clauvelin
,
N.
,
Audoly
,
B.
, and
Neukirch
,
S.
,
2009
, “
Matched Asymptotic Expansions for Twisted Elastic Knots: A Self-Contact Problem With Non-Trivial Contact Topology
,”
J. Mech. Phys. Solids
,
57
(
9
), pp.
1623
1656
.
55.
Jawed
,
M. K.
,
Da
,
F.
,
Joo
,
J.
,
Grinspun
,
E.
, and
Reis
,
P. M.
,
2014
, “
Coiling of Elastic Rods on Rigid Substrates
,”
Proc. Natl. Acad. Sci. USA
,
111
(
41
), pp.
14663
14668
.
56.
Lumelsky
,
V. J.
,
1985
, “
On Fast Computation of Distance Between Line Segments
,”
Inf. Process. Lett.
,
21
(
2
), pp.
55
61
.

Supplementary data