Dynamic modeling of a sliding ring on an elastic rod with incremental potential formulation

Here, one provides a brief overview of the DDG-based numerical frameworks applied to both 1D rod structures [42, 43] and 2D plate objects [59, 60]. These frameworks have undergone validation through comparisons with either experimental data or analytical solutions in previous studies [50, 18, 47]. The simulation domain is discretized into numerous nodes and elements, effectively conceptualized as a mass-spring-damper system. Within this system, a lumped mass is positioned at each vertex, accompanied by its corresponding discrete energy contributions. Referring to Fig. 1(a), the basic rod element combines $3$ vertices, $\mathcal{R}_{ijk}:=\left(\mathbf{x}_{i},\mathbf{x}_{j},\mathbf{x}_{k}\right)\in\mathbb{R}^{9}$ with $2$ edges, which are defined by

	$\displaystyle\mathbf{e}_{ij}$	$\displaystyle=\mathbf{x}_{j}-\mathbf{x}_{i},$		(1a)
	$\displaystyle\mathbf{e}_{jk}$	$\displaystyle=\mathbf{x}_{k}-\mathbf{x}_{j},$		(1b)
	$\displaystyle\mathbf{t}_{ij}$	$\displaystyle=\frac{\mathbf{e}_{ij}}{|\mathbf{e}_{ij}|},$		(1c)
	$\displaystyle\mathbf{t}_{jk}$	$\displaystyle=\frac{\mathbf{e}_{jk}}{|\mathbf{e}_{jk}|},$		(1d)
	$\displaystyle\Delta l_{ijk}$	$\displaystyle=\frac{1}{2}(|\mathbf{e}_{ij}|+|\mathbf{e}_{jk}|).$		(1e)

For a 1D rod with Young’s modulus $E$, cross section area $A$, and moment of inertia $I$, the simplified discrete elastic energies are given by [18]

	$\displaystyle E^{s}_{\mathcal{R}}$	$\displaystyle=\sum_{ij}\frac{1}{2}{EA}\frac{\left({|\mathbf{e}_{ij}|}-{|\bar{\mathbf{e}}_{ij}|}\right)^{2}}{{|\bar{\mathbf{e}}_{ij}|}},$		(2a)
	$\displaystyle E^{b}_{\mathcal{R}}$	$\displaystyle=\sum_{ijk}\frac{1}{2}{EI}\frac{\left(|\mathbf{t}_{ij}-\mathbf{t}_{jk}|\right)^{2}}{{{\Delta\bar{l}_{ijk}}}},$		(2b)

where a bar on the top represents the quantity in reference configuration. In this context, we have omitted consideration of the twisting energy associated with the rod’s centerline. However, it is worth noting that if deemed necessary, this twisting energy can be readily incorporated by introducing an additional degree of freedom in the form of a twisting angle, as explained in prior works [42, 43].

The plate element is defined by $4$ vertices, $\mathcal{P}_{ijkl}:=\left(\mathbf{x}_{i},\mathbf{x}_{j},\mathbf{x}_{k},\mathbf{x}_{l}\right)\in\mathbb{R}^{12}$ with $2$ triangular faces, which are given by

	$\displaystyle\mathbf{e}_{ij}$	$\displaystyle=\mathbf{x}_{j}-\mathbf{x}_{i},$		(3a)
	$\displaystyle\mathbf{e}_{ik}$	$\displaystyle=\mathbf{x}_{k}-\mathbf{x}_{i},$		(3b)
	$\displaystyle\mathbf{e}_{lj}$	$\displaystyle=\mathbf{x}_{j}-\mathbf{x}_{l},$		(3c)
	$\displaystyle\mathbf{e}_{lk}$	$\displaystyle=\mathbf{x}_{k}-\mathbf{x}_{l},$		(3d)
	$\displaystyle\mathbf{t}_{ij}$	$\displaystyle=\frac{\mathbf{e}_{ij}}{|\mathbf{e}_{ij}|},$		(3e)
	$\displaystyle\mathbf{t}_{ik}$	$\displaystyle=\frac{\mathbf{e}_{ik}}{|\mathbf{e}_{ik}|},$		(3f)
	$\displaystyle\mathbf{t}_{lj}$	$\displaystyle=\frac{\mathbf{e}_{lj}}{|\mathbf{e}_{lj}|},$		(3g)
	$\displaystyle\mathbf{t}_{lk}$	$\displaystyle=\frac{\mathbf{e}_{lk}}{|\mathbf{e}_{lk}|},$		(3h)
	$\displaystyle\mathbf{n}_{\alpha}$	$\displaystyle=\mathbf{t}_{ik}\times\mathbf{t}_{ij},$		(3i)
	$\displaystyle\mathbf{n}_{\beta}$	$\displaystyle=\mathbf{t}_{lj}\times\mathbf{t}_{lk},$		(3j)

