Large-Dimensional Multibody Dynamics Simulation Using Contact Nodalization and Diagonalization

We formulate discrete-time domain dynamics at the velocity and impulse level, as is generally the case to avoid inconsistency [45]. Dynamics integration methods can be broadly classified according to explicit/implicit type and linear/non-linear type. In this work, we construct the dynamic equation in linearized form, considering the constraints in an implicit manner. This is because: 1) while targeting complex systems including flexible bodies, the explicit method has clear limitations in its stability, and 2) non-linear integration can actually be expressed as an iteration of a linear integration. The equation of motion of mechanical system under contacts in the continuous-time domain can be written as follow [14]:

where $q$ is the generalized coordinate variable¹¹1can involve representation in Euclidean space as well as orientation representation such as SO(3). Note that generalized velocity $\dot{q}$ and acceleration $\ddot{q}$ can still be expressed as $\mathbb{R}^{n}$. of system, $M(q),C(q,\dot{q})\in\mathbb{R}^{n\times n}$ are the mass, Coriolis matrix, $d\psi^{T}\in\mathbb{R}^{n}$ is the potential action, $f_{ext}\in\mathbb{R}^{n}$ is the external force, $\lambda_{c}=[\lambda_{c,1};\cdots;\lambda_{c,{n_{c}}}]\in\mathbb{R}^{3n_{c}}$ is the contact impulse with $\lambda_{c,m}=[\lambda_{n,m};\lambda_{t_{1},m};\lambda_{t_{2},m}]\in\mathbb{R}^{3}$, $J_{c}=[J_{c,1};\cdots;J_{c,{n_{c}}}]\in\mathbb{R}^{3n_{c}\times n}$ is the contact Jacobian, and subscripts $n,t_{1},t_{2}$ denote the normal and tangential directions. Then we perform the discretization of the dynamics as

where $k$ denotes the time step index, $M_{k}=M(q_{k})$, $C_{k}=C(q_{k},v_{k})$, $t_{k}$ is the step size, $v_{k}\in\mathbb{R}^{n}$ is the velocity, and $\hat{v}_{k}$ is the representative velocity of each time step. We derive potential action from the passivity relation presented in [46] i.e.,

with the second-order approximation of exact potential energy deviation. Here, potential function may be non-convex therefore the Hessian term may not be symmetric positive definite. Yet, some common approximations can be adopted to solve the issue in a compact way. Consider the following constraint potential form:

where $e(q)\in\mathbb{R}^{n_{e}}$ is the constraint error and $K(q)\in\mathbb{R}^{n_{e}\times n_{e}}$ is the symmetric positive definite gain matrix, and $n_{s}$ is the constraint dimension. This form (7) is very versatile, as it can represent almost all types of constraint including simple spring, co-rotational finite element model [47] and even hyper-elastic material with generalized compliance model [37]. Then we can write as follow using outer product approximation similar to [48, 37, 43]:

where $J_{e}(q)\in\mathbb{R}^{n_{e}\times n}$ is the constraint Jacobian (i.e., $\frac{\partial e}{\partial q}$) and $E(q)\in\mathbb{R}^{n\times n}$ is the symmetric positive definite damping matrix. o maintain the exact energy relation along these approximations (or at least, stability preserved), determining $E(q)$ is another meaningful subject, as it is necessary to find an appropriate energy dissipation to preserve the passivity of the system. Yet in this paper, we do not delve deep into this issue and apply the following two simple policies: 1) user-defined constant damping matrix or 2) symmetric positive definite projection of geometric stiffness matrix [48] i.e.,

where $\text{proj}_{\mathbb{S}_{n}^{+}}$ is the projection to $n\times n$ symmetric positive definite matrix manifold [49]. One of classical ways to compute $\text{proj}_{\mathbb{S}_{n}^{+}}$ uses singular value decomposition [50], however, to reduce the computation time for this, we construct a diagonal matrix where each element is the sum of the absolute values of the elements in each column of the original matrix.

One of the widely used potentials that cannot be represented by (7) is gravity potential, however, its Hessian can be ignored by dropping out the derivative of the $J_{e}(q)$ as above. Finally, substituting (8) to (4), dynamics can be represented as in form of (1) with

where $\hat{M}_{k}=2t_{k}^{-1}M_{k}$. The structure of (9) is very similar to the one used in [48, 46]. Although it is based on the linearization of non-linear potential action, constraints are considered in implicit manner, therefore applicable to complex multibody systems including rigid and deformable bodies. Also here, note that $A_{k}$ is always a symmetric positive definite matrix.

The form of (1) is also applicable to the general integration method. For example, coordinate transform using $M(q)^{\frac{1}{2}}$ introduced in [46] can be utilized to improve the passivity property of the dynamics, as it maintains the linearized form with a symmetric positive definite property of $A_{k}$. For the cases of non-linear integration, that takes constraints into account in a completely implicit manner without linearization, we can write the equation in the form of

where $g:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is the non-linear function. For instance, for fully-implicit Euler integration [51], $g(\hat{v}_{k})$ is

