Cone complementarity approach for dynamic analysis of multiple pendulums

Abstract

The multibody system dynamics approach allows describing equations of motion for a dynamic system in a straightforward manner. This approach can be applied to a wide variety of applications that consist of interconnected components which may be rigid or deformable. Even though there are a number of applications in multibody dynamics, the contact description within multibody dynamics still remains challenging. A user of the multibody approach may face the problem of thousands or millions of contacts between particles and bodies. The objective of this article is to demonstrate a computationally straightforward approach for a planar case with multiple contacts. To this end, this article introduces a planar approach based on the cone complementarity problem and applies it to a practical problem of granular chains.

Introduction

A multibody system consists of a number of bodies which can interact via constraints or forces. The forces can be described conditionally if, for example, there is physical contact between the bodies. Accordingly, a number of individual solid bodies in a bulk of granular material ¹ can move freely until they establish contact with other bodies or solid walls during which the contact (collision) forces (impulses) alter the response of the bodies.

Granular chains ² are shown to be proper models for polymers driven far from equilibrium. ³ In the application of granular chains, it is important to obtain an accurate multiple impact models with a numerical method. Models which are governed by Newton’s and Poisson’s restitution law are widely used in describing the contact laws. Coulomb’s unilateral contact law with dry friction can be used to model the interaction between multiple particles. ⁴ Multibody dynamics and collision dynamics can be simultaneously applied to describing the behavior of granular chains. Individual multiple pendulums have been studied in the literature solely as a physics problem^5,6 with different practical applications in the analysis of walking^7–9 or bearings absorbing earthquake shocks. ¹⁰

Regardless of the applications, the problem of interacting pendulums is theoretically challenging and fascinating. Techniques of multibody dynamics can be combined with the collision dynamics of colliding pairs to obtain the mutual interactions of pendulums in time. Obviously, the time integrals of interactions are dependent on a number of geometric and physical parameters in the system under investigation. If one wants to simulate the granular with many contacts successfully, small time steps which can achieve numerical stability may be needed. Therefore, these issues motivate researchers to investigate the innovative time integration method to deal with the multiple contacts. In this field, Pang ¹¹ has investigated in-depth time integration methods for calculating multiple frictional contacts with local compliance.

Using unilateral constraints, complementarity formulas can compute contact impulses to avoid penetration between rigid bodies. Simulations of multiple contacts can be performed with linear complementarity problem (LCP) and nonlinear complementarity problem (NCP) methods. When solving the time integration problem, LCP can unify linear and quadratic programmer solvers. ¹² The methods of Gauss–Seidel and Jacobi are widely used to solve LCP. For the numerical solution of a large-scale symmetric positive-definite LCP, Kocvara and Zowe ¹³ have proposed an algorithm. Yao et al. ¹⁴ have established a linear complementary model to describe the non-holonomic system with friction by LCP. However, while the friction is nonzero or minor, the time integration interpretation will be lost, which makes the LCP solver imprecise. ¹⁵ Therefore, a more complex NCP which can be seen as an extension of LCP was proposed. DE Stewart and Trinkle ¹⁶ have proposed an implicit time-stepping method for simulating rigid body contact with Coulomb friction and inelastic impacts. The method is based on the NCP, which calculates the generalized position at the end of each time step. In this article, a semi-implicit time-stepping method is utilized with the complementarity principle to calculate the generalized position at the end of the time step. Tasora and Anitescu ¹⁷ have implemented C++ into the NCP solver to solve multiple unilateral contacts with friction with more than 100,000 colliding rigid bodies, which shows remarkable performance compared with other algorithms. However, when dealing with a large number of contacts and polyhedral approximation in friction, LCP and NCP solvers remain limited. ¹⁸

An alternative approach for solving contact problems is the so-called penalty method (i.e. continuous method contact force model ¹⁹ or complaint contact force model), ²⁰ which can be used to analyze contact forces as a continuous function of indentation and compliance of the contact surface. ²¹ This method typically needs small time step in the time integration scheme because of the stability limitations. ¹⁶ It is also noteworthy that choosing parameters requires little effort, because the selection of parameters depends on the contact case. Moreover, the selection of parameters’ value may lead to high stiffness of the ordinary differential equation (ODE), which can make time differentiation slower due to small time step requirements. Therefore, the selection has to be done with care keeping in mind that the system does not lead to unnecessary high eigenfrequencies.

In the two-dimensional Coulomb case, this article introduces an optimization-based method of the cone complementarity problem (CCP) for the simulation of interacting multiple pendulums. This method is used for the simulation of non-smooth rigid multibody dynamics with collision, contact, and friction, solving a convex quadratic program based on a fixed time step. ²² Tasora et al. have chosen to use the interior point method with a CCP, which turned out to be more accurate than the Gauss-Jacobi. ²¹ M Anitescu and colleagues^23,24 have proposed a time integration formulation method and a fixed-point iteration algorithm to solve a large CCP with low calculation which can handle large-scale contacts and large granular flow problems. Alessandro Tasora et al. have found that the fixed-point iteration algorithm displays linear complexity when solving a large CCP. In other words, the simulation time could increase linearly when the number of bodies increases in the model. ²¹ The complementarity condition to build a Coulomb dry friction model combined with a non-penetration condition in a spatial case is introduced in Negrut et al. ¹⁸

This article provides an important opportunity to advance the understanding of multibody contact in planar cases. The differential complementarity approach is utilized to analyze the contact dynamics of the two flexible quadruple pendulums. Comparing with the three-dimensional (3D) Coulomb friction model,^18,23,24 the explicit two-dimensional (2D) expressions in this article are much simpler. The changes of kinetic energy and potential energy have been explained by comparing the introduced approach with penalty method.

This article reviews and analyzes basic methods and formulations used in contact dynamics with the CCP. The overall structure of the study is as follows. Section “Methods and formulations” begins by laying out the theoretical methods and formulations of the research, including the friction model, motion, and CCP. The application of multibody contacts is studied through the contact description of multiple pendulums in section “Dynamic simulation of multiple pendulum interactions.” The contact dynamics is separated into three parts: contact with the ground, contact with each other, and the rope inextensibility constraint with the explanation of force transformation matrices and gap functions. Section “Performance investigation” displays the numerical results of contacts two pendulum. The paper gives some conclusions of the method for solving the CCP and a short discussion of the future direction of the study in the final section.

Methods and formulations

This article uses a complementarity approach to express the Coulomb friction model. The approach employs a non-penetration condition, and it is based on a set of differential equations and algebraic inequalities, which characterize the dynamics of multibody systems with friction and contact. The resulting problem is expressed in the form of a CCP.

System state

In the planar multibody case, the position of a particle of a body with respect to the global coordinate system can be described with the help of a body reference coordinate system. Therefore, a set of generalized coordinates for a multibody system can be written as

$$\mathit{q}={\left[\begin{array}{ccccc}{\mathit{R}}_1^T& {\phi }_1& \cdots & {\mathit{R}}_{n_b}^T& {\phi }_{n_b}\end{array}\right]}^T$$

(1)

where $\mathit{R}$ describes the translation of a body reference coordinate system with respect to the global coordinate system and $\phi$ is the rotation angle of the body reference system and $n_b$ is the number of bodies. Accordingly, generalized velocities can be written as

$$\stackrel{.}{\mathit{q}}={\left[\begin{array}{ccccc}{\stackrel{.}{\mathit{R}}}_1^T& {\omega }_1& \cdots & {\stackrel{.}{\mathit{R}}}_{n_b}^T& {\omega }_{n_b}\end{array}\right]}^T$$

(2)

where $\stackrel{·}{\mathit{R}}$ describes the velocity of a body reference coordinate system with respect to the global coordinate system and $\omega$ is the angular velocity of the body reference system.

Non-penetration contact constraints

If two rigid bodies come into contact, they should not penetrate. This means that any two bodies that are closer than a prescribed distance are considered to produce an active contact event. Accordingly, it is assumed that gap function ${\Phi }_i$ exists and can be used to describe the non-penetration constraint as follows