where $\mathbf{n}_{\alpha}$ and $\mathbf{n}_{\beta}$ are the surface normal vectors of the two triangular faces, referring to Fig. 1(b). For a 2D plate with Young’s modulus $Y$ and thickness $b$, the discrete elastic energies are the sum of in-plane stretching and out-of-plane bending [60],

	$\displaystyle E^{s}_{\mathcal{P}}$	$\displaystyle=\sum_{ij}\frac{\sqrt{3}}{4}{Yb}{\left({|\mathbf{e}_{ij}|}-{|\bar{\mathbf{e}}_{ij}|}\right)^{2}},$		(4a)
	$\displaystyle E^{b}_{\mathcal{P}}$	$\displaystyle=\sum_{\alpha\beta}\frac{1}{12\sqrt{3}}{Yb^{3}}{\left(|\mathbf{n}_{\alpha}-\mathbf{n}_{\beta}|\right)^{2}}.$		(4b)

The incremental potential contact method represents an effective numerical model for addressing nonlinear elastodynamic systems that involve contact [53]. Only particle-particle contact is considered in present method, and other arbitrary dimension (surfaces and curves) contact can be derived in a similar approach [54].

When the Euclidean distance between two approaching nodes, $\mathcal{C}_{ij}:=\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)\in\mathbb{R}^{6}$, is smaller than a threshold, the repulsive force needs to be included into the discrete dynamic system thus the non-penetration condition can be achieved. The incremental potential model uses a nonlinear barrier-type potential with $C^{2}$ continuity to achieve the intersection free condition between two approaching nodes [53],

$$E_{\mathcal{C}}=\begin{cases}-K_{\mathcal{C}}\left[(d-\hat{d})^{2}\log(\frac{d}{\hat{d}})\right]&\;\mathrm{when}\;0<d<\hat{d}\\ 0&\;\mathrm{when}\;d\geqslant\hat{d},\end{cases}$$

(5)

where $d=|\mathbf{x}_{j}-\mathbf{x}_{i}|$ is the Euclidean distance between two considered nodes, $K_{\mathcal{C}}$ is the stiffness parameter, and $\hat{d}$ is a barrier parameter. It is assumed here that the particle size is zero, and the relative distance can be easily derived by subtracting the particle radius if the ball with finite size is considered. By taking the variation of nonlinear barrier potential, the contact force along the normal direction is given by [53],

$$|\mathbf{F}_{\mathcal{C}}^{n}|=\begin{cases}K_{\mathcal{C}}\left[2(d-\hat{d})\log(\frac{d}{\hat{d}})+\frac{(d-\hat{d})^{2}}{d}\right]&\;\mathrm{when}\;0<d<\hat{d}\\ 0&\;\mathrm{when}\;d\geqslant\hat{d}.\end{cases}$$

(6)

The contact force is zero when the distance is beyond the threshold, $\hat{d}$, and would gradually increase as the distance between two nodes is smaller. The dependency of normalized contact force, $|\mathbf{F}_{\mathcal{C}}^{n}|/K_{\mathcal{C}}$, on the relative distance, $d$, for different barrier parameter, $\hat{d}\in\{0.50,0.75,1.00\}$, is plotted in Fig. 2(a).

Then, the tangential component of contact force is applied through a dissipative potential form, which follows the maximal dissipation principle [58, 53],