with $\hat{v}_{k}=v_{k+1},~{}\hat{q}_{k}=\text{update}(q_{k},\hat{v}_{k},t_{k})$. Similarly, implicit midpoint integration [52], variational integration [53] can be represented. Then solving this equation using the Newton method is equivalent to

and (10) can be simplified as

which is same as (1). Here, $A_{k}$ is a symmetric positive definite matrix, with some quasi-Newton style approximation. Meanwhile, still there remains room for expansion to be directly applicable to the non-linear integration methods, especially those that can be expressed in an optimization form [51]. However, we will remain the part as future work.

For specific cases, algorithms such as [54] and [55] can be utilized for effective factorization or linear solving of $A_{k}$. They are not considered here as they do not apply to general multibody dynamics.

To reflect the physical properties of the contact robustly and precisely, we consider solving the constructed dynamics (1) under the following complementarity-based Signorini-Coulomb condition (SCC):

for all contact indices $m=\left\{1,\cdots,n_{c}\right\}$ where $\mu_{m}$ is the friction coefficient, $\phi_{n,m}\in\mathbb{R}^{1}$ is the additional term for penetration compensation and restitution coefficient, and $\delta_{m}\in\mathbb{R}^{1}$ is the auxiliary variable. The first line of (11) is known as the Signorini condition which prevents penetration in a collision. The rest parts correspond to the Coulomb friction condition that ensures the contact impulse contained within the friction cone set is defined as

and tangential impulse operates in the opposite direction of motion if sliding occurs. There are three possible behavior outcomes as the result of this condition, which are illustrated in Fig. 1 - open ($\lambda_{n,m}=0$), stick ($\lambda_{n,m}>0,\delta_{m}=0$), and slip ($\lambda_{n,m}>0,\delta_{m}>0$).

We start with the following definition that categorizes contacts into two types:

Let us first consider solving (15) with only Signorini condition in (11). Then we can find that it is equivalent to Karush–Kuhn–Tucker (KKT) conditions of the following velocity-level optimization problem:

where $\mathcal{C_{V}}=\{\hat{v}\mid J_{n}\hat{v}+\phi_{n}\geq 0\}$ is the feasible set of velocity with $\phi_{n}=[\phi_{n,1},\cdots,\phi_{n,m}]^{T}$. Here and hereafter, notations for $k$-th time step are omitted for simplification, however all the components are still time-varying. This means that contact problems without friction can be replaced with solving (19). One way to solve (19) is the projected gradient descent method [56] which takes the following steps:

where $l$ is the iteration loop index, $W\in\mathbb{R}^{n\times n}$ is the symmetric positive definite step size matrix for gradient descent. Here, $W$ can alternatively interpreted as an inverse of surrogate dynamics matrix, which will be a key component of our derivation. Also, $\Pi(\cdot)$ is defined as

which is the projection from $x$ on the convex set $\mathcal{C_{V}}$, with respect to weighted Euclidean distance. Since $A$ is positive definite and $\mathcal{C_{V}}$ is convex, the projected gradient descent method will well converge to the solution of (19) if the problem is feasible. However, since (11) also include the impulse constraint like friction direction constraint, it is not enough to deal with the generic contact problem. Instead, we modify the projection step (21) as follow:

Noticing the similarity between (22) and (1), (22) can be considered as a surrogate dynamics, which is constructed from $W$ and $\hat{v}^{*}$. Then the update step can be interpreted as solving the contact problem with respect to this surrogate dynamics. Now recalling the form of (2), $\lambda_{c}$ is equivalent to

where $\Gamma_{c}=J_{c}WJ_{c}^{T}\in\mathbb{R}^{3n_{c}\times 3n_{c}}$ can be interpreted as a surrogate Delassus operator, $\eta_{c}=J_{c}\hat{v}^{*}\in\mathbb{R}^{3n_{c}}$, and SOL denotes the solution that $\lambda_{c}$ and $J_{c}\hat{v}^{l+1}=\Gamma_{c}\lambda_{c}+\eta_{c}$ satisfies SCC. After solving this contact problem (23), update step is performed as

Unlike the projection step in projected gradient descent that only enforces velocity-level constraint, this new velocity fixed point iteration (V-FPI, i.e., iteration of (20), (23) and (24)) directly achieves the contact condition (11).

In our V-FPI process above, $\hat{v}^{*}$ computing step (20) and $\hat{v}^{l+1}$ updating step (24) are simple processes, yet the part that solves surrogate contact problem (23) may be complicated. However, we find that selecting an appropriate surrogate dynamics matrix (i.e., $W$) under contact nodalization can make our surrogate Delassus operator $\Gamma_{c}=J_{c}WJ_{c}^{T}$ to be a diagonal matrix for all situations and the process can be drastically simplified. The proposition below demonstrates how it works.

Prop. 3 demonstrates that under a certain structure of $W$, extracting the components of $W$ is sufficient to construct the surrogate Delassus operator, without any matrix-matrix multiplication or factorization. Therefore, we can obtain significant advantages in time and memory compared to the original Delassus operator (i.e. $J_{c}A^{-1}J_{c}^{T}$) assembly. Also, the simple structure of the Delassus operator is significantly advantageous for the contact solving process (see Sec. III-D).