$$\begin{array}{cc}{\Phi }_i>0& \mathrm{if}\mathrm{two}\mathrm{bodies}\mathrm{share}\mathrm{no}\mathrm{point}\\ {\Phi }_i=0& \mathrm{if}\mathrm{two}\mathrm{bodies}\mathrm{share}\mathrm{one}\mathrm{point}\\ {\Phi }_i<0& \mathrm{if}\mathrm{two}\mathrm{bodies}\mathrm{are}\mathrm{in}\mathrm{penetration}\end{array}$$

(3)

where i represents the ith contact.

Modern collision detection algorithms can be utilized to define the contact points of the bodies with arbitrary shapes. Different methods of contact detection for several shape descriptors in a 2D and 3D case are proposed in Hogue, ²⁵ such as polygons, ellipses, and discrete function representations. It is important to note that the description of contact points for bodies with arbitrary shapes is often a challenging task. For example, there may be multiple contact points or the shape of the body may be concave, making it impossible to define the gap function between contact points.

In this study, circular bodies are considered for the sake of simplicity, whose gap function description is shown in Figure 1. Accordingly, the non-penetration constraint becomes ${\Phi }_i\ge 0$ .

Figure 1.

Gap function descriptions between two circular bodies.

Coulomb friction model

The frictional contact model used in this article is based on the Coulomb dry friction model. For a contact event i, when the contact is active between bodies, that is ${\Phi }_i=0$ , normal contact force $f_2^{(i)}$ and tangential force $f_1^{(i)}$ may exist at the contact point. Both normal contact forces and tangential contact forces can be described by set-valued force laws. It is for this reason that for any contact event i, the friction force can assume any value between 0 and $\mu f_2^{(i)}$ . Therefore, the contact forces follow two conditions

$$f_2^{\left(i\right)}>0$$

(4a)

$$f_1^{\left(i\right)}\le \mu f_2^{\left(i\right)}$$

(4b)

As for the normal contact force, if there is no contact $({\Phi }_i>0)$ at the contact point, the normal contact force must be zero $(f_2^{(i)}=0)$ . On the other hand, if contact exists $({\Phi }_{\mathrm{i}}=0)$ , the normal contact force $(f_2^{(i)}>0)$ exists as well.

The frictional force at the contact point changes instantaneously from sticking to sliding. If the contact model is sliding^26,27—that is, the relative tangential velocity $v_T$ at the contact point is not zero—the frictional contact force and normal contact force fulfill the relationship $f_1^{(i)}=\mu f_2^{(i)}$ . While the stick happens at the contact point, relative tangential velocity $v_T$ is zero; the friction model follows $f_1^{(i)}<\mu f_2^{(i)}$ . As soon as the normal contact force $f_2^{(i)}$ happens at the contact point, which is $f_2^{(i)}\ge 0$ , the contact forces can be expressed as a Coulomb friction law.

For a contact event i, when assuming ${\mathit{n}}_i$ to be the normal vector at the contact point which points toward the exterior of the body, ${\mathbf{\tau }}_i$ is the tangential vectors at the contact point. Accordingly, ${\mathit{n}}_i$ , ${\mathbf{\tau }}_i$ are unit orthogonal vectors. The reaction force on the contact point can be expressed by means of multipliers ${\hat{\gamma }}_{i,n}\ge 0$ , whereas ${\hat{\gamma }}_{i,u}$ and ${\hat{\gamma }}_{i,w}$ can have an arbitrary value.

Here, ${\mathit{f}}_1^{(i)}={\mathbf{\tau }}_i{\hat{\gamma }}_{i,\tau }$ , ${\mathit{f}}_2^{(i)}={\mathit{n}}_i{\hat{\gamma }}_{i,n}$ . So, for each contact i, both situations can be displayed by two equalities and one complementarity condition ²⁸

$${\hat{\gamma }}_{i,n}\ge 0,{\Phi }_i\ge 0$$

(5a)

$${\Phi }_i{\hat{\gamma }}_{i,n}=0$$

(5b)

$$\mu {\hat{\gamma }}_{i,n}-{\hat{\gamma }}_{i,\tau }\ge 0,\left|v_T\right|\ge 0$$

(5c)

$$\left|v_T\right|(\mu {\hat{\gamma }}_{i,n}-{\hat{\gamma }}_{i,\tau })=0$$

(5d)

Notation problem setup

In each configuration at time $(t)$ , the collection of $N_K$ contacts is denoted by $\mathcal{A}(\mathit{q}(t),{\Phi }_i)$ . Here, $N_K$ describes the total contact frequency of the bodies. As shown in Figure 2, it can be assumed that two bodies of index A and B are in contact in $n_b$ bodies. So, for contact event i, a collision detection process produces the point of contact P, a signed distance function ${\Phi }_i$ , and a set of two orthonormal vectors ${\mathit{n}}_i$ and ${\mathbf{\tau }}_i$ at the contact plane. The vector ${\mathit{n}}_i$ is normal with respect to the contact point, which leads from the body of lower index A to the body of high index B, and vector ${\mathbf{\tau }}_i$ is tangential with respect to the contact point.

Figure 2.

Illustration of rigid body A and B in contact. An arbitrary contact point is located via local constant vectors ${\mathit{s}}_{i,A}$ and ${\mathit{s}}_{i,B}$ . A set of two orthonormal vectors ${\mathit{n}}_i$ and ${\mathbf{\tau }}_i$ is generated at the contact point.

If the contact is active, then ${\Phi }_i=0$ . So, at the contact point, the force acting on body A at point P is as follows

$${\mathit{F}}_{i,A}={\mathbf{\tau }}_i{\hat{\gamma }}_{i,\tau }+{\mathit{n}}_i{\hat{\gamma }}_{i,n}=\left[\begin{array}{cc}{\mathbf{\tau }}_i\hfill & {\mathit{n}}_i\hfill \end{array}\right]\left[\begin{array}{c}{\hat{\gamma }}_{i,\tau }\hfill \\ \multicolumn{1}{c}{{\hat{\gamma }}_{i,n}}\end{array}\right]={\mathit{A}}_i{\hat{\gamma }}_i$$

(6)

where ${\mathit{A}}_i=[\begin{array}{cc}{\mathbf{\tau }}_i\hfill & {\mathit{n}}_i\hfill \end{array}]$ is the orthogonal rotation matrix consisting of unit orthogonal tangential and normal vectors at the ith contact point. The reaction force is imposed on the system by means of multipliers ${\hat{\gamma }}_{i,n}\mathrm{and}{\hat{\gamma }}_{i,\tau }$ ; that is, the normal component of the force is ${\mathit{f}}_2={\mathit{n}}_i{\hat{\gamma }}_{i,n}$ and the tangential component of the force is ${\mathit{f}}_1={\mathbf{\tau }}_i{\hat{\gamma }}_{i,\tau }$ . Here

$${\hat{\gamma }}_i=\left[\begin{array}{c}{\hat{\gamma }}_{i,\tau }\hfill \\ \multicolumn{1}{c}{{\hat{\gamma }}_{i,n}}\end{array}\right]$$

is the multiplier vector and is not known at time $(t)$ .

The position of the contact point P on body A can be written as

$${\mathit{r}}_A^P={\mathit{R}}_A+{\mathit{A}}_A{\overline{\mathit{s}}}_{i,A}$$

(7)

where ${\overline{\mathit{s}}}_{i,A}$ is the position vector of contact point P with respect to the body reference coordinate system of body A. ${\mathit{A}}_A$ is the orientation matrix of body A and can be written as

$${\mathit{A}}_A=\left[\begin{array}{cc}\cos {\phi }_A\hfill & -\sin {\phi }_A\hfill \\ \multicolumn{1}{c}{\sin {\phi }_A}& \cos {\phi }_A\hfill \end{array}\right]$$

(8)

The virtual displacement of the contact point is

$$\delta {\mathit{r}}_A^P=\delta ({\mathit{R}}_A+{\mathit{A}}_A{\overline{\mathit{s}}}_{i,A})=\delta {\mathit{R}}_A+{\mathit{A}}_{,{\phi }_A}{\overline{\mathit{s}}}_{i,A}\delta {\phi }_A$$

(9)

where the vector $\delta {\phi }_A$ is the virtual rotation associated with body A, and ${\mathit{A}}_{,{\phi }_A}$ is the derivative of ${\mathit{A}}_A$ with respect to ${\phi }_A$ .