$$|\mathbf{F}_{\mathcal{C}}^{t}|=\begin{cases}\mu|\mathbf{F}_{\mathcal{C}}^{n}|\left(-\frac{|\mathbf{v}^{t}|^{2}}{\epsilon_{v}^{2}}+\frac{|\mathbf{v}^{t}|}{\epsilon_{v}}\right)&\;\mathrm{when}\;0\leqslant|\mathbf{v}^{t}|\leqslant\epsilon_{v},\\ \mu|\mathbf{F}_{\mathcal{C}}^{n}|&\;\mathrm{when}\;|\mathbf{v}^{t}|>\epsilon_{v},\end{cases}$$

(7)

where $\mu$ is the frictional coefficient, $\mathbf{v}^{t}$ is the relative velocity along the tangential direction between two reference nodes, and $\epsilon_{v}$ is the velocity magnitude bound, and below which the sliding velocity is treated as static. The dependency of normalized tangential frictional force, $|\mathbf{F}_{\mathcal{C}}^{t}|/\mu|\mathbf{F}_{\mathcal{C}}^{n}|$, on the relative tangential velocity, $|\mathbf{v}^{t}|$, for different velocity bound, $\epsilon_{v}\in\{0.25,0.50,0.75\}$, is plotted in Fig. 2(c). This numerical treatment can smooth the discontinuity between the static friction and the dynamic friction, such that the elastodynamic system can be solved variationally.

The conventional incremental potential formulation could handle the interaction between a ring and a rod if they are both discretized into nodes and edges. However, the redundant nodal number would slow down the overall computational efficiency. Thus, a novel barrier functional is introduced here to handle the interaction between a rigid ring and a flexible rod.

As shown in Fig. 3, an entangled ring-rod coupled system is considered, and the rod is represented by multiple nodes while the ring is only a single vertex. The contact element is constructed based on the minimum distance between the sliding ring and the deformable rod, and $3$ nodes are considered here, $\mathcal{J}_{ijr}:=\left(\mathbf{x}_{i},\mathbf{x}_{j},\mathbf{x}_{r}\right)\in\mathbb{R}^{9}$, where $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$ are the vertices on rod, and $\mathbf{x}_{r}$ represents the position of sliding ring. The minimum distance $d$ between the joint and the cable here is defined by

$$d=\mathrm{min}\{|\mathbf{x}_{r}-\mathbf{x}_{i}|,|\mathbf{x}_{r}-\mathbf{x}_{j}|,\frac{|(\mathbf{x}_{r}-\mathbf{x}_{i})\times(\mathbf{x}_{r}-\mathbf{x}_{j})|}{|\mathbf{x}_{j}-\mathbf{x}_{i}|}\}.$$

(8)

To ensure the non-deviation condition between the sliding joint and the flexible cable, the following barrier energy functional is constructed by

$$E_{\mathcal{J}}=\begin{cases}K_{\mathcal{J}}\left[d^{2}\log(1-\frac{d}{\tilde{d}})\right]&\;\mathrm{when}\;0<d<\tilde{d}\\ +\infty&\;\mathrm{when}\;d\geqslant\tilde{d}\end{cases},$$

(9)

where $K_{\mathcal{J}}$ is the stiffness parameter and $\tilde{d}$ is a barrier parameter similar to the $\hat{d}$ in Eq.(6). Correspondingly, the ring radius here is assumed to be zero, and the ring-to-rod distance can be simply modified by minus the ring radius. It should also be noted that the ring-rod coupled element $\mathcal{J}_{ijr}$ would update with time when the minimum distance between the moving ring and deformable rod is changed. Again, This newly introduced barrier function for ring-to-rod model is also with $C^{2}$ continuous, and the normal contact force is close to zero when the distance is small, while the force is almost infinite when the distance is close to the threshold, $\tilde{d}$,

$$|\mathbf{F}_{\mathcal{J}}^{n}|=\begin{cases}K_{\mathcal{J}}\left[2d\log(1-\frac{d}{\tilde{d}})+\frac{d^{2}}{d-\tilde{d}}\right]&\;\mathrm{when}\;0<d<\tilde{d}\\ +\infty&\;\mathrm{when}\;d\geqslant\tilde{d}\end{cases}.$$

(10)