With this contact “diagonalization” method, we have built the basic structure of our simulation algorithm - COND, which is summarized in Alg. 1.

Now the remaining part is how to solve (23). In typical solvers, although the contact conditions (11) are constraints that are independent with each contact, coupling between each contact still exists since the Delassus operator is dense. These coupling terms make the global iteration and relaxation process essential for contact solvers when dealing with multi-contact situations [18, 17]. However, in COND, the diagonalized property of $\Gamma_{c}$ established in Sec. III-C resolves this problem. Convenience induced from the property of $\Gamma_{c}$ in the contact solving process can be summarized as follows:

• •

  Coupling relaxation through global iteration (e.g., Gauss-Seidel) is unnecessary since each contact situation is completely decoupled. Furthermore, completely parallel computation is possible for each contact.

• •

  Even for each single contact problem, coupling between normal and tangential impulse does not exist and it leads the exact solution to be obtained more simply without the need of numerical methods such as bisection [18], Newton [19], etc.

Overall, our contact solver consists only of a single linear solving of $-\Gamma_{c}^{-1}$ (which is very simple since $\Gamma_{c}$ is a diagonal matrix) with parallelized local projection step on $\mathcal{C}_{\lambda}$. We summarize this contact solver in Alg. 2.

The friction cone projection step (denoted as ProjectFC in Alg. 2) is necessary to enforce the output of the contact solver to satisfy the SCC. If temporary value (i.e., $\lambda_{c,m}^{*}$ in Alg. 2) is inside $\mathcal{C}_{\lambda}$, ProjectFC simply yields the same value as input. Otherwise, $\lambda_{c,m}^{*}$ is projected on the surface of $\mathcal{C}_{\lambda}$. In this paper, we use two projection schemes - “strict” and “proximal” operator as illustrated in Fig 3. Each can be written as

where $\mathcal{C}_{\lambda_{n,m}}$ is the cross-section of $\mathcal{C}_{\lambda}$ cut vertically from $\lambda_{n,m}$. It can be easily verified that for both strict operator and proximal operator, we can obtain the solution in a very simple analytic form. Also note that $\phi_{c}$ is only related to normal component, therefore $\phi_{c,m}=[\phi_{n,m};0;0]$.

We figure out that the two projection schemes have the following trade-offs in our solver: 1) the result of the strict operator exactly satisfies the SCC condition (11), yet sufficient condition for convergence is not always feasible²²2From the fundamental limitation of SCC. In practice, we do not experience the convergence problem. For detailed discussions and results, see Sec. IV,V. and 2) proximal operator can always guarantee the convergence of the solver, yet the solution does not exactly satisfy (11). Also, the result from the proximal operator is equivalent to the convex contact model proposed in [29, 30], therefore having a unique solution under strong convexity.

Contact diagonalization allows COND to easily extend with various contact models. Here we present some examples.

Summarizing the preceding contents, it is established that all components of V-FPI in COND (Alg. 1) are made up of simple scalar algebraic operations and matrix-vector multiplication. More precisely, the first multiplication of $W$ in (22) and (24) has $\mathcal{O}(n)$ complexity for both time and space since $W$ is a diagonal matrix. One-shot contact solving process (Alg. 2) also possesses $\mathcal{O}(n)$ complexity since the construction of $\Gamma_{c}$ does not involve any multiplication, but only element extraction from $W$, while contact solving is interpreted as $n_{c}$ (usually proportional to $n$) number of simple parallelizable operations. Also, since $J_{c}$ can be treated as just a stack of SO(3) matrices, matrix-vector multiplication on $J_{c}$ has also $\mathcal{O}(n)$ complexity. Therefore, the only part that has over linear complexity is a computation of $A\hat{v}^{l}$. Matrix-vector multiplication requires $\mathcal{O}(n^{2})$ complexity if $A$ is dense, but in many cases, especially for deformable body parts, $A$ contains many $0$ components. As a result, COND shows the complexity near $\mathcal{O}(n)$ - see Sec. V. Note that despite $A$ is sparse, $A^{-1}$ is generally fully dense, therefore computation of Delassus operator $J_{c}A^{-1}J_{c}^{T}$ is still has a complexity near $O(n^{3})$.

We find that Chebyshev acceleration [61] can be utilized as an efficient plug-in to accelerate the V-FPI in COND. For iterative linear solver to solve $A\hat{v}=b$ defined as

where $A=A_{1}-A_{2}$, Chebyshev acceleration method proposes iteration scheme as

with

where $P$ denotes the Chebyshev polynomial and $\varrho$ is the spectral radius of $A_{1}^{-1}A_{2}$. Looking closely at Alg. 1, we can find that our solver is quite similar to the linear solving scheme with $A_{1}=-W^{-1},A_{2}=W^{-1}-A$. Therefore, it can be expected that applying Chebyshev acceleration to COND will be effective. The algorithm is clarified in Alg. 3.