Here, ${\mathit{A}}_{,{\phi }_A}{\overline{\mathit{s}}}_{i,A}=-{\mathit{A}}_A{\widetilde{\overline{\mathit{s}}}}_{i,A}$ , where ~ is the tilde operator, which can be explained as follows

$${\overline{\mathit{s}}}_{i,A}={\left[\begin{array}{cc}{\overline{s}}_{i,A}^x\hfill & {\overline{s}}_{i,A}^y\hfill \end{array}\right]}^T{\widetilde{\overline{\mathit{s}}}}_{i,A}={\left[\begin{array}{cc}-{\overline{s}}_{i,A}^y\hfill & {\overline{s}}_{i,A}^x\hfill \end{array}\right]}^T$$

(10)

So, the virtual displacement can be rewritten as

$$\delta {\mathit{r}}_A^P=\delta {\mathit{R}}_A-{\mathit{A}}_A{\widetilde{\overline{\mathit{s}}}}_{i,A}\delta {\phi }_A$$

(11)

The virtual work associated with the contact force ${\mathit{F}}_{i,A}$ can now be expressed as

$$\delta {\mathit{W}}_{i,A}={\left[\delta {\mathit{r}}_A^P\right]}^T{\mathit{F}}_{i,A}={[\delta {\mathit{R}}_A-{\mathit{A}}_A{\widetilde{\overline{\mathit{s}}}}_{i,A}\delta {\phi }_A]}^T{\mathit{A}}_i{\hat{\gamma }}_i=\delta {\mathit{R}}_A^T{\mathit{A}}_i{\hat{\gamma }}_i+\delta {\phi }_A^T{\widetilde{\overline{\mathit{s}}}}_{i,A}{\mathit{A}}_A^T{\mathit{A}}_i{\hat{\gamma }}_i$$

(12)

Correspondingly, the virtual work for body B is

$$\begin{array}{c}\delta {\mathit{W}}_{i,B}={\left[\delta {\mathit{r}}_B^P\right]}^T(-{\mathit{F}}_{i,A})\hfill \\ ={[\delta {\mathit{R}}_B-{\mathit{A}}_B{\widetilde{\overline{\mathit{s}}}}_{i,B}\delta {\phi }_B]}^T(-{\mathit{A}}_i{\hat{\gamma }}_i)\hfill \end{array}$$

(13)

Then, the virtual work that the presence of the frictional contact force ${\mathit{F}}_{i,A}$ imparts is

$$\begin{array}{c}\delta {\mathit{W}}_i=\delta {\mathit{W}}_{i,A}+\delta {\mathit{W}}_{i,B}=\delta {\mathit{R}}_A^T{\mathit{A}}_i{\hat{\gamma }}_i+\delta {\phi }_A^T{\widetilde{\overline{\mathit{s}}}}_{i,A}{\mathit{A}}_A^T{\mathit{A}}_i{\hat{\gamma }}_i-\delta {\mathit{R}}_B^T{\mathit{A}}_i{\hat{\gamma }}_i-\delta {\phi }_B^T{\widetilde{\overline{\mathit{s}}}}_{i,B}{\mathit{A}}_B^T{\mathit{A}}_i{\hat{\gamma }}_i\hfill \\ =(\delta {\mathit{R}}_A^T+\delta {\phi }_A^T{\widetilde{\overline{\mathit{s}}}}_{i,A}{\mathit{A}}_A^T-\delta {\mathit{R}}_B^T-\delta {\phi }_B^T{\widetilde{\overline{\mathit{s}}}}_{i,B}{\mathit{A}}_B^T){\mathit{A}}_i{\hat{\gamma }}_i\hfill \\ =\delta {\mathit{q}}^T{\mathit{D}}_i{\hat{\gamma }}_i\hfill \end{array}$$

(14)

where

$$\delta \mathit{q}=\left[\begin{array}{c}\delta {\mathit{R}}_1\hfill \\ \multicolumn{1}{c}{\delta {\phi }_1}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\delta {\mathit{R}}_{n_b}}\\ \multicolumn{1}{c}{\delta {\phi }_{n_b}}\end{array}\right]{\mathit{D}}_i=\left[\begin{array}{c}{0}_{3\times 2}\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{0_{3\times 2}}\\ \multicolumn{1}{c}{+{\mathit{A}}_i}\\ \multicolumn{1}{c}{+{\widetilde{\overline{\mathit{s}}}}_{i,A}{\mathit{A}}_A^T{\mathit{A}}_i}\\ \multicolumn{1}{c}{0_{3\times 2}}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{0_{3\times 2}}\\ \multicolumn{1}{c}{-{\mathit{A}}_i}\\ \multicolumn{1}{c}{-{\widetilde{\overline{\mathit{s}}}}_{i,B}{\mathit{A}}_B^T{\mathit{A}}_i}\\ \multicolumn{1}{c}{0_{3\times 2}}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{0_{3\times 2}}\end{array}\right]$$

(15)

Therefore, the frictional contact force associated with the presence of $N_K$ contact events is ${\mathit{D}}_i{\hat{\gamma }}_i$ , which means that ${\mathit{D}}_i$ is the contact transformation matrix associated with contact $\mathcal{A}(\mathit{q}(t),{\Phi }_i)$ , which can convert the general force to contact force

$$\mathit{D}=\left[\begin{array}{ccc}{\mathit{D}}_1\hfill & \cdots \hfill & {\mathit{D}}_{N_K}\hfill \end{array}\right]$$

(16)

Finally, note that

$$\begin{array}{c}{\mathit{D}}_i^T\stackrel{·}{\mathit{q}}={\left[\begin{array}{c}{0}_{3\times 2}\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{0_{3\times 2}}\\ \multicolumn{1}{c}{+{\mathit{A}}_i}\\ \multicolumn{1}{c}{+{\widetilde{\overline{\mathit{s}}}}_{i,A}{\mathit{A}}_A^T{\mathit{A}}_i}\\ \multicolumn{1}{c}{0_{3\times 2}}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{0_{3\times 2}}\\ \multicolumn{1}{c}{-{\mathit{A}}_i}\\ \multicolumn{1}{c}{-{\widetilde{\overline{\mathit{s}}}}_{i,B}{\mathit{A}}_B^T{\mathit{A}}_i}\\ \multicolumn{1}{c}{0_{3\times 2}}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{0_{3\times 2}}\end{array}\right]}^T\left[\begin{array}{c}{\stackrel{·}{\mathit{R}}}_1^T\hfill \\ \multicolumn{1}{c}{{\omega }_1}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\stackrel{·}{\mathit{R}}}_A^T}\\ \multicolumn{1}{c}{{\omega }_A}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\stackrel{·}{\mathit{R}}}_B^T}\\ \multicolumn{1}{c}{{\omega }_B}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\stackrel{·}{\mathit{R}}}_{n_b}^T}\\ \multicolumn{1}{c}{{\omega }_{n_b}}\end{array}\right]\hfill \\ ={\mathit{A}}_i^T{\stackrel{·}{\mathit{R}}}_A+{\mathit{A}}_i^T{\mathit{A}}_A{\widetilde{\overline{\mathit{s}}}}_{i,A}{\omega }_A-{\mathit{A}}_i^T{\stackrel{·}{\mathit{R}}}_B-{\mathit{A}}_i^T{\mathit{A}}_B{\widetilde{\overline{\mathit{s}}}}_{i,B}{\omega }_B\hfill \\ ={\mathit{A}}_i^T({\stackrel{·}{\mathit{R}}}_A+{\mathit{A}}_{\phi ,A}{\overline{\mathit{s}}}_{i,A}{\omega }_A-{\stackrel{·}{\mathit{R}}}_B-{\mathit{A}}_{\phi ,B}{\overline{\mathit{s}}}_{i,B}{\omega }_B)\hfill \\ ={\mathit{A}}_i^T({\stackrel{·}{\mathit{r}}}_{i,A}-{\stackrel{·}{\mathit{r}}}_{i,B})=\left[\begin{array}{c}{v}_{i,\tau }\hfill \\ \multicolumn{1}{c}{v_{i,n}}\end{array}\right]\hfill \end{array}$$

(17)

where

$$\left[\begin{array}{c}{v}_{i,\tau }\hfill \\ \multicolumn{1}{c}{v_{i,n}}\end{array}\right]$$