With this barrier functional, the non-deviation condition between sliding ring joint and elastic cable is successfully achieved. The dependency of normalized interaction force, $|\mathbf{F}_{\mathcal{J}}^{n}|/K_{\mathcal{J}}$, on the relative distance, $d$, for different barrier parameter, $\tilde{d}\in\{0.50,0.75,1.00\}$, is plotted in Fig. 2(b). The formulation of frictional force for ring-to-rod contact is identical to Eq.(7).

This section discusses the time marching scheme for elastodynamic system with non-penetration contact and ring-to-rod interaction. An elastic structure is discretized into $N$ nodes and treated as a multibody dynamic system, with degrees of freedom $\mathbf{q}\in\mathbb{R}^{3N}$. Its second order dynamic growing equations of motion is,

$$\mathbb{M}\ddot{\mathbf{q}}(t)+\mathbb{C}\dot{\mathbf{q}}(t)+\mathbb{K}{\mathbf{q}}(t)=\mathbf{F}_{\mathrm{ext}}(t),$$

(11)

where $\mathbb{M}$ is the time-invariant $3N\times 3N$ diagonal mass matrix, $\mathbb{K}=\nabla^{2}E$ (where $E$ is the total elastic potential) is the tangential stiffness matrix, $\mathbb{C}=\alpha\mathbb{M}+\beta\mathbb{K}$ is the Rayleigh damping matrix, and $\mathbf{F}_{\mathrm{ext}}$ is the external force vector (such as the gravity). Notice that the system is nonlinear, i.e., $\mathbb{C}$ and $\mathbb{K}$ are time variant.

To achieve the non-penetration condition between two approaching nodes as well as the non-deviation condition between sliding joint and elastic rod, the incremental potential and the dissipative friction need to be included. Geometric detect is updated at each time step to compute $\mathbf{F}_{\mathcal{C}}$ (for particle-particle contact) and $\mathbf{F}_{\mathcal{J}}$ (for ring-to-rod interaction). Here, the normal contact force, $\mathbf{F}_{\mathcal{C}}^{n}$ (and $\mathbf{F}_{\mathcal{J}}^{n}$), is treated as a position-dependent elastic force, while the tangential friction force, $\mathbf{F}_{\mathcal{C}}^{t}$ (and $\mathbf{F}_{\mathcal{J}}^{t}$), is considered as a velocity-dependent damping force.

The implicit Euler method is used to numerically update the DOF vector and its velocity from time step $t_{k}$ to $t_{k+1}=t_{k}+h$ ($h$ is the time step size) [61]:

	$\displaystyle{\mathbb{M}}\ddot{\mathbf{q}}(t_{k+1})$	$\displaystyle={\mathbf{F}}_{\textrm{int}}(t_{k+1})+\mathbf{F}_{\mathcal{C}}(t_{k+1})+\mathbf{F}_{\mathcal{J}}(t_{k+1})+\mathbf{F}_{d}(t_{k+1})+\mathbf{F}_{{g}}(t_{k+1})$		(12a)
	$\displaystyle{\mathbf{q}}(t_{k+1})$	$\displaystyle={\mathbf{q}}(t_{k})+h\dot{\mathbf{q}}(t_{k+1})$		(12b)
	$\displaystyle\dot{\mathbf{q}}(t_{k+1})$	$\displaystyle=\dot{\mathbf{q}}(t_{k})+h\ddot{\mathbf{q}}(t_{k+1})$		(12c)

where ${\mathbf{F}}_{\textrm{int}}=-\nabla E$ is the internal elastic force, $\mathbf{F}_{d}=\mathbb{C}\dot{\mathbf{q}}$ is the damping force vector, and $\mathbf{F}_{{g}}$ is the external gravity force vector.

As the interaction potential is based on a $C^{2}$ continuous formulation, Newton-Raphson is adopted to solve the nonlinear equations of motion, e.g., the Jacobian matrix associated with Eq.(12) is available. It is important to note that the damping force and frictional force are treated semi-implicit, i.e., to get the viscoelastic damping force at $t_{k+1}$, the component of Rayleigh damping matrix, $\mathbb{K}$ is evaluated at time $t=t_{k}$, while the velocity, $\dot{\mathbf{q}}$, is evaluated at time $t=t_{k+1}$; to get the tangential frictional force at $t_{k+1}$, the normal contact, $\mathbf{F}_{\mathcal{C}}^{n}$ (as well as $\mathbf{F}_{\mathcal{J}}^{n}$), is derived at time $t=t_{k}$, while the relative velocity $\mathbf{v}_{t}$ is evaluated at time $t=t_{k+1}$; all other forces are treated fully implicitly and updated through a classical gradient descent algorithm.

The most basic scenario, which involves a ring sliding on an inclined rigid bar, is examined and compared to an analytical solution. In this particular case, our focus is solely on the interaction between the ring and the rod, as particle-particle contact does not come into play. The problem setup is illustrated in Fig. 4, and all relevant geometric and physical parameters are detailed in Table 1. The bar is discretized into $200$ nodes, resulting overall DOF vector (rod together with ring) $\mathbf{q}\in\mathbb{R}^{3\times 201}$. The time step size is set to be $h=1$ ms, but a larger time step size is also acceptable. The inclined angle is selected as $\theta=30^{\circ}$. The first node and the last node of the bar are constrained with Dirichlet boundary condition to avoid the rigid body motion. Due to the relatively large bending stiffness, the rod deformation is negligible and the structure is assumed to be rigid.

The point mass would slide along the tangential direction of the inclined bar if its frictional coefficient is smaller than the threshold, e.g., $\mu<\tan\theta$, and the trajectory of the ring is given by,

$$s=\frac{1}{2}gt^{2}(\sin\theta-\mu\cos\theta).$$

(13)

Fig. 4(b) shows the change of position as a function of time from both numerical data and analytical result. The agreement of lines (analytical solution) and symbols (numerical data) validates the proposed discrete numerical model.

For validation purposes, a more intricate scenario is examined, involving a single ring sliding on a rigid tautochrone curve. It’s noteworthy that an analytical solution for this case is also accessible. The tautochrone curve’s shape can be attained by subjecting a flexible rod to compression under the influence of gravity. In its deformable state, it adopts the form of a catenary [62]. Referring to Fig. 5, a circular ring entangled with a rigid curve is applied, and the function of which is given by,

	$\displaystyle x$	$\displaystyle=r(\theta-\sin\theta)$		(14a)
	$\displaystyle y$	$\displaystyle=r(\cos\theta-1)$		(14b)
	$\displaystyle\mathrm{with}\;\theta$	$\displaystyle\in[0,2\pi].$		(14c)

A tautochrone with $r=0.1$ m is adopted, resulting arc-length $L=0.8$ m. Again, similar to the previous case, the rod is cut into $200$ nodes and the first and last node of rod are fixed; The bending stiffness is also enormous to accomplish the structural rigidity.

In Fig. 5, we show the $x$ displacement of the ring evolves as a function of time from our numerical side. Intriguingly, the rings with different start heights, $y$ , reach the minimum with identical time, $T_{\mathrm{critical}}\sim 0.314$ s, as the name of Tautochrone curve. The analytical solution of drop time can be derived based on classical calculus of variations,

$$T_{\mathrm{critical}}=\pi\sqrt{\frac{r}{g}}.$$

(15)

A good agreement has been found for the critical time predicted by numerical simulation and analytical formulation.

Moving forward, the flexibility of rod and the nonlinear dynamic coupling between ring and rod are illustrated. Here, a rigid ring moving on a catenary curve is adopted, and the bending stiffness of the suspended flexible rod (which is known as catenary) is significantly reduced such that the structure changes from rigid to elastic. A circular ring entangled with a catenary curve is applied, and the function of which is given by,

$$y=A\cosh(\frac{x}{A}),$$

(16)

where its arc-length is $L=2.0$ m. With a compressive distance ratio $\Delta L/L=20\%$ (and, therefore, $x\in[-0.8,0.8]$ m), the control parameter would be $A=0.6764$ m. The geometric and physical parameters are in Table 1. Here, some material damping is added into dynamic system to avoid unrealistic vibratory motion of structure.