represents the relative velocity at the contact point between two bodies; that is, $v_{i,\tau }$ is the tangential relative velocity and $v_{i,n}$ is the normal relative velocity. The matrix ${\mathit{D}}_i$ can also convert the general velocity to contact velocity.

Equations of motion

Assume that there are $N_K$ potential contacts, so that the contact constraints are enforced by non-penetration constraints ${\Phi }_i\ge 0,i=1,2,\dots ,p$ . Here, superscript i is the number of contact event. Therefore, the equations of motion take the form

$$\mathit{M}\stackrel{··}{\mathit{q}}=\mathit{F}(\mathit{q},\stackrel{·}{\mathit{q}},t)+\mathit{D}\hat{\gamma }$$

(18)

where mass matrix $\mathit{M}=\mathrm{diag}\{{\mathit{M}}_1,\dots ,{\mathit{M}}_A,{\mathit{M}}_B,\dots ,{\mathit{M}}_{n_b}\}$ , where ${\mathit{M}}_i=\mathrm{diag}\{m_i,m_i,J_i\}$ and $J_i=m_i{r}_i^2/2$ . $\mathit{F}(\mathit{q},\stackrel{·}{\mathit{q}},t)$ is the generalized applied force, and $\mathit{D}\hat{\gamma }$ is the frictional contact force associated with the presence of $N_K$ contact events. Here, a matrix $\hat{\gamma }$ is built to contain all of the contact force as in equation (19)

$$\hat{\gamma }={\left[\begin{array}{cccc}{\hat{\gamma }}_1\hfill & {\hat{\gamma }}_2\hfill & \cdots \hfill & {\hat{\gamma }}_{N_K}\hfill \end{array}\right]}^T$$

(19)

Discretized equations of motion

Using a semi-implicit Euler numerical scheme ¹⁸ at time step $t^{(l+1)}=t^{(l)}+\mathrm{\Delta }t$ , the equations of motion can be expressed as follows

$$\stackrel{\mathrm{generalized}\mathrm{positions}}{\overbrace{\mathit{q}(l+1)}}={\mathit{q}}^{(l)}+\stackrel{\mathrm{step}\mathrm{size}}{\overbrace{\mathrm{\Delta }t}}{\stackrel{·}{\mathit{q}}}^{(l+1)}$$

(20a)

$$\mathit{M}(\stackrel{\mathrm{generalized}\mathrm{speeds}}{\overbrace{{\stackrel{·}{\mathit{q}}}^{(l+1)}}}-{\stackrel{·}{\mathit{q}}}^{(l)})=\underset{\mathrm{applied}\mathrm{impulse}}{\underbrace{{\mathit{f}}^{(l)}}}+\underset{\mathrm{frictional}\mathrm{contact}\mathrm{impulse}}{\underbrace{{\mathit{D}}^{(l)}{\mathbf{\gamma }}^{(l+1)}}}$$

(20b)

$$i\in {\mathcal{A}}^{(l)}(\mathit{q}\left(t\right),{\Phi }_i):\{\begin{array}{c}{\gamma }_{i,n}^{(l+1)}\ge 0\hfill \\ \multicolumn{1}{c}{({\mu }_i{\gamma }_{i,n}^{(l+1)}-{\gamma }_{i,\tau }^{(l+1)})\ge 0}\end{array}$$

(20c)

$${\Phi }_i\ge 0$$

(20d)

where ${\mathit{f}}^{(l)}=\mathrm{\Delta }t\mathit{F}({\mathit{q}}^{(l)},{\stackrel{·}{\mathit{q}}}^{(l)},t^{(l)})$ , ${\mathbf{\gamma }}^{(l+1)}=\mathrm{\Delta }t{\hat{\gamma }}^{(l+1)}$ and ${\mathit{D}}^{(l)}=\mathit{D}({\mathit{q}}^{(l)},t^{(l)})$ , ${\stackrel{·}{\mathit{q}}}^{(l)}$ is the general known velocity and ${\stackrel{·}{\mathit{q}}}^{(l+1)}$ is the unknown velocity with time step $\mathrm{\Delta }t$ . In equation (20c), ${\mathcal{A}}^{(l)}(\mathit{q}(t),{\Phi }_i)$ is the set of active contact events produced by the collision detection step carried out at $t^{(l)}$ in the configuration ${\mathit{q}}^{(l)}$ . According to equation (4), the Coulomb friction model used in this article is presented in equation (20c).

An approximation of the signed gap function at $t^{(l+1)}$ is utilized to reflect the complementarity

$$\exists i\in {\mathcal{A}}^{(l)}(\mathit{q}\left(t\right),{\Phi }_i):{\Phi }_i^{(l+1)}\approx {\Phi }_i^{(l)}+\mathrm{\Delta }t{\nu }_{i,n}^{(l+1)}$$

(21)

Reformulation as a CCP

Using the force balance condition in equation (19b), the velocity within time step ${\stackrel{·}{\mathit{q}}}^{(l+1)}$ can be calculated as

$${\stackrel{·}{\mathit{q}}}^{(l+1)}={\stackrel{·}{\mathit{q}}}^{(l)}+{\mathit{M}}^{-1}{\mathit{f}}^{(l)}+{\mathit{M}}^{-1}{\mathit{D}}^{(l)}{\mathbf{\gamma }}^{(l+1)}$$

(22)

Next, ${\mathit{d}}_i$ can be defined as

$${\mathit{d}}_i=\left[\begin{array}{c}{v}_{i,\tau }^{(l+1)}\hfill \\ \multicolumn{1}{c}{\frac{1}{\mathrm{\Delta }t}{\Phi }_i^{(l)}+v_{i,n}^{(l+1)}}\end{array}\right]$$

(23)

At the initial time, the velocity at the contact is all zero, so the initial term of ${\mathit{d}}_i$ is

$${\mathit{d}}_{i,0}=\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{\frac{1}{\mathrm{\Delta }t}{\Phi }_i^{(l)}}\end{array}\right]$$

(24)

According to equation (17), the matrix $\mathit{D}$ can convert the general velocity to contact velocity, and ${\mathit{d}}_i$ can also be presented as

$${\mathit{d}}_i=\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{\frac{1}{\mathrm{\Delta }t}{\Phi }_i^{(l)}}\end{array}\right]+\left[\begin{array}{c}{v}_{i,\tau }^{(l+1)}\hfill \\ \multicolumn{1}{c}{v_{i,n}^{(l+1)}}\end{array}\right]={\mathit{d}}_{i,0}+{\mathit{D}}^{\left(l\right),T}{\stackrel{·}{\mathit{q}}}^{(l+1)}$$

(25)

Then, equation (22) can be submitted into equation (25) as

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{d}}_i={\mathit{d}}_{i,0}+{\mathit{D}}^{\left(l\right),T}({\stackrel{·}{\mathit{q}}}^{(l)}+{\mathit{M}}^{-1}{\mathit{f}}^{(l)}+{\mathit{M}}^{-1}{\mathit{D}}^{(l)}{\mathbf{\gamma }}^{(l+1)})}\\ \multicolumn{1}{c}{={\mathit{d}}_{i,0}+{\mathit{D}}^{\left(l\right),T}({\stackrel{·}{\mathit{q}}}^{(l)}+{\mathit{M}}^{-1}{\mathit{f}}^{(l)})+{\mathit{D}}^{\left(l\right),T}{\mathit{M}}^{-1}{\mathit{D}}^{(l)}{\mathbf{\gamma }}^{(l+1)}}\end{array}$$

(26)

Here, ${\mathit{d}}_{i,1}$ can be defined as

$${\mathit{d}}_{i,1}={\mathit{d}}_{i,0}+{\mathit{D}}^{\left(l\right),T}({\stackrel{·}{\mathit{q}}}^{(l)}+{\mathit{M}}^{-1}{\mathit{f}}^{(l)})$$

(27)

Therefore, ${\mathit{d}}_i$ can be rewritten as

$${\mathit{d}}_i={\mathit{d}}_{i,1}+{\mathit{D}}^{\left(l\right),T}{\mathit{M}}^{-1}{\mathit{D}}^{(l)}{\mathbf{\gamma }}^{(l+1)}$$

(28)

Below, the vector p will be used as a union of all of the particular vectors ${\mathit{d}}_{i,1}$ .

$$\mathit{p}=\left[\begin{array}{c}{\mathit{d}}_{1,1}\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\mathit{d}}_{N_K,1}}\end{array}\right]$$

(29)

Analogously, also matrix $\mathit{N}$ can be introduced

$$\mathit{N}={\mathit{D}}^T{\mathit{M}}^{-1}\mathit{D}$$

(30)

Quadratic optimization problem

The CCP represents the first-order optimality conditions ¹¹ for the convex quadratic optimization problem with conic constraints

$$\begin{array}{c}{\mathbf{\gamma }}^{(l+1)}=\underset{\gamma }{min}\frac{1}{2}{\mathbf{\gamma }}^T\mathit{N}\mathbf{\gamma }+{\mathit{p}}^T\mathbf{\gamma }\hfill \\ \mathrm{subject}\mathrm{to}{\mathbf{\gamma }}_k\in {\mathit{C}}_k\hfill \end{array}$$

(31)

where $\mathbf{\gamma }$ is the contact impulse when the bodies collide, ${\mathbf{\gamma }}^{(l+1)}$ is the unknown contact impulse at the next time step $\mathrm{\Delta }t$ , and ${\mathit{C}}_k$ is from the Coulomb friction model equation (4).

With time integration, the new contact impulse ${\mathbf{\gamma }}^{(l+1)}$ will be computed using equation (31). Meanwhile, the new velocity ${\stackrel{·}{\mathit{q}}}^{(l+1)}$ can be obtained via equation (22). The new position ${\mathit{q}}^{(l+1)}$ will finally be calculated through equation (20a).

Dynamic simulation of multiple pendulum interactions

In this section, the contact description explained above is applied to study multiple pendulums. Multiple bodies may come into contact with each other in three possible situations: (1) the bodies hit the ground; (2) the bodies collide with each other; (3) adjacent bodies are restricted by a rope. Constraints and contact scenarios used are explained in the following sections of the paper.

Case of bodies hitting the ground

The first type of multibody contact is between bodies and the ground, as shown in Figure 3.

Figure 3.

Case of body hitting the ground. The gap function ${\Phi }^{G_{(i)}}$ is the distance between the contact point and the ground. Two orthogonal forces $f_1^{G_{(i)}}$ , $f_2^{G_{(i)}}$ are produced at the contact point. The gap function and the normal contact force follow the complementarity condition ${\Phi }^{G_{(i)}}{f}_2^{G_{(i)}}=0$ .

Each contact point results in two orthogonal forces $f_1^{G_{(i)}}\mathrm{and}{f}_2^{G_{(i)}}$ and one torque $T^{G_{(i)}}$ . The torque can be calculated as

$$T^{G_{\left(i\right)}}=r_i{f}_1^{G_{\left(i\right)}}+0f_2^{G_{\left(i\right)}}=\left[\begin{array}{cc}{r}_i\hfill & 0\hfill \end{array}\right]\left[\begin{array}{c}{f}_1^{G_{\left(i\right)}}\hfill \\ \multicolumn{1}{c}{f_2^{G_{\left(i\right)}}}\end{array}\right]$$

(32)

Therefore, the external force and torque can be calculated as

$$\left[\begin{array}{c}{f}_x^{G_{\left(i\right)}}\hfill \\ \multicolumn{1}{c}{f_y^{G_{\left(i\right)}}}\\ \multicolumn{1}{c}{T^{G_{\left(i\right)}}}\end{array}\right]=\left[\begin{array}{cc}1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \\ \multicolumn{1}{c}{r_i}& 0\hfill \end{array}\right]\left[\begin{array}{c}{f}_1^{G_{\left(i\right)}}\hfill \\ \multicolumn{1}{c}{f_2^{G_{\left(i\right)}}}\end{array}\right]$$

(33)

Here, the matrix ${\mathit{D}}^{G_{(i)}}$ is

$${\mathit{D}}^{G_{\left(i\right)}}=\left[\begin{array}{cc}1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \\ \multicolumn{1}{c}{r_i}& 0\hfill \end{array}\right]$$

(34)

Each body can potentially make contact with the ground. Therefore, matrix ${\mathit{D}}_i^G$ takes the following form

$${\mathit{D}}^G=\left[\underset{n_b\mathrm{times}}{\underbrace{\begin{array}{cccc}{\mathit{D}}^{G_{\left(1\right)}}\hfill & 0_{3\times 2}\hfill & \cdots \hfill & 0_{3\times 2}\hfill \\ \multicolumn{1}{c}{0_{3\times 2}}& {\mathit{D}}^{G_{\left(2\right)}}\hfill & \cdots \hfill & 0_{3\times 2}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{0_{3\times 2}}& 0_{3\times 2}\hfill & \cdots \hfill & {\mathit{D}}^{G_{\left(n_b\right)}}\hfill \end{array}}}\right]$$

(35)

The gap function ${\Phi }^{G_{(i)}}$ is the distance between the contact point and the ground; that is, the number of gap functions is the number of bodies $n_b$

$${\Phi }^{G_{\left(i\right)}}=y_i-r_i$$

(36)

where $y_i$ is the distance between the center of the body and the ground, and $r_i$ is the radius of the body i.

Case of bodies colliding with each other

As shown in Figure 4, when bodies collide with each other, two contact forces appear the normal force $f_2^{S_{(\mathrm{ij})}}$ and the tangential force $f_1^{S_{(\mathrm{ij})}}$ . The normal force $f_2^{S_{(\mathrm{ij})}}$ is directed along the line between two body centers.

Figure 4.

Case of bodies colliding with each other. The gap function ${\Phi }^{S_{(\mathrm{ij})}}$ is the distance between the contact point of the bodies. Two orthogonal forces $f_1^{S_{(\mathrm{ij})}}\mathrm{and}{f}_2^{S_{(\mathrm{ij})}}$ are produced at the contact point. The gap function and the normal contact force follow the complementarity condition ${\Phi }^{S_{(\mathrm{ij})}}{f}_2^{S_{(\mathrm{ij})}}=0$ .

The force acting on body i at point P is

$${\mathit{F}}^{S_{\left(\mathrm{ij}\right)}}={\mathbf{\tau }}_{\mathrm{ij}}{\hat{\gamma }}_1^{\left(\mathrm{ij}\right)}+{\mathit{n}}_{\mathrm{ij}}{\hat{\gamma }}_2^{\left(\mathrm{ij}\right)}={\mathit{A}}_{\mathrm{ij}}{\hat{\gamma }}_{\mathrm{ij}}$$

(37)

The distance of the center mass of two contact bodies can be written as $L=\sqrt{{(x_i-x_j)}^2+{(y_i-y_j)}^2}$ . For this situation, the orthogonal matrix ${\mathit{A}}_{\mathrm{ij}}$ is composed with the normal vector ${\mathit{n}}_{\mathrm{ij}}$ and the tangential vector ${\mathbf{\tau }}_{\mathrm{ij}}$ . The normal vector ${\mathit{n}}_{\mathrm{ij}}$ is the unit vector which is directed from the body j toward body i. Thus, it is calculated as

$${\mathit{n}}_{\mathrm{ij}}=\left[\begin{array}{c}\frac{x_i-x_j}{L}\hfill \\ \multicolumn{1}{c}{\frac{y_i-y_j}{L}}\end{array}\right]$$

The tangential vector ${\mathbf{\tau }}_{\mathrm{ij}}$ is orthogonal to ${\mathit{n}}_{\mathrm{ij}}$ , so it should be calculated as

$${\mathbf{\tau }}_{\mathrm{ij}}=\left[\begin{array}{c}\frac{y_i-y_j}{L}\hfill \\ \multicolumn{1}{c}{\frac{-x_i+x_j}{L}}\end{array}\right]$$

Accordingly, matrix ${\mathit{A}}_{\mathrm{ij}}$ takes the form

$${\mathit{A}}_{\mathrm{ij}}=\left[\begin{array}{cc}{\mathbf{\tau }}_{\mathrm{ij}}\hfill & {\mathit{n}}_{\mathrm{ij}}\hfill \end{array}\right]=\left[\begin{array}{cc}\frac{y_i-y_j}{L}\hfill & \frac{x_i-x_j}{L}\hfill \\ \multicolumn{1}{c}{\frac{-x_i+x_j}{L}}& \frac{y_i-y_j}{L}\hfill \end{array}\right]$$