The rigid ring is yoked onto the left side of flexible rod and would move to the right due to the existence of gravity, which is similar to the previous scenario. On the other side, the rod would also be bent as the movement of ring, because the ring can be regarded as a point load onto the curved beam. Snapshots are provided in Fig. 6. In contrast to the previous example, the ring will not move to the same height on the right side if the rod is relatively soft, as some of the gravitational potential is switched to the elastic energy of structure. Moverover, the energy would dissipate faster and the vibratory motion of ring would gradually disappear as the enlarge of frictional coefficient, as expected. Seeing Fig. 7 for more details. The convergent study and the performance of the numerical parameters ($\tilde{d}$ and $\epsilon_{v}$) can be found in Appendix A and Appendix B, separately.

Subsequently, the numerical model is employed to tackle a more intricate scenario - the simulation of hanging clothes. This particular scenario has garnered significant attention within the computer graphics community, primarily for generating visually appealing animations [46, 63, 64]. In this simulation, the 2D surface is divided into multiple triangular meshes, while the 1D rod is discretized into nodes and edges. All the relevant physical and geometric parameters for this numerical test can be found in Table 2. The elastic rod is represented by $100$ nodes and the flexible cloth is divided into $1402$ nodes (and, therefore, $4090$ connected edges and $2688$ triangular meshes). The first and last two nodes are manually fixed to establish the clamped-clamped boundary condition; all nodes on cloth are free to move and $5$ massless rings are used to connect thin cloth and slender rod. It is worth noting that the nodal positions on plate are used directly to represent the sliding rings and no extra vertex is employed to account for coupled joint. The discrete time step size is $h=0.1$ms.

Initially, the hanging rod would deform from a straight configuration into a curved shape; next, the sliding joint would move towards to the middle as the rod is inclined and the compressive loads are applied back to cloth; Finally, some vivid wrinkles would appear due the high flexibility of elastic cloth under compressive load, seeing Fig. 7(a)-(d) for details. To show the effectiveness of our ring-to-rod model, the constraints on the right side are released at $t=1.0$ s, such that the clamped-clamped beam would switch to the clamped-free boundary condition (a cantilever beam model). The hanging cloth would then slide together with the coupled joints from middle to right due to the asymmetric constraints of elastic rod and the external gravitational force; the rings would gradually shed one by one from the slender rods and loss all connections completely after $t=1.4$ s, referring to Fig. 7(e)-(h). The total numerical simulation can be done within $1$ minute.

Lastly, our numerical model is put to use in simulating the dynamic interaction between a tethered rod and a flexible net, involving the use of joint rings. The objective here is to emulate the mechanism by which the net closes during a debris capture mission [12, 16, 15]. Referring to Fig. 9(a), a rigid object with spherical shape is on top of an elastic net, and the tether rods (coloured by green) for close mechanism are connected with a hexagonal net by multiple rigid rings (coloured by red). The net is discretized into $6931$ nodes and $8316$ edges; each tether rod (overall $6$) is divided by $64$ vertex and $63$ edges; $8$ joint rings are setup for each tether rod. Again, the ring vertex is the same one with the nodes on net and no extra degrees of freedom are required to mimic the mechanics of rings. Here, both ring-to-rod interaction and ball-to-net contact are included into the numerical simulation. The time step size is $h=1$ ms, which is identical to the previous numerical examples. All other geometric and physical parameters for this numerical experiments are listed in Table 3.

The overall capture progress can be divided into $3$ phases. (i) During $0\mathrm{s}\leqslant t<5\mathrm{s}$, the rigid ball and the flexible net would drop under gravity and contact with each other. (ii) When $5\mathrm{s}\leqslant t<25\mathrm{s}$, the six corner nodes are manually moved towards center with speed $1.0\mathrm{m}/\mathrm{s}$, until the distance between adjacent corners is smaller than $5.0$m, e.g., configurations in Fig. 9(d). (iii) After $t=25$s, the six tether rods are stretched out, and, therefore, the net can be pulled up and partially closed due to the connection through joint rings, referring to Fig. 9(e)-(i). The stretching speed of green tether rods is also set to be $1.0\mathrm{m}/\mathrm{s}$ for convenience. The open size of tether-net system after closing can be tuned during phase (ii). The total computational time for this case is about $3$ minutes.