(38)

From equation (15), it can be concluded that

$${\mathit{D}}^{S_{\left(\mathrm{ij}\right)}}=\left[\begin{array}{c}+{\mathit{A}}_{\mathrm{ij}}\hfill \\ \multicolumn{1}{c}{+r_i{\mathbf{\tau }}_{\mathrm{ij}}^T{\mathit{A}}_{\mathrm{ij}}}\end{array}\right]{\mathit{D}}^{S_{\left(\mathrm{ji}\right)}}=\left[\begin{array}{c}-{\mathit{A}}_{\mathrm{ij}}\hfill \\ \multicolumn{1}{c}{-r_j{\mathbf{\tau }}_{\mathrm{ij}}^T{\mathit{A}}_{\mathrm{ij}}}\end{array}\right]$$

(39)

Considering that each body can potentially collide with all others, the number of contacts is $n_c=n_b(n_b-1)/2$ . That is why the full Jacobian matrix for this case is

$${\mathit{D}}^S=\left[\underset{\frac{n_b(n_b-1)}{2}\mathrm{times}}{\underbrace{\begin{array}{cccc}{\mathit{D}}^{S_{\left(12\right)}}\hfill & {\mathit{D}}^{S_{\left(13\right)}}\hfill & \cdots \hfill & 0_{3\times 2}\hfill \\ \multicolumn{1}{c}{{\mathit{D}}^{S_{\left(21\right)}}}& 0_{3\times 2}\hfill & \cdots \hfill & 0_{3\times 2}\hfill \\ \multicolumn{1}{c}{0_{3\times 2}}& {\mathit{D}}^{S_{\left(31\right)}}\hfill & \cdots \hfill & 0_{3\times 2}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{0_{3\times 2}}& 0_{3\times 2}\hfill & \cdots \hfill & {\mathit{D}}^{S_{(n_b-1)n_b}}\hfill \\ \multicolumn{1}{c}{0_{3\times 2}}& 0_{3\times 2}\hfill & \cdots \hfill & {\mathit{D}}^{S_{n_b(n_b-1)}}\hfill \end{array}}}\right]$$

(40)

The gap function ${\Phi }^{S_{(\mathrm{ij})}}$ expresses the distance between the two contact bodies, and the number of such functions is also $n_c$

$${\Phi }^{S_{\left(\mathrm{ij}\right)}}=L-r_i-r_j$$

(41)

where $r_i$ and $r_j$ are the radii of the body i and body j, respectively.

Rope inextensibility constraint

The adjacent bodies are connected pair-wise by a rope. Therefore, the adjacent bodies are restricted by the rope inextensibility constraint, as shown in Figure 5. The normal force $f_2^{L_{(\mathrm{ij})}}$ is directed along the line between the two adjacent bodies, which is also the direction of the rope.

Figure 5.

Case of rope inextensibility constraint. The gap function ${\Phi }^{L_{(\mathrm{ij})}}$ is the distance of the fully extended rope. The normal force $f_2^{L_{(\mathrm{ij})}}$ is produced along the rope. The gap function and the normal contact force follow the complementarity condition ${\Phi }^{L_{(\mathrm{ij})}}{f}_2^{L_{(\mathrm{ij})}}=0$ .

The force acting on body i at the contact point is

$${\mathit{F}}^{L_{\left(\mathrm{ij}\right)}}={\mathit{n}}_{\mathrm{ij}}{\hat{\gamma }}_2^{\left(\mathrm{ij}\right)}={\mathit{A}}_{\mathrm{ij}}{\hat{\gamma }}_{\mathrm{ij}}$$

(42)

where ${\mathit{A}}_{\mathrm{ij}}=[{\mathit{n}}_{\mathrm{ij}}]$ is the unit vector of the direction of the rope.

The normal vector ${\mathit{n}}_{\mathrm{ij}}$ is the unit vector which is oriented from body j to body i and it should be calculated as

$${\mathit{n}}_{\mathrm{ij}}=\left[\begin{array}{c}\frac{x_i-x_j}{L}\hfill \\ \multicolumn{1}{c}{\frac{y_i-y_j}{L}}\end{array}\right]$$

Thus, matrix ${\mathit{A}}_{\mathrm{ij}}$ is

$${\mathit{A}}_{\mathrm{ij}}=\left[{\mathit{n}}_{\mathrm{ij}}\right]=\left[\begin{array}{c}\frac{x_i-x_j}{L}\hfill \\ \multicolumn{1}{c}{\frac{y_i-y_j}{L}}\end{array}\right]$$

(43)

It can be deduced that

$${\mathit{D}}^{L_{\left(\mathrm{ij}\right)}}=\left[\begin{array}{c}+{\mathit{A}}_{\mathrm{ij}}\hfill \\ \multicolumn{1}{c}{0}\end{array}\right]{\mathit{D}}^{L_{\left(\mathrm{ji}\right)}}=\left[\begin{array}{c}-{\mathit{A}}_{\mathrm{ij}}\hfill \\ \multicolumn{1}{c}{0}\end{array}\right]$$

(44)

The number of contacts of this kind in one branch is equal to $n_b/2$ , so the total number is $n_b$ . ${\mathit{D}}^{L_{(01)}}$ and ${\mathit{D}}^{L_{(n_b/2(n_b/2+1))}}$ represent the constraints between the top body and the ceiling. The general form of the Jacobian matrix for this type of constraint is as follows

$$D^L=\left[\underset{n_b\mathrm{times}}{\underbrace{\begin{array}{cccccccc}{\mathbf{D}}^{L_{\left(01\right)}}& {\mathbf{D}}^{L_{\left(12\right)}}& \cdots & {\mathbf{0}}_{3\times 1}& {\mathbf{0}}_{3\times 1}& {\mathbf{0}}_{3\times 1}& \cdots & {\mathbf{0}}_{3\times 1}\\ {\mathbf{0}}_{3\times 1}& {\mathbf{D}}^{L_{\left(21\right)}}& \cdots & {\mathbf{0}}_{3\times 1}& {\mathbf{0}}_{3\times 1}& {\mathbf{0}}_{3\times 1}& \cdots & {\mathbf{0}}_{3\times 1}\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮\\ {\mathbf{0}}_{3\times 1}& {\mathbf{0}}_{3\times 1}& \cdots & {\mathbf{D}}^{L_{\left(\left(\frac{n_b}{2}-1\right)\frac{n_b}{2}\right)}}& {\mathbf{0}}_{3\times 1}& {\mathbf{0}}_{3\times 1}& \cdots & {\mathbf{0}}_{3\times 1}\\ {\mathbf{0}}_{3\times 1}& {\mathbf{0}}_{3\times 1}& \cdots & {\mathbf{D}}^{L_{\left(\frac{n_b}{2}\left(\frac{n_b}{2}-1\right)\right)}}& {\mathbf{0}}_{3\times 1}& {\mathbf{0}}_{3\times 1}& \cdots & {\mathbf{0}}_{3\times 1}\\ {\mathbf{0}}_{3\times 1}& {\mathbf{0}}_{3\times 1}& \cdots & {\mathbf{0}}_{3\times 1}& {\mathbf{D}}^{L_{\left(\frac{n_b}{2}\left(\frac{n_b}{2}+1\right)\right)}}& {\mathbf{D}}^{L_{\left(\left(\frac{n_b}{2}+1\right)\left(\frac{n_b}{2}+2\right)\right)}}& \cdots & {\mathbf{0}}_{3\times 1}\\ {\mathbf{0}}_{3\times 1}& {\mathbf{0}}_{3\times 1}& \cdots & {\mathbf{0}}_{3\times 1}& {\mathbf{0}}_{3\times 1}& {\mathbf{D}}^{L_{\left(\left(\frac{n_b}{2}+2\right)\left(\frac{n_b}{2}+1\right)\right)}}& \cdots & {\mathbf{0}}_{3\times 1}\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮\\ {\mathbf{0}}_{3\times 1}& {\mathbf{0}}_{3\times 1}& \cdots & {\mathbf{0}}_{3\times 1}& {\mathbf{0}}_{3\times 1}& {\mathbf{0}}_{3\times 1}& \cdots & {\mathbf{D}}^{L_{\left(\left(n_b-1\right)n_b\right)}}\\ {\mathbf{0}}_{3\times 1}& {\mathbf{0}}_{3\times 1}& \cdots & {\mathbf{0}}_{3\times 1}& {\mathbf{0}}_{3\times 1}& {\mathbf{0}}_{3\times 1}& \cdots & {\mathbf{D}}^{L_{\left(n_b\left(n_b-1\right)\right)}}\end{array}}}\right]$$

(45)

The gap function is computed as follows

$${\Phi }^{L_{\left(\mathrm{ij}\right)}}=L_{\max }-L$$

(46)

Then, put all D matrices together from equations (35), (40), and (45)

$$\mathit{D}=\left[\begin{array}{ccc}{\mathit{D}}^G\hfill & {\mathit{D}}^S\hfill & {\mathit{D}}^L\hfill \end{array}\right]$$

(47)

Solving the CCP

The overall solution scheme can be expressed with Algorithm 1.

Algorithm 1

1.  For $i=1,2,\dots n$ , the generalized position ${\mathit{q}}^{\left(l\right)}$ and generalized velocity ${\stackrel{·}{\mathit{q}}}^{\left(l\right)}$ of the system are known at the current time step as Eq.(1) and Eq. (2) show. 2.  For $i=1,2,\dots n$ , compute mass matrix $\mathit{M}$ . 3. Use Eq. (35), (40), (45), and (47) to obtain the D matrix and Eq. (36), (41) and (46) to obtain gap function $\Phi$ . 4. Use Eq. (48) and (27) to calculate ${\mathit{d}}_{i,0}$ and ${\mathit{d}}_{i,1}$ . 5. Use Eq. (30) and (29) to calculate matrix $\mathit{N}$ and $\mathit{p}$ , where ${\mathit{f}}^{\left(l\right)}$ is the applied impulse. In this situation, only the gravity force should be considered. So, ${\mathit{f}}^{\left(l\right)}={\left[\begin{array}{ccccccc}0\hfill & -m_{1}g\hfill & 0\hfill & \cdots \hfill & 0\hfill & -m_{n}g\hfill & 0\hfill \end{array}\right]}^T\mathrm{\Delta }\mathrm{t}$ . 6. Use the quadratic programming method to calculate the unknown force impulse ${\mathbf{\gamma }}^{(l+1)}$ , where matrix $\mathit{A}$ and $\mathit{b}$ can be obtained through Eq. (53) and (54) below. 7. Use Eq. (22) to update the speed ${\stackrel{·}{\mathit{q}}}^{(l+1)}$ and Eq. (20a) to update the position ${\mathit{q}}^{(l+1)}$ . Increment $t^{(l+1)}=t^{\left(l\right)}+\mathrm{\Delta }t$ , and repeat from step 3 until $t>t_{\mathrm{end}}$ .

According to equation (23), inserting the gap function which can be calculated in equation (36), (41), and (46), one can obtain ${\mathit{d}}_{i,0}$

$$\begin{array}{c}{\mathit{d}}_{i,0}=\hfill \\ \left[\underset{2*n_b\mathrm{times}}{\underbrace{\begin{array}{cccc}0\hfill & \frac{1}{\mathrm{\Delta }t}{\Phi }_1^G\hfill & \cdots \hfill & \frac{1}{\mathrm{\Delta }t}{\Phi }_n^G\hfill \end{array}}}\underset{2*\frac{n_b(n_b-1)}{2}\mathrm{times}}{\underbrace{\begin{array}{cccc}0\hfill & \frac{1}{\mathrm{\Delta }t}{\Phi }_{12}^S\hfill & \cdots \hfill & \frac{1}{\mathrm{\Delta }t}{\Phi }_{(n-1)n}^S\hfill \end{array}}}\underset{n_b\mathrm{times}}{\underbrace{\begin{array}{cccc}0\hfill & \frac{1}{\mathrm{\Delta }t}{\Phi }_{12}^L\hfill & \cdots \hfill & \frac{1}{\mathrm{\Delta }t}{\Phi }_{(n-1)n}^L\hfill \end{array}}}\right]\hfill \end{array}$$

(48)

This article employs the Coulomb friction model. When bodies hit the ground or collide with each other, there are two contact forces, $f_1$ and $f_2$ . Consequently, according to equation (4)

$$\begin{array}{c}-f_2\le 0\hfill \\ f_1-\mu f_2\le 0\hfill \end{array}$$

(49)

$${\mathit{A}}^{\left(i\right)}=\left[\begin{array}{cc}0\hfill & -1\hfill \\ \multicolumn{1}{c}{1}& -\mu \hfill \end{array}\right],{\mathit{b}}^{\left(i\right)}=\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\end{array}\right]$$

(50)

When adjacent bodies are loaded with the rope inextensibility constraint, only one contact force $f_2$ exists, so according to Coulomb’s law

$$-f_2\le 0$$

(51)

$${\mathit{A}}^{\left(j\right)}=[-1],{\mathit{b}}^{\left(j\right)}=\left[0\right]$$

(52)

Combining the A matrix and b matrix yields

$$\mathit{A}=\underset{n_p\mathrm{times}}{\left[\underbrace{\begin{array}{cccccccc}0\hfill & -1\hfill & \cdots \hfill & 0\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{1}& -\mu \hfill & \cdots \hfill & 0\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill & ⋮\hfill & ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & -1\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill & -\mu \hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & -1\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & ⋮\hfill & ⋮\hfill & ⋮\hfill & ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & \cdots \hfill & -1\hfill \end{array}}\right]}$$

(53)

$$b={\underset{n_p\mathrm{times}}{\left[\underbrace{\begin{array}{ccccc}0\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill \end{array}}\right]}}^T$$

(54)

where $n_p$ is the potential contact number

$$n_p=2*n_b+2*\frac{n_b(n_b-1)}{2}+n_b=n_b^2+2*n_b$$

(55)

Performance investigation

The numerical method proposed in this article can be utilized to simulate granular pendulum contact. This problem occurs in plenty of engineering applications; for example, a tree harvester truck can be assumed and simplified into a chain structure concept, which can be divided into several arbitrary bodies. ²⁸

In this section, two types of examples have been discussed. The CCP method developed in this article has been compared with one optimization-based method proposed in Kleinert et al. ²² in the first example. The second example investigates the simulation of multiple pendulum contact with CCP solvers.

Comparison between the CCP method and the optimization-based method

The dynamic example studies one particle fall on the horizontal ground. The initial position of the particle follows $x=0\mathrm{m},y=3\mathrm{m}$ , and the initial velocity is $\stackrel{·}{x}=3\mathrm{m}/\mathrm{s},\stackrel{·}{y}=0\mathrm{m}/\mathrm{s}$ . The particle is affected by gravity force, whose acceleration is $g=9.8\mathrm{N}/\mathrm{kg}$ , and the friction coefficient is $\mu =0.3$ . The mass of the particle is 1 kg and the time step $\mathrm{\Delta }t=0.001\mathrm{s}$ .

Figure 6 presents the displacements of the x and y directions for one particle with frictional contact. Some points at $x=[\begin{array}{ccccc}0.5\hfill & 1\hfill & 1.5\hfill & 2\hfill & 2.34\hfill \end{array}]$ have been separately measured from the figure of the study by M. Anitescu and added to Figure 6. Good agreement can be observed between the proposed optimization method in Anitescu ²³ and CCP. As the figure displays, the particle stops on the ground after the contact. The coefficient of restitution is 0, which agrees with Stewart and Trinkle. ¹⁶

Figure 6.

Displacements of X and Y directions for one particle with frictional contact. The coefficient of restitution is 0. Some points at $\mathrm{x}=[\begin{array}{ccccc}0.5\hfill & 1\hfill & 1.5\hfill & 2\hfill & 2.34\hfill \end{array}]$ have been measured from the figure of the study by M Anitescu. ²³

To show how the value of the time step affects the precision of the solution, Figure 7 shows numerical solutions of the Y coordinate of one particle contact with an initial condition of $x=0\mathrm{m},y=3\mathrm{m},\stackrel{·}{x}=3\mathrm{m}/\mathrm{s},\stackrel{·}{y}=0\mathrm{m}/\mathrm{s}$ . The simulation has been calculated with different time steps from $\mathrm{\Delta }t=0.1\mathrm{s}$ to $\mathrm{\Delta }t=0.00125\mathrm{s}$ . Here, the Y coordinates converge with the decrease of the value of $\mathrm{\Delta }t$ . Figure 8 gives an error estimate of the numerical results of Y coordinates at each time step from Figure 7. The solution for $\mathrm{\Delta }t=0.00125\mathrm{s}$ has been used as a converged value. For different values of $\mathrm{\Delta }t$ , the error measures used are $\mathrm{error}=|y(\mathrm{\Delta }t)-y(\mathrm{exact})|$ . The smaller the time step is, the closer the error converge to zero.

Figure 7.

Numerical stability of the CCP approach. Y coordinate of one particle contact with different time steps from $\mathrm{\Delta }t=0.1\mathrm{s}$ to $\mathrm{\Delta }t=0.00125\mathrm{s}$ .

Figure 8.

Errors for numerical results of Y coordinates associated with Figure 7. The solution of $\mathrm{\Delta }t=0.00125\mathrm{s}$ is used as a converged value.

Comparison between the CCP method and penalty method

To compare the CCP method with the penalty method, a simple case of two pendulums with two balls illustrated in Figure 9 has been analyzed. Aim of this simple contact problem is to clarify importance of the contact parameters of the penalty method in dynamic response after contact. The parameters for the problem are given in Figure 9. The gravity is $g=9.8\mathrm{kg}/{\mathrm{m}}^3$ , and the coefficient of friction is $\mu =0.3$ . The time step is $\mathrm{\Delta }t=0.0001\mathrm{s}$ and the total time is 3 s.

Figure 9.

Initial position of the two pendulum system with two rigid balls.

The penalty method is based on dissipative contact force models, combining a linear spring with a linear damper. The contact force can be calculated as

$$F_N=K{\delta }^{3/2}+D\stackrel{·}{\delta }\delta$$

(56)

where K is the stiffness coefficient and D is the damping coefficient. Indentation $\delta$ is the distance between the contact points of the two balls. In equation (55), $\stackrel{·}{\delta }$ is the relative normal velocity between two contact points. The contact parameters of the penalty method have been chosen so that it gives agreement for vertical position solved with the CCP approach. For the penalty method, contact stiffness coefficient $K=3.6e^5\mathrm{N}/{\mathrm{m}}^{3/2}$ and damping coefficient $D=1.45{\mathrm{e}}^5\mathrm{N}\mathrm{s}/{\mathrm{m}}^2$ .

The difference of vertical position Y between the contact approaches is shown in Figure 10. An explicit Runge–Kutta method of order 4 has been used as numerical time integrator scheme for solving penalty method. In case of the CCP approach, the semi-implicit Euler has been used. The time step is $\mathrm{\Delta }t=0.0001\mathrm{s}$ is used for both time integration methods. For this simple case, the computation time for CCP is around six times slower than for the penalty method. However, when simulating tens of thousands of contact problems, a complementarity method seems to be a better choice.^21,24 Furthermore, the CCP approach does not produce stability problems due to unnecessary high eigenfrequencies which maybe be problematic in the case of the penalty method.

Figure 10.

Y coordinates and the error of the pendulum ball from CCP approach and penalty method.

Potential and kinetic energy associated with the CCP approach and penalty method have been shown separately in Figures 11 and 12. As the two figures illustrate, the coefficient of restitution is 0. A change in the sum of kinetic and potential energy exists after inelastic contact. The material damping of the contact bodies is assumed to be the reason which causes energy loss. ¹⁹ In Figure 11, the energy changes instantly because of the nonlinear behavior of the CCP approach. However, in Figure 12, the energy changes continuously with the penalty method.

Figure 11.

Energies of two-pendulum contact with CCP approach.

Figure 12.

Two-pendulum contact with penalty method. Contact stiffness coefficient $K=3.6e^5\mathrm{N}/{\mathrm{m}}^{3/2}$ , damping coefficient $D=1.45e^5\mathrm{N}\mathrm{s}/{\mathrm{m}}^2$ .

Analysis of multiple pendulums

As shown in Figure 13, two pendulums of 20 bodies are fixed in the ceiling. The whole system is affected by the gravity force whose acceleration is $g=9.8\mathrm{kg}/{\mathrm{m}}^3$ , and the coefficient of friction is $\mu =0.3$ . The time step is $\mathrm{\Delta }t=0.001\mathrm{s}$ and the total time is 8 s.

Figure 13.

Initial position of the two pendulum system with multiple rigid balls.

Figure 14 displays the energies of the two pendulums when they fall down with zero initial velocity. The figure also shows the number of contacts during the process of two-pendulum contact. The sum of kinetic energy and potential energy decreases when the contact occurs between two pendulums. The change in the sum of kinetic and potential energy exists after each inelastic contact because the coefficient of the restitution is effectively zero with the cone complementarity approach.

Figure 14.

Energies and the number of contact of the whole system. The total number of the contact is 63 during an 8 s simulation.

Figure 15 shows the frames of the two pendulums from $t=1\mathrm{s}$ to $t=1.3\mathrm{s}$ . The contact bodies have been marked in the figure. In this example, neighbor bodies are falling down with gravity force and restricted by the rope inextensibility constraint.

Figure 15.

Frames of the simulation of two pendulums with inelastic collisions at $t=\left[\begin{array}{cccccc}1.804\hfill & 1.121\hfill & 1.142\hfill & 1.185\hfill & 1.230\hfill & 1.306\hfill \end{array}\right]\mathrm{s}$ .

Contact impulse $\mathbf{\gamma }\mathbf{=}\hat{\mathbf{\gamma }}\mathbf{\Delta }\mathbf{t}$ , which is the product of contact force $\hat{\gamma }$ and time step $\mathrm{\Delta }t$ is calculated from equation (31) at each time step. Figure 16 shows the normal contact impulses ${\hat{\gamma }}_n\mathrm{\Delta }t$ on body 10 and the frames of the simulation at corresponding times. Frames of the simulation at corresponding time have been shown in Figure 16.

Figure 16.

Normal contact impulse happened on body 10 during 8 s. Contact bodies and neighbor bodies at $t=[\begin{array}{cccccc}1.804\hfill & 1.099\hfill & 3.298\hfill & 3.489\hfill & 4.074\hfill & 4.122\hfill \end{array}]\mathrm{s}$ have been shown in the frames.

In Figure 16, the contact impulse between body 10 and body 20 at $t=1.084\mathrm{s}$ is much greater than the contact impulses with other bodies. The contact at $t=1.084\mathrm{s}$ is the first contact in the system. The generalized velocities of body 10 and body 20 are greater than the rest of bodies before contact.

Figure 17 shows the inertial force of body 10. The inertial force increases significantly when the contact takes place at the first contact when $t=1.084\mathrm{s}$ during the simulation. The inertial force is calculated as follows

$$F_{\mathrm{inertial}}=m\Vert \begin{array}{cc}{\stackrel{··}{R}}_x\hfill & {\stackrel{··}{R}}_y\hfill \end{array}\Vert$$

(57)

where $\stackrel{··}{R}=({\stackrel{·}{R}}^{(l+1)}-{\stackrel{·}{R}}^{(l)})\mathrm{\Delta }t$ .

Figure 17.

Inertial force of body 10.

The numerical simulations presented have been performed with MATLAB. These simulations show the convergence of the introduced algorithm and formulations in the case of multiple unilateral contacts.

Conclusion

The method studied in this article can solve non-smooth rigid multibody frictional contacts. This article aims for a solution to a planar system with thousands of dynamical contacts and presents a novel method for solving the CCP that appears in the time integration approach. The introduced method is able to simulate colliding rigid bodies on a large scale. In addition, this article proposed and analyzed a method to analyze the contact dynamics of two pendulums. The Coulomb friction model utilized in this article is simplified into two equations, which is simpler than Coulomb friction model for the 3D case. ¹⁸ The coefficient of restitution after contact is 0, which agrees with the conclusion made in Stewart and Trinkle. ¹⁶ It is also the reason why CCP is widely used in contact applications of powder composites, granular flows. ¹⁸ The energy change has been explained by comparing the introduced approach with penalty method.

In the future, this work can be combined with the implicit time-stepping method and elastic contacts which could be utilized in a virtual interactive environment for training, clinical therapies, military purposes, and video games.